[cig-commits] r19754 - in short/3D/PyLith/trunk: examples/3d/hex8 libsrc/pylith/materials libsrc/pylith/meshio modulesrc/meshio pylith/meshio
brad at geodynamics.org
brad at geodynamics.org
Thu Mar 8 17:22:29 PST 2012
Author: brad
Date: 2012-03-08 17:22:28 -0800 (Thu, 08 Mar 2012)
New Revision: 19754
Added:
short/3D/PyLith/trunk/examples/3d/hex8/output_points.txt
short/3D/PyLith/trunk/modulesrc/meshio/OutputSolnPoints.i
Modified:
short/3D/PyLith/trunk/examples/3d/hex8/step19.cfg
short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc
short/3D/PyLith/trunk/modulesrc/meshio/Makefile.am
short/3D/PyLith/trunk/modulesrc/meshio/meshio.i
short/3D/PyLith/trunk/pylith/meshio/OutputSolnPoints.py
Log:
Fixed small bugs in Python implementation of OutputSolnPoints.
Added: short/3D/PyLith/trunk/examples/3d/hex8/output_points.txt
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/output_points.txt (rev 0)
+++ short/3D/PyLith/trunk/examples/3d/hex8/output_points.txt 2012-03-09 01:22:28 UTC (rev 19754)
@@ -0,0 +1,7 @@
+# This file contains a list of points where we want the solution. The
+# solution will be interpolated to these points, so they need not
+# coincide with vertices of the mesh.
+10.0 10.0 -500.0
+10.0 10.0 -1500.0
+10.0 10.0 -2500.0
+10.0 10.0 -3500.0
Modified: short/3D/PyLith/trunk/examples/3d/hex8/step19.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/step19.cfg 2012-03-08 22:07:32 UTC (rev 19753)
+++ short/3D/PyLith/trunk/examples/3d/hex8/step19.cfg 2012-03-09 01:22:28 UTC (rev 19754)
@@ -36,11 +36,14 @@
[pylithapp.timedependent.implicit]
# Set the output to an array of 2 output managers.
# We will output the solution over the domain and the ground surface.
-output = [domain,subdomain]
+output = [domain,subdomain,points]
# Set subdomain component to OutputSolnSubset (subset of domain).
output.subdomain = pylith.meshio.OutputSolnSubset
+# Set points component to OutputSolnPoints (arbitrary points in domain).
+output.points = pylith.meshio.OutputSolnPoints
+
# Change the total simulation time to 700 years, and set the time
# step to 10 years.
[pylithapp.timedependent.implicit.time_step]
@@ -142,6 +145,13 @@
writer.time_constant = 1.0*year
writer.time_format = %04.0f
+# Give basename for VTK domain output of solution at points.
+[pylithapp.problem.formulation.output.points]
+reader.filename = output_points.txt
+writer.filename = output/step19-points.vtk
+writer.time_constant = 1.0*year
+writer.time_format = %04.0f
+
# Give basename for VTK output of upper_crust state variables.
[pylithapp.timedependent.materials.upper_crust.output]
cell_filter = pylith.meshio.CellFilterAvgMesh
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc 2012-03-08 22:07:32 UTC (rev 19753)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc 2012-03-09 01:22:28 UTC (rev 19754)
@@ -906,9 +906,10 @@
const PylithScalar dFac = 1.0/(sqrt(2.0) * ae);
const PylithScalar testMult = plasticFac *
(meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
- const bool tensileYield = (sqrt(2.0) * d < testMult) ? true: false;
- const PylithScalar plasticMult = tensileYield ? sqrt(2.0) * d: testMult;
+ const bool tensileYield = (sqrt(2.0)*d < testMult) ? true : false;
+ const PylithScalar plasticMult = (tensileYield) ? sqrt(2.0)*d : testMult;
+
// Define some constants, vectors, and matrices.
const PylithScalar third = 1.0/3.0;
const PylithScalar dEdEpsilon[6][6] = {
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc 2012-03-08 22:07:32 UTC (rev 19753)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnPoints.cc 2012-03-09 01:22:28 UTC (rev 19754)
@@ -101,6 +101,16 @@
assert(csMesh);
assert(csMesh->spaceDim() == spaceDim);
+#if 1 // DEBUGGING
+ std::cout << "OUTPUT SOLN POINTS" << std::endl;
+ for (int i=0; i < numPoints; ++i) {
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ std::cout << " " << points[i*spaceDim+iDim];
+ } // for
+ std::cout << "\n";
+ } // for
+#endif
+
scalar_array pointsArray(points, numPoints*spaceDim);
int_array cells(numPoints);
for (int i=0; i < numPoints; ++i)
Modified: short/3D/PyLith/trunk/modulesrc/meshio/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/modulesrc/meshio/Makefile.am 2012-03-08 22:07:32 UTC (rev 19753)
+++ short/3D/PyLith/trunk/modulesrc/meshio/Makefile.am 2012-03-09 01:22:28 UTC (rev 19754)
@@ -37,7 +37,8 @@
DataWriter.i \
DataWriterVTK.i \
OutputManager.i \
- OutputSolnSubset.i
+ OutputSolnSubset.i \
+ OutputSolnPoints.i
swig_generated = \
meshio_wrap.cxx \
Added: short/3D/PyLith/trunk/modulesrc/meshio/OutputSolnPoints.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/meshio/OutputSolnPoints.i (rev 0)
+++ short/3D/PyLith/trunk/modulesrc/meshio/OutputSolnPoints.i 2012-03-09 01:22:28 UTC (rev 19754)
@@ -0,0 +1,125 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2012 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+/**
+ * @file modulesrc/meshio/OutputSolnPoints.i
+ *
+ * @brief Python interface to C++ OutputSolnPoints object.
+ */
+
+namespace pylith {
+ namespace meshio {
+
+%template(_MeshOutputManager) pylith::meshio::OutputManager<pylith::topology::Mesh, pylith::topology::Field<pylith::topology::Mesh> >;
+
+
+ class pylith::meshio::OutputSolnPoints : public OutputManager<pylith::topology::Mesh, pylith::topology::Field<pylith::topology::Mesh> >
+ { // OutputSolnPoints
+
+ // PUBLIC METHODS ///////////////////////////////////////////////////
+ public :
+
+ /// Constructor
+ OutputSolnPoints(void);
+
+ /// Destructor
+ ~OutputSolnPoints(void);
+
+ /// Deallocate PETSc and local data structures.
+ void deallocate(void);
+
+ /** Get mesh associated with points.
+ *
+ * @returns Mesh associated with points.
+ */
+ const pylith::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.
+ */
+ %apply(PylithScalar* IN_ARRAY2, int DIM1, int DIM2) {
+ (const PylithScalar* points,
+ const int numPoints,
+ const int spaceDim)
+ };
+ void setupInterpolator(pylith::topology::Mesh* mesh,
+ const PylithScalar* points,
+ const int numPoints,
+ const int spaceDim);
+ %clear(const PylithScalar* points, const int numPoints, const int spaceDim);
+
+ /** Prepare for output.
+ *
+ * @param mesh Finite-element mesh object.
+ * @param numTimeSteps Expected number of time steps.
+ * @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 open(const pylith::topology::Mesh& mesh,
+ const int numTimeSteps,
+ const char* label =0,
+ const int labelId =0);
+
+ /** Setup file for writing fields at time step.
+ *
+ * @param t Time of time step.
+ * @param mesh Finite-element mesh object.
+ * @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 openTimeStep(const PylithScalar t,
+ const pylith::topology::Mesh& mesh,
+ const char* label =0,
+ const int labelId =0);
+
+ /** 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,
+ pylith::topology::Field<topology::Mesh>& field,
+ const pylith::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,
+ pylith::topology::Field<topology::Mesh>& field,
+ const char* label =0,
+ const int labelId =0);
+
+ }; // OutputSolnPoints
+
+ } // meshio
+} // pylith
+
+
+// End of file
Modified: short/3D/PyLith/trunk/modulesrc/meshio/meshio.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/meshio/meshio.i 2012-03-08 22:07:32 UTC (rev 19753)
+++ short/3D/PyLith/trunk/modulesrc/meshio/meshio.i 2012-03-09 01:22:28 UTC (rev 19754)
@@ -37,6 +37,7 @@
#include "pylith/meshio/DataWriterVTK.hh"
#include "pylith/meshio/OutputManager.hh"
#include "pylith/meshio/OutputSolnSubset.hh"
+#include "pylith/meshio/OutputSolnPoints.hh"
#if defined(ENABLE_HDF5)
#include "pylith/meshio/DataWriterHDF5.hh"
#include "pylith/meshio/DataWriterHDF5Ext.hh"
@@ -60,6 +61,15 @@
%include "typemaps.i"
%include "../include/scalartypemaps.i"
+// Numpy interface stuff
+%{
+#define SWIG_FILE_WITH_INIT
+%}
+%include "../include/numpy.i"
+%init %{
+import_array();
+%}
+
// Interfaces
%include "MeshIOObj.i"
%include "MeshIOAscii.i"
@@ -77,6 +87,7 @@
%include "DataWriterVTK.i"
%include "OutputManager.i"
%include "OutputSolnSubset.i"
+%include "OutputSolnPoints.i"
#if defined(ENABLE_HDF5)
%include "DataWriterHDF5.i"
%include "DataWriterHDF5Ext.i"
@@ -101,6 +112,7 @@
%template(MeshDataWriterVTK) pylith::meshio::DataWriterVTK<pylith::topology::Mesh, pylith::topology::Field<pylith::topology::Mesh> >;
%template(SubMeshDataWriterVTK) pylith::meshio::DataWriterVTK<pylith::topology::SubMesh, pylith::topology::Field<pylith::topology::Mesh> >;
%template(SubSubMeshDataWriterVTK) pylith::meshio::DataWriterVTK<pylith::topology::SubMesh, pylith::topology::Field<pylith::topology::SubMesh> >;
+%template(PointsDataWriterVTK) pylith::meshio::DataWriterVTK<pylith::topology::Mesh, pylith::topology::Field<pylith::topology::Mesh> >;
#if defined(ENABLE_HDF5)
%template(MeshDataWriterHDF5) pylith::meshio::DataWriterHDF5<pylith::topology::Mesh, pylith::topology::Field<pylith::topology::Mesh> >;
Modified: short/3D/PyLith/trunk/pylith/meshio/OutputSolnPoints.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/OutputSolnPoints.py 2012-03-08 22:07:32 UTC (rev 19753)
+++ short/3D/PyLith/trunk/pylith/meshio/OutputSolnPoints.py 2012-03-09 01:22:28 UTC (rev 19754)
@@ -98,15 +98,6 @@
return
- def verifyConfiguration(self, mesh):
- """
- Verify compatibility of configuration.
- """
- OutputManager.verifyConfiguration(self, mesh)
- ModuleOutputSolnPoints.verifyConfiguration(self, mesh)
- return
-
-
def initialize(self, mesh, normalizer):
"""
Initialize output manager.
@@ -114,10 +105,17 @@
logEvent = "%sinit" % self._loggingPrefix
self._eventLogger.eventBegin(logEvent)
+ OutputManager.initialize(self, normalizer)
+
+ # Read points
points = self.reader.read()
+
+ # Convert to mesh coordinate system
+ from spatialdata.geocoords.Converter import convert
+ convert(points, mesh.coordsys(), self.coordsys)
+
ModuleOutputSolnPoints.setupInterpolator(self, mesh, points)
self.mesh = ModuleOutputSolnPoints.createPointsMesh(self)
- OutputManager.initialize(self, normalizer)
self._eventLogger.eventEnd(logEvent)
return
@@ -153,8 +151,6 @@
"""
try:
OutputManager._configure(self)
- ModuleOutputSolnPoints.label(self, self.label)
- ModuleOutputSolnPoints.coordsys(self, self.inventory.coordsys)
ModuleOutputSolnPoints.writer(self, self.inventory.writer)
from pylith.utils.NullComponent import NullComponent
if not isinstance(self.inventory.vertexFilter, NullComponent):
More information about the CIG-COMMITS
mailing list