[cig-commits] r15152 - in short/3D/PyLith/trunk/libsrc: . bc
brad at geodynamics.org
brad at geodynamics.org
Mon Jun 8 22:29:58 PDT 2009
Author: brad
Date: 2009-06-08 22:29:57 -0700 (Mon, 08 Jun 2009)
New Revision: 15152
Added:
short/3D/PyLith/trunk/libsrc/bc/BCIntegratorSubMesh.cc
short/3D/PyLith/trunk/libsrc/bc/BCIntegratorSubMesh.hh
short/3D/PyLith/trunk/libsrc/bc/BCIntegratorSubMesh.icc
short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.cc
short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.hh
short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.icc
Modified:
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/bc/BoundaryConditionPoints.hh
short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.icc
short/3D/PyLith/trunk/libsrc/bc/Makefile.am
short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
short/3D/PyLith/trunk/libsrc/bc/TimeDependentPoints.hh
short/3D/PyLith/trunk/libsrc/bc/bcfwd.hh
Log:
Preparing for Neumann time-dependent BC.
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2009-06-09 00:19:01 UTC (rev 15151)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2009-06-09 05:29:57 UTC (rev 15152)
@@ -25,6 +25,7 @@
libpylith_la_SOURCES = \
bc/BoundaryCondition.cc \
bc/BoundaryConditionPoints.cc \
+ bc/BCIntegratorSubMesh.cc \
bc/TimeDependent.cc \
bc/TimeDependentPoints.cc \
bc/DirichletBC.cc \
Added: short/3D/PyLith/trunk/libsrc/bc/BCIntegratorSubMesh.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/BCIntegratorSubMesh.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/BCIntegratorSubMesh.cc 2009-06-09 05:29:57 UTC (rev 15152)
@@ -0,0 +1,102 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "BCIntegratorSubMesh.hh" // implementation of object methods
+
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::bc::BCIntegratorSubMesh::BCIntegratorSubMesh(void) :
+ _boundaryMesh(0),
+ _parameters(0)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::bc::BCIntegratorSubMesh::~BCIntegratorSubMesh(void)
+{ // destructor
+ delete _boundaryMesh; _boundaryMesh = 0;
+ delete _parameters; _parameters = 0;
+} // destructor
+
+// ----------------------------------------------------------------------
+// Get submesh associated with boundary condition.
+void
+pylith::bc::BCIntegratorSubMesh::createSubMesh(const topology::Mesh& mesh)
+{ // _createSubMesh
+ delete _boundaryMesh; _boundaryMesh = 0;
+ delete _parameters; _parameters = 0;
+
+ _boundaryMesh = new topology::SubMesh(mesh, _label.c_str());
+ assert(0 != _boundaryMesh);
+} // _createSubMesh
+
+// ----------------------------------------------------------------------
+// Verify configuration is acceptable.
+void
+pylith::bc::BCIntegratorSubMesh::verifyConfiguration(const topology::Mesh& mesh) const
+{ // verifyConfiguration
+ assert(0 != _quadrature);
+ assert(0 != _boundaryMesh);
+
+ BoundaryCondition::verifyConfiguration(mesh);
+
+ // check compatibility of quadrature and boundary mesh
+ if (_quadrature->cellDim() != _boundaryMesh->dimension()) {
+ std::ostringstream msg;
+ msg << "Quadrature is incompatible with cells for boundary condition '"
+ << _label << "'.\n"
+ << "Dimension of boundary mesh: " << _boundaryMesh->dimension()
+ << ", dimension of quadrature: " << _quadrature->cellDim()
+ << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ const int numCorners = _quadrature->numBasis();
+
+ // Get 'surface' cells (1 dimension lower than top-level cells)
+ const ALE::Obj<SieveSubMesh>& sieveSubMesh = _boundaryMesh->sieveMesh();
+ assert(!sieveSubMesh.isNull());
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells =
+ sieveSubMesh->heightStratum(1);
+ assert(!cells.isNull());
+ const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Make sure surface cells are compatible with quadrature.
+ const int boundaryDepth = sieveSubMesh->depth()-1; // depth of bndry cells
+ for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ const int cellNumCorners =
+ sieveSubMesh->getNumCellCorners(*c_iter, boundaryDepth);
+ if (numCorners != cellNumCorners) {
+ std::ostringstream msg;
+ msg << "Quadrature is incompatible with cell for boundary condition '"
+ << _label << "'.\n"
+ << "Cell " << *c_iter << " has " << cellNumCorners
+ << " vertices but quadrature reference cell has "
+ << numCorners << " vertices.";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // for
+
+} // verifyConfiguration
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/bc/BCIntegratorSubMesh.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/BCIntegratorSubMesh.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/BCIntegratorSubMesh.hh 2009-06-09 05:29:57 UTC (rev 15152)
@@ -0,0 +1,97 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/bc/BCIntegratorSubMesh.hh
+ *
+ * @brief C++ abstract base class for BoundaryCondition object with
+ * boundary condition applied at a simply-connected boundary (submesh).
+ *
+ * Interface definition for boundary conditions applied to a
+ * simply-connected boundary (submesh).
+ */
+
+#if !defined(pylith_bc_bcintegratorsubmesh_hh)
+#define pylith_bc_bcintegratorsubmesh_hh
+
+// Include directives ---------------------------------------------------
+#include "BoundaryCondition.hh" // ISA BoundaryCondition
+
+#include "pylith/topology/topologyfwd.hh" // HOLDSA Fields, SubMesh
+
+#include "pylith/topology/SubMesh.hh" // ISA Quadrature<SubMesh>
+#include "pylith/feassemble/Quadrature.hh" // ISA Integrator<Quadrature>
+#include "pylith/feassemble/Integrator.hh" // ISA Integrator
+
+// BCIntegratorSubMesh ----------------------------------------------
+class pylith::bc::BCIntegratorSubMesh : public BoundaryCondition,
+ public feassemble::Integrator<feassemble::Quadrature<topology::SubMesh> >
+{ // class BoundaryCondition
+ friend class TestBCIntegratorSubMesh; // unit testing
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Default constructor.
+ BCIntegratorSubMesh(void);
+
+ /// Destructor.
+ virtual
+ ~BCIntegratorSubMesh(void);
+
+ /** Get parameter fields.
+ *
+ * @returns Parameter fields.
+ */
+ const topology::Fields<topology::Field<topology::Mesh> >*
+ parameterFields(void) const;
+
+ /** Get boundary mesh.
+ *
+ * @return Boundary mesh.
+ */
+ const topology::SubMesh& boundaryMesh(void) const;
+
+ /** Get mesh labels for submesh associated with applied forces.
+ *
+ * @param mesh Finite-element mesh.
+ */
+ void createSubMesh(const topology::Mesh& mesh);
+
+ /** Verify configuration is acceptable.
+ *
+ * @param mesh Finite-element mesh
+ */
+ void verifyConfiguration(const topology::Mesh& mesh) const;
+
+ // PROTECTED MEMBERS //////////////////////////////////////////////////
+protected :
+
+ topology::SubMesh* _boundaryMesh; ///< Boundary mesh.
+
+ /// Parameters for boundary condition.
+ topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
+
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ /// Not implemented
+ BCIntegratorSubMesh(const BCIntegratorSubMesh&);
+
+ /// Not implemented
+ const BCIntegratorSubMesh& operator=(const BCIntegratorSubMesh&);
+
+}; // class BCIntegratorSubMesh
+
+#endif // pylith_bc_bcintegratorsubmesh_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/bc/BCIntegratorSubMesh.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/BCIntegratorSubMesh.icc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/BCIntegratorSubMesh.icc 2009-06-09 05:29:57 UTC (rev 15152)
@@ -0,0 +1,36 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#if !defined(pylith_bc_bcintegratorsubmesh_hh)
+#error "BCIntegratorSubMesh.icc can only be included from BCIntegratorSubMesh.hh"
+#endif
+
+#include <cassert> // USES assert()
+
+// ----------------------------------------------------------------------
+// Get parameter fields.
+const pylith::topology::Fields<pylith::topology::Field<pylith::topology::Mesh> >*
+pylith::bc::BCIntegratorSubMesh::parameterFields(void) const
+{ // parameterFields
+ return _parameters;
+} // parameterFields
+
+// ----------------------------------------------------------------------
+// Get boundary mesh.
+const pylith::topology::SubMesh&
+pylith::bc::BCIntegratorSubMesh::boundaryMesh(void) const {
+ assert(0 != _boundaryMesh);
+ return *_boundaryMesh;
+}
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/bc/BoundaryConditionPoints.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/BoundaryConditionPoints.hh 2009-06-09 00:19:01 UTC (rev 15151)
+++ short/3D/PyLith/trunk/libsrc/bc/BoundaryConditionPoints.hh 2009-06-09 05:29:57 UTC (rev 15152)
@@ -51,7 +51,7 @@
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
- /** Get mesh labels for points associated with applied forces.
+ /** Get mesh labels for points associated with boundary condition.
*
* @param mesh Finite-element mesh.
*/
@@ -63,7 +63,7 @@
/// Parameters for boundary condition.
topology::Fields<topology::Field<topology::Mesh> >* _parameters;
- int_array _points; ///< Points for forces.
+ int_array _points; ///< Points for boundary condition.
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.icc 2009-06-09 00:19:01 UTC (rev 15151)
+++ short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.icc 2009-06-09 05:29:57 UTC (rev 15152)
@@ -14,10 +14,13 @@
#error "DirichletBoundary.icc can only be included from DirichletBoundary.hh"
#endif
+#include <cassert> // USES assert()
+
// Get boundary mesh.
inline
const pylith::topology::SubMesh&
pylith::bc::DirichletBoundary::boundaryMesh(void) const {
+ assert(0 != _boundaryMesh);
return *_boundaryMesh;
} // boundaryMesh
Modified: short/3D/PyLith/trunk/libsrc/bc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Makefile.am 2009-06-09 00:19:01 UTC (rev 15151)
+++ short/3D/PyLith/trunk/libsrc/bc/Makefile.am 2009-06-09 05:29:57 UTC (rev 15152)
@@ -17,6 +17,7 @@
BoundaryCondition.hh \
BoundaryCondition.icc \
BoundaryConditionPoints.hh \
+ BCIntegratorSubMesh.hh \
TimeDependent.hh \
TimeDependent.icc \
TimeDependentPoints.hh \
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2009-06-09 00:19:01 UTC (rev 15151)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2009-06-09 05:29:57 UTC (rev 15152)
@@ -390,7 +390,7 @@
assert(0 != name);
if (0 == strcasecmp(name, "tractions")) {
- return _parameters->get("traction");;
+ return _parameters->get("traction");
} else {
std::ostringstream msg;
msg << "Unknown field '" << name << "' requested for Neumann BC '"
Added: short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.cc 2009-06-09 05:29:57 UTC (rev 15152)
@@ -0,0 +1,618 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "Neumann.hh" // implementation of object methods
+
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/Fields.hh" // USES Fields
+#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
+#include "spatialdata/spatialdb/TimeHistory.hh" // USES TimeHistory
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
+#include <cstring> // USES strcpy()
+#include <cassert> // USES assert()
+#include <stdexcept> // USES std::runtime_error
+#include <sstream> // USES std::ostringstream
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::SubMesh::SieveSubMesh SieveSubMesh;
+typedef pylith::topology::SubMesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::bc::Neumann::Neumann(void)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::bc::Neumann::~Neumann(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Initialize boundary condition. Determine orienation and compute traction
+// vector at integration points.
+void
+pylith::bc::Neumann::initialize(const topology::Mesh& mesh,
+ const double upDir[3])
+{ // initialize
+ double_array up(upDir, 3);
+ _queryDatabases(up);
+ _paramsLocalToGlobal();
+
+ // Get 'surface' cells (1 dimension lower than top-level cells)
+ const ALE::Obj<SieveSubMesh>& subSieveMesh = _boundaryMesh->sieveMesh();
+ assert(!subSieveMesh.isNull());
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells =
+ subSieveMesh->heightStratum(1);
+ assert(!cells.isNull());
+ const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Create section for traction vector in global coordinates
+ const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
+ const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
+ const int numBasis = _quadrature->numBasis();
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int spaceDim = cellGeometry.spaceDim();
+ const int fiberDim = spaceDim * numQuadPts;
+
+ _parameters =
+ new topology::Fields<topology::Field<topology::SubMesh> >(*_boundaryMesh);
+ assert(0 != _parameters);
+ _parameters->add("traction", "traction");
+ topology::Field<topology::SubMesh>& traction = _parameters->get("traction");
+ traction.newSection(cells, fiberDim);
+ traction.allocate();
+
+ // Containers for orientation information
+ const int orientationSize = spaceDim * spaceDim;
+ const int jacobianSize = spaceDim * cellDim;
+ double_array jacobian(jacobianSize);
+ double jacobianDet = 0;
+ double_array orientation(orientationSize);
+
+ // Set names based on dimension of problem.
+ // 1-D problem = {'normal-traction'}
+ // 2-D problem = {'shear-traction', 'normal-traction'}
+ // 3-D problem = {'horiz-shear-traction', 'vert-shear-traction',
+ // 'normal-traction'}
+ _db->open();
+ switch (spaceDim)
+ { // switch
+ case 1 : {
+ const char* valueNames[] = {"normal-traction"};
+ _db->queryVals(valueNames, 1);
+ break;
+ } // case 1
+ case 2 : {
+ const char* valueNames[] = {"shear-traction", "normal-traction"};
+ _db->queryVals(valueNames, 2);
+ break;
+ } // case 2
+ case 3 : {
+ const char* valueNames[] = {"horiz-shear-traction",
+ "vert-shear-traction",
+ "normal-traction"};
+ _db->queryVals(valueNames, 3);
+ break;
+ } // case 3
+ default :
+ assert(0);
+ } // switch
+
+ // Containers for database query results and quadrature coordinates in
+ // reference geometry.
+ double_array tractionDataLocal(spaceDim);
+ double_array quadPtRef(cellDim);
+ double_array quadPtsGlobal(numQuadPts*spaceDim);
+ const double_array& quadPtsRef = _quadrature->quadPtsRef();
+
+ // Container for cell tractions rotated to global coordinates.
+ double_array cellTractionsGlobal(fiberDim);
+
+ // Get sections.
+ double_array coordinatesCell(numCorners*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ subSieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
+ const ALE::Obj<SubRealSection>& tractionSection =
+ _parameters->get("traction").section();
+ assert(!tractionSection.isNull());
+
+ const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
+
+ assert(0 != _normalizer);
+ const double lengthScale = _normalizer->lengthScale();
+ const double pressureScale = _normalizer->pressureScale();
+
+ // Compute quadrature information
+ _quadrature->initializeGeometry();
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->computeGeometry(*_boundaryMesh, cells);
+#endif
+
+ // Loop over cells in boundary mesh, compute orientations, and then
+ // compute corresponding traction vector in global coordinates
+ // (store values in _tractionGlobal).
+ for(SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ subSieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
+ const double_array& quadPtsNondim = _quadrature->quadPts();
+ quadPtsGlobal = quadPtsNondim;
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
+ lengthScale);
+
+ cellTractionsGlobal = 0.0;
+ for(int iQuad=0, iRef=0, iSpace=0; iQuad < numQuadPts;
+ ++iQuad, iRef+=cellDim, iSpace+=spaceDim) {
+ // Get traction vector in local coordinate system at quadrature point
+ const int err = _db->query(&tractionDataLocal[0], spaceDim,
+ &quadPtsGlobal[iSpace], spaceDim, cs);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Could not find traction values at (";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << quadPtsGlobal[i+iSpace];
+ msg << ") for traction boundary condition " << _label << "\n"
+ << "using spatial database " << _db->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ _normalizer->nondimensionalize(&tractionDataLocal[0], spaceDim,
+ pressureScale);
+
+ // Compute Jacobian and determinant at quadrature point, then get
+ // orientation.
+ memcpy(&quadPtRef[0], &quadPtsRef[iRef], cellDim*sizeof(double));
+#if defined(PRECOMPUTE_GEOMETRY)
+ coordsVisitor.clear();
+ subSieveMesh->restrictClosure(*c_iter, coordsVisitor);
+#endif
+ cellGeometry.jacobian(&jacobian, &jacobianDet,
+ coordinatesCell, quadPtRef);
+ cellGeometry.orientation(&orientation, jacobian, jacobianDet, up);
+ assert(jacobianDet > 0.0);
+ orientation /= jacobianDet;
+
+ // Rotate traction vector from local coordinate system to global
+ // coordinate system
+ for(int iDim = 0; iDim < spaceDim; ++iDim) {
+ for(int jDim = 0; jDim < spaceDim; ++jDim)
+ cellTractionsGlobal[iDim+iSpace] +=
+ orientation[jDim*spaceDim+iDim] * tractionDataLocal[jDim];
+ } // for
+ } // for
+
+ // Update tractionsGlobal
+ tractionSection->updatePoint(*c_iter, &cellTractionsGlobal[0]);
+ } // for
+
+ _db->close();
+} // initialize
+
+// ----------------------------------------------------------------------
+// Integrate contributions to residual term (r) for operator.
+void
+pylith::bc::Neumann::integrateResidual(
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields)
+{ // integrateResidual
+ assert(0 != _quadrature);
+ assert(0 != _boundaryMesh);
+ assert(0 != _parameters);
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+
+ // Allocate vectors for cell values.
+ _initCellVector();
+ double_array tractionsCell(numQuadPts*spaceDim);
+
+ // Get cell information
+ const ALE::Obj<SieveSubMesh>& subSieveMesh = _boundaryMesh->sieveMesh();
+ assert(!subSieveMesh.isNull());
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells =
+ subSieveMesh->heightStratum(1);
+ assert(!cells.isNull());
+ const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<SubRealSection>& tractionSection =
+ _parameters->get("traction").section();
+ assert(!tractionSection.isNull());
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ topology::SubMesh::UpdateAddVisitor residualVisitor(*residualSection,
+ &_cellVector[0]);
+
+#if !defined(PRECOMPUTE_GEOMETRY)
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ subSieveMesh->getRealSection("coordinates");
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
+#endif
+
+ // Loop over faces and integrate contribution from each face
+ for (SieveSubMesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ subSieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
+
+ // Reset element vector to zero
+ _resetCellVector();
+
+ // Restrict tractions to cell
+ tractionSection->restrictPoint(*c_iter,
+ &tractionsCell[0], tractionsCell.size());
+
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
+ // Compute action for traction bc terms
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const double valI = wt*basis[iQuad*numBasis+iBasis];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const double valIJ = valI * basis[iQuad*numBasis+jBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _cellVector[iBasis*spaceDim+iDim] +=
+ tractionsCell[iQuad*spaceDim+iDim] * valIJ;
+ } // for
+ } // for
+ } // for
+ // Assemble cell contribution into field
+ residualVisitor.clear();
+ subSieveMesh->updateAdd(*c_iter, residualVisitor);
+
+ PetscLogFlops(numQuadPts*(1+numBasis*(1+numBasis*(1+2*spaceDim))));
+ } // for
+} // integrateResidual
+
+// ----------------------------------------------------------------------
+// Integrate contributions to Jacobian matrix (A) associated with
+void
+pylith::bc::Neumann::integrateJacobian(topology::Jacobian* jacobian,
+ const double t,
+ topology::SolutionFields* const fields)
+{ // integrateJacobian
+ _needNewJacobian = false;
+} // integrateJacobian
+
+// ----------------------------------------------------------------------
+// Verify configuration is acceptable.
+void
+pylith::bc::Neumann::verifyConfiguration(const topology::Mesh& mesh) const
+{ // verifyConfiguration
+ BCIntegratorSubMesh::verifyConfiguration(mesh);
+} // verifyConfiguration
+
+// ----------------------------------------------------------------------
+// Get cell field for tractions.
+const pylith::topology::Field<pylith::topology::SubMesh>&
+pylith::bc::Neumann::cellField(const char* name,
+ topology::SolutionFields* const fields)
+{ // cellField
+ assert(0 != _parameters);
+ assert(0 != name);
+
+ if (0 == strcasecmp(name, "tractions")) {
+ return _parameters->get("traction");
+ } else {
+ std::ostringstream msg;
+ msg << "Unknown field '" << name << "' requested for Neumann BC '"
+ << _label << "'.";
+ throw std::runtime_error(msg.str());
+ } // else
+
+ return _parameters->get("traction"); // Satisfy method definition
+} // cellField
+
+// ----------------------------------------------------------------------
+// Query databases for parameters.
+void
+pylith::bc::Neumann::_queryDatabases(void)
+{ // _queryDatabases
+ const double timeScale = _normalizer().timeScale();
+ const double rateScale = valueScale / timeScale;
+
+ const int numBCDOF = _bcDOF.size();
+ char** valueNames = (numBCDOF > 0) ? new char*[numBCDOF] : 0;
+ char** rateNames = (numBCDOF > 0) ? new char*[numBCDOF] : 0;
+ const std::string& valuePrefix = std::string(fieldName) + "-";
+ const std::string& ratePrefix = std::string(fieldName) + "-rate-";
+ const string_vector& components = _valueComponents();
+ for (int i=0; i < numBCDOF; ++i) {
+ std::string name = valuePrefix + components[_bcDOF[i]];
+ int size = 1 + name.length();
+ valueNames[i] = new char[size];
+ strcpy(valueNames[i], name.c_str());
+
+ name = ratePrefix + components[_bcDOF[i]];
+ size = 1 + name.length();
+ rateNames[i] = new char[size];
+ strcpy(rateNames[i], name.c_str());
+ } // for
+
+ delete _parameters;
+ _parameters =
+ new topology::Fields<topology::Field<topology::SubMesh> >(*_boundaryMesh);
+
+ // Create section to hold time dependent values
+ _parameters->add("value", fieldName);
+ topology::Field<topology::SubMesh>& value = _parameters->get("value");
+ value.scale(valueScale);
+ value.vectorFieldType(topology::FieldBase::OTHER);
+ value.newSection(topology::FieldBase::CELLS_FIELD, fiberDim);
+ value.allocate();
+
+ if (0 != _dbInitial) { // Setup initial values, if provided.
+ std::string fieldLabel = std::string("initial_") + std::string(fieldName);
+ _parameters->add("initial", fieldLabel.c_str());
+ topology::Field<topology::SubMesh>& initial =
+ _parameters->get("initial");
+ initial.cloneSection(value);
+ initial.scale(pressureScale);
+ initial.vectorFieldType(topology::FieldBase::OTHER);
+
+ _dbInitial->open();
+ _dbInitial->queryVals(valueNames, spaceDim);
+ _queryDB(&initial, _dbInitial, spaceDim, pressureScale);
+ _dbInitial->close();
+ } // if
+
+ if (0 != _dbRate) { // Setup rate of change of values, if provided.
+ std::string fieldLabel = std::string("rate_") + std::string(fieldName);
+ _parameters->add("rate", fieldLabel.c_str());
+ topology::Field<topology::SubMesh>& rate =
+ _parameters->get("rate");
+ rate.cloneSection(value);
+ rate.scale(rateScale);
+ rate.vectorFieldType(topology::FieldBase::OTHER);
+ const ALE::Obj<RealSection>& rateSection = rate.section();
+ assert(!rateSection.isNull());
+
+ _dbRate->open();
+ _dbRate->queryVals(rateNames, numBCDOF);
+ _queryDB(&rate, _dbRate, numBCDOF, rateScale);
+
+ std::string timeLabel =
+ std::string("rate_time_") + std::string(fieldName);
+ _parameters->add("rate time", timeLabel.c_str());
+ topology::Field<topology::SubMesh>& rateTime =
+ _parameters->get("rate time");
+ rateTime.newSection(rate, 1);
+ rateTime.allocate();
+ rateTime.scale(timeScale);
+ rateTime.vectorFieldType(topology::FieldBase::SCALAR);
+
+ const char* timeNames[1] = { "rate-start-time" };
+ _dbRate->queryVals(timeNames, 1);
+ _queryDB(&rateTime, _dbRate, 1, timeScale);
+ _dbRate->close();
+ } // if
+
+ if (0 != _dbChange) { // Setup change of values, if provided.
+ std::string fieldLabel = std::string("change_") + std::string(fieldName);
+ _parameters->add("change", fieldLabel.c_str());
+ topology::Field<topology::SubMesh>& change =
+ _parameters->get("change");
+ change.cloneSection(value);
+ change.scale(valueScale);
+ change.vectorFieldType(topology::FieldBase::OTHER);
+ const ALE::Obj<RealSection>& changeSection = change.section();
+ assert(!changeSection.isNull());
+
+ _dbChange->open();
+ _dbChange->queryVals(valueNames, numBCDOF);
+ _queryDB(&change, _dbChange, numBCDOF, valueScale);
+
+ std::string timeLabel =
+ std::string("change_time_") + std::string(fieldName);
+ _parameters->add("change time", timeLabel.c_str());
+ topology::Field<topology::SubMesh>& changeTime =
+ _parameters->get("change time");
+ changeTime.newSection(change, 1);
+ changeTime.allocate();
+ changeTime.scale(timeScale);
+ changeTime.vectorFieldType(topology::FieldBase::SCALAR);
+
+ const char* timeNames[1] = { "change-start-time" };
+ _dbChange->queryVals(timeNames, 1);
+ _queryDB(&changeTime, _dbChange, 1, timeScale);
+ _dbChange->close();
+
+ if (0 != _dbTimeHistory)
+ _dbTimeHistory->open();
+ } // if
+
+} // _queryDatabases
+
+// ----------------------------------------------------------------------
+// Query database for values.
+void
+pylith::bc::Neumann::_queryDB(topology::Field<topology::SubMesh>* field,
+ spatialdata::spatialdb::SpatialDB* const db,
+ const int querySize,
+ const double scale)
+{ // _queryDB
+ assert(0 != field);
+ assert(0 != db);
+ assert(0 != _boundaryMesh);
+
+ const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
+ assert(0 != cs);
+ const int spaceDim = cs->spaceDim();
+
+ const double lengthScale = _getNormalizer().lengthScale();
+
+ double_array coordsVertex(spaceDim);
+ const ALE::Obj<SieveMesh>& sieveMesh = _boundaryMesh->sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+
+ const ALE::Obj<RealSection>& section = field->section();
+ assert(!section.isNull());
+
+ double_array valuesCell(querySize);
+
+ // Get 'surface' cells (1 dimension lower than top-level cells)
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells =
+ sieveMesh->heightStratum(1);
+ assert(!cells.isNull());
+ const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Get dimensionalized coordinates of vertex
+ coordinates->restrictPoint(_Submesh[iPoint],
+ &coordsVertex[0], coordsVertex.size());
+ _getNormalizer().dimensionalize(&coordsVertex[0], coordsVertex.size(),
+ lengthScale);
+ int err = db->query(&valuesVertex[0], valuesVertex.size(),
+ &coordsVertex[0], coordsVertex.size(), cs);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Error querying for '" << field->label() << "' at (";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << coordsVertex[i];
+ msg << ") using spatial database " << db->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ _getNormalizer().nondimensionalize(&valuesVertex[0], valuesVertex.size(),
+ scale);
+ section->updatePoint(_Submesh[iPoint], &valuesVertex[0]);
+ } // for
+} // _queryDB
+
+// ----------------------------------------------------------------------
+// Calculate temporal and spatial variation of value over the list of Submesh.
+void
+pylith::bc::Neumann::_calculateValue(const double t)
+{ // _calculateValue
+ assert(0 != _parameters);
+
+ const ALE::Obj<RealSection>& valueSection =
+ _parameters->get("value").section();
+ assert(!valueSection.isNull());
+ valueSection->zero();
+
+ const int numSubmesh = _Submesh.size();
+ const int numBCDOF = _bcDOF.size();
+ const double timeScale = _getNormalizer().timeScale();
+
+ const ALE::Obj<RealSection>& initialSection = (0 != _dbInitial) ?
+ _parameters->get("initial").section() : 0;
+ const ALE::Obj<RealSection>& rateSection = ( 0 != _dbRate) ?
+ _parameters->get("rate").section() : 0;
+ const ALE::Obj<RealSection>& rateTimeSection = (0 != _dbRate) ?
+ _parameters->get("rate time").section() : 0;
+ const ALE::Obj<RealSection>& changeSection = ( 0 != _dbChange) ?
+ _parameters->get("change").section() : 0;
+ const ALE::Obj<RealSection>& changeTimeSection = ( 0 != _dbChange) ?
+ _parameters->get("change time").section() : 0;
+
+ double_array valuesVertex(numBCDOF);
+ double_array bufferVertex(numBCDOF);
+ for (int iPoint=0; iPoint < numSubmesh; ++iPoint) {
+ const int p_bc = _Submesh[iPoint]; // Get point label
+
+ valuesVertex = 0.0;
+
+ // Contribution from initial value
+ if (0 != _dbInitial) {
+ assert(!initialSection.isNull());
+ initialSection->restrictPoint(p_bc,
+ &bufferVertex[0], bufferVertex.size());
+ valuesVertex += bufferVertex;
+ } // if
+
+ // Contribution from rate of change of value
+ if (0 != _dbRate) {
+ assert(!rateSection.isNull());
+ assert(!rateTimeSection.isNull());
+ double tRate = 0.0;
+
+ rateSection->restrictPoint(p_bc, &bufferVertex[0], bufferVertex.size());
+ rateTimeSection->restrictPoint(p_bc, &tRate, 1);
+ if (t > tRate) { // rate of change integrated over time
+ bufferVertex *= (t - tRate);
+ valuesVertex += bufferVertex;
+ } // if
+ } // if
+
+ // Contribution from change of value
+ if (0 != _dbChange) {
+ assert(!changeSection.isNull());
+ assert(!changeTimeSection.isNull());
+ double tChange = 0.0;
+
+ changeSection->restrictPoint(p_bc, &bufferVertex[0], bufferVertex.size());
+ changeTimeSection->restrictPoint(p_bc, &tChange, 1);
+ if (t >= tChange) { // change in value over time
+ double scale = 1.0;
+ if (0 != _dbTimeHistory) {
+ double tDim = t - tChange;
+ _getNormalizer().dimensionalize(&tDim, 1, timeScale);
+ const int err = _dbTimeHistory->query(&scale, tDim);
+ if (0 != err) {
+ std::ostringstream msg;
+ msg << "Error querying for time '" << tDim
+ << "' in time history database "
+ << _dbTimeHistory->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // if
+ bufferVertex *= scale;
+ valuesVertex += bufferVertex;
+ } // if
+ } // if
+
+ valueSection->updateAddPoint(p_bc, &valuesVertex[0]);
+ } // for
+} // _calculateValue
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.hh 2009-06-09 05:29:57 UTC (rev 15152)
@@ -0,0 +1,142 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/bc/Neumann.hh
+ *
+ * @brief C++ implementation of time dependent Neumann (traction)
+ * boundary conditions applied to a simply-connected boundary.
+ */
+
+#if !defined(pylith_bc_neumann_hh)
+#define pylith_bc_neumann_hh
+
+// Include directives ---------------------------------------------------
+#include "BCIntegratorSubmesh.hh" // ISA BCIntegratorSubmesh
+#include "TimeDependent.hh" // ISA TimeDependent
+
+// Neumann ------------------------------------------------------
+class pylith::bc::Neumann : public BCIntegratorSubmesh,
+ public TimeDependent
+{ // class Neumann
+ friend class TestNeumann; // unit testing
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Default constructor.
+ Neumann(void);
+
+ /// Destructor.
+ ~Neumann(void);
+
+ /** Initialize boundary condition.
+ *
+ * @param mesh Finite-element mesh.
+ * @param upDir Direction perpendicular to horizontal surface tangent
+ * direction that is not collinear with surface normal.
+ */
+ void initialize(const topology::Mesh& mesh,
+ const double upDir[3]);
+
+ /** Integrate contributions to residual term (r) for operator.
+ *
+ * @param residual Field containing values for residual.
+ * @param t Current time.
+ * @param fields Solution fields.
+ */
+ void integrateResidual(const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields);
+
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
+ *
+ * @param jacobian Sparse matrix for Jacobian of system.
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ void integrateJacobian(topology::Jacobian* jacobian,
+ const double t,
+ topology::SolutionFields* const fields);
+
+ /** Verify configuration is acceptable.
+ *
+ * @param mesh Finite-element mesh
+ */
+ void verifyConfiguration(const topology::Mesh& mesh) const;
+
+ /** Get cell field with BC information.
+ *
+ * @param fieldType Type of field.
+ * @param name Name of field.
+ * @param mesh Finite-element mesh.
+ * @param fields Solution fields.
+ *
+ * @returns Traction vector at integration points.
+ */
+ const topology::Field<topology::SubMesh>&
+ cellField(const char* name,
+ topology::SolutionFields* const fields);
+
+ // PROTECTED METHODS //////////////////////////////////////////////////
+protected :
+
+ /** Get label of boundary condition surface.
+ *
+ * @returns Label of surface (from mesh generator).
+ */
+ const char* _getLabel(void) const;
+
+ /** Query databases for time dependent parameters.
+ *
+ * @param valueScale Dimension scale for value.
+ * @param fieldName Name of field associated with value.
+ */
+ void _queryDatabases(const double valueScale,
+ const char* fieldName);
+
+ /** Query database for values.
+ *
+ * @param field Field in which to store values.
+ * @param db Spatial database with values.
+ * @param querySize Number of values at each location.
+ * @param scale Dimension scale associated with values.
+ */
+ void _queryDB(topology::Field<topology::SubMesh>* field,
+ spatialdata::spatialdb::SpatialDB* const db,
+ const int querySize,
+ const double scale);
+
+ /// Convert parameters in local coordinates to global coordinates.
+ void _paramsLocalToGlobal(void);
+
+ /** Calculate spatial and temporal variation of value over the list
+ * of submesh.
+ *
+ * @param t Current time.
+ */
+ void _calculateValue(const double t);
+
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ Neumann(const Neumann&); ///< Not implemented.
+ const Neumann& operator=(const Neumann&); ///< Not implemented.
+
+}; // class Neumann
+
+#include "Neumann.icc"
+
+#endif // pylith_bc_neumann_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.icc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.icc 2009-06-09 05:29:57 UTC (rev 15152)
@@ -0,0 +1,25 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#if !defined(pylith_bc_neumann_hh)
+#error "Neumann.icc can only be included from Neumann.hh"
+#endif
+
+// Get label of boundary condition surface.
+inline
+const char*
+pylith::bc::Neumann::_getLabel(void) const {
+ return label();
+}
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/bc/TimeDependentPoints.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/TimeDependentPoints.hh 2009-06-09 00:19:01 UTC (rev 15151)
+++ short/3D/PyLith/trunk/libsrc/bc/TimeDependentPoints.hh 2009-06-09 05:29:57 UTC (rev 15152)
@@ -59,7 +59,7 @@
const double valueScale,
const char* fieldName);
- /** Wuery database for values.
+ /** Query database for values.
*
* @param field Field in which to store values.
* @param db Spatial database with values.
Modified: short/3D/PyLith/trunk/libsrc/bc/bcfwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/bcfwd.hh 2009-06-09 00:19:01 UTC (rev 15151)
+++ short/3D/PyLith/trunk/libsrc/bc/bcfwd.hh 2009-06-09 05:29:57 UTC (rev 15152)
@@ -26,6 +26,7 @@
class BoundaryCondition;
class BoundaryConditionPoints;
+ class BCIntegratorSubMesh;
class TimeDependent;
class TimeDependentPoints;
class TimeDependentSubMesh;
More information about the CIG-COMMITS
mailing list