[cig-commits] r5110 - in short/3D/PyLith/trunk: libsrc/feassemble
unittests/libtests/feassemble unittests/libtests/feassemble/data
brad at geodynamics.org
brad at geodynamics.org
Fri Oct 27 17:15:26 PDT 2006
Author: brad
Date: 2006-10-27 17:15:25 -0700 (Fri, 27 Oct 2006)
New Revision: 5110
Added:
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.hh
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorInertia1D.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorInertia1D.hh
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorApp.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.hh
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DLinear.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DLinear.hh
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DQuadratic.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DQuadratic.hh
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia1DLinear.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia1DQuadratic.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/feutils.py
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature1DLinear.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature1DQuadratic.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature2DLinear.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature2DQuadratic.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature3DLinear.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature3DQuadratic.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureApp.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXZ.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/generate.sh
Log:
Started work on unit tests for integration of inertia term(1-D tests are in place). Basic testing infrastructure for integrators is in place, except for checking matrices.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc 2006-10-28 00:15:25 UTC (rev 5110)
@@ -19,7 +19,9 @@
// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::Integrator::Integrator(void) :
- _quadrature(0)
+ _quadrature(0),
+ _cellVector(0),
+ _cellMatrix(0)
{ // constructor
} // constructor
@@ -28,11 +30,15 @@
pylith::feassemble::Integrator::~Integrator(void)
{ // destructor
delete _quadrature; _quadrature = 0;
+ delete[] _cellVector; _cellVector = 0;
+ delete[] _cellMatrix; _cellMatrix = 0;
} // destructor
// ----------------------------------------------------------------------
// Copy constructor
-pylith::feassemble::Integrator::Integrator(const Integrator& i)
+pylith::feassemble::Integrator::Integrator(const Integrator& i) :
+ _cellVector(0),
+ _cellMatrix(0)
{ // copy constructor
if (0 != i._quadrature)
_quadrature = i._quadrature->clone();
@@ -45,6 +51,10 @@
{ // 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
// ----------------------------------------------------------------------
@@ -54,7 +64,8 @@
{ // _initCellVector
assert(0 != _quadrature);
const int size = _quadrature->spaceDim() * _quadrature->numCorners();
- _cellVector = (size > 0) ? new real_section_type::value_type[size] : 0;
+ 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
@@ -79,7 +90,8 @@
const int size =
_quadrature->spaceDim() * _quadrature->numCorners() *
_quadrature->spaceDim() * _quadrature->numCorners();
- _cellMatrix = (size > 0) ? new real_section_type::value_type[size] : 0;
+ 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
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc 2006-10-28 00:15:25 UTC (rev 5110)
@@ -82,8 +82,8 @@
const int numBasis = _quadrature->numCorners();
const int spaceDim = _quadrature->spaceDim();
- // FIX THIS
- // Hardwire mass density
+ // :TODO: Get mass density at quadrature points from material database
+ // For now, hardwire mass density
const double density = 1.0;
// Compute action for cell
@@ -91,18 +91,19 @@
const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density;
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const int iBlock = iBasis * spaceDim;
- double val = wt*basis[iQ+iBasis];
+ const double valI = wt*basis[iQ+iBasis];
for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- val *= basis[iQ+jBasis];
+ const int jBlock = jBasis * spaceDim;
+ const double val = valI * basis[iQ+jBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
- _cellVector[iBlock+iDim] += val * fieldInCell[iBlock+iDim];
+ _cellVector[iBlock+iDim] += val * fieldInCell[jBlock+iDim];
} // for
} // for
} // for
PetscErrorCode err =
PetscLogFlops(numQuadPts*(2+numBasis*(1+numBasis*(1+2*spaceDim))));
if (err)
- throw std::runtime_error("Couldn't log PETSc flops.");
+ throw std::runtime_error("Logging PETSc flops failed.");
// Assemble cell contribution into field
fieldOut->updateAdd(patch, *cellIter, _cellVector);
@@ -119,6 +120,9 @@
assert(0 != mat);
assert(0 != _quadrature);
+ // Setup symmetric, sparse matrix
+ // :TODO: ADD STUFF HERE
+
// Get information about section
const topology_type::patch_type patch = 0;
const ALE::Obj<topology_type>& topology = coordinates->getTopology();
@@ -147,8 +151,8 @@
const int numBasis = _quadrature->numCorners();
const int spaceDim = _quadrature->spaceDim();
- // FIX THIS
- // Hardwire mass density
+ // :TODO: Get mass density at quadrature points from material database
+ // For now, hardwire mass density
const double density = 1.0;
// Integrate cell
@@ -156,10 +160,10 @@
const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density;
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const int iBlock = iBasis * spaceDim;
- double val = wt*basis[iQ+iBasis];
+ const double valI = wt*basis[iQ+iBasis];
for (int jBasis=0; jBasis < numBasis; ++jBasis) {
const int jBlock = jBasis * spaceDim;
- val *= basis[iQ+jBasis];
+ const double val = valI * basis[iQ+jBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
_cellMatrix[(iBlock+iDim)*(numBasis*spaceDim)+jBlock+iDim] += val;
} // for
@@ -168,10 +172,10 @@
PetscErrorCode err =
PetscLogFlops(numQuadPts*(2+numBasis*(1+numBasis*(1+spaceDim))));
if (err)
- throw std::runtime_error("Couldn't log PETSc flops.");
+ throw std::runtime_error("Logging PETSc flops failed.");
// Assemble cell contribution into sparse matrix
- // STUFF GOES HERE
+ // :TODO: ADD STUFF HERE
} // for
} // integrate
@@ -238,9 +242,9 @@
_cellVector[iBlock+iDim] *= diagScale;
} // for
PetscErrorCode err =
- PetscLogFlops(numQuadPts*(2+numBasis*2) + numBasis);
+ PetscLogFlops(numQuadPts*(3+numBasis*(3+spaceDim)) + numBasis*spaceDim);
if (err)
- throw std::runtime_error("Couldn't log PETSc flops.");
+ throw std::runtime_error("Logging PETSc flops failed.");
// Assemble cell contribution into field
fieldOut->updateAdd(patch, *cellIter, _cellVector);
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh 2006-10-28 00:15:25 UTC (rev 5110)
@@ -14,6 +14,21 @@
* @file pylith/feassemble/IntegratorInertia.hh
*
* @brief Integrate inertial term for 3-D finite elements.
+ *
+ * This object performs integrals over the domain of a finite-element
+ * associated with translational inertia.
+ * - Integrate action over cell
+ * \f[
+ * \int_{V^e} \rho N^p \sum_q N^q u_i^q \, dV
+ * \f]
+ * - Integrate to form matrix
+ * \f[
+ * \int_{V^e} (\rho N^q N^q)_i \, dV
+ * \f]
+ * - Integrate and lump to form lumped matrix (field)
+ *
+ * See governing equations section of user manual for more
+ * information.
*/
#if !defined(pylith_feassemble_integratorinertia_hh)
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc 2006-10-28 00:15:25 UTC (rev 5110)
@@ -54,10 +54,10 @@
_resetGeometry();
// Get coordinates of cell's vertices
- const topology_type::patch_type patch = 0;
+ const topology_type::patch_type patch = 0;
const real_section_type::value_type* vertCoords =
coordinates->restrict(patch, cell);
- //assert(3 == coordinates.GetFiberDimensionByDepth(patch,
+ //assert(3 == coordinates->getFiberDimensionByDepth(patch,
//*vertices->begin(), 0));
// Loop over quadrature points
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am 2006-10-28 00:15:25 UTC (rev 5110)
@@ -17,14 +17,33 @@
check_PROGRAMS = testfeassemble
+# Primary source files
testfeassemble_SOURCES = \
+ TestIntegrator.cc \
+ TestIntegratorInertia1D.cc \
TestQuadrature.cc \
TestQuadrature1D.cc \
TestQuadrature1Din2D.cc \
TestQuadrature2D.cc \
TestQuadrature2Din3D.cc \
TestQuadrature3D.cc \
- test_feassemble.cc \
+ test_feassemble.cc
+
+noinst_HEADERS = \
+ TestIntegrator.hh \
+ TestIntegratorInertia1D.hh \
+ TestQuadrature.hh \
+ TestQuadrature1D.hh \
+ TestQuadrature1Din2D.hh \
+ TestQuadrature1Din3D.hh \
+ TestQuadrature2D.hh \
+ TestQuadrature2Din3D.hh
+
+# Source files associated with testing data
+testfeassemble_SOURCES += \
+ data/IntegratorData.cc \
+ data/IntegratorDataInertia1DLinear.cc \
+ data/IntegratorDataInertia1DQuadratic.cc \
data/QuadratureData.cc \
data/QuadratureData1DLinear.cc \
data/QuadratureData1DQuadratic.cc \
@@ -42,13 +61,10 @@
data/QuadratureData3DLinear.cc \
data/QuadratureData3DQuadratic.cc
-noinst_HEADERS = \
- TestQuadrature.hh \
- TestQuadrature1D.hh \
- TestQuadrature1Din2D.hh \
- TestQuadrature1Din3D.hh \
- TestQuadrature2D.hh \
- TestQuadrature2Din3D.hh \
+noinst_HEADERS += \
+ data/IntegratorData.hh \
+ data/IntegratorDataInertia1DLinear.hh \
+ data/IntegratorDataInertia1DQuadratic.hh \
data/QuadratureData.hh \
data/QuadratureData1DLinear.hh \
data/QuadratureData1DQuadratic.hh \
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.cc 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.cc 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,206 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestIntegrator.hh" // Implementation of class methods
+
+#include "pylith/feassemble/IntegratorInertia.hh" // USES IntegratorInertia
+#include "pylith/feassemble/Quadrature1D.hh" // USES Quadrature1D
+#include "data/IntegratorData.hh" // USES IntegratorData
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::feassemble::TestIntegrator );
+
+// ----------------------------------------------------------------------
+namespace pylith {
+ namespace feassemble {
+ class _TestIntegrator;
+ } // feassemble
+} // pylith
+
+/// Helper class for TestIntegrator
+class pylith::feassemble::_TestIntegrator {
+
+public :
+ /** Setup mesh.
+ *
+ * @param data Integrator data
+ */
+ static
+ ALE::Obj<ALE::Mesh>
+ _setupMesh(const IntegratorData& data);
+}; // _TestIntegrator
+
+// ----------------------------------------------------------------------
+// Test clone().
+void
+pylith::feassemble::TestIntegrator::testClone(void)
+{ // testClone
+ // Test cloning by testing value of minJacobian value in quadrature
+
+ Quadrature1D quadrature;
+ const double minJacobian = 4.0;
+ quadrature.minJacobian(minJacobian);
+
+ IntegratorInertia iOrig;
+ iOrig.quadrature(&quadrature);
+
+ Integrator* iCopy = iOrig.clone();
+ CPPUNIT_ASSERT_EQUAL(minJacobian, iCopy->_quadrature->minJacobian());
+
+ delete iCopy;
+} // testClone
+
+// ----------------------------------------------------------------------
+// Test quadrature().
+void
+pylith::feassemble::TestIntegrator::testQuadrature(void)
+{ // testQuadrature
+ // Since quadrature is cloned, test setting quadrature by testing
+ // value of minJacobian
+
+ Quadrature1D quadrature;
+ const double minJacobian = 4.0;
+ quadrature.minJacobian(minJacobian);
+
+ IntegratorInertia integrator;
+ integrator.quadrature(&quadrature);
+
+ CPPUNIT_ASSERT_EQUAL(minJacobian, integrator._quadrature->minJacobian());
+} // testQuadrature
+
+
+// ----------------------------------------------------------------------
+// Test integrateAction()
+void
+pylith::feassemble::TestIntegrator::_testIntegrateAction(Integrator* integrator,
+ const IntegratorData& data) const
+{ // _testIntegrateAction
+ typedef ALE::Mesh::real_section_type real_section_type;
+ typedef ALE::Mesh::topology_type topology_type;
+
+ ALE::Obj<ALE::Mesh> mesh = _TestIntegrator::_setupMesh(data);
+ const ALE::Mesh::int_section_type::patch_type patch = 0;
+
+ // Fiber dimension (number of values in field per vertex) for fields
+ const int fiberDim = data.fiberDim;
+
+ // Setup input field for action
+ const ALE::Obj<real_section_type>& fieldIn =
+ mesh->getRealSection("fieldIn");
+ fieldIn->setName("fieldIn");
+ fieldIn->setFiberDimensionByDepth(patch, 0, fiberDim);
+ fieldIn->allocate();
+ int iVertex = 0;
+ const ALE::Obj<topology_type::label_sequence>& vertices =
+ mesh->getTopology()->depthStratum(patch, 0);
+ const topology_type::label_sequence::iterator verticesEnd =
+ vertices->end();
+ for (topology_type::label_sequence::iterator vIter=vertices->begin();
+ vIter != verticesEnd;
+ ++vIter, ++iVertex)
+ fieldIn->update(patch, *vIter, &data.fieldIn[iVertex*fiberDim]);
+
+ // Setup field for action result
+ const ALE::Obj<real_section_type>& fieldOut =
+ mesh->getRealSection("fieldOut");
+ fieldOut->setName("fieldOut");
+ fieldOut->setFiberDimensionByDepth(patch, 0, fiberDim);
+ fieldOut->allocate();
+
+ // Integrate action
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ integrator->integrateAction(fieldOut, fieldIn, coordinates);
+
+ // Check values in output field
+ iVertex = 0;
+ const double tolerance = 1.0e-06;
+ for (topology_type::label_sequence::iterator vIter=vertices->begin();
+ vIter != verticesEnd;
+ ++vIter, ++iVertex) {
+ const real_section_type::value_type* vals =
+ fieldOut->restrict(patch, *vIter);
+ const double* valsE = &data.valsAction[iVertex*fiberDim];
+ const int dim = fieldOut->getFiberDimension(patch, *vIter);
+ CPPUNIT_ASSERT_EQUAL(fiberDim, dim);
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[iDim]/valsE[iDim], tolerance);
+ } // for
+} // _testIntegrateAction
+
+// ----------------------------------------------------------------------
+// Test integrate()
+void
+pylith::feassemble::TestIntegrator::_testIntegrate(Integrator* integrator,
+ const IntegratorData& data) const
+{ // _testIntegrate
+ typedef ALE::Mesh::real_section_type real_section_type;
+
+ ALE::Obj<ALE::Mesh> mesh = _TestIntegrator::_setupMesh(data);
+
+ // Integrate
+ PetscMat mat;
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ integrator->integrate(&mat, coordinates);
+
+ // Crate dense matrix
+ // :TODO: ADD STUFF HERE
+
+ // Get values associated with dense matrix
+ // :TODO: ADD STUFF HERE
+
+ // :TODO: Check values from dense matrix
+ // ADD STUFF HERE
+
+ CPPUNIT_ASSERT(false);
+} // _testIntegrate
+
+// ----------------------------------------------------------------------
+// Setup mesh.
+ALE::Obj<ALE::Mesh>
+pylith::feassemble::_TestIntegrator::_setupMesh(const IntegratorData& data)
+{ // _setupMesh
+ typedef ALE::Mesh::topology_type topology_type;
+ typedef topology_type::sieve_type sieve_type;
+
+ const int cellDim = data.cellDim;
+ const int numCorners = data.numCorners;
+ const int spaceDim = data.spaceDim;
+ const int numVertices = data.numVertices;
+ const int numCells = data.numCells;
+ const double* vertCoords = data.vertices;
+ const int* cells = data.cells;
+ CPPUNIT_ASSERT(0 != vertCoords);
+ CPPUNIT_ASSERT(0 != cells);
+
+ ALE::Obj<ALE::Mesh> mesh = new ALE::Mesh(PETSC_COMM_WORLD, cellDim);
+ ALE::Obj<sieve_type> sieve = new sieve_type(mesh->comm());
+ ALE::Obj<topology_type> topology = new topology_type(mesh->comm());
+
+ const bool interpolate = false;
+ ALE::New::SieveBuilder<sieve_type>::buildTopology(sieve, cellDim, numCells,
+ const_cast<int*>(cells), numVertices, interpolate, numCorners);
+ sieve->stratify();
+ topology->setPatch(0, sieve);
+ topology->stratify();
+ mesh->setTopology(topology);
+ ALE::New::SieveBuilder<sieve_type>::buildCoordinates(
+ mesh->getRealSection("coordinates"), spaceDim, vertCoords);
+
+ return mesh;
+} // _setupMesh
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.hh 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegrator.hh 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,78 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/feassemble/TestIntegrator.hh
+ *
+ * @brief C++ TestIntegrator object
+ *
+ * C++ unit testing for Integrator.
+ */
+
+#if !defined(pylith_feassemble_testintegrator_hh)
+#define pylith_feassemble_testintegrator_hh
+
+#include <cppunit/extensions/HelperMacros.h>
+
+/// Namespace for spatialdata package
+namespace pylith {
+ namespace feassemble {
+ class Integrator;
+ class TestIntegrator;
+ class IntegratorData;
+ } // feassemble
+} // pylith
+
+/// C++ unit testing for Integrator
+class pylith::feassemble::TestIntegrator : public CppUnit::TestFixture
+{ // class TestIntegrator
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestIntegrator );
+ CPPUNIT_TEST( testClone );
+ CPPUNIT_TEST( testQuadrature );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Test clone()
+ void testClone(void);
+
+ /// Test quadrature()
+ void testQuadrature(void);
+
+ // PROTECTED METHODS //////////////////////////////////////////////////
+protected :
+
+ /** Test integrateAction()
+ *
+ * @param integrator Pointer to integrator
+ * @param data Data for testing integrator
+ */
+ void _testIntegrateAction(Integrator* integrator,
+ const IntegratorData& data) const;
+
+ /** Test integrate()
+ *
+ * @param integrator Pointer to integrator
+ * @param data Data for testing integrator
+ */
+ void _testIntegrate(Integrator* integrator,
+ const IntegratorData& data) const;
+
+}; // class TestIntegrator
+
+#endif // pylith_feassemble_testintegrator_hh
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorInertia1D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorInertia1D.cc 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorInertia1D.cc 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,84 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestIntegratorInertia1D.hh" // Implementation of class methods
+
+#include "pylith/feassemble/IntegratorInertia.hh"
+#include "pylith/feassemble/Quadrature1D.hh"
+
+#include "data/IntegratorDataInertia1DLinear.hh"
+#include "data/IntegratorDataInertia1DQuadratic.hh"
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::feassemble::TestIntegratorInertia1D );
+
+// ----------------------------------------------------------------------
+// Test constructor
+void
+pylith::feassemble::TestIntegratorInertia1D::testConstructor(void)
+{ // testConstructor
+ IntegratorInertia integrator;
+} // testConstructor
+
+// ----------------------------------------------------------------------
+// Test integrate() & integrateAction() w/linear basis fns
+void
+pylith::feassemble::TestIntegratorInertia1D::testLinear(void)
+{ // testLinear
+ IntegratorDataInertia1DLinear data;
+
+ Quadrature1D q;
+ q.initialize(data.basis,
+ data.basisDeriv,
+ data.quadPts,
+ data.quadWts,
+ data.cellDim,
+ data.numCorners,
+ data.numQuadPts,
+ data.spaceDim);
+
+ IntegratorInertia integrator;
+ integrator.quadrature(&q);
+
+ _testIntegrateAction(&integrator, data);
+
+ _testIntegrate(&integrator, data);
+} // testLinear
+
+// ----------------------------------------------------------------------
+// Test integrate() & integrateAction() w/quadratic basis fns
+void
+pylith::feassemble::TestIntegratorInertia1D::testQuadratic(void)
+{ // testQuadratic
+ IntegratorDataInertia1DQuadratic data;
+
+ Quadrature1D q;
+ q.initialize(data.basis,
+ data.basisDeriv,
+ data.quadPts,
+ data.quadWts,
+ data.cellDim,
+ data.numCorners,
+ data.numQuadPts,
+ data.spaceDim);
+
+ IntegratorInertia integrator;
+ integrator.quadrature(&q);
+
+ _testIntegrateAction(&integrator, data);
+
+ _testIntegrate(&integrator, data);
+} // testQuadratic
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorInertia1D.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorInertia1D.hh 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorInertia1D.hh 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,60 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/feassemble/TestIntegratorInertia1D.hh
+ *
+ * @brief C++ TestIntegratorInertia1D object
+ *
+ * C++ unit testing for Quadrature1D.
+ */
+
+#if !defined(pylith_feassemble_testintegratorinertia1d_hh)
+#define pylith_feassemble_testintegratorinertia1d_hh
+
+#include "TestIntegrator.hh"
+
+/// Namespace for spatialdata package
+namespace pylith {
+ namespace feassemble {
+ class TestIntegratorInertia1D;
+ } // feassemble
+} // pylith
+
+/// C++ unit testing for IntegratorInertia1D
+class pylith::feassemble::TestIntegratorInertia1D : public TestIntegrator
+{ // class TestIntegratorInertia1D
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestIntegratorInertia1D );
+ CPPUNIT_TEST( testConstructor );
+ CPPUNIT_TEST( testLinear );
+ CPPUNIT_TEST( testQuadratic );
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Test constructor
+ void testConstructor(void);
+
+ /// Test integrate() & integrateAction() w/linear basis fns
+ void testLinear(void);
+
+ /// Test integrate() & integrateAction() w/quadratic basis fns
+ void testQuadratic(void);
+
+}; // class TestIntegratorInertia1D
+
+#endif // pylith_feassemble_testintegratorinertia1d_hh
+
+// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc 2006-10-28 00:15:25 UTC (rev 5110)
@@ -17,8 +17,6 @@
#include "pylith/feassemble/Quadrature1D.hh" // USES Quadrature1D
#include "data/QuadratureData.hh" // USES QuadratureData
-#include <sstream> // USES std::stringstream
-
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::feassemble::TestQuadrature );
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorApp.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorApp.py 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorApp.py 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,185 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/feassemble/data/IntegratorApp.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ integrator objects.
+
+from pyre.applications.Script import Script
+
+import numpy
+
+# IntegratorApp class
+class IntegratorApp(Script):
+ """
+ Python application for generating C++ data files for testing C++
+ integrator objects.
+ """
+
+ # INVENTORY //////////////////////////////////////////////////////////
+
+ class Inventory(Script.Inventory):
+ """Python object for managing IntegratorApp facilities and properties."""
+
+ ## @class Inventory
+ ## Python object for managing IntegratorApp facilities and properties.
+ ##
+ ## \b Properties
+ ## @li None
+ ##
+ ## \b Facilities
+ ## @li \b data Data manager.
+
+ import pyre.inventory
+
+ from pylith.utils.CppData import CppData
+ data = pyre.inventory.facility("data", factory=CppData)
+ data.meta['tip'] = "Data manager."
+
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="integratorapp"):
+ """
+ Constructor.
+ """
+ Script.__init__(self, name)
+
+ self.numVertices = None
+ self.spaceDim = None
+ self.numCells = None
+ self.cellDim = None
+ self.numCorners = None
+ self.numQuadPts = None
+ self.fiberDim = None
+
+ self.quadPts = None
+ self.quadWts = None
+ self.vertices = None
+ self.cells = None
+
+ self.basis = None
+ self.basisDeriv = None
+ self.fieldIn = None
+ self.valsActions = None
+ self.valsMatrix = None
+ return
+
+
+ def main(self):
+ """
+ Run the application.
+ """
+ self._initialize()
+ self._calculateMatrix()
+ self._calculateAction()
+ self._initData()
+ self.data.write(self.name)
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Set members using inventory.
+ """
+ Script._configure(self)
+ self.data = self.inventory.data
+ return
+
+
+ def _initialize(self):
+ """
+ Get quadrature information.
+ """
+ q = self.quadrature
+ q.calculateBasis()
+ self.spaceDim = q.spaceDim
+ self.numCorners = q.numCorners
+ self.cellDim = q.cellDim
+ self.numQuadPts = q.numQuadPts
+ self.basis = q.basis
+ self.basisDeriv = q.basisDeriv
+ self.quadPts = q.quadPtsRef
+ self.quadWts = q.quadWts
+
+ return
+
+
+ def _calculateMatrix(self):
+ """
+ Calculate matrix associated with integration.
+ """
+ raise NotImplementedError("Not implemented in abstract base class.")
+
+
+ def _calculateAction(self):
+ """
+ Calculate integration action using matrix.
+ """
+ self.valsAction = numpy.dot(self.valsMatrix, self.fieldIn)[:]
+ return
+
+
+ def _initData(self):
+ self.data.addScalar(vtype="int", name="_numVertices",
+ value=self.numVertices,
+ format="%d")
+ self.data.addScalar(vtype="int", name="_spaceDim", value=self.spaceDim,
+ format="%d")
+ self.data.addScalar(vtype="int", name="_numCells", value=self.numCells,
+ format="%d")
+ self.data.addScalar(vtype="int", name="_cellDim", value=self.cellDim,
+ format="%d")
+ self.data.addScalar(vtype="int", name="_numCorners", value=
+ self.numCorners,
+ format="%d")
+ self.data.addScalar(vtype="int", name="_numQuadPts",
+ value=self.numQuadPts,
+ format="%d")
+ self.data.addScalar(vtype="int", name="_fiberDim",
+ value=self.fiberDim,
+ format="%d")
+
+ self.data.addArray(vtype="double", name="_vertices", values=self.vertices,
+ format="%16.8e", ncols=self.spaceDim)
+ self.data.addArray(vtype="int", name="_cells", values=self.cells,
+ format="%8d", ncols=self.numVertices)
+
+ self.data.addArray(vtype="double", name="_quadPts",
+ values=self.quadPts,
+ format="%16.8e", ncols=self.cellDim)
+ self.data.addArray(vtype="double", name="_quadWts", values=self.quadWts,
+ format="%16.8e", ncols=self.numQuadPts)
+
+ self.data.addArray(vtype="double", name="_basis", values=self.basis,
+ format="%16.8e", ncols=self.cellDim)
+ self.data.addArray(vtype="double", name="_basisDeriv",
+ values=self.basisDeriv,
+ format="%16.8e", ncols=self.cellDim)
+
+ self.data.addArray(vtype="double", name="_fieldIn",
+ values=self.fieldIn,
+ format="%16.8e", ncols=self.spaceDim)
+ self.data.addArray(vtype="double", name="_valsAction",
+ values=self.valsAction,
+ format="%16.8e", ncols=self.spaceDim)
+ self.data.addArray(vtype="double", name="_valsMatrix",
+ values=self.valsMatrix,
+ format="%16.8e", ncols=self.spaceDim)
+
+ return
+
+
+# End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.cc 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.cc 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,43 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include "IntegratorData.hh"
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::feassemble::IntegratorData::IntegratorData(void) :
+ numVertices(0),
+ spaceDim(0),
+ numCells(0),
+ cellDim(0),
+ numCorners(0),
+ numQuadPts(0),
+ fiberDim(0),
+ vertices(0),
+ cells(0),
+ quadPts(0),
+ quadWts(0),
+ basis(0),
+ basisDeriv(0),
+ fieldIn(0),
+ valsAction(0),
+ valsMatrix(0)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::IntegratorData::~IntegratorData(void)
+{ // destructor
+} // destructor
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.hh 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.hh 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,69 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_feassemble_integratordata_hh)
+#define pylith_feassemble_integratordata_hh
+
+namespace pylith {
+ namespace feassemble {
+ class IntegratorData;
+ } // pylith
+} // feassemble
+
+class pylith::feassemble::IntegratorData
+{
+
+// PUBLIC METHODS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ IntegratorData(void);
+
+ /// Destructor
+ ~IntegratorData(void);
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public:
+
+ int numVertices; ///< Number of vertices
+ int spaceDim; ///< Number of dimensions in vertex coordinates
+ int numCells; ///< Number of cells
+ int cellDim; ///< Number of dimensions associated with cell
+ int numCorners; ///< Number of vertices in cell
+ int numQuadPts; ///< Number of quadrature points
+ int fiberDim; ///< Number of values per vertex in field
+
+ /// @name Mesh information
+ //@{
+ double* vertices; ///< Pointer to coordinates of vertices
+ int* cells; ///< Pointer to indices of vertices in cells
+ //@}
+
+ /// @name Discretization information
+ //@{
+ double* quadPts; ///< Coordinates of quad pts in ref cell
+ double* quadWts; ///< Weights of quadrature points
+ double* basis; ///< Basis fns at quadrature points
+ double* basisDeriv; ///< Derivatives of basis fns at quad pts
+ //@}
+
+ /// @name Integration information
+ //@{
+ double* fieldIn; ///< Input field for integration action
+ double* valsAction; ///< Expected output for integration action
+ double* valsMatrix; ///< Expected output for integration
+ //@}
+};
+
+#endif // pylith_feassemble_integratordata_hh
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DLinear.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DLinear.cc 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DLinear.cc 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,100 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application itnegratorinertia1dlinear.
+
+#include "IntegratorDataInertia1DLinear.hh"
+
+const int pylith::feassemble::IntegratorDataInertia1DLinear::_numVertices = 2;
+
+const int pylith::feassemble::IntegratorDataInertia1DLinear::_spaceDim = 1;
+
+const int pylith::feassemble::IntegratorDataInertia1DLinear::_numCells = 1;
+
+const int pylith::feassemble::IntegratorDataInertia1DLinear::_cellDim = 1;
+
+const int pylith::feassemble::IntegratorDataInertia1DLinear::_numCorners = 2;
+
+const int pylith::feassemble::IntegratorDataInertia1DLinear::_numQuadPts = 1;
+
+const int pylith::feassemble::IntegratorDataInertia1DLinear::_fiberDim = 1;
+
+const double pylith::feassemble::IntegratorDataInertia1DLinear::_vertices[] = {
+ -2.50000000e-01,
+ 2.00000000e+00,
+};
+
+const int pylith::feassemble::IntegratorDataInertia1DLinear::_cells[] = {
+ 0, 1,
+};
+
+const double pylith::feassemble::IntegratorDataInertia1DLinear::_quadPts[] = {
+ 0.00000000e+00,
+};
+
+const double pylith::feassemble::IntegratorDataInertia1DLinear::_quadWts[] = {
+ 2.00000000e+00,
+};
+
+const double pylith::feassemble::IntegratorDataInertia1DLinear::_basis[] = {
+ 5.00000000e-01,
+ 5.00000000e-01,
+};
+
+const double pylith::feassemble::IntegratorDataInertia1DLinear::_basisDeriv[] = {
+ -5.00000000e-01,
+ 5.00000000e-01,
+};
+
+const double pylith::feassemble::IntegratorDataInertia1DLinear::_fieldIn[] = {
+ 1.20000000e+00,
+ 1.50000000e+00,
+};
+
+const double pylith::feassemble::IntegratorDataInertia1DLinear::_valsAction[] = {
+ 1.51875000e+00,
+ 1.51875000e+00,
+};
+
+const double pylith::feassemble::IntegratorDataInertia1DLinear::_valsMatrix[] = {
+ 5.62500000e-01,
+ 5.62500000e-01,
+ 5.62500000e-01,
+ 5.62500000e-01,
+};
+
+pylith::feassemble::IntegratorDataInertia1DLinear::IntegratorDataInertia1DLinear(void)
+{ // constructor
+ numVertices = _numVertices;
+ spaceDim = _spaceDim;
+ numCells = _numCells;
+ cellDim = _cellDim;
+ numCorners = _numCorners;
+ numQuadPts = _numQuadPts;
+ fiberDim = _fiberDim;
+ vertices = const_cast<double*>(_vertices);
+ cells = const_cast<int*>(_cells);
+ quadPts = const_cast<double*>(_quadPts);
+ quadWts = const_cast<double*>(_quadWts);
+ basis = const_cast<double*>(_basis);
+ basisDeriv = const_cast<double*>(_basisDeriv);
+ fieldIn = const_cast<double*>(_fieldIn);
+ valsAction = const_cast<double*>(_valsAction);
+ valsMatrix = const_cast<double*>(_valsMatrix);
+} // constructor
+
+pylith::feassemble::IntegratorDataInertia1DLinear::~IntegratorDataInertia1DLinear(void)
+{}
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DLinear.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DLinear.hh 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DLinear.hh 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,76 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application itnegratorinertia1dlinear.
+
+#if !defined(pylith_feassemble_integratordatainertia1dlinear_hh)
+#define pylith_feassemble_integratordatainertia1dlinear_hh
+
+#include "IntegratorData.hh"
+
+namespace pylith {
+ namespace feassemble {
+ class IntegratorDataInertia1DLinear;
+ } // pylith
+} // feassemble
+
+class pylith::feassemble::IntegratorDataInertia1DLinear : public IntegratorData
+{
+
+public:
+
+ /// Constructor
+ IntegratorDataInertia1DLinear(void);
+
+ /// Destructor
+ ~IntegratorDataInertia1DLinear(void);
+
+private:
+
+ static const int _numVertices;
+
+ static const int _spaceDim;
+
+ static const int _numCells;
+
+ static const int _cellDim;
+
+ static const int _numCorners;
+
+ static const int _numQuadPts;
+
+ static const int _fiberDim;
+
+ static const double _vertices[];
+
+ static const int _cells[];
+
+ static const double _quadPts[];
+
+ static const double _quadWts[];
+
+ static const double _basis[];
+
+ static const double _basisDeriv[];
+
+ static const double _fieldIn[];
+
+ static const double _valsAction[];
+
+ static const double _valsMatrix[];
+
+};
+
+#endif // pylith_feassemble_integratordatainertia1dlinear_hh
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DQuadratic.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DQuadratic.cc 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DQuadratic.cc 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,117 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application itnegratorinertia1dlinear.
+
+#include "IntegratorDataInertia1DQuadratic.hh"
+
+const int pylith::feassemble::IntegratorDataInertia1DQuadratic::_numVertices = 3;
+
+const int pylith::feassemble::IntegratorDataInertia1DQuadratic::_spaceDim = 1;
+
+const int pylith::feassemble::IntegratorDataInertia1DQuadratic::_numCells = 1;
+
+const int pylith::feassemble::IntegratorDataInertia1DQuadratic::_cellDim = 1;
+
+const int pylith::feassemble::IntegratorDataInertia1DQuadratic::_numCorners = 3;
+
+const int pylith::feassemble::IntegratorDataInertia1DQuadratic::_numQuadPts = 2;
+
+const int pylith::feassemble::IntegratorDataInertia1DQuadratic::_fiberDim = 1;
+
+const double pylith::feassemble::IntegratorDataInertia1DQuadratic::_vertices[] = {
+ -2.50000000e-01,
+ 8.75000000e-01,
+ 2.00000000e+00,
+};
+
+const int pylith::feassemble::IntegratorDataInertia1DQuadratic::_cells[] = {
+ 0, 1, 2,
+};
+
+const double pylith::feassemble::IntegratorDataInertia1DQuadratic::_quadPts[] = {
+ -5.77350269e-01,
+ 5.77350269e-01,
+};
+
+const double pylith::feassemble::IntegratorDataInertia1DQuadratic::_quadWts[] = {
+ 1.00000000e+00, 1.00000000e+00,
+};
+
+const double pylith::feassemble::IntegratorDataInertia1DQuadratic::_basis[] = {
+ 4.55341801e-01,
+ 6.66666667e-01,
+ -1.22008468e-01,
+ -1.22008468e-01,
+ 6.66666667e-01,
+ 4.55341801e-01,
+};
+
+const double pylith::feassemble::IntegratorDataInertia1DQuadratic::_basisDeriv[] = {
+ -1.07735027e+00,
+ 1.15470054e+00,
+ -7.73502692e-02,
+ 7.73502692e-02,
+ -1.15470054e+00,
+ 1.07735027e+00,
+};
+
+const double pylith::feassemble::IntegratorDataInertia1DQuadratic::_fieldIn[] = {
+ 1.20000000e+00,
+ 1.50000000e+00,
+ -8.00000000e-01,
+};
+
+const double pylith::feassemble::IntegratorDataInertia1DQuadratic::_valsAction[] = {
+ 7.75000000e-01,
+ 1.60000000e+00,
+ 2.50000000e-02,
+};
+
+const double pylith::feassemble::IntegratorDataInertia1DQuadratic::_valsMatrix[] = {
+ 2.50000000e-01,
+ 2.50000000e-01,
+ -1.25000000e-01,
+ 2.50000000e-01,
+ 1.00000000e-00,
+ 2.50000000e-01,
+ -1.25000000e-01,
+ 2.50000000e-01,
+ 2.50000000e-01,
+};
+
+pylith::feassemble::IntegratorDataInertia1DQuadratic::IntegratorDataInertia1DQuadratic(void)
+{ // constructor
+ numVertices = _numVertices;
+ spaceDim = _spaceDim;
+ numCells = _numCells;
+ cellDim = _cellDim;
+ numCorners = _numCorners;
+ numQuadPts = _numQuadPts;
+ fiberDim = _fiberDim;
+ vertices = const_cast<double*>(_vertices);
+ cells = const_cast<int*>(_cells);
+ quadPts = const_cast<double*>(_quadPts);
+ quadWts = const_cast<double*>(_quadWts);
+ basis = const_cast<double*>(_basis);
+ basisDeriv = const_cast<double*>(_basisDeriv);
+ fieldIn = const_cast<double*>(_fieldIn);
+ valsAction = const_cast<double*>(_valsAction);
+ valsMatrix = const_cast<double*>(_valsMatrix);
+} // constructor
+
+pylith::feassemble::IntegratorDataInertia1DQuadratic::~IntegratorDataInertia1DQuadratic(void)
+{}
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DQuadratic.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DQuadratic.hh 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorDataInertia1DQuadratic.hh 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,76 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application itnegratorinertia1dlinear.
+
+#if !defined(pylith_feassemble_integratordatainertia1dquadratic_hh)
+#define pylith_feassemble_integratordatainertia1dquadratic_hh
+
+#include "IntegratorData.hh"
+
+namespace pylith {
+ namespace feassemble {
+ class IntegratorDataInertia1DQuadratic;
+ } // pylith
+} // feassemble
+
+class pylith::feassemble::IntegratorDataInertia1DQuadratic : public IntegratorData
+{
+
+public:
+
+ /// Constructor
+ IntegratorDataInertia1DQuadratic(void);
+
+ /// Destructor
+ ~IntegratorDataInertia1DQuadratic(void);
+
+private:
+
+ static const int _numVertices;
+
+ static const int _spaceDim;
+
+ static const int _numCells;
+
+ static const int _cellDim;
+
+ static const int _numCorners;
+
+ static const int _numQuadPts;
+
+ static const int _fiberDim;
+
+ static const double _vertices[];
+
+ static const int _cells[];
+
+ static const double _quadPts[];
+
+ static const double _quadWts[];
+
+ static const double _basis[];
+
+ static const double _basisDeriv[];
+
+ static const double _fieldIn[];
+
+ static const double _valsAction[];
+
+ static const double _valsMatrix[];
+
+};
+
+#endif // pylith_feassemble_integratordatainertia1dquadratic_hh
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia.py 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia.py 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,77 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/feassemble/data/IntegratorInertia.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ IntegratorInertia object.
+
+from IntegratorApp import IntegratorApp
+
+import numpy
+
+# ----------------------------------------------------------------------
+
+# IntegratorInertia class
+class IntegratorInertia(IntegratorApp):
+ """
+ Python application for generating C++ data files for testing C++
+ IntegratorInertia object.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="itnegratorinertia"):
+ """
+ Constructor.
+ """
+ IntegratorApp.__init__(self, name)
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _calculateMatrix(self):
+ """
+ Calculate matrix associated with integration.
+ """
+ import feutils
+
+ self.valsMatrix = numpy.zeros( (self.spaceDim*self.numVertices,
+ self.spaceDim*self.numVertices),
+ dtype=numpy.float64)
+
+ n = numpy.zeros( (self.spaceDim, self.spaceDim*self.numCorners),
+ dtype=numpy.float64)
+
+ for cell in self.cells:
+ cellMatrix = numpy.zeros( (self.spaceDim*self.numCorners,
+ self.spaceDim*self.numCorners),
+ dtype=numpy.float64)
+
+ vertices = self.vertices[cell, :]
+ (jacobian, jacobianInv, jacobianDet) = \
+ feutils.calculateJacobian(self.quadrature, vertices)
+ density = 1.0
+ for iQuad in xrange(self.numQuadPts):
+ n *= 0.0
+ for iBasis in xrange(self.numCorners):
+ for iDim in xrange(self.spaceDim):
+ n[iDim, iBasis*self.spaceDim+iDim] = self.basis[iQuad, iBasis]
+
+ wt = density * self.quadWts[iQuad] * jacobianDet[iQuad]
+ cellMatrix[:] += wt * numpy.dot(n.transpose(), n)
+ feutils.assembleMat(self.valsMatrix, cellMatrix, cell)
+ return
+
+
+# End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia1DLinear.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia1DLinear.py 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia1DLinear.py 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,60 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/feassemble/data/IntegratorInertia1DLinear.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ IntegratorInertia object with 1-D cell and linear basis
+## functions.
+
+from IntegratorInertia import IntegratorInertia
+
+import numpy
+
+# ----------------------------------------------------------------------
+
+# IntegratorInertia1DLinear class
+class IntegratorInertia1DLinear(IntegratorInertia):
+ """
+ Python application for generating C++ data files for testing C++
+ IntegratorInertia object with 1-D cell and linear basis functions.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="itnegratorinertia1dlinear"):
+ """
+ Constructor.
+ """
+ IntegratorInertia.__init__(self, name)
+
+ from Quadrature1DLinear import Quadrature1DLinear
+ self.quadrature = Quadrature1DLinear()
+
+ self.numVertices = 2
+ self.numCells = 1
+ self.fiberDim = 1
+
+ self.vertices = numpy.array( [[-0.25], [2.0]], dtype=numpy.float64)
+ self.cells = numpy.array( [[0, 1]], dtype=numpy.int32)
+ self.fieldIn = numpy.array( [[1.2], [1.5]], dtype=numpy.float64)
+ return
+
+
+# MAIN /////////////////////////////////////////////////////////////////
+if __name__ == "__main__":
+
+ app = IntegratorInertia1DLinear()
+ app.run()
+
+
+# End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia1DQuadratic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia1DQuadratic.py 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia1DQuadratic.py 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,61 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/feassemble/data/IntegratorInertia1DQuadratic.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ IntegratorInertia object with 1-D cell and linear basis
+## functions.
+
+from IntegratorInertia import IntegratorInertia
+
+import numpy
+
+# ----------------------------------------------------------------------
+
+# IntegratorInertia1DQuadratic class
+class IntegratorInertia1DQuadratic(IntegratorInertia):
+ """
+ Python application for generating C++ data files for testing C++
+ IntegratorInertia object with 1-D cell and linear basis functions.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="itnegratorinertia1dlinear"):
+ """
+ Constructor.
+ """
+ IntegratorInertia.__init__(self, name)
+
+ from Quadrature1DQuadratic import Quadrature1DQuadratic
+ self.quadrature = Quadrature1DQuadratic()
+
+ self.numVertices = 3
+ self.numCells = 1
+ self.fiberDim = 1
+
+ self.vertices = numpy.array( [[-0.25], [0.875], [2.0]],
+ dtype=numpy.Float64)
+ self.cells = numpy.array( [[0, 1, 2]], dtype=numpy.Int32)
+ self.fieldIn = numpy.array( [[1.2], [1.5], [-0.8]], dtype=numpy.float64)
+ return
+
+
+# MAIN /////////////////////////////////////////////////////////////////
+if __name__ == "__main__":
+
+ app = IntegratorInertia1DQuadratic()
+ app.run()
+
+
+# End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature1DLinear.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature1DLinear.py 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature1DLinear.py 2006-10-28 00:15:25 UTC (rev 5110)
@@ -65,7 +65,7 @@
# PRIVATE METHODS ////////////////////////////////////////////////////
- def _calculateBasis(self):
+ def calculateBasis(self):
"""
Calculate basis functions and their derivatives at quadrature points.
"""
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature1DQuadratic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature1DQuadratic.py 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature1DQuadratic.py 2006-10-28 00:15:25 UTC (rev 5110)
@@ -75,7 +75,7 @@
# PRIVATE METHODS ////////////////////////////////////////////////////
- def _calculateBasis(self):
+ def calculateBasis(self):
"""
Calculate basis functions and their derivatives at quadrature points.
"""
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature2DLinear.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature2DLinear.py 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature2DLinear.py 2006-10-28 00:15:25 UTC (rev 5110)
@@ -80,9 +80,7 @@
return
- # PRIVATE METHODS ////////////////////////////////////////////////////
-
- def _calculateBasis(self):
+ def calculateBasis(self):
"""
Calculate basis functions and their derivatives at quadrature points.
"""
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature2DQuadratic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature2DQuadratic.py 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature2DQuadratic.py 2006-10-28 00:15:25 UTC (rev 5110)
@@ -116,9 +116,7 @@
return
- # PRIVATE METHODS ////////////////////////////////////////////////////
-
- def _calculateBasis(self):
+ def calculateBasis(self):
"""
Calculate basis functions and their derivatives at quadrature points.
"""
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature3DLinear.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature3DLinear.py 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature3DLinear.py 2006-10-28 00:15:25 UTC (rev 5110)
@@ -103,9 +103,7 @@
return
- # PRIVATE METHODS ////////////////////////////////////////////////////
-
- def _calculateBasis(self):
+ def calculateBasis(self):
"""
Calculate basis functions and their derivatives at quadrature points.
"""
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature3DQuadratic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature3DQuadratic.py 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Quadrature3DQuadratic.py 2006-10-28 00:15:25 UTC (rev 5110)
@@ -201,9 +201,7 @@
return
- # PRIVATE METHODS ////////////////////////////////////////////////////
-
- def _calculateBasis(self):
+ def calculateBasis(self):
"""
Calculate basis functions and their derivatives at quadrature points.
"""
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureApp.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureApp.py 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureApp.py 2006-10-28 00:15:25 UTC (rev 5110)
@@ -29,10 +29,10 @@
# INVENTORY //////////////////////////////////////////////////////////
class Inventory(Script.Inventory):
- """Python object for managing PyLithApp facilities and properties."""
+ """Python object for managing QuadratureApp facilities and properties."""
## @class Inventory
- ## Python object for managing PyLithApp facilities and properties.
+ ## Python object for managing QuadratureApp facilities and properties.
##
## \b Properties
## @li None
@@ -55,20 +55,23 @@
"""
Script.__init__(self, name)
+ # Mesh information
self.numVertices = None
self.spaceDim = None
self.numCells = None
+ self.vertices = None
+ self.cells = None
+
+ # Reference cell information
self.cellDim = None
self.numCorners = None
self.numQuadPts = None
-
self.quadPtsRef = None
self.quadWts = None
- self.vertices = None
- self.cells = None
-
self.basis = None
self.basisDeriv = None
+
+ # Computed quadrature information
self.quadPts = None
self.jacobian = None
self.jacobianDet = None
@@ -80,110 +83,46 @@
"""
Run the application.
"""
- self._calculateBasis()
- self._calculateJacobian()
+ self.calculateBasis()
+ self.calculateJacobian()
+
+ # Quadrature points in cell
+ self.quadPts = numpy.dot(self.basis, self.vertices)
+
self._initData()
self.data.write(self.name)
return
- # PRIVATE METHODS ////////////////////////////////////////////////////
-
- def _configure(self):
+ def calculateBasis(self):
"""
- Set members using inventory.
+ Calculate basis functions and derivatives at quadrature points.
"""
- Script._configure(self)
- self.data = self.inventory.data
- return
+ raise NotImplementedError
- def _calculateJacobian(self):
+ def calculateJacobian(self):
"""
- Calculate jacobian, its determinant, and its inverse at quadrature
+ Calculate Jacobian, its determinant, and its inverse at quadrature
pts plus coordinates of quadrature points in the cell.
"""
- self.jacobian = numpy.zeros( (self.numQuadPts,
- self.cellDim, self.spaceDim),
- dtype=numpy.Float64)
- self.jacobianInv = numpy.zeros( (self.numQuadPts,
- self.spaceDim, self.cellDim),
- dtype=numpy.Float64)
- self.jacobianDet = numpy.zeros( (self.numQuadPts,),
- dtype=numpy.Float64)
-
- iQuad = 0
- for q in self.quadPtsRef:
- # Jacobian at quadrature points
- deriv = self.basisDeriv[iQuad]
- jacobian = numpy.dot(deriv.transpose(), self.vertices)
- self.jacobian[iQuad] = jacobian
+ import feutils
+ (self.jacobian, self.jacobianInv, self.jacobianDet) = \
+ feutils.calculateJacobian(self, self.vertices)
+ return
- # Determinant of Jacobian and Jacobian inverse at quadrature points
- if self.spaceDim == self.cellDim:
- self.jacobianDet[iQuad] = numpy.linalg.det(jacobian)
- self.jacobianInv[iQuad] = numpy.linalg.inv(jacobian)
- else:
- det = numpy.linalg.det(numpy.dot(jacobian, jacobian.transpose()))**0.5
- self.jacobianDet[iQuad] = det
- if 1 == self.cellDim:
- self.jacobianInv[iQuad] = 1.0 / jacobian
- elif 2 == self.cellDim:
- minJacobian = 1.0e-06
- jj01 = jacobian[:,[0,1]]
- jj12 = jacobian[:,[1,2]]
- jj02 = jacobian[:,[0,2]]
- det01 = numpy.linalg.det(jj01)
- det12 = numpy.linalg.det(jj12)
- det02 = numpy.linalg.det(jj02)
- if abs(det01) > minJacobian:
- ij01 = numpy.linalg.inv(jj01)
- if abs(det12) > minJacobian:
- ij12 = numpy.linalg.inv(jj12)
- self.jacobianInv[iQuad] = numpy.array([ [ij01[0,0], ij01[0,1]],
- [ij01[1,0], ij01[1,1]],
- [ij12[1,0], ij12[1,1]] ],
- dtype=numpy.Float64)
- elif abs(det02) > minJacobian:
- ij02 = numpy.linalg.inv(jj02)
- self.jacobianInv[iQuad] = numpy.array([ [ij01[0,0], ij01[0,1]],
- [ij01[1,0], ij01[1,1]],
- [ij02[1,0], ij02[1,1]] ],
- dtype=numpy.Float64)
- else:
- self.jacobianInv[iQuad] = numpy.array([ [ij01[0,0], ij01[0,1]],
- [ij01[1,0], ij01[1,1]],
- [ 0.0, 0.0] ],
- dtype=numpy.Float64)
- elif abs(det02) > minJacobian:
- ij02 = numpy.linalg.inv(jj02)
- if abs(det12) > minJacobian:
- ij12 = numpy.linalg.inv(jj12)
- self.jacobianInv[iQuad] = numpy.array([ [ij02[0,0], ij02[0,1]],
- [ij12[0,0], ij12[0,1]],
- [ij02[1,0], ij02[1,1]] ],
- dtype=numpy.Float64)
- else:
- self.jacobianInv[iQuad] = numpy.array([ [ij02[0,0], ij02[0,1]],
- [ 0.0, 0.0],
- [ij02[1,0], ij02[1,1]] ],
- dtype=numpy.Float64)
- elif abs(det12) > minJacobian:
- ij12 = numpy.linalg.inv(jj12)
- self.jacobianInv[iQuad] = numpy.array([ [ 0.0, 0.0],
- [ij12[0,0], ij12[0,1]],
- [ij12[1,0], ij12[1,1]] ],
- dtype=numpy.Float64)
- else:
- raise ValueError("Could not find inverse of Jacobian.")
- iQuad += 1
+ # PRIVATE METHODS ////////////////////////////////////////////////////
- # Quadrature points in cell
- self.quadPts = numpy.dot(self.basis, self.vertices)
+ def _configure(self):
+ """
+ Set members using inventory.
+ """
+ Script._configure(self)
+ self.data = self.inventory.data
return
-
+
def _initData(self):
self.data.addScalar(vtype="int", name="_numVertices",
value=self.numVertices,
@@ -234,12 +173,4 @@
return
- def _calculate(self):
- """
- Calculate basis functions, derivatives, and Jacobian information
- at quadrature points.
- """
- raise NotImplementedError
-
-
# End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXZ.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXZ.cc 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/QuadratureData2Din3DLinearXZ.cc 2006-10-28 00:15:25 UTC (rev 5110)
@@ -69,7 +69,7 @@
};
const double pylith::feassemble::QuadratureData2Din3DLinearXZ::_jacobianInv[] = {
- -1.00000000e+00, -0.00000000e+00,
+ -1.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 1.00000000e+00,
};
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/feutils.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/feutils.py 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/feutils.py 2006-10-28 00:15:25 UTC (rev 5110)
@@ -0,0 +1,118 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/feassemble/data/feutils.py
+
+## @brief Python routines for doing simple finite-element related
+## calculations.
+
+import numpy
+
+# ----------------------------------------------------------------------
+def calculateJacobian(quadrature, vertices):
+ """
+ Calculate jacobian, its determinant, and its inverse at quadrature
+ points for a given cell.
+
+ @param quadrature Quadrature information
+ @param vertices Coordinates of cell's vertices
+ """
+ jacobian = numpy.zeros( (quadrature.numQuadPts,
+ quadrature.cellDim, quadrature.spaceDim),
+ dtype=numpy.Float64)
+ jacobianInv = numpy.zeros( (quadrature.numQuadPts,
+ quadrature.spaceDim, quadrature.cellDim),
+ dtype=numpy.Float64)
+ jacobianDet = numpy.zeros( (quadrature.numQuadPts,), dtype=numpy.Float64)
+
+ iQuad = 0
+ for q in quadrature.quadPtsRef:
+ # Jacobian at quadrature points
+ deriv = quadrature.basisDeriv[iQuad]
+ j = numpy.dot(deriv.transpose(), vertices)
+ jacobian[iQuad] = j
+
+ # Determinant of Jacobian and Jacobian inverse at quadrature points
+ if quadrature.spaceDim == quadrature.cellDim:
+ jacobianDet[iQuad] = numpy.linalg.det(j)
+ jacobianInv[iQuad] = numpy.linalg.inv(j)
+ else:
+ det = numpy.linalg.det(numpy.dot(j, j.transpose()))**0.5
+ jacobianDet[iQuad] = det
+
+ if 1 == quadrature.cellDim:
+ jacobianInv[iQuad] = 1.0 / j
+ elif 2 == quadrature.cellDim:
+ minJacobian = 1.0e-06
+ jj01 = j[:,[0,1]]
+ jj12 = j[:,[1,2]]
+ jj02 = j[:,[0,2]]
+ det01 = numpy.linalg.det(jj01)
+ det12 = numpy.linalg.det(jj12)
+ det02 = numpy.linalg.det(jj02)
+ if abs(det01) > minJacobian:
+ ij01 = numpy.linalg.inv(jj01)
+ if abs(det12) > minJacobian:
+ ij12 = numpy.linalg.inv(jj12)
+ jacobianInv[iQuad] = numpy.array([ [ij01[0,0], ij01[0,1]],
+ [ij01[1,0], ij01[1,1]],
+ [ij12[1,0], ij12[1,1]] ],
+ dtype=numpy.Float64)
+ elif abs(det02) > minJacobian:
+ ij02 = numpy.linalg.inv(jj02)
+ jacobianInv[iQuad] = numpy.array([ [ij01[0,0], ij01[0,1]],
+ [ij01[1,0], ij01[1,1]],
+ [ij02[1,0], ij02[1,1]] ],
+ dtype=numpy.Float64)
+ else:
+ jacobianInv[iQuad] = numpy.array([ [ij01[0,0], ij01[0,1]],
+ [ij01[1,0], ij01[1,1]],
+ [ 0.0, 0.0] ],
+ dtype=numpy.Float64)
+ elif abs(det02) > minJacobian:
+ ij02 = numpy.linalg.inv(jj02)
+ if abs(det12) > minJacobian:
+ ij12 = numpy.linalg.inv(jj12)
+ jacobianInv[iQuad] = numpy.array([ [ij02[0,0], ij02[0,1]],
+ [ij12[0,0], ij12[0,1]],
+ [ij02[1,0], ij02[1,1]] ],
+ dtype=numpy.Float64)
+ else:
+ jacobianInv[iQuad] = numpy.array([ [ij02[0,0], ij02[0,1]],
+ [ 0.0, 0.0],
+ [ij02[1,0], ij02[1,1]] ],
+ dtype=numpy.Float64)
+ elif abs(det12) > minJacobian:
+ ij12 = numpy.linalg.inv(jj12)
+ jacobianInv[iQuad] = numpy.array([ [ 0.0, 0.0],
+ [ij12[0,0], ij12[0,1]],
+ [ij12[1,0], ij12[1,1]] ],
+ dtype=numpy.Float64)
+ else:
+ raise ValueError("Could not find inverse of Jacobian.")
+ else:
+ raise ValueError("Could not find inverse of Jacobian.")
+ iQuad += 1
+ return (jacobian, jacobianInv, jacobianDet)
+
+
+# ----------------------------------------------------------------------
+def assembleMat(globalMat, cellMat, cell):
+ """
+ Assemble cell matrix into global matrix.
+ """
+ # :KLUDGE: Assume only 1 cell in problem so assembly is trivial
+ globalMat += cellMat
+ return
+
+
+# End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/generate.sh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/generate.sh 2006-10-27 21:26:23 UTC (rev 5109)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/generate.sh 2006-10-28 00:15:25 UTC (rev 5110)
@@ -10,85 +10,113 @@
# ======================================================================
#
-# 1-D ------------------------------------------------------------------
+if (( $# != 1 )); then
+ echo "usage: generate.sh quadrature|integrator|all"
+ exit 1
+fi
-python Quadrature1DLinear.py \
+# //////////////////////////////////////////////////////////////////////
+if [ $1 == "quadrature" ] || [ $1 == "all" ]; then
+
+ # 1-D ----------------------------------------------------------------
+
+ python Quadrature1DLinear.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData1DLinear \
--data.parent=QuadratureData
-python Quadrature1DQuadratic.py \
+ python Quadrature1DQuadratic.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData1DQuadratic \
--data.parent=QuadratureData
-python Quadrature1Din2DLinear.py \
+ python Quadrature1Din2DLinear.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData1Din2DLinear \
--data.parent=QuadratureData
-python Quadrature1Din2DQuadratic.py \
+ python Quadrature1Din2DQuadratic.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData1Din2DQuadratic \
--data.parent=QuadratureData
-python Quadrature1Din3DLinear.py \
+ python Quadrature1Din3DLinear.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData1Din3DLinear \
--data.parent=QuadratureData
-python Quadrature1Din3DQuadratic.py \
+ python Quadrature1Din3DQuadratic.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData1Din3DQuadratic \
--data.parent=QuadratureData
-# 2-D ------------------------------------------------------------------
+ # 2-D ----------------------------------------------------------------
-python Quadrature2DLinear.py \
+ python Quadrature2DLinear.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData2DLinear \
--data.parent=QuadratureData
-python Quadrature2DQuadratic.py \
+ python Quadrature2DQuadratic.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData2DQuadratic \
--data.parent=QuadratureData
-python Quadrature2Din3DLinearXYZ.py \
+ python Quadrature2Din3DLinearXYZ.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData2Din3DLinearXYZ \
--data.parent=QuadratureData
-python Quadrature2Din3DLinearXY.py \
+ python Quadrature2Din3DLinearXY.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData2Din3DLinearXY \
--data.parent=QuadratureData
-python Quadrature2Din3DLinearYZ.py \
+ python Quadrature2Din3DLinearYZ.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData2Din3DLinearYZ \
--data.parent=QuadratureData
-python Quadrature2Din3DLinearXZ.py \
+ python Quadrature2Din3DLinearXZ.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData2Din3DLinearXZ \
--data.parent=QuadratureData
-python Quadrature2Din3DQuadratic.py \
+ python Quadrature2Din3DQuadratic.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData2Din3DQuadratic \
--data.parent=QuadratureData
-# 3-D ------------------------------------------------------------------
+ # 3-D ----------------------------------------------------------------
-python Quadrature3DLinear.py \
+ python Quadrature3DLinear.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData3DLinear \
--data.parent=QuadratureData
-python Quadrature3DQuadratic.py \
+ python Quadrature3DQuadratic.py \
--data.namespace=pylith,feassemble \
--data.object=QuadratureData3DQuadratic \
--data.parent=QuadratureData
+fi
+
+# //////////////////////////////////////////////////////////////////////
+if [ $1 == "integrator" ] || [ $1 == "all" ]; then
+
+ # 1-D ----------------------------------------------------------------
+
+ python IntegratorInertia1DLinear.py \
+ --data.namespace=pylith,feassemble \
+ --data.object=IntegratorDataInertia1DLinear \
+ --data.parent=IntegratorData
+
+ python IntegratorInertia1DQuadratic.py \
+ --data.namespace=pylith,feassemble \
+ --data.object=IntegratorDataInertia1DQuadratic \
+ --data.parent=IntegratorData
+
+fi
+
+
# End of file
More information about the cig-commits
mailing list