[cig-commits] r5111 - in short/3D/PyLith/trunk: libsrc/feassemble
playpen/integrate/src
knepley at geodynamics.org
knepley at geodynamics.org
Fri Oct 27 18:21:54 PDT 2006
Author: knepley
Date: 2006-10-27 18:21:54 -0700 (Fri, 27 Oct 2006)
New Revision: 5111
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity3D.hh
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh
short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc
Log:
Added code for updating a global matrix
- This is untested
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2006-10-28 00:15:25 UTC (rev 5110)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2006-10-28 01:21:54 UTC (rev 5111)
@@ -19,7 +19,7 @@
#if !defined(pylith_feassemble_integrator_hh)
#define pylith_feassemble_integrator_hh
-#include <Mesh.hh> // USES Mesh
+#include <petscmesh.h> // USES Mesh
#include "pylith/utils/petscfwd.h" // USES PetscMat
namespace pylith {
@@ -75,6 +75,7 @@
*/
virtual
void integrate(PetscMat* mat,
+ const ALE::Obj<real_section_type>& fieldIn,
const ALE::Obj<real_section_type>& coordinates) = 0;
/** Set quadrature for integrating finite-element quantities.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity3D.cc 2006-10-28 00:15:25 UTC (rev 5110)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity3D.cc 2006-10-28 01:21:54 UTC (rev 5111)
@@ -50,6 +50,7 @@
// Compute matrix associated with operator.
void
pylith::feassemble::IntegratorElasticity3D::integrate(PetscMat* mat,
+ const ALE::Obj<ALE::Mesh::real_section_type>& fieldIn,
const ALE::Obj<ALE::Mesh::real_section_type>& coordinates)
{ // integrate
} // integrate
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity3D.hh 2006-10-28 00:15:25 UTC (rev 5110)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity3D.hh 2006-10-28 01:21:54 UTC (rev 5111)
@@ -60,6 +60,7 @@
* @param coordinates Field of cell vertex coordinates
*/
void integrate(PetscMat* mat,
+ const ALE::Obj<ALE::Mesh::real_section_type>& fieldIn,
const ALE::Obj<ALE::Mesh::real_section_type>& coordinates);
// PROTECTED METHODS ////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc 2006-10-28 00:15:25 UTC (rev 5110)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc 2006-10-28 01:21:54 UTC (rev 5111)
@@ -115,21 +115,29 @@
void
pylith::feassemble::IntegratorInertia::integrate(
PetscMat* mat,
+ const ALE::Obj<real_section_type>& fieldIn,
const ALE::Obj<real_section_type>& coordinates)
{ // integrate
assert(0 != mat);
assert(0 != _quadrature);
+ PetscErrorCode ierr;
- // Setup symmetric, sparse matrix
- // :TODO: ADD STUFF HERE
-
// Get information about section
const topology_type::patch_type patch = 0;
const ALE::Obj<topology_type>& topology = coordinates->getTopology();
const ALE::Obj<topology_type::label_sequence>& cells =
topology->heightStratum(patch, 0);
const topology_type::label_sequence::iterator cellsEnd = cells->end();
+ const ALE::Obj<ALE::Mesh::order_type>& globalOrder = ALE::New::NumberingFactory<topology_type>::singleton(topology->debug())->getGlobalOrder(topology, patch, "default", fieldIn->getAtlas());
+ // Setup symmetric, sparse matrix
+ int localSize = globalOrder->getLocalSize();
+ int globalSize = globalOrder->getGlobalSize();
+ ierr = MatCreate(topology->comm(), mat);
+ ierr = MatSetSizes(*mat, localSize, localSize, globalSize, globalSize);
+ ierr = MatSetFromOptions(*mat);
+ ierr = preallocateMatrix(topology, fieldIn->getAtlas(), globalOrder, *mat);
+
// Allocate matrix for cell values (if necessary)
_initCellMatrix();
@@ -175,7 +183,7 @@
throw std::runtime_error("Logging PETSc flops failed.");
// Assemble cell contribution into sparse matrix
- // :TODO: ADD STUFF HERE
+ ierr = updateOperator(*mat, fieldIn, globalOrder, *cellIter, _cellMatrix, ADD_VALUES);
} // for
} // integrate
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh 2006-10-28 00:15:25 UTC (rev 5110)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh 2006-10-28 01:21:54 UTC (rev 5111)
@@ -66,8 +66,8 @@
* @param coordinates Field of cell vertex coordinates
*/
void integrateAction(const ALE::Obj<ALE::Mesh::real_section_type>& fieldOut,
- const ALE::Obj<ALE::Mesh::real_section_type>& fieldIn,
- const ALE::Obj<ALE::Mesh::real_section_type>& coordinates);
+ const ALE::Obj<real_section_type>& fieldIn,
+ const ALE::Obj<real_section_type>& coordinates);
/** Compute matrix associated with operator.
*
@@ -75,7 +75,8 @@
* @param coordinates Field of cell vertex coordinates
*/
void integrate(PetscMat* mat,
- const ALE::Obj<real_section_type>& coordinates);
+ const ALE::Obj<real_section_type>& fieldIn,
+ const ALE::Obj<real_section_type>& coordinates);
/** Compute lumped mass matrix.
*
Modified: short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc
===================================================================
--- short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc 2006-10-28 00:15:25 UTC (rev 5110)
+++ short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc 2006-10-28 01:21:54 UTC (rev 5111)
@@ -220,7 +220,7 @@
const Obj<topology_type>& topology = X->getTopology();
const Obj<topology_type::label_sequence>& elements = topology->heightStratum(patch, 0);
const topology_type::label_sequence::iterator end = elements->end();
- const double C = 1.0;
+ const double *C;
value_type detJ;
for(topology_type::label_sequence::iterator e_iter = elements->begin(); e_iter != end; ++e_iter) {
@@ -241,10 +241,12 @@
for(int d = 0; d < 3; d++) {
for(int e = 0; e < 3; e++) {
for(int f = 0; f < 3; f++) {
- _elemMatrix[((i*3+c)*NUM_BASIS_FUNCTIONS + j)*3+d] += 0.25*C*(_testWork[d] + _testWork[c])*(_basisWork[f] +_basisWork[e])*_weights[q]*detJ;
+ _elemMatrix[((i*3+c)*NUM_BASIS_FUNCTIONS + j)*3+d] += 0.25*C[q*81+(((c*3+d)*3)+e)*3)+f]*(_testWork[d] + _testWork[c])*(_basisWork[f] +_basisWork[e])*_weights[q]*detJ;
}
}
- _elemMatrix[((i*3+c)*NUM_BASIS_FUNCTIONS + j)*3+d] += _basis[q*NUM_BASIS_FUNCTIONS+i]*_basis[q*NUM_BASIS_FUNCTIONS+j]*_weights[q]*detJ;
+ if (c == d) {
+ _elemMatrix[((i*3+c)*NUM_BASIS_FUNCTIONS + j)*3+d] += _basis[q*NUM_BASIS_FUNCTIONS+i]*_basis[q*NUM_BASIS_FUNCTIONS+j]*_weights[q]*detJ;
+ }
}
}
}
More information about the cig-commits
mailing list