[cig-commits] r15038 - in short/3D/PyLith/trunk: libsrc/faults libsrc/meshio libsrc/topology unittests/libtests/meshio
brad at geodynamics.org
brad at geodynamics.org
Sun May 24 13:23:46 PDT 2009
Author: brad
Date: 2009-05-24 13:23:46 -0700 (Sun, 24 May 2009)
New Revision: 15038
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc
short/3D/PyLith/trunk/libsrc/meshio/OutputManager.cc
short/3D/PyLith/trunk/libsrc/meshio/OutputManager.hh
short/3D/PyLith/trunk/libsrc/meshio/VertexFilterVecNorm.cc
short/3D/PyLith/trunk/libsrc/topology/Field.cc
short/3D/PyLith/trunk/libsrc/topology/Field.hh
short/3D/PyLith/trunk/unittests/libtests/meshio/TestOutputManager.cc
Log:
Added dimensionalization of output and fixed some related bugs.
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2009-05-22 17:37:36 UTC (rev 15037)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2009-05-24 20:23:46 UTC (rev 15038)
@@ -154,7 +154,7 @@
_calcArea();
// Create empty tractions field for change in fault tractions.
- _fields->add("tractions", "tractions_change");
+ _fields->add("tractions", "traction_change");
topology::Field<topology::SubMesh>& tractions = _fields->get("tractions");
tractions.vectorFieldType(topology::FieldBase::VECTOR);
tractions.scale(_normalizer->pressureScale());
@@ -1166,7 +1166,7 @@
assert(0 != _normalizer);
tractions->scale(_normalizer->pressureScale());
- tractions->label("traction");
+ tractions->label("traction_change");
// Get vertices from mesh of domain.
const ALE::Obj<SieveMesh>& sieveMesh = dispT.mesh().sieveMesh();
@@ -1230,7 +1230,6 @@
// Allocate buffer for tractions field (if nec.).
const ALE::Obj<RealSection>& tractionsSection = tractions->section();
- tractions->label("tractions");
if (tractionsSection.isNull()) {
tractions->newSection(topology::FieldBase::VERTICES_FIELD, fiberDim);
tractions->allocate();
Modified: short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc 2009-05-22 17:37:36 UTC (rev 15037)
+++ short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc 2009-05-24 20:23:46 UTC (rev 15038)
@@ -91,7 +91,6 @@
_fieldAvg = new field_type(fieldIn.mesh());
assert(0 != _fieldAvg);
_fieldAvg->newSection(sectionIn->getChart(), fiberDim);
- _fieldAvg->label("average");
_fieldAvg->allocate();
switch (fieldIn.vectorFieldType())
@@ -120,6 +119,8 @@
} // if
assert(0 != _fieldAvg);
const ALE::Obj<RealSection>& sectionAvg = _fieldAvg->section();
+ _fieldAvg->label(fieldIn.label());
+ _fieldAvg->scale(fieldIn.scale());
double_array fieldAvgCell(fiberDim);
double scalar = 0.0;
Modified: short/3D/PyLith/trunk/libsrc/meshio/OutputManager.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/OutputManager.cc 2009-05-22 17:37:36 UTC (rev 15037)
+++ short/3D/PyLith/trunk/libsrc/meshio/OutputManager.cc 2009-05-24 20:23:46 UTC (rev 15038)
@@ -16,6 +16,8 @@
#include "VertexFilter.hh" // USES VertexFilter
#include "CellFilter.hh" // USES CellFilter
+#include "pylith/topology/Fields.hh" // USES Fields
+
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
// ----------------------------------------------------------------------
@@ -25,7 +27,8 @@
_coordsys(0),
_writer(0),
_vertexFilter(0),
- _cellFilter(0)
+ _cellFilter(0),
+ _fields(0)
{ // constructor
} // constructor
@@ -38,6 +41,7 @@
_vertexFilter = 0; // :TODO: Use shared pointer
_cellFilter = 0; // :TODO: Use shared pointer
delete _coordsys; _coordsys = 0;
+ delete _fields; _fields = 0;
} // destructor
// ----------------------------------------------------------------------
@@ -140,8 +144,9 @@
{ // appendVertexField
const field_type& fieldFiltered =
(0 == _vertexFilter) ? field : _vertexFilter->filter(field);
-
- _writer->writeVertexField(t, fieldFiltered, mesh);
+ const field_type& fieldDimensioned = _dimension(fieldFiltered);
+
+ _writer->writeVertexField(t, fieldDimensioned, mesh);
} // appendVertexField
// ----------------------------------------------------------------------
@@ -156,9 +161,85 @@
{ // appendCellField
const field_type& fieldFiltered =
(0 == _cellFilter) ? field : _cellFilter->filter(field, label, labelId);
+ const field_type& fieldDimensioned = _dimension(fieldFiltered);
- _writer->writeCellField(t, fieldFiltered, label, labelId);
+ _writer->writeCellField(t, fieldDimensioned, label, labelId);
} // appendCellField
+// ----------------------------------------------------------------------
+// Dimension field.
+template<typename mesh_type, typename field_type>
+const field_type&
+pylith::meshio::OutputManager<mesh_type, field_type>::_dimension(const field_type& fieldIn)
+{ // _dimension
+ if (fieldIn.addDimensionOkay()) {
+ std::cout << "Adding dimensions to field '" << fieldIn.label()
+ << "'." << std::endl;
+ fieldIn.dimensionalize();
+
+ return fieldIn;
+ } else {
+ std::cout << "Creating temporary field to add dimensions to field '"
+ << fieldIn.label() << "'." << std::endl;
+
+ std::string fieldName = "buffer (other)";
+ switch (fieldIn.vectorFieldType())
+ { // switch
+ case topology::FieldBase::SCALAR :
+ fieldName = "buffer (scalar)";
+ break;
+ case topology::FieldBase::VECTOR :
+ fieldName = "buffer (vector)";
+ break;
+ case topology::FieldBase::TENSOR :
+ fieldName = "buffer (tensor)";
+ break;
+ case topology::FieldBase::OTHER :
+ fieldName = "buffer (other)";
+ break;
+ case topology::FieldBase::MULTI_SCALAR :
+ fieldName = "buffer (multiple scalars)";
+ break;
+ case topology::FieldBase::MULTI_VECTOR :
+ fieldName = "buffer (multiple vectors)";
+ break;
+ case topology::FieldBase::MULTI_TENSOR :
+ fieldName = "buffer (multiple tensors)";
+ break;
+ case topology::FieldBase::MULTI_OTHER :
+ fieldName = "buffer (multiple others)";
+ break;
+ default :
+ // Spit out usefule error message and stop via assert. If
+ // optimized, throw exception.
+ std::cerr << "Unknown field type '" << fieldIn.vectorFieldType()
+ << "'";
+ assert(0);
+ throw std::logic_error("Unknown field type");
+ } // switch
+
+ if (0 == _fields)
+ _fields = new topology::Fields<field_type>(fieldIn.mesh());
+
+ if (!_fields->hasField(fieldName.c_str())) {
+ _fields->add(fieldName.c_str(), fieldIn.label());
+ field_type& fieldOut = _fields->get(fieldName.c_str());
+ fieldOut.newSection(fieldIn);
+ fieldOut.allocate();
+ } // if
+ field_type& fieldOut = _fields->get(fieldName.c_str());
+ fieldOut.copy(fieldIn);
+ fieldOut.addDimensionOkay(true);
+ fieldOut.dimensionalize();
+
+ return fieldOut;
+ } // if/else
+
+ // Satisfy return value. Should never get this far.
+ assert(0 != _fields);
+ return _fields->get("buffer (dimensioned)");
+} // _dimension
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/meshio/OutputManager.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/OutputManager.hh 2009-05-22 17:37:36 UTC (rev 15037)
+++ short/3D/PyLith/trunk/libsrc/meshio/OutputManager.hh 2009-05-24 20:23:46 UTC (rev 15038)
@@ -121,6 +121,15 @@
const char* label =0,
const int labelId =0);
+// PROTECTED MEMBERS ////////////////////////////////////////////////////
+protected :
+
+ /** Dimension field.
+ *
+ * @param fieldIn Field to dimensionalize.
+ */
+ const field_type& _dimension(const field_type& fieldIn);
+
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
@@ -137,6 +146,8 @@
VertexFilter<field_type>* _vertexFilter; ///< Filter applied to vertex data.
CellFilter<mesh_type, field_type>* _cellFilter; ///< Filter applied to cell data.
+ topology::Fields<field_type>* _fields; ///< Buffer fields.
+
}; // OutputManager
#include "OutputManager.cc" // template methods
Modified: short/3D/PyLith/trunk/libsrc/meshio/VertexFilterVecNorm.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/VertexFilterVecNorm.cc 2009-05-22 17:37:36 UTC (rev 15037)
+++ short/3D/PyLith/trunk/libsrc/meshio/VertexFilterVecNorm.cc 2009-05-24 20:23:46 UTC (rev 15038)
@@ -79,7 +79,8 @@
_fieldVecNorm->newSection(sectionIn->getChart(), fiberDimNorm);
_fieldVecNorm->allocate();
- _fieldVecNorm->label(fieldIn.label());
+ _fieldVecNorm->label(fieldIn.label());
+ _fieldVecNorm->scale(fieldIn.scale());
switch (fieldIn.vectorFieldType())
{ // switch
case topology::FieldBase::SCALAR:
Modified: short/3D/PyLith/trunk/libsrc/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Field.cc 2009-05-22 17:37:36 UTC (rev 15037)
+++ short/3D/PyLith/trunk/libsrc/topology/Field.cc 2009-05-24 20:23:46 UTC (rev 15038)
@@ -305,7 +305,6 @@
const int dstSize = (!_section.isNull()) ? _section->size() : 0;
if (field.spaceDim() != spaceDim() ||
field._vecFieldType != _vecFieldType ||
- field._scale != _scale ||
srcSize != dstSize) {
std::ostringstream msg;
@@ -315,7 +314,7 @@
<< " space dim: " << field.spaceDim() << "\n"
<< " vector field type: " << field._vecFieldType << "\n"
<< " scale: " << field._scale << "\n"
- << " size: " << srcSize
+ << " size: " << srcSize << "\n"
<< " Destination section:\n"
<< " space dim: " << spaceDim() << "\n"
<< " vector field type: " << _vecFieldType << "\n"
@@ -340,6 +339,9 @@
_section->updatePointAll(*c_iter, field._section->restrictPoint(*c_iter));
} // for
} // if
+
+ _label = field._label;
+ _scale = field._scale;
} // copy
// ----------------------------------------------------------------------
@@ -357,7 +359,7 @@
msg << "Cannot copy values from Sieve section "
<< _label << "'. Sections are incompatible.\n"
<< " Source section:\n"
- << " size: " << srcSize
+ << " size: " << srcSize << "\n"
<< " Destination section:\n"
<< " space dim: " << spaceDim() << "\n"
<< " vector field type: " << _vecFieldType << "\n"
@@ -404,7 +406,7 @@
<< " space dim: " << field.spaceDim() << "\n"
<< " vector field type: " << field._vecFieldType << "\n"
<< " scale: " << field._scale << "\n"
- << " size: " << srcSize
+ << " size: " << srcSize << "\n"
<< " Destination section:\n"
<< " space dim: " << spaceDim() << "\n"
<< " vector field type: " << _vecFieldType << "\n"
@@ -443,7 +445,7 @@
// Dimensionalize field.
template<typename mesh_type>
void
-pylith::topology::Field<mesh_type>::dimensionalize(void)
+pylith::topology::Field<mesh_type>::dimensionalize(void) const
{ // dimensionalize
if (!_dimensionsOkay) {
std::ostringstream msg;
@@ -470,7 +472,7 @@
_section->restrictPoint(*c_iter, &values[0], values.size());
normalizer.dimensionalize(&values[0], values.size(), _scale);
- _section->updatePoint(*c_iter, &values[0]);
+ _section->updatePointAll(*c_iter, &values[0]);
} // for
} // if
} // dimensionalize
@@ -479,7 +481,7 @@
// Print field to standard out.
template<typename mesh_type>
void
-pylith::topology::Field<mesh_type>::view(const char* label)
+pylith::topology::Field<mesh_type>::view(const char* label) const
{ // view
std::string vecFieldString;
switch(_vecFieldType)
Modified: short/3D/PyLith/trunk/libsrc/topology/Field.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Field.hh 2009-05-22 17:37:36 UTC (rev 15037)
+++ short/3D/PyLith/trunk/libsrc/topology/Field.hh 2009-05-24 20:23:46 UTC (rev 15038)
@@ -216,13 +216,13 @@
/** Dimensionalize field. Throws runtime_error if field is not
* allowed to be dimensionalized.
*/
- void dimensionalize(void);
+ void dimensionalize(void) const;
/** Print field to standard out.
*
* @param label Label for output.
*/
- void view(const char* label);
+ void view(const char* label) const;
/// Create PETSc vector for field.
void createVector(void);
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestOutputManager.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestOutputManager.cc 2009-05-22 17:37:36 UTC (rev 15037)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestOutputManager.cc 2009-05-24 20:23:46 UTC (rev 15038)
@@ -179,6 +179,7 @@
3.1, 3.2,
4.1, 4.2
};
+ const double scale = 2.0;
topology::Mesh mesh;
MeshIOAscii iohandler;
@@ -198,17 +199,20 @@
field.allocate();
field.label(label);
field.vectorFieldType(fieldType);
+ field.scale(scale);
const ALE::Obj<RealSection>& section = field.section();
CPPUNIT_ASSERT(!section.isNull());
CPPUNIT_ASSERT_EQUAL(nvertices, int(vertices->size()));
+ double_array values(nvertices*fiberDim);
+ for (int i=0; i < nvertices*fiberDim; ++i)
+ values[i] = fieldValues[i];
+ values /= scale;
int ipt = 0;
for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
v_iter != verticesEnd;
- ++v_iter, ++ipt) {
- const double* values = &fieldValues[ipt*fiberDim];
- section->updatePoint(*v_iter, values);
- } // for
+ ++v_iter, ++ipt)
+ section->updatePoint(*v_iter, &values[ipt*fiberDim]);
spatialdata::geocoords::CSCart cs;
const int numTimeSteps = 1;
@@ -261,6 +265,7 @@
1.1, 1.2,
2.1, 2.2,
};
+ const double scale = 4.0;
topology::Mesh mesh;
MeshIOAscii iohandler;
@@ -280,17 +285,20 @@
field.allocate();
field.label(label);
field.vectorFieldType(fieldType);
+ field.scale(scale);
const ALE::Obj<RealSection>& section = field.section();
CPPUNIT_ASSERT(!section.isNull());
CPPUNIT_ASSERT_EQUAL(ncells, int(cells->size()));
+ double_array values(ncells*fiberDim);
+ for (int i=0; i < ncells*fiberDim; ++i)
+ values[i] = fieldValues[i];
+ values /= scale;
int ipt = 0;
for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
- ++c_iter, ++ipt) {
- const double* values = &fieldValues[ipt*fiberDim];
- section->updatePoint(*c_iter, values);
- } // for
+ ++c_iter, ++ipt)
+ section->updatePoint(*c_iter, &values[ipt*fiberDim]);
spatialdata::geocoords::CSCart cs;
const int numTimeSteps = 1;
More information about the CIG-COMMITS
mailing list