[cig-commits] r14049 - in short/3D/PyLith/branches/pylith-swig: libsrc/feassemble libsrc/topology unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Fri Feb 13 15:05:03 PST 2009
Author: brad
Date: 2009-02-13 15:05:03 -0800 (Fri, 13 Feb 2009)
New Revision: 14049
Modified:
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/topology/Mesh.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
short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadratureEngine.cc
Log:
Updated testing of basic quadrature object.
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.cc 2009-02-13 19:52:03 UTC (rev 14048)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.cc 2009-02-13 23:05:03 UTC (rev 14049)
@@ -12,11 +12,9 @@
#include <portinfo>
-//#include "Quadrature.hh" // implementation of class methods
-
#include "CellGeometry.hh" // USES CellGeometry
-#if 0
+#include "QuadratureEngine.hh" // USES QuadratureEngine
#include "Quadrature0D.hh"
#include "Quadrature1D.hh"
#include "Quadrature1Din2D.hh"
@@ -24,7 +22,6 @@
#include "Quadrature2D.hh"
#include "Quadrature2Din3D.hh"
#include "Quadrature3D.hh"
-#endif
#include "pylith/topology/Field.hh" // HOLDSA Field
@@ -62,7 +59,7 @@
// Copy constructor
template<typename mesh_type>
pylith::feassemble::Quadrature<mesh_type>::Quadrature(const Quadrature& q) :
- QuadratureBase(q),
+ QuadratureRefCell(q),
_engine(0),
_quadPtsField(0),
_jacobianField(0),
@@ -82,14 +79,16 @@
const ALE::Obj<typename mesh_type::SieveMesh::label_sequence>& cells)
{ // precomputeGeometry
typedef typename mesh_type::RealSection RealSection;
+ typedef typename mesh_type::SieveMesh::label_sequence label_sequence;
+ typedef typename mesh_type::RestrictVisitor RestrictVisitor;
- _setupEngine();
-
const char* loggingStage = "QuadratureCreation";
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush(loggingStage);
clear();
+ _setupEngine();
+ assert(0 != _engine);
// Allocate field and cell buffer for quadrature points
int fiberDim = _numQuadPts * _spaceDim;
@@ -140,9 +139,6 @@
<< std::endl;
#endif
- typedef typename mesh_type::SieveMesh::label_sequence label_sequence;
- typedef ALE::ISieveVisitor::RestrictVisitor<RealSection> RealSectionVisitor;
-
const typename label_sequence::iterator cellsEnd = cells->end();
assert(0 != _geometry);
const int numCorners = _geometry->numCorners();
@@ -150,7 +146,7 @@
assert(!sieveMesh.isNull());
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
- RealSectionVisitor coordsVisitor(coordinates, numCorners*_spaceDim);
+ RestrictVisitor coordsVisitor(*coordinates, numCorners*_spaceDim);
const ALE::Obj<RealSection>& quadPtsSection = _quadPtsField->section();
const ALE::Obj<RealSection>& jacobianSection = _jacobianField->section();
@@ -158,20 +154,24 @@
_jacobianDetField->section();
const ALE::Obj<RealSection>& basisDerivSection = _basisDerivField->section();
+ const double_array& quadPts = _engine->quadPts();
+ const double_array& jacobian = _engine->jacobian();
+ const double_array& jacobianDet = _engine->jacobianDet();
+ const double_array& basisDeriv = _engine->basisDeriv();
+
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);
- _resetGeometry();
- computeGeometry(cellVertexCoords, _spaceDim, *c_iter);
+ _engine->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]);
+ 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
} // computeGeometry
@@ -186,20 +186,29 @@
assert(0 != _jacobianField);
assert(0 != _jacobianDetField);
assert(0 != _basisDerivField);
+ assert(0 != _engine);
+ const double_array& quadPts = _engine->quadPts();
+ const double_array& jacobian = _engine->jacobian();
+ const double_array& jacobianDet = _engine->jacobianDet();
+ const double_array& basisDeriv = _engine->basisDeriv();
+
const ALE::Obj<RealSection>& quadPtsSection = _quadPtsField->section();
- quadPtsSection->restrictPoint(cell, &_quadPts[0], _quadPts.size());
+ quadPtsSection->restrictPoint(cell, const_cast<double*>(&quadPts[0]),
+ quadPts.size());
const ALE::Obj<RealSection>& jacobianSection = _jacobianField->section();
- jacobianSection->restrictPoint(cell, &_jacobian[0], _jacobian.size());
+ jacobianSection->restrictPoint(cell, const_cast<double*>(&jacobian[0]),
+ jacobian.size());
const ALE::Obj<RealSection>& jacobianDetSection =
_jacobianDetField->section();
- jacobianDetSection->restrictPoint(cell,
- &_jacobianDet[0], _jacobianDet.size());
+ jacobianDetSection->restrictPoint(cell, const_cast<double*>(&jacobianDet[0]),
+ jacobianDet.size());
const ALE::Obj<RealSection>& basisDerivSection = _basisDerivField->section();
- basisDerivSection->restrictPoint(cell, &_basisDeriv[0], _basisDeriv.size());
+ basisDerivSection->restrictPoint(cell, const_cast<double*>(&basisDeriv[0]),
+ basisDeriv.size());
} // retrieveGeometry
// ----------------------------------------------------------------------
@@ -228,7 +237,6 @@
const int cellDim = _cellDim;
const int spaceDim = _spaceDim;
-#if 0
if (1 == spaceDim)
if (1 == cellDim)
_engine = new Quadrature1D(*this);
@@ -240,7 +248,7 @@
<< std::endl;
assert(0);
} // if/else
- else if (2 == spaceDim) {
+ else if (2 == spaceDim)
if (2 == cellDim)
_engine = new Quadrature2D(*this);
else if (1 == cellDim)
@@ -253,7 +261,7 @@
<< std::endl;
assert(0);
} // if/else
- else if (3 == spaceDim) {
+ else if (3 == spaceDim)
if (3 == cellDim)
_engine = new Quadrature3D(*this);
else if (2 == cellDim)
@@ -268,7 +276,15 @@
<< std::endl;
assert(0);
} // if/else
-#endif
+ else {
+ std::cerr << "Unknown quadrature case with cellDim '"
+ << cellDim << "' and spaceDim '" << spaceDim << "'"
+ << std::endl;
+ assert(0);
+ } // if/else
+
+ assert(0 != _engine);
+ _engine->initialize();
} // _setupEngine
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.hh 2009-02-13 19:52:03 UTC (rev 14048)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.hh 2009-02-13 23:05:03 UTC (rev 14049)
@@ -18,7 +18,7 @@
*
* This object contains the informatio needed to perform numerical
* quadrature over a finite-element cell. It inherits quadrature
- * information over the reference cell from the QuadratureBase object.
+ * information over the reference cell from the QuadratureRefCell object.
* Given a cell this object will compute the cell's Jacobian, the
* determinant of the Jacobian, the inverse of the Jacobian, and the
@@ -121,6 +121,8 @@
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
+ QuadratureEngine* _engine; ///< Quadrature geometry engine.
+
/** Fields and visitors for precomputing geometry information for
* cells associated with this quadrature.
*/
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.icc 2009-02-13 19:52:03 UTC (rev 14048)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.icc 2009-02-13 23:05:03 UTC (rev 14049)
@@ -15,6 +15,7 @@
#else
#include <cassert> // USES assert()
+#include "QuadratureEngine.hh" // USES QuadratureEngine
// Set flag for checking ill-conditioning.
template<typename mesh_type>
@@ -38,7 +39,7 @@
const pylith::double_array&
pylith::feassemble::Quadrature<mesh_type>::quadPts(void) const {
assert(0 != _engine);
- return _engine->_quadPts;
+ return _engine->quadPts();
}
// Get derivatives of basis fns evaluated at quadrature points.
@@ -47,7 +48,7 @@
const pylith::double_array&
pylith::feassemble::Quadrature<mesh_type>::basisDeriv(void) const {
assert(0 != _engine);
- return _engine->_basisDeriv;
+ return _engine->basisDeriv();
}
// Get Jacobians evaluated at quadrature points.
@@ -56,7 +57,7 @@
const pylith::double_array&
pylith::feassemble::Quadrature<mesh_type>::jacobian(void) const {
assert(0 != _engine);
- return _engine->_jacobian;
+ return _engine->jacobian();
}
// Get determinants of Jacobian evaluated at quadrature points.
@@ -65,7 +66,7 @@
const pylith::double_array&
pylith::feassemble::Quadrature<mesh_type>::jacobianDet(void) const {
assert(0 != _engine);
- return _engine->_jacobianDet;
+ return _engine->jacobianDet();
}
#endif
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.hh 2009-02-13 19:52:03 UTC (rev 14048)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.hh 2009-02-13 23:05:03 UTC (rev 14049)
@@ -49,7 +49,9 @@
typedef SieveMesh::real_section_type RealSection;
typedef SieveMesh::int_section_type IntSection;
typedef ALE::IMesh<ALE::LabelSifter<int, SieveMesh::point_type> > SieveSubMesh;
+ typedef ALE::ISieveVisitor::RestrictVisitor<RealSection> RestrictVisitor;
+
// PUBLIC METHODS ///////////////////////////////////////////////////////
public :
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-13 19:52:03 UTC (rev 14048)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/Makefile.am 2009-02-13 23:05:03 UTC (rev 14049)
@@ -43,6 +43,7 @@
TestQuadrature2D.cc \
TestQuadrature2Din3D.cc \
TestQuadrature3D.cc \
+ TestQuadrature.cc \
test_feassemble.cc
# TestIntegrator.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-13 19:52:03 UTC (rev 14048)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadrature.cc 2009-02-13 23:05:03 UTC (rev 14049)
@@ -12,13 +12,15 @@
#include <portinfo>
-#include "pylith/topology/Mesh.hh" // USES Mesh
#include "TestQuadrature.hh" // Implementation of class methods
-#include "pylith/feassemble/Quadrature1D.hh" // USES Quadrature1D
+#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+
#include "pylith/feassemble/GeometryLine1D.hh" // USES GeometryLine1D
+#include "pylith/feassemble/GeometryTri2D.hh" // USES GeometryTri2D
-#include "data/QuadratureData.hh" // USES QuadratureData
+#include "data/QuadratureData2DLinear.hh" // USES QuadratureData2DLinear
#include <string.h> // USES memcpy()
@@ -26,9 +28,9 @@
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::feassemble::TestQuadrature );
// ----------------------------------------------------------------------
-// Test clone
+// Test copy constuctor.
void
-pylith::feassemble::TestQuadrature::testClone(void)
+pylith::feassemble::TestQuadrature::testCopyConstructor(void)
{ // testClone
// Semi-random values manually set to check cloning
const double minJacobianE = 1.0;
@@ -48,123 +50,61 @@
GeometryLine1D geometry;
// Set values
- Quadrature1D<topology::Mesh> qOrig;
- qOrig._minJacobian = minJacobianE;
- qOrig._checkConditioning = checkConditioning;
- qOrig._cellDim = cellDimE;
- qOrig._numBasis = numBasisE;
- qOrig._numQuadPts = numQuadPtsE;
- qOrig._spaceDim = spaceDimE;
+ Quadrature<topology::Mesh> qOrig;
+ qOrig.refGeometry(&geometry);
+ qOrig.minJacobian(minJacobianE);
+ qOrig.checkConditioning(checkConditioning);
+ qOrig.initialize(basisE, basisDerivE, quadPtsRefE, quadWtsE,
+ cellDimE, numBasisE, numQuadPtsE, spaceDimE);
- size_t size = 2;
- qOrig._basis.resize(size);
- memcpy(&qOrig._basis[0], basisE, size*sizeof(double));
+ // Copy
+ Quadrature<topology::Mesh> qCopy(qOrig);
- size = 2;
- qOrig._basisDerivRef.resize(size);
- memcpy(&qOrig._basisDerivRef[0], basisDerivE, size*sizeof(double));
+ // Check copy
+ CPPUNIT_ASSERT(0 == qCopy._engine);
+ CPPUNIT_ASSERT_EQUAL(minJacobianE, qCopy._minJacobian);
+ CPPUNIT_ASSERT_EQUAL(checkConditioning, qCopy._checkConditioning);
+ CPPUNIT_ASSERT_EQUAL(cellDimE, qCopy.cellDim());
+ CPPUNIT_ASSERT_EQUAL(numBasisE, qCopy.numBasis());
+ CPPUNIT_ASSERT_EQUAL(numQuadPtsE, qCopy.numQuadPts());
+ CPPUNIT_ASSERT_EQUAL(spaceDimE, qCopy.spaceDim());
- size = 1;
- qOrig._quadPtsRef.resize(size);
- memcpy(&qOrig._quadPtsRef[0], quadPtsRefE, size*sizeof(double));
-
- size = 1;
- qOrig._quadWts.resize(size);
- memcpy(&qOrig._quadWts[0], quadWtsE, size*sizeof(double));
-
- size = 1;
- qOrig._quadPts.resize(size);
- memcpy(&qOrig._quadPts[0], quadPtsE, size*sizeof(double));
-
- size = 1;
- qOrig._jacobian.resize(size);
- memcpy(&qOrig._jacobian[0], jacobianE, size*sizeof(double));
-
- size = 1;
- qOrig._jacobianInv.resize(size);
- memcpy(&qOrig._jacobianInv[0], jacobianInvE, size*sizeof(double));
-
- size = 1;
- qOrig._jacobianDet.resize(size);
- memcpy(&qOrig._jacobianDet[0], jacobianDetE, size*sizeof(double));
-
- qOrig._geometry = geometry.clone();
-
- // Clone
- const Quadrature<topology::Mesh>* qCopy = qOrig.clone();
-
- // Check clone
- CPPUNIT_ASSERT(0 != qCopy);
-
- CPPUNIT_ASSERT_EQUAL(minJacobianE, qCopy->_minJacobian);
- CPPUNIT_ASSERT_EQUAL(checkConditioning, qCopy->_checkConditioning);
- CPPUNIT_ASSERT_EQUAL(cellDimE, qCopy->cellDim());
- CPPUNIT_ASSERT_EQUAL(numBasisE, qCopy->numBasis());
- CPPUNIT_ASSERT_EQUAL(numQuadPtsE, qCopy->numQuadPts());
- CPPUNIT_ASSERT_EQUAL(spaceDimE, qCopy->spaceDim());
-
- const double_array& basis = qCopy->basis();
- size = numBasisE * numQuadPtsE;
+ const double_array& basis = qCopy.basis();
+ size_t size = numBasisE * numQuadPtsE;
CPPUNIT_ASSERT_EQUAL(size, basis.size());
for (int i=0; i < size; ++i)
CPPUNIT_ASSERT_EQUAL(basisE[i], basis[i]);
- const double_array& basisDerivRef = qCopy->_basisDerivRef;
+ const double_array& basisDerivRef = qCopy._basisDerivRef;
size = numBasisE * numQuadPtsE * spaceDimE;
CPPUNIT_ASSERT_EQUAL(size, basisDerivRef.size());
for (int i=0; i < size; ++i)
CPPUNIT_ASSERT_EQUAL(basisDerivE[i], basisDerivRef[i]);
- const double_array& quadPtsRef = qCopy->_quadPtsRef;
+ const double_array& quadPtsRef = qCopy._quadPtsRef;
size = numQuadPtsE * cellDimE;
CPPUNIT_ASSERT_EQUAL(size, quadPtsRef.size());
for (int i=0; i < size; ++i)
CPPUNIT_ASSERT_EQUAL(quadPtsRefE[i], quadPtsRef[i]);
- const double_array& quadWts = qCopy->quadWts();
+ const double_array& quadWts = qCopy.quadWts();
size = numQuadPtsE;
CPPUNIT_ASSERT_EQUAL(size, quadWts.size());
for (int i=0; i < size; ++i)
CPPUNIT_ASSERT_EQUAL(quadWtsE[i], quadWts[i]);
- const double_array& quadPts = qCopy->quadPts();
- size = 1;
- CPPUNIT_ASSERT_EQUAL(size, quadPts.size());
- for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_EQUAL(quadPtsE[i], quadPts[i]);
+ CPPUNIT_ASSERT_EQUAL(geometry.cellDim(), qCopy.refGeometry().cellDim());
+ CPPUNIT_ASSERT_EQUAL(geometry.spaceDim(), qCopy.refGeometry().spaceDim());
+ CPPUNIT_ASSERT_EQUAL(geometry.numCorners(), qCopy.refGeometry().numCorners());
+} // testCopyConstructor
- const double_array& jacobian = qCopy->_jacobian;
- size = 1;
- CPPUNIT_ASSERT_EQUAL(size, jacobian.size());
- for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_EQUAL(jacobianE[i], jacobian[i]);
-
- const double_array& jacobianInv = qCopy->jacobianInv();
- size = 1;
- CPPUNIT_ASSERT_EQUAL(size, jacobianInv.size());
- for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_EQUAL(jacobianInvE[i], jacobianInv[i]);
-
- const double_array& jacobianDet = qCopy->jacobianDet();
- size = 1;
- CPPUNIT_ASSERT_EQUAL(size, jacobianDet.size());
- for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_EQUAL(jacobianDetE[i], jacobianDet[i]);
-
- CPPUNIT_ASSERT(0 != qCopy->_geometry);
- CPPUNIT_ASSERT_EQUAL(geometry.cellDim(), qCopy->_geometry->cellDim());
- CPPUNIT_ASSERT_EQUAL(geometry.spaceDim(), qCopy->_geometry->spaceDim());
- CPPUNIT_ASSERT_EQUAL(geometry.numCorners(), qCopy->_geometry->numCorners());
-
- delete qCopy; qCopy = 0;
-} // testClone
-
// ----------------------------------------------------------------------
// Test checkConditioning()
void
pylith::feassemble::TestQuadrature::testCheckConditioning(void)
{ // testCheckConditioning
- Quadrature1D<topology::Mesh> q;
+ Quadrature<topology::Mesh> q;
+
CPPUNIT_ASSERT_EQUAL(false, q.checkConditioning());
q.checkConditioning(true);
CPPUNIT_ASSERT_EQUAL(true, q.checkConditioning());
@@ -172,92 +112,158 @@
CPPUNIT_ASSERT_EQUAL(false, q.checkConditioning());
} // testCheckConditioning
-#if 0
// ----------------------------------------------------------------------
-// Test computeGeometry() and retrieveGeometry() for meshes.
+// Test quadPts(), basisDeriv(), jacobian(), and jacobianDet().
void
-pylith::feassemble::TestQuadrature::_testComputeGeometry(Quadrature<topology::Mesh>* pQuad,
- const QuadratureData& data) const
+pylith::feassemble::TestQuadrature::testEngineAccessors(void)
+{ // testEngineAccessors
+ const int cellDim = 2;
+ const int numBasis = 5;
+ const int numQuadPts = 1;
+ const int spaceDim = 3;
+ const double basis[] = {
+ 1.1, 1.2, 1.3, 1.4, 1.5
+ };
+ const double basisDerivRef[] = {
+ 2.1, 2.2, 2.3,
+ 2.4, 2.5, 2.6,
+ 2.7, 2.8, 2.9,
+ 2.10, 2.11, 2.12,
+ 2.13, 2.14, 2.15,
+ };
+ const double quadPtsRef[] = { 3.1, 3.2, 3.3 };
+ const double quadWts[] = { 4.0 };
+
+ QuadratureRefCell refCell;
+ refCell.initialize(basis, basisDerivRef, quadPtsRef, quadWts,
+ cellDim, numBasis, numQuadPts, spaceDim);
+
+ Quadrature1D engine(refCell);
+ engine.initialize();
+
+ Quadrature<topology::Mesh> q;
+ q._engine = engine.clone();
+
+ size_t size = 0;
+
+ size = numQuadPts * spaceDim;
+ CPPUNIT_ASSERT_EQUAL(size, q.quadPts().size());
+
+ size = numQuadPts * cellDim * spaceDim;
+ CPPUNIT_ASSERT_EQUAL(size, q.jacobian().size());
+
+ size = numQuadPts;
+ CPPUNIT_ASSERT_EQUAL(size, q.jacobianDet().size());
+
+ size = numQuadPts * numBasis * spaceDim;
+ CPPUNIT_ASSERT_EQUAL(size, q.basisDeriv().size());
+} // testEngineAccessors
+
+// ----------------------------------------------------------------------
+// Test computeGeometry() and retrieveGeometry()
+void
+pylith::feassemble::TestQuadrature::testComputeGeometry(void)
{ // testComputeGeometry
+ typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+
+ QuadratureData2DLinear data;
const int cellDim = data.cellDim;
const int numBasis = data.numBasis;
const int numQuadPts = data.numQuadPts;
const int spaceDim = data.spaceDim;
- const double* basis = data.basis;
- const double* basisDerivRef = data.basisDerivRef;
- const double* basisDeriv = data.basisDeriv;
- const double* quadPtsRef = data.quadPtsRef;
- const double* quadWts = data.quadWts;
- const int numVertices = data.numVertices;
const int numCells = data.numCells;
const double* vertCoords = data.vertices;
- const int* cells = data.cells;
- const double* quadPts = data.quadPts;
- const double* jacobian = data.jacobian;
- const double* jacobianInv = data.jacobianInv;
- const double* jacobianDet = data.jacobianDet;
+ const double* quadPtsE = data.quadPts;
+ const double* jacobianE = data.jacobian;
+ const double* jacobianDetE = data.jacobianDet;
+ const double* basisDerivE = data.basisDeriv;
const double minJacobian = 1.0e-06;
- pQuad->minJacobian(minJacobian);
- pQuad->initialize(basis, basisDerivRef, quadPtsRef, quadWts,
- cellDim, numBasis, numQuadPts, spaceDim);
-
// Create mesh with test cell
- ALE::Obj<Mesh> mesh = new Mesh(PETSC_COMM_WORLD, cellDim);
- CPPUNIT_ASSERT(!mesh.isNull());
- ALE::Obj<sieve_type> sieve = new sieve_type(mesh->comm());
+ topology::Mesh mesh(data.cellDim);
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ ALE::Obj<SieveMesh::sieve_type> sieve =
+ new SieveMesh::sieve_type(mesh.comm());
CPPUNIT_ASSERT(!sieve.isNull());
+ // Cells and vertices
const bool interpolate = false;
- ALE::Obj<ALE::Mesh::sieve_type> s = new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
-
+ ALE::Obj<ALE::Mesh::sieve_type> s =
+ new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
+
ALE::SieveBuilder<ALE::Mesh>::buildTopology(s, cellDim, numCells,
- (int*) cells, numVertices, interpolate, numBasis);
- std::map<Mesh::point_type,Mesh::point_type> renumbering;
+ const_cast<int*>(data.cells),
+ data.numVertices,
+ interpolate, numBasis);
+ std::map<ALE::Mesh::point_type,ALE::Mesh::point_type> renumbering;
ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
- mesh->setSieve(sieve);
- mesh->stratify();
- ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, spaceDim, vertCoords);
+ sieveMesh->setSieve(sieve);
+ sieveMesh->stratify();
+ ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, spaceDim,
+ data.vertices);
+
+ // Setup quadrature and compute geometry
+ GeometryTri2D geometry;
+ Quadrature<topology::Mesh> quadrature;
+ quadrature.refGeometry(&geometry);
+ quadrature.minJacobian(minJacobian);
+ quadrature.initialize(data.basis, data.basisDerivRef,
+ data.quadPtsRef, data.quadWts,
+ cellDim, numBasis, numQuadPts, spaceDim);
+
+ const ALE::Obj<SieveMesh::label_sequence>& cells = sieveMesh->heightStratum(0);
+ CPPUNIT_ASSERT(!cells.isNull());
+ quadrature.computeGeometry(mesh, cells);
+ size_t size = 0;
+
// Check values from computeGeometry()
- const ALE::Obj<Mesh::label_sequence>& cellsMesh = mesh->heightStratum(0);
- CPPUNIT_ASSERT(!cellsMesh.isNull());
- const Mesh::label_sequence::iterator e_iter = cellsMesh->begin();
- const ALE::Obj<Mesh::real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- CPPUNIT_ASSERT(!coordinates.isNull());
- pQuad->computeGeometry(mesh, coordinates, *e_iter);
-
- CPPUNIT_ASSERT(1 == numCells);
-
const double tolerance = 1.0e-06;
- int size = numQuadPts * spaceDim;
- for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(quadPts[i], pQuad->_quadPts[i], tolerance);
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ quadrature.retrieveGeometry(*c_iter);
- size = numQuadPts * cellDim * spaceDim;
- for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(jacobian[i], pQuad->_jacobian[i],
- tolerance);
+ const double_array& quadPts = quadrature.quadPts();
+ size = numQuadPts * spaceDim;
+ CPPUNIT_ASSERT_EQUAL(size, quadPts.size());
+ for (size_t i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(quadPtsE[i], quadPts[i], tolerance);
- size = numQuadPts * spaceDim * cellDim;
- for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(jacobianInv[i], pQuad->_jacobianInv[i],
- tolerance);
+ const double_array& jacobian = quadrature.jacobian();
+ size = numQuadPts * cellDim * spaceDim;
+ CPPUNIT_ASSERT_EQUAL(size, jacobian.size());
+ for (size_t i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(jacobianE[i], jacobian[i], tolerance);
- size = numQuadPts;
- for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(jacobianDet[i], pQuad->_jacobianDet[i],
- tolerance);
+ const double_array& jacobianDet = quadrature.jacobianDet();
+ size = numQuadPts;
+ CPPUNIT_ASSERT_EQUAL(size, jacobianDet.size());
+ for (size_t i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(jacobianDetE[i], jacobianDet[i], tolerance);
+
+ const double_array& basisDeriv = quadrature.basisDeriv();
+ size = numQuadPts * numBasis * spaceDim;
+ CPPUNIT_ASSERT_EQUAL(size, basisDeriv.size());
+ for (size_t i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(basisDerivE[i], basisDeriv[i], tolerance);
+ } // for
- size = numQuadPts * numBasis * spaceDim;
- for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(basisDeriv[i], pQuad->_basisDeriv[i],
- tolerance);
+ // Check clear()
+ quadrature.clear();
+
+ CPPUNIT_ASSERT(0 == quadrature._quadPtsField);
+ CPPUNIT_ASSERT(0 == quadrature._jacobianField);
+ CPPUNIT_ASSERT(0 == quadrature._jacobianDetField);
+ CPPUNIT_ASSERT(0 == quadrature._basisDerivField);
+ CPPUNIT_ASSERT(0 == quadrature._engine);
+ quadrature.clear();
} // 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-13 19:52:03 UTC (rev 14048)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadrature.hh 2009-02-13 23:05:03 UTC (rev 14049)
@@ -23,14 +23,10 @@
#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 TestQuadrature;
- class QuadratureData;
} // feassemble
} // pylith
@@ -41,30 +37,27 @@
// CPPUNIT TEST SUITE /////////////////////////////////////////////////
CPPUNIT_TEST_SUITE( TestQuadrature );
- CPPUNIT_TEST( testClone );
+ CPPUNIT_TEST( testCopyConstructor );
CPPUNIT_TEST( testCheckConditioning );
+ CPPUNIT_TEST( testEngineAccessors );
+ CPPUNIT_TEST( testComputeGeometry );
CPPUNIT_TEST_SUITE_END();
// PUBLIC METHODS /////////////////////////////////////////////////////
public :
- /// Test clone()
- void testClone(void);
+ /// Test copy constructor.
+ void testCopyConstructor(void);
- void testCheckConditioning(void);
/// Test checkConditioning()
+ void testCheckConditioning(void);
- // PROTECTED METHODS //////////////////////////////////////////////////
-protected :
+ /// Test quadPts(), basisDeriv(), jacobian(), and jacobianDet().
+ void testEngineAccessors(void);
- /** Test computeGeometry() and retrieveGeometry().
- *
- * @param pQuad Pointer to quadrature
- * @param data Data for testing quadrature
- */
- void _testComputeGeometry(Quadrature<topology::Mesh>* pQuad,
- const QuadratureData& data) const;
+ /// Test computeGeometry(), retrieveGeometry(), and clear().
+ void testComputeGeometry(void);
}; // class TestQuadrature
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadratureEngine.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadratureEngine.cc 2009-02-13 19:52:03 UTC (rev 14048)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestQuadratureEngine.cc 2009-02-13 23:05:03 UTC (rev 14049)
@@ -129,6 +129,9 @@
size = numQuadPts;
CPPUNIT_ASSERT_EQUAL(size, engine.jacobianDet().size());
+
+ size = numQuadPts * numBasis * spaceDim;
+ CPPUNIT_ASSERT_EQUAL(size, engine.basisDeriv().size());
} // testInitialize
// ----------------------------------------------------------------------
More information about the CIG-COMMITS
mailing list