[cig-commits] r20620 - in short/3D/PyLith/branches/v1.7-trunk: libsrc/pylith/feassemble libsrc/pylith/materials unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Wed Aug 22 12:06:21 PDT 2012
Author: brad
Date: 2012-08-22 12:06:20 -0700 (Wed, 22 Aug 2012)
New Revision: 20620
Modified:
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/CellGeometry.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicit.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryHex3D.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryHex3D.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine1D.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine1D.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine2D.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine2D.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine3D.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine3D.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint1D.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint1D.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint2D.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint2D.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint3D.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint3D.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad2D.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad3D.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad3D.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri2D.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri2D.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri3D.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri3D.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/Quadrature.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/QuadratureRefCell.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/QuadratureRefCell.icc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/materials/ElasticMaterial.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/materials/ElasticMaterial.hh
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/TestIntegrator.cc
Log:
Started work on calculating stable time step for explicit time stepping. Created routines to get smallest cell width.
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/CellGeometry.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/CellGeometry.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/CellGeometry.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -151,6 +151,14 @@
const int dim,
const int npts) const = 0;
+ /** Compute minimum width across cell.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Minimum width across cell.
+ */
+ virtual
+ PylithScalar minCellWidth(const scalar_array& coordinatesCell) const = 0;
+
/** Compute orientation of cell at location.
*
* The orientation is returned as an array of direction cosines
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -83,6 +83,15 @@
} // timeStep
// ----------------------------------------------------------------------
+// Get stable time step for advancing from time t to time t+dt.
+PylithScalar
+pylith::feassemble::ElasticityExplicit::stableTimeStep(const topology::Mesh& mesh) const
+{ // stableTimeStep
+ assert(_material);
+ return _material->stableTimeStepExplicit(mesh, _quadrature);
+} // stableTimeStep
+
+// ----------------------------------------------------------------------
// Set normalized viscosity for numerical damping.
void
pylith::feassemble::ElasticityExplicit::normViscosity(const PylithScalar viscosity)
@@ -108,10 +117,10 @@
typedef void (pylith::feassemble::ElasticityExplicit::*elasticityResidual_fn_type)
(const scalar_array&);
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(0 != _logger);
- assert(0 != fields);
+ assert(_quadrature);
+ assert(_material);
+ assert(_logger);
+ assert(fields);
const int setupEvent = _logger->eventId("ElIR setup");
const int geometryEvent = _logger->eventId("ElIR geometry");
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicit.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicit.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicit.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -90,6 +90,15 @@
*/
void timeStep(const PylithScalar dt);
+ /** Get stable time step for advancing from time t to time t+dt.
+ *
+ * Default is current time step.
+ *
+ * @param mesh Finite-element mesh.
+ * @returns Time step
+ */
+ PylithScalar stableTimeStep(const topology::Mesh& mesh) const;
+
/** Set normalized viscosity for numerical damping.
*
* @param viscosity Normalized viscosity (viscosity / elastic modulus).
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -81,6 +81,15 @@
} // timeStep
// ----------------------------------------------------------------------
+// Get stable time step for advancing from time t to time t+dt.
+PylithScalar
+pylith::feassemble::ElasticityExplicitLgDeform::stableTimeStep(const topology::Mesh& mesh) const
+{ // stableTimeStep
+ assert(_material);
+ return _material->stableTimeStepExplicit(mesh, _quadrature);
+} // stableTimeStep
+
+// ----------------------------------------------------------------------
// Set normalized viscosity for numerical damping.
void
pylith::feassemble::ElasticityExplicitLgDeform::normViscosity(const PylithScalar viscosity)
@@ -106,10 +115,10 @@
typedef void (pylith::feassemble::ElasticityExplicitLgDeform::*elasticityResidual_fn_type)
(const scalar_array&, const scalar_array&);
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(0 != _logger);
- assert(0 != fields);
+ assert(_quadrature);
+ assert(_material);
+ assert(_logger);
+ assert(fields);
const int setupEvent = _logger->eventId("ElIR setup");
const int geometryEvent = _logger->eventId("ElIR geometry");
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -92,6 +92,15 @@
*/
void timeStep(const PylithScalar dt);
+ /** Get stable time step for advancing from time t to time t+dt.
+ *
+ * Default is current time step.
+ *
+ * @param mesh Finite-element mesh.
+ * @returns Time step
+ */
+ PylithScalar stableTimeStep(const topology::Mesh& mesh) const;
+
/** Set normalized viscosity for numerical damping.
*
* @param viscosity Normalized viscosity (viscosity / elastic modulus).
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryHex3D.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryHex3D.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryHex3D.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -25,6 +25,7 @@
#include "petsc.h" // USES PetscLogFlops
#include "pylith/utils/array.hh" // USES scalar_array
+#include "pylith/utils/constdefs.h" // USES PYLITH_MAXSCALAR
#include <cassert> // USES assert()
@@ -366,5 +367,43 @@
PetscLogFlops(78 + npts*69);
} // jacobian
+// ----------------------------------------------------------------------
+// Compute minimum width across cell.
+PylithScalar
+pylith::feassemble::GeometryHex3D::minCellWidth(const scalar_array& coordinatesCell) const
+{ // minCellWidth
+ const int numCorners = 8;
+ const int spaceDim = 3;
+ assert(numCorners*spaceDim == coordinatesCell.size());
+ const int numEdges = 12;
+ const int edges[numEdges][2] = {
+ {0, 1}, {1, 2}, {2, 3}, {3, 0},
+ {4, 5}, {5, 6}, {6, 7}, {7, 0},
+ {0, 4}, {1, 5}, {2, 6}, {3, 7},
+ };
+
+ PylithScalar minWidth = PYLITH_MAXSCALAR;
+ for (int iedge=0; iedge < numEdges; ++iedge) {
+ const int iA = edges[iedge][0];
+ const int iB = edges[iedge][1];
+ const PylithScalar xA = coordinatesCell[spaceDim*iA ];
+ const PylithScalar yA = coordinatesCell[spaceDim*iA+1];
+ const PylithScalar zA = coordinatesCell[spaceDim*iA+2];
+ const PylithScalar xB = coordinatesCell[spaceDim*iB ];
+ const PylithScalar yB = coordinatesCell[spaceDim*iB+1];
+ const PylithScalar zB = coordinatesCell[spaceDim*iB+2];
+
+ const PylithScalar edgeLen = sqrt(pow(xB-xA,2) + pow(yB-yA,2) + pow(zB-zA,2));
+ if (edgeLen < minWidth) {
+ minWidth = edgeLen;
+ } // if
+ } // for
+
+ PetscLogFlops(numEdges*9);
+
+ return minWidth;
+} // minCellWidth
+
+
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryHex3D.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryHex3D.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryHex3D.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -114,6 +114,13 @@
const int dim,
const int npts =1) const;
+ /** Compute minimum width across cell.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Minimum width across cell.
+ */
+ PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine1D.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine1D.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine1D.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -23,6 +23,7 @@
#include "GeometryPoint1D.hh" // USES GeometryPoint
#include "pylith/utils/array.hh" // USES scalar_array
+#include "pylith/utils/constdefs.h" // USES scalar_array
#include "petsc.h" // USES PetscLogFlops
@@ -140,4 +141,25 @@
} // jacobian
+// ----------------------------------------------------------------------
+// Compute minimum width across cell.
+PylithScalar
+pylith::feassemble::GeometryLine1D::minCellWidth(const scalar_array& coordinatesCell) const
+{ // minCellWidth
+ const int numCorners = 2;
+ const int spaceDim = 1;
+ assert(2*spaceDim == coordinatesCell.size() ||
+ 3*spaceDim == coordinatesCell.size()); // :KLUDGE: allow quadratic
+
+ const PylithScalar xA = coordinatesCell[0];
+ const PylithScalar xB = coordinatesCell[1];
+
+ const PylithScalar minWidth = fabs(xB-xA);
+
+ PetscLogFlops(2);
+
+ return minWidth;
+} // minCellWidth
+
+
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine1D.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine1D.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine1D.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -98,6 +98,13 @@
const int dim,
const int npts =1) const;
+ /** Compute minimum width across cell.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Minimum width across cell.
+ */
+ PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine2D.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine2D.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine2D.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -23,6 +23,7 @@
#include "GeometryPoint2D.hh" // USES GeometryPoint
#include "pylith/utils/array.hh" // USES scalar_array
+#include "pylith/utils/constdefs.h" // USES scalar_array
#include "petsc.h" // USES PetscLogFlops
@@ -161,4 +162,26 @@
} // jacobian
+// ----------------------------------------------------------------------
+// Compute minimum width across cell.
+PylithScalar
+pylith::feassemble::GeometryLine2D::minCellWidth(const scalar_array& coordinatesCell) const
+{ // minCellWidth
+ const int numCorners = 2;
+ const int spaceDim = 2;
+ assert(numCorners*spaceDim == coordinatesCell.size());
+
+ const PylithScalar xA = coordinatesCell[0];
+ const PylithScalar yA = coordinatesCell[1];
+ const PylithScalar xB = coordinatesCell[2];
+ const PylithScalar yB = coordinatesCell[3];
+
+ const PylithScalar minWidth = sqrt(pow(xB-xA,2) + pow(yB-yA,2));
+
+ PetscLogFlops(6);
+
+ return minWidth;
+} // minCellWidth
+
+
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine2D.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine2D.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine2D.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -98,6 +98,13 @@
const int dim,
const int npts =1) const;
+ /** Compute minimum width across cell.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Minimum width across cell.
+ */
+ PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine3D.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine3D.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine3D.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -25,6 +25,7 @@
#include "petsc.h" // USES PetscLogFlops
#include "pylith/utils/array.hh" // USES scalar_array
+#include "pylith/utils/constdefs.h" // USES scalar_array
#include <cassert> // USES assert()
@@ -172,4 +173,28 @@
} // jacobian
+// ----------------------------------------------------------------------
+// Compute minimum width across cell.
+PylithScalar
+pylith::feassemble::GeometryLine3D::minCellWidth(const scalar_array& coordinatesCell) const
+{ // minCellWidth
+ const int numCorners = 2;
+ const int spaceDim = 3;
+ assert(numCorners*spaceDim == coordinatesCell.size());
+
+ const PylithScalar xA = coordinatesCell[0];
+ const PylithScalar yA = coordinatesCell[1];
+ const PylithScalar zA = coordinatesCell[2];
+ const PylithScalar xB = coordinatesCell[3];
+ const PylithScalar yB = coordinatesCell[4];
+ const PylithScalar zB = coordinatesCell[5];
+
+ const PylithScalar minWidth = sqrt(pow(xB-xA,2) + pow(yB-yA,2) + pow(zB-zA,2));
+
+ PetscLogFlops(9);
+
+ return minWidth;
+} // minCellWidth
+
+
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine3D.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine3D.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryLine3D.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -98,6 +98,13 @@
const int dim,
const int npts =1) const;
+ /** Compute minimum width across cell.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Minimum width across cell.
+ */
+ PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint1D.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint1D.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint1D.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -21,6 +21,7 @@
#include "GeometryPoint1D.hh" // implementation of class methods
#include "pylith/utils/array.hh" // USES scalar_array
+#include "pylith/utils/constdefs.h" // USES scalar_array
#include <cassert> // USES assert()
@@ -115,4 +116,13 @@
} // jacobian
+// ----------------------------------------------------------------------
+// Compute minimum width across cell.
+PylithScalar
+pylith::feassemble::GeometryPoint1D::minCellWidth(const scalar_array& coordinatesCell) const
+{ // minCellWidth
+ return PYLITH_MAXSCALAR;
+} // minCellWidth
+
+
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint1D.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint1D.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint1D.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -98,6 +98,13 @@
const int dim,
const int npts =1) const;
+ /** Compute minimum width across cell.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Minimum width across cell.
+ */
+ PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint2D.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint2D.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint2D.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -21,6 +21,7 @@
#include "GeometryPoint2D.hh" // implementation of class methods
#include "pylith/utils/array.hh" // USES scalar_array
+#include "pylith/utils/constdefs.h" // USES scalar_array
#include <cassert> // USES assert()
@@ -116,4 +117,13 @@
} // jacobian
+// ----------------------------------------------------------------------
+// Compute minimum width across cell.
+PylithScalar
+pylith::feassemble::GeometryPoint2D::minCellWidth(const scalar_array& coordinatesCell) const
+{ // minCellWidth
+ return PYLITH_MAXSCALAR;
+} // minCellWidth
+
+
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint2D.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint2D.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint2D.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -98,6 +98,13 @@
const int dim,
const int npts =1) const;
+ /** Compute minimum width across cell.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Minimum width across cell.
+ */
+ PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint3D.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint3D.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint3D.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -21,6 +21,7 @@
#include "GeometryPoint3D.hh" // implementation of class methods
#include "pylith/utils/array.hh" // USES scalar_array
+#include "pylith/utils/constdefs.h" // USES scalar_array
#include <cassert> // USES assert()
@@ -116,4 +117,13 @@
} // jacobian
+// ----------------------------------------------------------------------
+// Compute minimum width across cell.
+PylithScalar
+pylith::feassemble::GeometryPoint3D::minCellWidth(const scalar_array& coordinatesCell) const
+{ // minCellWidth
+ return PYLITH_MAXSCALAR;
+} // minCellWidth
+
+
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint3D.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint3D.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryPoint3D.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -98,6 +98,13 @@
const int dim,
const int npts =1) const;
+ /** Compute minimum width across cell.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Minimum width across cell.
+ */
+ PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -25,6 +25,7 @@
#include "petsc.h" // USES PetscLogFlops
#include "pylith/utils/array.hh" // USES scalar_array
+#include "pylith/utils/constdefs.h" // USES scalar_array
#include <cassert> // USES assert()
@@ -217,4 +218,39 @@
} // jacobian
+// ----------------------------------------------------------------------
+// Compute minimum width across cell.
+PylithScalar
+pylith::feassemble::GeometryQuad2D::minCellWidth(const scalar_array& coordinatesCell) const
+{ // minCellWidth
+ const int numCorners = 4;
+ const int spaceDim = 2;
+ assert(numCorners*spaceDim == coordinatesCell.size());
+
+ const int numEdges = 4;
+ const int edges[numEdges][2] = {
+ {0, 1}, {1, 2}, {2, 3}, {3, 0},
+ };
+
+ PylithScalar minWidth = PYLITH_MAXSCALAR;
+ for (int iedge=0; iedge < numEdges; ++iedge) {
+ const int iA = edges[iedge][0];
+ const int iB = edges[iedge][1];
+ const PylithScalar xA = coordinatesCell[spaceDim*iA ];
+ const PylithScalar yA = coordinatesCell[spaceDim*iA+1];
+ const PylithScalar xB = coordinatesCell[spaceDim*iB ];
+ const PylithScalar yB = coordinatesCell[spaceDim*iB+1];
+
+ const PylithScalar edgeLen = sqrt(pow(xB-xA,2) + pow(yB-yA,2));
+ if (edgeLen < minWidth) {
+ minWidth = edgeLen;
+ } // if
+ } // for
+
+ PetscLogFlops(numEdges*6);
+
+ return minWidth;
+} // minCellWidth
+
+
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad2D.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad2D.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad2D.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -116,6 +116,13 @@
const int dim,
const int npts =1) const;
+ /** Compute minimum width across cell.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Minimum width across cell.
+ */
+ PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad3D.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad3D.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad3D.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -25,6 +25,7 @@
#include "petsc.h" // USES PetscLogFlops
#include "pylith/utils/array.hh" // USES scalar_array
+#include "pylith/utils/constdefs.h" // USES scalar_array
#include <cassert> // USES assert()
@@ -262,4 +263,41 @@
} // jacobian
+// ----------------------------------------------------------------------
+// Compute minimum width across cell.
+PylithScalar
+pylith::feassemble::GeometryQuad3D::minCellWidth(const scalar_array& coordinatesCell) const
+{ // minCellWidth
+ const int numCorners = 4;
+ const int spaceDim = 3;
+ assert(numCorners*spaceDim == coordinatesCell.size());
+
+ const int numEdges = 4;
+ const int edges[numEdges][2] = {
+ {0, 1}, {1, 2}, {2, 3}, {3, 0},
+ };
+
+ PylithScalar minWidth = PYLITH_MAXSCALAR;
+ for (int iedge=0; iedge < numEdges; ++iedge) {
+ const int iA = edges[iedge][0];
+ const int iB = edges[iedge][1];
+ const PylithScalar xA = coordinatesCell[spaceDim*iA ];
+ const PylithScalar yA = coordinatesCell[spaceDim*iA+1];
+ const PylithScalar zA = coordinatesCell[spaceDim*iA+2];
+ const PylithScalar xB = coordinatesCell[spaceDim*iB ];
+ const PylithScalar yB = coordinatesCell[spaceDim*iB+1];
+ const PylithScalar zB = coordinatesCell[spaceDim*iB+2];
+
+ const PylithScalar edgeLen = sqrt(pow(xB-xA,2) + pow(yB-yA,2) + pow(zB-zA,2));
+ if (edgeLen < minWidth) {
+ minWidth = edgeLen;
+ } // if
+ } // for
+
+ PetscLogFlops(numEdges*9);
+
+ return minWidth;
+} // minCellWidth
+
+
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad3D.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad3D.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryQuad3D.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -114,6 +114,13 @@
const int dim,
const int npts =1) const;
+ /** Compute minimum width across cell.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Minimum width across cell.
+ */
+ PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -25,6 +25,7 @@
#include "petsc.h" // USES PetscLogFlops
#include "pylith/utils/array.hh" // USES scalar_array
+#include "pylith/utils/constdefs.h" // USES scalar_array
#include <cassert> // USES assert()
@@ -241,4 +242,46 @@
} // jacobian
+// ----------------------------------------------------------------------
+// Compute minimum width across cell.
+PylithScalar
+pylith::feassemble::GeometryTet3D::minCellWidth(const scalar_array& coordinatesCell) const
+{ // minCellWidth
+ const int numCorners = 4;
+ const int spaceDim = 3;
+ assert(4*spaceDim == coordinatesCell.size() ||
+ 10*spaceDim == coordinatesCell.size()); // :KLUDGE: allow quadratic
+
+ const int numEdges = 6;
+ const int edges[numEdges][2] = {
+ {0, 1}, {1, 2}, {2, 0},
+ {0, 3}, {1, 3}, {2, 3},
+ };
+
+ PylithScalar minWidth = PYLITH_MAXSCALAR;
+ for (int iedge=0; iedge < numEdges; ++iedge) {
+ const int iA = edges[iedge][0];
+ const int iB = edges[iedge][1];
+ const PylithScalar xA = coordinatesCell[spaceDim*iA ];
+ const PylithScalar yA = coordinatesCell[spaceDim*iA+1];
+ const PylithScalar zA = coordinatesCell[spaceDim*iA+2];
+ const PylithScalar xB = coordinatesCell[spaceDim*iB ];
+ const PylithScalar yB = coordinatesCell[spaceDim*iB+1];
+ const PylithScalar zB = coordinatesCell[spaceDim*iB+2];
+
+ const PylithScalar edgeLen = sqrt(pow(xB-xA,2) + pow(yB-yA,2) + pow(zB-zA,2));
+ if (edgeLen < minWidth) {
+ minWidth = edgeLen;
+ } // if
+ } // for
+
+ PetscLogFlops(numEdges*9);
+
+ // TODO: Add radius of inscribed sphere
+ // r = V / 3*A
+
+ return minWidth;
+} // minCellWidth
+
+
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTet3D.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -109,6 +109,13 @@
const int dim,
const int npts =1) const;
+ /** Compute minimum width across cell.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Minimum width across cell.
+ */
+ PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri2D.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri2D.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri2D.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -25,6 +25,7 @@
#include "petsc.h" // USES PetscLogFlops
#include "pylith/utils/array.hh" // USES scalar_array
+#include "pylith/utils/constdefs.h" // USES scalar_array
#include <cassert> // USES assert()
@@ -107,9 +108,9 @@
// Compute Jacobian at location in cell.
void
pylith::feassemble::GeometryTri2D::jacobian(scalar_array* jacobian,
- PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const
+ PylithScalar* det,
+ const scalar_array& vertices,
+ const scalar_array& location) const
{ // jacobian
assert(0 != jacobian);
@@ -182,4 +183,61 @@
} // jacobian
+// ----------------------------------------------------------------------
+// Compute minimum width across cell.
+PylithScalar
+pylith::feassemble::GeometryTri2D::minCellWidth(const scalar_array& coordinatesCell) const
+{ // minCellWidth
+ const int numCorners = 3;
+ const int spaceDim = 2;
+ assert(3*spaceDim == coordinatesCell.size() ||
+ 6*spaceDim == coordinatesCell.size()); // :KLUDGE: allow quadratic
+
+ const int numEdges = 3;
+ const int edges[numEdges][2] = {
+ {0, 1}, {1, 2}, {2, 0},
+ };
+
+ PylithScalar minWidth = PYLITH_MAXSCALAR;
+ for (int iedge=0; iedge < numEdges; ++iedge) {
+ const int iA = edges[iedge][0];
+ const int iB = edges[iedge][1];
+ const PylithScalar xA = coordinatesCell[spaceDim*iA ];
+ const PylithScalar yA = coordinatesCell[spaceDim*iA+1];
+ const PylithScalar xB = coordinatesCell[spaceDim*iB ];
+ const PylithScalar yB = coordinatesCell[spaceDim*iB+1];
+
+ const PylithScalar edgeLen = sqrt(pow(xB-xA,2) + pow(yB-yA,2));
+ if (edgeLen < minWidth) {
+ minWidth = edgeLen;
+ } // if
+ } // for
+
+ PetscLogFlops(numEdges*6);
+
+#if 1
+ // Ad-hoc to account for distorted cells.
+ // Radius of inscribed circle.
+ const PylithScalar xA = coordinatesCell[0];
+ const PylithScalar yA = coordinatesCell[1];
+ const PylithScalar xB = coordinatesCell[2];
+ const PylithScalar yB = coordinatesCell[3];
+ const PylithScalar xC = coordinatesCell[4];
+ const PylithScalar yC = coordinatesCell[5];
+
+ const PylithScalar c = sqrt(pow(xB-xA,2) + pow(yB-yA, 2));
+ const PylithScalar a = sqrt(pow(xC-xB,2) + pow(yC-yB, 2));
+ const PylithScalar b = sqrt(pow(xA-xC,2) + pow(yA-yC, 2));
+
+ const PylithScalar k = (a + b + c) / 3.0;
+ const PylithScalar r = sqrt(k*(k-a)*(k-b)*(k-c)) / k;
+ minWidth = r;
+
+ PetscLogFlops(3*6 + 3 + 8);
+#endif
+
+ return minWidth;
+} // minCellWidth
+
+
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri2D.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri2D.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri2D.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -115,6 +115,13 @@
const int dim,
const int npts =1) const;
+ /** Compute minimum width across cell.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Minimum width across cell.
+ */
+ PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri3D.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri3D.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri3D.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -25,6 +25,7 @@
#include "petsc.h" // USES PetscLogFlops
#include "pylith/utils/array.hh" // USES scalar_array
+#include "pylith/utils/constdefs.h" // USES scalar_array
#include <cassert> // USES assert()
@@ -218,4 +219,64 @@
} // jacobian
+// ----------------------------------------------------------------------
+// Compute minimum width across cell.
+PylithScalar
+pylith::feassemble::GeometryTri3D::minCellWidth(const scalar_array& coordinatesCell) const
+{ // minCellWidth
+ const int numCorners = 2;
+ const int spaceDim = 3;
+ assert(numCorners*spaceDim == coordinatesCell.size());
+
+ const int numEdges = 3;
+ const int edges[numEdges][2] = {
+ {0, 1}, {1, 2}, {2, 0},
+ };
+
+ PylithScalar minWidth = PYLITH_MAXSCALAR;
+ for (int iedge=0; iedge < numEdges; ++iedge) {
+ const int iA = edges[iedge][0];
+ const int iB = edges[iedge][1];
+ const PylithScalar xA = coordinatesCell[spaceDim*iA ];
+ const PylithScalar yA = coordinatesCell[spaceDim*iA+1];
+ const PylithScalar zA = coordinatesCell[spaceDim*iA+2];
+ const PylithScalar xB = coordinatesCell[spaceDim*iB ];
+ const PylithScalar yB = coordinatesCell[spaceDim*iB+1];
+ const PylithScalar zB = coordinatesCell[spaceDim*iB+2];
+
+ const PylithScalar edgeLen = sqrt(pow(xB-xA,2) + pow(yB-yA,2) + pow(zB-zA,2));
+ if (edgeLen < minWidth) {
+ minWidth = edgeLen;
+ } // if
+ } // for
+
+ PetscLogFlops(numEdges*9);
+
+#if 1
+ // Ad-hoc to account for distorted cells.
+ // Radius of inscribed circle.
+ const PylithScalar xA = coordinatesCell[0];
+ const PylithScalar yA = coordinatesCell[1];
+ const PylithScalar zA = coordinatesCell[2];
+ const PylithScalar xB = coordinatesCell[3];
+ const PylithScalar yB = coordinatesCell[4];
+ const PylithScalar zB = coordinatesCell[5];
+ const PylithScalar xC = coordinatesCell[6];
+ const PylithScalar yC = coordinatesCell[7];
+ const PylithScalar zC = coordinatesCell[8];
+
+ const PylithScalar c = sqrt(pow(xB-xA,2) + pow(yB-yA, 2) + pow(zB-zA, 2));
+ const PylithScalar a = sqrt(pow(xC-xB,2) + pow(yC-yB, 2) + pow(zC-zB, 2));
+ const PylithScalar b = sqrt(pow(xA-xC,2) + pow(yA-yC, 2) + pow(zA-zC, 2));
+
+ const PylithScalar k = (a + b + c) / 3.0;
+ const PylithScalar r = sqrt(k*(k-a)*(k-b)*(k-c)) / k;
+ minWidth = r;
+
+#endif
+
+ return minWidth;
+} // minCellWidth
+
+
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri3D.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri3D.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/GeometryTri3D.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -115,6 +115,13 @@
const int dim,
const int npts =1) const;
+ /** Compute minimum width across cell.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Minimum width across cell.
+ */
+ PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/Quadrature.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/Quadrature.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -135,7 +135,7 @@
/// Deallocate temporary storage.
void clear(void);
- /** Precompute geometric quantities for each cell.
+ /** Compute geometric quantities for cell.
*
* @param coordinatesCell Coordinates of vertices in cell.
* @param cell Finite-element cell
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/QuadratureRefCell.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/QuadratureRefCell.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/QuadratureRefCell.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -29,6 +29,7 @@
// Include directives ---------------------------------------------------
#include "feassemblefwd.hh" // forward declarations
+#include "CellGeometry.hh" // USES CellGeometry
#include "pylith/utils/array.hh" // HASA scalar_array
// Quadrature -----------------------------------------------------------
@@ -186,6 +187,13 @@
*/
int spaceDim(void) const;
+ /** Compute minimum width across cell.
+ *
+ * @param coordinatesCell Coordinates of vertices in cell.
+ * @returns Minimum width across cell.
+ */
+ PylithScalar minCellWidth(const scalar_array& coordinatesCell) const;
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/QuadratureRefCell.icc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/QuadratureRefCell.icc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/feassemble/QuadratureRefCell.icc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -90,6 +90,14 @@
return _spaceDim;
}
+// Compute minimum width across cell.
+inline
+PylithScalar
+pylith::feassemble::QuadratureRefCell::minCellWidth(const scalar_array& coordinatesCell) const {
+ assert(_geometry);
+ return _geometry->minCellWidth(coordinatesCell);
+}
+
#endif
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/materials/ElasticMaterial.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/materials/ElasticMaterial.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -274,6 +274,70 @@
} // stableTimeStepImplicit
// ----------------------------------------------------------------------
+// Get stable time step for explicit time integration.
+PylithScalar
+pylith::materials::ElasticMaterial::stableTimeStepExplicit(const topology::Mesh& mesh,
+ feassemble::Quadrature<topology::Mesh>* quadrature)
+{ // stableTimeStepImplicit
+ assert(quadrature);
+
+ const int numQuadPts = _numQuadPts;
+ const int numPropsQuadPt = _numPropsQuadPt;
+ const int tensorSize = _tensorSize;
+ const int numVarsQuadPt = _numVarsQuadPt;
+ assert(_propertiesCell.size() == numQuadPts*numPropsQuadPt);
+ assert(_stateVarsCell.size() == numQuadPts*numVarsQuadPt);
+ assert(_elasticConstsCell.size() == numQuadPts*_numElasticConsts);
+ assert(_initialStressCell.size() == numQuadPts*_tensorSize);
+ assert(_initialStrainCell.size() == numQuadPts*_tensorSize);
+
+ PylithScalar dtStable = pylith::PYLITH_MAXSCALAR;
+
+ // Get cells associated with material
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", id());
+ assert(!cells.isNull());
+ const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ const int spaceDim = quadrature->spaceDim();
+ const int numBasis = quadrature->numBasis();
+ scalar_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
+
+
+ for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ retrievePropsAndVars(*c_iter);
+
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ const double minCellWidth = quadrature->minCellWidth(coordinatesCell);
+
+#if 0
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const PylithScalar dt =
+ _stableTimeStepExplicit(&_propertiesCell[iQuad*numPropsQuadPt],
+ numPropsQuadPt,
+ &_stateVarsCell[iQuad*numVarsQuadPt],
+ numVarsQuadPt,
+ minCellWidth);
+ if (dt < dtStable)
+ dtStable = dt;
+ } // for
+#endif
+ } // for
+
+ return dtStable;
+} // stableTimeStepExplicit
+
+// ----------------------------------------------------------------------
// Allocate cell arrays.
void
pylith::materials::ElasticMaterial::_allocateCellArrays(void)
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/materials/ElasticMaterial.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/materials/ElasticMaterial.hh 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/materials/ElasticMaterial.hh 2012-08-22 19:06:20 UTC (rev 20620)
@@ -183,6 +183,21 @@
virtual
PylithScalar stableTimeStepImplicit(const topology::Mesh& mesh);
+ /** Get stable time step for explicit time integration.
+ *
+ * @pre Must call retrievePropsAndVars for cell before calling
+ * stableTimeStep().
+ *
+ * Default is MAXFLOAT (or 1.0e+30 if MAXFLOAT is not defined in math.h).
+ *
+ * @param mesh Finite-element mesh.
+ * @param quadrature Quadrature for finite-element integration
+ * @returns Time step
+ */
+ virtual
+ PylithScalar stableTimeStepExplicit(const topology::Mesh& mesh,
+ feassemble::Quadrature<topology::Mesh>* quadrature);
+
/** Get initial stress/strain fields.
*
* @returns Initial stress field.
@@ -308,10 +323,29 @@
*/
virtual
PylithScalar _stableTimeStepImplicit(const PylithScalar* properties,
- const int numProperties,
- const PylithScalar* stateVars,
- const int numStateVars) const = 0;
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars) const = 0;
+ /** Get stable time step for explicit time integration.
+ *
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param minCellWidth Minimum width across cell.
+ *
+ * @returns Time step
+ */
+#if 0
+ virtual
+ PylithScalar _stableTimeStepExplicit(const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars,
+ const double minCellWidth) const = 0;
+#endif
+
/** Compute 2D deviatoric stress/strain from vector and mean value.
*
* @param deviatoric Array for deviatoric tensor.
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/TestIntegrator.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/TestIntegrator.cc 2012-08-22 01:09:35 UTC (rev 20619)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/feassemble/TestIntegrator.cc 2012-08-22 19:06:20 UTC (rev 20620)
@@ -22,6 +22,7 @@
#include "pylith/feassemble/ElasticityExplicit.hh" // USES ElasticityExplicit
#include "pylith/feassemble/ElasticityImplicit.hh" // USES ElasticityImplicit
+#include "pylith/bc/Neumann.hh" // USES Neumann
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
#include "pylith/utils/constdefs.h" // USES MAXSCALAR
@@ -48,7 +49,7 @@
void
pylith::feassemble::TestIntegrator::testStableTimeStep(void)
{ // testStableTimeStep
- ElasticityExplicit integrator;
+ bc::Neumann integrator;
topology::Mesh mesh;
CPPUNIT_ASSERT_EQUAL(pylith::PYLITH_MAXSCALAR,
More information about the CIG-COMMITS
mailing list