[cig-commits] r14035 - in short/3D/PyLith/branches/pylith-swig: libsrc/feassemble libsrc/topology unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Tue Feb 10 16:06:18 PST 2009
Author: brad
Date: 2009-02-10 16:06:18 -0800 (Tue, 10 Feb 2009)
New Revision: 14035
Added:
short/3D/PyLith/branches/pylith-swig/libsrc/topology/RestrictVisitor.hh
short/3D/PyLith/branches/pylith-swig/libsrc/topology/RestrictVisitor.icc
short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadratureBase.cc
short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadratureBase.hh
Modified:
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.cc
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.hh
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.icc
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature1D.cc
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature1D.hh
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature1D.icc
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature3D.cc
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature3D.hh
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.cc
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.hh
short/3D/PyLith/branches/pylith-swig/libsrc/topology/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/topology/topologyfwd.hh
short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/Makefile.am
short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadrature.cc
short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadrature.hh
Log:
Worked on quadrature.
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Makefile.am 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Makefile.am 2009-02-11 00:06:18 UTC (rev 14035)
@@ -39,10 +39,12 @@
QuadratureBase.icc \
Quadrature.hh \
Quadrature.icc \
+ Quadrature.cc \
Quadrature0D.hh \
Quadrature0D.icc \
Quadrature1D.hh \
Quadrature1D.icc \
+ Quadrature1D.cc \
Quadrature1Din2D.hh \
Quadrature1Din2D.icc \
Quadrature1Din3D.hh \
@@ -53,6 +55,7 @@
Quadrature2Din3D.icc \
Quadrature3D.hh \
Quadrature3D.icc \
+ Quadrature3D.cc \
feassemblefwd.hh
noinst_HEADERS =
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.cc 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.cc 2009-02-11 00:06:18 UTC (rev 14035)
@@ -15,6 +15,7 @@
#include "Quadrature.hh" // implementation of class methods
#include "CellGeometry.hh" // USES CellGeometry
+#include "pylith/topology/Field.hh" // HOLDSA Field
#include <cstring> // USES memcpy()
#include <cassert> // USES assert()
@@ -23,354 +24,195 @@
// ----------------------------------------------------------------------
// Constructor
-pylith::feassemble::Quadrature::Quadrature(void) :
- _minJacobian(0),
- _cellDim(0),
- _numBasis(0),
- _numQuadPts(0),
- _spaceDim(0),
- _geometry(0),
- _precomputed(false),
+template<typename mesh_type>
+pylith::feassemble::Quadrature<mesh_type>::Quadrature(void) :
+ _quadPtsField(0),
+ _jacobianField(0),
+ _jacobianDetField(0),
+ _basisDerivField(0),
_checkConditioning(false)
{ // constructor
- _quadPtsPre = new real_section_type(PETSC_COMM_WORLD);
- _jacobianPre = new real_section_type(PETSC_COMM_WORLD);
- _jacobianDetPre = new real_section_type(PETSC_COMM_WORLD);
- _jacobianInvPre = new real_section_type(PETSC_COMM_WORLD);
- _basisDerivPre = new real_section_type(PETSC_COMM_WORLD);
} // constructor
// ----------------------------------------------------------------------
// Destructor
-pylith::feassemble::Quadrature::~Quadrature(void)
+template<typename mesh_type>
+pylith::feassemble::Quadrature<mesh_type>::~Quadrature(void)
{ // destructor
- delete _geometry; _geometry = 0;
+ delete _quadPtsField; _quadPtsField = 0;
+ delete _jacobianField; _jacobianField = 0;
+ delete _jacobianDetField; _jacobianDetField = 0;
+ delete _basisDerivField; _basisDerivField = 0;
} // destructor
// ----------------------------------------------------------------------
// Copy constructor
-pylith::feassemble::Quadrature::Quadrature(const Quadrature& q) :
- _minJacobian(q._minJacobian),
- _quadPtsRef(q._quadPtsRef),
- _quadPts(q._quadPts),
- _quadWts(q._quadWts),
- _basis(q._basis),
- _basisDerivRef(q._basisDerivRef),
- _basisDeriv(q._basisDeriv),
- _jacobian(q._jacobian),
- _jacobianInv(q._jacobianInv),
- _jacobianDet(q._jacobianDet),
- _cellDim(q._cellDim),
- _numBasis(q._numBasis),
- _numQuadPts(q._numQuadPts),
- _spaceDim(q._spaceDim),
- _geometry(0),
- _precomputed(q._precomputed),
+template<typename mesh_type>
+pylith::feassemble::Quadrature<mesh_type>::Quadrature(const Quadrature& q) :
+ QuadratureBase(q),
+ _quadPtsField(0),
+ _jacobianField(0),
+ _jacobianDetField(0),
+ _basisDerivField(0),
_checkConditioning(q._checkConditioning)
{ // copy constructor
- if (0 != q._geometry)
- _geometry = q._geometry->clone();
- _quadPtsPre = q._quadPtsPre;
- _jacobianPre = q._jacobianPre;
- _jacobianDetPre = q._jacobianDetPre;
- _jacobianInvPre = q._jacobianInvPre;
- _basisDerivPre = q._basisDerivPre;
} // copy constructor
// ----------------------------------------------------------------------
-// Set basis functions and their derivatives and coordinates and
-// weights of the quadrature points.
+// Deallocate temporary storage;
+template<typename mesh_type>
void
-pylith::feassemble::Quadrature::initialize(const double* basis,
- const double* basisDerivRef,
- const double* quadPtsRef,
- const double* quadWts,
- const int cellDim,
- const int numBasis,
- const int numQuadPts,
- const int spaceDim)
-{ // initialize
- if (0 == basis ||
- 0 == basisDerivRef ||
- 0 == quadPtsRef ||
- 0 == quadWts ||
- cellDim < 0 || cellDim > 3 ||
- numBasis < 1 ||
- numQuadPts < 1 ||
- spaceDim < 1 || spaceDim > 3) {
- std::ostringstream msg;
- msg << "Incompatible values for quadrature information. Basis functions,\n"
- << "their derivatives, and coordinates and weights of quadrature\n"
- << "points must all be specified.\n"
- << "Values:\n"
- << " basis pointer: " << basis << "\n"
- << " basis derivatites pointer: " << basisDerivRef << "\n"
- << " quadrature points pointer: " << quadPtsRef << "\n"
- << " quadrature weights pointer: " << quadWts << "\n"
- << " space dimension: " << spaceDim << "\n"
- << " # basis functions: " << numBasis << "\n"
- << " # quadrature points: " << numQuadPts << "\n"
- << " dimension of coordinate space: " << spaceDim << "\n";
- throw std::runtime_error(msg.str());
- } // if
+pylith::feassemble::Quadrature<mesh_type>::clear(void)
+{ // clear
+ // Clear storage for fields
+ delete _quadPtsField; _quadPtsField = 0;
+ delete _jacobianField; _jacobianField = 0;
+ delete _jacobianDetField; _jacobianDetField = 0;
+ delete _basisDerivField; _basisDerivField = 0;
+} // clear
- if (cellDim > 0) {
- int size = numBasis * numQuadPts; assert(size > 0);
- _basis.resize(size);
- for (int i=0; i < size; ++i)
- _basis[i] = basis[i];
-
- size = numQuadPts * numBasis * cellDim; assert(size > 0);
- _basisDerivRef.resize(size);
- for (int i=0; i < size; ++i)
- _basisDerivRef[i] = basisDerivRef[i];
-
- size = numQuadPts * cellDim; assert(size > 0);
- _quadPtsRef.resize(size);
- for (int i=0; i < size; ++i)
- _quadPtsRef[i] = quadPtsRef[i];
-
- size = numQuadPts; assert(size > 0);
- _quadWts.resize(size);
- for (int i=0; i < size; ++i)
- _quadWts[i] = quadWts[i];
-
- _cellDim = cellDim;
- _numBasis = numBasis;
- _numQuadPts = numQuadPts;
- _spaceDim = spaceDim;
-
- // Allocate for Jacobian and its inverse
- size = numQuadPts * cellDim * spaceDim; assert(size > 0);
- _jacobian.resize(size);
- _jacobianInv.resize(size);
-
- // Allocate for Jacobian determinant
- size = numQuadPts; assert(size > 0);
- _jacobianDet.resize(size);
-
- // Allocate for basis derivatives (in global coordinates)
- size = numQuadPts * numBasis * spaceDim; assert(size > 0);
- _basisDeriv.resize(size);
-
- // Allocate for quad pts
- size = numQuadPts*spaceDim; assert(size > 0);
- _quadPts.resize(size);
- } else {
- if (1 != numBasis ||
- 1 != numQuadPts ||
- 1 != spaceDim) {
- std::ostringstream msg;
- msg << "0-D quadrature only works in 1-D and is limited to 1 basis "
- << "function and 1 quadrature point.\n"
- << "Values:\n"
- << " cell dimension: " << cellDim << "\n"
- << " spatial dimension: " << spaceDim << "\n"
- << " # basis functions: " << numBasis << "\n"
- << " # quadrature points: " << numQuadPts << "\n";
- throw std::runtime_error(msg.str());
- } // if
-
- int size = 1;
- _basis.resize(size);
- for (int i=0; i < size; ++i)
- _basis[i] = basis[i];
-
- size = 1;
- _basisDerivRef.resize(size);
- for (int i=0; i < size; ++i)
- _basisDerivRef[i] = basisDerivRef[i];
-
- size = 1;
- _quadPtsRef.resize(size);
- for (int i=0; i < size; ++i)
- _quadPtsRef[i] = quadPtsRef[i];
-
- size = 1;
- _quadWts.resize(size);
- for (int i=0; i < size; ++i)
- _quadWts[i] = quadWts[i];
-
- _cellDim = cellDim;
- _numBasis = numBasis;
- _numQuadPts = numQuadPts;
- _spaceDim = spaceDim;
-
- // Allocate for Jacobian and its inverse
- size = 1;
- _jacobian.resize(size);
- _jacobianInv.resize(size);
-
- // Allocate for Jacobian determinant
- size = 1;
- _jacobianDet.resize(size);
-
- // Allocate for basis derivatives (in global coordinates)
- size = numQuadPts * numBasis * spaceDim; assert(size > 0);
- _basisDeriv.resize(size);
-
- // Allocate for quad pts
- size = spaceDim; assert(size > 0);
- _quadPts.resize(size);
- } // else
-} // initialize
-
// ----------------------------------------------------------------------
-// Set geometry associated with reference cell.
-void
-pylith::feassemble::Quadrature::refGeometry(CellGeometry* const geometry)
-{ // refGeometry
- delete _geometry; _geometry = (0 != geometry) ? geometry->clone() : 0;
-} // refGeometry
-
-// ----------------------------------------------------------------------
-// Get geometry associated with reference cell.
-const pylith::feassemble::CellGeometry&
-pylith::feassemble::Quadrature::refGeometry(void) const
-{ // refGeometry
- assert(0 != _geometry);
- return *_geometry;
-} // refGeometry
-
-// ----------------------------------------------------------------------
// Set entries in geometry arrays to zero.
+template<typename mesh_type>
void
-pylith::feassemble::Quadrature::_resetGeometry(void)
+pylith::feassemble::Quadrature<mesh_type>::_resetGeometry(void)
{ // _resetGeometry
+ _quadPts = 0.0;
_jacobian = 0.0;
_jacobianDet = 0.0;
_jacobianInv = 0.0;
_basisDeriv = 0.0;
- _quadPts = 0.0;
} // _resetGeometry
// ----------------------------------------------------------------------
-// Check determinant of Jacobian against minimum allowable value
+template<typename mesh_type>
void
-pylith::feassemble::Quadrature::_checkJacobianDet(const double det,
- const int cell) const
-{ // _checkJacobianDet
- if (det < _minJacobian) {
- std::ostringstream msg;
- msg << "Determinant of Jacobian (" << det << ") for cell " << cell
- << " is smaller than minimum permissible value (" << _minJacobian
- << ")!\n";
- throw std::runtime_error(msg.str());
- } // if
-} // _checkJacobianDet
-
-// ----------------------------------------------------------------------
-void
-pylith::feassemble::Quadrature::resetPrecomputation()
-{ // resetPrecomputation
- _precomputed = false;
- _quadPtsPre->clear();
- _jacobianPre->clear();
- _jacobianDetPre->clear();
- _jacobianInvPre->clear();
- _basisDerivPre->clear();
-} // resetPrecomputation
-
-// ----------------------------------------------------------------------
-void
-pylith::feassemble::Quadrature::precomputeGeometry(
- const ALE::Obj<Mesh>& mesh,
- const ALE::Obj<real_section_type>& coordinates,
- const ALE::Obj<Mesh::label_sequence>& cells)
+pylith::feassemble::Quadrature<mesh_type>::computeGeometry(
+ const mesh_type& mesh,
+ const ALE::Obj<typename mesh_type::SieveMesh::label_sequence>& cells)
{ // precomputeGeometry
- if (_precomputed)
- return;
+ typedef typename mesh_type::RealSection RealSection;
- _quadPtsPre->setChart(real_section_type::chart_type(
- *std::min_element(cells->begin(), cells->end()),
- *std::max_element(cells->begin(), cells->end())+1));
- _quadPtsPre->setFiberDimension(cells, _numQuadPts*_spaceDim);
- _quadPtsPre->allocatePoint();
- _quadPtsPreV = new ALE::ISieveVisitor::RestrictVisitor<real_section_type>(
- *_quadPtsPre, _numQuadPts*_spaceDim, &_quadPts[0]);
+ const char* loggingStage = "QuadratureCreation";
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush(loggingStage);
- _jacobianPre->getAtlas()->setAtlas(_quadPtsPre->getAtlas()->getAtlas());
- _jacobianPre->getAtlas()->allocatePoint();
- _jacobianPre->setFiberDimension(cells, _numQuadPts*_cellDim*_spaceDim);
- _jacobianPre->allocatePoint();
- _jacobianPreV = new ALE::ISieveVisitor::RestrictVisitor<real_section_type>(
- *_jacobianPre, _numQuadPts*_cellDim*_spaceDim, &_jacobian[0]);
+ clear();
- _jacobianDetPre->getAtlas()->setAtlas(_quadPtsPre->getAtlas()->getAtlas());
- _jacobianDetPre->getAtlas()->allocatePoint();
- _jacobianDetPre->setFiberDimension(cells, _numQuadPts);
- _jacobianDetPre->allocatePoint();
- _jacobianDetPreV = new ALE::ISieveVisitor::RestrictVisitor<real_section_type>(*_jacobianDetPre, _numQuadPts, &_jacobianDet[0]);
+ // Allocate field and cell buffer for quadrature points
+ int fiberDim = _numQuadPts * _spaceDim;
+ _quadPts.resize(fiberDim);
+ _quadPtsField = new topology::Field<mesh_type>(mesh);
+ assert(0 != _quadPtsField);
+ _quadPtsField->newSection(cells, fiberDim);
+ _quadPtsField->allocate();
- _jacobianInvPre->setAtlas(_jacobianPre->getAtlas());
- _jacobianInvPre->setFiberDimension(cells, _numQuadPts*_cellDim*_spaceDim);
- _jacobianInvPre->allocatePoint();
- _jacobianInvPreV = new ALE::ISieveVisitor::RestrictVisitor<real_section_type>(*_jacobianInvPre, _numQuadPts*_cellDim*_spaceDim, &_jacobianInv[0]);
- //_jITag = _jacobianInvPre->copyCustomAtlas(_jacobianPre, _jTag);
+ // Get chart for reuse in other fields
+ const ALE::Obj<RealSection>& section = _quadPtsField->section();
+ assert(!section.isNull());
+ const typename RealSection::chart_type& chart = section->getChart();
- _basisDerivPre->getAtlas()->setAtlas(_quadPtsPre->getAtlas()->getAtlas());
- _basisDerivPre->getAtlas()->allocatePoint();
- _basisDerivPre->setFiberDimension(cells, _numQuadPts*_numBasis*_spaceDim);
- _basisDerivPre->allocatePoint();
- _basisDerivPreV = new ALE::ISieveVisitor::RestrictVisitor<real_section_type>(*_basisDerivPre, _numQuadPts*_numBasis*_spaceDim, &_basisDeriv[0]);
-
-#if 0
- const int ncells = cells->size();
- const int nbytes = (_numQuadPts*_spaceDim + // quadPts
- _numQuadPts*_cellDim*_spaceDim + // jacobian
- _numQuadPts*_cellDim*_spaceDim + // jacobianInv
- _numQuadPts + // jacobianDet
- _numQuadPts*_numBasis*_spaceDim // basisDeriv
- ) * ncells * sizeof(double);
+ // Allocate field and cell buffer for Jacobian at quadrature points
+ fiberDim = _numQuadPts * _cellDim * _spaceDim;
+ _jacobian.resize(fiberDim);
+ _jacobianInv.resize(fiberDim);
+ _jacobianField = new topology::Field<mesh_type>(mesh);
+ assert(0 != _jacobianField);
+ _jacobianField->newSection(chart, fiberDim);
+ _jacobianField->allocate();
- std::cout << "Quadrature::precomputeGeometry() allocating "
- << nbytes/(1024*1024) << " MB."
- << std::endl;
-#endif
+ // Allocate field and cell buffer for determinant of Jacobian at quad pts
+ fiberDim = _numQuadPts;
+ _jacobianDet.resize(fiberDim);
+ _jacobianDetField = new topology::Field<mesh_type>(mesh);
+ assert(0 != _jacobianDetField);
+ _jacobianDetField->newSection(chart, fiberDim);
+ _jacobianDetField->allocate();
+
+ // Allocate field for derivatives of basis functions at quad pts
+ fiberDim = _numQuadPts * _numBasis * _spaceDim;
+ _basisDeriv.resize(fiberDim);
+ _basisDerivField = new topology::Field<mesh_type>(mesh);
+ assert(0 != _basisDerivField);
+ _basisDerivField->newSection(chart, fiberDim);
+ _basisDerivField->allocate();
- const Mesh::label_sequence::iterator cellsEnd = cells->end();
- for(Mesh::label_sequence::iterator c_iter = cells->begin();
- c_iter != cellsEnd;
- ++c_iter) {
- this->computeGeometry(mesh, coordinates, *c_iter);
+ logger.stagePop();
- // Set coordinates of quadrature points in cell
- _quadPtsPre->updatePoint(*c_iter, &_quadPts[0]);
+#if defined(ALE_MEM_LOGGING)
+ std::cout
+ << loggingStage << ": "
+ << logger.getNumAllocations(loggingStage)
+ << " allocations " << logger.getAllocationTotal(loggingStage)
+ << " bytes"
+ << std::endl
+ << loggingStage << ": "
+ << logger.getNumDeallocations(loggingStage)
+ << " deallocations " << logger.getDeallocationTotal(loggingStage)
+ << " bytes"
+ << std::endl;
+#endif
- // Set Jacobian at quadrature points in cell
- _jacobianPre->updatePoint(*c_iter, &_jacobian[0]);
+ typedef typename mesh_type::SieveMesh::label_sequence label_sequence;
+ typedef ALE::ISieveVisitor::RestrictVisitor<RealSection> RealSectionVisitor;
- // Set determinant of Jacobian at quadrature points in cell
- _jacobianDetPre->updatePoint(*c_iter, &_jacobianDet[0]);
+ const typename label_sequence::iterator cellsEnd = cells->end();
+ assert(0 != _geometry);
+ const int numCorners = _geometry->numCorners();
+ const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ RealSectionVisitor coordsVisitor(coordinates, numCorners*_spaceDim);
- // Set inverse of Jacobian at quadrature points in cell
- _jacobianInvPre->updatePoint(*c_iter, &_jacobianInv[0]);
+ const ALE::Obj<RealSection>& quadPtsSection = _quadPtsField->section();
+ const ALE::Obj<RealSection>& jacobianSection = _jacobianField->section();
+ const ALE::Obj<RealSection>& jacobianDetSection =
+ _jacobianDetField->section();
+ const ALE::Obj<RealSection>& basisDerivSection = _basisDerivField->section();
- // Set derivatives of basis functions with respect to global
- _basisDerivPre->updatePoint(*c_iter, &_basisDeriv[0]);
+ for(typename label_sequence::iterator c_iter = cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ const double* cellVertexCoords = coordsVisitor.getValues();
+ assert(0 != cellVertexCoords);
+ computeGeometry(cellVertexCoords, _spaceDim, *c_iter);
+
+ // Update fields with cell data
+ quadPtsSection->updatePoint(*c_iter, &_quadPts[0]);
+ jacobianSection->updatePoint(*c_iter, &_jacobian[0]);
+ jacobianDetSection->updatePoint(*c_iter, &_jacobianDet[0]);
+ basisDerivSection->updatePoint(*c_iter, &_basisDeriv[0]);
} // for
- _precomputed = true;
-} // precomputeGeometry
+} // computeGeometry
// ----------------------------------------------------------------------
+template<typename mesh_type>
void
-pylith::feassemble::Quadrature::retrieveGeometry(const ALE::Obj<SieveMesh>& mesh,
- const SieveMesh::point_type& cell)
+pylith::feassemble::Quadrature<mesh_type>::retrieveGeometry(const typename mesh_type::SieveMesh::point_type& cell)
{ // retrieveGeometry
- _quadPtsPreV->clear();
- mesh->restrictClosure(cell, *_quadPtsPreV);
+ typedef typename mesh_type::RealSection RealSection;
- _jacobianPreV->clear();
- mesh->restrictClosure(cell, *_jacobianPreV);
+ assert(0 != _quadPtsField);
+ assert(0 != _jacobianField);
+ assert(0 != _jacobianDetField);
+ assert(0 != _basisDerivField);
- _jacobianDetPreV->clear();
- mesh->restrictClosure(cell, *_jacobianDetPreV);
+ const ALE::Obj<RealSection>& quadPtsSection = _quadPtsField->section();
+ quadPtsSection->restrictPoint(cell, &_quadPts[0], _quadPts.size());
- _jacobianInvPreV->clear();
- mesh->restrictClosure(cell, *_jacobianInvPreV);
+ const ALE::Obj<RealSection>& jacobianSection = _jacobianField->section();
+ jacobianSection->restrictPoint(cell, &_jacobian[0], _jacobian.size());
- _basisDerivPreV->clear();
- mesh->restrictClosure(cell, *_basisDerivPreV);
+ const ALE::Obj<RealSection>& jacobianDetSection =
+ _jacobianDetField->section();
+ jacobianDetSection->restrictPoint(cell,
+ &_jacobianDet[0], _jacobianDet.size());
+
+ const ALE::Obj<RealSection>& basisDerivSection = _basisDerivField->section();
+ basisDerivSection->restrictPoint(cell, &_basisDeriv[0], _basisDeriv.size());
} // retrieveGeometry
+
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.hh 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.hh 2009-02-11 00:06:18 UTC (rev 14035)
@@ -33,12 +33,13 @@
#include "QuadratureBase.hh" // ISA QuadratureBase
#include "pylith/topology/topologyfwd.hh" // forward declarations
+#include "pylith/topology/topologyfwd.hh" // forward declarations
#include "pylith/utils/array.hh" // HASA double_array
// Quadrature -----------------------------------------------------------
template<typename mesh_type>
-class pylith::feassemble::Quadrature
+class pylith::feassemble::Quadrature : public QuadratureBase
{ // Quadrature
friend class TestQuadrature; // unit testing
@@ -93,25 +94,24 @@
*/
const double_array& jacobianDet(void) const;
- /// Reset the precomputation structures.
- void resetPrecomputation(void);
-
/** Precompute geometric quantities for each cell.
*
* @param mesh Finite-element mesh
* @param cells Finite-element cells for geometry.
*/
- void precomputeGeometry(const mesh_type& mesh,
- const ALE::Obj<mesh_type::SieveMesh::label_sequence>& cells);
+ void computeGeometry(const mesh_type& mesh,
+ const ALE::Obj<typename mesh_type::SieveMesh::label_sequence>& cells);
/** Retrieve precomputed geometric quantities for a cell.
*
* @param mesh Finite-element mesh
* @param cell Finite-element cell
*/
- void retrieveGeometry(const ALE::Obj<mesh_type::SieveMesh>& mesh,
- const mesh_type::SieveMesh::point_type& cell);
+ void retrieveGeometry(const typename mesh_type::SieveMesh::point_type& cell);
+ /// Deallocate temporary storage.
+ void clear(void);
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
@@ -121,14 +121,6 @@
*/
Quadrature(const Quadrature& q);
- /* Check determinant of Jacobian against minimum allowable value.
- *
- * @param det Value of determinant of Jacobian
- * @param cell Label of finite-element cell
- */
- void _checkJacobianDet(const double det,
- const mesh_type::SieveMesh::point_type& cell) const;
-
/// Set entries in geometry arrays to zero.
void _resetGeometry(void);
@@ -141,30 +133,28 @@
virtual
void _computeGeometry(const double* coordsVertices,
const int spaceDim,
- const mesh_type::SieveMesh::point_type& cell) = 0;
+ const int cell) = 0;
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
+ /** Buffers for cell data */
+ double_array _quadPts; ///< Coordinates of quad pts.
+ double_array _jacobian; ///< Jacobian at quad pts;
+ double_array _jacobianDet; ///< |J| at quad pts.
+ double_array _jacobianInv; /// Inverse of Jacobian at quad pts.
+ double_array _basisDeriv; ///< Deriv. of basis fns at quad pts.
+
/** Fields and visitors for precomputing geometry information for
* cells associated with this quadrature.
*/
- /** :OPTIMIZATION:
- * Need to check speed of retrieving geometry using Fields and Visitors.
- * Can switch to Sieve sections and visitors if too slow.
- */
- Field<mesh_type>* _quadPtsField; ///< Coordinates of quad pts.
- Visitor<Field<mesh_type> >* _quadPtsVisitor; ///< Visitor for quad pts.
+ topology::Field<mesh_type>* _quadPtsField; ///< Coordinates of quad pts.
+ topology::Field<mesh_type>* _jacobianField; ///< Jacobian at quad pts.
+ topology::Field<mesh_type>* _jacobianDetField; ///< |J| at quad pts.
- Field<mesh_type>* _jacobianField; ///< Jacobian at quad pts.
- Visitor<Field<mesh_type> >* _jacobianVisitor; ///< Visitor for Jacobian.
+ /// Derivatives of basis fns at quad pts.
+ topology::Field<mesh_type>* _basisDerivField;
- Field<mesh_type>* _jacobianDetField; ///< |J| at quad pts.
- Visitor<Field<mesh_type> >* _jacobianDetVisitor; ///< Visitor for |J|.
-
- Field<mesh_type>* _basisDerivField; ///< Deriv. of basis fns at quad pts.
- Visitor<Field<mesh_type> >* _basisDerivVisitor; ///< Visitor for derivatives.
-
bool _checkConditioning; ///< True if checking for ill-conditioning.
// NOT IMPLEMENTED //////////////////////////////////////////////////////
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.icc 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.icc 2009-02-11 00:06:18 UTC (rev 14035)
@@ -14,118 +14,54 @@
#error "Quadrature.icc must be included only from Quadrature.hh"
#else
-// Get minimum allowable Jacobian.
-inline
-double
-pylith::feassemble::Quadrature::minJacobian(void) const {
- return _minJacobian;
-}
-
-// Set minimum allowable Jacobian.
-inline
-void
-pylith::feassemble::Quadrature::minJacobian(const double min) {
- _minJacobian = min;
-}
-
// Set flag for checking ill-conditioning.
+template<typename mesh_type>
inline
void
-pylith::feassemble::Quadrature::checkConditioning(const bool flag) {
+pylith::feassemble::Quadrature<mesh_type>::checkConditioning(const bool flag) {
_checkConditioning = flag;
}
// Get flag for checking ill-conditioning.
+template<typename mesh_type>
inline
bool
-pylith::feassemble::Quadrature::checkConditioning(void) const {
+pylith::feassemble::Quadrature<mesh_type>::checkConditioning(void) const {
return _checkConditioning;
}
-// Get coordinates of quadrature points in reference cell.
-inline
-const pylith::double_array&
-pylith::feassemble::Quadrature::quadPtsRef(void) const {
- return _quadPtsRef;
-}
-
// Get coordinates of quadrature points in cell (NOT reference cell).
+template<typename mesh_type>
inline
const pylith::double_array&
-pylith::feassemble::Quadrature::quadPts(void) const {
+pylith::feassemble::Quadrature<mesh_type>::quadPts(void) const {
return _quadPts;
}
-// Get weights of quadrature points.
-inline
-const pylith::double_array&
-pylith::feassemble::Quadrature::quadWts(void) const {
- return _quadWts;
-}
-
-// Get basis fns evaluated at quadrature points.
-inline
-const pylith::double_array&
-pylith::feassemble::Quadrature::basis(void) const {
- return _basis;
-}
-
// Get derivatives of basis fns evaluated at quadrature points.
+template<typename mesh_type>
inline
const pylith::double_array&
-pylith::feassemble::Quadrature::basisDeriv(void) const {
+pylith::feassemble::Quadrature<mesh_type>::basisDeriv(void) const {
return _basisDeriv;
}
// Get Jacobians evaluated at quadrature points.
+template<typename mesh_type>
inline
const pylith::double_array&
-pylith::feassemble::Quadrature::jacobian(void) const {
+pylith::feassemble::Quadrature<mesh_type>::jacobian(void) const {
return _jacobian;
}
// Get determinants of Jacobian evaluated at quadrature points.
+template<typename mesh_type>
inline
const pylith::double_array&
-pylith::feassemble::Quadrature::jacobianDet(void) const {
+pylith::feassemble::Quadrature<mesh_type>::jacobianDet(void) const {
return _jacobianDet;
}
-// Get Jacobian inverses evaluated at quadrature points.
-inline
-const pylith::double_array&
-pylith::feassemble::Quadrature::jacobianInv(void) const {
- return _jacobianInv;
-}
-
-// Get number of dimensions in reference cell.
-inline
-int
-pylith::feassemble::Quadrature::cellDim(void) const {
- return _cellDim;
-}
-
-// Get number of basis functions for cell.
-inline
-int
-pylith::feassemble::Quadrature::numBasis(void) const {
- return _numBasis;
-}
-
-// Get number of quadrature points.
-inline
-int
-pylith::feassemble::Quadrature::numQuadPts(void) const {
- return _numQuadPts;
-}
-
-// Get number of dimensions in coordinates of cell vertices.
-inline
-int
-pylith::feassemble::Quadrature::spaceDim(void) const {
- return _spaceDim;
-}
-
#endif
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature1D.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature1D.cc 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature1D.cc 2009-02-11 00:06:18 UTC (rev 14035)
@@ -12,9 +12,8 @@
#include <portinfo>
-#include "Quadrature1D.hh" // implementation of class methods
+//#include "Quadrature1D.hh" // implementation of class methods
-#include "pylith/utils/array.hh" // USES double_array
#include "CellGeometry.hh" // USES CellGeometry
#include "petsc.h" // USES PetscLogFlops
@@ -25,30 +24,35 @@
// ----------------------------------------------------------------------
// Constructor
-pylith::feassemble::Quadrature1D::Quadrature1D(void) : Quadrature()
+template<typename mesh_type>
+pylith::feassemble::Quadrature1D<mesh_type>::Quadrature1D(void)
{ // constructor
} // constructor
// ----------------------------------------------------------------------
// Destructor
-pylith::feassemble::Quadrature1D::~Quadrature1D(void)
+template<typename mesh_type>
+pylith::feassemble::Quadrature1D<mesh_type>::~Quadrature1D(void)
{ // destructor
} // destructor
// ----------------------------------------------------------------------
// Copy constructor.
-pylith::feassemble::Quadrature1D::Quadrature1D(const Quadrature1D& q) :
- Quadrature(q)
+template<typename mesh_type>
+pylith::feassemble::Quadrature1D<mesh_type>::Quadrature1D(const Quadrature1D& q) :
+ Quadrature<mesh_type>(q)
{ // copy constructor
} // copy constructor
+#if 0
// ----------------------------------------------------------------------
// Compute geometric quantities for a cell at quadrature points.
+template<typename mesh_type>
void
-pylith::feassemble::Quadrature1D::computeGeometry(
- const real_section_type::value_type* vertCoords,
- const int coordDim,
- const Mesh::point_type& cell)
+pylith::feassemble::Quadrature1D<mesh_type>::computeGeometry(
+ const double* vertCoords,
+ const int coordDim,
+ const int cell)
{ // computeGeometry
assert(1 == _cellDim);
assert(1 == _spaceDim);
@@ -110,6 +114,6 @@
PetscLogFlops(_numQuadPts * (1 + _numBasis * 4));
} // computeGeometry
+#endif
-
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature1D.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature1D.hh 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature1D.hh 2009-02-11 00:06:18 UTC (rev 14035)
@@ -19,16 +19,10 @@
#if !defined(pylith_feassemble_quadrature1d_hh)
#define pylith_feassemble_quadrature1d_hh
-#include "Quadrature.hh"
+#include "Quadrature.hh" // ISA Quadrature
-namespace pylith {
- namespace feassemble {
- class Quadrature1D;
- class TestQuadrature1D;
- } // feassemble
-} // pylith
-
-class pylith::feassemble::Quadrature1D : public Quadrature
+template<typename mesh_type>
+class pylith::feassemble::Quadrature1D : public Quadrature<mesh_type>
{ // Quadrature1D
friend class TestQuadrature1D; // unit testing
@@ -42,7 +36,7 @@
~Quadrature1D(void);
/// Create a copy of this object.
- Quadrature* clone(void) const;
+ Quadrature<mesh_type>* clone(void) const;
/** Compute geometric quantities for a cell at quadrature points.
*
@@ -50,9 +44,9 @@
* @param coordinates Section containing vertex coordinates
* @param cell Finite-element cell
*/
- void computeGeometry(const real_section_type::value_type* vertCoords,
+ void computeGeometry(const double* vertCoords,
const int coordDim,
- const Mesh::point_type& cell);
+ const int cell);
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
@@ -71,6 +65,7 @@
}; // Quadrature1D
#include "Quadrature1D.icc" // inline methods
+#include "Quadrature1D.cc" // template methods
#endif // pylith_feassemble_quadrature1d_hh
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature1D.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature1D.icc 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature1D.icc 2009-02-11 00:06:18 UTC (rev 14035)
@@ -15,10 +15,11 @@
#else
// Create a copy of this object.
+template<typename mesh_type>
inline
-pylith::feassemble::Quadrature*
-pylith::feassemble::Quadrature1D::clone(void) const {
- return new Quadrature1D(*this);
+pylith::feassemble::Quadrature<mesh_type>*
+pylith::feassemble::Quadrature1D<mesh_type>::clone(void) const {
+ return new Quadrature1D<mesh_type>(*this);
} // clone
#endif
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature3D.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature3D.cc 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature3D.cc 2009-02-11 00:06:18 UTC (rev 14035)
@@ -12,7 +12,7 @@
#include <portinfo>
-#include "Quadrature3D.hh" // implementation of class methods
+//#include "Quadrature3D.hh" // implementation of class methods
#include "CellGeometry.hh" // USES CellGeometry
@@ -24,30 +24,34 @@
// ----------------------------------------------------------------------
// Constructor
-pylith::feassemble::Quadrature3D::Quadrature3D(void) : pylith::feassemble::Quadrature::Quadrature()
+template<typename mesh_type>
+pylith::feassemble::Quadrature3D<mesh_type>::Quadrature3D(void) : pylith::feassemble::Quadrature::Quadrature()
{ // constructor
} // constructor
// ----------------------------------------------------------------------
// Destructor
-pylith::feassemble::Quadrature3D::~Quadrature3D(void)
+template<typename mesh_type>
+pylith::feassemble::Quadrature3D<mesh_type>::~Quadrature3D(void)
{ // destructor
} // destructor
// ----------------------------------------------------------------------
// Copy constructor.
-pylith::feassemble::Quadrature3D::Quadrature3D(const Quadrature3D& q) :
+template<typename mesh_type>
+pylith::feassemble::Quadrature3D<mesh_type>::Quadrature3D(const Quadrature3D& q) :
Quadrature(q)
{ // copy constructor
} // copy constructor
// ----------------------------------------------------------------------
// Compute geometric quantities for a cell at quadrature points.
+template<typename mesh_type>
void
-pylith::feassemble::Quadrature3D::computeGeometry(
- const real_section_type::value_type* vertCoords,
- const int coordDim,
- const Mesh::point_type& cell)
+pylith::feassemble::Quadrature3D<mesh_type>::computeGeometry(
+ const double* vertCoords,
+ const int coordDim,
+ const int cell)
{ // computeGeometry
assert(3 == _cellDim);
assert(3 == _spaceDim);
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature3D.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature3D.hh 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature3D.hh 2009-02-11 00:06:18 UTC (rev 14035)
@@ -19,16 +19,10 @@
#if !defined(pylith_feassemble_quadrature3d_hh)
#define pylith_feassemble_quadrature3d_hh
-#include "Quadrature.hh"
+#include "Quadrature.hh" // ISA Quadrature
-namespace pylith {
- namespace feassemble {
- class Quadrature3D;
- class TestQuadrature3D;
- } // feassemble
-} // pylith
-
-class pylith::feassemble::Quadrature3D : public Quadrature
+template<typename mesh_type>
+class pylith::feassemble::Quadrature3D : public Quadrature<mesh_type>
{ // Quadrature3D
friend class TestQuadrature3D; // unit testing
@@ -50,9 +44,9 @@
* @param coordinates Section containing vertex coordinates
* @param cell Finite-element cell
*/
- void computeGeometry(const real_section_type::value_type* vertCoords,
+ void computeGeometry(const double* vertCoords,
const int coordDim,
- const Mesh::point_type& cell);
+ const int cell);
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
@@ -72,6 +66,7 @@
}; // Quadrature3D
#include "Quadrature3D.icc" // inline methods
+#include "Quadrature3D.cc" // template methods
#endif // pylith_feassemble_quadrature3d_hh
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.cc 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.cc 2009-02-11 00:06:18 UTC (rev 14035)
@@ -181,5 +181,20 @@
return *_geometry;
} // refGeometry
+// ----------------------------------------------------------------------
+// Check determinant of Jacobian against minimum allowable value
+void
+pylith::feassemble::QuadratureBase::_checkJacobianDet(const double det,
+ const int cell) const
+{ // _checkJacobianDet
+ if (det < _minJacobian) {
+ std::ostringstream msg;
+ msg << "Determinant of Jacobian (" << det << ") for cell " << cell
+ << " is smaller than minimum permissible value (" << _minJacobian
+ << ")!\n";
+ throw std::runtime_error(msg.str());
+ } // if
+} // _checkJacobianDet
+
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.hh 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.hh 2009-02-11 00:06:18 UTC (rev 14035)
@@ -164,6 +164,14 @@
*/
QuadratureBase(const QuadratureBase& q);
+ /* Check determinant of Jacobian against minimum allowable value.
+ *
+ * @param det Value of determinant of Jacobian
+ * @param cell Label of finite-element cell
+ */
+ void _checkJacobianDet(const double det,
+ const int cell) const;
+
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Makefile.am 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Makefile.am 2009-02-11 00:06:18 UTC (rev 14035)
@@ -28,6 +28,8 @@
SubMesh.hh \
SubMesh.icc \
SubMesh.cc \
+ RestrictVisitor.hh \
+ RestrictVisitor.icc \
topologyfwd.hh
noinst_HEADERS =
Added: short/3D/PyLith/branches/pylith-swig/libsrc/topology/RestrictVisitor.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/RestrictVisitor.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/RestrictVisitor.hh 2009-02-11 00:06:18 UTC (rev 14035)
@@ -0,0 +1,98 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/topology/Visitor.hh
+ *
+ * @brief Visitor for managing values extracted from restriction of a
+ * field to a cell.
+ *
+ * Wraps a Sieve RestrictVisitor.
+ */
+
+#if !defined(pylith_topology_restrictvisitor_hh)
+#define pylith_topology_restrictvisitor_hh
+
+// Include directives ---------------------------------------------------
+#include "topologyfwd.hh" // forward declarations
+
+// RestrictVisitor ------------------------------------------------------
+template<typename field_type>
+class pylith::topology::RestrictVisitor
+{ // Field
+ friend class TestRestrictVisitor; // unit testing
+
+// PRIVATE TYPEDEFS /////////////////////////////////////////////////////
+private:
+
+ // Convenience typedefs
+ typedef typename field_type::RealSection RealSection;
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /** Constructor with field and size.
+ *
+ * @param field Field over finite-element mesh.
+ * @param size Fiber dimension for field.
+ */
+ RestrictVisitor(const field_type& field,
+ const int size);
+
+ /** Constructor with field and array.
+ *
+ * @param field Field over finite-element mesh.
+ * @param values Array of values to use for storage.
+ */
+ RestrictVisitor(const field_type& field,
+ const double_array& values);
+
+ /// Destructor.
+ ~RestrictVisitor(void);
+
+ /** Visit field for cell and retrieve values.
+ *
+ * @param cell Cell in finite-element mesh.
+ */
+ void visit(const typename field_type::Mesh::SieveMesh::point_type& cell);
+
+ /** Get field values previously retrieved for cell.
+ *
+ * @returns Values for field.
+ */
+ const double* values(void) const;
+
+ /// Clear values associated with cell.
+ void clear(void);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ const field_type& _field; ///< Field associated with visitor
+
+ /// Sieve visitor that manages storage.
+ ALE::ISieveVisitor::RestrictVisitor<RealSection> _visitor;
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ RestrictVisitor(const RestrictVisitor&); ///< Not implemented
+ const RestrictVisitor& operator=(const RestrictVisitor&); ///< Not implemented
+
+}; // Field
+
+#include "RestrictVisitor.icc"
+
+#endif // pylith_topology_restrictvisitor_hh
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/topology/RestrictVisitor.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/RestrictVisitor.icc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/RestrictVisitor.icc 2009-02-11 00:06:18 UTC (rev 14035)
@@ -0,0 +1,73 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_topology_restrictvisitor_hh)
+#error "RestrictVisitor.icc must be included only from RestrictVisitor.hh"
+#else
+
+// Constructor with field and size.
+template<typename field_type>
+pylith::topology::RestrictVisitor<field_type>::RestrictVisitor(const field_type& field,
+ const int size) :
+ _field(field),
+ _visitor(field.section(), size)
+{ // constructor
+} // constructor
+
+// Constructor with field and array.
+template<typename field_type>
+pylith::topology::RestrictVisitor<field_type>::RestrictVisitor(const field_type& field,
+ const double_array& values) :
+ _field(field),
+ _visitor(field.section(), &values[0], values.size()
+{ // constructor
+} // constructor
+
+// Destructor.
+template<typename field_type>
+~RestrictVisitor(void)
+{ // destructor
+} // destructor
+
+// Visit field for cell and retrieve values.
+template<typename field_type>
+inline
+void
+pylith::topology::RestrictVisitor<field_type>::visit(
+ const typename field_type::Mesh::SieveMesh::point_type& cell)
+{ // visit
+ const ALE::Obj<SieveMesh>& sieveMesh = field.mesh().sieveMesh();
+ mesh->restrictClosure(cell, _visitor);
+} // visit
+
+// Get field values previously retrieved for cell.
+template<typename field_type>
+inline
+const double*
+pylith::topology::RestrictVisitor<field_type>::values(void) const
+{ // values
+ return _visitor->getValues();
+} // values
+
+// Clear values associated with cell.
+template<typename field_type>
+void
+pylith::topology::RestrictVisitor<field_type>::clear(void)
+{ // clear
+ return _visitor->clear();
+} // clear
+
+
+#endif
+
+
+// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/topologyfwd.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/topologyfwd.hh 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/topologyfwd.hh 2009-02-11 00:06:18 UTC (rev 14035)
@@ -31,6 +31,7 @@
typedef Fields<Field<Mesh> > FieldsMesh;
typedef Fields<Field<SubMesh> > FieldsSubMesh;
class SolutionFields;
+ template<typename field_type> class RestrictVisitor;
class MeshOps;
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/Makefile.am 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/Makefile.am 2009-02-11 00:06:18 UTC (rev 14035)
@@ -34,17 +34,18 @@
TestGeometryQuad3D.cc \
TestGeometryTet3D.cc \
TestGeometryHex3D.cc \
+ TestQuadratureBase.cc \
+ TestQuadrature.cc \
+ TestQuadrature3D.cc \
test_feassemble.cc
# TestIntegrator.cc \
# TestIntegratorElasticity.cc \
-# TestQuadrature.cc \
# TestQuadrature0D.cc \
# TestQuadrature1D.cc \
# TestQuadrature1Din2D.cc \
# TestQuadrature2D.cc \
# TestQuadrature2Din3D.cc \
-# TestQuadrature3D.cc \
# TestElasticityExplicit.cc \
# TestElasticityExplicit1DLinear.cc \
# TestElasticityExplicit1DQuadratic.cc \
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadrature.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadrature.cc 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadrature.cc 2009-02-11 00:06:18 UTC (rev 14035)
@@ -12,6 +12,7 @@
#include <portinfo>
+#include "pylith/topology/Mesh.hh" // USES Mesh
#include "TestQuadrature.hh" // Implementation of class methods
#include "pylith/feassemble/Quadrature1D.hh" // USES Quadrature1D
@@ -47,7 +48,7 @@
GeometryLine1D geometry;
// Set values
- Quadrature1D qOrig;
+ Quadrature1D<topology::Mesh> qOrig;
qOrig._minJacobian = minJacobianE;
qOrig._checkConditioning = checkConditioning;
qOrig._cellDim = cellDimE;
@@ -90,7 +91,7 @@
qOrig._geometry = geometry.clone();
// Clone
- const Quadrature* qCopy = qOrig.clone();
+ const Quadrature<topology::Mesh>* qCopy = qOrig.clone();
// Check clone
CPPUNIT_ASSERT(0 != qCopy);
@@ -156,25 +157,14 @@
CPPUNIT_ASSERT_EQUAL(geometry.numCorners(), qCopy->_geometry->numCorners());
delete qCopy; qCopy = 0;
-} // testCopy
+} // testClone
// ----------------------------------------------------------------------
-// Test minJacobian()
-void
-pylith::feassemble::TestQuadrature::testMinJacobian(void)
-{ // testMinJacobian
- Quadrature1D q;
- const double min = 1.0;
- q.minJacobian(min);
- CPPUNIT_ASSERT_EQUAL(min, q._minJacobian);
-} // testMinJacobian
-
-// ----------------------------------------------------------------------
// Test checkConditioning()
void
pylith::feassemble::TestQuadrature::testCheckConditioning(void)
{ // testCheckConditioning
- Quadrature1D q;
+ Quadrature1D<topology::Mesh> q;
CPPUNIT_ASSERT_EQUAL(false, q.checkConditioning());
q.checkConditioning(true);
CPPUNIT_ASSERT_EQUAL(true, q.checkConditioning());
@@ -182,80 +172,11 @@
CPPUNIT_ASSERT_EQUAL(false, q.checkConditioning());
} // testCheckConditioning
+#if 0
// ----------------------------------------------------------------------
-// Test refGeometry()
+// Test computeGeometry() and retrieveGeometry() for meshes.
void
-pylith::feassemble::TestQuadrature::testRefGeometry(void)
-{ // testRefGeometry
- GeometryLine1D geometry;
- Quadrature1D quadrature;
- quadrature.refGeometry(&geometry);
- const CellGeometry& test = quadrature.refGeometry();
-
- CPPUNIT_ASSERT_EQUAL(geometry.cellDim(), test.cellDim());
- CPPUNIT_ASSERT_EQUAL(geometry.spaceDim(), test.spaceDim());
- CPPUNIT_ASSERT_EQUAL(geometry.numCorners(), test.numCorners());
-} // testRefGeometry
-
-// ----------------------------------------------------------------------
-// Test initialize()
-void
-pylith::feassemble::TestQuadrature::testInitialize(void)
-{ // initialize
-
- const int cellDim = 1;
- const int numBasis = 2;
- const int numQuadPts = 1;
- const int spaceDim = 1;
- const double basis[] = { 0.5, 0.5 };
- const double basisDerivRef[] = { -0.5, 0.5 };
- const double quadPtsRef[] = { 0.0 };
- const double quadWts[] = { 2.0 };
- const double minJacobian = 1.0;
-
- Quadrature1D q;
- q.initialize(basis, basisDerivRef, quadPtsRef, quadWts,
- cellDim, numBasis, numQuadPts, spaceDim);
-
- CPPUNIT_ASSERT_EQUAL(cellDim, q._cellDim);
- CPPUNIT_ASSERT_EQUAL(numBasis, q._numBasis);
- CPPUNIT_ASSERT_EQUAL(numQuadPts, q._numQuadPts);
- CPPUNIT_ASSERT_EQUAL(spaceDim, q._spaceDim);
-
- size_t size = numBasis * numQuadPts;
- for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_EQUAL(basis[i], q._basis[i]);
-
- size = numBasis * numQuadPts * spaceDim;
- for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_EQUAL(basisDerivRef[i], q._basisDerivRef[i]);
-
- size = numQuadPts * cellDim;
- for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_EQUAL(quadPtsRef[i], q._quadPtsRef[i]);
-
- size = numQuadPts;
- for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_EQUAL(quadWts[i], q._quadWts[i]);
-
- // Make sure Jacobian stuff has been allocated
- size = numQuadPts*cellDim*spaceDim;
- CPPUNIT_ASSERT_EQUAL(size, q._jacobian.size());
-
- size = numQuadPts*spaceDim*cellDim;
- CPPUNIT_ASSERT_EQUAL(size, q._jacobianInv.size());
-
- size = numQuadPts;
- CPPUNIT_ASSERT_EQUAL(size, q._jacobianDet.size());
-
- size = numQuadPts*spaceDim;
- CPPUNIT_ASSERT_EQUAL(size, q._quadPts.size());
-} // initialize
-
-// ----------------------------------------------------------------------
-// Test initialize() & computeGeometry()
-void
-pylith::feassemble::TestQuadrature::_testComputeGeometry(Quadrature* pQuad,
+pylith::feassemble::TestQuadrature::_testComputeGeometry(Quadrature<topology::Mesh>* pQuad,
const QuadratureData& data) const
{ // testComputeGeometry
const int cellDim = data.cellDim;
@@ -337,6 +258,6 @@
tolerance);
} // testComputeGeometry
+#endif
-
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadrature.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadrature.hh 2009-02-10 23:21:26 UTC (rev 14034)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadrature.hh 2009-02-11 00:06:18 UTC (rev 14035)
@@ -23,10 +23,12 @@
#include <cppunit/extensions/HelperMacros.h>
+#include "pylith/feassemble/feassemblefwd.hh" // forward declarations
+#include "pylith/topology/topologyfwd.hh" // forward declarations
+
/// Namespace for pylith package
namespace pylith {
namespace feassemble {
- class Quadrature;
class TestQuadrature;
class QuadratureData;
} // feassemble
@@ -40,10 +42,7 @@
CPPUNIT_TEST_SUITE( TestQuadrature );
CPPUNIT_TEST( testClone );
- CPPUNIT_TEST( testMinJacobian );
CPPUNIT_TEST( testCheckConditioning );
- CPPUNIT_TEST( testRefGeometry );
- CPPUNIT_TEST( testInitialize );
CPPUNIT_TEST_SUITE_END();
@@ -53,27 +52,18 @@
/// Test clone()
void testClone(void);
- /// Test minJacobian()
- void testMinJacobian(void);
-
void testCheckConditioning(void);
/// Test checkConditioning()
- /// Test refGeometry()
- void testRefGeometry(void);
-
- /// Test initialize()
- void testInitialize(void);
-
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
- /** Test initialize() & computeGeometry()
+ /** Test computeGeometry() and retrieveGeometry().
*
* @param pQuad Pointer to quadrature
* @param data Data for testing quadrature
*/
- void _testComputeGeometry(Quadrature* pQuad,
+ void _testComputeGeometry(Quadrature<topology::Mesh>* pQuad,
const QuadratureData& data) const;
}; // class TestQuadrature
Added: short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadratureBase.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadratureBase.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadratureBase.cc 2009-02-11 00:06:18 UTC (rev 14035)
@@ -0,0 +1,106 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestQuadratureBase.hh" // Implementation of class methods
+
+#include "pylith/feassemble/QuadratureBase.hh" // USES QuadratureBase
+#include "pylith/feassemble/GeometryLine1D.hh" // USES GeometryLine1D
+
+#include "data/QuadratureData.hh" // USES QuadratureData
+
+#include <string.h> // USES memcpy()
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::feassemble::TestQuadratureBase );
+
+// ----------------------------------------------------------------------
+// Test constructor.
+void
+pylith::feassemble::TestQuadratureBase::testConstructor(void)
+{ // testConstructor
+ QuadratureBase q;
+} // testMinJacobian
+
+// ----------------------------------------------------------------------
+// Test minJacobian()
+void
+pylith::feassemble::TestQuadratureBase::testMinJacobian(void)
+{ // testMinJacobian
+ QuadratureBase q;
+
+ const double min = 1.0;
+ q.minJacobian(min);
+ CPPUNIT_ASSERT_EQUAL(min, q._minJacobian);
+} // testMinJacobian
+
+// ----------------------------------------------------------------------
+// Test refGeometry()
+void
+pylith::feassemble::TestQuadratureBase::testRefGeometry(void)
+{ // testRefGeometry
+ GeometryLine1D geometry;
+ QuadratureBase quadrature;
+
+ quadrature.refGeometry(&geometry);
+ const CellGeometry& test = quadrature.refGeometry();
+
+ CPPUNIT_ASSERT_EQUAL(geometry.cellDim(), test.cellDim());
+ CPPUNIT_ASSERT_EQUAL(geometry.spaceDim(), test.spaceDim());
+ CPPUNIT_ASSERT_EQUAL(geometry.numCorners(), test.numCorners());
+} // testRefGeometry
+
+// ----------------------------------------------------------------------
+// Test initialize()
+void
+pylith::feassemble::TestQuadratureBase::testInitialize(void)
+{ // initialize
+
+ const int cellDim = 1;
+ const int numBasis = 2;
+ const int numQuadPts = 1;
+ const int spaceDim = 1;
+ const double basis[] = { 0.5, 0.5 };
+ const double basisDerivRef[] = { -0.5, 0.5 };
+ const double quadPtsRef[] = { 0.0 };
+ const double quadWts[] = { 2.0 };
+ const double minJacobian = 1.0;
+
+ QuadratureBase q;
+ q.initialize(basis, basisDerivRef, quadPtsRef, quadWts,
+ cellDim, numBasis, numQuadPts, spaceDim);
+
+ CPPUNIT_ASSERT_EQUAL(cellDim, q._cellDim);
+ CPPUNIT_ASSERT_EQUAL(numBasis, q._numBasis);
+ CPPUNIT_ASSERT_EQUAL(numQuadPts, q._numQuadPts);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, q._spaceDim);
+
+ size_t size = numBasis * numQuadPts;
+ for (int i=0; i < size; ++i)
+ CPPUNIT_ASSERT_EQUAL(basis[i], q._basis[i]);
+
+ size = numBasis * numQuadPts * spaceDim;
+ for (int i=0; i < size; ++i)
+ CPPUNIT_ASSERT_EQUAL(basisDerivRef[i], q._basisDerivRef[i]);
+
+ size = numQuadPts * cellDim;
+ for (int i=0; i < size; ++i)
+ CPPUNIT_ASSERT_EQUAL(quadPtsRef[i], q._quadPtsRef[i]);
+
+ size = numQuadPts;
+ for (int i=0; i < size; ++i)
+ CPPUNIT_ASSERT_EQUAL(quadWts[i], q._quadWts[i]);
+} // initialize
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadratureBase.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadratureBase.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadratureBase.hh 2009-02-11 00:06:18 UTC (rev 14035)
@@ -0,0 +1,67 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/feassemble/TestQuadratureBase.hh
+ *
+ * @brief C++ TestQuadratureBase object
+ *
+ * C++ unit testing for QuadratureBase.
+ */
+
+#if !defined(pylith_feassemble_testquadraturebase_hh)
+#define pylith_feassemble_testquadraturebase_hh
+
+#include <cppunit/extensions/HelperMacros.h>
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace feassemble {
+ class TestQuadratureBase;
+ class QuadratureBaseData;
+ } // feassemble
+} // pylith
+
+/// C++ unit testing for QuadratureBase
+class pylith::feassemble::TestQuadratureBase : public CppUnit::TestFixture
+{ // class TestQuadratureBase
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestQuadratureBase );
+
+ CPPUNIT_TEST( testConstructor );
+ CPPUNIT_TEST( testMinJacobian );
+ CPPUNIT_TEST( testRefGeometry );
+ CPPUNIT_TEST( testInitialize );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Test constructor.
+ void testConstructor(void);
+
+ /// Test minJacobian()
+ void testMinJacobian(void);
+
+ /// Test refGeometry()
+ void testRefGeometry(void);
+
+ /// Test initialize()
+ void testInitialize(void);
+
+}; // class TestQuadratureBase
+
+#endif // pylith_feassemble_testquadraturebase_hh
+
+// End of file
More information about the CIG-COMMITS
mailing list