[cig-commits] r12058 - in short/3D/PyLith/trunk/libsrc: bc faults feassemble
willic3 at geodynamics.org
willic3 at geodynamics.org
Thu May 29 13:13:34 PDT 2008
Author: willic3
Date: 2008-05-29 13:13:34 -0700 (Thu, 29 May 2008)
New Revision: 12058
Modified:
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh
short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
short/3D/PyLith/trunk/libsrc/bc/Neumann.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/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
Log:
Added GravityField and made associated changes.
Gravity should now be included in integrateResidual for implicit elasticity
problems. Also, coordinate system is now passed along with mesh object for
anything that uses integrateResidual.
Body forces have not yet been tested, but they do not appear to break any
of the unit tests. New unit tests and examples will need to be added.
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2008-05-29 19:50:01 UTC (rev 12057)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2008-05-29 20:13:34 UTC (rev 12058)
@@ -223,7 +223,8 @@
const ALE::Obj<real_section_type>& residual,
const double t,
topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh)
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs)
{ // integrateResidual
assert(0 != _quadrature);
assert(!_boundaryMesh.isNull());
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh 2008-05-29 19:50:01 UTC (rev 12057)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh 2008-05-29 20:13:34 UTC (rev 12058)
@@ -64,7 +64,15 @@
} // bc
} // pylith
+/*
+namespace spatialdata {
+ namespace geocoords {
+ class CoordSys; // USES CoordSys
+ } // geocoords
+} // spatialdata
+*/
+
/// C++ implementation of AbsorbingDampers boundary conditions.
class pylith::bc::AbsorbingDampers : public BoundaryCondition,
public feassemble::Integrator
@@ -97,11 +105,13 @@
* @param t Current time
* @param fields Solution fields
* @param mesh Finite-element mesh
+ * @param cs Mesh coordinate system
*/
void integrateResidual(const ALE::Obj<real_section_type>& residual,
const double t,
topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh);
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs);
/** Integrate contributions to Jacobian matrix (A) associated with
* operator.
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2008-05-29 19:50:01 UTC (rev 12057)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2008-05-29 20:13:34 UTC (rev 12058)
@@ -254,7 +254,8 @@
const ALE::Obj<real_section_type>& residual,
const double t,
topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh)
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs)
{ // integrateResidual
assert(0 != _quadrature);
assert(!_boundaryMesh.isNull());
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.hh 2008-05-29 19:50:01 UTC (rev 12057)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.hh 2008-05-29 20:13:34 UTC (rev 12058)
@@ -34,7 +34,15 @@
} // bc
} // pylith
+/*
+namespace spatialdata {
+ namespace geocoords {
+ class CoordSys; // USES CoordSys
+ } // geocoords
+} // spatialdata
+*/
+
/// C++ implementation of Neumann boundary conditions.
class pylith::bc::Neumann : public BoundaryCondition,
public feassemble::Integrator
@@ -67,11 +75,13 @@
* @param t Current time
* @param fields Solution fields
* @param mesh Finite-element mesh
+ * @param cs Mesh coordinate system
*/
void integrateResidual(const ALE::Obj<real_section_type>& residual,
const double t,
topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh);
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs);
/** Integrate contributions to Jacobian matrix (A) associated with
* operator.
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2008-05-29 19:50:01 UTC (rev 12057)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2008-05-29 20:13:34 UTC (rev 12058)
@@ -18,6 +18,8 @@
#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
#include "pylith/utils/array.hh" // USES double_array
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+
#include <assert.h> // USES assert()
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
@@ -53,7 +55,8 @@
const ALE::Obj<real_section_type>& residual,
const double t,
topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh)
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs)
{ // integrateResidual
throw std::logic_error("FaultCohesiveDyn::integrateResidual() not implemented.");
} // integrateResidual
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2008-05-29 19:50:01 UTC (rev 12057)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2008-05-29 20:13:34 UTC (rev 12058)
@@ -34,6 +34,15 @@
} // faults
} // pylith
+/*
+namespace spatialdata {
+ namespace geocoords {
+ class CoordSys; // USES CoordSys
+ } // geocoords
+} // spatialdata
+*/
+
+
/// @brief C++ implementation for a fault surface with spontaneous
/// (dynamic) slip implemented with cohesive elements.
class pylith::faults::FaultCohesiveDyn : public FaultCohesive,
@@ -77,11 +86,13 @@
* @param t Current time
* @param fields Solution fields
* @param mesh Finite-element mesh
+ * @param cs Mesh coordinate system
*/
void integrateResidual(const ALE::Obj<real_section_type>& residual,
const double t,
topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh);
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs);
/** Compute Jacobian matrix (A) associated with operator.
*
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-05-29 19:50:01 UTC (rev 12057)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-05-29 20:13:34 UTC (rev 12058)
@@ -105,7 +105,8 @@
const ALE::Obj<real_section_type>& residual,
const double t,
topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh)
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs)
{ // integrateResidual
assert(!residual.isNull());
assert(0 != fields);
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2008-05-29 19:50:01 UTC (rev 12057)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2008-05-29 20:13:34 UTC (rev 12058)
@@ -52,6 +52,15 @@
} // faults
} // pylith
+/*
+namespace spatialdata {
+ namespace geocoords {
+ class CoordSys; // USES CoordSys
+ } // geocoords
+} // spatialdata
+*/
+
+
/// C++ implementation for a fault surface with kinematic (prescribed)
/// slip implemented with cohesive elements.
class pylith::faults::FaultCohesiveKin : public FaultCohesive,
@@ -105,7 +114,8 @@
void integrateResidual(const ALE::Obj<real_section_type>& residual,
const double t,
topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh);
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs);
/** Integrate contributions to Jacobian matrix (A) associated with
* operator.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2008-05-29 19:50:01 UTC (rev 12057)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2008-05-29 20:13:34 UTC (rev 12058)
@@ -23,7 +23,7 @@
#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
#include "petscmat.h" // USES PetscMat
-#include "spatialdata/spatialdb/SpatialDB.hh"
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include <assert.h> // USES assert()
#include <stdexcept> // USES std::runtime_error
@@ -75,7 +75,8 @@
const ALE::Obj<real_section_type>& residual,
const double t,
topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh)
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs)
{ // integrateResidual
/// Member prototype for _elasticityResidualXD()
typedef void (pylith::feassemble::ElasticityExplicit::*elasticityResidual_fn_type)
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh 2008-05-29 19:50:01 UTC (rev 12057)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh 2008-05-29 20:13:34 UTC (rev 12058)
@@ -110,7 +110,8 @@
void integrateResidual(const ALE::Obj<real_section_type>& residual,
const double t,
topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh);
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs);
/** Integrate contributions to Jacobian matrix (A) associated with
* operator.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2008-05-29 19:50:01 UTC (rev 12057)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2008-05-29 20:13:34 UTC (rev 12058)
@@ -23,7 +23,8 @@
#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
#include "petscmat.h" // USES PetscMat
-#include "spatialdata/spatialdb/SpatialDB.hh"
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/spatialdb/GravityField.hh" // USES GravityField
#include <assert.h> // USES assert()
#include <stdexcept> // USES std::runtime_error
@@ -73,7 +74,8 @@
const ALE::Obj<real_section_type>& residual,
const double t,
topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh)
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs)
{ // integrateResidual
/// Member prototype for _elasticityResidualXD()
typedef void (pylith::feassemble::ElasticityImplicit::*elasticityResidual_fn_type)
@@ -160,7 +162,6 @@
_initCellVector();
const int cellVecSize = numBasis*spaceDim;
double_array dispTBctpdtCell(cellVecSize);
- //double_array gravCell(cellVecSize);
// Allocate vector for total strain
double_array totalStrain(numQuadPts*tensorSize);
@@ -194,33 +195,54 @@
const double_array& basis = _quadrature->basis();
const double_array& basisDeriv = _quadrature->basisDeriv();
const double_array& jacobianDet = _quadrature->jacobianDet();
+ const double_array& quadPts = _quadrature->quadPts();
if (cellDim != spaceDim)
throw std::logic_error("Not implemented yet.");
-#if 0
- // Comment out gravity section for now, until we figure out how to deal
- // with gravity vector.
- // Get density at quadrature points for this cell
- const double_array& density = _material->calcDensity();
+ // Compute body force vector if gravity is being used.
+ if (0 != _gravityField) {
- // Compute action for element body forces
- if (!grav.isNull()) {
- mesh->restrictClosure(grav, *c_iter, &gravCell[0], cellVecSize);
+ // Make sure coordinate names exist in gravity field.
+ _gravityField->open();
+ if (1 == spaceDim){
+ const char* queryNames[] = { "x"};
+ _gravityField->queryVals(queryNames, spaceDim);
+ } else if (2 == spaceDim){
+ const char* queryNames[] = { "x", "y"};
+ _gravityField->queryVals(queryNames, spaceDim);
+ } else if (3 == spaceDim){
+ const char* queryNames[] = { "x", "y", "z"};
+ _gravityField->queryVals(queryNames, spaceDim);
+ } else
+ assert(0);
+
+ // Get density at quadrature points for this cell
+ const double_array& density = _material->calcDensity();
+
+ // Compute action for element body forces
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ double gravVec[spaceDim];
+ double coords[spaceDim];
+ memcpy(coords, &quadPts[iQuad * spaceDim], sizeof(double)*spaceDim);
+
+ const int err = _gravityField->query(gravVec, spaceDim,
+ coords, spaceDim, cs);
+ if (err)
+ throw std::runtime_error("Unable to get gravity vector for point.");
const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad];
for (int iBasis=0, iQ=iQuad*numBasis*cellDim;
iBasis < numBasis; ++iBasis) {
const double valI = wt*basis[iQ+iBasis];
for (int iDim=0; iDim < spaceDim; ++iDim) {
- _cellVector[iBasis*spaceDim+iDim] += valI*gravCell[iDim];
+ _cellVector[iBasis*spaceDim+iDim] += valI*gravVec[iDim];
} // for
} // for
} // for
- PetscLogFlops(numQuadPts*(2+numBasis*(2+2*spaceDim)));
+ PetscLogFlops(numQuadPts*(2+numBasis*(1+2*spaceDim)));
+ _gravityField->close();
} // if
-#endif
// Compute B(transpose) * sigma, first computing strains
PetscLogEventBegin(stressEvent,0,0,0,0);
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh 2008-05-29 19:50:01 UTC (rev 12057)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh 2008-05-29 20:13:34 UTC (rev 12058)
@@ -99,18 +99,20 @@
* included in the residual (tractions, concentrated nodal forces,
* and contributions to internal force vector due to
* displacement/velocity BC). This routine computes the additional
- * external loads due to body forces (not yet implemented) plus the
+ * external loads due to body forces plus the
* element internal forces for the current stress state.
*
* @param residual Residual field (output)
* @param t Current time
* @param fields Solution fields
* @param mesh Mesh object
+ * @param cs Coordinate system
*/
void integrateResidual(const ALE::Obj<real_section_type>& residual,
const double t,
topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh);
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs);
/** Compute Jacobian matrix (A) associated with operator.
*
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc 2008-05-29 19:50:01 UTC (rev 12057)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc 2008-05-29 20:13:34 UTC (rev 12058)
@@ -16,6 +16,8 @@
#include "Quadrature.hh" // USES Quadrature
+#include "spatialdata/spatialdb/GravityField.hh" // USES GravityField
+
#include <assert.h> // USES assert()
// ----------------------------------------------------------------------
@@ -23,6 +25,7 @@
pylith::feassemble::Integrator::Integrator(void) :
_dt(-1.0),
_quadrature(0),
+ _gravityField(0),
_cellVector(0),
_cellMatrix(0),
_needNewJacobian(true),
@@ -35,6 +38,7 @@
pylith::feassemble::Integrator::~Integrator(void)
{ // destructor
delete _quadrature; _quadrature = 0;
+ delete _gravityField; _gravityField = 0;
delete[] _cellVector; _cellVector = 0;
delete[] _cellMatrix; _cellMatrix = 0;
} // destructor
@@ -53,6 +57,15 @@
} // quadrature
// ----------------------------------------------------------------------
+// Set gravity field.
+void
+pylith::feassemble::Integrator::gravityField(spatialdata::spatialdb::GravityField* const gravityField)
+{ // gravityField
+ _gravityField = gravityField;
+
+} // gravityField
+
+// ----------------------------------------------------------------------
// Initialize vector containing result of integration action for cell.
void
pylith::feassemble::Integrator::_initCellVector(void)
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2008-05-29 19:50:01 UTC (rev 12057)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2008-05-29 20:13:34 UTC (rev 12058)
@@ -40,6 +40,15 @@
} // topology
} // pylith
+namespace spatialdata {
+ namespace spatialdb {
+ class GravityField; // HOLDSA GravityField
+ } // spatialdb
+ namespace geocoords {
+ class CoordSys; // USES CoordSys
+ } // geocoords
+} // spatialdata
+
class pylith::feassemble::Integrator
{ // Integrator
friend class TestIntegrator; // unit testing
@@ -61,6 +70,12 @@
*/
void quadrature(const Quadrature* q);
+ /** Set gravity field. Gravity Field should already be initialized.
+ *
+ * @param g Gravity field.
+ */
+ void gravityField(spatialdata::spatialdb::GravityField* const gravityField);
+
/** Set time step for advancing from time t to time t+dt.
*
* @param dt Time step
@@ -98,12 +113,14 @@
* @param t Current time
* @param fields Solution fields
* @param mesh Finite-element mesh
+ * @param cs Mesh coordinate system
*/
virtual
void integrateResidual(const ALE::Obj<real_section_type>& residual,
const double t,
topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh) = 0;
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs) = 0;
/** Integrate contributions to Jacobian matrix (A) associated with
* operator.
@@ -168,6 +185,8 @@
Quadrature* _quadrature; ///< Quadrature for integrating finite-element
+ spatialdata::spatialdb::GravityField* _gravityField; ///< Gravity field.
+
/// Vector local to cell containing result of integration action
real_section_type::value_type* _cellVector;
More information about the cig-commits
mailing list