[cig-commits] r15041 - in short/3D/PyLith/trunk: . libsrc/feassemble libsrc/meshio tests/1d/line2
brad at geodynamics.org
brad at geodynamics.org
Sun May 24 20:41:47 PDT 2009
Author: brad
Date: 2009-05-24 20:41:47 -0700 (Sun, 24 May 2009)
New Revision: 15041
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc
short/3D/PyLith/trunk/tests/1d/line2/TestDislocation.py
Log:
Fixed some nondimensionalization bugs.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2009-05-25 02:40:17 UTC (rev 15040)
+++ short/3D/PyLith/trunk/TODO 2009-05-25 03:41:47 UTC (rev 15041)
@@ -9,6 +9,8 @@
Brad
Test precompute geometry [DONE]
Nondimensionalization
+ Fault coordinates - need coordinates_dimensioned
+ Material state vars - use dimstatevars and dimproperties to get scales
cleanup
full-scale testing
test uniform refinement
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2009-05-25 02:40:17 UTC (rev 15040)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2009-05-25 03:41:47 UTC (rev 15041)
@@ -314,6 +314,7 @@
topology::SolutionFields* fields)
{ // cellField
assert(0 != _material);
+ assert(0 != _normalizer);
// We assume the material stores the total_strain field if
// hasStateVars() is TRUE.
@@ -322,16 +323,28 @@
_outputFields =
new topology::Fields<topology::Field<topology::Mesh> >(mesh);
- if (!_material->hasStateVars() &&
- (0 == strcasecmp(name, "total_strain") ||
- 0 == strcasecmp(name, "stress") )) {
+ if (!_material->hasStateVars() && 0 == strcasecmp(name, "total_strain")) {
assert(0 != fields);
_allocateTensorField(mesh);
topology::Field<topology::Mesh>& buffer =
_outputFields->get("buffer (tensor)");
_calcStrainStressField(&buffer, name, fields);
- buffer.label(name);
+ buffer.label("total_strain");
+ buffer.scale(1.0);
+ buffer.addDimensionOkay(true);
return buffer;
+
+ } else if (!_material->hasStateVars() && 0 == strcasecmp(name, "stress")) {
+ assert(0 != fields);
+ _allocateTensorField(mesh);
+ topology::Field<topology::Mesh>& buffer =
+ _outputFields->get("buffer (tensor)");
+ _calcStrainStressField(&buffer, name, fields);
+ buffer.label("stress");
+ buffer.scale(_normalizer->pressureScale());
+ buffer.addDimensionOkay(true);
+ return buffer;
+
} else if (0 == strcasecmp(name, "stress")) {
assert(0 != fields);
_allocateTensorField(mesh);
@@ -340,13 +353,17 @@
_material->getField(&buffer, "total_strain");
_calcStressFromStrain(&buffer);
buffer.label(name);
+ buffer.scale(_normalizer->pressureScale());
+ buffer.addDimensionOkay(true);
return buffer;
+
} else {
if (!_outputFields->hasField("buffer (other)"))
_outputFields->add("buffer (other)", "buffer");
topology::Field<topology::Mesh>& buffer =
_outputFields->get("buffer (other)");
_material->getField(&buffer, name);
+ buffer.addDimensionOkay(true);
return buffer;
} // if/else
Modified: short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc 2009-05-25 02:40:17 UTC (rev 15040)
+++ short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc 2009-05-25 03:41:47 UTC (rev 15041)
@@ -14,6 +14,8 @@
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+#include "pylith/topology/Field.hh" // USES Field
+
// ----------------------------------------------------------------------
// Constructor
template<typename mesh_type, typename field_type>
@@ -82,7 +84,8 @@
// Only processors with cells for output get the correct fiber dimension.
const ALE::Obj<RealSection>& sectionIn = fieldIn.section();
assert(!sectionIn.isNull());
- const int totalFiberDim = sectionIn->getFiberDimension(*cellsBegin);
+ const int totalFiberDim = (cellsBegin != cellsEnd) ?
+ sectionIn->getFiberDimension(*cellsBegin) : 0;
const int fiberDim = totalFiberDim / numQuadPts;
assert(fiberDim * numQuadPts == totalFiberDim);
@@ -92,35 +95,45 @@
assert(0 != _fieldAvg);
_fieldAvg->newSection(sectionIn->getChart(), fiberDim);
_fieldAvg->allocate();
+ } else if (_fieldAvg->size() != cells->size()*fiberDim) {
+ _fieldAvg->clear();
+ _fieldAvg->newSection(sectionIn->getChart(), fiberDim);
+ _fieldAvg->allocate();
+ } // else
- switch (fieldIn.vectorFieldType())
- { // switch
- case topology::FieldBase::MULTI_SCALAR:
- _fieldAvg->vectorFieldType(topology::FieldBase::SCALAR);
- break;
- case topology::FieldBase::MULTI_VECTOR:
- _fieldAvg->vectorFieldType(topology::FieldBase::VECTOR);
- break;
- case topology::FieldBase::MULTI_TENSOR:
- _fieldAvg->vectorFieldType(topology::FieldBase::TENSOR);
- break;
- case topology::FieldBase::MULTI_OTHER:
- _fieldAvg->vectorFieldType(topology::FieldBase::OTHER);
- break;
- case topology::FieldBase::SCALAR:
- case topology::FieldBase::VECTOR:
- case topology::FieldBase::TENSOR:
- case topology::FieldBase::OTHER:
- default :
- std::cerr << "Bad vector field type for CellFilterAvg." << std::endl;
- assert(0);
- throw std::logic_error("Bad vector field type for CellFilterAvg.");
- } // switch
- } // if
+ std::cout << "FIELD " << fieldIn.label()
+ << ", fiberDim: " << totalFiberDim << "\n"
+ << "AVG FIELD: " << _fieldAvg->label()
+ << ", fiberDim: " << fiberDim << std::endl;
+
assert(0 != _fieldAvg);
+ switch (fieldIn.vectorFieldType())
+ { // switch
+ case topology::FieldBase::MULTI_SCALAR:
+ _fieldAvg->vectorFieldType(topology::FieldBase::SCALAR);
+ break;
+ case topology::FieldBase::MULTI_VECTOR:
+ _fieldAvg->vectorFieldType(topology::FieldBase::VECTOR);
+ break;
+ case topology::FieldBase::MULTI_TENSOR:
+ _fieldAvg->vectorFieldType(topology::FieldBase::TENSOR);
+ break;
+ case topology::FieldBase::MULTI_OTHER:
+ _fieldAvg->vectorFieldType(topology::FieldBase::OTHER);
+ break;
+ case topology::FieldBase::SCALAR:
+ case topology::FieldBase::VECTOR:
+ case topology::FieldBase::TENSOR:
+ case topology::FieldBase::OTHER:
+ default :
+ std::cerr << "Bad vector field type for CellFilterAvg." << std::endl;
+ assert(0);
+ throw std::logic_error("Bad vector field type for CellFilterAvg.");
+ } // switch
const ALE::Obj<RealSection>& sectionAvg = _fieldAvg->section();
_fieldAvg->label(fieldIn.label());
_fieldAvg->scale(fieldIn.scale());
+ _fieldAvg->addDimensionOkay(true);
double_array fieldAvgCell(fiberDim);
double scalar = 0.0;
@@ -142,8 +155,6 @@
} // for
PetscLogFlops( cells->size() * numQuadPts*fiberDim*3 );
- _fieldAvg->label(fieldIn.label());
-
return *_fieldAvg;
} // filter
Modified: short/3D/PyLith/trunk/tests/1d/line2/TestDislocation.py
===================================================================
--- short/3D/PyLith/trunk/tests/1d/line2/TestDislocation.py 2009-05-25 02:40:17 UTC (rev 15040)
+++ short/3D/PyLith/trunk/tests/1d/line2/TestDislocation.py 2009-05-25 03:41:47 UTC (rev 15041)
@@ -111,7 +111,6 @@
maskN*(-0.20 - 0.025*vertices[:,0]) + \
maskP*(+0.30 - 0.025*vertices[:,0])
- print "VERTICES",vertices
return disp
More information about the CIG-COMMITS
mailing list