[cig-commits] r19108 - in short/3D/PyLith/trunk: libsrc/pylith libsrc/pylith/meshio libsrc/pylith/utils pylith/meshio
brad at geodynamics.org
brad at geodynamics.org
Thu Oct 20 21:56:38 PDT 2011
Author: brad
Date: 2011-10-20 21:56:37 -0700 (Thu, 20 Oct 2011)
New Revision: 19108
Modified:
short/3D/PyLith/trunk/libsrc/pylith/Makefile.am
short/3D/PyLith/trunk/libsrc/pylith/meshio/Makefile.am
short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputManager.hh
short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.hh
short/3D/PyLith/trunk/libsrc/pylith/meshio/meshiofwd.hh
short/3D/PyLith/trunk/libsrc/pylith/utils/petscfwd.h
short/3D/PyLith/trunk/pylith/meshio/OutputSolnPoints.py
Log:
More work on output at arbitrary points.
Modified: short/3D/PyLith/trunk/libsrc/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/Makefile.am 2011-10-20 23:26:58 UTC (rev 19107)
+++ short/3D/PyLith/trunk/libsrc/pylith/Makefile.am 2011-10-21 04:56:37 UTC (rev 19108)
@@ -153,6 +153,7 @@
utils/EventLogger.cc \
utils/TestArray.cc
+#meshio/OutputSolnPoints.cc
# TEMPORARY
libpylith_la_SOURCES += \
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/Makefile.am 2011-10-20 23:26:58 UTC (rev 19107)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/Makefile.am 2011-10-21 04:56:37 UTC (rev 19108)
@@ -41,6 +41,7 @@
OutputManager.hh \
OutputManager.cc \
OutputSolnSubset.hh \
+ OutputSolnPoints.hh \
UCDFaultFile.hh \
VertexFilter.hh \
VertexFilter.cc \
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputManager.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputManager.hh 2011-10-20 23:26:58 UTC (rev 19107)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputManager.hh 2011-10-21 04:56:37 UTC (rev 19108)
@@ -144,7 +144,7 @@
const char* label =0,
const int labelId =0);
-// PROTECTED MEMBERS ////////////////////////////////////////////////////
+// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
/** Dimension field.
@@ -153,6 +153,11 @@
*/
field_type& _dimension(field_type& fieldIn);
+// PROTECTED MEMBERS ////////////////////////////////////////////////////
+protected :
+
+ topology::Fields<field_type>* _fields; ///< Buffer fields.
+
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
@@ -169,8 +174,6 @@
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/pylith/meshio/OutputSolnPoints.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc 2011-10-20 23:26:58 UTC (rev 19107)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc 2011-10-21 04:56:37 UTC (rev 19108)
@@ -21,8 +21,15 @@
#include "OutputSolnPoints.hh" // implementation of class methods
#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "MeshBuilder.hh" // USES MeshBuilder
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+
// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
// Constructor
pylith::meshio::OutputSolnPoints::OutputSolnPoints(void) :
_mesh(0),
@@ -45,10 +52,17 @@
{ // deallocate
OutputManager<topology::Mesh, topology::Field<topology::Mesh> >::deallocate();
+
+ err = DMMeshCreateGlobalScatter(sieveMesh, _section, order, false,
+ &sinfo.scatter);
+
+
if (_interpolator) {
assert(_mesh);
- DMMeshInterpolationDestroy(_mesh->sieveMesh(), _interpolator);
- delete _interpolater; _interpolator = 0;
+ const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
+ PetscErrorCode err =
+ DMMeshInterpolationDestroy(sieveMesh, &_interpolator);
+ PETSC_CHECK_ERROR(err);
} // if
delete _mesh; _mesh = 0;
@@ -58,14 +72,132 @@
// ----------------------------------------------------------------------
// Get mesh associated with points.
const pylith::topology::Mesh&
-pylith::meshio::OutputSolnPoints::createPointsMesh(const PylithScalar* points,
- const int numPoints,
- const int spaceDim)
+pylith::meshio::OutputSolnPoints::pointsMesh(void)
+{ // pointsMesh
+ assert(_pointsMesh);
+ return *_pointsMesh;
+} // pointsMesh
+
+
+// ----------------------------------------------------------------------
+// Setup interpolator.
+void
+pylith::meshio::OutputSolnPoints::setupInterpolator(topology::Mesh* mesh,
+ const PylithScalar* points,
+ const int numPoints,
+ const int spaceDim)
{ // createPointsMesh
- // delete _pointsMesh; _pointsMesh = new topology::Mesh();
+ assert(mesh);
+ assert(points);
+
+ _mesh = mesh;
+
+ // Create mesh without cells for points.
+ const int meshDim = 0;
+ delete _pointsMesh; _pointsMesh = new topology::Mesh(meshDim);
+ _pointsMesh->createSieveMesh(0);
assert(_pointsMesh);
- return *_pointsMesh;
+
+ scalar_array pointsArray(points, numPoints*spaceDim);
+ int_array cells;
+ const int numCells = 0;
+ const int numCorners = 0;
+ const bool interpolate = false;
+ MeshBuilder::buildMesh(_pointsMesh,
+ &pointsArray, numPoints, spaceDim,
+ cells, numCells, numCorners, meshDim,
+ interpolate);
+ _pointsMesh->coordsys(_mesh->coordsys());
+
+ // Setup interpolator object
+ PetscErrorCode err = 0;
+
+ err = DMMeshInterpolationCreate(_mesh->sieveMesh(), &_interpolator);
+ CHECK_PETSC_ERROR(err);
+
+ const spatialdata::geocoords::CoordSys* cs = _pointsMesh().coordsys();
+ assert(0 != cs);
+ assert(cs->spaceDim() == spaceDim);
+
+ err = DMMeshInterpolationSetDim(_mesh->sieveMesh(), spaceDim,
+ _interpolator); CHECK_PETSC_ERROR(err);
+
+ err = DMMeshInterpolationAddPoints(_mesh->sieveMesh(), numPoints, points,
+ _interpolator); CHECK_PETSC_ERROR(err);
+
+ err = DMMeshInterpolationSetUp(_mesh->sieveMesh(), _interpolator);
+ CHECK_PETSC_ERROR(err);
+
} // createPointsMesh
+// ----------------------------------------------------------------------
+// Append finite-element vertex field to file.
+void
+pylith::meshio::OutputSolnPoints::appendVertexField(const PylithScalar t,
+ topology::Field<topology::Mesh>& field,
+ const topology::Mesh& mesh)
+{ // appendVertexField
+ assert(_mesh);
+ assert(_fields);
+
+ const ALE::Obj<SieveMesh>& pointsSieveMesh = _pointsMesh->sieveMesh();
+ assert(!pointsSieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ pointsSieveMesh->depthStratum(0);
+ assert(!vertices.isNull());
+
+ const int fiberDim = (vertices->begin() != vertices->end()) ?
+ field.section()->getFiberDimension(*vertices->begin()) : 0;
+
+ // Create field if necessary for interpolated values or recreate
+ // field if mismatch in size between buffer and field.
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("Output");
+ if (!_fields->hasField("buffer (interpolated)")) {
+ _fields->add("buffer (interpolated)", field.label());
+ } // if
+
+ topology::Field<topology::Mesh>& fieldInterp =
+ _fields->get("buffer (interpolated)");
+ if (vertices->size()*fiberDim != fieldInterp.sectionSize()) {
+ fieldInterp.newSection(vertices, fiberDim);
+ fieldInterp.allocate();
+ } // if
+ logger.stagePop();
+
+ fieldInterp.zero();
+ fieldInterp.vectorFieldType(field.vectorFieldType());
+ fieldInterp.scale(field.scale());
+
+ PetscVec fieldInterpVec = fieldInterp.vector();
+ assert(fieldInterpVec);
+
+ const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
+ assert(!sieveMesh.isNull());
+
+ PetscErrorCode err = 0;
+
+ err = DMMeshInterpolationSetDof(sieveMesh, fiberDim, _interpolator);
+ CHECK_PETSC_ERROR(err);
+
+ err = DMMeshInterpolationEvaluate(sieveMesh, field.section(), fieldInterpVec,
+ _interpolator); CHECK_PETSC_ERROR(err);
+
+ OutputManager<topology::Mesh, topology::Field<topology::Mesh> >::appendVertexField(t, fieldInterp, mesh);
+} // appendVertexField
+
+
+// ----------------------------------------------------------------------
+// Append finite-element cell field to file.
+void
+pylith::meshio::OutputSolnPoints::appendCellField(const PylithScalar t,
+ topology::Field<topology::Mesh>& field,
+ const char* label,
+ const int labelId)
+{ // appendCellField
+ throw std::logic_error("OutputSolnPoints::appendCellField() not implemented.");
+} // appendCellField
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.hh 2011-10-20 23:26:58 UTC (rev 19107)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.hh 2011-10-21 04:56:37 UTC (rev 19108)
@@ -28,13 +28,12 @@
// Include directives ---------------------------------------------------
#include "meshiofwd.hh" // forward declarations
+#include "pylith/utils/petscfwd.h"
-#include "pylith/topology/SubMesh.hh" // ISA OutputManager<Mesh>
+#include "pylith/topology/Mesh.hh" // ISA OutputManager<Mesh>
#include "pylith/topology/Field.hh" // ISA OutputManager<Field<Mesh>>
#include "OutputManager.hh" // ISA OutputManager
-#include <string> // HASA std::string
-
// OutputSolnPoints -----------------------------------------------------
/** @brief C++ object for managing output of finite-element data over
* a subdomain.
@@ -60,27 +59,55 @@
*
* @returns Mesh associated with points.
*/
- const topology::Mesh& createPointsMesh(const PylithScalar* points,
- const int numPoints,
- const int spaceDim);
+ const topology::Mesh& pointsMesh(void);
+
+ /** Setup interpolator.
+ *
+ * @param mesh Domain mesh.
+ * @param points Array of coordinates for points [numPoints*spaceDim].
+ * @param numPoints Number of points.
+ * @param spaceDim Spatial dimension for coordinates.
+ */
+ void setupInterpolator(topology::Mesh* mesh,
+ const PylithScalar* points,
+ const int numPoints,
+ const int spaceDim);
+ /** Append finite-element vertex field to file.
+ *
+ * @param t Time associated with field.
+ * @param field Vertex field.
+ * @param mesh Mesh for output.
+ */
+ void appendVertexField(const PylithScalar t,
+ topology::Field<topology::Mesh>& field,
+ const topology::Mesh& mesh);
+
+ /** Append finite-element cell field to file.
+ *
+ * @param t Time associated with field.
+ * @param field Cell field.
+ * @param label Name of label defining cells to include in output
+ * (=0 means use all cells in mesh).
+ * @param labelId Value of label defining which cells to include.
+ */
+ void appendCellField(const PylithScalar t,
+ topology::Field<topology::Mesh>& field,
+ const char* label =0,
+ const int labelId =0);
+
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
OutputSolnPoints(const OutputSolnPoints&); ///< Not implemented.
const OutputSolnPoints& operator=(const OutputSolnPoints&); ///< Not implemented
-// PRIVATE TYPEDEFS /////////////////////////////////////////////////////
-private :
-
- typedef DMMeshInterpolationInfo PetscDMMeshInterpolationInfo;
-
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
topology::Mesh* _mesh; ///< Domain mesh.
topology::Mesh* _pointsMesh; ///< Mesh for points (no cells).
- PetscDMMeshInterpolationInfo* _interpolator;
+ PetscDMMeshInterpolationInfo _interpolator; ///< Field interpolator.
}; // OutputSolnPoints
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/meshiofwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/meshiofwd.hh 2011-10-20 23:26:58 UTC (rev 19107)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/meshiofwd.hh 2011-10-21 04:56:37 UTC (rev 19108)
@@ -57,6 +57,7 @@
template<typename field_type> class VertexFilter;
template<typename field_type> class VertexFilterVecNorm;
class OutputSolnSubset;
+ class OutputSolnPoints;
class UCDFaultFile;
Modified: short/3D/PyLith/trunk/libsrc/pylith/utils/petscfwd.h
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/utils/petscfwd.h 2011-10-20 23:26:58 UTC (rev 19107)
+++ short/3D/PyLith/trunk/libsrc/pylith/utils/petscfwd.h 2011-10-21 04:56:37 UTC (rev 19108)
@@ -25,6 +25,8 @@
#if !defined(pylith_utils_petscfwd_h)
#define pylith_utils_petscfwd_h
+/// forward declaration for PETSc PetscErrorCode
+typedef int PetscErrorCode;
/// forward declaration for PETSc Mat
typedef struct _p_Mat* PetscMat;
@@ -35,9 +37,6 @@
/// forward declaration for PETSc VecScatter
typedef struct _p_VecScatter* PetscVecScatter;
-/// forward declaration for PETSc ISLocalToGlobalMapping
-typedef struct _p_ISLocalToGlobalMapping* PetscISLocalToGlobalMapping;
-
/// forward declaration for PETSc KSP
typedef struct _p_KSP* PetscKSP;
@@ -47,9 +46,13 @@
/// forward declaration for PETSc PC
typedef struct _p_PC* PetscPC;
-/// forward declaration for PETSc PetscErrorCode
-typedef int PetscErrorCode;
+/// forward declaration for PETSc ISLocalToGlobalMapping
+typedef struct _p_ISLocalToGlobalMapping* PetscISLocalToGlobalMapping;
+/// forward declaration for PETSc DMMeshInterpolationInfo
+typedef struct _DMMeshInterpolationInfo* PetscDMMeshInterpolationInfo;
+
+
#endif // pylith_utils_petscfwd_h
// End of file
Modified: short/3D/PyLith/trunk/pylith/meshio/OutputSolnPoints.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/OutputSolnPoints.py 2011-10-20 23:26:58 UTC (rev 19107)
+++ short/3D/PyLith/trunk/pylith/meshio/OutputSolnPoints.py 2011-10-21 04:56:37 UTC (rev 19108)
@@ -115,7 +115,8 @@
self._eventLogger.eventBegin(logEvent)
points = self.reader.read()
- self.mesh = ModuleOutputSolnPoints.createPointsMesh(self, points)
+ ModuleOutputSolnPoints.setupInterpolator(self, mesh, points)
+ self.mesh = ModuleOutputSolnPoints.createPointsMesh(self)
OutputManager.initialize(self, normalizer)
self._eventLogger.eventEnd(logEvent)
More information about the CIG-COMMITS
mailing list