[cig-commits] r7000 - in short/3D/PyLith/trunk: libsrc/feassemble
unittests/libtests/feassemble unittests/libtests/feassemble/data
brad at geodynamics.org
brad at geodynamics.org
Wed May 30 10:10:57 PDT 2007
Author: brad
Date: 2007-05-30 10:10:56 -0700 (Wed, 30 May 2007)
New Revision: 7000
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.hh
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.hh
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit1DLinear.hh
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicit.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData1DLinear.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorApp.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorElasticity.py
Log:
Finished implementing initial C++ unit tests for ElasticityExplicit. Fixed a few bugs in ElasticityExplicit.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.cc 2007-05-30 10:33:02 UTC (rev 6999)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.cc 2007-05-30 17:10:56 UTC (rev 7000)
@@ -22,16 +22,16 @@
pylith::feassemble::Elasticity::calcTotalStrain1D(
std::vector<double_array>* strain,
const double_array& basisDeriv,
- const double* disp,
+ const double_array& disp,
const int numBasis)
{ // calcTotalStrain1D
assert(0 != strain);
- assert(0 != disp);
const int dimension = 1;
const int numQuadPts = strain->size();
assert(basisDeriv.size() == numQuadPts*numBasis*dimension);
+ assert(disp.size() == numBasis*dimension);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
assert(1 == (*strain)[iQuad].size());
@@ -46,16 +46,16 @@
pylith::feassemble::Elasticity::calcTotalStrain2D(
std::vector<double_array>* strain,
const double_array& basisDeriv,
- const double* disp,
+ const double_array& disp,
const int numBasis)
{ // calcTotalStrain2D
assert(0 != strain);
- assert(0 != disp);
const int dimension = 2;
const int numQuadPts = strain->size();
assert(basisDeriv.size() == numQuadPts*numBasis*dimension);
+ assert(disp.size() == numBasis*dimension);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
assert(3 == (*strain)[iQuad].size());
@@ -76,16 +76,16 @@
pylith::feassemble::Elasticity::calcTotalStrain3D(
std::vector<double_array>* strain,
const double_array& basisDeriv,
- const double* disp,
+ const double_array& disp,
const int numBasis)
{ // calcTotalStrain3D
assert(0 != strain);
- assert(0 != disp);
const int dimension = 3;
const int numQuadPts = strain->size();
assert(basisDeriv.size() == numQuadPts*numBasis*dimension);
+ assert(disp.size() == numBasis*dimension);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
assert(6 == (*strain)[iQuad].size());
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.hh 2007-05-30 10:33:02 UTC (rev 6999)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.hh 2007-05-30 17:10:56 UTC (rev 7000)
@@ -48,7 +48,7 @@
static
void calcTotalStrain1D(std::vector<double_array>* strain,
const double_array& basisDeriv,
- const double* disp,
+ const double_array& disp,
const int numBasis);
/** Compute total strain in at quadrature points of a cell.
@@ -62,7 +62,7 @@
static
void calcTotalStrain2D(std::vector<double_array>* strain,
const double_array& basisDeriv,
- const double* disp,
+ const double_array& disp,
const int numBasis);
/** Compute total strain in at quadrature points of a cell.
@@ -76,7 +76,7 @@
static
void calcTotalStrain3D(std::vector<double_array>* strain,
const double_array& basisDeriv,
- const double* disp,
+ const double_array& disp,
const int numBasis);
}; // Elasticity
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2007-05-30 10:33:02 UTC (rev 6999)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2007-05-30 17:10:56 UTC (rev 7000)
@@ -150,12 +150,10 @@
_resetCellVector();
// Restrict input fields to cell
- const real_section_type::value_type* dispTpdtCell =
- mesh->restrict(dispTpdt, *c_iter);
- const real_section_type::value_type* dispTCell =
- mesh->restrict(dispT, *c_iter);
- const real_section_type::value_type* dispTmdtCell =
- mesh->restrict(dispTmdt, *c_iter);
+ const int vecSize = spaceDim*numBasis;
+ double_array dispTpdtCell(mesh->restrict(dispTpdt, *c_iter), vecSize);
+ double_array dispTCell(mesh->restrict(dispT, *c_iter), vecSize);
+ double_array dispTmdtCell(mesh->restrict(dispTmdt, *c_iter), vecSize);
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -184,12 +182,12 @@
PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(6*spaceDim))));
if (err)
throw std::runtime_error("Logging PETSc flops failed.");
-
+
// Compute action for elastic terms
if (1 == cellDim) {
// Compute stresses
Elasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
+ dispTCell, numBasis);
const std::vector<double_array>& stress =
_material->calcStress(totalStrain);
@@ -210,7 +208,7 @@
} else if (2 == cellDim) {
// Compute stresses
Elasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
+ dispTCell, numBasis);
const std::vector<double_array>& stress =
_material->calcStress(totalStrain);
@@ -276,13 +274,13 @@
// Compute matrix associated with operator.
void
pylith::feassemble::ElasticityExplicit::integrateJacobian(
- PetscMat* mat,
+ PetscMat* jacobian,
topology::FieldsManager* fields,
const ALE::Obj<Mesh>& mesh)
{ // integrateJacobian
assert(0 != _quadrature);
assert(0 != _material);
- assert(0 != mat);
+ assert(0 != jacobian);
assert(0 != fields);
assert(!mesh.isNull());
@@ -356,7 +354,7 @@
const ALE::Obj<Mesh::order_type>& globalOrder =
mesh->getFactory()->getGlobalOrder(mesh, "default", dispT);
- err = updateOperator(*mat, mesh, dispT, globalOrder,
+ err = updateOperator(*jacobian, mesh, dispT, globalOrder,
*c_iter, _cellMatrix, ADD_VALUES);
if (err)
throw std::runtime_error("Update to PETSc Mat failed.");
@@ -421,8 +419,8 @@
_material->initCellData(*c_iter, numQuadPts);
// Restrict input fields to cell
- const real_section_type::value_type* dispCell =
- mesh->restrict(disp, *c_iter);
+ const int vecSize = spaceDim*numBasis;
+ double_array dispCell(mesh->restrict(disp, *c_iter), vecSize);
// Get cell geometry information that depends on cell
const double_array& basisDeriv = _quadrature->basisDeriv();
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh 2007-05-30 10:33:02 UTC (rev 6999)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh 2007-05-30 17:10:56 UTC (rev 7000)
@@ -115,11 +115,11 @@
/** Integrate contributions to Jacobian matrix (A) associated with
* operator.
*
- * @param mat Sparse matrix
- * @param fields Solution fields
- * @param mesh Finite-element mesh
+ * @param jacobian Sparse matrix to hold Jacobian of operator.
+ * @param fields Solution fields.
+ * @param mesh Finite-element mesh.
*/
- void integrateJacobian(PetscMat* mat,
+ void integrateJacobian(PetscMat* jacobian,
topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh);
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-05-30 10:33:02 UTC (rev 6999)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-05-30 17:10:56 UTC (rev 7000)
@@ -130,8 +130,8 @@
_resetCellVector();
// Restrict input fields to cell
- const real_section_type::value_type* dispTCell =
- mesh->restrict(dispT, *c_iter);
+ const int vecSize = spaceDim*numBasis;
+ double_array dispTCell(mesh->restrict(dispT, *c_iter), vecSize);
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -192,7 +192,7 @@
} else if (2 == cellDim) {
// Compute total strains and then use these to compute stresses
Elasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
+ dispTCell, numBasis);
const std::vector<double_array>& stress =
_material->calcStress(totalStrain);
@@ -217,7 +217,7 @@
} else if (3 == cellDim) {
// Compute total strains and then use these to compute stresses
Elasticity::calcTotalStrain3D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
+ dispTCell, numBasis);
const std::vector<double_array>& stress =
_material->calcStress(totalStrain);
@@ -331,8 +331,8 @@
_resetCellMatrix();
// Restrict input fields to cell
- const real_section_type::value_type* dispTCell =
- mesh->restrict(dispT, *c_iter);
+ const int vecSize = spaceDim*numBasis;
+ double_array dispTCell(mesh->restrict(dispT, *c_iter), vecSize);
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -569,8 +569,8 @@
_material->initCellData(*c_iter, numQuadPts);
// Restrict input fields to cell
- const real_section_type::value_type* dispCell =
- mesh->restrict(disp, *c_iter);
+ const int vecSize = spaceDim*numBasis;
+ double_array dispCell(mesh->restrict(disp, *c_iter), vecSize);
// Get cell geometry information that depends on cell
const double_array& basisDeriv = _quadrature->basisDeriv();
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2007-05-30 10:33:02 UTC (rev 6999)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2007-05-30 17:10:56 UTC (rev 7000)
@@ -25,6 +25,8 @@
#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
+#include <math.h> // USES fabs()
+
#include <stdexcept> // TEMPORARY
// ----------------------------------------------------------------------
@@ -88,7 +90,16 @@
void
pylith::feassemble::TestElasticityExplicit::testUpdateState(void)
{ // testUpdateState
- throw std::logic_error("Unit test not implemented.");
+ CPPUNIT_ASSERT(0 != _data);
+
+ ALE::Obj<ALE::Mesh> mesh;
+ ElasticityExplicit integrator;
+ topology::FieldsManager fields(mesh);
+ _initialize(&mesh, &integrator, &fields);
+
+ const ALE::Obj<real_section_type>& dispT = fields.getReal("dispT");
+ CPPUNIT_ASSERT(!dispT.isNull());
+ integrator.updateState(dispT, mesh);
} // testUpdateState
// ----------------------------------------------------------------------
@@ -116,7 +127,7 @@
const double tolerance = 1.0e-06;
for (int i=0; i < size; ++i)
- if (valsE[i] > 1.0)
+ if (fabs(valsE[i]) > 1.0)
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
@@ -127,10 +138,63 @@
void
pylith::feassemble::TestElasticityExplicit::testIntegrateJacobian(void)
{ // testIntegrateJacobian
- throw std::logic_error("Unit test not implemented.");
+ CPPUNIT_ASSERT(0 != _data);
+
+ ALE::Obj<ALE::Mesh> mesh;
+ ElasticityExplicit integrator;
+ topology::FieldsManager fields(mesh);
+ _initialize(&mesh, &integrator, &fields);
+
+ const ALE::Obj<pylith::real_section_type>& dispTpdt =
+ fields.getReal("dispTpdt");
+ CPPUNIT_ASSERT(!dispTpdt.isNull());
+
+ PetscMat jacobian;
+ PetscErrorCode err = MeshCreateMatrix(mesh, dispTpdt, MATMPIBAIJ, &jacobian);
+ CPPUNIT_ASSERT(0 == err);
+
+ integrator.integrateJacobian(&jacobian, &fields, mesh);
+ err = MatAssemblyBegin(jacobian, MAT_FINAL_ASSEMBLY);
+ CPPUNIT_ASSERT(0 == err);
+ err = MatAssemblyEnd(jacobian, MAT_FINAL_ASSEMBLY);
+ CPPUNIT_ASSERT(0 == err);
+
+ const double* valsE = _data->valsJacobian;
+ const int nrowsE = _data->numVertices * _data->spaceDim;
+ const int ncolsE = _data->numVertices * _data->spaceDim;
+
+ int nrows = 0;
+ int ncols = 0;
+ MatGetSize(jacobian, &nrows, &ncols);
+ CPPUNIT_ASSERT_EQUAL(nrowsE, nrows);
+ CPPUNIT_ASSERT_EQUAL(ncolsE, ncols);
+
+ PetscMat jDense;
+ PetscMat jSparseAIJ;
+ MatConvert(jacobian, MATSEQAIJ, MAT_INITIAL_MATRIX, &jSparseAIJ);
+ MatConvert(jSparseAIJ, MATSEQDENSE, MAT_INITIAL_MATRIX, &jDense);
+
+ double_array vals(nrows*ncols);
+ int_array rows(nrows);
+ int_array cols(ncols);
+ for (int iRow=0; iRow < nrows; ++iRow)
+ rows[iRow] = iRow;
+ for (int iCol=0; iCol < ncols; ++iCol)
+ cols[iCol] = iCol;
+ MatGetValues(jDense, nrows, &rows[0], ncols, &cols[0], &vals[0]);
+ const double tolerance = 1.0e-06;
+ for (int iRow=0; iRow < nrows; ++iRow)
+ for (int iCol=0; iCol < ncols; ++iCol) {
+ const int index = ncols*iRow+iCol;
+ if (fabs(valsE[index]) > 1.0)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[index]/valsE[index], tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[index], vals[index], tolerance);
+ } // for
+ MatDestroy(jDense);
+ MatDestroy(jSparseAIJ);
} // testIntegrateJacobian
-
// ----------------------------------------------------------------------
// Initialize elasticity integrator.
void
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.hh 2007-05-30 10:33:02 UTC (rev 6999)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.hh 2007-05-30 17:10:56 UTC (rev 7000)
@@ -55,7 +55,6 @@
CPPUNIT_TEST( testTimeStep );
CPPUNIT_TEST( testStableTimeStep );
CPPUNIT_TEST( testMaterial );
- CPPUNIT_TEST( testUpdateState );
CPPUNIT_TEST_SUITE_END();
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit1DLinear.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit1DLinear.hh 2007-05-30 10:33:02 UTC (rev 6999)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit1DLinear.hh 2007-05-30 17:10:56 UTC (rev 7000)
@@ -38,6 +38,7 @@
// CPPUNIT TEST SUITE /////////////////////////////////////////////////
CPPUNIT_TEST_SUITE( TestElasticityExplicit1DLinear );
+ CPPUNIT_TEST( testUpdateState );
CPPUNIT_TEST( testIntegrateResidual );
CPPUNIT_TEST( testIntegrateJacobian );
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicit.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicit.py 2007-05-30 10:33:02 UTC (rev 6999)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicit.py 2007-05-30 17:10:56 UTC (rev 7000)
@@ -48,9 +48,11 @@
[K]{u(t)}
"""
K = self._calculateStiffnessMat()
+ print K
+
M = self._calculateMassMat()
+ print M
-
dispResult = -self.fieldTpdt + 2.0*self.fieldT - self.fieldTmdt
self.valsResidual = 1.0/self.dt**2 * numpy.dot(M, dispResult) - \
numpy.dot(K, self.fieldT)
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData1DLinear.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData1DLinear.cc 2007-05-30 10:33:02 UTC (rev 6999)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData1DLinear.cc 2007-05-30 17:10:56 UTC (rev 7000)
@@ -47,7 +47,7 @@
};
const double pylith::feassemble::ElasticityExplicitData1DLinear::_quadPts[] = {
- 8.75000000e-01,
+ 0.00000000e+00,
};
const double pylith::feassemble::ElasticityExplicitData1DLinear::_quadWts[] = {
@@ -70,18 +70,18 @@
};
const double pylith::feassemble::ElasticityExplicitData1DLinear::_fieldT[] = {
- 1.20000000e+00,
- 1.70000000e+00,
+ 1.10000000e+00,
+ 1.50000000e+00,
};
const double pylith::feassemble::ElasticityExplicitData1DLinear::_fieldTmdt[] = {
- 1.20000000e+00,
- 1.70000000e+00,
+ 1.00000000e+00,
+ 1.30000000e+00,
};
const double pylith::feassemble::ElasticityExplicitData1DLinear::_valsResidual[] = {
- -5.56875000e+10,
- 3.12250226e-09,
+ 2.02500000e+10,
+ -2.02500000e+10,
};
const double pylith::feassemble::ElasticityExplicitData1DLinear::_valsJacobian[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorApp.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorApp.py 2007-05-30 10:33:02 UTC (rev 6999)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorApp.py 2007-05-30 17:10:56 UTC (rev 7000)
@@ -150,7 +150,6 @@
self.quadPts = self.quadrature.quadPtsRef
self.quadWts = self.quadrature.quadWts
(self.basis, self.basisDeriv) = self.quadrature.calculateBasis()
- self.quadPts = numpy.dot(self.basis, self.vertices)
# Material information
self.matType = self.material.type
@@ -233,10 +232,10 @@
values=self.fieldTpdt,
format="%16.8e", ncols=self.spaceDim)
self.data.addArray(vtype="double", name="_fieldT",
- values=self.fieldTpdt,
+ values=self.fieldT,
format="%16.8e", ncols=self.spaceDim)
self.data.addArray(vtype="double", name="_fieldTmdt",
- values=self.fieldTpdt,
+ values=self.fieldTmdt,
format="%16.8e", ncols=self.spaceDim)
# Calculated values
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorElasticity.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorElasticity.py 2007-05-30 10:33:02 UTC (rev 6999)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorElasticity.py 2007-05-30 17:10:56 UTC (rev 7000)
@@ -59,6 +59,7 @@
# Matrix of elasticity values
D = self._calculateElasticityMat()
+ print "D: ",D
for cell in self.cells:
cellK = numpy.zeros( (self.spaceDim*self.numBasis,
@@ -70,6 +71,7 @@
for iQuad in xrange(self.numQuadPts):
wt = self.quadWts[iQuad] * jacobianDet[iQuad]
B = self._calculateBasisDerivMat(iQuad)
+ print "B: ",B
cellK[:] += wt * numpy.dot(numpy.dot(B.transpose(), D), B)
feutils.assembleMat(K, cellK, cell, self.spaceDim)
return K
@@ -94,6 +96,7 @@
for iQuad in xrange(self.numQuadPts):
wt = self.quadWts[iQuad] * jacobianDet[iQuad]
N = self._calculateBasisMat(iQuad)
+ print "N: ",N
cellM[:] += self.density * wt * numpy.dot(N.transpose(), N)
feutils.assembleMat(M, cellM, cell, self.spaceDim)
return M
@@ -191,10 +194,10 @@
B[2, iBasis*self.spaceDim+0] = self.basisDeriv[iQuad, iBasis, 1]
B[2, iBasis*self.spaceDim+1] = self.basisDeriv[iQuad, iBasis, 0]
elif 1 == self.spaceDim:
- iBasis = 0
B = numpy.zeros( (1, self.spaceDim*self.numBasis),
dtype=numpy.float64)
- B[0, iBasis*self.spaceDim+0] = self.basisDeriv[iQuad, iBasis, 0]
+ for iBasis in xrange(self.numBasis):
+ B[0, iBasis*self.spaceDim+0] = self.basisDeriv[iQuad, iBasis, 0]
else:
raise ValueError("Unknown spatial dimension '%d'." % self.spaceDim)
return B
More information about the cig-commits
mailing list