[cig-commits] r5973 - in short/3D/PyLith/trunk: libsrc
libsrc/feassemble unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Tue Feb 6 22:31:09 PST 2007
Author: brad
Date: 2007-02-06 22:31:09 -0800 (Tue, 06 Feb 2007)
New Revision: 5973
Added:
short/3D/PyLith/trunk/libsrc/feassemble/DynExplicitElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/DynExplicitElasticity.hh
short/3D/PyLith/trunk/libsrc/feassemble/DynExplicitElasticity.icc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorDynExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorDynExplicit.hh
Modified:
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
Log:
Started reworking integrator layout. Created skeletal IntegratorDynExplicit and DynExplicitElasticity. Removed integrator from unit tests.
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2007-02-07 03:49:38 UTC (rev 5972)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2007-02-07 06:31:09 UTC (rev 5973)
@@ -20,9 +20,8 @@
lib_LTLIBRARIES = libpylith.la
libpylith_la_SOURCES = \
- feassemble/Integrator.cc \
- feassemble/IntegratorElasticity3D.cc \
- feassemble/IntegratorInertia.cc \
+ feassemble/DynExplicitElasticity.cc \
+ feassemble/IntegratorDynExplicit.cc \
feassemble/Quadrature.cc \
feassemble/Quadrature1D.cc \
feassemble/Quadrature1Din2D.cc \
Added: short/3D/PyLith/trunk/libsrc/feassemble/DynExplicitElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/DynExplicitElasticity.cc 2007-02-07 03:49:38 UTC (rev 5972)
+++ short/3D/PyLith/trunk/libsrc/feassemble/DynExplicitElasticity.cc 2007-02-07 06:31:09 UTC (rev 5973)
@@ -0,0 +1,359 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "DynExplicitElasticity.hh" // implementation of class methods
+
+#include "Quadrature.hh" // USES Quadrature
+
+#include "petscmat.h" // USES PetscMat
+#include "spatialdata/spatialdb/SpatialDB.hh"
+
+#include <assert.h> // USES assert()
+#include <stdexcept> // USES std::runtime_error
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::feassemble::DynExplicitElasticity::DynExplicitElasticity(void)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::DynExplicitElasticity::~DynExplicitElasticity(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Copy constructor.
+pylith::feassemble::DynExplicitElasticity::DynExplicitElasticity(const DynExplicitElasticity& i) :
+ IntegratorDynExplicit(i)
+{ // copy constructor
+} // copy constructor
+
+// ----------------------------------------------------------------------
+// Integrate residual term (b) for dynamic elasticity term for 3-D
+// finite elements.
+void
+pylith::feassemble::DynExplicitElasticity::integrateResidual(
+ const ALE::Obj<real_section_type>& fieldOut,
+ const ALE::Obj<real_section_type>& fieldInT,
+ const ALE::Obj<real_section_type>& fieldInTmdt,
+ const ALE::Obj<real_section_type>& coordinates)
+{ // integrateResidual
+ assert(0 != _quadrature);
+
+ // Get information about section
+ const topology_type::patch_type patch = 0;
+ const ALE::Obj<topology_type>& topology = fieldInT->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();
+ cellIter != cellsEnd;
+ ++cellIter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(coordinates, *cellIter);
+
+ // Reset element vector to zero
+ _resetCellVector();
+
+ // Restrict input fields to cell
+ const real_section_type::value_type* fieldInTCell =
+ fieldInT->restrict(patch, *cellIter);
+ const real_section_type::value_type* fieldInTmdtCell =
+ fieldInTmdt->restrict(patch, *cellIter);
+
+ // 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();
+
+ // Restrict material properties material database to quadrature
+ // points for this cell
+#if 0
+ // :QUESTION: Not sure where we will store the density section
+ const real_section_type::value_type* density =
+ _density->restrict(patch, *cellIter);
+#else
+ const double density = 1.0;
+#endif
+
+ // Compute action for cell
+ // :TODO: Start with inertia terms, then add elastic term
+ 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 valI = wt*basis[iQ+iBasis];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const int jBlock = jBasis * spaceDim;
+ const double val = valI * basis[iQ+jBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _cellVector[iBlock+iDim] += val * fieldInTCell[jBlock+iDim];
+ } // for
+ } // for
+ } // for
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(2+numBasis*(1+numBasis*(1+2*spaceDim))));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ // Assemble cell contribution into field
+ fieldOut->updateAdd(patch, *cellIter, _cellVector);
+ } // for
+} // integrateResidual
+
+// ----------------------------------------------------------------------
+// Compute matrix associated with operator.
+void
+pylith::feassemble::DynExplicitElasticity::integrateJacobian(
+ PetscMat* mat,
+ const ALE::Obj<real_section_type>& fieldIn,
+ const ALE::Obj<real_section_type>& coordinates)
+{ // integrateJacobian
+ assert(0 != mat);
+ assert(0 != _quadrature);
+ PetscErrorCode err;
+
+ // 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,
+ fieldIn->getName(),
+ fieldIn->getAtlas());
+
+ // Setup symmetric, sparse matrix
+ // :TODO: This needs to be moved outside Integrator object, because
+ // integrator object will be specific to cell type and material type
+ int localSize = globalOrder->getLocalSize();
+ int globalSize = globalOrder->getGlobalSize();
+ err = MatCreate(topology->comm(), mat);
+ err = MatSetSizes(*mat, localSize, localSize, globalSize, globalSize);
+ err = MatSetFromOptions(*mat);
+ err = preallocateMatrix(topology, fieldIn->getAtlas(), globalOrder, *mat);
+
+ // 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();
+
+ // Restrict material properties material database to quadrature
+ // points for this cell
+#if 0
+ // :QUESTION: Not sure where we will store the density section
+ const real_section_type::value_type* density =
+ _density->restrict(patch, *cellIter);
+#else
+ const double density = 1.0;
+#endif
+
+ // 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;
+ const double valI = wt*basis[iQ+iBasis];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const int jBlock = jBasis * spaceDim;
+ const double val = valI * basis[iQ+jBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _cellMatrix[(iBlock+iDim)*(numBasis*spaceDim)+jBlock+iDim] += val;
+ } // for
+ } // for
+ } // for
+ err = PetscLogFlops(numQuadPts*(2+numBasis*(1+numBasis*(1+spaceDim))));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ // Assemble cell contribution into sparse matrix
+ err = updateOperator(*mat, fieldIn, globalOrder, *cellIter, _cellMatrix,
+ ADD_VALUES);
+ } // for
+} // integrateResidual
+
+// ----------------------------------------------------------------------
+// Compute lumped matrix associated with operator.
+void
+pylith::feassemble::DynExplicitElasticity::integrateJacobian(
+ const ALE::Obj<real_section_type>& fieldOut,
+ const ALE::Obj<real_section_type>& fieldIn,
+ const ALE::Obj<real_section_type>& coordinates)
+{ // integrateJacobian
+ 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();
+
+ // Restrict material properties material database to quadrature
+ // points for this cell
+#if 0
+ // :QUESTION: Not sure where we will store the density section
+ const real_section_type::value_type* density =
+ _density->restrict(patch, *cellIter);
+#else
+ const double density = 1.0;
+#endif
+
+ // Compute lumped mass matrix for 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;
+ const double valI = wt*basis[iQ+iBasis];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const int jBlock = jBasis * spaceDim;
+ const double val = valI*basis[iQ+jBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _cellVector[iBlock+iDim] += val;
+ } // for
+ } // for
+ } // for
+
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(2+numBasis*(1+numBasis*(1+spaceDim))));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ // Assemble cell contribution into field
+ fieldOut->updateAdd(patch, *cellIter, _cellVector);
+ } // for
+} // integrateLumped
+
+// ----------------------------------------------------------------------
+// Initialize, get material property parameters from database.
+void
+pylith::feassemble::DynExplicitElasticity::initialize(
+ ALE::Obj<ALE::Mesh>& mesh,
+ spatialdata::geocoords::CoordSys* cs,
+ spatialdata::spatialdb::SpatialDB* db)
+{ // initialize
+ assert(0 != cs);
+ assert(0 != db);
+
+ typedef ALE::Mesh::real_section_type real_section_type;
+ typedef ALE::Mesh::topology_type topology_type;
+
+ // :QUESTION: Where should we store and
+ // create the density section. Perhaps one level up.
+
+ // Create density section
+ const int numQuadPts = _quadrature->numQuadPts();
+ const ALE::Mesh::int_section_type::patch_type patch = 0;
+#if 0
+ _density = mesh->getRealSection("density");
+ const int fiberDim = numQuadPts; // number of values in field per cell
+ _density->setName("density");
+ _density->setFiberDimensionByDepth(patch, 0, fiberDim);
+ _density->allocate();
+#endif
+
+ // Open database
+ db->open();
+ const int numVals = 1;
+ const char* names[numVals];
+ names[0] = "density";
+ db->queryVals(names, numVals);
+
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ 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();
+
+ // Loop over cells
+ double* cellDensity = (numQuadPts > 0) ? new double[numQuadPts] : 0;
+ for (topology_type::label_sequence::iterator cellIter=cells->begin();
+ cellIter != cellsEnd;
+ ++cellIter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(coordinates, *cellIter);
+
+ const double* quadPts = _quadrature->quadPts();
+ const int spaceDim = _quadrature->spaceDim();
+
+ // Loop over quadrature points in cell and query database
+ for (int iQuadPt=0, index=0;
+ iQuadPt < numQuadPts;
+ ++iQuadPt, index+=spaceDim)
+ // account for differences in spaceDim
+ const int err = db->query(&cellDensity[iQuadPt], numVals,
+ &quadPts[index], spaceDim, cs);
+ // Assemble cell contribution into field
+#if 0
+ _density->updateAdd(patch, *cellIter, cellDensity);
+#endif
+ } // for
+ delete[] cellDensity; cellDensity = 0;
+
+ // Close database
+ db->close();
+} // initialize
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/feassemble/DynExplicitElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/DynExplicitElasticity.hh 2007-02-07 03:49:38 UTC (rev 5972)
+++ short/3D/PyLith/trunk/libsrc/feassemble/DynExplicitElasticity.hh 2007-02-07 06:31:09 UTC (rev 5973)
@@ -0,0 +1,136 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/feassemble/DynExplicitElasticity.hh
+ *
+ * @brief Explicit time integration of dynamic elasticity equation
+ * using finite-elements.
+ *
+ * Computes terms A and b in A(t) u(t+dt) = b(u(t), u(t-dt)), where
+ * A(t) is a sparse matrix or vector, u(t+dt) is the field we want to
+ * compute at time t+dt and b is a vector that depends on the field at
+ * time t and t-dt.
+ *
+ * A = 1/(dt*dt) [M]
+ *
+ * b = 2/(dt*dt)[M]{u(t)} - 1/(dt*dt)[M]{u(t-dt)} - 1/s[K]{u(t)} + {f(t)}
+ */
+
+#if !defined(pylith_feassemble_dynexplicitelasticity_hh)
+#define pylith_feassemble_dynexplicitelasticity_hh
+
+#include <petscmesh.h> // USES Mesh
+#include "pylith/utils/petscfwd.h" // USES PetscMat
+
+#include "IntegratorDynExplicit.hh"
+
+namespace pylith {
+ namespace feassemble {
+ class DynExplicitElasticity;
+ class TestDynExplicitElasticity;
+
+ class Quadrature; // HOLDSA Quadrature
+ } // feassemble
+} // pylith
+
+class pylith::feassemble::DynExplicitElasticity :
+ public IntegratorDynExplicit
+{ // DynExplicitElasticity
+ friend class TestDynExplicitElasticity; // unit testing
+
+// PUBLIC TYPEDEFS //////////////////////////////////////////////////////
+public :
+
+ typedef ALE::Mesh Mesh;
+ typedef Mesh::topology_type topology_type;
+ typedef topology_type::point_type point_type;
+ typedef Mesh::real_section_type real_section_type;
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ DynExplicitElasticity(void);
+
+ /// Destructor
+ ~DynExplicitElasticity(void);
+
+ /// Create a copy of this object.
+ IntegratorDynExplicit* clone(void) const;
+
+ /** Integrate residual term (b) for dynamic elasticity term
+ * for 3-D finite elements.
+ *
+ * @param fieldOut Output field
+ * @param fieldInT Input field at time t
+ * @param fieldInTmdt Input field at time t-dt
+ * @param coordinates Field of cell vertex coordinates
+ */
+ void integrateResidual(const ALE::Obj<real_section_type>& fieldOut,
+ const ALE::Obj<real_section_type>& fieldInT,
+ const ALE::Obj<real_section_type>& fieldInTmdt,
+ const ALE::Obj<real_section_type>& coordinates);
+
+ /** Compute matrix (A) associated with operator.
+ *
+ * @param mat Sparse matrix
+ * @param fieldIn Input field at time t
+ * @param coordinates Field of cell vertex coordinates
+ */
+ void integrateJacobian(PetscMat* mat,
+ const ALE::Obj<real_section_type>& fieldIn,
+ const ALE::Obj<real_section_type>& coordinates);
+
+ /** Compute field (A) associated with lumped operator.
+ *
+ * @param fieldOut Output Jacobian field
+ * @param fieldIn Input field at time t
+ * @param coordinates Field of cell vertex coordinates
+ */
+ void integrateJacobian(const ALE::Obj<real_section_type>& fieldOut,
+ const ALE::Obj<real_section_type>& fieldIn,
+ const ALE::Obj<real_section_type>& coordinates);
+
+ /** Initialize, get material property parameters from database.
+ *
+ * @param mesh PETSc mesh
+ * @param cs Pointer to coordinate system of vertices
+ * @param db Pointer to spatial database with material property parameters
+ */
+ void initialize(ALE::Obj<ALE::Mesh>& mesh,
+ spatialdata::geocoords::CoordSys* cs,
+ spatialdata::spatialdb::SpatialDB* db);
+
+// PROTECTED METHODS ////////////////////////////////////////////////////
+protected :
+
+ /** Copy constructor.
+ *
+ * @param i Integrator to copy
+ */
+ DynExplicitElasticity(const DynExplicitElasticity& i);
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+ /// Not implemented
+ const DynExplicitElasticity& operator=(const DynExplicitElasticity&);
+
+}; // DynExplicitElasticity
+
+#include "DynExplicitElasticity.icc" // inline methods
+
+#endif // pylith_feassemble_dynexplicitelasticity_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/feassemble/DynExplicitElasticity.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/DynExplicitElasticity.icc 2007-02-07 03:49:38 UTC (rev 5972)
+++ short/3D/PyLith/trunk/libsrc/feassemble/DynExplicitElasticity.icc 2007-02-07 06:31:09 UTC (rev 5973)
@@ -0,0 +1,27 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_feassemble_dynexplicitelasticity_hh)
+#error "DynExplicitElasticity.icc must be included only from DynExplicitElasticity.hh"
+#else
+
+// Create a copy of this object.
+inline
+pylith::feassemble::IntegratorDynExplicit*
+pylith::feassemble::DynExplicitElasticity::clone(void) const {
+ return new DynExplicitElasticity(*this);
+} // clone
+
+#endif
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorDynExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorDynExplicit.cc 2007-02-07 03:49:38 UTC (rev 5972)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorDynExplicit.cc 2007-02-07 06:31:09 UTC (rev 5973)
@@ -0,0 +1,116 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "IntegratorDynExplicit.hh" // implementation of class methods
+
+#include "Quadrature.hh"
+#include "spatialdata/spatialdb/SpatialDB.hh"
+
+#include <assert.h> // USES assert()
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::feassemble::IntegratorDynExplicit::IntegratorDynExplicit(void) :
+ _quadrature(0),
+ _cellVector(0),
+ _cellMatrix(0)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::IntegratorDynExplicit::~IntegratorDynExplicit(void)
+{ // destructor
+ delete _quadrature; _quadrature = 0;
+ delete[] _cellVector; _cellVector = 0;
+ delete[] _cellMatrix; _cellMatrix = 0;
+} // destructor
+
+// ----------------------------------------------------------------------
+// Copy constructor
+pylith::feassemble::IntegratorDynExplicit::IntegratorDynExplicit(const IntegratorDynExplicit& i) :
+ _cellVector(0),
+ _cellMatrix(0)
+{ // copy constructor
+ if (0 != i._quadrature)
+ _quadrature = i._quadrature->clone();
+} // copy constructor
+
+// ----------------------------------------------------------------------
+// Set quadrature for integrating finite-element quantities.
+void
+pylith::feassemble::IntegratorDynExplicit::quadrature(const Quadrature* q)
+{ // quadrature
+ delete _quadrature;
+ _quadrature = (0 != q) ? q->clone() : 0;
+
+ // Deallocate cell vector and matrix since size may change
+ delete[] _cellVector; _cellVector = 0;
+ delete[] _cellMatrix; _cellMatrix = 0;
+} // quadrature
+
+// ----------------------------------------------------------------------
+// Initialize vector containing result of integration action for cell.
+void
+pylith::feassemble::IntegratorDynExplicit::_initCellVector(void)
+{ // _initCellVector
+ assert(0 != _quadrature);
+ const int size = _quadrature->spaceDim() * _quadrature->numCorners();
+ if (0 == _cellVector)
+ _cellVector = (size > 0) ? new real_section_type::value_type[size] : 0;
+ for (int i=0; i < size; ++i)
+ _cellVector[i] = 0.0;
+} // _initCellVector
+
+// ----------------------------------------------------------------------
+// Zero out vector containing result of integration actions for cell.
+void
+pylith::feassemble::IntegratorDynExplicit::_resetCellVector(void)
+{ // _resetCellVector
+ assert(0 != _quadrature);
+ const int size = _quadrature->spaceDim() * _quadrature->numCorners();
+ for (int i=0; i < size; ++i)
+ _cellVector[i] = 0.0;
+} // _resetCellVector
+
+// ----------------------------------------------------------------------
+// Initialize matrix containing result of integration for cell.
+void
+pylith::feassemble::IntegratorDynExplicit::_initCellMatrix(void)
+{ // _initCellMatrix
+ assert(0 != _quadrature);
+ const int size =
+ _quadrature->spaceDim() * _quadrature->numCorners() *
+ _quadrature->spaceDim() * _quadrature->numCorners();
+ if (0 == _cellMatrix)
+ _cellMatrix = (size > 0) ? new real_section_type::value_type[size] : 0;
+ for (int i=0; i < size; ++i)
+ _cellMatrix[i] = 0.0;
+} // _initCellMatrix
+
+// ----------------------------------------------------------------------
+// Zero out matrix containing result of integration for cell.
+void
+pylith::feassemble::IntegratorDynExplicit::_resetCellMatrix(void)
+{ // _resetCellMatrix
+ assert(0 != _quadrature);
+ const int size =
+ _quadrature->spaceDim() * _quadrature->numCorners() *
+ _quadrature->spaceDim() * _quadrature->numCorners();
+ for (int i=0; i < size; ++i)
+ _cellMatrix[i] = 0.0;
+} // _resetCellMatrix
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorDynExplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorDynExplicit.hh 2007-02-07 03:49:38 UTC (rev 5972)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorDynExplicit.hh 2007-02-07 06:31:09 UTC (rev 5973)
@@ -0,0 +1,171 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/feassemble/IntegratorDynExplicit.hh
+ *
+ * @brief Abstract base class for explicit time integration of
+ * finite-element actions.
+ *
+ * Computes terms A and b in A(t) u(t+dt) = b(u(t), u(t-dt)), where
+ * A(t) is a sparse matrix or vector, u(t+dt) is the field we want to
+ * compute at time t+dt and b is a vector that depends on the field at
+ * time t and t-dt.
+ */
+
+#if !defined(pylith_feassemble_integratordynexplicit_hh)
+#define pylith_feassemble_integratordynexplicit_hh
+
+#include <petscmesh.h> // USES Mesh
+#include "pylith/utils/petscfwd.h" // USES PetscMat
+
+namespace pylith {
+ namespace feassemble {
+ class IntegratorDynExplicit;
+ class TestIntegratorDynExplicit;
+
+ class Quadrature; // HOLDSA Quadrature
+ } // feassemble
+} // pylith
+
+namespace spatialdata {
+ namespace spatialdb {
+ class SpatialDB; // USES SpatialDB
+ } // spatialdb
+ namespace geocoords {
+ class CoordSys; // USES CoordSys
+ } // geocoords
+} // spatialdata
+
+class pylith::feassemble::IntegratorDynExplicit
+{ // Integrator
+ friend class TestIntegratorDynExplicit; // unit testing
+
+// PUBLIC TYPEDEFS //////////////////////////////////////////////////////
+public :
+
+ typedef ALE::Mesh Mesh;
+ typedef Mesh::topology_type topology_type;
+ typedef topology_type::point_type point_type;
+ typedef Mesh::real_section_type real_section_type;
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ IntegratorDynExplicit(void);
+
+ /// Destructor
+ virtual
+ ~IntegratorDynExplicit(void);
+
+ /// Create a copy of this object.
+ virtual
+ IntegratorDynExplicit* clone(void) const = 0;
+
+ /** Integrate residual term (b) for dynamic elasticity term
+ * for 3-D finite elements.
+ *
+ * @param fieldOut Output field
+ * @param fieldInT Input field at time t
+ * @param fieldInTmdt Input field at time t-dt
+ * @param coordinates Field of cell vertex coordinates
+ */
+ virtual
+ void integrateResidual(const ALE::Obj<real_section_type>& fieldOut,
+ const ALE::Obj<real_section_type>& fieldInT,
+ const ALE::Obj<real_section_type>& fieldInTmdt,
+ const ALE::Obj<real_section_type>& coordinates) = 0;
+
+ /** Compute matrix (A) associated with operator.
+ *
+ * @param mat Sparse matrix
+ * @param fieldIn Input field at time t
+ * @param coordinates Field of cell vertex coordinates
+ */
+ virtual
+ void integrateJacobian(PetscMat* mat,
+ const ALE::Obj<real_section_type>& fieldIn,
+ const ALE::Obj<real_section_type>& coordinates) = 0;
+
+ /** Compute field (A) associated with lumped operator.
+ *
+ * @param fieldOut Output Jacobian field
+ * @param fieldIn Input field at time t
+ * @param coordinates Field of cell vertex coordinates
+ */
+ virtual
+ void integrateJacobian(const ALE::Obj<real_section_type>& fieldOut,
+ const ALE::Obj<real_section_type>& fieldIn,
+ const ALE::Obj<real_section_type>& coordinates) = 0;
+
+ /** Set quadrature for integrating finite-element quantities.
+ *
+ * @param q Quadrature for integrating.
+ */
+ void quadrature(const Quadrature* q);
+
+ /** Initialize, get material property parameters from database.
+ *
+ * @param mesh PETSc mesh
+ * @param cs Pointer to coordinate system of vertices
+ * @param db Pointer to spatial database with material property parameters
+ */
+ virtual
+ void initialize(ALE::Obj<ALE::Mesh>& mesh,
+ spatialdata::geocoords::CoordSys* cs,
+ spatialdata::spatialdb::SpatialDB* db) = 0;
+
+// PROTECTED METHODS ////////////////////////////////////////////////////
+protected :
+
+ /** Copy constructor.
+ *
+ * @param i Integrator to copy
+ */
+ IntegratorDynExplicit(const IntegratorDynExplicit& i);
+
+ /// Initialize vector containing result of integration action for cell.
+ void _initCellVector(void);
+
+ /// Zero out vector containing result of integration actions for cell.
+ void _resetCellVector(void);
+
+ /// Initialize matrix containing result of integration for cell.
+ void _initCellMatrix(void);
+
+ /// Zero out matrix containing result of integration for cell.
+ void _resetCellMatrix(void);
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+ /// Not implemented
+ const IntegratorDynExplicit& operator=(const IntegratorDynExplicit&);
+
+// PROTECTED MEMBERS ////////////////////////////////////////////////////
+protected :
+
+ Quadrature* _quadrature; ///< Quadrature for integrating finite-element
+
+ /// Vector local to cell containing result of integration action
+ real_section_type::value_type* _cellVector;
+
+ /// Matrix local to cell containing result of integration
+ real_section_type::value_type* _cellMatrix;
+
+}; // IntegratorDynExplicit
+
+#endif // pylith_feassemble_integratordynexplicit_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2007-02-07 03:49:38 UTC (rev 5972)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2007-02-07 06:31:09 UTC (rev 5973)
@@ -14,11 +14,8 @@
include $(top_srcdir)/subpackage.am
subpkginclude_HEADERS = \
- Integrator.hh \
- IntegratorElasticity3D.hh \
- IntegratorElasticity3D.icc \
- IntegratorInertia.hh \
- IntegratorInertia.icc \
+ DynExplicitElasticity.hh \
+ IntegratorDynExplicit.hh \
Quadrature.hh \
Quadrature.icc \
Quadrature1D.hh \
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am 2007-02-07 03:49:38 UTC (rev 5972)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am 2007-02-07 06:31:09 UTC (rev 5973)
@@ -19,11 +19,6 @@
# Primary source files
testfeassemble_SOURCES = \
- TestIntegrator.cc \
- TestIntegratorInertia.cc \
- TestIntegratorInertia1D.cc \
- TestIntegratorInertia2D.cc \
- TestIntegratorInertia3D.cc \
TestQuadrature.cc \
TestQuadrature1D.cc \
TestQuadrature1Din2D.cc \
@@ -32,6 +27,12 @@
TestQuadrature3D.cc \
test_feassemble.cc
+#TestIntegrator.cc \
+#TestIntegratorInertia.cc \
+#TestIntegratorInertia1D.cc \
+#TestIntegratorInertia2D.cc \
+#TestIntegratorInertia3D.cc
+
noinst_HEADERS = \
TestIntegrator.hh \
TestIntegratorInertia.hh \
More information about the cig-commits
mailing list