[cig-commits] r9146 - in short/3D/PyLith/trunk: libsrc/faults modulesrc/meshio pylith/faults pylith/meshio pylith/problems
brad at geodynamics.org
brad at geodynamics.org
Sun Jan 27 14:26:52 PST 2008
Author: brad
Date: 2008-01-27 14:26:52 -0800 (Sun, 27 Jan 2008)
New Revision: 9146
Modified:
short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/modulesrc/meshio/meshio.pyxe.src
short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py
short/3D/PyLith/trunk/pylith/meshio/OutputFaultKin.py
short/3D/PyLith/trunk/pylith/meshio/OutputManager.py
short/3D/PyLith/trunk/pylith/problems/Formulation.py
Log:
More work on output of fault information. Output of orientation directions now works.
Modified: short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc 2008-01-27 19:07:29 UTC (rev 9145)
+++ short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc 2008-01-27 22:26:52 UTC (rev 9146)
@@ -676,12 +676,26 @@
int color = 0;
const int coneSize = cone->size();
- const int faceSize = (constraintCell) ? coneSize / 3 : coneSize / 2;
- assert(0 == coneSize % faceSize);
+ if (!constraintCell) {
+ const int faceSize = coneSize / 2;
+ assert(0 == coneSize % faceSize);
- sieve_type::traits::coneSequence::iterator v_iter = cone->begin();
- for(int i=0; i < faceSize; ++i, ++v_iter)
- faultSieve->addArrow(*v_iter, face, color++);
+ // Use first vertices (negative side of the fault) for fault mesh
+ sieve_type::traits::coneSequence::iterator v_iter = cone->begin();
+ for(int i=0; i < faceSize; ++i, ++v_iter)
+ faultSieve->addArrow(*v_iter, face, color++);
+ } else {
+ const int faceSize = coneSize / 3;
+ assert(0 == coneSize % faceSize);
+
+ // Use last vertices (contraints) for fault mesh
+ sieve_type::traits::coneSequence::iterator v_iter = cone->begin();
+ for(int i=0; i < 2*faceSize; ++i)
+ ++v_iter;
+ for(int i=0; i < faceSize; ++i, ++v_iter)
+ faultSieve->addArrow(*v_iter, face, color++);
+ } // if/else
+
++face;
} // for
(*fault)->setSieve(faultSieve);
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-01-27 19:07:29 UTC (rev 9145)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-01-27 22:26:52 UTC (rev 9146)
@@ -119,13 +119,14 @@
} // for
// Create orientation section for constraint vertices
+ const int cohesiveDim = _faultMesh->getDimension();
const int spaceDim = cs->spaceDim();
const int orientationSize = spaceDim*spaceDim;
_orientation = new real_section_type(mesh->comm(), mesh->debug());
assert(!_orientation.isNull());
- for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int iDim=0; iDim <= cohesiveDim; ++iDim)
_orientation->addSpace();
- assert(spaceDim == _orientation->getNumSpaces());
+ assert(cohesiveDim+1 == _orientation->getNumSpaces());
const std::set<Mesh::point_type>::const_iterator vertConstraintBegin =
_constraintVert.begin();
const std::set<Mesh::point_type>::const_iterator vertConstraintEnd =
@@ -137,7 +138,7 @@
_orientation->setFiberDimension(*v_iter, orientationSize);
// Create fibrations, one for each direction
- for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int iDim=0; iDim <= cohesiveDim; ++iDim)
_orientation->setFiberDimension(*v_iter, spaceDim, iDim);
} // for
mesh->allocate(_orientation);
@@ -150,7 +151,6 @@
assert(!coordinates.isNull());
// Set orientation function
- const int cohesiveDim = mesh->getDimension()-1;
assert(cohesiveDim == _quadrature->cellDim());
assert(spaceDim == _quadrature->spaceDim());
@@ -669,14 +669,44 @@
const char* name,
const ALE::Obj<Mesh>& mesh)
{ // vertexField
- if (0 == strcasecmp("strike_dir", name)) {
- _allocateOutputVertexScalar();
+ const int cohesiveDim = _faultMesh->getDimension();
+
+ if (cohesiveDim > 0 && 0 == strcasecmp("strike_dir", name)) {
+ _allocateOutputVertexVector();
const ALE::Obj<real_section_type>& strikeDir =
_orientation->getFibration(0);
- _projectCohesiveVertexField(&_outputVertexScalar, strikeDir, mesh);
+ _projectCohesiveVertexField(&_outputVertexVector, strikeDir, mesh);
+ *fieldType = meshio::DataWriter::VECTOR_FIELD;
+ return _outputVertexVector;
+ } else if (2 == cohesiveDim && 0 == strcasecmp("dip_dir", name)) {
+ _allocateOutputVertexVector();
+ const ALE::Obj<real_section_type>& dipDir =
+ _orientation->getFibration(1);
+ _projectCohesiveVertexField(&_outputVertexVector, dipDir, mesh);
+ *fieldType = meshio::DataWriter::VECTOR_FIELD;
+ return _outputVertexVector;
+ } else if (0 == strcasecmp("normal_dir", name)) {
+ _allocateOutputVertexVector();
+ const int space =
+ (0 == cohesiveDim) ? 0 : (1 == cohesiveDim) ? 1 : 2;
+ const ALE::Obj<real_section_type>& normalDir =
+ _orientation->getFibration(space);
+ _projectCohesiveVertexField(&_outputVertexVector, normalDir, mesh);
+ *fieldType = meshio::DataWriter::VECTOR_FIELD;
+ return _outputVertexVector;
+ } else if (0 == strcasecmp("final_slip", name)) {
+ _allocateOutputVertexVector();
+ // ADD STUFF HERE
+ //_projectCohesiveVertexField(&_outputVertexVector, finalSlip, mesh);
+ *fieldType = meshio::DataWriter::VECTOR_FIELD;
+ return _outputVertexVector;
+ } else if (0 == strcasecmp("slip_time", name)) {
+ _allocateOutputVertexScalar();
+ // ADD STUFF HERE
+ //_projectCohesiveVertexField(&_outputVertexScalar, slipTime, mesh);
*fieldType = meshio::DataWriter::SCALAR_FIELD;
return _outputVertexScalar;
- } // if
+ } // if/else
// Should not reach this point if requested field was found
std::ostringstream msg;
@@ -718,6 +748,7 @@
const ALE::Obj<Mesh::label_sequence>& vertices =
_faultMesh->depthStratum(0);
_outputVertexScalar->setFiberDimension(vertices, fiberDim);
+ _faultMesh->allocate(_outputVertexScalar);
} // if
} // _allocateOutputVertexScalar
@@ -734,7 +765,8 @@
const ALE::Obj<Mesh::label_sequence>& vertices =
_faultMesh->depthStratum(0);
_outputVertexVector->setFiberDimension(vertices, fiberDim);
- } // if
+ _faultMesh->allocate(_outputVertexVector);
+ } // if
} // _allocateOutputVertexVector
// ----------------------------------------------------------------------
@@ -745,10 +777,9 @@
const ALE::Obj<real_section_type>& fieldCohesive,
const ALE::Obj<Mesh>& mesh)
{ // _projectCohesiveVertexField
+ assert(0 != _quadrature);
// Get cohesive cells
- const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
- assert(!sieve.isNull());
const ALE::Obj<ALE::Mesh::label_sequence>& cellsCohesive =
mesh->getLabelStratum("material-id", id());
assert(!cellsCohesive.isNull());
@@ -756,53 +787,31 @@
cellsCohesive->begin();
const Mesh::label_sequence::iterator cellsCohesiveEnd =
cellsCohesive->end();
+ assert(!fieldCohesive.isNull());
- assert(!fieldCohesive.isNull());
+ // Get fault cells
+ const ALE::Obj<Mesh::label_sequence>& cellsFault =
+ _faultMesh->heightStratum(0);
+ assert(!cellsFault.isNull());
+ const Mesh::label_sequence::iterator cellsFaultBegin =
+ cellsFault->begin();
+ const Mesh::label_sequence::iterator cellsFaultEnd =
+ cellsFault->end();
assert(!fieldFault->isNull());
- // Replace this with simultaneous loop over cohesive cells and fault
- // cells. Restrict on cohesive cells and update on fault cells. Only
- // works if sections are limited to lagrange vertices (include
- // assert on fiber dimensions).
- for (Mesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
+ // Project field at constraint (Lagrange multiplier) vertices
+ // defined over cohesive cells to fault mesh
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ _faultMesh->depthStratum(0);
+ const int fiberDim = fieldCohesive->getFiberDimension(*vertices->begin());
+ const int numBasis = _quadrature->numBasis();
+ double_array cellData(numBasis*fiberDim);
+ for (Mesh::label_sequence::iterator c_iter=cellsCohesive->begin(),
+ f_iter=cellsFault->begin();
c_iter != cellsCohesiveEnd;
- ++c_iter) {
- // Vertices for each cohesive cell are in three groups of N.
- // 0 to N-1: vertices on negative side of the fault
- // N-1 to 2N-1: vertices on positive side of the fault
- // 2N to 3N-1: vertices associated with constraint forces
- const ALE::Obj<sieve_type::traits::coneSequence>& cone =
- sieve->cone(*c_iter);
- assert(!cone.isNull());
- const int coneSize = cone->size();
- assert(coneSize % 3 == 0);
- sieve_type::traits::coneSequence::iterator v_iter = cone->begin();
-
- // Loop over each group of 3 vertices (negative, positive,
- // Lagrange) in cohesive cells.
- for (int i=0; i < coneSize; i += 3) {
- // Get vertex on negative side of the fault
- const sieve_type::traits::coneSequence::iterator v_fault = v_iter;
-
- // Get vertex associated with Lagrange multiplier (constraint)
- ++v_iter;
- ++v_iter;
- const sieve_type::traits::coneSequence::iterator v_cohesive = v_iter;
-
- // Copy values from Lagrange multplier vertex in field over
- // cohesive cells to vertex (on negative side of the fault) in
- // field over fault.
- //
- // NOTE: THIS IS INEFFICIENT BECAUSE WE COPY VALUES MULTIPLE
- // TIMES. If the fields are only defined at the Lagrange
- // multiplier vertices and the fault mesh is ordered the same as
- // the cohseive cells, we could do a direct copy of the values
- // for the entire section.
- assert((*fieldFault)->getFiberDimension(*v_fault) ==
- fieldCohesive->getFiberDimension(*v_cohesive));
- (*fieldFault)->updatePoint(*v_fault,
- fieldCohesive->restrictPoint(*v_cohesive));
- } // for
+ ++c_iter, ++f_iter) {
+ mesh->restrict(fieldCohesive, *c_iter, &cellData[0], cellData.size());
+ _faultMesh->update(*fieldFault, *f_iter, &cellData[0]);
} // for
} // _projectCohesiveVertexField
Modified: short/3D/PyLith/trunk/modulesrc/meshio/meshio.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/meshio/meshio.pyxe.src 2008-01-27 19:07:29 UTC (rev 9145)
+++ short/3D/PyLith/trunk/modulesrc/meshio/meshio.pyxe.src 2008-01-27 22:26:52 UTC (rev 9146)
@@ -722,7 +722,7 @@
Append field over vertices to file.
"""
# create shim for method 'appendVertexField'
- #embed{ void OutputManager_appendVertexField(void* objVptr, double t, char* name, void* fieldVptr, char* fieldTypeName, void* meshVptr)
+ #embed{ void OutputManager_appendVertexField(void* objVptr, double t, char* name, void* fieldVptr, int fieldTypeInt, void* meshVptr)
try {
assert(0 != objVptr);
assert(0 != fieldVptr);
@@ -731,14 +731,8 @@
ALE::Obj<pylith::real_section_type>* field =
(ALE::Obj<pylith::real_section_type>*) fieldVptr;
ALE::Obj<ALE::Mesh>* mesh = (ALE::Obj<ALE::Mesh>*) meshVptr;
- pylith::meshio::DataWriter::FieldEnum fieldType =
- pylith::meshio::DataWriter::OTHER_FIELD;
- if (0 == strcasecmp("scalar field", fieldTypeName))
- fieldType = pylith::meshio::DataWriter::SCALAR_FIELD;
- else if (0 == strcasecmp("vector field", fieldTypeName))
- fieldType = pylith::meshio::DataWriter::VECTOR_FIELD;
- else if (0 == strcasecmp("tensor field", fieldTypeName))
- fieldType = pylith::meshio::DataWriter::TENSOR_FIELD;
+ const pylith::meshio::DataWriter::FieldEnum fieldType =
+ (pylith::meshio::DataWriter::FieldEnum) fieldTypeInt;
((pylith::meshio::OutputManager*) objVptr)->appendVertexField(t, name,
*field, fieldType, *mesh);
} catch (const std::exception& err) {
@@ -767,7 +761,7 @@
Append field over cells to file.
"""
# create shim for method 'appendCellField'
- #embed{ void OutputManager_appendCellField(void* objVptr, double t, char* name, void* fieldVptr, char* fieldTypeName, void* meshVptr)
+ #embed{ void OutputManager_appendCellField(void* objVptr, double t, char* name, void* fieldVptr, int fieldTypeInt, void* meshVptr)
try {
assert(0 != objVptr);
assert(0 != fieldVptr);
@@ -776,14 +770,8 @@
ALE::Obj<pylith::real_section_type>* field =
(ALE::Obj<pylith::real_section_type>*) fieldVptr;
ALE::Obj<ALE::Mesh>* mesh = (ALE::Obj<ALE::Mesh>*) meshVptr;
- pylith::meshio::DataWriter::FieldEnum fieldType =
- pylith::meshio::DataWriter::OTHER_FIELD;
- if (0 == strcasecmp("scalar field", fieldTypeName))
- fieldType = pylith::meshio::DataWriter::SCALAR_FIELD;
- else if (0 == strcasecmp("vector field", fieldTypeName))
- fieldType = pylith::meshio::DataWriter::VECTOR_FIELD;
- else if (0 == strcasecmp("tensor field", fieldTypeName))
- fieldType = pylith::meshio::DataWriter::TENSOR_FIELD;
+ const pylith::meshio::DataWriter::FieldEnum fieldType =
+ (pylith::meshio::DataWriter::FieldEnum) fieldTypeInt;
((pylith::meshio::OutputManager*) objVptr)->appendCellField(t, name,
*field,
fieldType,
Modified: short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py 2008-01-27 19:07:29 UTC (rev 9145)
+++ short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py 2008-01-27 22:26:52 UTC (rev 9146)
@@ -73,9 +73,7 @@
self.availableFields = \
{'vertex': \
- {'info': ["strike_dir",
- "dip_dir",
- "normal_dir",
+ {'info': ["normal_dir",
"final_slip",
"slip_time"],
'data': ["slip"]},
@@ -96,6 +94,12 @@
self.eqsrc.preinitialize()
self.cppHandle.eqsrc = self.eqsrc.cppHandle
self.output.preinitialize(self)
+
+ if mesh.dimension() == 2:
+ self.availableFields['vertex']['info'] += ["strike_dir"]
+ elif mesh.dimension() == 3:
+ self.availableFields['vertex']['info'] += ["strike_dir",
+ "dip_dir"]
return
@@ -141,36 +145,18 @@
return
- def getVertexField(self, name, mesh):
+ def getVertexField(self, name):
"""
Get vertex field.
"""
- field = self.cppHandle.getVertexField(name, mesh.cppHandle)
- fieldType = None
- if name in ["strike_dir",
- "dip_dir",
- "normal_dir",
- "final_slip",
- "traction_change"]:
- fieldType = "vector field"
- elif name in ["slip_time"]:
- fieldType = "scalar field"
- else:
- raise ValueError, "Vertex field '%s' not available." % name
- return (field, fieldType)
+ return self.cppHandle.vertexField(name, self.mesh.cppHandle)
- def getCellField(self, name, mesh):
+ def getCellField(self, name):
"""
Get cell field.
"""
- field = self.cppHandle.getVertexField(name, mesh.cppHandle)
- fieldType = None
- if name in ["traction_change"]:
- fieldType = "vector field"
- else:
- raise ValueError, "Cell field '%s' not available." % name
- return (field, fieldType)
+ return self.cppHandle.cellField(name, self.mesh.cppHandle)
# PRIVATE METHODS ////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/pylith/meshio/OutputFaultKin.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/OutputFaultKin.py 2008-01-27 19:07:29 UTC (rev 9145)
+++ short/3D/PyLith/trunk/pylith/meshio/OutputFaultKin.py 2008-01-27 22:26:52 UTC (rev 9146)
@@ -50,9 +50,7 @@
import pyre.inventory
vertexInfoFields = pyre.inventory.list("vertex_info_fields",
- default=["strike_dir",
- "dip_dir",
- "normal_dir",
+ default=["normal_dir",
"final_slip",
"slip_time"])
vertexInfoFields.meta['tip'] = "Names of vertex info fields to output."
Modified: short/3D/PyLith/trunk/pylith/meshio/OutputManager.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/OutputManager.py 2008-01-27 19:07:29 UTC (rev 9145)
+++ short/3D/PyLith/trunk/pylith/meshio/OutputManager.py 2008-01-27 22:26:52 UTC (rev 9146)
@@ -201,15 +201,17 @@
self.cppHandle.openTimeStep(t.value,
mesh.cppHandle, mesh.coordsys.cppHandle)
- #for name in self.vertexInfoFields:
- # (field, fieldType) = self.dataProvider.getVertexField(name)
- # self.cppHandle.appendVertexField(t.value, name, field, fieldType,
- # mesh.cppHandle)
+ for name in self.vertexInfoFields:
+ print "Getting field '%s'." % name
+ (field, fieldType) = self.dataProvider.getVertexField(name)
+ print "Writing field '%s'." % name
+ self.cppHandle.appendVertexField(t.value, name, field, fieldType,
+ mesh.cppHandle)
- #for name in self.cellInfoFields:
- # (field, fieldType) = self.dataProvider.getCellField(name)
- # self.cppHandle.appendCellField(t.value, name, field, fieldType,
- # mesh.cppHandle)
+ for name in self.cellInfoFields:
+ (field, fieldType) = self.dataProvider.getCellField(name)
+ self.cppHandle.appendCellField(t.value, name, field, fieldType,
+ mesh.cppHandle)
self.cppHandle.closeTimeStep()
self.close()
@@ -328,9 +330,9 @@
"Field type: '%s'\n" \
"Data type: '%s'\n" % (fieldCategory, dataCategory)
msg += "Available fields: "
- for name in available:
+ for name in available[fieldCategory][dataCategory]:
msg += " '%s'" % name
- msg += "\n"
+ msg += "\n"
msg += "Fields not available: "
for name in notavailable:
msg += " '%s'" % name
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2008-01-27 19:07:29 UTC (rev 9145)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2008-01-27 22:26:52 UTC (rev 9146)
@@ -313,7 +313,7 @@
fieldType = None
if name == "displacements":
field = self.fields.getSolution()
- fieldType = "vector field"
+ fieldType = 1 # vector field
else:
raise ValueError, "Vertex field '%s' not available." % name
return (field, fieldType)
More information about the cig-commits
mailing list