[cig-commits] r9144 - in short/3D/PyLith/trunk: libsrc/faults libsrc/meshio modulesrc/faults modulesrc/meshio pylith/faults pylith/meshio
brad at geodynamics.org
brad at geodynamics.org
Sun Jan 27 10:18:29 PST 2008
Author: brad
Date: 2008-01-27 10:18:29 -0800 (Sun, 27 Jan 2008)
New Revision: 9144
Modified:
short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.hh
short/3D/PyLith/trunk/libsrc/faults/Fault.cc
short/3D/PyLith/trunk/libsrc/faults/Fault.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
short/3D/PyLith/trunk/libsrc/meshio/DataWriterVTK.cc
short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src
short/3D/PyLith/trunk/modulesrc/meshio/meshio.pyxe.src
short/3D/PyLith/trunk/pylith/faults/Fault.py
short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py
short/3D/PyLith/trunk/pylith/meshio/OutputManager.py
Log:
Worked on output of fault information. Fault mesh is missing coordinates of vertices.
Modified: short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc 2008-01-27 18:18:29 UTC (rev 9144)
@@ -640,7 +640,56 @@
if (debug) coordinates->view("Coordinates with shadow vertices");
} // createCohesiveCells
+// ----------------------------------------------------------------------
+// Form a parallel fault mesh using the cohesive cell information
void
+pylith::faults::CohesiveTopology::createParallel(ALE::Obj<Mesh>* fault,
+ const ALE::Obj<Mesh>& mesh,
+ const int materialId,
+ const bool constraintCell)
+{
+ assert(0 != fault);
+ *fault = new Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
+
+ const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
+ const ALE::Obj<sieve_type> faultSieve =
+ new sieve_type(sieve->comm(), sieve->debug());
+ const ALE::Obj<Mesh::label_sequence>& cohesiveCells =
+ mesh->getLabelStratum("material-id", materialId);
+ const Mesh::label_sequence::iterator cBegin = cohesiveCells->begin();
+ const Mesh::label_sequence::iterator cEnd = cohesiveCells->end();
+ const int sieveEnd = sieve->base()->size() + sieve->cap()->size();
+ const int numFaces = cohesiveCells->size();
+ int globalSieveEnd = 0;
+ int globalFaceOffset = 0;
+
+ MPI_Allreduce((void *) &sieveEnd, (void *) &globalSieveEnd, 1,
+ MPI_INT, MPI_SUM, sieve->comm());
+ MPI_Scan((void *) &numFaces, (void *) &globalFaceOffset, 1,
+ MPI_INT, MPI_SUM, sieve->comm());
+ int face = globalSieveEnd + globalFaceOffset;
+ for(Mesh::label_sequence::iterator c_iter = cBegin;
+ c_iter != cEnd;
+ ++c_iter) {
+ const ALE::Obj<sieve_type::traits::coneSequence>& cone =
+ sieve->cone(*c_iter);
+ int color = 0;
+
+ const int coneSize = cone->size();
+ const int faceSize = (constraintCell) ? coneSize / 3 : 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++);
+ ++face;
+ } // for
+ (*fault)->setSieve(faultSieve);
+ (*fault)->stratify();
+}
+
+// ----------------------------------------------------------------------
+void
pylith::faults::CohesiveTopology::classifyCells(const ALE::Obj<Mesh::sieve_type>& sieve,
const Mesh::point_type& vertex,
const int depth,
@@ -729,6 +778,7 @@
noReplaceCells.insert(vNoReplaceCells.begin(), vNoReplaceCells.end());
}
+// ----------------------------------------------------------------------
template<class InputPoints>
bool
pylith::faults::CohesiveTopology::_compatibleOrientation(const ALE::Obj<Mesh>& mesh,
@@ -766,6 +816,7 @@
return compatible;
}
+// ----------------------------------------------------------------------
void
pylith::faults::CohesiveTopology::_replaceCell(const Obj<sieve_type>& sieve,
const Mesh::point_type cell,
@@ -798,6 +849,7 @@
}
}
+// ----------------------------------------------------------------------
template<class InputPoints>
void
pylith::faults::CohesiveTopology::_computeCensoredDepth(const ALE::Obj<Mesh>& mesh,
@@ -824,43 +876,7 @@
if(modifiedPoints->size() > 0) {
_computeCensoredDepth(mesh, depth, sieve, sieve->support(modifiedPoints), firstCohesiveCell, modifiedPoints);
}
-};
+}
-// Form a parallel fault mesh using the cohesive cell information
-void
-pylith::faults::CohesiveTopology::createParallel(ALE::Obj<Mesh>* fault,
- const ALE::Obj<Mesh>& mesh,
- const int materialId)
-{
- assert(0 != fault);
- *fault = new Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
- const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
- const ALE::Obj<sieve_type> faultSieve = new sieve_type(sieve->comm(), sieve->debug());
- const ALE::Obj<Mesh::label_sequence>& cohesiveCells = mesh->getLabelStratum("material-id", materialId);
- const Mesh::label_sequence::iterator cBegin = cohesiveCells->begin();
- const Mesh::label_sequence::iterator cEnd = cohesiveCells->end();
- const int sieveEnd = sieve->base()->size() + sieve->cap()->size();
- const int numFaces = cohesiveCells->size();
- int globalSieveEnd;
- int globalFaceOffset;
- int face;
-
- MPI_Allreduce((void *) &sieveEnd, (void *) &globalSieveEnd, 1, MPI_INT, MPI_SUM, sieve->comm());
- MPI_Scan((void *) &numFaces, (void *) &globalFaceOffset, 1, MPI_INT, MPI_SUM, sieve->comm());
- face = globalSieveEnd + globalFaceOffset;
- for(Mesh::label_sequence::iterator c_iter = cBegin; c_iter != cEnd; ++c_iter) {
- const ALE::Obj<sieve_type::traits::coneSequence>& cone = sieve->cone(*c_iter);
- const sieve_type::traits::coneSequence::iterator begin = cone->begin();
- const sieve_type::traits::coneSequence::iterator end = cone->end();
- int color = 0;
-
- for(sieve_type::traits::coneSequence::iterator v_iter = begin; v_iter != end; ++v_iter) {
- faultSieve->addArrow(*v_iter, face, color++);
- }
- ++face;
- }
- (*fault)->setSieve(faultSieve);
-};
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.hh 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.hh 2008-01-27 18:18:29 UTC (rev 9144)
@@ -38,6 +38,7 @@
// PUBLIC METHODS /////////////////////////////////////////////////////
public :
+
/** Create cohesive cells.
*
* @param fault Finite-element mesh of fault (output)
@@ -53,8 +54,22 @@
const ALE::Obj<Mesh>& mesh,
const ALE::Obj<Mesh::int_section_type>& groupField,
const int materialId,
- const bool constraintCell = false);
+ const bool constraintCell =false);
+ /** Create (distributed) fault mesh from cohesive cells.
+ *
+ * @param fault Finite-element mesh of fault (output)
+ * @param mesh Finite-element mesh
+ * @param materialId Material id for cohesive elements.
+ * @param constraintCell True if creating cells constrained with
+ * Lagrange multipliers that require extra vertices, false otherwise
+ */
+ static
+ void createParallel(ALE::Obj<Mesh>* fault,
+ const ALE::Obj<Mesh>& mesh,
+ const int materialId,
+ const bool constraintCell =false);
+
// PRIVATE METHODS ////////////////////////////////////////////////////
private :
/** Get number of vertices on face.
@@ -136,9 +151,6 @@
PointSet& noReplaceCells,
const int debug);
- static void createParallel(ALE::Obj<Mesh>* fault,
- const ALE::Obj<Mesh>& mesh,
- const int materialId);
}; // class CohesiveTopology
#endif // pylith_faults_cohesivetopology_hh
Modified: short/3D/PyLith/trunk/libsrc/faults/Fault.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/Fault.cc 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/libsrc/faults/Fault.cc 2008-01-27 18:18:29 UTC (rev 9144)
@@ -18,7 +18,8 @@
// Default constructor.
pylith::faults::Fault::Fault(void) :
_id(0),
- _label("")
+ _label(""),
+ _faultMesh(0)
{ // constructor
} // constructor
@@ -28,5 +29,13 @@
{ // destructor
} // destructor
+// ----------------------------------------------------------------------
+// Get mesh associated with fault fields.
+const ALE::Obj<ALE::Mesh>&
+pylith::faults::Fault:: faultMesh(void) const
+{ // faultMesh
+ return _faultMesh;
+} // faultMesh
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/Fault.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/Fault.hh 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/libsrc/faults/Fault.hh 2008-01-27 18:18:29 UTC (rev 9144)
@@ -20,8 +20,9 @@
#if !defined(pylith_faults_fault_hh)
#define pylith_faults_fault_hh
-#include "pylith/utils/sievefwd.hh" // USES PETSc Mesh
+#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh, real_section_type
#include "pylith/utils/arrayfwd.hh" // USES double_array
+#include "pylith/meshio/DataWriter.hh" // USES FieldEnum
#include <string> // HASA std::string
@@ -111,9 +112,38 @@
const double_array& normalDir,
spatialdata::spatialdb::SpatialDB* matDB) = 0;
- // PROTECTED METHODS //////////////////////////////////////////////////
-protected :
+ /** Get mesh associated with fault fields.
+ *
+ * @returns PETSc mesh object
+ */
+ const ALE::Obj<ALE::Mesh>& faultMesh(void) const;
+ /** Get vertex field associated with integrator.
+ *
+ * @param fieldType Type of field.
+ * @param name Name of vertex field.
+ * @param mesh PETSc mesh for problem.
+ * @returns Vertex field.
+ */
+ virtual
+ const ALE::Obj<real_section_type>&
+ vertexField(meshio::DataWriter::FieldEnum* fieldType,
+ const char* name,
+ const ALE::Obj<Mesh>& mesh) = 0;
+
+ /** Get cell field associated with integrator.
+ *
+ * @param fieldType Type of field.
+ * @param name Name of vertex field.
+ * @param mesh PETSc mesh for problem.
+ * @returns Cell field.
+ */
+ virtual
+ const ALE::Obj<real_section_type>&
+ cellField(meshio::DataWriter::FieldEnum* fieldType,
+ const char* name,
+ const ALE::Obj<Mesh>& mesh) = 0;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
@@ -123,6 +153,11 @@
/// Not implemented
const Fault& operator=(const Fault& m);
+// PROTECTED MEMBERS ////////////////////////////////////////////////////
+protected :
+
+ ALE::Obj<ALE::Mesh> _faultMesh; ///< Mesh over fault surface
+
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2008-01-27 18:18:29 UTC (rev 9144)
@@ -23,8 +23,7 @@
// ----------------------------------------------------------------------
// Default constructor.
-pylith::faults::FaultCohesive::FaultCohesive(void) :
- _faultMesh(new ALE::Obj<ALE::Mesh>)
+pylith::faults::FaultCohesive::FaultCohesive(void)
{ // constructor
} // constructor
@@ -32,7 +31,6 @@
// Destructor.
pylith::faults::FaultCohesive::~FaultCohesive(void)
{ // destructor
- delete _faultMesh; _faultMesh = 0;
} // destructor
// ----------------------------------------------------------------------
@@ -47,7 +45,8 @@
mesh->getIntSection(label());
assert(!groupField.isNull());
- CohesiveTopology::create(_faultMesh, mesh, groupField, id(),
+ ALE::Obj<Mesh> faultMesh;
+ CohesiveTopology::create(&faultMesh, mesh, groupField, id(),
_useLagrangeConstraints());
} // adjustTopology
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2008-01-27 18:18:29 UTC (rev 9144)
@@ -74,11 +74,6 @@
/// Not implemented
const FaultCohesive& operator=(const FaultCohesive& m);
-// PROTECTED MEMBERS ////////////////////////////////////////////////////
-protected :
-
- ALE::Obj<ALE::Mesh>* _faultMesh;
-
}; // class FaultCohesive
#endif // pylith_faults_faultcohesive_hh
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2008-01-27 18:18:29 UTC (rev 9144)
@@ -112,4 +112,42 @@
} // if
} // for
} // verifyConfiguration
+
+// ----------------------------------------------------------------------
+// Get vertex field associated with integrator.
+const ALE::Obj<pylith::real_section_type>&
+pylith::faults::FaultCohesiveDyn::vertexField(
+ meshio::DataWriter::FieldEnum* fieldType,
+ const char* name,
+ const ALE::Obj<Mesh>& mesh)
+{ // vertexField
+ // Should not reach this point if requested field was found
+ std::ostringstream msg;
+ msg << "Request for unknown vertex field '" << name
+ << "' for fault '" << label() << ".";
+ throw std::runtime_error(msg.str());
+
+ // Return generic section to satisfy member function definition.
+ //return _outputVertxScalar;
+} // vertexField
+
+// ----------------------------------------------------------------------
+// Get cell field associated with integrator.
+const ALE::Obj<pylith::real_section_type>&
+pylith::faults::FaultCohesiveDyn::cellField(
+ meshio::DataWriter::FieldEnum* fieldType,
+ const char* name,
+ const ALE::Obj<Mesh>& mesh)
+{ // cellField
+ // Should not reach this point if requested field was found
+ std::ostringstream msg;
+ msg << "Request for unknown vertex field '" << name
+ << "' for fault '" << label() << ".";
+ throw std::runtime_error(msg.str());
+
+ // Return generic section to satisfy member function definition.
+ //return _outputCellVector;
+} // cellField
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2008-01-27 18:18:29 UTC (rev 9144)
@@ -101,6 +101,30 @@
*/
void verifyConfiguration(const ALE::Obj<Mesh>& mesh);
+ /** Get vertex field associated with integrator.
+ *
+ * @param fieldType Type of field.
+ * @param name Name of vertex field.
+ * @param mesh PETSc mesh for problem.
+ * @returns Vertex field.
+ */
+ const ALE::Obj<real_section_type>&
+ vertexField(meshio::DataWriter::FieldEnum* fieldType,
+ const char* name,
+ 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.
+ * @returns Cell field.
+ */
+ const ALE::Obj<real_section_type>&
+ cellField(meshio::DataWriter::FieldEnum* fieldType,
+ const char* name,
+ const ALE::Obj<Mesh>& mesh);
+
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-01-27 18:18:29 UTC (rev 9144)
@@ -15,6 +15,7 @@
#include "FaultCohesiveKin.hh" // implementation of object methods
#include "EqKinSrc.hh" // USES EqKinSrc
+#include "CohesiveTopology.hh" // USES CohesiveTopology
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
#include "pylith/feassemble/CellGeometry.hh" // USES CellGeometry
@@ -73,6 +74,9 @@
throw std::runtime_error("Normal direction for fault orientation must be "
"a vector with 3 components.");
+ CohesiveTopology::createParallel(&_faultMesh, mesh, id(),
+ _useLagrangeConstraints());
+
/* First find vertices associated with Lagrange multiplier
* constraints in cohesive cells and compute the orientation of the
* fault at these locations.
@@ -101,7 +105,6 @@
sieve->cone(*c_iter);
assert(!cone.isNull());
const sieve_type::traits::coneSequence::iterator vBegin = cone->begin();
- const sieve_type::traits::coneSequence::iterator vEnd = cone->end();
const int coneSize = cone->size();
assert(coneSize % 3 == 0);
sieve_type::traits::coneSequence::iterator v_iter = vBegin;
@@ -120,6 +123,9 @@
const int orientationSize = spaceDim*spaceDim;
_orientation = new real_section_type(mesh->comm(), mesh->debug());
assert(!_orientation.isNull());
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _orientation->addSpace();
+ assert(spaceDim == _orientation->getNumSpaces());
const std::set<Mesh::point_type>::const_iterator vertConstraintBegin =
_constraintVert.begin();
const std::set<Mesh::point_type>::const_iterator vertConstraintEnd =
@@ -127,8 +133,13 @@
const int vertConstraintSize = _constraintVert.size();
for (std::set<Mesh::point_type>::const_iterator v_iter=vertConstraintBegin;
v_iter != vertConstraintEnd;
- ++v_iter)
+ ++v_iter) {
_orientation->setFiberDimension(*v_iter, orientationSize);
+
+ // Create fibrations, one for each direction
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _orientation->setFiberDimension(*v_iter, spaceDim, iDim);
+ } // for
mesh->allocate(_orientation);
// Compute orientation of fault at constraint vertices
@@ -258,7 +269,7 @@
PetscLogFlopsNoCheck(5 + vertConstraintSize * 3);
} // if
- _eqsrc->initialize(mesh, *_faultMesh, _constraintVert, cs);
+ _eqsrc->initialize(mesh, _faultMesh, _constraintVert, cs);
// Establish pairing between constraint vertices and first cell they
// appear in to prevent overlap in integrating Jacobian
@@ -650,5 +661,150 @@
} // for
} // verifyConfiguration
+// ----------------------------------------------------------------------
+// Get vertex field associated with integrator.
+const ALE::Obj<pylith::real_section_type>&
+pylith::faults::FaultCohesiveKin::vertexField(
+ meshio::DataWriter::FieldEnum* fieldType,
+ const char* name,
+ const ALE::Obj<Mesh>& mesh)
+{ // vertexField
+ if (0 == strcasecmp("strike_dir", name)) {
+ _allocateOutputVertexScalar();
+ const ALE::Obj<real_section_type>& strikeDir =
+ _orientation->getFibration(0);
+ _projectCohesiveVertexField(&_outputVertexScalar, strikeDir, mesh);
+ *fieldType = meshio::DataWriter::SCALAR_FIELD;
+ return _outputVertexScalar;
+ } // if
+ // Should not reach this point if requested field was found
+ std::ostringstream msg;
+ msg << "Request for unknown vertex field '" << name
+ << "' for fault '" << label() << ".";
+ throw std::runtime_error(msg.str());
+
+ // Return generic section to satisfy member function definition.
+ return _outputVertexScalar;
+} // vertexField
+
+// ----------------------------------------------------------------------
+// Get cell field associated with integrator.
+const ALE::Obj<pylith::real_section_type>&
+pylith::faults::FaultCohesiveKin::cellField(
+ meshio::DataWriter::FieldEnum* fieldType,
+ const char* name,
+ const ALE::Obj<Mesh>& mesh)
+{ // cellField
+ // Should not reach this point if requested field was found
+ std::ostringstream msg;
+ msg << "Request for unknown vertex field '" << name
+ << "' for fault '" << label() << ".";
+ throw std::runtime_error(msg.str());
+
+ // Return generic section to satisfy member function definition.
+ //return _outputCellVector;
+} // cellField
+
+// ----------------------------------------------------------------------
+// Allocate scalar field for output.
+void
+pylith::faults::FaultCohesiveKin::_allocateOutputVertexScalar(void)
+{ // _allocateOutputVertexScalar
+ const int fiberDim = 1;
+ if (_outputVertexScalar.isNull()) {
+ _outputVertexScalar = new real_section_type(_faultMesh->comm(),
+ _faultMesh->debug());
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ _faultMesh->depthStratum(0);
+ _outputVertexScalar->setFiberDimension(vertices, fiberDim);
+ } // if
+} // _allocateOutputVertexScalar
+
+// ----------------------------------------------------------------------
+// Allocate vector field for output.
+void
+pylith::faults::FaultCohesiveKin::_allocateOutputVertexVector(void)
+{ // _allocateOutputVertexVector
+ assert(0 != _quadrature);
+ const int fiberDim = _quadrature->spaceDim();
+ if (_outputVertexVector.isNull()) {
+ _outputVertexVector = new real_section_type(_faultMesh->comm(),
+ _faultMesh->debug());
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ _faultMesh->depthStratum(0);
+ _outputVertexVector->setFiberDimension(vertices, fiberDim);
+ } // if
+} // _allocateOutputVertexVector
+
+// ----------------------------------------------------------------------
+// Project field defined over cohesive cells to fault mesh.
+void
+pylith::faults::FaultCohesiveKin::_projectCohesiveVertexField(
+ ALE::Obj<real_section_type>* fieldFault,
+ const ALE::Obj<real_section_type>& fieldCohesive,
+ const ALE::Obj<Mesh>& mesh)
+{ // _projectCohesiveVertexField
+
+ // 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());
+ const Mesh::label_sequence::iterator cellsCohesiveBegin =
+ cellsCohesive->begin();
+ const Mesh::label_sequence::iterator cellsCohesiveEnd =
+ cellsCohesive->end();
+
+ assert(!fieldCohesive.isNull());
+ 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;
+ 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
+ } // for
+} // _projectCohesiveVertexField
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2008-01-27 18:18:29 UTC (rev 9144)
@@ -124,6 +124,30 @@
*/
void verifyConfiguration(const ALE::Obj<Mesh>& mesh);
+ /** Get vertex field associated with integrator.
+ *
+ * @param fieldType Type of field.
+ * @param name Name of vertex field.
+ * @param mesh PETSc mesh for problem.
+ * @returns Vertex field.
+ */
+ const ALE::Obj<real_section_type>&
+ vertexField(meshio::DataWriter::FieldEnum* fieldType,
+ const char* name,
+ 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.
+ * @returns Cell field.
+ */
+ const ALE::Obj<real_section_type>&
+ cellField(meshio::DataWriter::FieldEnum* fieldType,
+ const char* name,
+ const ALE::Obj<Mesh>& mesh);
+
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -134,6 +158,27 @@
*/
bool _useLagrangeConstraints(void) const;
+ // PRIVATE METHODS ////////////////////////////////////////////////////
+private :
+
+ /// Allocate scalar field for output of vertex information.
+ void _allocateOutputVertexScalar(void);
+
+ /// Allocate vector field for output of vertex information.
+ void _allocateOutputVertexVector(void);
+
+ /** Project field defined over cohesive cells to fault mesh.
+ *
+ * @param fieldFault Field defined over fault mesh.
+ * @param fieldCohesive Field defined over cohesive cells.
+ * @param mesh PETSc mesh for problem.
+ */
+ void _projectCohesiveVertexField(
+ ALE::Obj<real_section_type>* fieldFault,
+ const ALE::Obj<real_section_type>& fieldCohesive,
+ const ALE::Obj<Mesh>& mesh);
+
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
@@ -164,6 +209,12 @@
/// each vertex).
ALE::Obj<int_section_type> _constraintCell;
+ /// Scalar field for output of vertex information.
+ ALE::Obj<real_section_type> _outputVertexScalar;
+
+ /// Vector field for output of vertex information.
+ ALE::Obj<real_section_type> _outputVertexVector;
+
}; // class FaultCohesiveKin
#include "FaultCohesiveKin.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/meshio/DataWriterVTK.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/DataWriterVTK.cc 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/libsrc/meshio/DataWriterVTK.cc 2008-01-27 18:18:29 UTC (rev 9144)
@@ -56,6 +56,9 @@
const ALE::Obj<ALE::Mesh>& mesh,
const spatialdata::geocoords::CoordSys* csMesh)
{ // openTimeStep
+ assert(!mesh.isNull());
+ assert(0 != csMesh);
+
try {
PetscErrorCode err;
@@ -92,6 +95,11 @@
msg << "Error while preparing for writing data to VTK file "
<< _filename << " at time " << t << ".\n" << err.what();
throw std::runtime_error(msg.str());
+ } catch (const ALE::Exception& err) {
+ std::ostringstream msg;
+ msg << "Error while preparing for writing data to VTK file "
+ << _filename << " at time " << t << ".\n" << err.msg();
+ throw std::runtime_error(msg.str());
} catch (...) {
std::ostringstream msg;
msg << "Unknown error while preparing for writing data to VTK file "
Modified: short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src 2008-01-27 18:18:29 UTC (rev 9144)
@@ -39,6 +39,7 @@
void* malloc(size_t size)
void free(void* mem)
+# ----------------------------------------------------------------------
cdef void Fault_destructor(void* obj):
"""
Destroy Fault object.
@@ -51,33 +52,6 @@
Fault_destructor_cpp(obj)
return
-cdef void EqKinSrc_destructor(void* obj):
- """
- Destroy EqKinSrc object.
- """
- # create shim for destructor
- #embed{ void EqKinSrc_destructor_cpp(void* objVptr)
- pylith::faults::EqKinSrc* src = (pylith::faults::EqKinSrc*) objVptr;
- delete src;
- #}embed
- EqKinSrc_destructor_cpp(obj)
- return
-
-
-cdef void SlipTimeFn_destructor(void* obj):
- """
- Destroy SlipTimeFn object.
- """
- # create shim for destructor
- #embed{ void SlipTimeFn_destructor_cpp(void* objVptr)
- pylith::faults::SlipTimeFn* f = (pylith::faults::SlipTimeFn*) objVptr;
- delete f;
- #}embed
- SlipTimeFn_destructor_cpp(obj)
- return
-
-
-# ----------------------------------------------------------------------
cdef class Fault:
cdef void* thisptr # Pointer to C++ object
@@ -187,6 +161,79 @@
return
+ def faultMesh(self, mesh):
+ """
+ Get mesh associated with fault fields.
+ """
+ # create shim for method 'faultMesh'
+ #embed{ void Fault_faultMesh(void* objVptr, void* meshVptr)
+ try {
+ assert(0 != objVptr);
+ assert(0 != meshVptr);
+ ALE::Obj<ALE::Mesh>* mesh = (ALE::Obj<ALE::Mesh>*) meshVptr;
+ *mesh = ((pylith::faults::Fault*) objVptr)->faultMesh();
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+
+ if not mesh.name == "pylith_topology_Mesh":
+ raise TypeError, \
+ "Argument 'mesh' must be extension module type " \
+ "'pylith::topology::Mesh'."
+ Fault_faultMesh(self.thisptr, ptrFromHandle(mesh))
+ return
+
+
+ def vertexField(self, name, mesh):
+ """
+ Get vertex field.
+ """
+ # create shim for method 'vertexField'
+ #embed{ void* Fault_vertexField(void* objVptr, int* fieldPtr, char* name, void* meshVptr)
+ void* result = 0;
+ try {
+ assert(0 != objVptr);
+ assert(0 != fieldPtr);
+ assert(0 != name);
+ assert(0 != meshVptr);
+ pylith::faults::Fault* fault = (pylith::faults::Fault*) objVptr;
+ ALE::Obj<ALE::Mesh>* mesh = (ALE::Obj<ALE::Mesh>*) meshVptr;
+ pylith::meshio::DataWriter::FieldEnum fieldType;
+ const ALE::Obj<pylith::real_section_type>& field =
+ fault->vertexField(&fieldType, name, *mesh);
+ *fieldPtr = fieldType;
+ result = (void*) &field;
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ return result;
+ #}embed
+ if not mesh.name == "pylith_topology_Mesh":
+ raise TypeError, \
+ "Argument 'mesh' must be extension module type " \
+ "'pylith::topology::Mesh'."
+ cdef void* ptr
+ cdef int fieldType
+ ptr = Fault_vertexField(self.thisptr, &fieldType, name, ptrFromHandle(mesh))
+ return (PyCObject_FromVoidPtr(ptr, NULL), fieldType)
+
+
+
def _createHandle(self):
"""
Wrap pointer to C++ object in PyCObject.
@@ -621,6 +668,19 @@
# ----------------------------------------------------------------------
+cdef void EqKinSrc_destructor(void* obj):
+ """
+ Destroy EqKinSrc object.
+ """
+ # create shim for destructor
+ #embed{ void EqKinSrc_destructor_cpp(void* objVptr)
+ pylith::faults::EqKinSrc* src = (pylith::faults::EqKinSrc*) objVptr;
+ delete src;
+ #}embed
+ EqKinSrc_destructor_cpp(obj)
+ return
+
+
cdef class EqKinSrc:
cdef void* thisptr # Pointer to C++ object
@@ -695,6 +755,19 @@
# ----------------------------------------------------------------------
+cdef void SlipTimeFn_destructor(void* obj):
+ """
+ Destroy SlipTimeFn object.
+ """
+ # create shim for destructor
+ #embed{ void SlipTimeFn_destructor_cpp(void* objVptr)
+ pylith::faults::SlipTimeFn* f = (pylith::faults::SlipTimeFn*) objVptr;
+ delete f;
+ #}embed
+ SlipTimeFn_destructor_cpp(obj)
+ return
+
+
cdef class SlipTimeFn:
cdef void* thisptr # Pointer to C++ object
Modified: short/3D/PyLith/trunk/modulesrc/meshio/meshio.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/meshio/meshio.pyxe.src 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/modulesrc/meshio/meshio.pyxe.src 2008-01-27 18:18:29 UTC (rev 9144)
@@ -689,7 +689,7 @@
raise TypeError, \
"Argument must be extension module type 'CoordSys'."
OutputManager_openTimeStep(self.thisptr, t,
- ptrFromHandle(mesh), ptrFromHandle(cs))
+ ptrFromHandle(mesh), ptrFromHandle(cs))
return
Modified: short/3D/PyLith/trunk/pylith/faults/Fault.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/Fault.py 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/pylith/faults/Fault.py 2008-01-27 18:18:29 UTC (rev 9144)
@@ -177,7 +177,7 @@
from pylith.topology.Mesh import Mesh
self.faultMesh = Mesh()
self.faultMesh.initialize(self.mesh.coordsys)
- #self.faultMesh.cppHandle = self.cppHandle.faultMesh() # TODO
+ self.cppHandle.faultMesh(self.faultMesh.cppHandle)
return
Modified: short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py 2008-01-27 18:18:29 UTC (rev 9144)
@@ -121,6 +121,7 @@
self.eqsrc.initialize()
FaultCohesive.initialize(self)
self.output.initialize(self.quadrature.cppHandle)
+ self.output.writeInfo()
self._logger.eventEnd(logEvent)
return
@@ -134,17 +135,17 @@
self._logger.eventBegin(logEvent)
self._info.log("Writing fault data.")
- #self.cppHandle.writeData(t+dt)
+ #self.output.writeData(t+dt)
self._logger.eventEnd(logEvent)
return
- def getVertexField(self, name):
+ def getVertexField(self, name, mesh):
"""
Get vertex field.
"""
- field = self.cppHandle.getVertexField(name) # TODO
+ field = self.cppHandle.getVertexField(name, mesh.cppHandle)
fieldType = None
if name in ["strike_dir",
"dip_dir",
@@ -159,11 +160,11 @@
return (field, fieldType)
- def getCellField(self):
+ def getCellField(self, name, mesh):
"""
Get cell field.
"""
- field = self.cppHandle.getVertexField(name) # TODO
+ field = self.cppHandle.getVertexField(name, mesh.cppHandle)
fieldType = None
if name in ["traction_change"]:
fieldType = "vector field"
Modified: short/3D/PyLith/trunk/pylith/meshio/OutputManager.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/OutputManager.py 2008-01-25 22:38:36 UTC (rev 9143)
+++ short/3D/PyLith/trunk/pylith/meshio/OutputManager.py 2008-01-27 18:18:29 UTC (rev 9144)
@@ -201,18 +201,18 @@
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:
+ # (field, fieldType) = self.dataProvider.getVertexField(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()
+ self.cppHandle.closeTimeStep()
+ self.close()
self._logger.eventEnd(logEvent)
return
More information about the cig-commits
mailing list