[cig-commits] r15337 - in short/3D/PyLith/trunk: examples/twocells/twoquad4 libsrc libsrc/bc modulesrc/bc pylith/bc pylith/meshio unittests/libtests/bc unittests/libtests/bc/data unittests/pytests/bc unittests/pytests/bc/data
brad at geodynamics.org
brad at geodynamics.org
Thu Jun 18 16:43:22 PDT 2009
Author: brad
Date: 2009-06-18 16:43:10 -0700 (Thu, 18 Jun 2009)
New Revision: 15337
Added:
short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
short/3D/PyLith/trunk/libsrc/bc/Neumann.icc
short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.cc
short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.hh
short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.icc
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.hh
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_OLD.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_OLD.hh
short/3D/PyLith/trunk/unittests/pytests/bc/data/tri3_vel.spatialdb
Removed:
short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
short/3D/PyLith/trunk/libsrc/bc/Neumann.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
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.hh
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.hh
Modified:
short/3D/PyLith/trunk/examples/twocells/twoquad4/axialtract.cfg
short/3D/PyLith/trunk/examples/twocells/twoquad4/axialtract.spatialdb
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/bc/Makefile.am
short/3D/PyLith/trunk/modulesrc/bc/Neumann.i
short/3D/PyLith/trunk/pylith/bc/Neumann.py
short/3D/PyLith/trunk/pylith/meshio/OutputDirichlet.py
short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannHex8.hh
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannLine2.hh
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannQuad4.hh
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannTet4.hh
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannTri3.hh
short/3D/PyLith/trunk/unittests/libtests/bc/data/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/bc/data/hex8b_traction.spatialdb
short/3D/PyLith/trunk/unittests/libtests/bc/data/line2_tractions.spatialdb
short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_tractions.spatialdb
short/3D/PyLith/trunk/unittests/libtests/bc/data/tet4_tractions.spatialdb
short/3D/PyLith/trunk/unittests/libtests/bc/data/tri3_tractions.spatialdb
short/3D/PyLith/trunk/unittests/pytests/bc/TestDirichletBoundary.py
short/3D/PyLith/trunk/unittests/pytests/bc/TestNeumann.py
short/3D/PyLith/trunk/unittests/pytests/bc/data/Makefile.am
short/3D/PyLith/trunk/unittests/pytests/bc/data/tri3_tractions.spatialdb
short/3D/PyLith/trunk/unittests/pytests/bc/testbc.py
Log:
Replaced old Neumann BC with new time-dependent implementation.
Modified: short/3D/PyLith/trunk/examples/twocells/twoquad4/axialtract.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/twocells/twoquad4/axialtract.cfg 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/examples/twocells/twoquad4/axialtract.cfg 2009-06-18 23:43:10 UTC (rev 15337)
@@ -94,15 +94,15 @@
# condition have the name 'x_pos'.
label = x_pos
-db = spatialdata.spatialdb.SimpleDB
+db_initial = spatialdata.spatialdb.SimpleDB
# We are assigning the label 'Dirichlet BC +x edge' to the database.
-db.label = Neumann BC +x edge
+db_initial.label = Neumann BC +x edge
# The name of the file containing the spatial database for the BC
# specification.
-db.iohandler.filename = axialtract.spatialdb
-db.query_type = nearest
+db_initial.iohandler.filename = axialtract.spatialdb
+db_initial.query_type = nearest
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 1
@@ -168,7 +168,7 @@
# Give basename for VTK output of traction BC information.
[pylithapp.timedependent.bc.x_pos.output]
-cell_info_fields = [tractions]
+cell_info_fields = [initial-value]
writer.filename = axialtract-tractions.vtk
# Give basename for VTK fault output.
Modified: short/3D/PyLith/trunk/examples/twocells/twoquad4/axialtract.spatialdb
===================================================================
--- short/3D/PyLith/trunk/examples/twocells/twoquad4/axialtract.spatialdb 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/examples/twocells/twoquad4/axialtract.spatialdb 2009-06-18 23:43:10 UTC (rev 15337)
@@ -6,7 +6,7 @@
// There are two values specified in the database, corresponding to the
// two components of the traction vector.
num-values = 2
- value-names = shear-traction normal-traction
+ value-names = traction-shear traction-normal
// The tractions have units of Pa.
value-units = Pa Pa
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2009-06-18 23:43:10 UTC (rev 15337)
@@ -31,7 +31,6 @@
bc/DirichletBC.cc \
bc/DirichletBoundary.cc \
bc/Neumann.cc \
- bc/Neumann_NEW.cc \
bc/AbsorbingDampers.cc \
bc/PointForce.cc \
faults/Fault.cc \
Modified: short/3D/PyLith/trunk/libsrc/bc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Makefile.am 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/libsrc/bc/Makefile.am 2009-06-18 23:43:10 UTC (rev 15337)
@@ -31,8 +31,6 @@
DirichletBoundary.icc \
Neumann.hh \
Neumann.icc \
- Neumann_NEW.hh \
- Neumann_NEW.icc \
PointForce.hh \
PointForce.icc \
bcfwd.hh
Deleted: short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,376 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-#include <portinfo>
-
-#include "Neumann.hh" // implementation of object methods
-
-#include "pylith/topology/Field.hh" // HOLDSA Field
-#include "pylith/feassemble/CellGeometry.hh" // USES CellGeometry
-
-#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
-#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
-#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
-
-#include <cstring> // USES memcpy()
-#include <strings.h> // USES strcasecmp()
-#include <cassert> // USES assert()
-#include <stdexcept> // USES std::runtime_error
-#include <sstream> // USES std::ostringstream
-
-//#define PRECOMPUTE_GEOMETRY
-
-// ----------------------------------------------------------------------
-typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-typedef pylith::topology::SubMesh::RealSection SubRealSection;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
-
-// ----------------------------------------------------------------------
-// Default constructor.
-pylith::bc::Neumann::Neumann(void) :
- _db(0)
-{ // constructor
-} // constructor
-
-// ----------------------------------------------------------------------
-// Destructor.
-pylith::bc::Neumann::~Neumann(void)
-{ // destructor
- deallocate();
-} // destructor
-
-// ----------------------------------------------------------------------
-// Deallocate PETSc and local data structures.
-void
-pylith::bc::Neumann::deallocate(void)
-{ // deallocate
- _db = 0; // :TODO: Use shared pointer
-} // deallocate
-
-// ----------------------------------------------------------------------
-// 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
- assert(0 != _boundaryMesh);
- assert(0 != _quadrature);
- assert(0 != _db);
-
- double_array up(upDir, 3);
- const int numCorners = _quadrature->numBasis();
-
- assert(0 != _normalizer);
- const double lengthScale = _normalizer->lengthScale();
- const double pressureScale = _normalizer->pressureScale();
-
- // 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;
- ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("BoundaryConditions");
-
- delete _parameters;
- _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();
- traction.scale(pressureScale);
- traction.vectorFieldType(topology::FieldBase::VECTOR);
-
- logger.stagePop();
-
- // 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 :
- std::cerr << "Bad spatial dimension '" << spaceDim << "'." << std::endl;
- assert(0);
- throw std::logic_error("Bad spatial dimension in Neumann.");
- } // 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();
-
- // 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
-
-
-// End of file
Copied: short/3D/PyLith/trunk/libsrc/bc/Neumann.cc (from rev 15335, short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.cc)
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2009-06-18 23:43:10 UTC (rev 15337)
@@ -0,0 +1,759 @@
+// -*- 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 <strings.h> // USES strcasecmp()
+#include <cassert> // USES assert()
+#include <stdexcept> // USES std::runtime_error
+#include <sstream> // USES std::ostringstream
+
+//#define PRECOMPUTE_GEOMETRY
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
+typedef pylith::topology::SubMesh::RealSection SubRealSection;
+typedef pylith::topology::Mesh::RealSection RealSection;
+typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::bc::Neumann::Neumann(void)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::bc::Neumann::~Neumann(void)
+{ // destructor
+ deallocate();
+} // destructor
+
+// ----------------------------------------------------------------------
+// Deallocate PETSc and local data structures.
+void
+pylith::bc::Neumann::deallocate(void)
+{ // deallocate
+} // deallocate
+
+// ----------------------------------------------------------------------
+// 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
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("BoundaryConditions");
+
+ _queryDatabases();
+ _paramsLocalToGlobal(upDir);
+
+ logger.stagePop();
+} // 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
+ _calculateValue(t);
+ const ALE::Obj<SubRealSection>& tractionSection =
+ _parameters->get("value").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
+ if (1 == mesh.dimension())
+ throw std::runtime_error("Neumann boundary conditions are not "
+ "implemented for a 1-D mesh.");
+
+ 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, "initial-value"))
+ return _parameters->get("initial");
+
+ else if (0 == strcasecmp(name, "rate-of-change"))
+ return _parameters->get("rate");
+
+ else if (0 == strcasecmp(name, "change-in-value"))
+ return _parameters->get("change");
+
+ else if (0 == strcasecmp(name, "rate-start-time"))
+ return _parameters->get("rate time");
+
+ else if (0 == strcasecmp(name, "change-start-time"))
+ return _parameters->get("change time");
+
+ 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
+ assert(0 != _quadrature);
+ assert(0 != _boundaryMesh);
+
+ const double pressureScale = _normalizer->pressureScale();
+ const double timeScale = _normalizer->timeScale();
+ const double rateScale = pressureScale / timeScale;
+
+ const int spaceDim = _quadrature->spaceDim();
+ const int numQuadPts = _quadrature->numQuadPts();
+
+ delete _parameters;
+ _parameters =
+ new topology::Fields<topology::Field<topology::SubMesh> >(*_boundaryMesh);
+
+ // Create section to hold time dependent values
+ _parameters->add("value", "traction");
+ topology::Field<topology::SubMesh>& value = _parameters->get("value");
+ value.scale(pressureScale);
+ value.vectorFieldType(topology::FieldBase::OTHER);
+ value.newSection(topology::FieldBase::CELLS_FIELD, numQuadPts*spaceDim, 1);
+ value.allocate();
+
+ if (0 != _dbInitial) { // Setup initial values, if provided.
+ _parameters->add("initial", "initial_traction");
+ topology::Field<topology::SubMesh>& initial =
+ _parameters->get("initial");
+ initial.cloneSection(value);
+ initial.scale(pressureScale);
+ initial.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
+
+ _dbInitial->open();
+ switch (spaceDim)
+ { // switch
+ case 1 : {
+ const char* valueNames[] = {"traction-normal"};
+ _dbInitial->queryVals(valueNames, 1);
+ break;
+ } // case 1
+ case 2 : {
+ const char* valueNames[] = {"traction-shear", "traction-normal"};
+ _dbInitial->queryVals(valueNames, 2);
+ break;
+ } // case 2
+ case 3 : {
+ const char* valueNames[] = {"traction-shear-horiz",
+ "traction-shear-vert",
+ "traction-normal"};
+ _dbInitial->queryVals(valueNames, 3);
+ break;
+ } // case 3
+ default :
+ std::cerr << "Bad spatial dimension '" << spaceDim << "'." << std::endl;
+ assert(0);
+ throw std::logic_error("Bad spatial dimension in Neumann.");
+ } // switch
+ _queryDB(&initial, _dbInitial, spaceDim, pressureScale);
+ _dbInitial->close();
+ } // if
+
+ if (0 != _dbRate) { // Setup rate of change of values, if provided.
+ _parameters->add("rate", "traction_rate");
+ topology::Field<topology::SubMesh>& rate =
+ _parameters->get("rate");
+ rate.cloneSection(value);
+ rate.scale(rateScale);
+ rate.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
+ const ALE::Obj<RealSection>& rateSection = rate.section();
+ assert(!rateSection.isNull());
+
+ _dbRate->open();
+ switch (spaceDim)
+ { // switch
+ case 1 : {
+ const char* valueNames[] = {"traction-rate-normal"};
+ _dbRate->queryVals(valueNames, 1);
+ break;
+ } // case 1
+ case 2 : {
+ const char* valueNames[] = {"traction-rate-shear",
+ "traction-rate-normal"};
+ _dbRate->queryVals(valueNames, 2);
+ break;
+ } // case 2
+ case 3 : {
+ const char* valueNames[] = {"traction-rate-shear-horiz",
+ "traction-rate-shear-vert",
+ "traction-rate-normal"};
+ _dbRate->queryVals(valueNames, 3);
+ break;
+ } // case 3
+ default :
+ std::cerr << "Bad spatial dimension '" << spaceDim << "'." << std::endl;
+ assert(0);
+ throw std::logic_error("Bad spatial dimension in Neumann.");
+ } // switch
+ _queryDB(&rate, _dbRate, spaceDim, rateScale);
+
+ _parameters->add("rate time", "rate_traction_time");
+ topology::Field<topology::SubMesh>& rateTime =
+ _parameters->get("rate time");
+ rateTime.newSection(rate, numQuadPts);
+ rateTime.allocate();
+ rateTime.scale(timeScale);
+ rateTime.vectorFieldType(topology::FieldBase::MULTI_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.
+ _parameters->add("change", "change_traction");
+ topology::Field<topology::SubMesh>& change =
+ _parameters->get("change");
+ change.cloneSection(value);
+ change.scale(pressureScale);
+ change.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
+ const ALE::Obj<RealSection>& changeSection = change.section();
+ assert(!changeSection.isNull());
+
+ _dbChange->open();
+ switch (spaceDim)
+ { // switch
+ case 1 : {
+ const char* valueNames[] = {"traction-normal"};
+ _dbChange->queryVals(valueNames, 1);
+ break;
+ } // case 1
+ case 2 : {
+ const char* valueNames[] = {"traction-shear", "traction-normal"};
+ _dbChange->queryVals(valueNames, 2);
+ break;
+ } // case 2
+ case 3 : {
+ const char* valueNames[] = {"traction-shear-horiz",
+ "traction-shear-vert",
+ "traction-normal"};
+ _dbChange->queryVals(valueNames, 3);
+ break;
+ } // case 3
+ default :
+ std::cerr << "Bad spatial dimension '" << spaceDim << "'." << std::endl;
+ assert(0);
+ throw std::logic_error("Bad spatial dimension in Neumann.");
+ } // switch
+ _queryDB(&change, _dbChange, spaceDim, pressureScale);
+
+ _parameters->add("change time", "change_traction_time");
+ topology::Field<topology::SubMesh>& changeTime =
+ _parameters->get("change time");
+ changeTime.newSection(change, numQuadPts);
+ changeTime.allocate();
+ changeTime.scale(timeScale);
+ changeTime.vectorFieldType(topology::FieldBase::MULTI_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);
+ assert(0 != _quadrature);
+
+ // 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();
+
+ const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
+ const int numBasis = _quadrature->numBasis();
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int spaceDim = _quadrature->spaceDim();
+
+ // Containers for database query results and quadrature coordinates in
+ // reference geometry.
+ double_array valuesCell(numQuadPts*querySize);
+ double_array quadPtsGlobal(numQuadPts*spaceDim);
+
+ // Get sections.
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ subSieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
+ const ALE::Obj<RealSection>& section = field->section();
+ assert(!section.isNull());
+
+ const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
+
+ assert(0 != _normalizer);
+ const double lengthScale = _normalizer->lengthScale();
+
+ // Compute quadrature information
+ _quadrature->initializeGeometry();
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->computeGeometry(*_boundaryMesh, cells);
+#endif
+
+ // Loop over cells in boundary mesh and perform queries.
+ 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);
+
+ valuesCell = 0.0;
+ for (int iQuad=0, iSpace=0;
+ iQuad < numQuadPts;
+ ++iQuad, iSpace+=spaceDim) {
+ const int err = db->query(&valuesCell[iQuad*querySize], querySize,
+ &quadPtsGlobal[iSpace], spaceDim, cs);
+
+ if (err) {
+ std::ostringstream msg;
+ msg << "Could not find 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
+
+ } // for
+ _normalizer->nondimensionalize(&valuesCell[0], valuesCell.size(),
+ scale);
+
+ // Update section
+ assert(valuesCell.size() == section->getFiberDimension(*c_iter));
+ section->updatePoint(*c_iter, &valuesCell[0]);
+ } // for
+} // _queryDB
+
+// ----------------------------------------------------------------------
+// Initialize boundary condition. Determine orienation and compute traction
+// vector at integration points.
+void
+ pylith::bc::Neumann::_paramsLocalToGlobal(const double upDir[3])
+{ // _paramsLocalToGlobal
+ assert(0 != _boundaryMesh);
+ assert(0 != _quadrature);
+
+ double_array up(upDir, 3);
+
+
+ // 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();
+
+ // Quadrature related values.
+ 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 = numQuadPts * spaceDim;
+ double_array quadPtRef(cellDim);
+ const double_array& quadPtsRef = _quadrature->quadPtsRef();
+
+ // 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);
+
+ // Container for cell tractions rotated to global coordinates.
+ double_array initialCellLocal(fiberDim);
+ double_array initialCellGlobal(fiberDim);
+ double_array rateCellLocal(fiberDim);
+ double_array rateCellGlobal(fiberDim);
+ double_array changeCellLocal(fiberDim);
+ double_array changeCellGlobal(fiberDim);
+
+ // Get sections.
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ subSieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
+ 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>& changeSection = ( 0 != _dbChange) ?
+ _parameters->get("change").section() : 0;
+
+ // Loop over cells in boundary mesh, compute orientations, and then
+ // rotate corresponding traction vector from local to global coordinates.
+ 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
+ // Reset traction vectors
+ initialCellLocal = 0.0;
+ initialCellGlobal = 0.0;
+ rateCellLocal = 0.0;
+ rateCellGlobal = 0.0;
+ changeCellLocal = 0.0;
+ changeCellGlobal = 0.0;
+
+ // Get values for cell from each of the sections
+ if (0 != _dbInitial) {
+ assert(!initialSection.isNull());
+ initialSection->restrictPoint(*c_iter, &initialCellLocal[0],
+ initialCellLocal.size());
+ } // if
+ if (0 != _dbRate) {
+ assert(!rateSection.isNull());
+ rateSection->restrictPoint(*c_iter, &rateCellLocal[0],
+ rateCellLocal.size());
+ } // if
+ if (0 != _dbChange) {
+ assert(!changeSection.isNull());
+ changeSection->restrictPoint(*c_iter, &changeCellLocal[0],
+ changeCellLocal.size());
+ } // if
+
+ for(int iQuad=0, iRef=0, iSpace=0; iQuad < numQuadPts;
+ ++iQuad, iRef+=cellDim, iSpace+=spaceDim) {
+ // 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;
+
+ if (0 != _dbInitial) {
+ assert(!initialSection.isNull());
+ // 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)
+ initialCellGlobal[iSpace+iDim] +=
+ orientation[jDim*spaceDim+iDim] * initialCellLocal[iSpace+jDim];
+ } // for
+ } // if
+
+ if (0 != _dbRate) {
+ assert(!rateSection.isNull());
+ // 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)
+ rateCellGlobal[iDim+iSpace] +=
+ orientation[jDim*spaceDim+iDim] * rateCellLocal[iSpace+jDim];
+ } // for
+ } // if
+
+ if (0 != _dbChange) {
+ assert(!changeSection.isNull());
+ // 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)
+ changeCellGlobal[iDim+iSpace] +=
+ orientation[jDim*spaceDim+iDim] * changeCellLocal[iSpace+jDim];
+ } // for
+ } // if
+
+ } // for
+
+ // Update sections
+ if (0 != _dbInitial)
+ initialSection->updatePoint(*c_iter, &initialCellGlobal[0]);
+ if (0 != _dbRate)
+ rateSection->updatePoint(*c_iter, &rateCellGlobal[0]);
+ if (0 != _dbChange)
+ changeSection->updatePoint(*c_iter, &changeCellGlobal[0]);
+ } // for
+} // paramsLocalToGlobal
+
+// ----------------------------------------------------------------------
+// Calculate temporal and spatial variation of value over the list of Submesh.
+void
+pylith::bc::Neumann::_calculateValue(const double t)
+{ // _calculateValue
+ assert(0 != _parameters);
+ assert(0 != _boundaryMesh);
+ assert(0 != _quadrature);
+
+ const ALE::Obj<RealSection>& valueSection =
+ _parameters->get("value").section();
+ assert(!valueSection.isNull());
+ valueSection->zero();
+
+ 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;
+
+ // 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();
+
+ const int spaceDim = _quadrature->spaceDim();
+ const int numBasis = _quadrature->numBasis();
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int fiberDim = numQuadPts * spaceDim;
+
+ double_array valuesCell(fiberDim);
+ double_array bufferCell(fiberDim);
+ double_array timeCell(numQuadPts);
+ for(SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+
+ valuesCell = 0.0;
+ timeCell = 0.0;
+
+ // Contribution from initial value
+ if (0 != _dbInitial) {
+ assert(!initialSection.isNull());
+ initialSection->restrictPoint(*c_iter,
+ &bufferCell[0], bufferCell.size());
+ valuesCell += bufferCell;
+ } // if
+
+ // Contribution from rate of change of value
+ if (0 != _dbRate) {
+ assert(!rateSection.isNull());
+ assert(!rateTimeSection.isNull());
+
+ rateSection->restrictPoint(*c_iter, &bufferCell[0], bufferCell.size());
+ rateTimeSection->restrictPoint(*c_iter, &timeCell[0], timeCell.size());
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
+ if (t > timeCell[iQuad]) // rate of change integrated over time
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ valuesCell[iQuad*spaceDim+iDim] +=
+ bufferCell[iQuad*spaceDim+iDim] * (t - timeCell[iQuad]);
+ } // if
+
+ // Contribution from change of value
+ if (0 != _dbChange) {
+ assert(!changeSection.isNull());
+ assert(!changeTimeSection.isNull());
+
+ changeSection->restrictPoint(*c_iter, &bufferCell[0], bufferCell.size());
+ changeTimeSection->restrictPoint(*c_iter, &timeCell[0], timeCell.size());
+ for (int i=0; i < numQuadPts; ++i)
+ if (t >= timeCell[i]) { // change in value over time
+ double scale = 1.0;
+ if (0 != _dbTimeHistory) {
+ double tDim = t - timeCell[i];
+ _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
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ valuesCell[i*spaceDim+iDim] += bufferCell[i*spaceDim+iDim] * scale;
+ } // if
+ } // if
+
+ valueSection->updateAddPoint(*c_iter, &valuesCell[0]);
+ } // for
+} // _calculateValue
+
+
+// End of file
Property changes on: short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
___________________________________________________________________
Name: svn:mergeinfo
+
Deleted: short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.hh 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.hh 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,115 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-/** @file libsrc/bc/Neumann.hh
- *
- * @brief C++ implementation of Neumann (prescribed tractions
- * on a surface) boundary conditions.
- */
-
-#if !defined(pylith_bc_neumann_hh)
-#define pylith_bc_neumann_hh
-
-// Include directives ---------------------------------------------------
-#include "BCIntegratorSubMesh.hh" // ISA BCIntegratorSubMesh
-
-// Neumann --------------------------------------------------------------
-class pylith::bc::Neumann : public BCIntegratorSubMesh
-{ // class Neumann
- friend class TestNeumann; // unit testing
-
- // PUBLIC METHODS /////////////////////////////////////////////////////
-public :
-
- /// Default constructor.
- Neumann(void);
-
- /// Destructor.
- ~Neumann(void);
-
- /// Deallocate PETSc and local data structures.
- void deallocate(void);
-
- /** Set database for boundary condition parameters.
- *
- * @param db Spatial database
- */
- void db(spatialdata::spatialdb::SpatialDB* const db);
-
- /** 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 =0);
-
- // PRIVATE MEMBERS ////////////////////////////////////////////////////
-private :
-
- spatialdata::spatialdb::SpatialDB* _db; ///< Spatial database w/parameters
-
- // NOT IMPLEMENTED ////////////////////////////////////////////////////
-private :
-
- Neumann(const Neumann&); ///< Not implemented
- const Neumann& operator=(const Neumann&); ///< Not implemented
-
-}; // class Neumann
-
-#include "Neumann.icc" // inline methods
-
-#endif // pylith_bc_neumann_hh
-
-
-// End of file
Copied: short/3D/PyLith/trunk/libsrc/bc/Neumann.hh (from rev 15333, short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.hh)
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.hh 2009-06-18 23:43:10 UTC (rev 15337)
@@ -0,0 +1,150 @@
+// -*- 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);
+
+ /// Deallocate PETSc and local data structures.
+ void deallocate(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 =0);
+
+ // PROTECTED METHODS //////////////////////////////////////////////////
+protected :
+
+ /** Get label of boundary condition surface.
+ *
+ * @returns Label of surface (from mesh generator).
+ */
+ const char* _getLabel(void) const;
+
+ /** Get manager of scales used to nondimensionalize problem.
+ *
+ * @returns Nondimensionalizer.
+ */
+ const spatialdata::units::Nondimensional& _getNormalizer(void) const;
+
+ /// Query databases for time dependent parameters.
+ void _queryDatabases(void);
+
+ /** 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.
+ *
+ * @param upDir Direction perpendicular to horizontal surface tangent
+ * direction that is not collinear with surface normal.
+ */
+ void _paramsLocalToGlobal(const double upDir[3]);
+
+ /** 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
Property changes on: short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
___________________________________________________________________
Name: svn:mergeinfo
+
Deleted: short/3D/PyLith/trunk/libsrc/bc/Neumann.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.icc 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.icc 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,25 +0,0 @@
-// -*- 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
-
-// Set database for boundary condition parameters.
-inline
-void
-pylith::bc::Neumann::db(spatialdata::spatialdb::SpatialDB* const db) {
- _db = db;
-}
-
-
-// End of file
Copied: short/3D/PyLith/trunk/libsrc/bc/Neumann.icc (from rev 15333, short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.icc)
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.icc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.icc 2009-06-18 23:43:10 UTC (rev 15337)
@@ -0,0 +1,35 @@
+// -*- 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
+
+#include <cassert> // USES assert()
+
+// Get label of boundary condition surface.
+inline
+const char*
+pylith::bc::Neumann::_getLabel(void) const {
+ return label();
+}
+
+// Get manager of scales used to nondimensionalize problem.
+inline
+const spatialdata::units::Nondimensional&
+pylith::bc::Neumann::_getNormalizer(void) const {
+ assert(0 != _normalizer);
+ return *_normalizer;
+}
+
+
+// End of file
Property changes on: short/3D/PyLith/trunk/libsrc/bc/Neumann.icc
___________________________________________________________________
Name: svn:mergeinfo
+
Deleted: short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.cc 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.cc 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,760 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-#include <portinfo>
-
-#include "Neumann_NEW.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 <strings.h> // USES strcasecmp()
-#include <cassert> // USES assert()
-#include <stdexcept> // USES std::runtime_error
-#include <sstream> // USES std::ostringstream
-
-//#define PRECOMPUTE_GEOMETRY
-
-// ----------------------------------------------------------------------
-typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-typedef pylith::topology::SubMesh::RealSection SubRealSection;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
-
-// ----------------------------------------------------------------------
-// Default constructor.
-pylith::bc::Neumann_NEW::Neumann_NEW(void)
-{ // constructor
-} // constructor
-
-// ----------------------------------------------------------------------
-// Destructor.
-pylith::bc::Neumann_NEW::~Neumann_NEW(void)
-{ // destructor
- deallocate();
-} // destructor
-
-// ----------------------------------------------------------------------
-// Deallocate PETSc and local data structures.
-void
-pylith::bc::Neumann_NEW::deallocate(void)
-{ // deallocate
-} // deallocate
-
-// ----------------------------------------------------------------------
-// Initialize boundary condition. Determine orienation and compute traction
-// vector at integration points.
-void
-pylith::bc::Neumann_NEW::initialize(const topology::Mesh& mesh,
- const double upDir[3])
-{ // initialize
- ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("BoundaryConditions");
-
- _queryDatabases();
-
- _paramsLocalToGlobal(upDir);
-
- logger.stagePop();
-} // initialize
-
-// ----------------------------------------------------------------------
-// Integrate contributions to residual term (r) for operator.
-void
-pylith::bc::Neumann_NEW::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
- _calculateValue(t);
- const ALE::Obj<SubRealSection>& tractionSection =
- _parameters->get("value").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_NEW::integrateJacobian(topology::Jacobian* jacobian,
- const double t,
- topology::SolutionFields* const fields)
-{ // integrateJacobian
- _needNewJacobian = false;
-} // integrateJacobian
-
-// ----------------------------------------------------------------------
-// Verify configuration is acceptable.
-void
-pylith::bc::Neumann_NEW::verifyConfiguration(const topology::Mesh& mesh) const
-{ // verifyConfiguration
- if (1 == mesh.dimension())
- throw std::runtime_error("Neumann_NEW boundary conditions are not "
- "implemented for a 1-D mesh.");
-
- BCIntegratorSubMesh::verifyConfiguration(mesh);
-} // verifyConfiguration
-
-// ----------------------------------------------------------------------
-// Get cell field for tractions.
-const pylith::topology::Field<pylith::topology::SubMesh>&
-pylith::bc::Neumann_NEW::cellField(const char* name,
- topology::SolutionFields* const fields)
-{ // cellField
- assert(0 != _parameters);
- assert(0 != name);
-
- if (0 == strcasecmp(name, "initial-value"))
- return _parameters->get("initial");
-
- else if (0 == strcasecmp(name, "rate-of-change"))
- return _parameters->get("rate");
-
- else if (0 == strcasecmp(name, "change-in-value"))
- return _parameters->get("change");
-
- else if (0 == strcasecmp(name, "rate-start-time"))
- return _parameters->get("rate time");
-
- else if (0 == strcasecmp(name, "change-start-time"))
- return _parameters->get("change time");
-
- else {
- std::ostringstream msg;
- msg << "Unknown field '" << name << "' requested for Neumann_NEW 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_NEW::_queryDatabases(void)
-{ // _queryDatabases
- assert(0 != _quadrature);
- assert(0 != _boundaryMesh);
-
- const double pressureScale = _normalizer->pressureScale();
- const double timeScale = _normalizer->timeScale();
- const double rateScale = pressureScale / timeScale;
-
- const int spaceDim = _quadrature->spaceDim();
- const int numQuadPts = _quadrature->numQuadPts();
-
- delete _parameters;
- _parameters =
- new topology::Fields<topology::Field<topology::SubMesh> >(*_boundaryMesh);
-
- // Create section to hold time dependent values
- _parameters->add("value", "traction");
- topology::Field<topology::SubMesh>& value = _parameters->get("value");
- value.scale(pressureScale);
- value.vectorFieldType(topology::FieldBase::OTHER);
- value.newSection(topology::FieldBase::CELLS_FIELD, numQuadPts*spaceDim, 1);
- value.allocate();
-
- if (0 != _dbInitial) { // Setup initial values, if provided.
- _parameters->add("initial", "initial_traction");
- topology::Field<topology::SubMesh>& initial =
- _parameters->get("initial");
- initial.cloneSection(value);
- initial.scale(pressureScale);
- initial.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
-
- _dbInitial->open();
- switch (spaceDim)
- { // switch
- case 1 : {
- const char* valueNames[] = {"traction-normal"};
- _dbInitial->queryVals(valueNames, 1);
- break;
- } // case 1
- case 2 : {
- const char* valueNames[] = {"traction-shear", "traction-normal"};
- _dbInitial->queryVals(valueNames, 2);
- break;
- } // case 2
- case 3 : {
- const char* valueNames[] = {"traction-shear-horiz",
- "traction-shear-vert",
- "traction-normal"};
- _dbInitial->queryVals(valueNames, 3);
- break;
- } // case 3
- default :
- std::cerr << "Bad spatial dimension '" << spaceDim << "'." << std::endl;
- assert(0);
- throw std::logic_error("Bad spatial dimension in Neumann_NEW.");
- } // switch
- _queryDB(&initial, _dbInitial, spaceDim, pressureScale);
- _dbInitial->close();
- } // if
-
- if (0 != _dbRate) { // Setup rate of change of values, if provided.
- _parameters->add("rate", "traction_rate");
- topology::Field<topology::SubMesh>& rate =
- _parameters->get("rate");
- rate.cloneSection(value);
- rate.scale(rateScale);
- rate.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
- const ALE::Obj<RealSection>& rateSection = rate.section();
- assert(!rateSection.isNull());
-
- _dbRate->open();
- switch (spaceDim)
- { // switch
- case 1 : {
- const char* valueNames[] = {"traction-rate-normal"};
- _dbRate->queryVals(valueNames, 1);
- break;
- } // case 1
- case 2 : {
- const char* valueNames[] = {"traction-rate-shear",
- "traction-rate-normal"};
- _dbRate->queryVals(valueNames, 2);
- break;
- } // case 2
- case 3 : {
- const char* valueNames[] = {"traction-rate-shear-horiz",
- "traction-rate-shear-vert",
- "traction-rate-normal"};
- _dbRate->queryVals(valueNames, 3);
- break;
- } // case 3
- default :
- std::cerr << "Bad spatial dimension '" << spaceDim << "'." << std::endl;
- assert(0);
- throw std::logic_error("Bad spatial dimension in Neumann_NEW.");
- } // switch
- _queryDB(&rate, _dbRate, spaceDim, rateScale);
-
- _parameters->add("rate time", "rate_traction_time");
- topology::Field<topology::SubMesh>& rateTime =
- _parameters->get("rate time");
- rateTime.newSection(rate, numQuadPts);
- rateTime.allocate();
- rateTime.scale(timeScale);
- rateTime.vectorFieldType(topology::FieldBase::MULTI_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.
- _parameters->add("change", "change_traction");
- topology::Field<topology::SubMesh>& change =
- _parameters->get("change");
- change.cloneSection(value);
- change.scale(pressureScale);
- change.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
- const ALE::Obj<RealSection>& changeSection = change.section();
- assert(!changeSection.isNull());
-
- _dbChange->open();
- switch (spaceDim)
- { // switch
- case 1 : {
- const char* valueNames[] = {"traction-normal"};
- _dbChange->queryVals(valueNames, 1);
- break;
- } // case 1
- case 2 : {
- const char* valueNames[] = {"traction-shear", "traction-normal"};
- _dbChange->queryVals(valueNames, 2);
- break;
- } // case 2
- case 3 : {
- const char* valueNames[] = {"traction-shear-horiz",
- "traction-shear-vert",
- "traction-normal"};
- _dbChange->queryVals(valueNames, 3);
- break;
- } // case 3
- default :
- std::cerr << "Bad spatial dimension '" << spaceDim << "'." << std::endl;
- assert(0);
- throw std::logic_error("Bad spatial dimension in Neumann_NEW.");
- } // switch
- _queryDB(&change, _dbChange, spaceDim, pressureScale);
-
- _parameters->add("change time", "change_traction_time");
- topology::Field<topology::SubMesh>& changeTime =
- _parameters->get("change time");
- changeTime.newSection(change, numQuadPts);
- changeTime.allocate();
- changeTime.scale(timeScale);
- changeTime.vectorFieldType(topology::FieldBase::MULTI_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_NEW::_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);
- assert(0 != _quadrature);
-
- // 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();
-
- const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
- const int numBasis = _quadrature->numBasis();
- const int numQuadPts = _quadrature->numQuadPts();
- const int spaceDim = _quadrature->spaceDim();
-
- // Containers for database query results and quadrature coordinates in
- // reference geometry.
- double_array valuesCell(numQuadPts*querySize);
- double_array quadPtsGlobal(numQuadPts*spaceDim);
-
- // Get sections.
- double_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- subSieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
-
- const ALE::Obj<RealSection>& section = field->section();
- assert(!section.isNull());
-
- const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
-
- assert(0 != _normalizer);
- const double lengthScale = _normalizer->lengthScale();
-
- // Compute quadrature information
- _quadrature->initializeGeometry();
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->computeGeometry(*_boundaryMesh, cells);
-#endif
-
- // Loop over cells in boundary mesh and perform queries.
- 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);
-
- valuesCell = 0.0;
- for (int iQuad=0, iSpace=0;
- iQuad < numQuadPts;
- ++iQuad, iSpace+=spaceDim) {
- const int err = db->query(&valuesCell[iQuad*querySize], querySize,
- &quadPtsGlobal[iSpace], spaceDim, cs);
-
- if (err) {
- std::ostringstream msg;
- msg << "Could not find 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
-
- } // for
- _normalizer->nondimensionalize(&valuesCell[0], valuesCell.size(),
- scale);
-
- // Update section
- assert(valuesCell.size() == section->getFiberDimension(*c_iter));
- section->updatePoint(*c_iter, &valuesCell[0]);
- } // for
-} // _queryDB
-
-// ----------------------------------------------------------------------
-// Initialize boundary condition. Determine orienation and compute traction
-// vector at integration points.
-void
- pylith::bc::Neumann_NEW::_paramsLocalToGlobal(const double upDir[3])
-{ // _paramsLocalToGlobal
- assert(0 != _boundaryMesh);
- assert(0 != _quadrature);
-
- double_array up(upDir, 3);
-
-
- // 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();
-
- // Quadrature related values.
- 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 = numQuadPts * spaceDim;
- double_array quadPtRef(cellDim);
- const double_array& quadPtsRef = _quadrature->quadPtsRef();
-
- // 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);
-
- // Container for cell tractions rotated to global coordinates.
- double_array initialCellLocal(fiberDim);
- double_array initialCellGlobal(fiberDim);
- double_array rateCellLocal(fiberDim);
- double_array rateCellGlobal(fiberDim);
- double_array changeCellLocal(fiberDim);
- double_array changeCellGlobal(fiberDim);
-
- // Get sections.
- double_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- subSieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
-
- 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>& changeSection = ( 0 != _dbChange) ?
- _parameters->get("change").section() : 0;
-
- // Loop over cells in boundary mesh, compute orientations, and then
- // rotate corresponding traction vector from local to global coordinates.
- 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
- // Reset traction vectors
- initialCellLocal = 0.0;
- initialCellGlobal = 0.0;
- rateCellLocal = 0.0;
- rateCellGlobal = 0.0;
- changeCellLocal = 0.0;
- changeCellGlobal = 0.0;
-
- // Get values for cell from each of the sections
- if (0 != _dbInitial) {
- assert(!initialSection.isNull());
- initialSection->restrictPoint(*c_iter, &initialCellLocal[0],
- initialCellLocal.size());
- } // if
- if (0 != _dbRate) {
- assert(!rateSection.isNull());
- rateSection->restrictPoint(*c_iter, &rateCellLocal[0],
- rateCellLocal.size());
- } // if
- if (0 != _dbChange) {
- assert(!changeSection.isNull());
- changeSection->restrictPoint(*c_iter, &changeCellLocal[0],
- changeCellLocal.size());
- } // if
-
- for(int iQuad=0, iRef=0, iSpace=0; iQuad < numQuadPts;
- ++iQuad, iRef+=cellDim, iSpace+=spaceDim) {
- // 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;
-
- if (0 != _dbInitial) {
- assert(!initialSection.isNull());
- // 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)
- initialCellGlobal[iSpace+iDim] +=
- orientation[jDim*spaceDim+iDim] * initialCellLocal[iSpace+jDim];
- } // for
- } // if
-
- if (0 != _dbRate) {
- assert(!rateSection.isNull());
- // 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)
- rateCellGlobal[iDim+iSpace] +=
- orientation[jDim*spaceDim+iDim] * rateCellLocal[iSpace+jDim];
- } // for
- } // if
-
- if (0 != _dbChange) {
- assert(!changeSection.isNull());
- // 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)
- changeCellGlobal[iDim+iSpace] +=
- orientation[jDim*spaceDim+iDim] * changeCellLocal[iSpace+jDim];
- } // for
- } // if
-
- } // for
-
- // Update sections
- if (0 != _dbInitial)
- initialSection->updatePoint(*c_iter, &initialCellGlobal[0]);
- if (0 != _dbRate)
- rateSection->updatePoint(*c_iter, &rateCellGlobal[0]);
- if (0 != _dbChange)
- changeSection->updatePoint(*c_iter, &changeCellGlobal[0]);
- } // for
-} // paramsLocalToGlobal
-
-// ----------------------------------------------------------------------
-// Calculate temporal and spatial variation of value over the list of Submesh.
-void
-pylith::bc::Neumann_NEW::_calculateValue(const double t)
-{ // _calculateValue
- assert(0 != _parameters);
- assert(0 != _boundaryMesh);
- assert(0 != _quadrature);
-
- const ALE::Obj<RealSection>& valueSection =
- _parameters->get("value").section();
- assert(!valueSection.isNull());
- valueSection->zero();
-
- 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;
-
- // 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();
-
- const int spaceDim = _quadrature->spaceDim();
- const int numBasis = _quadrature->numBasis();
- const int numQuadPts = _quadrature->numQuadPts();
- const int fiberDim = numQuadPts * spaceDim;
-
- double_array valuesCell(fiberDim);
- double_array bufferCell(fiberDim);
- double_array timeCell(numQuadPts);
- for(SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
-
- valuesCell = 0.0;
- timeCell = 0.0;
-
- // Contribution from initial value
- if (0 != _dbInitial) {
- assert(!initialSection.isNull());
- initialSection->restrictPoint(*c_iter,
- &bufferCell[0], bufferCell.size());
- valuesCell += bufferCell;
- } // if
-
- // Contribution from rate of change of value
- if (0 != _dbRate) {
- assert(!rateSection.isNull());
- assert(!rateTimeSection.isNull());
-
- rateSection->restrictPoint(*c_iter, &bufferCell[0], bufferCell.size());
- rateTimeSection->restrictPoint(*c_iter, &timeCell[0], timeCell.size());
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- if (t > timeCell[iQuad]) // rate of change integrated over time
- for (int iDim=0; iDim < spaceDim; ++iDim)
- valuesCell[iQuad*spaceDim+iDim] +=
- bufferCell[iQuad*spaceDim+iDim] * (t - timeCell[iQuad]);
- } // if
-
- // Contribution from change of value
- if (0 != _dbChange) {
- assert(!changeSection.isNull());
- assert(!changeTimeSection.isNull());
-
- changeSection->restrictPoint(*c_iter, &bufferCell[0], bufferCell.size());
- changeTimeSection->restrictPoint(*c_iter, &timeCell[0], timeCell.size());
- for (int i=0; i < numQuadPts; ++i)
- if (t >= timeCell[i]) { // change in value over time
- double scale = 1.0;
- if (0 != _dbTimeHistory) {
- double tDim = t - timeCell[i];
- _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
- for (int iDim=0; iDim < spaceDim; ++iDim)
- valuesCell[i*spaceDim+iDim] += bufferCell[i*spaceDim+iDim] * scale;
- } // if
- } // if
-
- valueSection->updateAddPoint(*c_iter, &valuesCell[0]);
- } // for
-} // _calculateValue
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.hh 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.hh 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,150 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-/** @file libsrc/bc/Neumann_NEW.hh
- *
- * @brief C++ implementation of time dependent Neumann_NEW (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_NEW ------------------------------------------------------
-class pylith::bc::Neumann_NEW : public BCIntegratorSubMesh,
- public TimeDependent
-{ // class Neumann_NEW
- friend class TestNeumann_NEW; // unit testing
-
- // PUBLIC METHODS /////////////////////////////////////////////////////
-public :
-
- /// Default constructor.
- Neumann_NEW(void);
-
- /// Destructor.
- ~Neumann_NEW(void);
-
- /// Deallocate PETSc and local data structures.
- void deallocate(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;
-
- /** Get manager of scales used to nondimensionalize problem.
- *
- * @returns Nondimensionalizer.
- */
- const spatialdata::units::Nondimensional& _getNormalizer(void) const;
-
- /// Query databases for time dependent parameters.
- void _queryDatabases(void);
-
- /** 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.
- *
- * @param upDir Direction perpendicular to horizontal surface tangent
- * direction that is not collinear with surface normal.
- */
- void _paramsLocalToGlobal(const double upDir[3]);
-
- /** 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_NEW(const Neumann_NEW&); ///< Not implemented.
- const Neumann_NEW& operator=(const Neumann_NEW&); ///< Not implemented.
-
-}; // class Neumann_NEW
-
-#include "Neumann_NEW.icc"
-
-#endif // pylith_bc_neumann_hh
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.icc 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.icc 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,35 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-#if !defined(pylith_bc_neumann_hh)
-#error "Neumann_NEW.icc can only be included from Neumann_NEW.hh"
-#endif
-
-#include <cassert> // USES assert()
-
-// Get label of boundary condition surface.
-inline
-const char*
-pylith::bc::Neumann_NEW::_getLabel(void) const {
- return label();
-}
-
-// Get manager of scales used to nondimensionalize problem.
-inline
-const spatialdata::units::Nondimensional&
-pylith::bc::Neumann_NEW::_getNormalizer(void) const {
- assert(0 != _normalizer);
- return *_normalizer;
-}
-
-
-// End of file
Copied: short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.cc (from rev 15331, short/3D/PyLith/trunk/libsrc/bc/Neumann.cc)
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.cc 2009-06-18 23:43:10 UTC (rev 15337)
@@ -0,0 +1,376 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "Neumann.hh" // implementation of object methods
+
+#include "pylith/topology/Field.hh" // HOLDSA Field
+#include "pylith/feassemble/CellGeometry.hh" // USES CellGeometry
+
+#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
+#include <cstring> // USES memcpy()
+#include <strings.h> // USES strcasecmp()
+#include <cassert> // USES assert()
+#include <stdexcept> // USES std::runtime_error
+#include <sstream> // USES std::ostringstream
+
+//#define PRECOMPUTE_GEOMETRY
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
+typedef pylith::topology::SubMesh::RealSection SubRealSection;
+typedef pylith::topology::Mesh::RealSection RealSection;
+typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::bc::Neumann::Neumann(void) :
+ _db(0)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::bc::Neumann::~Neumann(void)
+{ // destructor
+ deallocate();
+} // destructor
+
+// ----------------------------------------------------------------------
+// Deallocate PETSc and local data structures.
+void
+pylith::bc::Neumann::deallocate(void)
+{ // deallocate
+ _db = 0; // :TODO: Use shared pointer
+} // deallocate
+
+// ----------------------------------------------------------------------
+// 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
+ assert(0 != _boundaryMesh);
+ assert(0 != _quadrature);
+ assert(0 != _db);
+
+ double_array up(upDir, 3);
+ const int numCorners = _quadrature->numBasis();
+
+ assert(0 != _normalizer);
+ const double lengthScale = _normalizer->lengthScale();
+ const double pressureScale = _normalizer->pressureScale();
+
+ // 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;
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("BoundaryConditions");
+
+ delete _parameters;
+ _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();
+ traction.scale(pressureScale);
+ traction.vectorFieldType(topology::FieldBase::VECTOR);
+
+ logger.stagePop();
+
+ // 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 :
+ std::cerr << "Bad spatial dimension '" << spaceDim << "'." << std::endl;
+ assert(0);
+ throw std::logic_error("Bad spatial dimension in Neumann.");
+ } // 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();
+
+ // 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
+
+
+// End of file
Property changes on: short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.cc
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.hh (from rev 15331, short/3D/PyLith/trunk/libsrc/bc/Neumann.hh)
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.hh 2009-06-18 23:43:10 UTC (rev 15337)
@@ -0,0 +1,115 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/bc/Neumann.hh
+ *
+ * @brief C++ implementation of Neumann (prescribed tractions
+ * on a surface) boundary conditions.
+ */
+
+#if !defined(pylith_bc_neumann_hh)
+#define pylith_bc_neumann_hh
+
+// Include directives ---------------------------------------------------
+#include "BCIntegratorSubMesh.hh" // ISA BCIntegratorSubMesh
+
+// Neumann --------------------------------------------------------------
+class pylith::bc::Neumann : public BCIntegratorSubMesh
+{ // class Neumann
+ friend class TestNeumann; // unit testing
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Default constructor.
+ Neumann(void);
+
+ /// Destructor.
+ ~Neumann(void);
+
+ /// Deallocate PETSc and local data structures.
+ void deallocate(void);
+
+ /** Set database for boundary condition parameters.
+ *
+ * @param db Spatial database
+ */
+ void db(spatialdata::spatialdb::SpatialDB* const db);
+
+ /** 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 =0);
+
+ // PRIVATE MEMBERS ////////////////////////////////////////////////////
+private :
+
+ spatialdata::spatialdb::SpatialDB* _db; ///< Spatial database w/parameters
+
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ Neumann(const Neumann&); ///< Not implemented
+ const Neumann& operator=(const Neumann&); ///< Not implemented
+
+}; // class Neumann
+
+#include "Neumann.icc" // inline methods
+
+#endif // pylith_bc_neumann_hh
+
+
+// End of file
Property changes on: short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.hh
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.icc (from rev 15331, short/3D/PyLith/trunk/libsrc/bc/Neumann.icc)
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.icc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.icc 2009-06-18 23:43:10 UTC (rev 15337)
@@ -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
+
+// Set database for boundary condition parameters.
+inline
+void
+pylith::bc::Neumann::db(spatialdata::spatialdb::SpatialDB* const db) {
+ _db = db;
+}
+
+
+// End of file
Property changes on: short/3D/PyLith/trunk/libsrc/bc/Neumann_OLD.icc
___________________________________________________________________
Name: svn:mergeinfo
+
Modified: short/3D/PyLith/trunk/modulesrc/bc/Neumann.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/bc/Neumann.i 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/modulesrc/bc/Neumann.i 2009-06-18 23:43:10 UTC (rev 15337)
@@ -18,7 +18,8 @@
namespace pylith {
namespace bc {
- class Neumann : public BCIntegratorSubMesh
+ class Neumann : public BCIntegratorSubMesh,
+ public TimeDependent
{ // class Neumann
// PUBLIC METHODS /////////////////////////////////////////////////
@@ -33,12 +34,6 @@
/// Deallocate PETSc and local data structures.
void deallocate(void);
- /** Set database for boundary condition parameters.
- *
- * @param db Spatial database
- */
- void db(spatialdata::spatialdb::SpatialDB* const db);
-
/** Initialize boundary condition.
*
* @param mesh Finite-element mesh.
@@ -88,6 +83,21 @@
cellField(const char* name,
pylith::topology::SolutionFields* const fields =0);
+ // PROTECTED METHODS //////////////////////////////////////////////
+ protected :
+
+ /** Get label of boundary condition surface.
+ *
+ * @returns Label of surface (from mesh generator).
+ */
+ const char* _getLabel(void) const;
+
+ /** Get manager of scales used to nondimensionalize problem.
+ *
+ * @returns Nondimensionalizer.
+ */
+ const spatialdata::units::Nondimensional& _getNormalizer(void) const;
+
}; // class Neumann
} // bc
Modified: short/3D/PyLith/trunk/pylith/bc/Neumann.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/Neumann.py 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/pylith/bc/Neumann.py 2009-06-18 23:43:10 UTC (rev 15337)
@@ -17,11 +17,15 @@
## Factory: boundary_condition
from BoundaryCondition import BoundaryCondition
+from TimeDependent import TimeDependent
from pylith.feassemble.Integrator import Integrator
from bc import Neumann as ModuleNeumann
# Neumann class
-class Neumann(BoundaryCondition, Integrator, ModuleNeumann):
+class Neumann(BoundaryCondition,
+ TimeDependent,
+ Integrator,
+ ModuleNeumann):
"""
Python object for managing traction boundary conditions.
@@ -30,40 +34,19 @@
# INVENTORY //////////////////////////////////////////////////////////
- class Inventory(BoundaryCondition.Inventory):
- """
- Python object for managing BoundaryCondition facilities and properties.
- """
+ import pyre.inventory
+
+ from pylith.feassemble.Quadrature import SubMeshQuadrature
+ bcQuadrature = pyre.inventory.facility("quadrature",
+ factory=SubMeshQuadrature)
+ bcQuadrature.meta['tip'] = "Quadrature object for numerical integration."
+
+ from pylith.meshio.OutputNeumann import OutputNeumann
+ output = pyre.inventory.facility("output", family="output_manager",
+ factory=OutputNeumann)
+ output.meta['tip'] = "Output manager associated with diagnostic output."
- ## @class Inventory
- ## Python object for managing BoundaryCondition facilities and properties.
- ##
- ## \b Properties
- ## @li None
- ##
- ## \b Facilities
- ## @li \b quadrature Quadrature object for numerical integration
- ## @li \b db Database of boundary condition parameters
- ## @li \b output Output manager associated with diagnostic output.
- import pyre.inventory
-
- from pylith.feassemble.Quadrature import SubMeshQuadrature
- quadrature = pyre.inventory.facility("quadrature",
- factory=SubMeshQuadrature)
- quadrature.meta['tip'] = "Quadrature object for numerical integration."
-
- from spatialdata.spatialdb.SimpleDB import SimpleDB
- db = pyre.inventory.facility("db", factory=SimpleDB,
- family="spatial_database")
- db.meta['tip'] = "Database of boundary condition parameters."
-
- from pylith.meshio.OutputNeumann import OutputNeumann
- output = pyre.inventory.facility("output", family="output_manager",
- factory=OutputNeumann)
- output.meta['tip'] = "Output manager associated with diagnostic output."
-
-
# PUBLIC METHODS /////////////////////////////////////////////////////
def __init__(self, name="neumann"):
@@ -72,13 +55,19 @@
"""
BoundaryCondition.__init__(self, name)
Integrator.__init__(self)
+ TimeDependent.__init__(self)
self._loggingPrefix = "NeBC "
self.availableFields = \
{'vertex': \
{'info': [],
'data': []},
'cell': \
- {'info': ["tractions"],
+ {'info': ["initial-value",
+ "rate-of-change",
+ "rate-start-time",
+ "change-in-value",
+ "change-start-time",
+ ],
'data': []}}
return
@@ -161,8 +150,8 @@
Setup members using inventory.
"""
BoundaryCondition._configure(self)
- self.bcQuadrature = self.inventory.quadrature
- ModuleNeumann.db(self, self.inventory.db)
+ TimeDependent._configure(self)
+ self.bcQuadrature = self.inventory.bcQuadrature
self.output = self.inventory.output
return
Modified: short/3D/PyLith/trunk/pylith/meshio/OutputDirichlet.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/OutputDirichlet.py 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/pylith/meshio/OutputDirichlet.py 2009-06-18 23:43:10 UTC (rev 15337)
@@ -17,7 +17,7 @@
##
## Factory: output_manager
-from OutputManager import OutputManagerSubMesh
+from OutputManagerSubMesh import OutputManagerSubMesh
# OutputDirichlet class
class OutputDirichlet(OutputManagerSubMesh):
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am 2009-06-18 23:43:10 UTC (rev 15337)
@@ -52,7 +52,6 @@
TestDirichletBoundaryTet4.cc \
TestDirichletBoundaryHex8.cc \
TestNeumann.cc \
- TestNeumann_NEW.cc \
TestNeumannLine2.cc \
TestNeumannTri3.cc \
TestNeumannQuad4.cc \
@@ -99,7 +98,6 @@
TestDirichletBoundaryTet4.hh \
TestDirichletBoundaryHex8.hh \
TestNeumann.hh \
- TestNeumann_NEW.hh \
TestNeumannLine2.hh \
TestNeumannTri3.hh \
TestNeumannQuad4.hh \
Deleted: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,295 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-#include <portinfo>
-
-#include "TestNeumann.hh" // Implementation of class methods
-
-#include "pylith/bc/Neumann.hh" // USES Neumann
-
-#include "data/NeumannData.hh" // USES NeumannData
-
-#include "pylith/topology/Mesh.hh" // USES Mesh
-#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
-#include "pylith/topology/SubMesh.hh" // USES SubMesh
-#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
-#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
-
-#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
-#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
-#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
-
-#include <stdexcept> // USES std::runtime_erro
-
-// ----------------------------------------------------------------------
-CPPUNIT_TEST_SUITE_REGISTRATION( pylith::bc::TestNeumann );
-
-// ----------------------------------------------------------------------
-typedef pylith::topology::SubMesh::SieveMesh SieveMesh;
-typedef pylith::topology::SubMesh::RealSection RealSection;
-typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-typedef pylith::topology::SubMesh::RealSection SubRealSection;
-typedef pylith::topology::SubMesh::RestrictVisitor RestrictVisitor;
-
-// ----------------------------------------------------------------------
-// Setup testing data.
-void
-pylith::bc::TestNeumann::setUp(void)
-{ // setUp
- _data = 0;
- _quadrature = new feassemble::Quadrature<topology::SubMesh>();
- CPPUNIT_ASSERT(0 != _quadrature);
-} // setUp
-
-// ----------------------------------------------------------------------
-// Tear down testing data.
-void
-pylith::bc::TestNeumann::tearDown(void)
-{ // tearDown
- delete _data; _data = 0;
- delete _quadrature; _quadrature = 0;
-} // tearDown
-
-// ----------------------------------------------------------------------
-// Test constructor.
-void
-pylith::bc::TestNeumann::testConstructor(void)
-{ // testConstructor
- Neumann bc;
-} // testConstructor
-
-// ----------------------------------------------------------------------
-// Test db().
-void
-pylith::bc::TestNeumann::testDB(void)
-{ // testDB
- const std::string& label = "my db";
- spatialdata::spatialdb::SimpleDB db(label.c_str());
- Neumann bc;
- bc.db(&db);
-
- CPPUNIT_ASSERT(0 != bc._db);
- CPPUNIT_ASSERT_EQUAL(label, std::string(bc._db->label()));
-} // testDB
-
-// ----------------------------------------------------------------------
-// Test initialize().
-void
-pylith::bc::TestNeumann::testInitialize(void)
-{ // testInitialize
- topology::Mesh mesh;
- Neumann bc;
- topology::SolutionFields fields(mesh);
- _initialize(&mesh, &bc, &fields);
-
- CPPUNIT_ASSERT(0 != _data);
-
- const topology::SubMesh& boundaryMesh = *bc._boundaryMesh;
- const ALE::Obj<SieveSubMesh>& submesh = boundaryMesh.sieveMesh();
-
- // Check boundary mesh
- CPPUNIT_ASSERT(!submesh.isNull());
-
- const int cellDim = boundaryMesh.dimension();
- const int numCorners = _data->numCorners;
- const int spaceDim = _data->spaceDim;
- const ALE::Obj<SieveSubMesh::label_sequence>& cells = submesh->heightStratum(1);
- const int numBoundaryVertices = submesh->depthStratum(0)->size();
- const int numBoundaryCells = cells->size();
-
- CPPUNIT_ASSERT_EQUAL(_data->cellDim, cellDim);
- CPPUNIT_ASSERT_EQUAL(_data->numBoundaryVertices, numBoundaryVertices);
- CPPUNIT_ASSERT_EQUAL(_data->numBoundaryCells, numBoundaryCells);
-
- const int boundaryDepth = submesh->depth()-1; // depth of boundary cells
- const ALE::Obj<SubRealSection>& coordinates =
- submesh->getRealSection("coordinates");
- RestrictVisitor coordsVisitor(*coordinates, numCorners*spaceDim);
- // coordinates->view("Mesh coordinates from TestNeumann::testInitialize");
-
- const int numBasis = bc._quadrature->numBasis();
- const int cellVertSize = _data->numCorners * spaceDim;
- double_array cellVertices(cellVertSize);
-
- const double tolerance = 1.0e-06;
-
- // check cell vertices
- int iCell = 0;
- for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
- c_iter != cells->end();
- ++c_iter) {
- const int numCorners = submesh->getNumCellCorners(*c_iter, boundaryDepth);
- CPPUNIT_ASSERT_EQUAL(_data->numCorners, numCorners);
-
- coordsVisitor.clear();
- submesh->restrictClosure(*c_iter, coordsVisitor);
- double vert =0.0;
- double vertE =0.0;
- const double* cellVertices = coordsVisitor.getValues();
- // std::cout << "c_iter " << *c_iter << " vertex coords:" << std::endl;
- for(int iVert = 0; iVert < numCorners; ++iVert) {
- for(int iDim = 0; iDim < spaceDim; ++iDim) {
- vertE = _data->cellVertices[iDim+spaceDim*iVert+iCell*cellVertSize];
- vert = cellVertices[iDim+spaceDim*iVert];
- // std::cout << " " << cellVertices[iDim+spaceDim*iVert];
- if (fabs(vertE) > 1.0)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vert/vertE, tolerance);
- else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(vert, vertE, tolerance);
- } // for
- // std::cout << std::endl;
- } // for
- iCell++;
- } // for
-
- // Check traction values
- const int numQuadPts = _data->numQuadPts;
- const int fiberDim = numQuadPts * spaceDim;
- double_array tractionsCell(fiberDim);
- int index = 0;
- CPPUNIT_ASSERT(0 != bc._parameters);
- const ALE::Obj<SubRealSection>& tractionSection =
- bc._parameters->get("traction").section();
-
- for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
- c_iter != cells->end();
- ++c_iter) {
- tractionSection->restrictPoint(*c_iter,
- &tractionsCell[0], tractionsCell.size());
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- for (int iDim =0; iDim < spaceDim; ++iDim) {
- const double tractionsCellData = _data->tractionsCell[index];
- CPPUNIT_ASSERT_DOUBLES_EQUAL(tractionsCellData,
- tractionsCell[iQuad*spaceDim+iDim],
- tolerance);
- ++index;
- } // for
- } // for
-
-} // testInitialize
-
-// ----------------------------------------------------------------------
-// Test integrateResidual().
-void
-pylith::bc::TestNeumann::testIntegrateResidual(void)
-{ // testIntegrateResidual
- CPPUNIT_ASSERT(0 != _data);
-
- topology::Mesh mesh;
- Neumann bc;
- topology::SolutionFields fields(mesh);
- _initialize(&mesh, &bc, &fields);
-
- topology::Field<topology::Mesh>& residual = fields.get("residual");
- const double t = 0.0;
- bc.integrateResidual(residual, t, &fields);
-
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- CPPUNIT_ASSERT(!sieveMesh->depthStratum(0).isNull());
-
- const double* valsE = _data->valsResidual;
- const int totalNumVertices = sieveMesh->depthStratum(0)->size();
- const int sizeE = _data->spaceDim * totalNumVertices;
-
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
-
- const double* vals = residualSection->restrictSpace();
- const int size = residualSection->sizeWithBC();
- CPPUNIT_ASSERT_EQUAL(sizeE, size);
-
- //residual.view("RESIDUAL");
-
- const double tolerance = 1.0e-06;
- // std::cout << "computed residuals: " << std::endl;
- for (int i=0; i < size; ++i)
- // std::cout << " " << vals[i] << std::endl;
- if (fabs(valsE[i]) > 1.0)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
- else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
-} // testIntegrateResidual
-
-// ----------------------------------------------------------------------
-void
-pylith::bc::TestNeumann::_initialize(topology::Mesh* mesh,
- Neumann* const bc,
- topology::SolutionFields* fields) const
-{ // _initialize
- CPPUNIT_ASSERT(0 != _data);
- CPPUNIT_ASSERT(0 != mesh);
- CPPUNIT_ASSERT(0 != bc);
- CPPUNIT_ASSERT(0 != fields);
- CPPUNIT_ASSERT(0 != _quadrature);
-
- try {
- // Set up mesh
- meshio::MeshIOAscii iohandler;
- iohandler.filename(_data->meshFilename);
- iohandler.read(mesh);
-
- // Set up coordinates
- spatialdata::geocoords::CSCart cs;
- spatialdata::units::Nondimensional normalizer;
- cs.setSpaceDim(mesh->dimension());
- cs.initialize();
- mesh->coordsys(&cs);
- mesh->nondimensionalize(normalizer);
-
- // Set up quadrature
- _quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
- _data->basisDerivRef, _data->numQuadPts,
- _data->numBasis, _data->cellDim,
- _data->quadPts, _data->numQuadPts, _data->cellDim,
- _data->quadWts, _data->numQuadPts,
- _data->spaceDim);
-
- // Set up database
- spatialdata::spatialdb::SimpleDB db("TestNeumann");
- spatialdata::spatialdb::SimpleIOAscii dbIO;
- dbIO.filename(_data->spatialDBFilename);
- db.ioHandler(&dbIO);
- db.queryType(spatialdata::spatialdb::SimpleDB::LINEAR);
-
- const double upDir[] = { 0.0, 0.0, 1.0 };
-
- bc->quadrature(_quadrature);
- bc->label(_data->label);
- bc->db(&db);
- bc->createSubMesh(*mesh);
- bc->initialize(*mesh, upDir);
-
- // Set up fields
- CPPUNIT_ASSERT(0 != fields);
- fields->add("residual", "residual");
- fields->add("disp(t), bc(t+dt)", "displacement");
- fields->solutionName("disp(t), bc(t+dt)");
-
- topology::Field<topology::Mesh>& residual = fields->get("residual");
- const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- residual.newSection(vertices, _data->spaceDim);
- residual.allocate();
- residual.zero();
-
- fields->copyLayout("residual");
- } catch (const ALE::Exception& err) {
- throw std::runtime_error(err.msg());
- } // catch
-} // _initialize
-
-
-// End of file
Copied: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc (from rev 15333, short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.cc)
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2009-06-18 23:43:10 UTC (rev 15337)
@@ -0,0 +1,839 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestNeumann.hh" // Implementation of class methods
+
+#include "pylith/bc/Neumann.hh" // USES Neumann
+
+#include "data/NeumannData.hh" // USES NeumannData
+
+#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
+#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+
+#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
+#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
+#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
+#include "spatialdata/spatialdb/TimeHistory.hh" // USES TimeHistory
+
+#include "data/NeumannDataQuad4.hh" // USES NeumannDataQuad4
+#include "pylith/feassemble/GeometryLine2D.hh" // USES GeometryLine2D
+
+#include <stdexcept> // USES std::runtime_erro
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::bc::TestNeumann );
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::SubMesh::SieveMesh SieveMesh;
+typedef pylith::topology::SubMesh::RealSection RealSection;
+typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
+typedef pylith::topology::SubMesh::RealSection SubRealSection;
+typedef pylith::topology::SubMesh::RestrictVisitor RestrictVisitor;
+
+// ----------------------------------------------------------------------
+namespace pylith {
+ namespace bc {
+ namespace _TestNeumann {
+ const double pressureScale = 4.0;
+ const double lengthScale = 1.0; // Mesh coordinates have scale=1.0
+ const double timeScale = 0.5;
+ const int ncells = 2;
+ const int numQuadPts = 2;
+ const int spaceDim = 2;
+
+ const double initial[ncells*numQuadPts*spaceDim] = {
+ 0.3, 0.4, 0.7, 0.6,
+ 1.3, 1.4, 1.7, 1.6,
+ };
+ const double rate[ncells*numQuadPts*spaceDim] = {
+ -0.2, -0.1, 0.4, 0.3,
+ -1.2, -1.1, 1.4, 1.3,
+ };
+ const double rateTime[ncells*numQuadPts] = {
+ 0.5, 0.8,
+ 0.6, 0.9,
+ };
+ const double change[ncells*numQuadPts*spaceDim] = {
+ 1.3, 1.4, 1.7, 1.6,
+ 2.3, 2.4, 2.7, 2.6,
+ };
+ const double changeTime[ncells*numQuadPts] = {
+ 2.0, 2.4,
+ 2.1, 2.5,
+ };
+
+ const double tValue = 2.2;
+ const double valuesRate[ncells*numQuadPts*spaceDim] = {
+ -0.34, -0.17, 0.56, 0.42,
+ -1.92, -1.76, 1.82, 1.69,
+ };
+ const double valuesChange[ncells*numQuadPts*spaceDim] = {
+ 1.3, 1.4, 0.0, 0.0,
+ 2.3, 2.4, 0.0, 0.0,
+ };
+ const double valuesChangeTH[ncells*numQuadPts*spaceDim] = {
+ 1.3*0.98, 1.4*0.98, 0.0, 0.0,
+ 2.3*0.99, 2.4*0.99, 0.0, 0.0,
+ };
+
+ // Check values in section against expected values.
+ static
+ void _checkValues(const double* valuesE,
+ const int fiberDimE,
+ const topology::Field<topology::SubMesh>& field);
+ } // _TestNeumann
+ } // bc
+} // pylith
+
+// ----------------------------------------------------------------------
+// Setup testing data.
+void
+pylith::bc::TestNeumann::setUp(void)
+{ // setUp
+ _data = 0;
+ _quadrature = new feassemble::Quadrature<topology::SubMesh>();
+ CPPUNIT_ASSERT(0 != _quadrature);
+} // setUp
+
+// ----------------------------------------------------------------------
+// Tear down testing data.
+void
+pylith::bc::TestNeumann::tearDown(void)
+{ // tearDown
+ delete _data; _data = 0;
+ delete _quadrature; _quadrature = 0;
+} // tearDown
+
+// ----------------------------------------------------------------------
+// Test constructor.
+void
+pylith::bc::TestNeumann::testConstructor(void)
+{ // testConstructor
+ Neumann bc;
+} // testConstructor
+
+// ----------------------------------------------------------------------
+// Test _getLabel().
+void
+pylith::bc::TestNeumann::test_getLabel(void)
+{ // test_getLabel
+ Neumann bc;
+
+ const std::string& label = "traction bc";
+ bc.label(label.c_str());
+ CPPUNIT_ASSERT_EQUAL(label, std::string(bc._getLabel()));
+} // test_getLabel
+
+// ----------------------------------------------------------------------
+// Test initialize().
+void
+pylith::bc::TestNeumann::testInitialize(void)
+{ // testInitialize
+ topology::Mesh mesh;
+ Neumann bc;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &bc, &fields);
+
+ CPPUNIT_ASSERT(0 != _data);
+
+ const topology::SubMesh& boundaryMesh = *bc._boundaryMesh;
+ const ALE::Obj<SieveSubMesh>& submesh = boundaryMesh.sieveMesh();
+
+ // Check boundary mesh
+ CPPUNIT_ASSERT(!submesh.isNull());
+
+ const int cellDim = boundaryMesh.dimension();
+ const int numCorners = _data->numCorners;
+ const int spaceDim = _data->spaceDim;
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells = submesh->heightStratum(1);
+ const int numBoundaryVertices = submesh->depthStratum(0)->size();
+ const int numBoundaryCells = cells->size();
+
+ CPPUNIT_ASSERT_EQUAL(_data->cellDim, cellDim);
+ CPPUNIT_ASSERT_EQUAL(_data->numBoundaryVertices, numBoundaryVertices);
+ CPPUNIT_ASSERT_EQUAL(_data->numBoundaryCells, numBoundaryCells);
+
+ const int boundaryDepth = submesh->depth()-1; // depth of boundary cells
+ const ALE::Obj<SubRealSection>& coordinates =
+ submesh->getRealSection("coordinates");
+ RestrictVisitor coordsVisitor(*coordinates, numCorners*spaceDim);
+ // coordinates->view("Mesh coordinates from TestNeumann::testInitialize");
+
+ const int numBasis = bc._quadrature->numBasis();
+ const int cellVertSize = _data->numCorners * spaceDim;
+ double_array cellVertices(cellVertSize);
+
+ const double tolerance = 1.0e-06;
+
+ // check cell vertices
+ int iCell = 0;
+ for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
+ c_iter != cells->end();
+ ++c_iter) {
+ const int numCorners = submesh->getNumCellCorners(*c_iter, boundaryDepth);
+ CPPUNIT_ASSERT_EQUAL(_data->numCorners, numCorners);
+
+ coordsVisitor.clear();
+ submesh->restrictClosure(*c_iter, coordsVisitor);
+ double vert =0.0;
+ double vertE =0.0;
+ const double* cellVertices = coordsVisitor.getValues();
+ // std::cout << "c_iter " << *c_iter << " vertex coords:" << std::endl;
+ for(int iVert = 0; iVert < numCorners; ++iVert) {
+ for(int iDim = 0; iDim < spaceDim; ++iDim) {
+ vertE = _data->cellVertices[iDim+spaceDim*iVert+iCell*cellVertSize];
+ vert = cellVertices[iDim+spaceDim*iVert];
+ // std::cout << " " << cellVertices[iDim+spaceDim*iVert];
+ if (fabs(vertE) > 1.0)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vert/vertE, tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(vert, vertE, tolerance);
+ } // for
+ // std::cout << std::endl;
+ } // for
+ iCell++;
+ } // for
+
+ // Check traction values
+ const int numQuadPts = _data->numQuadPts;
+ const int fiberDim = numQuadPts * spaceDim;
+ double_array tractionsCell(fiberDim);
+ int index = 0;
+ CPPUNIT_ASSERT(0 != bc._parameters);
+ const ALE::Obj<SubRealSection>& tractionSection =
+ bc._parameters->get("initial").section();
+
+ for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
+ c_iter != cells->end();
+ ++c_iter) {
+ tractionSection->restrictPoint(*c_iter,
+ &tractionsCell[0], tractionsCell.size());
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
+ for (int iDim =0; iDim < spaceDim; ++iDim) {
+ const double tractionsCellData = _data->tractionsCell[index];
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(tractionsCellData,
+ tractionsCell[iQuad*spaceDim+iDim],
+ tolerance);
+ ++index;
+ } // for
+ } // for
+
+} // testInitialize
+
+// ----------------------------------------------------------------------
+// Test integrateResidual().
+void
+pylith::bc::TestNeumann::testIntegrateResidual(void)
+{ // testIntegrateResidual
+ CPPUNIT_ASSERT(0 != _data);
+
+ topology::Mesh mesh;
+ Neumann bc;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &bc, &fields);
+
+ topology::Field<topology::Mesh>& residual = fields.get("residual");
+ const double t = 0.0;
+ bc.integrateResidual(residual, t, &fields);
+
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ CPPUNIT_ASSERT(!sieveMesh->depthStratum(0).isNull());
+
+ const double* valsE = _data->valsResidual;
+ const int totalNumVertices = sieveMesh->depthStratum(0)->size();
+ const int sizeE = _data->spaceDim * totalNumVertices;
+
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ CPPUNIT_ASSERT(!residualSection.isNull());
+
+ const double* vals = residualSection->restrictSpace();
+ const int size = residualSection->sizeWithBC();
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
+
+ //residual.view("RESIDUAL");
+
+ const double tolerance = 1.0e-06;
+ // std::cout << "computed residuals: " << std::endl;
+ for (int i=0; i < size; ++i)
+ // std::cout << " " << vals[i] << std::endl;
+ if (fabs(valsE[i]) > 1.0)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
+} // testIntegrateResidual
+
+// ----------------------------------------------------------------------
+// Test _queryDB().
+void
+pylith::bc::TestNeumann::test_queryDB(void)
+{ // testQueryDB
+ _data = new NeumannDataQuad4();
+ feassemble::GeometryLine2D geometry;
+ CPPUNIT_ASSERT(0 != _quadrature);
+ _quadrature->refGeometry(&geometry);
+
+ topology::Mesh mesh;
+ Neumann bc;
+ _preinitialize(&mesh, &bc, true);
+
+ spatialdata::spatialdb::SimpleDB dbInitial("_TestNeumann _queryDB");
+ spatialdata::spatialdb::SimpleIOAscii dbInitialIO;
+ dbInitialIO.filename("data/quad4_traction_initial.spatialdb");
+ dbInitial.ioHandler(&dbInitialIO);
+ dbInitial.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ const double scale = _TestNeumann::pressureScale;
+ const int spaceDim = _TestNeumann::spaceDim;
+ const int numQuadPts = _TestNeumann::numQuadPts;
+ const int fiberDim = numQuadPts*spaceDim;
+ const char* queryVals[spaceDim] = { "traction-shear", "traction-normal" };
+
+ topology::Field<topology::SubMesh> initial(*bc._boundaryMesh);
+ initial.newSection(topology::FieldBase::CELLS_FIELD, fiberDim, 1);
+ initial.allocate();
+ initial.zero();
+ initial.scale(_TestNeumann::pressureScale);
+
+ dbInitial.open();
+ dbInitial.queryVals(queryVals, spaceDim);
+ bc._queryDB(&initial, &dbInitial, spaceDim, scale);
+ dbInitial.close();
+
+ const ALE::Obj<RealSection>& initialSection = initial.section();
+ CPPUNIT_ASSERT(!initialSection.isNull());
+ _TestNeumann::_checkValues(_TestNeumann::initial,
+ fiberDim, initial);
+} // testQueryDB
+
+// ----------------------------------------------------------------------
+// Test _queryDatabases().
+void
+pylith::bc::TestNeumann::test_queryDatabases(void)
+{ // test_queryDatabases
+ _data = new NeumannDataQuad4();
+ feassemble::GeometryLine2D geometry;
+ CPPUNIT_ASSERT(0 != _quadrature);
+ _quadrature->refGeometry(&geometry);
+
+ topology::Mesh mesh;
+ Neumann bc;
+ _preinitialize(&mesh, &bc, true);
+
+ spatialdata::spatialdb::SimpleDB dbInitial("_TestNeumann _queryDatabases");
+ spatialdata::spatialdb::SimpleIOAscii dbInitialIO;
+ dbInitialIO.filename("data/quad4_traction_initial.spatialdb");
+ dbInitial.ioHandler(&dbInitialIO);
+ dbInitial.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ spatialdata::spatialdb::SimpleDB dbRate("_TestNeumann _queryDatabases");
+ spatialdata::spatialdb::SimpleIOAscii dbRateIO;
+ dbRateIO.filename("data/quad4_traction_rate.spatialdb");
+ dbRate.ioHandler(&dbRateIO);
+ dbRate.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ spatialdata::spatialdb::SimpleDB dbChange("_TestNeumann _queryDatabases");
+ spatialdata::spatialdb::SimpleIOAscii dbChangeIO;
+ dbChangeIO.filename("data/quad4_traction_change.spatialdb");
+ dbChange.ioHandler(&dbChangeIO);
+ dbChange.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ spatialdata::spatialdb::TimeHistory th("_TestNeumann _queryDatabases");
+ th.filename("data/quad4_traction.timedb");
+
+ bc.dbInitial(&dbInitial);
+ bc.dbRate(&dbRate);
+ bc.dbChange(&dbChange);
+ bc.dbTimeHistory(&th);
+
+ const double pressureScale = _TestNeumann::pressureScale;
+ const double timeScale = _TestNeumann::timeScale;
+ bc._queryDatabases();
+
+ const double tolerance = 1.0e-06;
+ const int spaceDim = _TestNeumann::spaceDim;
+ const int numQuadPts = _TestNeumann::numQuadPts;
+ CPPUNIT_ASSERT(0 != bc._parameters);
+
+ // Check initial values.
+ const topology::Field<topology::SubMesh>& initial =
+ bc._parameters->get("initial");
+ _TestNeumann::_checkValues(_TestNeumann::initial,
+ numQuadPts*spaceDim, initial);
+
+ // Check rate values.
+ const topology::Field<topology::SubMesh>& rate = bc._parameters->get("rate");
+ _TestNeumann::_checkValues(_TestNeumann::rate,
+ numQuadPts*spaceDim, rate);
+
+ // Check rate start time.
+ const topology::Field<topology::SubMesh>& rateTime =
+ bc._parameters->get("rate time");
+ _TestNeumann::_checkValues(_TestNeumann::rateTime,
+ numQuadPts, rateTime);
+
+ // Check change values.
+ const topology::Field<topology::SubMesh>& change =
+ bc._parameters->get("change");
+ _TestNeumann::_checkValues(_TestNeumann::change,
+ numQuadPts*spaceDim, change);
+
+ // Check change start time.
+ const topology::Field<topology::SubMesh>& changeTime =
+ bc._parameters->get("change time");
+ _TestNeumann::_checkValues(_TestNeumann::changeTime,
+ numQuadPts*1, changeTime);
+ th.close();
+} // test_queryDatabases
+
+// ----------------------------------------------------------------------
+// Test _paramsLocalToGlobal().
+void
+pylith::bc::TestNeumann::test_paramsLocalToGlobal(void)
+{ // test_paramsLocalToGlobal
+ _data = new NeumannDataQuad4();
+ feassemble::GeometryLine2D geometry;
+ CPPUNIT_ASSERT(0 != _quadrature);
+ _quadrature->refGeometry(&geometry);
+
+ topology::Mesh mesh;
+ Neumann bc;
+ _preinitialize(&mesh, &bc, true);
+
+ spatialdata::spatialdb::SimpleDB dbInitial("_TestNeumann _queryDatabases");
+ spatialdata::spatialdb::SimpleIOAscii dbInitialIO;
+ dbInitialIO.filename("data/quad4_traction_initial.spatialdb");
+ dbInitial.ioHandler(&dbInitialIO);
+ dbInitial.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ spatialdata::spatialdb::SimpleDB dbRate("_TestNeumann _queryDatabases");
+ spatialdata::spatialdb::SimpleIOAscii dbRateIO;
+ dbRateIO.filename("data/quad4_traction_rate.spatialdb");
+ dbRate.ioHandler(&dbRateIO);
+ dbRate.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ spatialdata::spatialdb::SimpleDB dbChange("_TestNeumann _queryDatabases");
+ spatialdata::spatialdb::SimpleIOAscii dbChangeIO;
+ dbChangeIO.filename("data/quad4_traction_change.spatialdb");
+ dbChange.ioHandler(&dbChangeIO);
+ dbChange.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ bc.dbInitial(&dbInitial);
+ bc.dbRate(&dbRate);
+ bc.dbChange(&dbChange);
+
+ const double pressureScale = _TestNeumann::pressureScale;
+ const double timeScale = _TestNeumann::timeScale;
+ bc._queryDatabases();
+ const double upDir[3] = { 0.0, 0.0, 1.0 };
+ bc._paramsLocalToGlobal(upDir);
+
+ const double tolerance = 1.0e-06;
+ const int spaceDim = _TestNeumann::spaceDim;
+ const int numQuadPts = _TestNeumann::numQuadPts;
+ CPPUNIT_ASSERT(0 != bc._parameters);
+
+ // Orientation for quad4 is +x, -y for shear and normal tractions.
+ CPPUNIT_ASSERT_EQUAL(2, spaceDim);
+ const int ncells = _TestNeumann::ncells;
+ double_array valuesE(ncells*numQuadPts*spaceDim);
+
+ // Check initial values.
+ for (int i=0; i < valuesE.size(); i+=spaceDim) {
+ valuesE[i+0] = _TestNeumann::initial[i+0]; // x
+ valuesE[i+1] = -_TestNeumann::initial[i+1]; // y
+ } // for
+ const topology::Field<topology::SubMesh>& initial =
+ bc._parameters->get("initial");
+ _TestNeumann::_checkValues(&valuesE[0],
+ numQuadPts*spaceDim, initial);
+
+ // Check rate values.
+ for (int i=0; i < valuesE.size(); i+=spaceDim) {
+ valuesE[i+0] = _TestNeumann::rate[i+0]; // x
+ valuesE[i+1] = -_TestNeumann::rate[i+1]; // y
+ } // for
+ const topology::Field<topology::SubMesh>& rate = bc._parameters->get("rate");
+ _TestNeumann::_checkValues(&valuesE[0],
+ numQuadPts*spaceDim, rate);
+
+ // Check change values.
+ for (int i=0; i < valuesE.size(); i+=spaceDim) {
+ valuesE[i+0] = _TestNeumann::change[i+0]; // x
+ valuesE[i+1] = -_TestNeumann::change[i+1]; // y
+ } // for
+ const topology::Field<topology::SubMesh>& change =
+ bc._parameters->get("change");
+ _TestNeumann::_checkValues(&valuesE[0],
+ numQuadPts*spaceDim, change);
+} // test_paramsLocalToGlobal
+
+// ----------------------------------------------------------------------
+// Test _calculateValue() with initial value.
+void
+pylith::bc::TestNeumann::test_calculateValueInitial(void)
+{ // test_calculateValueInitial
+ _data = new NeumannDataQuad4();
+ feassemble::GeometryLine2D geometry;
+ CPPUNIT_ASSERT(0 != _quadrature);
+ _quadrature->refGeometry(&geometry);
+
+ topology::Mesh mesh;
+ Neumann bc;
+ _preinitialize(&mesh, &bc, true);
+
+ spatialdata::spatialdb::SimpleDB dbInitial("_TestNeumann _queryDatabases");
+ spatialdata::spatialdb::SimpleIOAscii dbInitialIO;
+ dbInitialIO.filename("data/quad4_traction_initial.spatialdb");
+ dbInitial.ioHandler(&dbInitialIO);
+ dbInitial.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ bc.dbInitial(&dbInitial);
+
+ const double timeScale = _TestNeumann::timeScale;
+ bc._queryDatabases();
+ bc._calculateValue(_TestNeumann::tValue/timeScale);
+
+ const double tolerance = 1.0e-06;
+ const int spaceDim = _TestNeumann::spaceDim;
+ const int numQuadPts = _TestNeumann::numQuadPts;
+ CPPUNIT_ASSERT(0 != bc._parameters);
+
+ // Check values.
+ const topology::Field<topology::SubMesh>& value =
+ bc._parameters->get("value");
+ _TestNeumann::_checkValues(_TestNeumann::initial,
+ numQuadPts*spaceDim, value);
+} // test_calculateValueInitial
+
+// ----------------------------------------------------------------------
+// Test _calculateValue() with rate.
+void
+pylith::bc::TestNeumann::test_calculateValueRate(void)
+{ // test_calculateValueRate
+ _data = new NeumannDataQuad4();
+ feassemble::GeometryLine2D geometry;
+ CPPUNIT_ASSERT(0 != _quadrature);
+ _quadrature->refGeometry(&geometry);
+
+ topology::Mesh mesh;
+ Neumann bc;
+ _preinitialize(&mesh, &bc, true);
+
+ spatialdata::spatialdb::SimpleDB dbRate("_TestNeumann _queryDatabases");
+ spatialdata::spatialdb::SimpleIOAscii dbRateIO;
+ dbRateIO.filename("data/quad4_traction_rate.spatialdb");
+ dbRate.ioHandler(&dbRateIO);
+ dbRate.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ bc.dbRate(&dbRate);
+
+ const double timeScale = _TestNeumann::timeScale;
+ bc._queryDatabases();
+ bc._calculateValue(_TestNeumann::tValue/timeScale);
+
+ const double tolerance = 1.0e-06;
+ const int spaceDim = _TestNeumann::spaceDim;
+ const int numQuadPts = _TestNeumann::numQuadPts;
+ CPPUNIT_ASSERT(0 != bc._parameters);
+
+ // Check values.
+ const topology::Field<topology::SubMesh>& value =
+ bc._parameters->get("value");
+ _TestNeumann::_checkValues(_TestNeumann::valuesRate,
+ numQuadPts*spaceDim, value);
+} // test_calculateValueRate
+
+// ----------------------------------------------------------------------
+// Test _calculateValue() with temporal change.
+void
+pylith::bc::TestNeumann::test_calculateValueChange(void)
+{ // test_calculateValueChange
+ _data = new NeumannDataQuad4();
+ feassemble::GeometryLine2D geometry;
+ CPPUNIT_ASSERT(0 != _quadrature);
+ _quadrature->refGeometry(&geometry);
+
+ topology::Mesh mesh;
+ Neumann bc;
+ _preinitialize(&mesh, &bc, true);
+
+ spatialdata::spatialdb::SimpleDB dbChange("_TestNeumann _queryDatabases");
+ spatialdata::spatialdb::SimpleIOAscii dbChangeIO;
+ dbChangeIO.filename("data/quad4_traction_change.spatialdb");
+ dbChange.ioHandler(&dbChangeIO);
+ dbChange.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ bc.dbChange(&dbChange);
+
+ const double timeScale = _TestNeumann::timeScale;
+ bc._queryDatabases();
+ bc._calculateValue(_TestNeumann::tValue/timeScale);
+
+ const double tolerance = 1.0e-06;
+ const int spaceDim = _TestNeumann::spaceDim;
+ const int numQuadPts = _TestNeumann::numQuadPts;
+ CPPUNIT_ASSERT(0 != bc._parameters);
+
+ // Check values.
+ const topology::Field<topology::SubMesh>& value =
+ bc._parameters->get("value");
+ _TestNeumann::_checkValues(_TestNeumann::valuesChange,
+ numQuadPts*spaceDim, value);
+} // test_calculateValueChange
+
+// ----------------------------------------------------------------------
+// Test _calculateValue() with temporal change w/time history.
+void
+pylith::bc::TestNeumann::test_calculateValueChangeTH(void)
+{ // test_calculateValueChangeTH
+ _data = new NeumannDataQuad4();
+ feassemble::GeometryLine2D geometry;
+ CPPUNIT_ASSERT(0 != _quadrature);
+ _quadrature->refGeometry(&geometry);
+
+ topology::Mesh mesh;
+ Neumann bc;
+ _preinitialize(&mesh, &bc, true);
+
+
+ spatialdata::spatialdb::SimpleDB dbChange("_TestNeumann _queryDatabases");
+ spatialdata::spatialdb::SimpleIOAscii dbChangeIO;
+ dbChangeIO.filename("data/quad4_traction_change.spatialdb");
+ dbChange.ioHandler(&dbChangeIO);
+ dbChange.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ spatialdata::spatialdb::TimeHistory th("_TestNeumann _queryDatabases");
+ th.filename("data/quad4_traction.timedb");
+
+ bc.dbChange(&dbChange);
+ bc.dbTimeHistory(&th);
+
+ const double timeScale = _TestNeumann::timeScale;
+ bc._queryDatabases();
+ bc._calculateValue(_TestNeumann::tValue/timeScale);
+
+ const double tolerance = 1.0e-06;
+ const int spaceDim = _TestNeumann::spaceDim;
+ const int numQuadPts = _TestNeumann::numQuadPts;
+ CPPUNIT_ASSERT(0 != bc._parameters);
+
+ // Check values.
+ const topology::Field<topology::SubMesh>& value =
+ bc._parameters->get("value");
+ _TestNeumann::_checkValues(_TestNeumann::valuesChangeTH,
+ numQuadPts*spaceDim, value);
+} // test_calculateValueChangeTH
+
+// ----------------------------------------------------------------------
+// Test _calculateValue() with initial, rate, and temporal change w/time history.
+void
+pylith::bc::TestNeumann::test_calculateValueAll(void)
+{ // test_calculateValueAll
+ _data = new NeumannDataQuad4();
+ feassemble::GeometryLine2D geometry;
+ CPPUNIT_ASSERT(0 != _quadrature);
+ _quadrature->refGeometry(&geometry);
+
+ topology::Mesh mesh;
+ Neumann bc;
+ _preinitialize(&mesh, &bc, true);
+
+
+ spatialdata::spatialdb::SimpleDB dbInitial("_TestNeumann _queryDatabases");
+ spatialdata::spatialdb::SimpleIOAscii dbInitialIO;
+ dbInitialIO.filename("data/quad4_traction_initial.spatialdb");
+ dbInitial.ioHandler(&dbInitialIO);
+ dbInitial.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ spatialdata::spatialdb::SimpleDB dbRate("_TestNeumann _queryDatabases");
+ spatialdata::spatialdb::SimpleIOAscii dbRateIO;
+ dbRateIO.filename("data/quad4_traction_rate.spatialdb");
+ dbRate.ioHandler(&dbRateIO);
+ dbRate.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ spatialdata::spatialdb::SimpleDB dbChange("_TestNeumann _queryDatabases");
+ spatialdata::spatialdb::SimpleIOAscii dbChangeIO;
+ dbChangeIO.filename("data/quad4_traction_change.spatialdb");
+ dbChange.ioHandler(&dbChangeIO);
+ dbChange.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
+
+ spatialdata::spatialdb::TimeHistory th("_TestNeumann _queryDatabases");
+ th.filename("data/quad4_traction.timedb");
+
+ bc.dbInitial(&dbInitial);
+ bc.dbRate(&dbRate);
+ bc.dbChange(&dbChange);
+ bc.dbTimeHistory(&th);
+
+ const double timeScale = _TestNeumann::timeScale;
+ bc._queryDatabases();
+ bc._calculateValue(_TestNeumann::tValue/timeScale);
+
+ const double tolerance = 1.0e-06;
+ const int spaceDim = _TestNeumann::spaceDim;
+ const int numQuadPts = _TestNeumann::numQuadPts;
+ CPPUNIT_ASSERT(0 != bc._parameters);
+
+ // Check values.
+ const int ncells = _TestNeumann::ncells;
+ double_array valuesE(ncells*numQuadPts*spaceDim);
+ for (int i=0; i < valuesE.size(); ++i)
+ valuesE[i] =
+ _TestNeumann::initial[i] +
+ _TestNeumann::valuesRate[i] +
+ _TestNeumann::valuesChangeTH[i];
+
+ const topology::Field<topology::SubMesh>& value =
+ bc._parameters->get("value");
+ _TestNeumann::_checkValues(&valuesE[0], numQuadPts*spaceDim, value);
+} // test_calculateValueAll
+
+// ----------------------------------------------------------------------
+void
+pylith::bc::TestNeumann::_preinitialize(topology::Mesh* mesh,
+ Neumann* const bc,
+ const bool useScales) const
+{ // _initialize
+ CPPUNIT_ASSERT(0 != _data);
+ CPPUNIT_ASSERT(0 != _quadrature);
+ CPPUNIT_ASSERT(0 != mesh);
+ CPPUNIT_ASSERT(0 != bc);
+
+ try {
+ // Set up mesh
+ meshio::MeshIOAscii iohandler;
+ iohandler.filename(_data->meshFilename);
+ iohandler.read(mesh);
+
+ // Set up coordinates
+ spatialdata::geocoords::CSCart cs;
+ cs.setSpaceDim(mesh->dimension());
+ cs.initialize();
+
+ spatialdata::units::Nondimensional normalizer;
+ if (useScales) {
+ normalizer.lengthScale(_TestNeumann::lengthScale);
+ normalizer.pressureScale(_TestNeumann::pressureScale);
+ normalizer.timeScale(_TestNeumann::timeScale);
+ } // if
+
+ mesh->coordsys(&cs);
+ mesh->nondimensionalize(normalizer);
+
+ // Set up quadrature
+ _quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
+ _data->basisDerivRef, _data->numQuadPts,
+ _data->numBasis, _data->cellDim,
+ _data->quadPts, _data->numQuadPts, _data->cellDim,
+ _data->quadWts, _data->numQuadPts,
+ _data->spaceDim);
+
+ bc->quadrature(_quadrature);
+ bc->label(_data->label);
+ bc->normalizer(normalizer);
+ bc->createSubMesh(*mesh);
+ } catch (const ALE::Exception& err) {
+ throw std::runtime_error(err.msg());
+ } // catch
+} // _preinitialize
+
+// ----------------------------------------------------------------------
+void
+pylith::bc::TestNeumann::_initialize(topology::Mesh* mesh,
+ Neumann* const bc,
+ topology::SolutionFields* fields) const
+{ // _initialize
+ CPPUNIT_ASSERT(0 != _data);
+ CPPUNIT_ASSERT(0 != mesh);
+ CPPUNIT_ASSERT(0 != bc);
+ CPPUNIT_ASSERT(0 != fields);
+ CPPUNIT_ASSERT(0 != _quadrature);
+
+ try {
+ _preinitialize(mesh, bc);
+
+ // Set up database
+ spatialdata::spatialdb::SimpleDB db("TestNeumann");
+ spatialdata::spatialdb::SimpleIOAscii dbIO;
+ dbIO.filename(_data->spatialDBFilename);
+ db.ioHandler(&dbIO);
+ db.queryType(spatialdata::spatialdb::SimpleDB::LINEAR);
+
+ const double upDir[] = { 0.0, 0.0, 1.0 };
+
+ bc->dbInitial(&db);
+ bc->initialize(*mesh, upDir);
+
+ // Set up fields
+ CPPUNIT_ASSERT(0 != fields);
+ fields->add("residual", "residual");
+ fields->add("disp(t), bc(t+dt)", "displacement");
+ fields->solutionName("disp(t), bc(t+dt)");
+
+ topology::Field<topology::Mesh>& residual = fields->get("residual");
+ residual.newSection(topology::FieldBase::VERTICES_FIELD, _data->spaceDim);
+ residual.allocate();
+ residual.zero();
+
+ fields->copyLayout("residual");
+ } catch (const ALE::Exception& err) {
+ throw std::runtime_error(err.msg());
+ } // catch
+} // _initialize
+
+// ----------------------------------------------------------------------
+// Check values in section against expected values.
+void
+pylith::bc::_TestNeumann::_checkValues(const double* valuesE,
+ const int fiberDimE,
+ const topology::Field<topology::SubMesh>& field)
+{ // _checkValues
+ assert(0 != valuesE);
+
+ const topology::SubMesh& boundaryMesh = field.mesh();
+ const ALE::Obj<SieveSubMesh>& submesh = boundaryMesh.sieveMesh();
+ CPPUNIT_ASSERT(!submesh.isNull());
+ const ALE::Obj<RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells =
+ submesh->heightStratum(1);
+
+ const double scale = field.scale();
+
+ const size_t ncells = _TestNeumann::ncells;
+ CPPUNIT_ASSERT_EQUAL(ncells, cells->size());
+
+ // Check values associated with BC.
+ int icell = 0;
+ const double tolerance = 1.0e-06;
+ for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
+ c_iter != cells->end();
+ ++c_iter, ++icell) {
+
+ const int fiberDim = section->getFiberDimension(*c_iter);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+
+ const double* values = section->restrictPoint(*c_iter);
+ for (int iDim=0; iDim < fiberDimE; ++iDim)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[icell*fiberDimE+iDim]/scale,
+ values[iDim], tolerance);
+ } // for
+} // _checkValues
+
+
+// End of file
Property changes on: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
___________________________________________________________________
Name: svn:mergeinfo
+
Deleted: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.hh 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.hh 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,95 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-/**
- * @file unittests/libtests/bc/TestNeumann.hh
- *
- * @brief C++ TestNeumann object.
- *
- * C++ unit testing for Neumann.
- */
-
-#if !defined(pylith_bc_testneumann_hh)
-#define pylith_bc_testneumann_hh
-
-#include <cppunit/extensions/HelperMacros.h>
-
-#include "pylith/bc/bcfwd.hh" // forward declarations
-#include "pylith/topology/topologyfwd.hh" // forward declarations
-#include "pylith/feassemble/feassemblefwd.hh" // forward declarations
-
-/// Namespace for pylith package
-namespace pylith {
- namespace bc {
- class TestNeumann;
- class NeumannData; // HOLDSA NeumannData
- } // bc
-} // pylith
-
-/// C++ unit testing for Neumann.
-class pylith::bc::TestNeumann : public CppUnit::TestFixture
-{ // class TestNeumann
-
- // CPPUNIT TEST SUITE /////////////////////////////////////////////////
- CPPUNIT_TEST_SUITE( TestNeumann );
-
- CPPUNIT_TEST( testConstructor );
- CPPUNIT_TEST( testDB );
-
- CPPUNIT_TEST_SUITE_END();
-
- // PUBLIC METHODS /////////////////////////////////////////////////////
-public :
-
- /// Setup testing data.
- void setUp(void);
-
- /// Tear down testing data.
- void tearDown(void);
-
- /// Test constructor.
- void testConstructor(void);
-
- /// Test db()
- void testDB(void);
-
- /// Test initialize().
- void testInitialize(void);
-
- /// Test integrateResidual().
- void testIntegrateResidual(void);
-
- // PROTECTED MEMBERS //////////////////////////////////////////////////
-protected :
-
- NeumannData* _data; ///< Data for testing
- feassemble::Quadrature<topology::SubMesh>* _quadrature; ///< Used in testing.
-
- // PRIVATE METHODS ////////////////////////////////////////////////////
-private :
-
- /** Initialize Neumann boundary condition.
- *
- * @param mesh Finite-element mesh to initialize
- * @param bc Neumann boundary condition to initialize.
- * @param fields Solution fields.
- */
- void _initialize(topology::Mesh* mesh,
- Neumann* const bc,
- topology::SolutionFields* fields) const;
-
-}; // class TestNeumann
-
-#endif // pylith_bc_neumann_hh
-
-
-// End of file
Copied: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.hh (from rev 15333, short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.hh)
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.hh 2009-06-18 23:43:10 UTC (rev 15337)
@@ -0,0 +1,141 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/bc/TestNeumann.hh
+ *
+ * @brief C++ TestNeumann object.
+ *
+ * C++ unit testing for Neumann.
+ */
+
+#if !defined(pylith_bc_testneumann_hh)
+#define pylith_bc_testneumann_hh
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "pylith/bc/bcfwd.hh" // forward declarations
+#include "pylith/topology/topologyfwd.hh" // forward declarations
+#include "pylith/feassemble/feassemblefwd.hh" // forward declarations
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace bc {
+ class TestNeumann;
+ class NeumannData; // HOLDSA NeumannData
+ } // bc
+} // pylith
+
+/// C++ unit testing for Neumann.
+class pylith::bc::TestNeumann : public CppUnit::TestFixture
+{ // class TestNeumann
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestNeumann );
+
+ CPPUNIT_TEST( testConstructor );
+ CPPUNIT_TEST( test_getLabel );
+ CPPUNIT_TEST( test_queryDB );
+ CPPUNIT_TEST( test_queryDatabases );
+ CPPUNIT_TEST( test_paramsLocalToGlobal );
+ CPPUNIT_TEST( test_calculateValueInitial );
+ CPPUNIT_TEST( test_calculateValueRate );
+ CPPUNIT_TEST( test_calculateValueChange );
+ CPPUNIT_TEST( test_calculateValueChangeTH );
+ CPPUNIT_TEST( test_calculateValueAll );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Setup testing data.
+ void setUp(void);
+
+ /// Tear down testing data.
+ void tearDown(void);
+
+ /// Test constructor.
+ void testConstructor(void);
+
+ /// Test db()
+ void testDB(void);
+
+ /// Test initialize().
+ void testInitialize(void);
+
+ /// Test integrateResidual().
+ void testIntegrateResidual(void);
+
+ /// Test _getLabel().
+ void test_getLabel(void);
+
+ /// Test _queryDB().
+ void test_queryDB(void);
+
+ /// Test _queryDatabases().
+ void test_queryDatabases(void);
+
+ /// Test _paramsLocalToGlobal().
+ void test_paramsLocalToGlobal(void);
+
+ /// Test _calculateValue() with initial value.
+ void test_calculateValueInitial(void);
+
+ /// Test _calculateValue() with rate.
+ void test_calculateValueRate(void);
+
+ /// Test _calculateValue() with temporal change.
+ void test_calculateValueChange(void);
+
+ /// Test _calculateValue() with temporal change w/time history.
+ void test_calculateValueChangeTH(void);
+
+ /// Test _calculateValue() with initial, rate, and temporal change
+ /// w/time history.
+ void test_calculateValueAll(void);
+
+ // PROTECTED MEMBERS //////////////////////////////////////////////////
+protected :
+
+ NeumannData* _data; ///< Data for testing
+ feassemble::Quadrature<topology::SubMesh>* _quadrature; ///< Used in testing.
+
+ // PRIVATE METHODS ////////////////////////////////////////////////////
+private :
+
+ /** Do minimal initialization of Neumann boundary condition.
+ *
+ * @param mesh Finite-element mesh to initialize
+ * @param bc Neumann boundary condition to initialize.
+ * @param useScales Use scales provided by local constants.
+ */
+ void _preinitialize(topology::Mesh* mesh,
+ Neumann* const bc,
+ const bool useScales =false) const;
+
+ /** Initialize Neumann boundary condition.
+ *
+ * @param mesh Finite-element mesh to initialize
+ * @param bc Neumann boundary condition to initialize.
+ * @param fields Solution fields.
+ */
+ void _initialize(topology::Mesh* mesh,
+ Neumann* const bc,
+ topology::SolutionFields* fields) const;
+
+}; // class TestNeumann
+
+#endif // pylith_bc_neumann_hh
+
+
+// End of file
Property changes on: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.hh
___________________________________________________________________
Name: svn:mergeinfo
+
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannHex8.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannHex8.hh 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannHex8.hh 2009-06-18 23:43:10 UTC (rev 15337)
@@ -35,9 +35,11 @@
{ // class TestNeumann
// CPPUNIT TEST SUITE /////////////////////////////////////////////////
- CPPUNIT_TEST_SUB_SUITE( TestNeumannHex8, TestNeumann );
+ CPPUNIT_TEST_SUITE( TestNeumannHex8 );
+
CPPUNIT_TEST( testInitialize );
CPPUNIT_TEST( testIntegrateResidual );
+
CPPUNIT_TEST_SUITE_END();
// PUBLIC METHODS /////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannLine2.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannLine2.hh 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannLine2.hh 2009-06-18 23:43:10 UTC (rev 15337)
@@ -35,9 +35,11 @@
{ // class TestNeumann
// CPPUNIT TEST SUITE /////////////////////////////////////////////////
- CPPUNIT_TEST_SUB_SUITE( TestNeumannLine2, TestNeumann );
+ CPPUNIT_TEST_SUITE( TestNeumannLine2 );
+
CPPUNIT_TEST( testInitialize );
CPPUNIT_TEST( testIntegrateResidual );
+
CPPUNIT_TEST_SUITE_END();
// PUBLIC METHODS /////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannQuad4.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannQuad4.hh 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannQuad4.hh 2009-06-18 23:43:10 UTC (rev 15337)
@@ -35,9 +35,11 @@
{ // class TestNeumann
// CPPUNIT TEST SUITE /////////////////////////////////////////////////
- CPPUNIT_TEST_SUB_SUITE( TestNeumannQuad4, TestNeumann );
+ CPPUNIT_TEST_SUITE( TestNeumannQuad4 );
+
CPPUNIT_TEST( testInitialize );
CPPUNIT_TEST( testIntegrateResidual );
+
CPPUNIT_TEST_SUITE_END();
// PUBLIC METHODS /////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannTet4.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannTet4.hh 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannTet4.hh 2009-06-18 23:43:10 UTC (rev 15337)
@@ -35,9 +35,11 @@
{ // class TestNeumann
// CPPUNIT TEST SUITE /////////////////////////////////////////////////
- CPPUNIT_TEST_SUB_SUITE( TestNeumannTet4, TestNeumann );
+ CPPUNIT_TEST_SUITE( TestNeumannTet4 );
+
CPPUNIT_TEST( testInitialize );
CPPUNIT_TEST( testIntegrateResidual );
+
CPPUNIT_TEST_SUITE_END();
// PUBLIC METHODS /////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannTri3.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannTri3.hh 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannTri3.hh 2009-06-18 23:43:10 UTC (rev 15337)
@@ -35,9 +35,11 @@
{ // class TestNeumann
// CPPUNIT TEST SUITE /////////////////////////////////////////////////
- CPPUNIT_TEST_SUB_SUITE( TestNeumannTri3, TestNeumann );
+ CPPUNIT_TEST_SUITE( TestNeumannTri3 );
+
CPPUNIT_TEST( testInitialize );
CPPUNIT_TEST( testIntegrateResidual );
+
CPPUNIT_TEST_SUITE_END();
// PUBLIC METHODS /////////////////////////////////////////////////////
Deleted: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.cc 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.cc 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,839 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-#include <portinfo>
-
-#include "TestNeumann_NEW.hh" // Implementation of class methods
-
-#include "pylith/bc/Neumann_NEW.hh" // USES Neumann_NEW
-
-#include "data/NeumannData.hh" // USES NeumannData
-
-#include "pylith/topology/Mesh.hh" // USES Mesh
-#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
-#include "pylith/topology/SubMesh.hh" // USES SubMesh
-#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
-#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
-
-#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
-#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
-#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
-#include "spatialdata/spatialdb/TimeHistory.hh" // USES TimeHistory
-
-#include "data/NeumannDataQuad4.hh" // USES NeumannDataQuad4
-#include "pylith/feassemble/GeometryLine2D.hh" // USES GeometryLine2D
-
-#include <stdexcept> // USES std::runtime_erro
-
-// ----------------------------------------------------------------------
-CPPUNIT_TEST_SUITE_REGISTRATION( pylith::bc::TestNeumann_NEW );
-
-// ----------------------------------------------------------------------
-typedef pylith::topology::SubMesh::SieveMesh SieveMesh;
-typedef pylith::topology::SubMesh::RealSection RealSection;
-typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-typedef pylith::topology::SubMesh::RealSection SubRealSection;
-typedef pylith::topology::SubMesh::RestrictVisitor RestrictVisitor;
-
-// ----------------------------------------------------------------------
-namespace pylith {
- namespace bc {
- namespace _TestNeumann_NEW {
- const double pressureScale = 4.0;
- const double lengthScale = 1.0; // Mesh coordinates have scale=1.0
- const double timeScale = 0.5;
- const int ncells = 2;
- const int numQuadPts = 2;
- const int spaceDim = 2;
-
- const double initial[ncells*numQuadPts*spaceDim] = {
- 0.3, 0.4, 0.7, 0.6,
- 1.3, 1.4, 1.7, 1.6,
- };
- const double rate[ncells*numQuadPts*spaceDim] = {
- -0.2, -0.1, 0.4, 0.3,
- -1.2, -1.1, 1.4, 1.3,
- };
- const double rateTime[ncells*numQuadPts] = {
- 0.5, 0.8,
- 0.6, 0.9,
- };
- const double change[ncells*numQuadPts*spaceDim] = {
- 1.3, 1.4, 1.7, 1.6,
- 2.3, 2.4, 2.7, 2.6,
- };
- const double changeTime[ncells*numQuadPts] = {
- 2.0, 2.4,
- 2.1, 2.5,
- };
-
- const double tValue = 2.2;
- const double valuesRate[ncells*numQuadPts*spaceDim] = {
- -0.34, -0.17, 0.56, 0.42,
- -1.92, -1.76, 1.82, 1.69,
- };
- const double valuesChange[ncells*numQuadPts*spaceDim] = {
- 1.3, 1.4, 0.0, 0.0,
- 2.3, 2.4, 0.0, 0.0,
- };
- const double valuesChangeTH[ncells*numQuadPts*spaceDim] = {
- 1.3*0.98, 1.4*0.98, 0.0, 0.0,
- 2.3*0.99, 2.4*0.99, 0.0, 0.0,
- };
-
- // Check values in section against expected values.
- static
- void _checkValues(const double* valuesE,
- const int fiberDimE,
- const topology::Field<topology::SubMesh>& field);
- } // _TestNeumann_NEW
- } // bc
-} // pylith
-
-// ----------------------------------------------------------------------
-// Setup testing data.
-void
-pylith::bc::TestNeumann_NEW::setUp(void)
-{ // setUp
- _data = 0;
- _quadrature = new feassemble::Quadrature<topology::SubMesh>();
- CPPUNIT_ASSERT(0 != _quadrature);
-} // setUp
-
-// ----------------------------------------------------------------------
-// Tear down testing data.
-void
-pylith::bc::TestNeumann_NEW::tearDown(void)
-{ // tearDown
- delete _data; _data = 0;
- delete _quadrature; _quadrature = 0;
-} // tearDown
-
-// ----------------------------------------------------------------------
-// Test constructor.
-void
-pylith::bc::TestNeumann_NEW::testConstructor(void)
-{ // testConstructor
- Neumann_NEW bc;
-} // testConstructor
-
-// ----------------------------------------------------------------------
-// Test _getLabel().
-void
-pylith::bc::TestNeumann_NEW::test_getLabel(void)
-{ // test_getLabel
- Neumann_NEW bc;
-
- const std::string& label = "traction bc";
- bc.label(label.c_str());
- CPPUNIT_ASSERT_EQUAL(label, std::string(bc._getLabel()));
-} // test_getLabel
-
-// ----------------------------------------------------------------------
-// Test initialize().
-void
-pylith::bc::TestNeumann_NEW::testInitialize(void)
-{ // testInitialize
- topology::Mesh mesh;
- Neumann_NEW bc;
- topology::SolutionFields fields(mesh);
- _initialize(&mesh, &bc, &fields);
-
- CPPUNIT_ASSERT(0 != _data);
-
- const topology::SubMesh& boundaryMesh = *bc._boundaryMesh;
- const ALE::Obj<SieveSubMesh>& submesh = boundaryMesh.sieveMesh();
-
- // Check boundary mesh
- CPPUNIT_ASSERT(!submesh.isNull());
-
- const int cellDim = boundaryMesh.dimension();
- const int numCorners = _data->numCorners;
- const int spaceDim = _data->spaceDim;
- const ALE::Obj<SieveSubMesh::label_sequence>& cells = submesh->heightStratum(1);
- const int numBoundaryVertices = submesh->depthStratum(0)->size();
- const int numBoundaryCells = cells->size();
-
- CPPUNIT_ASSERT_EQUAL(_data->cellDim, cellDim);
- CPPUNIT_ASSERT_EQUAL(_data->numBoundaryVertices, numBoundaryVertices);
- CPPUNIT_ASSERT_EQUAL(_data->numBoundaryCells, numBoundaryCells);
-
- const int boundaryDepth = submesh->depth()-1; // depth of boundary cells
- const ALE::Obj<SubRealSection>& coordinates =
- submesh->getRealSection("coordinates");
- RestrictVisitor coordsVisitor(*coordinates, numCorners*spaceDim);
- // coordinates->view("Mesh coordinates from TestNeumann_NEW::testInitialize");
-
- const int numBasis = bc._quadrature->numBasis();
- const int cellVertSize = _data->numCorners * spaceDim;
- double_array cellVertices(cellVertSize);
-
- const double tolerance = 1.0e-06;
-
- // check cell vertices
- int iCell = 0;
- for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
- c_iter != cells->end();
- ++c_iter) {
- const int numCorners = submesh->getNumCellCorners(*c_iter, boundaryDepth);
- CPPUNIT_ASSERT_EQUAL(_data->numCorners, numCorners);
-
- coordsVisitor.clear();
- submesh->restrictClosure(*c_iter, coordsVisitor);
- double vert =0.0;
- double vertE =0.0;
- const double* cellVertices = coordsVisitor.getValues();
- // std::cout << "c_iter " << *c_iter << " vertex coords:" << std::endl;
- for(int iVert = 0; iVert < numCorners; ++iVert) {
- for(int iDim = 0; iDim < spaceDim; ++iDim) {
- vertE = _data->cellVertices[iDim+spaceDim*iVert+iCell*cellVertSize];
- vert = cellVertices[iDim+spaceDim*iVert];
- // std::cout << " " << cellVertices[iDim+spaceDim*iVert];
- if (fabs(vertE) > 1.0)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vert/vertE, tolerance);
- else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(vert, vertE, tolerance);
- } // for
- // std::cout << std::endl;
- } // for
- iCell++;
- } // for
-
- // Check traction values
- const int numQuadPts = _data->numQuadPts;
- const int fiberDim = numQuadPts * spaceDim;
- double_array tractionsCell(fiberDim);
- int index = 0;
- CPPUNIT_ASSERT(0 != bc._parameters);
- const ALE::Obj<SubRealSection>& tractionSection =
- bc._parameters->get("initial").section();
-
- for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
- c_iter != cells->end();
- ++c_iter) {
- tractionSection->restrictPoint(*c_iter,
- &tractionsCell[0], tractionsCell.size());
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- for (int iDim =0; iDim < spaceDim; ++iDim) {
- const double tractionsCellData = _data->tractionsCell[index];
- CPPUNIT_ASSERT_DOUBLES_EQUAL(tractionsCellData,
- tractionsCell[iQuad*spaceDim+iDim],
- tolerance);
- ++index;
- } // for
- } // for
-
-} // testInitialize
-
-// ----------------------------------------------------------------------
-// Test integrateResidual().
-void
-pylith::bc::TestNeumann_NEW::testIntegrateResidual(void)
-{ // testIntegrateResidual
- CPPUNIT_ASSERT(0 != _data);
-
- topology::Mesh mesh;
- Neumann_NEW bc;
- topology::SolutionFields fields(mesh);
- _initialize(&mesh, &bc, &fields);
-
- topology::Field<topology::Mesh>& residual = fields.get("residual");
- const double t = 0.0;
- bc.integrateResidual(residual, t, &fields);
-
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- CPPUNIT_ASSERT(!sieveMesh->depthStratum(0).isNull());
-
- const double* valsE = _data->valsResidual;
- const int totalNumVertices = sieveMesh->depthStratum(0)->size();
- const int sizeE = _data->spaceDim * totalNumVertices;
-
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
-
- const double* vals = residualSection->restrictSpace();
- const int size = residualSection->sizeWithBC();
- CPPUNIT_ASSERT_EQUAL(sizeE, size);
-
- //residual.view("RESIDUAL");
-
- const double tolerance = 1.0e-06;
- // std::cout << "computed residuals: " << std::endl;
- for (int i=0; i < size; ++i)
- // std::cout << " " << vals[i] << std::endl;
- if (fabs(valsE[i]) > 1.0)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
- else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
-} // testIntegrateResidual
-
-// ----------------------------------------------------------------------
-// Test _queryDB().
-void
-pylith::bc::TestNeumann_NEW::test_queryDB(void)
-{ // testQueryDB
- _data = new NeumannDataQuad4();
- feassemble::GeometryLine2D geometry;
- CPPUNIT_ASSERT(0 != _quadrature);
- _quadrature->refGeometry(&geometry);
-
- topology::Mesh mesh;
- Neumann_NEW bc;
- _preinitialize(&mesh, &bc, true);
-
- spatialdata::spatialdb::SimpleDB dbInitial("_TestNeumann_NEW _queryDB");
- spatialdata::spatialdb::SimpleIOAscii dbInitialIO;
- dbInitialIO.filename("data/quad4_traction_initial.spatialdb");
- dbInitial.ioHandler(&dbInitialIO);
- dbInitial.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
-
- const double scale = _TestNeumann_NEW::pressureScale;
- const int spaceDim = _TestNeumann_NEW::spaceDim;
- const int numQuadPts = _TestNeumann_NEW::numQuadPts;
- const int fiberDim = numQuadPts*spaceDim;
- const char* queryVals[spaceDim] = { "traction-shear", "traction-normal" };
-
- topology::Field<topology::SubMesh> initial(*bc._boundaryMesh);
- initial.newSection(topology::FieldBase::CELLS_FIELD, fiberDim, 1);
- initial.allocate();
- initial.zero();
- initial.scale(_TestNeumann_NEW::pressureScale);
-
- dbInitial.open();
- dbInitial.queryVals(queryVals, spaceDim);
- bc._queryDB(&initial, &dbInitial, spaceDim, scale);
- dbInitial.close();
-
- const ALE::Obj<RealSection>& initialSection = initial.section();
- CPPUNIT_ASSERT(!initialSection.isNull());
- _TestNeumann_NEW::_checkValues(_TestNeumann_NEW::initial,
- fiberDim, initial);
-} // testQueryDB
-
-// ----------------------------------------------------------------------
-// Test _queryDatabases().
-void
-pylith::bc::TestNeumann_NEW::test_queryDatabases(void)
-{ // test_queryDatabases
- _data = new NeumannDataQuad4();
- feassemble::GeometryLine2D geometry;
- CPPUNIT_ASSERT(0 != _quadrature);
- _quadrature->refGeometry(&geometry);
-
- topology::Mesh mesh;
- Neumann_NEW bc;
- _preinitialize(&mesh, &bc, true);
-
- spatialdata::spatialdb::SimpleDB dbInitial("_TestNeumann_NEW _queryDatabases");
- spatialdata::spatialdb::SimpleIOAscii dbInitialIO;
- dbInitialIO.filename("data/quad4_traction_initial.spatialdb");
- dbInitial.ioHandler(&dbInitialIO);
- dbInitial.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
-
- spatialdata::spatialdb::SimpleDB dbRate("_TestNeumann_NEW _queryDatabases");
- spatialdata::spatialdb::SimpleIOAscii dbRateIO;
- dbRateIO.filename("data/quad4_traction_rate.spatialdb");
- dbRate.ioHandler(&dbRateIO);
- dbRate.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
-
- spatialdata::spatialdb::SimpleDB dbChange("_TestNeumann_NEW _queryDatabases");
- spatialdata::spatialdb::SimpleIOAscii dbChangeIO;
- dbChangeIO.filename("data/quad4_traction_change.spatialdb");
- dbChange.ioHandler(&dbChangeIO);
- dbChange.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
-
- spatialdata::spatialdb::TimeHistory th("_TestNeumann_NEW _queryDatabases");
- th.filename("data/quad4_traction.timedb");
-
- bc.dbInitial(&dbInitial);
- bc.dbRate(&dbRate);
- bc.dbChange(&dbChange);
- bc.dbTimeHistory(&th);
-
- const double pressureScale = _TestNeumann_NEW::pressureScale;
- const double timeScale = _TestNeumann_NEW::timeScale;
- bc._queryDatabases();
-
- const double tolerance = 1.0e-06;
- const int spaceDim = _TestNeumann_NEW::spaceDim;
- const int numQuadPts = _TestNeumann_NEW::numQuadPts;
- CPPUNIT_ASSERT(0 != bc._parameters);
-
- // Check initial values.
- const topology::Field<topology::SubMesh>& initial =
- bc._parameters->get("initial");
- _TestNeumann_NEW::_checkValues(_TestNeumann_NEW::initial,
- numQuadPts*spaceDim, initial);
-
- // Check rate values.
- const topology::Field<topology::SubMesh>& rate = bc._parameters->get("rate");
- _TestNeumann_NEW::_checkValues(_TestNeumann_NEW::rate,
- numQuadPts*spaceDim, rate);
-
- // Check rate start time.
- const topology::Field<topology::SubMesh>& rateTime =
- bc._parameters->get("rate time");
- _TestNeumann_NEW::_checkValues(_TestNeumann_NEW::rateTime,
- numQuadPts, rateTime);
-
- // Check change values.
- const topology::Field<topology::SubMesh>& change =
- bc._parameters->get("change");
- _TestNeumann_NEW::_checkValues(_TestNeumann_NEW::change,
- numQuadPts*spaceDim, change);
-
- // Check change start time.
- const topology::Field<topology::SubMesh>& changeTime =
- bc._parameters->get("change time");
- _TestNeumann_NEW::_checkValues(_TestNeumann_NEW::changeTime,
- numQuadPts*1, changeTime);
- th.close();
-} // test_queryDatabases
-
-// ----------------------------------------------------------------------
-// Test _paramsLocalToGlobal().
-void
-pylith::bc::TestNeumann_NEW::test_paramsLocalToGlobal(void)
-{ // test_paramsLocalToGlobal
- _data = new NeumannDataQuad4();
- feassemble::GeometryLine2D geometry;
- CPPUNIT_ASSERT(0 != _quadrature);
- _quadrature->refGeometry(&geometry);
-
- topology::Mesh mesh;
- Neumann_NEW bc;
- _preinitialize(&mesh, &bc, true);
-
- spatialdata::spatialdb::SimpleDB dbInitial("_TestNeumann_NEW _queryDatabases");
- spatialdata::spatialdb::SimpleIOAscii dbInitialIO;
- dbInitialIO.filename("data/quad4_traction_initial.spatialdb");
- dbInitial.ioHandler(&dbInitialIO);
- dbInitial.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
-
- spatialdata::spatialdb::SimpleDB dbRate("_TestNeumann_NEW _queryDatabases");
- spatialdata::spatialdb::SimpleIOAscii dbRateIO;
- dbRateIO.filename("data/quad4_traction_rate.spatialdb");
- dbRate.ioHandler(&dbRateIO);
- dbRate.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
-
- spatialdata::spatialdb::SimpleDB dbChange("_TestNeumann_NEW _queryDatabases");
- spatialdata::spatialdb::SimpleIOAscii dbChangeIO;
- dbChangeIO.filename("data/quad4_traction_change.spatialdb");
- dbChange.ioHandler(&dbChangeIO);
- dbChange.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
-
- bc.dbInitial(&dbInitial);
- bc.dbRate(&dbRate);
- bc.dbChange(&dbChange);
-
- const double pressureScale = _TestNeumann_NEW::pressureScale;
- const double timeScale = _TestNeumann_NEW::timeScale;
- bc._queryDatabases();
- const double upDir[3] = { 0.0, 0.0, 1.0 };
- bc._paramsLocalToGlobal(upDir);
-
- const double tolerance = 1.0e-06;
- const int spaceDim = _TestNeumann_NEW::spaceDim;
- const int numQuadPts = _TestNeumann_NEW::numQuadPts;
- CPPUNIT_ASSERT(0 != bc._parameters);
-
- // Orientation for quad4 is +x, -y for shear and normal tractions.
- CPPUNIT_ASSERT_EQUAL(2, spaceDim);
- const int ncells = _TestNeumann_NEW::ncells;
- double_array valuesE(ncells*numQuadPts*spaceDim);
-
- // Check initial values.
- for (int i=0; i < valuesE.size(); i+=spaceDim) {
- valuesE[i+0] = _TestNeumann_NEW::initial[i+0]; // x
- valuesE[i+1] = -_TestNeumann_NEW::initial[i+1]; // y
- } // for
- const topology::Field<topology::SubMesh>& initial =
- bc._parameters->get("initial");
- _TestNeumann_NEW::_checkValues(&valuesE[0],
- numQuadPts*spaceDim, initial);
-
- // Check rate values.
- for (int i=0; i < valuesE.size(); i+=spaceDim) {
- valuesE[i+0] = _TestNeumann_NEW::rate[i+0]; // x
- valuesE[i+1] = -_TestNeumann_NEW::rate[i+1]; // y
- } // for
- const topology::Field<topology::SubMesh>& rate = bc._parameters->get("rate");
- _TestNeumann_NEW::_checkValues(&valuesE[0],
- numQuadPts*spaceDim, rate);
-
- // Check change values.
- for (int i=0; i < valuesE.size(); i+=spaceDim) {
- valuesE[i+0] = _TestNeumann_NEW::change[i+0]; // x
- valuesE[i+1] = -_TestNeumann_NEW::change[i+1]; // y
- } // for
- const topology::Field<topology::SubMesh>& change =
- bc._parameters->get("change");
- _TestNeumann_NEW::_checkValues(&valuesE[0],
- numQuadPts*spaceDim, change);
-} // test_paramsLocalToGlobal
-
-// ----------------------------------------------------------------------
-// Test _calculateValue() with initial value.
-void
-pylith::bc::TestNeumann_NEW::test_calculateValueInitial(void)
-{ // test_calculateValueInitial
- _data = new NeumannDataQuad4();
- feassemble::GeometryLine2D geometry;
- CPPUNIT_ASSERT(0 != _quadrature);
- _quadrature->refGeometry(&geometry);
-
- topology::Mesh mesh;
- Neumann_NEW bc;
- _preinitialize(&mesh, &bc, true);
-
- spatialdata::spatialdb::SimpleDB dbInitial("_TestNeumann_NEW _queryDatabases");
- spatialdata::spatialdb::SimpleIOAscii dbInitialIO;
- dbInitialIO.filename("data/quad4_traction_initial.spatialdb");
- dbInitial.ioHandler(&dbInitialIO);
- dbInitial.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
-
- bc.dbInitial(&dbInitial);
-
- const double timeScale = _TestNeumann_NEW::timeScale;
- bc._queryDatabases();
- bc._calculateValue(_TestNeumann_NEW::tValue/timeScale);
-
- const double tolerance = 1.0e-06;
- const int spaceDim = _TestNeumann_NEW::spaceDim;
- const int numQuadPts = _TestNeumann_NEW::numQuadPts;
- CPPUNIT_ASSERT(0 != bc._parameters);
-
- // Check values.
- const topology::Field<topology::SubMesh>& value =
- bc._parameters->get("value");
- _TestNeumann_NEW::_checkValues(_TestNeumann_NEW::initial,
- numQuadPts*spaceDim, value);
-} // test_calculateValueInitial
-
-// ----------------------------------------------------------------------
-// Test _calculateValue() with rate.
-void
-pylith::bc::TestNeumann_NEW::test_calculateValueRate(void)
-{ // test_calculateValueRate
- _data = new NeumannDataQuad4();
- feassemble::GeometryLine2D geometry;
- CPPUNIT_ASSERT(0 != _quadrature);
- _quadrature->refGeometry(&geometry);
-
- topology::Mesh mesh;
- Neumann_NEW bc;
- _preinitialize(&mesh, &bc, true);
-
- spatialdata::spatialdb::SimpleDB dbRate("_TestNeumann_NEW _queryDatabases");
- spatialdata::spatialdb::SimpleIOAscii dbRateIO;
- dbRateIO.filename("data/quad4_traction_rate.spatialdb");
- dbRate.ioHandler(&dbRateIO);
- dbRate.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
-
- bc.dbRate(&dbRate);
-
- const double timeScale = _TestNeumann_NEW::timeScale;
- bc._queryDatabases();
- bc._calculateValue(_TestNeumann_NEW::tValue/timeScale);
-
- const double tolerance = 1.0e-06;
- const int spaceDim = _TestNeumann_NEW::spaceDim;
- const int numQuadPts = _TestNeumann_NEW::numQuadPts;
- CPPUNIT_ASSERT(0 != bc._parameters);
-
- // Check values.
- const topology::Field<topology::SubMesh>& value =
- bc._parameters->get("value");
- _TestNeumann_NEW::_checkValues(_TestNeumann_NEW::valuesRate,
- numQuadPts*spaceDim, value);
-} // test_calculateValueRate
-
-// ----------------------------------------------------------------------
-// Test _calculateValue() with temporal change.
-void
-pylith::bc::TestNeumann_NEW::test_calculateValueChange(void)
-{ // test_calculateValueChange
- _data = new NeumannDataQuad4();
- feassemble::GeometryLine2D geometry;
- CPPUNIT_ASSERT(0 != _quadrature);
- _quadrature->refGeometry(&geometry);
-
- topology::Mesh mesh;
- Neumann_NEW bc;
- _preinitialize(&mesh, &bc, true);
-
- spatialdata::spatialdb::SimpleDB dbChange("_TestNeumann_NEW _queryDatabases");
- spatialdata::spatialdb::SimpleIOAscii dbChangeIO;
- dbChangeIO.filename("data/quad4_traction_change.spatialdb");
- dbChange.ioHandler(&dbChangeIO);
- dbChange.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
-
- bc.dbChange(&dbChange);
-
- const double timeScale = _TestNeumann_NEW::timeScale;
- bc._queryDatabases();
- bc._calculateValue(_TestNeumann_NEW::tValue/timeScale);
-
- const double tolerance = 1.0e-06;
- const int spaceDim = _TestNeumann_NEW::spaceDim;
- const int numQuadPts = _TestNeumann_NEW::numQuadPts;
- CPPUNIT_ASSERT(0 != bc._parameters);
-
- // Check values.
- const topology::Field<topology::SubMesh>& value =
- bc._parameters->get("value");
- _TestNeumann_NEW::_checkValues(_TestNeumann_NEW::valuesChange,
- numQuadPts*spaceDim, value);
-} // test_calculateValueChange
-
-// ----------------------------------------------------------------------
-// Test _calculateValue() with temporal change w/time history.
-void
-pylith::bc::TestNeumann_NEW::test_calculateValueChangeTH(void)
-{ // test_calculateValueChangeTH
- _data = new NeumannDataQuad4();
- feassemble::GeometryLine2D geometry;
- CPPUNIT_ASSERT(0 != _quadrature);
- _quadrature->refGeometry(&geometry);
-
- topology::Mesh mesh;
- Neumann_NEW bc;
- _preinitialize(&mesh, &bc, true);
-
-
- spatialdata::spatialdb::SimpleDB dbChange("_TestNeumann_NEW _queryDatabases");
- spatialdata::spatialdb::SimpleIOAscii dbChangeIO;
- dbChangeIO.filename("data/quad4_traction_change.spatialdb");
- dbChange.ioHandler(&dbChangeIO);
- dbChange.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
-
- spatialdata::spatialdb::TimeHistory th("_TestNeumann_NEW _queryDatabases");
- th.filename("data/quad4_traction.timedb");
-
- bc.dbChange(&dbChange);
- bc.dbTimeHistory(&th);
-
- const double timeScale = _TestNeumann_NEW::timeScale;
- bc._queryDatabases();
- bc._calculateValue(_TestNeumann_NEW::tValue/timeScale);
-
- const double tolerance = 1.0e-06;
- const int spaceDim = _TestNeumann_NEW::spaceDim;
- const int numQuadPts = _TestNeumann_NEW::numQuadPts;
- CPPUNIT_ASSERT(0 != bc._parameters);
-
- // Check values.
- const topology::Field<topology::SubMesh>& value =
- bc._parameters->get("value");
- _TestNeumann_NEW::_checkValues(_TestNeumann_NEW::valuesChangeTH,
- numQuadPts*spaceDim, value);
-} // test_calculateValueChangeTH
-
-// ----------------------------------------------------------------------
-// Test _calculateValue() with initial, rate, and temporal change w/time history.
-void
-pylith::bc::TestNeumann_NEW::test_calculateValueAll(void)
-{ // test_calculateValueAll
- _data = new NeumannDataQuad4();
- feassemble::GeometryLine2D geometry;
- CPPUNIT_ASSERT(0 != _quadrature);
- _quadrature->refGeometry(&geometry);
-
- topology::Mesh mesh;
- Neumann_NEW bc;
- _preinitialize(&mesh, &bc, true);
-
-
- spatialdata::spatialdb::SimpleDB dbInitial("_TestNeumann_NEW _queryDatabases");
- spatialdata::spatialdb::SimpleIOAscii dbInitialIO;
- dbInitialIO.filename("data/quad4_traction_initial.spatialdb");
- dbInitial.ioHandler(&dbInitialIO);
- dbInitial.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
-
- spatialdata::spatialdb::SimpleDB dbRate("_TestNeumann_NEW _queryDatabases");
- spatialdata::spatialdb::SimpleIOAscii dbRateIO;
- dbRateIO.filename("data/quad4_traction_rate.spatialdb");
- dbRate.ioHandler(&dbRateIO);
- dbRate.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
-
- spatialdata::spatialdb::SimpleDB dbChange("_TestNeumann_NEW _queryDatabases");
- spatialdata::spatialdb::SimpleIOAscii dbChangeIO;
- dbChangeIO.filename("data/quad4_traction_change.spatialdb");
- dbChange.ioHandler(&dbChangeIO);
- dbChange.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
-
- spatialdata::spatialdb::TimeHistory th("_TestNeumann_NEW _queryDatabases");
- th.filename("data/quad4_traction.timedb");
-
- bc.dbInitial(&dbInitial);
- bc.dbRate(&dbRate);
- bc.dbChange(&dbChange);
- bc.dbTimeHistory(&th);
-
- const double timeScale = _TestNeumann_NEW::timeScale;
- bc._queryDatabases();
- bc._calculateValue(_TestNeumann_NEW::tValue/timeScale);
-
- const double tolerance = 1.0e-06;
- const int spaceDim = _TestNeumann_NEW::spaceDim;
- const int numQuadPts = _TestNeumann_NEW::numQuadPts;
- CPPUNIT_ASSERT(0 != bc._parameters);
-
- // Check values.
- const int ncells = _TestNeumann_NEW::ncells;
- double_array valuesE(ncells*numQuadPts*spaceDim);
- for (int i=0; i < valuesE.size(); ++i)
- valuesE[i] =
- _TestNeumann_NEW::initial[i] +
- _TestNeumann_NEW::valuesRate[i] +
- _TestNeumann_NEW::valuesChangeTH[i];
-
- const topology::Field<topology::SubMesh>& value =
- bc._parameters->get("value");
- _TestNeumann_NEW::_checkValues(&valuesE[0], numQuadPts*spaceDim, value);
-} // test_calculateValueAll
-
-// ----------------------------------------------------------------------
-void
-pylith::bc::TestNeumann_NEW::_preinitialize(topology::Mesh* mesh,
- Neumann_NEW* const bc,
- const bool useScales) const
-{ // _initialize
- CPPUNIT_ASSERT(0 != _data);
- CPPUNIT_ASSERT(0 != _quadrature);
- CPPUNIT_ASSERT(0 != mesh);
- CPPUNIT_ASSERT(0 != bc);
-
- try {
- // Set up mesh
- meshio::MeshIOAscii iohandler;
- iohandler.filename(_data->meshFilename);
- iohandler.read(mesh);
-
- // Set up coordinates
- spatialdata::geocoords::CSCart cs;
- cs.setSpaceDim(mesh->dimension());
- cs.initialize();
-
- spatialdata::units::Nondimensional normalizer;
- if (useScales) {
- normalizer.lengthScale(_TestNeumann_NEW::lengthScale);
- normalizer.pressureScale(_TestNeumann_NEW::pressureScale);
- normalizer.timeScale(_TestNeumann_NEW::timeScale);
- } // if
-
- mesh->coordsys(&cs);
- mesh->nondimensionalize(normalizer);
-
- // Set up quadrature
- _quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
- _data->basisDerivRef, _data->numQuadPts,
- _data->numBasis, _data->cellDim,
- _data->quadPts, _data->numQuadPts, _data->cellDim,
- _data->quadWts, _data->numQuadPts,
- _data->spaceDim);
-
- bc->quadrature(_quadrature);
- bc->label(_data->label);
- bc->normalizer(normalizer);
- bc->createSubMesh(*mesh);
- } catch (const ALE::Exception& err) {
- throw std::runtime_error(err.msg());
- } // catch
-} // _preinitialize
-
-// ----------------------------------------------------------------------
-void
-pylith::bc::TestNeumann_NEW::_initialize(topology::Mesh* mesh,
- Neumann_NEW* const bc,
- topology::SolutionFields* fields) const
-{ // _initialize
- CPPUNIT_ASSERT(0 != _data);
- CPPUNIT_ASSERT(0 != mesh);
- CPPUNIT_ASSERT(0 != bc);
- CPPUNIT_ASSERT(0 != fields);
- CPPUNIT_ASSERT(0 != _quadrature);
-
- try {
- _preinitialize(mesh, bc);
-
- // Set up database
- spatialdata::spatialdb::SimpleDB db("TestNeumann_NEW");
- spatialdata::spatialdb::SimpleIOAscii dbIO;
- dbIO.filename(_data->spatialDBFilename);
- db.ioHandler(&dbIO);
- db.queryType(spatialdata::spatialdb::SimpleDB::LINEAR);
-
- const double upDir[] = { 0.0, 0.0, 1.0 };
-
- bc->dbInitial(&db);
- bc->initialize(*mesh, upDir);
-
- // Set up fields
- CPPUNIT_ASSERT(0 != fields);
- fields->add("residual", "residual");
- fields->add("disp(t), bc(t+dt)", "displacement");
- fields->solutionName("disp(t), bc(t+dt)");
-
- topology::Field<topology::Mesh>& residual = fields->get("residual");
- residual.newSection(topology::FieldBase::VERTICES_FIELD, _data->spaceDim);
- residual.allocate();
- residual.zero();
-
- fields->copyLayout("residual");
- } catch (const ALE::Exception& err) {
- throw std::runtime_error(err.msg());
- } // catch
-} // _initialize
-
-// ----------------------------------------------------------------------
-// Check values in section against expected values.
-void
-pylith::bc::_TestNeumann_NEW::_checkValues(const double* valuesE,
- const int fiberDimE,
- const topology::Field<topology::SubMesh>& field)
-{ // _checkValues
- assert(0 != valuesE);
-
- const topology::SubMesh& boundaryMesh = field.mesh();
- const ALE::Obj<SieveSubMesh>& submesh = boundaryMesh.sieveMesh();
- CPPUNIT_ASSERT(!submesh.isNull());
- const ALE::Obj<RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- submesh->heightStratum(1);
-
- const double scale = field.scale();
-
- const size_t ncells = _TestNeumann_NEW::ncells;
- CPPUNIT_ASSERT_EQUAL(ncells, cells->size());
-
- // Check values associated with BC.
- int icell = 0;
- const double tolerance = 1.0e-06;
- for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
- c_iter != cells->end();
- ++c_iter, ++icell) {
-
- const int fiberDim = section->getFiberDimension(*c_iter);
- CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
-
- const double* values = section->restrictPoint(*c_iter);
- for (int iDim=0; iDim < fiberDimE; ++iDim)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[icell*fiberDimE+iDim]/scale,
- values[iDim], tolerance);
- } // for
-} // _checkValues
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.hh 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.hh 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,141 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-/**
- * @file unittests/libtests/bc/TestNeumann_NEW.hh
- *
- * @brief C++ TestNeumann_NEW object.
- *
- * C++ unit testing for Neumann_NEW.
- */
-
-#if !defined(pylith_bc_testneumann_hh)
-#define pylith_bc_testneumann_hh
-
-#include <cppunit/extensions/HelperMacros.h>
-
-#include "pylith/bc/bcfwd.hh" // forward declarations
-#include "pylith/topology/topologyfwd.hh" // forward declarations
-#include "pylith/feassemble/feassemblefwd.hh" // forward declarations
-
-/// Namespace for pylith package
-namespace pylith {
- namespace bc {
- class TestNeumann_NEW;
- class NeumannData; // HOLDSA NeumannData
- } // bc
-} // pylith
-
-/// C++ unit testing for Neumann_NEW.
-class pylith::bc::TestNeumann_NEW : public CppUnit::TestFixture
-{ // class TestNeumann_NEW
-
- // CPPUNIT TEST SUITE /////////////////////////////////////////////////
- CPPUNIT_TEST_SUITE( TestNeumann_NEW );
-
- CPPUNIT_TEST( testConstructor );
- CPPUNIT_TEST( test_getLabel );
- CPPUNIT_TEST( test_queryDB );
- CPPUNIT_TEST( test_queryDatabases );
- CPPUNIT_TEST( test_paramsLocalToGlobal );
- CPPUNIT_TEST( test_calculateValueInitial );
- CPPUNIT_TEST( test_calculateValueRate );
- CPPUNIT_TEST( test_calculateValueChange );
- CPPUNIT_TEST( test_calculateValueChangeTH );
- CPPUNIT_TEST( test_calculateValueAll );
-
- CPPUNIT_TEST_SUITE_END();
-
- // PUBLIC METHODS /////////////////////////////////////////////////////
-public :
-
- /// Setup testing data.
- void setUp(void);
-
- /// Tear down testing data.
- void tearDown(void);
-
- /// Test constructor.
- void testConstructor(void);
-
- /// Test db()
- void testDB(void);
-
- /// Test initialize().
- void testInitialize(void);
-
- /// Test integrateResidual().
- void testIntegrateResidual(void);
-
- /// Test _getLabel().
- void test_getLabel(void);
-
- /// Test _queryDB().
- void test_queryDB(void);
-
- /// Test _queryDatabases().
- void test_queryDatabases(void);
-
- /// Test _paramsLocalToGlobal().
- void test_paramsLocalToGlobal(void);
-
- /// Test _calculateValue() with initial value.
- void test_calculateValueInitial(void);
-
- /// Test _calculateValue() with rate.
- void test_calculateValueRate(void);
-
- /// Test _calculateValue() with temporal change.
- void test_calculateValueChange(void);
-
- /// Test _calculateValue() with temporal change w/time history.
- void test_calculateValueChangeTH(void);
-
- /// Test _calculateValue() with initial, rate, and temporal change
- /// w/time history.
- void test_calculateValueAll(void);
-
- // PROTECTED MEMBERS //////////////////////////////////////////////////
-protected :
-
- NeumannData* _data; ///< Data for testing
- feassemble::Quadrature<topology::SubMesh>* _quadrature; ///< Used in testing.
-
- // PRIVATE METHODS ////////////////////////////////////////////////////
-private :
-
- /** Do minimal initialization of Neumann_NEW boundary condition.
- *
- * @param mesh Finite-element mesh to initialize
- * @param bc Neumann_NEW boundary condition to initialize.
- * @param useScales Use scales provided by local constants.
- */
- void _preinitialize(topology::Mesh* mesh,
- Neumann_NEW* const bc,
- const bool useScales =false) const;
-
- /** Initialize Neumann_NEW boundary condition.
- *
- * @param mesh Finite-element mesh to initialize
- * @param bc Neumann_NEW boundary condition to initialize.
- * @param fields Solution fields.
- */
- void _initialize(topology::Mesh* mesh,
- Neumann_NEW* const bc,
- topology::SolutionFields* fields) const;
-
-}; // class TestNeumann_NEW
-
-#endif // pylith_bc_neumann_hh
-
-
-// End of file
Copied: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_OLD.cc (from rev 15331, short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc)
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_OLD.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_OLD.cc 2009-06-18 23:43:10 UTC (rev 15337)
@@ -0,0 +1,295 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestNeumann.hh" // Implementation of class methods
+
+#include "pylith/bc/Neumann.hh" // USES Neumann
+
+#include "data/NeumannData.hh" // USES NeumannData
+
+#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
+#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+
+#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
+#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
+#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
+
+#include <stdexcept> // USES std::runtime_erro
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::bc::TestNeumann );
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::SubMesh::SieveMesh SieveMesh;
+typedef pylith::topology::SubMesh::RealSection RealSection;
+typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
+typedef pylith::topology::SubMesh::RealSection SubRealSection;
+typedef pylith::topology::SubMesh::RestrictVisitor RestrictVisitor;
+
+// ----------------------------------------------------------------------
+// Setup testing data.
+void
+pylith::bc::TestNeumann::setUp(void)
+{ // setUp
+ _data = 0;
+ _quadrature = new feassemble::Quadrature<topology::SubMesh>();
+ CPPUNIT_ASSERT(0 != _quadrature);
+} // setUp
+
+// ----------------------------------------------------------------------
+// Tear down testing data.
+void
+pylith::bc::TestNeumann::tearDown(void)
+{ // tearDown
+ delete _data; _data = 0;
+ delete _quadrature; _quadrature = 0;
+} // tearDown
+
+// ----------------------------------------------------------------------
+// Test constructor.
+void
+pylith::bc::TestNeumann::testConstructor(void)
+{ // testConstructor
+ Neumann bc;
+} // testConstructor
+
+// ----------------------------------------------------------------------
+// Test db().
+void
+pylith::bc::TestNeumann::testDB(void)
+{ // testDB
+ const std::string& label = "my db";
+ spatialdata::spatialdb::SimpleDB db(label.c_str());
+ Neumann bc;
+ bc.db(&db);
+
+ CPPUNIT_ASSERT(0 != bc._db);
+ CPPUNIT_ASSERT_EQUAL(label, std::string(bc._db->label()));
+} // testDB
+
+// ----------------------------------------------------------------------
+// Test initialize().
+void
+pylith::bc::TestNeumann::testInitialize(void)
+{ // testInitialize
+ topology::Mesh mesh;
+ Neumann bc;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &bc, &fields);
+
+ CPPUNIT_ASSERT(0 != _data);
+
+ const topology::SubMesh& boundaryMesh = *bc._boundaryMesh;
+ const ALE::Obj<SieveSubMesh>& submesh = boundaryMesh.sieveMesh();
+
+ // Check boundary mesh
+ CPPUNIT_ASSERT(!submesh.isNull());
+
+ const int cellDim = boundaryMesh.dimension();
+ const int numCorners = _data->numCorners;
+ const int spaceDim = _data->spaceDim;
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells = submesh->heightStratum(1);
+ const int numBoundaryVertices = submesh->depthStratum(0)->size();
+ const int numBoundaryCells = cells->size();
+
+ CPPUNIT_ASSERT_EQUAL(_data->cellDim, cellDim);
+ CPPUNIT_ASSERT_EQUAL(_data->numBoundaryVertices, numBoundaryVertices);
+ CPPUNIT_ASSERT_EQUAL(_data->numBoundaryCells, numBoundaryCells);
+
+ const int boundaryDepth = submesh->depth()-1; // depth of boundary cells
+ const ALE::Obj<SubRealSection>& coordinates =
+ submesh->getRealSection("coordinates");
+ RestrictVisitor coordsVisitor(*coordinates, numCorners*spaceDim);
+ // coordinates->view("Mesh coordinates from TestNeumann::testInitialize");
+
+ const int numBasis = bc._quadrature->numBasis();
+ const int cellVertSize = _data->numCorners * spaceDim;
+ double_array cellVertices(cellVertSize);
+
+ const double tolerance = 1.0e-06;
+
+ // check cell vertices
+ int iCell = 0;
+ for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
+ c_iter != cells->end();
+ ++c_iter) {
+ const int numCorners = submesh->getNumCellCorners(*c_iter, boundaryDepth);
+ CPPUNIT_ASSERT_EQUAL(_data->numCorners, numCorners);
+
+ coordsVisitor.clear();
+ submesh->restrictClosure(*c_iter, coordsVisitor);
+ double vert =0.0;
+ double vertE =0.0;
+ const double* cellVertices = coordsVisitor.getValues();
+ // std::cout << "c_iter " << *c_iter << " vertex coords:" << std::endl;
+ for(int iVert = 0; iVert < numCorners; ++iVert) {
+ for(int iDim = 0; iDim < spaceDim; ++iDim) {
+ vertE = _data->cellVertices[iDim+spaceDim*iVert+iCell*cellVertSize];
+ vert = cellVertices[iDim+spaceDim*iVert];
+ // std::cout << " " << cellVertices[iDim+spaceDim*iVert];
+ if (fabs(vertE) > 1.0)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vert/vertE, tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(vert, vertE, tolerance);
+ } // for
+ // std::cout << std::endl;
+ } // for
+ iCell++;
+ } // for
+
+ // Check traction values
+ const int numQuadPts = _data->numQuadPts;
+ const int fiberDim = numQuadPts * spaceDim;
+ double_array tractionsCell(fiberDim);
+ int index = 0;
+ CPPUNIT_ASSERT(0 != bc._parameters);
+ const ALE::Obj<SubRealSection>& tractionSection =
+ bc._parameters->get("traction").section();
+
+ for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
+ c_iter != cells->end();
+ ++c_iter) {
+ tractionSection->restrictPoint(*c_iter,
+ &tractionsCell[0], tractionsCell.size());
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
+ for (int iDim =0; iDim < spaceDim; ++iDim) {
+ const double tractionsCellData = _data->tractionsCell[index];
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(tractionsCellData,
+ tractionsCell[iQuad*spaceDim+iDim],
+ tolerance);
+ ++index;
+ } // for
+ } // for
+
+} // testInitialize
+
+// ----------------------------------------------------------------------
+// Test integrateResidual().
+void
+pylith::bc::TestNeumann::testIntegrateResidual(void)
+{ // testIntegrateResidual
+ CPPUNIT_ASSERT(0 != _data);
+
+ topology::Mesh mesh;
+ Neumann bc;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &bc, &fields);
+
+ topology::Field<topology::Mesh>& residual = fields.get("residual");
+ const double t = 0.0;
+ bc.integrateResidual(residual, t, &fields);
+
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ CPPUNIT_ASSERT(!sieveMesh->depthStratum(0).isNull());
+
+ const double* valsE = _data->valsResidual;
+ const int totalNumVertices = sieveMesh->depthStratum(0)->size();
+ const int sizeE = _data->spaceDim * totalNumVertices;
+
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ CPPUNIT_ASSERT(!residualSection.isNull());
+
+ const double* vals = residualSection->restrictSpace();
+ const int size = residualSection->sizeWithBC();
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
+
+ //residual.view("RESIDUAL");
+
+ const double tolerance = 1.0e-06;
+ // std::cout << "computed residuals: " << std::endl;
+ for (int i=0; i < size; ++i)
+ // std::cout << " " << vals[i] << std::endl;
+ if (fabs(valsE[i]) > 1.0)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
+} // testIntegrateResidual
+
+// ----------------------------------------------------------------------
+void
+pylith::bc::TestNeumann::_initialize(topology::Mesh* mesh,
+ Neumann* const bc,
+ topology::SolutionFields* fields) const
+{ // _initialize
+ CPPUNIT_ASSERT(0 != _data);
+ CPPUNIT_ASSERT(0 != mesh);
+ CPPUNIT_ASSERT(0 != bc);
+ CPPUNIT_ASSERT(0 != fields);
+ CPPUNIT_ASSERT(0 != _quadrature);
+
+ try {
+ // Set up mesh
+ meshio::MeshIOAscii iohandler;
+ iohandler.filename(_data->meshFilename);
+ iohandler.read(mesh);
+
+ // Set up coordinates
+ spatialdata::geocoords::CSCart cs;
+ spatialdata::units::Nondimensional normalizer;
+ cs.setSpaceDim(mesh->dimension());
+ cs.initialize();
+ mesh->coordsys(&cs);
+ mesh->nondimensionalize(normalizer);
+
+ // Set up quadrature
+ _quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
+ _data->basisDerivRef, _data->numQuadPts,
+ _data->numBasis, _data->cellDim,
+ _data->quadPts, _data->numQuadPts, _data->cellDim,
+ _data->quadWts, _data->numQuadPts,
+ _data->spaceDim);
+
+ // Set up database
+ spatialdata::spatialdb::SimpleDB db("TestNeumann");
+ spatialdata::spatialdb::SimpleIOAscii dbIO;
+ dbIO.filename(_data->spatialDBFilename);
+ db.ioHandler(&dbIO);
+ db.queryType(spatialdata::spatialdb::SimpleDB::LINEAR);
+
+ const double upDir[] = { 0.0, 0.0, 1.0 };
+
+ bc->quadrature(_quadrature);
+ bc->label(_data->label);
+ bc->db(&db);
+ bc->createSubMesh(*mesh);
+ bc->initialize(*mesh, upDir);
+
+ // Set up fields
+ CPPUNIT_ASSERT(0 != fields);
+ fields->add("residual", "residual");
+ fields->add("disp(t), bc(t+dt)", "displacement");
+ fields->solutionName("disp(t), bc(t+dt)");
+
+ topology::Field<topology::Mesh>& residual = fields->get("residual");
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ residual.newSection(vertices, _data->spaceDim);
+ residual.allocate();
+ residual.zero();
+
+ fields->copyLayout("residual");
+ } catch (const ALE::Exception& err) {
+ throw std::runtime_error(err.msg());
+ } // catch
+} // _initialize
+
+
+// End of file
Property changes on: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_OLD.cc
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_OLD.hh (from rev 15331, short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.hh)
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_OLD.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_OLD.hh 2009-06-18 23:43:10 UTC (rev 15337)
@@ -0,0 +1,95 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/bc/TestNeumann.hh
+ *
+ * @brief C++ TestNeumann object.
+ *
+ * C++ unit testing for Neumann.
+ */
+
+#if !defined(pylith_bc_testneumann_hh)
+#define pylith_bc_testneumann_hh
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "pylith/bc/bcfwd.hh" // forward declarations
+#include "pylith/topology/topologyfwd.hh" // forward declarations
+#include "pylith/feassemble/feassemblefwd.hh" // forward declarations
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace bc {
+ class TestNeumann;
+ class NeumannData; // HOLDSA NeumannData
+ } // bc
+} // pylith
+
+/// C++ unit testing for Neumann.
+class pylith::bc::TestNeumann : public CppUnit::TestFixture
+{ // class TestNeumann
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestNeumann );
+
+ CPPUNIT_TEST( testConstructor );
+ CPPUNIT_TEST( testDB );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Setup testing data.
+ void setUp(void);
+
+ /// Tear down testing data.
+ void tearDown(void);
+
+ /// Test constructor.
+ void testConstructor(void);
+
+ /// Test db()
+ void testDB(void);
+
+ /// Test initialize().
+ void testInitialize(void);
+
+ /// Test integrateResidual().
+ void testIntegrateResidual(void);
+
+ // PROTECTED MEMBERS //////////////////////////////////////////////////
+protected :
+
+ NeumannData* _data; ///< Data for testing
+ feassemble::Quadrature<topology::SubMesh>* _quadrature; ///< Used in testing.
+
+ // PRIVATE METHODS ////////////////////////////////////////////////////
+private :
+
+ /** Initialize Neumann boundary condition.
+ *
+ * @param mesh Finite-element mesh to initialize
+ * @param bc Neumann boundary condition to initialize.
+ * @param fields Solution fields.
+ */
+ void _initialize(topology::Mesh* mesh,
+ Neumann* const bc,
+ topology::SolutionFields* fields) const;
+
+}; // class TestNeumann
+
+#endif // pylith_bc_neumann_hh
+
+
+// End of file
Property changes on: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_OLD.hh
___________________________________________________________________
Name: svn:mergeinfo
+
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/Makefile.am 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/Makefile.am 2009-06-18 23:43:10 UTC (rev 15337)
@@ -22,7 +22,6 @@
tri3_disp2.spatialdb \
tri3_vel2.spatialdb \
tri3_tractions.spatialdb \
- tri3_tractions_NEW.spatialdb \
tri3_force.spatialdb \
tri3_force_rate.spatialdb \
tri3_force_change.spatialdb \
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/hex8b_traction.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/hex8b_traction.spatialdb 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/hex8b_traction.spatialdb 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,7 +1,7 @@
#SPATIAL.ascii 1
SimpleDB {
num-values = 3
- value-names = horiz-shear-traction vert-shear-traction normal-traction
+ value-names = traction-shear-horiz traction-shear-vert traction-normal
value-units = Pa Pa Pa
num-locs = 1
data-dim = 0
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/line2_tractions.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/line2_tractions.spatialdb 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/line2_tractions.spatialdb 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,7 +1,7 @@
#SPATIAL.ascii 1
SimpleDB {
num-values = 1
- value-names = normal-traction
+ value-names = traction-normal
value-units = Pa
num-locs = 1
data-dim = 0
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_tractions.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_tractions.spatialdb 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_tractions.spatialdb 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,7 +1,7 @@
#SPATIAL.ascii 1
SimpleDB {
num-values = 2
- value-names = shear-traction normal-traction
+ value-names = traction-shear traction-normal
value-units = Pa Pa
num-locs = 2
data-dim = 1
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/tet4_tractions.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/tet4_tractions.spatialdb 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/tet4_tractions.spatialdb 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,7 +1,7 @@
#SPATIAL.ascii 1
SimpleDB {
num-values = 3
- value-names = horiz-shear-traction vert-shear-traction normal-traction
+ value-names = traction-shear-horiz traction-shear-vert traction-normal
value-units = Pa Pa Pa
num-locs = 1
data-dim = 0
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/tri3_tractions.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/tri3_tractions.spatialdb 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/tri3_tractions.spatialdb 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,7 +1,7 @@
#SPATIAL.ascii 1
SimpleDB {
num-values = 2
- value-names = shear-traction normal-traction
+ value-names = traction-shear traction-normal
value-units = Pa Pa
num-locs = 1
data-dim = 0
Modified: short/3D/PyLith/trunk/unittests/pytests/bc/TestDirichletBoundary.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/bc/TestDirichletBoundary.py 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/pytests/bc/TestDirichletBoundary.py 2009-06-18 23:43:10 UTC (rev 15337)
@@ -153,8 +153,13 @@
"""
from spatialdata.geocoords.CSCart import CSCart
cs = CSCart()
- cs.spaceDim = 2
+ cs.inventory.spaceDim = 2
+ cs._configure()
+ from spatialdata.units.Nondimensional import Nondimensional
+ normalizer = Nondimensional()
+ normalizer._configure()
+
from pylith.meshio.MeshIOAscii import MeshIOAscii
importer = MeshIOAscii()
importer.inventory.filename = "data/tri3.mesh"
@@ -165,46 +170,36 @@
from spatialdata.spatialdb.SimpleDB import SimpleDB
db = SimpleDB()
db.inventory.label = "TestDirichletBoundary tri3"
- db.inventory.iohandler.inventory.filename = "data/tri3.spatialdb"
+ db.inventory.iohandler.inventory.filename = "data/tri3_disp.spatialdb"
db.inventory.iohandler._configure()
db._configure()
- from pylith.bc.FixedDOFDB import FixedDOFDB
- dbRate = FixedDOFDB()
- dbRate.inventory.label = "TestDirichletBoundary rate tri3"
+ from spatialdata.spatialdb.SimpleDB import SimpleDB
+ dbRate = SimpleDB()
+ dbRate.inventory.label = "TestDirichletBoundary tri3"
+ dbRate.inventory.iohandler.inventory.filename = "data/tri3_vel.spatialdb"
+ dbRate.inventory.iohandler._configure()
dbRate._configure()
- bc.inventory.db = db
- bc.inventory.dbRate = dbRate
-
-
from pylith.bc.DirichletBoundary import DirichletBoundary
bc = DirichletBoundary()
bc.inventory.output._configure()
- bc.output.writer._configure()
- bc.label = "bc"
- bc.fixedDOF = [1]
+ bc.inventory.output.writer._configure()
+ bc.inventory.label = "bc"
+ bc.inventory.bcDOF = [1]
+ bc.inventory.dbInitial = db
+ bc.inventory.dbRate = dbRate
bc._configure()
- from pyre.units.time import second
- bc.tRef = -1.0*second
-
- from spatialdata.units.Nondimensional import Nondimensional
- normalizer = Nondimensional()
- normalizer.initialize()
-
bc.preinitialize(mesh)
bc.initialize(totalTime=0.0, numTimeSteps=1, normalizer=normalizer)
- # Setup fields
- from pylith.topology.FieldsManager import FieldsManager
- fields = FieldsManager(mesh)
- fields.addReal("field")
- fields.setFiberDimension("field", cs.spaceDim)
- fields.allocate("field")
-
- import pylith.topology.topology as bindings
- bindings.zeroRealSection(fields.getReal("field"))
+ # Setup field
+ from pylith.topology.Field import MeshField
+ field = MeshField(mesh)
+ field.newSection(field.VERTICES_FIELD, cs.spaceDim())
+ field.allocate()
+ field.zero()
return (mesh, bc, field)
Modified: short/3D/PyLith/trunk/unittests/pytests/bc/TestNeumann.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/bc/TestNeumann.py 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/pytests/bc/TestNeumann.py 2009-06-18 23:43:10 UTC (rev 15337)
@@ -201,8 +201,8 @@
from pylith.bc.Neumann import Neumann
bc = Neumann()
- bc.inventory.quadrature = quadrature
- bc.inventory.db = db
+ bc.inventory.bcQuadrature = quadrature
+ bc.inventory.dbInitial = db
bc.inventory.label = "bc"
bc.inventory.output.inventory.writer._configure()
bc.inventory.output._configure()
Modified: short/3D/PyLith/trunk/unittests/pytests/bc/data/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/bc/data/Makefile.am 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/pytests/bc/data/Makefile.am 2009-06-18 23:43:10 UTC (rev 15337)
@@ -12,6 +12,7 @@
dist_noinst_DATA = \
tri3_disp.spatialdb \
+ tri3_vel.spatialdb \
tri3_tractions.spatialdb \
tri3.mesh \
elasticplanestrain.spatialdb
Modified: short/3D/PyLith/trunk/unittests/pytests/bc/data/tri3_tractions.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/bc/data/tri3_tractions.spatialdb 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/pytests/bc/data/tri3_tractions.spatialdb 2009-06-18 23:43:10 UTC (rev 15337)
@@ -1,7 +1,7 @@
#SPATIAL.ascii 1
SimpleDB {
num-values = 2
- value-names = shear-traction normal-traction
+ value-names = traction-shear traction-normal
value-units = Pa Pa
num-locs = 1
data-dim = 0
Added: short/3D/PyLith/trunk/unittests/pytests/bc/data/tri3_vel.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/bc/data/tri3_vel.spatialdb (rev 0)
+++ short/3D/PyLith/trunk/unittests/pytests/bc/data/tri3_vel.spatialdb 2009-06-18 23:43:10 UTC (rev 15337)
@@ -0,0 +1,15 @@
+#SPATIAL.ascii 1
+SimpleDB {
+ num-values = 2
+ value-names = displacement-rate-y rate-start-time
+ value-units = m s
+ num-locs = 2
+ data-dim = 1
+ space-dim = 2
+ cs-data = cartesian {
+ to-meters = 1.0
+ space-dim = 2
+ }
+}
+ 0.0 -1.0 0.3 0.0
+ 1.0 0.0 0.7 0.1
Modified: short/3D/PyLith/trunk/unittests/pytests/bc/testbc.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/bc/testbc.py 2009-06-18 22:15:58 UTC (rev 15336)
+++ short/3D/PyLith/trunk/unittests/pytests/bc/testbc.py 2009-06-18 23:43:10 UTC (rev 15337)
@@ -59,8 +59,8 @@
from TestDirichletBC import TestDirichletBC
suite.addTest(unittest.makeSuite(TestDirichletBC))
- #from TestDirichletBoundary import TestDirichletBoundary
- #suite.addTest(unittest.makeSuite(TestDirichletBoundary))
+ from TestDirichletBoundary import TestDirichletBoundary
+ suite.addTest(unittest.makeSuite(TestDirichletBoundary))
from TestAbsorbingDampers import TestAbsorbingDampers
suite.addTest(unittest.makeSuite(TestAbsorbingDampers))
More information about the CIG-COMMITS
mailing list