[cig-commits] r15333 - in short/3D/PyLith/trunk: . libsrc/bc libsrc/topology unittests/libtests/bc unittests/libtests/bc/data unittests/libtests/faults
brad at geodynamics.org
brad at geodynamics.org
Thu Jun 18 13:54:34 PDT 2009
Author: brad
Date: 2009-06-18 13:54:31 -0700 (Thu, 18 Jun 2009)
New Revision: 15333
Added:
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.hh
short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction.timedb
short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction_change.spatialdb
short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction_initial.spatialdb
short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction_rate.spatialdb
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/bc/Makefile.am
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/libsrc/topology/SubMesh.hh
short/3D/PyLith/trunk/libsrc/topology/SubMesh.icc
short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/bc/data/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
Log:
Fixed bugs in new time-dependent Neumann BC. Added corresponding C++ unit tests.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2009-06-18 18:56:56 UTC (rev 15332)
+++ short/3D/PyLith/trunk/TODO 2009-06-18 20:54:31 UTC (rev 15333)
@@ -7,8 +7,6 @@
Need power-law full-scale "test" to insure that the power-law and
nonlinear solve give reasonable results.
- Need generalized Maxwell full-scale "test".
-
Time dependent Neumann BC
Run valgrind on the unit tests.
Modified: short/3D/PyLith/trunk/libsrc/bc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Makefile.am 2009-06-18 18:56:56 UTC (rev 15332)
+++ short/3D/PyLith/trunk/libsrc/bc/Makefile.am 2009-06-18 20:54:31 UTC (rev 15333)
@@ -31,6 +31,8 @@
DirichletBoundary.icc \
Neumann.hh \
Neumann.icc \
+ Neumann_NEW.hh \
+ Neumann_NEW.icc \
PointForce.hh \
PointForce.icc \
bcfwd.hh
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.cc 2009-06-18 18:56:56 UTC (rev 15332)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.cc 2009-06-18 20:54:31 UTC (rev 15333)
@@ -38,8 +38,7 @@
// ----------------------------------------------------------------------
// Default constructor.
-pylith::bc::Neumann_NEW::Neumann_NEW(void) :
- _db(0)
+pylith::bc::Neumann_NEW::Neumann_NEW(void)
{ // constructor
} // constructor
@@ -55,7 +54,6 @@
void
pylith::bc::Neumann_NEW::deallocate(void)
{ // deallocate
- _db = 0; // :TODO: Use shared pointer
} // deallocate
// ----------------------------------------------------------------------
@@ -198,12 +196,22 @@
assert(0 != _parameters);
assert(0 != name);
- if (0 == strcasecmp(name, "tractions")) {
- return _parameters->get("traction");
+ if (0 == strcasecmp(name, "initial-value"))
+ return _parameters->get("initial");
- // ADD STUFF HERE
+ else if (0 == strcasecmp(name, "rate-of-change"))
+ return _parameters->get("rate");
- } else {
+ 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 << "'.";
@@ -227,7 +235,6 @@
const int spaceDim = _quadrature->spaceDim();
const int numQuadPts = _quadrature->numQuadPts();
- const int fiberDim = numQuadPts * spaceDim;
delete _parameters;
_parameters =
@@ -238,7 +245,7 @@
topology::Field<topology::SubMesh>& value = _parameters->get("value");
value.scale(pressureScale);
value.vectorFieldType(topology::FieldBase::OTHER);
- value.newSection(topology::FieldBase::CELLS_FIELD, fiberDim);
+ value.newSection(topology::FieldBase::CELLS_FIELD, numQuadPts*spaceDim, 1);
value.allocate();
if (0 != _dbInitial) { // Setup initial values, if provided.
@@ -314,15 +321,15 @@
assert(0);
throw std::logic_error("Bad spatial dimension in Neumann_NEW.");
} // switch
- _queryDB(&rate, _dbRate, spaceDim, pressureScale);
+ _queryDB(&rate, _dbRate, spaceDim, rateScale);
_parameters->add("rate time", "rate_traction_time");
topology::Field<topology::SubMesh>& rateTime =
_parameters->get("rate time");
- rateTime.newSection(rate, 1);
+ rateTime.newSection(rate, numQuadPts);
rateTime.allocate();
rateTime.scale(timeScale);
- rateTime.vectorFieldType(topology::FieldBase::SCALAR);
+ rateTime.vectorFieldType(topology::FieldBase::MULTI_SCALAR);
const char* timeNames[1] = { "rate-start-time" };
_dbRate->queryVals(timeNames, 1);
@@ -370,10 +377,10 @@
_parameters->add("change time", "change_traction_time");
topology::Field<topology::SubMesh>& changeTime =
_parameters->get("change time");
- changeTime.newSection(change, 1);
+ changeTime.newSection(change, numQuadPts);
changeTime.allocate();
changeTime.scale(timeScale);
- changeTime.vectorFieldType(topology::FieldBase::SCALAR);
+ changeTime.vectorFieldType(topology::FieldBase::MULTI_SCALAR);
const char* timeNames[1] = { "change-start-time" };
_dbChange->queryVals(timeNames, 1);
@@ -412,11 +419,10 @@
const int numBasis = _quadrature->numBasis();
const int numQuadPts = _quadrature->numQuadPts();
const int spaceDim = _quadrature->spaceDim();
- const int fiberDim = numQuadPts * querySize;
// Containers for database query results and quadrature coordinates in
// reference geometry.
- double_array valuesCell(querySize);
+ double_array valuesCell(numQuadPts*querySize);
double_array quadPtRef(cellDim);
double_array quadPtsGlobal(numQuadPts*spaceDim);
const double_array& quadPtsRef = _quadrature->quadPtsRef();
@@ -462,25 +468,28 @@
lengthScale);
valuesCell = 0.0;
- for(int iQuad=0, iRef=0, iSpace=0; iQuad < numQuadPts;
- ++iQuad, iRef+=cellDim, iSpace+=spaceDim) {
- const int err = _db->query(&valuesCell[iQuad*querySize], querySize,
- &quadPtsGlobal[iSpace], spaceDim, cs);
+ 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() << ".";
+ << "using spatial database " << db->label() << ".";
throw std::runtime_error(msg.str());
} // if
- _normalizer->nondimensionalize(&valuesCell[0], valuesCell.size(),
- scale);
} // 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
@@ -512,7 +521,7 @@
const int numBasis = _quadrature->numBasis();
const int numQuadPts = _quadrature->numQuadPts();
const int spaceDim = cellGeometry.spaceDim();
- const int fiberDim = spaceDim * numQuadPts;
+ const int fiberDim = numQuadPts * spaceDim;
double_array quadPtRef(cellDim);
const double_array& quadPtsRef = _quadrature->quadPtsRef();
@@ -606,8 +615,8 @@
// coordinate system
for(int iDim = 0; iDim < spaceDim; ++iDim) {
for(int jDim = 0; jDim < spaceDim; ++jDim)
- initialCellGlobal[iDim+iSpace] +=
- orientation[jDim*spaceDim+iDim] * initialCellLocal[jDim];
+ initialCellGlobal[iSpace+iDim] +=
+ orientation[jDim*spaceDim+iDim] * initialCellLocal[iSpace+jDim];
} // for
} // if
@@ -618,7 +627,7 @@
for(int iDim = 0; iDim < spaceDim; ++iDim) {
for(int jDim = 0; jDim < spaceDim; ++jDim)
rateCellGlobal[iDim+iSpace] +=
- orientation[jDim*spaceDim+iDim] * rateCellLocal[jDim];
+ orientation[jDim*spaceDim+iDim] * rateCellLocal[iSpace+jDim];
} // for
} // if
@@ -629,7 +638,7 @@
for(int iDim = 0; iDim < spaceDim; ++iDim) {
for(int jDim = 0; jDim < spaceDim; ++jDim)
changeCellGlobal[iDim+iSpace] +=
- orientation[jDim*spaceDim+iDim] * changeCellLocal[jDim];
+ orientation[jDim*spaceDim+iDim] * changeCellLocal[iSpace+jDim];
} // for
} // if
@@ -684,15 +693,17 @@
const int spaceDim = _quadrature->spaceDim();
const int numBasis = _quadrature->numBasis();
const int numQuadPts = _quadrature->numQuadPts();
- const int fiberDim = spaceDim * 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) {
@@ -706,41 +717,41 @@
if (0 != _dbRate) {
assert(!rateSection.isNull());
assert(!rateTimeSection.isNull());
- double tRate = 0.0;
rateSection->restrictPoint(*c_iter, &bufferCell[0], bufferCell.size());
- rateTimeSection->restrictPoint(*c_iter, &tRate, 1);
- if (t > tRate) { // rate of change integrated over time
- bufferCell *= (t - tRate);
- valuesCell += bufferCell;
- } // if
+ 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());
- double tChange = 0.0;
changeSection->restrictPoint(*c_iter, &bufferCell[0], bufferCell.size());
- changeTimeSection->restrictPoint(*c_iter, &tChange, 1);
- if (t >= tChange) { // change in value over time
- double scale = 1.0;
- if (0 != _dbTimeHistory) {
- double tDim = t - tChange;
- _getNormalizer().dimensionalize(&tDim, 1, timeScale);
- const int err = _dbTimeHistory->query(&scale, tDim);
- if (0 != err) {
- std::ostringstream msg;
- msg << "Error querying for time '" << tDim
- << "' in time history database "
- << _dbTimeHistory->label() << ".";
- throw std::runtime_error(msg.str());
+ 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
- bufferCell *= scale;
- valuesCell += bufferCell;
- } // if
} // if
valueSection->updateAddPoint(*c_iter, &valuesCell[0]);
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.hh 2009-06-18 18:56:56 UTC (rev 15332)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.hh 2009-06-18 20:54:31 UTC (rev 15333)
@@ -41,12 +41,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.
@@ -96,11 +90,6 @@
cellField(const char* name,
topology::SolutionFields* const fields);
- // PRIVATE MEMBERS ////////////////////////////////////////////////////
-private :
-
- spatialdata::spatialdb::SpatialDB* _db; ///< Spatial database w/parameters
-
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -110,6 +99,12 @@
*/
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);
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.icc 2009-06-18 18:56:56 UTC (rev 15332)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann_NEW.icc 2009-06-18 20:54:31 UTC (rev 15333)
@@ -14,12 +14,7 @@
#error "Neumann_NEW.icc can only be included from Neumann_NEW.hh"
#endif
-// Set database for boundary condition parameters.
-inline
-void
-pylith::bc::Neumann_NEW::db(spatialdata::spatialdb::SpatialDB* const db) {
- _db = db;
-}
+#include <cassert> // USES assert()
// Get label of boundary condition surface.
inline
@@ -28,5 +23,13 @@
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
Modified: short/3D/PyLith/trunk/libsrc/topology/SubMesh.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/SubMesh.hh 2009-06-18 18:56:56 UTC (rev 15332)
+++ short/3D/PyLith/trunk/libsrc/topology/SubMesh.hh 2009-06-18 20:54:31 UTC (rev 15333)
@@ -130,7 +130,7 @@
*
* @param label Label for mesh.
*/
- void view(const char* label);
+ void view(const char* label) const;
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/topology/SubMesh.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/SubMesh.icc 2009-06-18 18:56:56 UTC (rev 15332)
+++ short/3D/PyLith/trunk/libsrc/topology/SubMesh.icc 2009-06-18 20:54:31 UTC (rev 15333)
@@ -68,7 +68,7 @@
// Print mesh to stdout.
inline
void
-pylith::topology::SubMesh::view(const char* label) {
+pylith::topology::SubMesh::view(const char* label) const {
_mesh->view(label);
}
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am 2009-06-18 18:56:56 UTC (rev 15332)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am 2009-06-18 20:54:31 UTC (rev 15333)
@@ -52,6 +52,7 @@
TestDirichletBoundaryTet4.cc \
TestDirichletBoundaryHex8.cc \
TestNeumann.cc \
+ TestNeumann_NEW.cc \
TestNeumannLine2.cc \
TestNeumannTri3.cc \
TestNeumannQuad4.cc \
@@ -98,6 +99,7 @@
TestDirichletBoundaryTet4.hh \
TestDirichletBoundaryHex8.hh \
TestNeumann.hh \
+ TestNeumann_NEW.hh \
TestNeumannLine2.hh \
TestNeumannTri3.hh \
TestNeumannQuad4.hh \
Added: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.cc 2009-06-18 20:54:31 UTC (rev 15333)
@@ -0,0 +1,839 @@
+// -*- 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
Added: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann_NEW.hh 2009-06-18 20:54:31 UTC (rev 15333)
@@ -0,0 +1,141 @@
+// -*- 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
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/Makefile.am 2009-06-18 18:56:56 UTC (rev 15332)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/Makefile.am 2009-06-18 20:54:31 UTC (rev 15333)
@@ -22,6 +22,7 @@
tri3_disp2.spatialdb \
tri3_vel2.spatialdb \
tri3_tractions.spatialdb \
+ tri3_tractions_NEW.spatialdb \
tri3_force.spatialdb \
tri3_force_rate.spatialdb \
tri3_force_change.spatialdb \
@@ -30,6 +31,10 @@
quad4.mesh \
quad4_disp.spatialdb \
quad4_tractions.spatialdb \
+ quad4_traction_initial.spatialdb \
+ quad4_traction_rate.spatialdb \
+ quad4_traction_change.spatialdb \
+ quad4_traction.timedb \
quad4_force.spatialdb \
tet4.mesh \
tet4_disp.spatialdb \
Added: short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction.timedb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction.timedb (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction.timedb 2009-06-18 20:54:31 UTC (rev 15333)
@@ -0,0 +1,7 @@
+#TIME HISTORY ascii
+TimeHistory {
+ num-points = 2
+ time-units = second
+}
+ 0.0 1.0
+ 10.0 0.0
Added: short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction_change.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction_change.spatialdb (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction_change.spatialdb 2009-06-18 20:54:31 UTC (rev 15333)
@@ -0,0 +1,17 @@
+#SPATIAL.ascii 1
+SimpleDB {
+ num-values = 3
+ value-names = traction-shear traction-normal change-start-time
+ value-units = Pa Pa s
+ num-locs = 4
+ data-dim = 1
+ space-dim = 2
+ cs-data = cartesian {
+ to-meters = 1.0
+ space-dim = 2
+ }
+}
+-0.667 0.0 1.3 1.4 2.0
+-0.333 0.0 1.7 1.6 2.4
++0.333 0.0 2.3 2.4 2.1
++0.667 0.0 2.7 2.6 2.5
Added: short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction_initial.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction_initial.spatialdb (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction_initial.spatialdb 2009-06-18 20:54:31 UTC (rev 15333)
@@ -0,0 +1,17 @@
+#SPATIAL.ascii 1
+SimpleDB {
+ num-values = 2
+ value-names = traction-shear traction-normal
+ value-units = Pa Pa
+ num-locs = 4
+ data-dim = 1
+ space-dim = 2
+ cs-data = cartesian {
+ to-meters = 1.0
+ space-dim = 2
+ }
+}
+-0.667 0.0 0.3 0.4
+-0.333 0.0 0.7 0.6
++0.333 0.0 1.3 1.4
++0.667 0.0 1.7 1.6
Added: short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction_rate.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction_rate.spatialdb (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4_traction_rate.spatialdb 2009-06-18 20:54:31 UTC (rev 15333)
@@ -0,0 +1,17 @@
+#SPATIAL.ascii 1
+SimpleDB {
+ num-values = 3
+ value-names = traction-rate-shear traction-rate-normal rate-start-time
+ value-units = Pa Pa s
+ num-locs = 4
+ data-dim = 1
+ space-dim = 2
+ cs-data = cartesian {
+ to-meters = 1.0
+ space-dim = 2
+ }
+}
+-0.667 0.0 -0.2 -0.1 0.5
+-0.333 0.0 0.4 0.3 0.8
++0.333 0.0 -1.2 -1.1 0.6
++0.667 0.0 1.4 1.3 0.9
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2009-06-18 18:56:56 UTC (rev 15332)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2009-06-18 20:54:31 UTC (rev 15333)
@@ -167,7 +167,6 @@
const ALE::Obj<RealSection>& orientationSection =
fault._fields->get("orientation").section();
CPPUNIT_ASSERT(!orientationSection.isNull());
- const int cellDim = _data->cellDim;
const int spaceDim = _data->spaceDim;
const int orientationSize = spaceDim*spaceDim;
iVertex = 0;
@@ -246,7 +245,7 @@
fault.useSolnIncr(false);
fault.integrateResidual(residual, t, &fields);
- residual.view("RESIDUAL"); // DEBUGGING
+ //residual.view("RESIDUAL"); // DEBUGGING
// Check values
const double* valsE = _data->valsResidual;
More information about the CIG-COMMITS
mailing list