[cig-commits] r5105 - short/3D/PyLith/trunk/libsrc/feassemble
brad at geodynamics.org
brad at geodynamics.org
Thu Oct 26 21:49:03 PDT 2006
Author: brad
Date: 2006-10-26 21:49:03 -0700 (Thu, 26 Oct 2006)
New Revision: 5105
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia3D.hh
Log:
More work on inertial operator.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia3D.cc 2006-10-26 23:37:16 UTC (rev 5104)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia3D.cc 2006-10-27 04:49:03 UTC (rev 5105)
@@ -19,6 +19,7 @@
#include "petscmat.h" // USES PetscMat
#include <assert.h> // USES assert()
+#include <stdexcept> // USES std::runtime_error
// ----------------------------------------------------------------------
// Constructor
@@ -49,12 +50,14 @@
{ // integrateAction
assert(0 != _quadrature);
+ // Get information about section
const topology_type::patch_type patch = 0;
const ALE::Obj<topology_type>& topology = fieldIn->getTopology();
const ALE::Obj<topology_type::label_sequence>& cells =
topology->heightStratum(patch, 0);
const topology_type::label_sequence::iterator cellsEnd = cells->end();
+ // Allocate vector for cell values (if necessary)
_initCellVector();
for (topology_type::label_sequence::iterator cellIter=cells->begin();
@@ -70,6 +73,7 @@
const real_section_type::value_type* fieldInCell =
fieldIn->restrict(patch, *cellIter);
+ // Get cell geometry information
const int numQuadPts = _quadrature->numQuadPts();
const double* basis = _quadrature->basis();
const double* quadPts = _quadrature->quadPts();
@@ -78,18 +82,27 @@
const int numBasis = _quadrature->numCorners();
const int spaceDim = _quadrature->spaceDim();
+ // FIX THIS
+ // Hardwire mass density
+ const double density = 1.0;
+
+ // Compute action for cell
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad]*jacobianDet[iQuad];
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density;
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const int iC = iBasis*spaceDim;
+ const int iBlock = iBasis * spaceDim;
double val = wt*basis[iQ+iBasis];
for (int jBasis=0; jBasis < numBasis; ++jBasis) {
val *= basis[iQ+jBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
- _cellVector[iC+iDim] += val*fieldInCell[iC+iDim];
+ _cellVector[iBlock+iDim] += val * fieldInCell[iBlock+iDim];
} // for
} // for
} // for
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(2+numBasis*(1+numBasis*(1+2*spaceDim))));
+ if (err)
+ throw std::runtime_error("Couldn't log PETSc flops.");
// Assemble cell contribution into field
fieldOut->updateAdd(patch, *cellIter, _cellVector);
@@ -99,10 +112,140 @@
// ----------------------------------------------------------------------
// Compute matrix associated with operator.
void
-pylith::feassemble::IntegratorInertia3D::integrate(PetscMat* mat,
- const ALE::Obj<real_section_type>& coordinates)
+pylith::feassemble::IntegratorInertia3D::integrate(
+ PetscMat* mat,
+ const ALE::Obj<real_section_type>& coordinates)
{ // integrate
+ assert(0 != mat);
+ assert(0 != _quadrature);
+
+ // 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();
+
+ // Allocate matrix for cell values (if necessary)
+ _initCellMatrix();
+
+ for (topology_type::label_sequence::iterator cellIter=cells->begin();
+ cellIter != cellsEnd;
+ ++cellIter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(coordinates, *cellIter);
+
+ // Reset element matrix to zero
+ _resetCellMatrix();
+
+ // Get cell geometry information
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double* basis = _quadrature->basis();
+ const double* quadPts = _quadrature->quadPts();
+ const double* quadWts = _quadrature->quadWts();
+ const double* jacobianDet = _quadrature->jacobianDet();
+ const int numBasis = _quadrature->numCorners();
+ const int spaceDim = _quadrature->spaceDim();
+
+ // FIX THIS
+ // Hardwire mass density
+ const double density = 1.0;
+
+ // Integrate cell
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density;
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ double val = wt*basis[iQ+iBasis];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const int jBlock = jBasis * spaceDim;
+ val *= basis[iQ+jBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _cellMatrix[(iBlock+iDim)*(numBasis*spaceDim)+jBlock+iDim] += val;
+ } // for
+ } // for
+ } // for
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(2+numBasis*(1+numBasis*(1+spaceDim))));
+ if (err)
+ throw std::runtime_error("Couldn't log PETSc flops.");
+
+ // Assemble cell contribution into sparse matrix
+ // STUFF GOES HERE
+ } // for
} // integrate
+// ----------------------------------------------------------------------
+// Compute lumped matrix associated with operator.
+void
+pylith::feassemble::IntegratorInertia3D::integrateLumped(
+ const ALE::Obj<real_section_type>& fieldOut,
+ const ALE::Obj<real_section_type>& coordinates)
+{ // integrateLumped
+ assert(0 != _quadrature);
+ // 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();
+
+ // Allocate matrix for cell values (if necessary)
+ _initCellVector();
+
+ for (topology_type::label_sequence::iterator cellIter=cells->begin();
+ cellIter != cellsEnd;
+ ++cellIter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(coordinates, *cellIter);
+
+ // Reset element matrix to zero
+ _resetCellVector();
+
+ // Get cell geometry information
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double* basis = _quadrature->basis();
+ const double* quadPts = _quadrature->quadPts();
+ const double* quadWts = _quadrature->quadWts();
+ const double* jacobianDet = _quadrature->jacobianDet();
+ const int numBasis = _quadrature->numCorners();
+ const int spaceDim = _quadrature->spaceDim();
+
+ // FIX THIS
+ // Hardwire mass density
+ const double density = 1.0;
+
+ // Compute lumped mass matrix for cell
+ double sumdiag = 0;
+ double diagScale = 0;
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density;
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ const double val = wt*basis[iQ+iBasis]*basis[iQ+iBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _cellVector[iBlock+iDim] += val;
+ sumdiag += val;
+ } // for
+ diagScale += numBasis*wt;
+ } // for
+ diagScale /= sumdiag*numBasis;
+
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ for (int iDim=0; iDim < numBasis; ++iDim)
+ _cellVector[iBlock+iDim] *= diagScale;
+ } // for
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(2+numBasis*2) + numBasis);
+ if (err)
+ throw std::runtime_error("Couldn't log PETSc flops.");
+
+ // Assemble cell contribution into field
+ fieldOut->updateAdd(patch, *cellIter, _cellVector);
+ } // for
+} // integrateLumped
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia3D.hh 2006-10-26 23:37:16 UTC (rev 5104)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia3D.hh 2006-10-27 04:49:03 UTC (rev 5105)
@@ -60,8 +60,16 @@
* @param coordinates Field of cell vertex coordinates
*/
void integrate(PetscMat* mat,
- const ALE::Obj<ALE::Mesh::real_section_type>& coordinates);
+ const ALE::Obj<real_section_type>& coordinates);
+ /** Compute lumped mass matrix.
+ *
+ * @param fieldOut Lumped mass matrix as field
+ * @param coordinates Field of cell vertex coordinates
+ */
+ void integrateLumped(const ALE::Obj<real_section_type>& fieldOut,
+ const ALE::Obj<real_section_type>& coordinates);
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
More information about the cig-commits
mailing list