[cig-commits] r14152 - in short/3D/PyLith/branches/pylith-swig: libsrc/materials unittests/libtests/materials unittests/libtests/materials/data
brad at geodynamics.org
brad at geodynamics.org
Wed Feb 25 17:06:43 PST 2009
Author: brad
Date: 2009-02-25 17:06:43 -0800 (Wed, 25 Feb 2009)
New Revision: 14152
Modified:
short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.cc
short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.hh
short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/TestElasticMaterial.cc
short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/TestElasticMaterial.hh
short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/TestMaterial.cc
short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/data/Makefile.am
short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/data/matinitialize.spatialdb
Log:
Updated more ElasticMaterial libtests.
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.cc 2009-02-25 19:54:52 UTC (rev 14151)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.cc 2009-02-26 01:06:43 UTC (rev 14152)
@@ -91,11 +91,30 @@
assert(!propertiesSection.isNull());
propertiesSection->restrictPoint(cell, &_propertiesCell[0],
_propertiesCell.size());
-
- const ALE::Obj<RealSection>& stateVarsSection = _stateVars->section();
- assert(!stateVarsSection.isNull());
- stateVarsSection->restrictPoint(cell, &_stateVarsCell[0],
- _stateVarsCell.size());
+
+ if (hasStateVars()) {
+ const ALE::Obj<RealSection>& stateVarsSection = _stateVars->section();
+ assert(!stateVarsSection.isNull());
+ stateVarsSection->restrictPoint(cell, &_stateVarsCell[0],
+ _stateVarsCell.size());
+ } // if
+
+ if (0 == _initialStress)
+ _initialStressCell = 0.0;
+ else {
+ const ALE::Obj<RealSection>& stressSection = _initialStress->section();
+ assert(!stressSection.isNull());
+ stressSection->restrictPoint(cell, &_initialStressCell[0],
+ _initialStressCell.size());
+ } // if/else
+ if (0 == _initialStrain)
+ _initialStrainCell = 0.0;
+ else {
+ const ALE::Obj<RealSection>& strainSection = _initialStrain->section();
+ assert(!strainSection.isNull());
+ strainSection->restrictPoint(cell, &_initialStrainCell[0],
+ _initialStrainCell.size());
+ } // if/else
} // retrievePropsAndVars
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.hh 2009-02-25 19:54:52 UTC (rev 14151)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/materials/ElasticMaterial.hh 2009-02-26 01:06:43 UTC (rev 14152)
@@ -77,6 +77,9 @@
/** Compute density for cell at quadrature points.
*
+ * @pre Must call retrievePropsAndVars for cell before calling
+ * calcDensity().
+ *
* @returns Array of density values at cell's quadrature points.
*/
const double_array& calcDensity(void);
@@ -86,6 +89,9 @@
* should be set to true so that the state variables are updated
* (but not stored) when computing the stresses.
*
+ * @pre Must call retrievePropsAndVars for cell before calling
+ * calcStress().
+ *
* Size of array of stress tensors = [numQuadPts][tensorSize].
*
* Order of stresses for 3-D:
@@ -109,6 +115,9 @@
/** Compute derivative of elasticity matrix for cell at quadrature points.
*
+ * @pre Must call retrievePropsAndVars for cell before calling
+ * calcDerivElastic().
+ *
* Size of array of elasticity constants = [numQuadPts][numElasticConsts]
*
* Order of elasticity constants for 3-D:
@@ -151,6 +160,9 @@
/** Get stable time step for implicit time integration.
*
+ * @pre Must call retrievePropsAndVars for cell before calling
+ * stableTimeStep().
+ *
* Default is MAXFLOAT (or 1.0e+30 if MAXFLOAT is not defined in math.h).
*
* @returns Time step
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/TestElasticMaterial.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/TestElasticMaterial.cc 2009-02-25 19:54:52 UTC (rev 14151)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/TestElasticMaterial.cc 2009-02-26 01:06:43 UTC (rev 14152)
@@ -15,9 +15,16 @@
#include "TestElasticMaterial.hh" // Implementation of class methods
#include "data/ElasticMaterialData.hh" // USES ElasticMaterialData
-//#include "data/ElasticIsotropic3DData.hh" // USES ElasticIsotropic3DData
+#include "data/ElasticStrain1DData.hh" // USES ElasticStrain1DData
+#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
#include "pylith/materials/ElasticIsotropic3D.hh" // USES ElasticIsotropic3D
+#include "pylith/materials/ElasticStrain1D.hh" // USES ElasticStrain1D
+#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+#include "pylith/feassemble/GeometryLine1D.hh" // USES GeometryLine1D
+
#include "pylith/utils/array.hh" // USES double_array
#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
@@ -31,6 +38,10 @@
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::materials::TestElasticMaterial );
// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
// Test dbInitialStress()
void
pylith::materials::TestElasticMaterial::testDBInitialStress(void)
@@ -62,45 +73,195 @@
CPPUNIT_ASSERT_EQUAL(label, std::string(material._dbInitialStrain->label()));
} // testDBInitialStrain
-#if 0
// ----------------------------------------------------------------------
// Test initialize()
void
pylith::materials::TestElasticMaterial::testInitialize(void)
{ // testInitialize
topology::Mesh mesh;
- Quadrature<topology::Mesh> quadrature;
ElasticStrain1D material;
ElasticStrain1DData data;
- _initialize(&mesh, &quadrature, &material, &data);
+ _initialize(&mesh, &material, &data);
- CPPUNIT_ASSERT(false);
+ // Get cells associated with material
+ const int materialId = 24;
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+ SieveMesh::point_type cell = *cells->begin();
+
+ const double tolerance = 1.0e-06;
+ const int tensorSize = material._tensorSize;
+ const int numQuadPts = data.numLocs;
+
+ // Test initialStress field
+ CPPUNIT_ASSERT(0 != material._initialStress);
+ const ALE::Obj<RealSection>& stressSection =
+ material._initialStress->section();
+ CPPUNIT_ASSERT(!stressSection.isNull());
+ int fiberDim = numQuadPts * tensorSize;
+ CPPUNIT_ASSERT_EQUAL(fiberDim, stressSection->getFiberDimension(cell));
+ const double* initialStress = stressSection->restrictPoint(cell);
+ CPPUNIT_ASSERT(0 != initialStress);
+ const double* initialStressE = data.initialStress;
+ CPPUNIT_ASSERT(0 != initialStressE);
+ for (int i=0; i < fiberDim; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, initialStress[i]/initialStressE[i],
+ tolerance);
+
+ // Test initialStrain field
+ CPPUNIT_ASSERT(0 != material._initialStrain);
+ const ALE::Obj<RealSection>& strainSection =
+ material._initialStrain->section();
+ CPPUNIT_ASSERT(!strainSection.isNull());
+ fiberDim = numQuadPts * tensorSize;
+ CPPUNIT_ASSERT_EQUAL(fiberDim, strainSection->getFiberDimension(cell));
+ const double* initialStrain = strainSection->restrictPoint(cell);
+ CPPUNIT_ASSERT(0 != initialStrain);
+ const double* initialStrainE = data.initialStrain;
+ CPPUNIT_ASSERT(0 != initialStrainE);
+ for (int i=0; i < fiberDim; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, initialStrain[i]/initialStrainE[i],
+ tolerance);
+
+ // Test cell arrays
+ size_t size = data.numLocs*data.numPropsQuadPt;
+ CPPUNIT_ASSERT_EQUAL(size, material._propertiesCell.size());
+
+ size = data.numLocs*data.numVarsQuadPt;
+ CPPUNIT_ASSERT_EQUAL(size, material._stateVarsCell.size());
+
+ size = data.numLocs*tensorSize;
+ CPPUNIT_ASSERT_EQUAL(size, material._initialStressCell.size());
+
+ size = data.numLocs*tensorSize;
+ CPPUNIT_ASSERT_EQUAL(size, material._initialStrainCell.size());
+
+ size = data.numLocs;
+ CPPUNIT_ASSERT_EQUAL(size, material._densityCell.size());
+
+ size = data.numLocs*tensorSize;
+ CPPUNIT_ASSERT_EQUAL(size, material._stressCell.size());
+
+ int numElasticConsts = 0;
+ switch (data.dimension)
+ { // switch
+ case 1 :
+ numElasticConsts = 1;
+ break;
+ case 2 :
+ numElasticConsts = 6;
+ break;
+ case 3 :
+ numElasticConsts = 21;
+ break;
+ default :
+ assert(0);
+ } // switch
+ size = data.numLocs*numElasticConsts;
+ CPPUNIT_ASSERT_EQUAL(size, material._elasticConstsCell.size());
} // testInitialize
// ----------------------------------------------------------------------
+// Test retrievePropsAndVars().
+void
+pylith::materials::TestElasticMaterial::testRetrievePropsAndVars(void)
+{ // testRetrievePropsAndVars
+ topology::Mesh mesh;
+ ElasticStrain1D material;
+ ElasticStrain1DData data;
+ _initialize(&mesh, &material, &data);
+
+ // Get cells associated with material
+ const int materialId = 24;
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+ SieveMesh::point_type cell = *cells->begin();
+
+ material.retrievePropsAndVars(cell);
+
+ const double tolerance = 1.0e-06;
+ const int tensorSize = material._tensorSize;
+ const int numQuadPts = data.numLocs;
+ const int numVarsQuadPt = data.numVarsQuadPt;
+
+ // Test cell arrays
+ const double* propertiesE = data.properties;
+ CPPUNIT_ASSERT(0 != propertiesE);
+ const double_array& properties = material._propertiesCell;
+ size_t size = data.numLocs*data.numPropsQuadPt;
+ CPPUNIT_ASSERT_EQUAL(size, properties.size());
+ for (size_t i=0; i < size; ++i)
+ {
+ std::cout << "propertyE: " << propertiesE[i]
+ << ", property: " << properties[i]
+ << std::endl;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, properties[i]/propertiesE[i],
+ tolerance);
+ }
+
+ const double* stateVarsE = data.stateVars;
+ CPPUNIT_ASSERT( (0 < numVarsQuadPt && 0 != stateVarsE) ||
+ (0 == numVarsQuadPt && 0 == stateVarsE) );
+ const double_array& stateVars = material._stateVarsCell;
+ size = data.numLocs*numVarsQuadPt;
+ CPPUNIT_ASSERT_EQUAL(size, stateVars.size());
+ for (size_t i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, stateVars[i]/stateVarsE[i],
+ tolerance);
+
+ const double* initialStressE = data.initialStress;
+ CPPUNIT_ASSERT(0 != initialStressE);
+ const double_array& initialStress = material._initialStressCell;
+ size = data.numLocs*tensorSize;
+ CPPUNIT_ASSERT_EQUAL(size, initialStress.size());
+ for (size_t i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, initialStress[i]/initialStressE[i],
+ tolerance);
+
+ const double* initialStrainE = data.initialStrain;
+ CPPUNIT_ASSERT(0 != initialStrainE);
+ const double_array& initialStrain = material._initialStrainCell;
+ size = data.numLocs*tensorSize;
+ CPPUNIT_ASSERT_EQUAL(size, initialStrain.size());
+ for (size_t i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, initialStrain[i]/initialStrainE[i],
+ tolerance);
+} // testRetrievePropsAndVars
+
+// ----------------------------------------------------------------------
// Test calcDensity()
void
pylith::materials::TestElasticMaterial::testCalcDensity(void)
{ // testCalcDensity
topology::Mesh mesh;
- Quadrature<topology::Mesh> quadrature;
ElasticStrain1D material;
ElasticStrain1DData data;
- _initialize(&mesh, &quadrature, &material, &data);
+ _initialize(&mesh, &material, &data);
// Get cells associated with material
- const ALE::Obj<SieveMesh::label_sequence>& cells = mesh->heightStratum(0);
- SieveMesh::label_sequence::iterator c_iter = cells->begin();
- material.retrievePropsAndVars(*c_iter);
+ const int materialId = 24;
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+ SieveMesh::point_type cell = *cells->begin();
+
+ material.retrievePropsAndVars(cell);
const double_array& density = material.calcDensity();
- CPPUNIT_ASSERT_EQUAL(data._numQuadPts, density.size());
+ const int tensorSize = material._tensorSize;
+ const int numQuadPts = data.numLocs;
- const double tolerance = 1.0e-06;
const double* densityE = data.density;
CPPUNIT_ASSERT(0 != densityE);
- const double size = data._numQuadPts;
- for (int i=0; i < size; ++i)
+ const size_t size = numQuadPts;
+ CPPUNIT_ASSERT_EQUAL(size, density.size());
+ const double tolerance = 1.0e-06;
+ for (size_t i=0; i < size; ++i)
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, density[i]/densityE[i], tolerance);
} // testCalcDensity
@@ -108,80 +269,36 @@
// Test calcStress()
void
pylith::materials::TestElasticMaterial::testCalcStress(void)
-{ // testCalcProperties
- ALE::Obj<Mesh> mesh;
- { // create mesh
- const int cellDim = 1;
- const int numCorners = 2;
- const int spaceDim = 1;
- const int numVertices = 2;
- const int numCells = 1;
- const double vertCoords[] = { -1.0, 1.0};
- const int cells[] = { 0, 1};
- CPPUNIT_ASSERT(0 != vertCoords);
- CPPUNIT_ASSERT(0 != cells);
+{ // testCalcStress
+ topology::Mesh mesh;
+ ElasticStrain1D material;
+ ElasticStrain1DData data;
+ _initialize(&mesh, &material, &data);
- mesh = new Mesh(PETSC_COMM_WORLD, cellDim);
- ALE::Obj<sieve_type> sieve = new sieve_type(mesh->comm());
-
- const bool interpolate = false;
- ALE::Obj<ALE::Mesh::sieve_type> s = new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
-
- ALE::SieveBuilder<ALE::Mesh>::buildTopology(s, cellDim, numCells,
- const_cast<int*>(cells), numVertices, interpolate, numCorners);
- std::map<Mesh::point_type,Mesh::point_type> renumbering;
- ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
- mesh->setSieve(sieve);
- mesh->stratify();
- ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, spaceDim, vertCoords);
- } // create mesh
-
// Get cells associated with material
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ const int materialId = 24;
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+ SieveMesh::point_type cell = *cells->begin();
- ElasticIsotropic3D material;
- ElasticIsotropic3DData data;
- const int numQuadPts = 2;
- const int numParams = data.numParameters;
- const int numParamsQuadPt = data.numParamsQuadPt;
+ const int tensorSize = material._tensorSize;
+ const int numQuadPts = data.numLocs;
- Mesh::label_sequence::iterator c_iter = cells->begin();
+ // Setup total strain
+ double_array strain(data.strain, numQuadPts*tensorSize);
- const int fiberDim = numQuadPts * numParamsQuadPt;
- material._properties = new real_section_type(mesh->comm(), mesh->debug());
- material._properties->setChart(mesh->getSieve()->getChart());
- material._properties->setFiberDimension(cells, fiberDim);
- mesh->allocate(material._properties);
-
- material._properties->updatePoint(*c_iter, data.parameterData);
- const int strainSize = material._tensorSize;
- double_array strain(data.strain, numQuadPts*strainSize);
-
- const int initialStateSize = material._tensorSize;
- material._initialStateSize = initialStateSize;
- const int initialStateFiberDim = numQuadPts * initialStateSize;
- material._initialState = new real_section_type(mesh->comm(), mesh->debug());
- material._initialState->setChart(mesh->getSieve()->getChart());
- material._initialState->setFiberDimension(cells, initialStateFiberDim);
- mesh->allocate(material._initialState);
- material._initialState->updatePoint(*c_iter, data.initialState);
-
- material.getPropertiesCell(*c_iter, numQuadPts);
+ material.retrievePropsAndVars(cell);
const double_array& stress = material.calcStress(strain);
- const double tolerance = 1.0e-06;
const double* stressE = data.stress;
CPPUNIT_ASSERT(0 != stressE);
-
- const int size = numQuadPts * strainSize;;
- CPPUNIT_ASSERT_EQUAL(size, int(stress.size()));
- for (int i=0; i < size; ++i)
- if (fabs(stressE[i]) > tolerance)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, stress[i]/stressE[i], tolerance);
- else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(stressE[i], stress[i], tolerance);
+ const size_t size = numQuadPts * tensorSize;
+ CPPUNIT_ASSERT_EQUAL(size, stress.size());
+ const double tolerance = 1.0e-06;
+ for (size_t i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, stress[i]/stressE[i], tolerance);
} // testCalcStress
// ----------------------------------------------------------------------
@@ -189,81 +306,52 @@
void
pylith::materials::TestElasticMaterial::testCalcDerivElastic(void)
{ // testCalcDerivElastic
- ALE::Obj<Mesh> mesh;
- { // create mesh
- const int cellDim = 1;
- const int numCorners = 2;
- const int spaceDim = 1;
- const int numVertices = 2;
- const int numCells = 1;
- const double vertCoords[] = { -1.0, 1.0};
- const int cells[] = { 0, 1};
- CPPUNIT_ASSERT(0 != vertCoords);
- CPPUNIT_ASSERT(0 != cells);
+ topology::Mesh mesh;
+ ElasticStrain1D material;
+ ElasticStrain1DData data;
+ _initialize(&mesh, &material, &data);
- mesh = new Mesh(PETSC_COMM_WORLD, cellDim);
- ALE::Obj<sieve_type> sieve = new sieve_type(mesh->comm());
-
- const bool interpolate = false;
- ALE::Obj<ALE::Mesh::sieve_type> s = new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
-
- ALE::SieveBuilder<ALE::Mesh>::buildTopology(s, cellDim, numCells,
- const_cast<int*>(cells), numVertices, interpolate, numCorners);
- std::map<Mesh::point_type,Mesh::point_type> renumbering;
- ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
- mesh->setSieve(sieve);
- mesh->stratify();
- ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, spaceDim, vertCoords);
- } // create mesh
-
// Get cells associated with material
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ const int materialId = 24;
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+ SieveMesh::point_type cell = *cells->begin();
- ElasticIsotropic3D material;
- ElasticIsotropic3DData data;
- const int numQuadPts = 2;
- const int numParams = data.numParameters;
- const int numParamsQuadPt = data.numParamsQuadPt;
+ const int tensorSize = material._tensorSize;
+ const int numQuadPts = data.numLocs;
- Mesh::label_sequence::iterator c_iter = cells->begin();
+ // Setup total strain
+ double_array strain(data.strain, numQuadPts*tensorSize);
- const int fiberDim = numQuadPts * numParamsQuadPt;
- material._properties = new real_section_type(mesh->comm(), mesh->debug());
- material._properties->setChart(mesh->getSieve()->getChart());
- material._properties->setFiberDimension(cells, fiberDim);
- mesh->allocate(material._properties);
+ material.retrievePropsAndVars(cell);
+ const double_array& elasticConsts = material.calcDerivElastic(strain);
- material._properties->updatePoint(*c_iter, data.parameterData);
- const int strainSize = material._tensorSize;
- double_array strain(data.strain, numQuadPts*strainSize);
+ int numElasticConsts = 0;
+ switch (data.dimension)
+ { // switch
+ case 1 :
+ numElasticConsts = 1;
+ break;
+ case 2 :
+ numElasticConsts = 6;
+ break;
+ case 3 :
+ numElasticConsts = 21;
+ break;
+ default :
+ assert(0);
+ } // switch
- const int initialStateSize = material._tensorSize;
- material._initialStateSize = initialStateSize;
- const int initialStateFiberDim = numQuadPts * initialStateSize;
- material._initialState = new real_section_type(mesh->comm(), mesh->debug());
- material._initialState->setChart(mesh->getSieve()->getChart());
- material._initialState->setFiberDimension(cells, initialStateFiberDim);
- mesh->allocate(material._initialState);
-
- material.getPropertiesCell(*c_iter, numQuadPts);
- const double_array& elasticConsts = material.calcDerivElastic(strain);
-
- const double tolerance = 1.0e-06;
const double* elasticConstsE = data.elasticConsts;
CPPUNIT_ASSERT(0 != elasticConstsE);
-
- const int size = elasticConsts.size();
- for (int i=0; i < size; ++i)
- if (fabs(elasticConstsE[i]) > tolerance)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0,
- elasticConsts[i]/elasticConstsE[i],
- tolerance);
- else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(elasticConstsE[i],
- elasticConsts[i],
- tolerance);
+ const size_t size = numQuadPts * numElasticConsts;
+ CPPUNIT_ASSERT_EQUAL(size, elasticConsts.size());
+ const double tolerance = 1.0e-06;
+ for (size_t i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, elasticConsts[i]/elasticConstsE[i],
+ tolerance);
} // testCalcDerivElastic
// ----------------------------------------------------------------------
@@ -279,71 +367,26 @@
void
pylith::materials::TestElasticMaterial::testStableTimeStepImplicit(void)
{ // testStableTimeStepImplicit
- ALE::Obj<Mesh> mesh;
- { // create mesh
- const int cellDim = 1;
- const int numCorners = 2;
- const int spaceDim = 1;
- const int numVertices = 2;
- const int numCells = 1;
- const double vertCoords[] = { -1.0, 1.0};
- const int cells[] = { 0, 1};
- CPPUNIT_ASSERT(0 != vertCoords);
- CPPUNIT_ASSERT(0 != cells);
+ topology::Mesh mesh;
+ ElasticStrain1D material;
+ ElasticStrain1DData data;
+ _initialize(&mesh, &material, &data);
- mesh = new Mesh(PETSC_COMM_WORLD, cellDim);
- ALE::Obj<sieve_type> sieve = new sieve_type(mesh->comm());
-
- const bool interpolate = false;
- ALE::Obj<ALE::Mesh::sieve_type> s =
- new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
-
- ALE::SieveBuilder<ALE::Mesh>::buildTopology(s, cellDim, numCells,
- const_cast<int*>(cells), numVertices, interpolate, numCorners);
- std::map<Mesh::point_type,Mesh::point_type> renumbering;
- ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
- mesh->setSieve(sieve);
- mesh->stratify();
- ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, spaceDim, vertCoords);
- } // create mesh
-
// Get cells associated with material
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ const int materialId = 24;
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+ SieveMesh::point_type cell = *cells->begin();
- ElasticIsotropic3D material;
- ElasticIsotropic3DData data;
- const int numQuadPts = 2;
- const int numParams = data.numParameters;
- const int numParamsQuadPt = data.numParamsQuadPt;
-
- const int initialStateSize = material._tensorSize;
- material._initialStateSize = initialStateSize;
- const int initialStateFiberDim = numQuadPts * initialStateSize;
- material._initialState = new real_section_type(mesh->comm(), mesh->debug());
- material._initialState->setChart(mesh->getSieve()->getChart());
- material._initialState->setFiberDimension(cells, initialStateFiberDim);
- mesh->allocate(material._initialState);
-
- Mesh::label_sequence::iterator c_iter = cells->begin();
-
- const int fiberDim = numQuadPts * numParamsQuadPt;
- material._properties = new real_section_type(mesh->comm(), mesh->debug());
- material._properties->setChart(mesh->getSieve()->getChart());
- material._properties->setFiberDimension(cells, fiberDim);
- mesh->allocate(material._properties);
-
- material._properties->updatePoint(*c_iter, data.parameterData);
-
- material.getPropertiesCell(*c_iter, numQuadPts);
+ material.retrievePropsAndVars(cell);
const double dt = material.stableTimeStepImplicit();
const double tolerance = 1.0e-06;
const double dtE = data.dtStableImplicit;
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, dt/dtE, tolerance);
} // testStableTimeStepImplicit
-#endif
// ----------------------------------------------------------------------
// Setup testing data.
@@ -615,5 +658,90 @@
_matElastic->normalizer(normalizer);
} // setupNormalizer
+// ----------------------------------------------------------------------
+// Setup mesh and material.
+void
+pylith::materials::TestElasticMaterial::_initialize(
+ topology::Mesh* mesh,
+ ElasticStrain1D* material,
+ const ElasticStrain1DData* data)
+{ // _initialize
+ CPPUNIT_ASSERT(0 != mesh);
+ CPPUNIT_ASSERT(0 != material);
+ CPPUNIT_ASSERT(0 != data);
+ meshio::MeshIOAscii iohandler;
+ iohandler.filename("data/line3.mesh");
+ iohandler.read(mesh);
+
+ // Set up coordinates
+ spatialdata::geocoords::CSCart cs;
+ cs.setSpaceDim(mesh->dimension());
+ cs.initialize();
+ mesh->coordsys(&cs);
+
+ // Setup quadrature
+ feassemble::Quadrature<topology::Mesh> quadrature;
+ feassemble::GeometryLine1D geometry;
+ quadrature.refGeometry(&geometry);
+ const int cellDim = 1;
+ const int numCorners = 3;
+ const int numQuadPts = 2;
+ const int spaceDim = 1;
+ const double basis[] = { 0.455, -0.122, 0.667, -0.122, 0.455, 0.667 };
+ const double basisDeriv[] = {
+ -1.07735027e+00,
+ -7.73502692e-02,
+ 1.15470054e+00,
+ 7.73502692e-02,
+ 1.07735027e+00,
+ -1.15470054e+00,
+ };
+ const double quadPtsRef[] = { -0.577350269, 0.577350269 };
+ const double quadWts[] = { 1.0, 1.0 };
+ quadrature.initialize(basis, basisDeriv, quadPtsRef, quadWts,
+ cellDim, numCorners, numQuadPts, spaceDim);
+
+
+ // Get cells associated with material
+ const int materialId = 24;
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+
+ // Compute geometry for cells
+ quadrature.computeGeometry(*mesh, cells);
+
+ spatialdata::spatialdb::SimpleDB db;
+ spatialdata::spatialdb::SimpleIOAscii dbIO;
+ dbIO.filename("data/matinitialize.spatialdb");
+ db.ioHandler(&dbIO);
+ db.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ spatialdata::spatialdb::SimpleDB dbStress;
+ spatialdata::spatialdb::SimpleIOAscii dbIOStress;
+ dbIOStress.filename("data/matstress.spatialdb");
+ dbStress.ioHandler(&dbIOStress);
+ dbStress.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ spatialdata::spatialdb::SimpleDB dbStrain;
+ spatialdata::spatialdb::SimpleIOAscii dbIOStrain;
+ dbIOStrain.filename("data/matstrain.spatialdb");
+ dbStrain.ioHandler(&dbIOStrain);
+ dbStrain.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ spatialdata::units::Nondimensional normalizer;
+
+ material->dbProperties(&db);
+ material->id(materialId);
+ material->label("my_material");
+ material->normalizer(normalizer);
+ material->dbInitialStress(&dbStress);
+ material->dbInitialStrain(&dbStrain);
+
+ material->initialize(*mesh, &quadrature);
+} // _initialize
+
+
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/TestElasticMaterial.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/TestElasticMaterial.hh 2009-02-25 19:54:52 UTC (rev 14151)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/TestElasticMaterial.hh 2009-02-26 01:06:43 UTC (rev 14152)
@@ -23,12 +23,14 @@
#include "TestMaterial.hh"
#include "pylith/materials/materialsfwd.hh" // forward declarations
+#include "pylith/topology/topologyfwd.hh" // forward declarations
/// Namespace for pylith package
namespace pylith {
namespace materials {
class TestElasticMaterial;
class ElasticMaterialData;
+ class ElasticStrain1DData;
} // materials
} // pylith
@@ -41,15 +43,13 @@
CPPUNIT_TEST( testDBInitialStress );
CPPUNIT_TEST( testDBInitialStrain );
-#if 0
CPPUNIT_TEST( testInitialize );
- CPPUNIT_TETS( testRetrievePropsAndVars );
+ CPPUNIT_TEST( testRetrievePropsAndVars );
CPPUNIT_TEST( testCalcDensity );
CPPUNIT_TEST( testCalcStress );
CPPUNIT_TEST( testCalcDerivElastic );
CPPUNIT_TEST( testUpdateStateVars );
CPPUNIT_TEST( testStableTimeStepImplicit );
-#endif
CPPUNIT_TEST_SUITE_END();
@@ -126,7 +126,6 @@
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
-#if 0
/** Setup mesh and material.
*
* @param mesh Finite-element mesh.
@@ -134,9 +133,9 @@
* @param data Data with properties for elastic material.
*/
void _initialize(topology::Mesh* mesh,
- ElasticMaterial* material,
- const ElasticMaterialData* data);
-#endif
+ ElasticStrain1D* material,
+ const ElasticStrain1DData* data);
+
}; // class TestElasticMaterial
#endif // pylith_materials_testelasticmaterial_hh
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/TestMaterial.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/TestMaterial.cc 2009-02-25 19:54:52 UTC (rev 14151)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/TestMaterial.cc 2009-02-26 01:06:43 UTC (rev 14152)
@@ -33,6 +33,7 @@
#include <cstring> // USES strcmp()
#include <cassert> // USES assert()
+#include <cmath> // USES sqrt()
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::materials::TestMaterial );
@@ -206,27 +207,27 @@
material.normalizer(normalizer);
material.initialize(mesh, &quadrature);
- const double densityA = 2000.0;
- const double vsA = 100.0;
- const double vpA = 180.0;
+ const double densityA = 2500.0;
+ const double vsA = 3000.0;
+ const double vpA = vsA*sqrt(3.0);
const double muA = vsA*vsA*densityA;
const double lambdaA = vpA*vpA*densityA - 2.0*muA;
- const double densityB = 3000.0;
- const double vsB = 200.0;
- const double vpB = 400.0;
+ const double densityB = 2000.0;
+ const double vsB = 1200.0;
+ const double vpB = vsB*sqrt(3.0);
const double muB = vsB*vsB*densityB;
const double lambdaB = vpB*vpB*densityB - 2.0*muB;
const double densityE[] = { densityA, densityB };
const double muE[] = { muA, muB };
const double lambdaE[] = { lambdaA, lambdaB };
- SieveMesh::label_sequence::iterator c_iter = cells->begin();
+ SieveMesh::point_type cell = *cells->begin();
const double tolerance = 1.0e-06;
CPPUNIT_ASSERT(0 != material._properties);
const Obj<RealSection>& propertiesSection = material._properties->section();
CPPUNIT_ASSERT(!propertiesSection.isNull());
- const double* propertiesCell = propertiesSection->restrictPoint(*c_iter);
+ const double* propertiesCell = propertiesSection->restrictPoint(cell);
CPPUNIT_ASSERT(0 != propertiesCell);
const int p_density = 0;
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/data/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/data/Makefile.am 2009-02-25 19:54:52 UTC (rev 14151)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/data/Makefile.am 2009-02-26 01:06:43 UTC (rev 14152)
@@ -12,6 +12,8 @@
dist_noinst_DATA = \
matinitialize.spatialdb \
+ matstress.spatialdb \
+ matstrain.spatialdb \
line3.mesh
noinst_TMP =
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/data/matinitialize.spatialdb
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/data/matinitialize.spatialdb 2009-02-25 19:54:52 UTC (rev 14151)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/materials/data/matinitialize.spatialdb 2009-02-26 01:06:43 UTC (rev 14152)
@@ -11,5 +11,5 @@
space-dim = 1
}
}
--0.5 2000.0 100.0 180.0
-+0.5 3000.0 200.0 400.0
+-0.5 2500.0 3000.0 5196.15242
++0.5 2000.0 1200.0 2078.46097
More information about the CIG-COMMITS
mailing list