[cig-commits] r9281 - in short/3D/PyLith/trunk/libsrc: feassemble materials
brad at geodynamics.org
brad at geodynamics.org
Tue Feb 12 17:04:44 PST 2008
Author: brad
Date: 2008-02-12 17:04:43 -0800 (Tue, 12 Feb 2008)
New Revision: 9281
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
short/3D/PyLith/trunk/libsrc/materials/Material.cc
short/3D/PyLith/trunk/libsrc/materials/Material.hh
Log:
Implemented getting fields from elasticity integrator.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2008-02-12 22:58:47 UTC (rev 9280)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2008-02-13 01:04:43 UTC (rev 9281)
@@ -24,8 +24,6 @@
#include <assert.h> // USES assert()
#include <stdexcept> // USES std::runtime_error
-#define FASTER
-
// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::IntegratorElasticity::IntegratorElasticity(void) :
@@ -121,9 +119,7 @@
totalStrain = 0.0;
const ALE::Obj<real_section_type>& disp = fields->getSolution();
-#ifdef FASTER
const int dispAtlasTag = fields->getSolutionAtlasTag(materialId);
-#endif
// Loop over cells
int c_index = 0;
@@ -134,11 +130,7 @@
_quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
// Restrict input fields to cell
-#ifdef FASTER
mesh->restrict(disp, dispAtlasTag, c_index, &dispCell[0], cellVecSize);
-#else
- mesh->restrict(disp, *c_iter, &dispCell[0], cellVecSize);
-#endif
// Get cell geometry information that depends on cell
const double_array& basisDeriv = _quadrature->basisDeriv();
@@ -209,6 +201,205 @@
} // verifyConfiguration
// ----------------------------------------------------------------------
+// Get cell field associated with integrator.
+const ALE::Obj<pylith::real_section_type>&
+pylith::feassemble::IntegratorElasticity::cellField(
+ VectorFieldEnum* fieldType,
+ const char* name,
+ const ALE::Obj<Mesh>& mesh,
+ topology::FieldsManager* const fields)
+{ // cellField
+ assert(0 != _material);
+
+ // We assume the material stores the total-strain field if
+ // usesUpdateProperties() is TRUE.
+
+ if (!_material->usesUpdateProperties() &&
+ (0 == strcasecmp(name, "total-strain") ||
+ 0 == strcasecmp(name, "stress")) ) {
+
+ _calcStrainStressField(&_bufferCellTensor, name, mesh, fields);
+ return _bufferCellTensor;
+ } else if (0 == strcasecmp(name, "stress")) {
+ int fiberDim = 0;
+ const ALE::Obj<real_section_type>& strain =
+ _material->propertyField(&fiberDim, fieldType, "total-strain");
+ _calcStressFromStrain(&_bufferCellTensor, mesh, strain);
+ return _bufferCellTensor;
+ } else {
+ int fiberDim = 0;
+ const ALE::Obj<real_section_type>& field =
+ _material->propertyField(&fiberDim, fieldType, name);
+ switch (*fieldType)
+ { // switch
+ case SCALAR_FIELD :
+ _bufferCellScalar = field;
+ return _bufferCellScalar;
+ break;
+ case TENSOR_FIELD :
+ _bufferCellTensor = field;
+ return _bufferCellTensor;
+ break;
+ case OTHER_FIELD :
+ default:
+ _bufferCellOther = field;
+ return _bufferCellOther;
+ } // switch
+ } // else
+
+ return _bufferCellScalar;
+} // cellField
+
+// ----------------------------------------------------------------------
+void
+pylith::feassemble::IntegratorElasticity::_calcStrainStressField(
+ ALE::Obj<real_section_type>* field,
+ const char* name,
+ const ALE::Obj<Mesh>& mesh,
+ topology::FieldsManager* const fields)
+{ // _calcStrainStressField
+ assert(0 != _quadrature);
+ assert(0 != _material);
+
+ bool calcStress = (0 == strcasecmp(name, "stress")) ? true : false;
+
+ const int cellDim = _quadrature->cellDim();
+ int tensorSize = 0;
+ totalStrain_fn_type calcTotalStrainFn;
+ if (1 == cellDim) {
+ tensorSize = 1;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D;
+ } else if (2 == cellDim) {
+ tensorSize = 3;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D;
+ } else if (3 == cellDim) {
+ tensorSize = 6;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain3D;
+ } else
+ assert(0);
+
+ // Get cell information
+ const int materialId = _material->id();
+ const ALE::Obj<Mesh::label_sequence>& cells =
+ mesh->getLabelStratum("material-id", materialId);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+
+ const int cellVecSize = numBasis*spaceDim;
+ double_array dispCell(cellVecSize);
+
+ // Allocate vector for total strain
+ double_array totalStrain(numQuadPts*tensorSize);
+ totalStrain = 0.0;
+
+ // Allocate buffer for tensor field.
+ if (field->isNull()) {
+ const int fiberDim = numQuadPts * tensorSize;
+ *field = new real_section_type(mesh->comm(), mesh->debug());
+ (*field)->setFiberDimension(cells, fiberDim);
+ mesh->allocate(*field);
+ } // if
+
+ const ALE::Obj<real_section_type>& disp = fields->getSolution();
+ const int dispAtlasTag = fields->getSolutionAtlasTag(materialId);
+
+ // Loop over cells
+ int c_index = 0;
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter, ++c_index) {
+ // Compute geometry information for current cell
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
+
+ // Restrict input fields to cell
+ mesh->restrict(disp, dispAtlasTag, c_index, &dispCell[0], cellVecSize);
+
+ // Get cell geometry information that depends on cell
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ // Compute strains
+ calcTotalStrainFn(&totalStrain, basisDeriv, dispCell, numBasis,
+ numQuadPts);
+
+ if (!calcStress) {
+ (*field)->updatePoint(*c_iter, &totalStrain[0]);
+ } else {
+ const double_array& stress = _material->calcStress(totalStrain);
+ (*field)->updatePoint(*c_iter, &stress[0]);
+ } // else
+ } // for
+} // _calcStrainStressField
+
+// ----------------------------------------------------------------------
+void
+pylith::feassemble::IntegratorElasticity::_calcStressFromStrain(
+ ALE::Obj<real_section_type>* field,
+ const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<real_section_type>& strain)
+{ // _calcStressFromStrain
+ assert(0 != _quadrature);
+ assert(0 != _material);
+
+ const int cellDim = _quadrature->cellDim();
+ int tensorSize = 0;
+ if (1 == cellDim) {
+ tensorSize = 1;
+ } else if (2 == cellDim) {
+ tensorSize = 3;
+ } else if (3 == cellDim) {
+ tensorSize = 6;
+ } else
+ assert(0);
+
+ // Get cell information
+ const int materialId = _material->id();
+ const ALE::Obj<Mesh::label_sequence>& cells =
+ mesh->getLabelStratum("material-id", materialId);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+
+ // Allocate vector for total strain
+ double_array totalStrain(numQuadPts*tensorSize);
+ totalStrain = 0.0;
+
+ // Allocate buffer for tensor field.
+ if (field->isNull()) {
+ const int fiberDim = numQuadPts * tensorSize;
+ *field = new real_section_type(mesh->comm(), mesh->debug());
+ (*field)->setFiberDimension(cells, fiberDim);
+ mesh->allocate(*field);
+ } // if
+
+ // Loop over cells
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ const real_section_type::value_type* strainVals =
+ strain->restrictPoint(*c_iter);
+ memcpy(&totalStrain[0], strainVals, sizeof(double)*totalStrain.size());
+
+ const double_array& stress = _material->calcStress(totalStrain);
+ (*field)->updatePoint(*c_iter, &stress[0]);
+ } // for
+} // _calcStressFromStrain
+
+// ----------------------------------------------------------------------
// Integrate elasticity term in residual for 1-D cells.
void
pylith::feassemble::IntegratorElasticity::_elasticityResidual1D(
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2008-02-12 22:58:47 UTC (rev 9280)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2008-02-13 01:04:43 UTC (rev 9281)
@@ -22,6 +22,7 @@
#include "Integrator.hh" // ISA Integrator
#include "pylith/utils/array.hh" // USES std::vector, double_array
+#include "pylith/utils/vectorfields.hh" // USES VectorFieldEnum
namespace pylith {
namespace feassemble {
@@ -85,9 +86,47 @@
*/
void verifyConfiguration(const ALE::Obj<Mesh>& mesh);
+ /** Get cell field associated with integrator.
+ *
+ * @param fieldType Type of field.
+ * @param name Name of vertex field.
+ * @param mesh PETSc mesh for problem.
+ * @param fields Fields manager.
+ * @returns Cell field.
+ */
+ const ALE::Obj<real_section_type>&
+ cellField(VectorFieldEnum* fieldType,
+ const char* name,
+ const ALE::Obj<Mesh>& mesh,
+ topology::FieldsManager* const fields);
+
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
+ /** Calculate stress or strain field from solution field.
+ *
+ * @param field Field in which to store stress or strain.
+ * @param name Name of field to compute ('total-strain' or 'stress')
+ * @param mesh PETSc mesh for problem.
+ * @param fields Fields manager with solution.
+ */
+ void _calcStrainStressField(ALE::Obj<real_section_type>* field,
+ const char* name,
+ const ALE::Obj<Mesh>& mesh,
+ topology::FieldsManager* const fields);
+
+ /** Calculate stress field from total strain field.
+ *
+ * @param field Field in which to store stress.
+ * @param mesh PETSc mesh for problem.
+ * @param strain Total strain field.
+ */
+ void _calcStressFromStrain(ALE::Obj<real_section_type>* field,
+ const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<real_section_type>& strain);
+
+
/** Integrate elasticity term in residual for 1-D cells.
*
* @param stress Stress tensor for cell at quadrature points.
@@ -185,6 +224,15 @@
/// Elastic material associated with integrator
materials::ElasticMaterial* _material;
+ /// Scalar field for output of cell information.
+ ALE::Obj<real_section_type> _bufferCellScalar;
+
+ /// Tensor field for output of cell information.
+ ALE::Obj<real_section_type> _bufferCellTensor;
+
+ /// Other field for output of cell information.
+ ALE::Obj<real_section_type> _bufferCellOther;
+
}; // IntegratorElasticity
#endif // pylith_feassemble_integratorelasticity_hh
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.cc 2008-02-12 22:58:47 UTC (rev 9280)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.cc 2008-02-13 01:04:43 UTC (rev 9281)
@@ -76,9 +76,12 @@
assert(!cells.isNull());
const ALE::Mesh::label_sequence::iterator cellsEnd = cells->end();
- // Create sections to hold parameters for physical properties
+ // Create sections to hold physical properties and state variables.
_properties = new real_section_type(mesh->comm(), mesh->debug());
assert(!_properties.isNull());
+ for (int i=0; i < _numProperties; ++i)
+ _properties->addSpace();
+
const int numQuadPts = quadrature->numQuadPts();
const int spaceDim = quadrature->spaceDim();
@@ -87,6 +90,9 @@
const int totalPropsQuadPt = _totalPropsQuadPt;
const int fiberDim = totalPropsQuadPt * numQuadPts;
_properties->setFiberDimension(cells, fiberDim);
+ for (int i=0; i < _numProperties; ++i)
+ _properties->setFiberDimension(cells, spaceDim,
+ _propMetaData[i].fiberDim*numQuadPts);
mesh->allocate(_properties);
// Setup database for querying
@@ -137,20 +143,17 @@
} // initialize
// ----------------------------------------------------------------------
-// Get metadata for physical property. Values are returned through the
-// arguments.
-void
-pylith::materials::Material::propertyInfo(int* space,
- int* fiberDim,
- VectorFieldEnum* fieldType,
- const char* name) const
-{ // propertyInfo
+// Get physical property field.
+ALE::Obj<pylith::real_section_type>
+pylith::materials::Material::propertyField(int* fiberDim,
+ VectorFieldEnum* fieldType,
+ const char* name) const
+{ // propertyField
int i=0;
while (i < _numProperties)
if (0 == strcasecmp(name, _propMetaData[i].name))
break;
if (i < _numProperties) {
- *space = i;
*fiberDim = _propMetaData[i].fiberDim;
*fieldType = _propMetaData[i].fieldType;
} else {
@@ -159,7 +162,10 @@
<< _label << "'.";
throw std::runtime_error(msg.str());
} // else
-} // propertyInfo
+
+ assert(i < _numProperties);
+ return _properties->getFibration(i);
+} // propertyField
// End of file
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.hh 2008-02-12 22:58:47 UTC (rev 9280)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.hh 2008-02-13 01:04:43 UTC (rev 9281)
@@ -154,18 +154,16 @@
/// current state.
void resetNeedNewJacobian(void);
- /** Get metadata for physical property. Values are returned through
- * the arguments.
+ /** Get physical property field. Meta data is returned via the
+ * arguments.
*
- * @param space Subspace index.
* @param fiberDim Fiber dimension.
* @param fieldType Vector field type.
* @param name Name of physical property.
*/
- void propertyInfo(int* space,
- int* fiberDim,
- VectorFieldEnum* fieldType,
- const char* name) const;
+ ALE::Obj<real_section_type> propertyField(int* fiberDim,
+ VectorFieldEnum* fieldType,
+ const char* name) const;
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
More information about the cig-commits
mailing list