[cig-commits] r15125 - in short/3D/PyLith/trunk: libsrc libsrc/bc libsrc/topology unittests/libtests/topology
brad at geodynamics.org
brad at geodynamics.org
Fri Jun 5 17:19:51 PDT 2009
Author: brad
Date: 2009-06-05 17:19:51 -0700 (Fri, 05 Jun 2009)
New Revision: 15125
Modified:
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/bc/DirichletBC.cc
short/3D/PyLith/trunk/libsrc/bc/PointForce.cc
short/3D/PyLith/trunk/libsrc/bc/PointForce.hh
short/3D/PyLith/trunk/libsrc/bc/PointForce.icc
short/3D/PyLith/trunk/libsrc/topology/Field.cc
short/3D/PyLith/trunk/libsrc/topology/Field.hh
short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.cc
short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.hh
Log:
Worked on setting up time dependent point force BC.
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2009-06-05 19:26:11 UTC (rev 15124)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2009-06-06 00:19:51 UTC (rev 15125)
@@ -28,6 +28,7 @@
bc/DirichletBoundary.cc \
bc/Neumann.cc \
bc/AbsorbingDampers.cc \
+ bc/PointForce.cc \
faults/Fault.cc \
faults/SlipTimeFn.cc \
faults/StepSlipFn.cc \
Modified: short/3D/PyLith/trunk/libsrc/bc/DirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/DirichletBC.cc 2009-06-05 19:26:11 UTC (rev 15124)
+++ short/3D/PyLith/trunk/libsrc/bc/DirichletBC.cc 2009-06-06 00:19:51 UTC (rev 15125)
@@ -25,6 +25,8 @@
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
+//#define FIELD_SPLIT
+
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
@@ -106,7 +108,7 @@
} // if
_offsetLocal[iPoint] = curNumConstraints;
section->addConstraintDimension(_points[iPoint], numFixedDOF);
-#if 0 // TODO: FIELD SPLIT
+#if defined(FIELD_SPLIT)
if (fibration >= 0) {
assert(fiberDim == section->getFiberDimension(_points[iPoint],
fibration));
@@ -178,7 +180,7 @@
// Update list of constrained DOF
section->setConstraintDof(point, &allFixedDOF[0]);
-#if 0 // TODO: FIELD SPLIT
+#if defined(FIELD_SPLIT)
if (fibration >= 0)
section->setConstraintDof(point, &allFixedDOF[0], fibration);
#endif
Modified: short/3D/PyLith/trunk/libsrc/bc/PointForce.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/PointForce.cc 2009-06-05 19:26:11 UTC (rev 15124)
+++ short/3D/PyLith/trunk/libsrc/bc/PointForce.cc 2009-06-06 00:19:51 UTC (rev 15125)
@@ -15,8 +15,8 @@
#include "PointForce.hh" // implementation of object methods
#include "pylith/topology/Field.hh" // USES Field
-#include "pylith/topology/Mesh.hh" // USES Mesh
#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
@@ -31,7 +31,14 @@
// ----------------------------------------------------------------------
// Default constructor.
-pylith::bc::PointForce::PointForce(void)
+pylith::bc::PointForce::PointForce(void) :
+ _parameters(0),
+ _dbInitial(0),
+ _dbRate(0),
+ _dbRateTime(0),
+ _dbChange(0),
+ _dbChangeTime(0),
+ _dbTimeHistory(0)
{ // constructor
} // constructor
@@ -39,6 +46,14 @@
// Destructor.
pylith::bc::PointForce::~PointForce(void)
{ // destructor
+ delete _parameters; _parameters = 0;
+
+ _dbInitial = 0; // TODO: Use shared pointers
+ _dbRate = 0; // TODO: Use shared pointers
+ _dbRateTime = 0; // TODO: Use shared pointers
+ _dbChange = 0; // TODO: Use shared pointers
+ _dbChangeTime = 0; // TODO: Use shared pointers
+ _dbTimeHistory = 0; // TODO: Use shared pointers
} // destructor
// ----------------------------------------------------------------------
@@ -50,9 +65,9 @@
if (size > 0)
assert(0 != flags);
- _forceDOF.resize(size);
+ _bcDOF.resize(size);
for (int i=0; i < size; ++i)
- _forceDOF[i] = flags[i];
+ _bcDOF[i] = flags[i];
} // forceDOF
// ----------------------------------------------------------------------
@@ -61,67 +76,200 @@
pylith::bc::PointForce::initialize(const topology::Mesh& mesh,
const double upDir[3])
{ // initialize
- const int numForceDOF = _forceDOF.size();
- if (0 == numForceDOF)
+ assert(0 != _normalizer);
+
+ if (0 == _bcDOF.size())
return;
_getPoints(mesh);
- _setupQueryDatabases();
- _queryDatabases(mesh);
+
+ const double lengthScale = _normalizer->lengthScale();
+ const double pressureScale = _normalizer->pressureScale();
+ const double forceScale = pressureScale * lengthScale * lengthScale;
+
+ _queryDatabases(mesh, forceScale, "force");
} // initialize
// ----------------------------------------------------------------------
// Integrate contributions to residual term (r) for operator that
// do not require assembly over cells, vertices, or processors.
void
-integrateResidualAssembled(const topology::Field<topology::Mesh>& residual,
+pylith::bc::PointForce::integrateResidualAssembled(
+ topology::Field<topology::Mesh>* residual,
const double t,
topology::SolutionFields* const fields)
{ // integrateResidualAssembled
+ assert(0 != residual);
+ assert(0 != _parameters);
+ assert(0 != _normalizer);
+
+ // Calculate spatial and temporal variation of value for BC.
+ _calculateValue(t);
+
const int numPoints = _points.size();
- const int numForceDOF = _forceDOF.size();
+ const int numBCDOF = _bcDOF.size();
- const ALE::Obj<RealSection>& amplitudeSection =
- _parameters->get("amplitude").section();
- assert(!amplitudeSection.isNull());
+ const topology::Mesh& mesh = residual->mesh();
+ const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
+ assert(0 != cs);
+ const int spaceDim = cs->spaceDim();
- const ALE::Obj<RealSection>& residualSection = residual.section();
+ double_array residualVertex(spaceDim);
+ const ALE::Obj<RealSection>& residualSection = residual->section();
assert(!residualSection.isNull());
- const ALE::Obj<RealSection>& startTimeSection = (0 == _dbStartTime) ?
- 0 : _parameters->get("start time").section();
-
- double_array forcesVertex(spaceDim);
- double scale = 1.0;
+ double_array valuesVertex(numBCDOF);
+ const ALE::Obj<RealSection>& valueSection =
+ _parameters->get("value").section();
+ assert(!valueSection.isNull());
+
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
- const int p_force = _points[iPoint];
- assert(numForceDOF == amplitudeSection->getFiberDimension(p_force));
- const double* amplitudeVertex = amplitudeSection->restrictPoint(p_force);
- assert(0 != amplitudeVertex);
+ const int p_bc = _points[iPoint]; // Get point label.
+ residualVertex *= 0.0; // Reset residual contribution to zero.
+
+ valueSection->restrictPoint(p_bc, &valuesVertex[0], valuesVertex.size());
+ for (int iDOF=0; iDOF < numBCDOF; ++iDOF)
+ residualVertex[_bcDOF[iDOF]] += valuesVertex[iDOF];
+ residualSection->updateAddPoint(p_bc, &residualVertex[0]);
+ } // for
- if (0 != _dbStartTime) {
- assert(!startTimeSection.isNull()); // Expect section with start times
- assert(0 != _dbTimeAmp); // Expect database with temporal evolution
+} // integrateResidualAssembled
- assert(1 == startTimeSection->getFiberDimension(p_force));
- const double tRef = *startTimeSection->restrictPoint(p_force);
- const double tRel = t - tRef;
- if (tRel > 0.0) {
- _normalizer->dimensionalize(&tRel, 1, timeScale);
- scale = _dbTimeAmp->query(tRel);
+// ----------------------------------------------------------------------
+// Calculate temporal and spatial variation of value over the list of points.
+void
+pylith::bc::PointForce::_calculateValue(const double t)
+{ // _calculateValue
+ assert(0 != _parameters);
+ assert(0 != _normalizer);
+
+ const ALE::Obj<RealSection>& valueSection =
+ _parameters->get("value").section();
+ valueSection->zero();
+
+ const int numPoints = _points.size();
+ const int numBCDOF = _bcDOF.size();
+ const double timeScale = _normalizer->timeScale();
+
+ 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;
+
+ // Contribution from initial value
+ if (0 != _dbInitial) {
+ double_array values(numBCDOF);
+ const ALE::Obj<RealSection>& initialSection =
+ _parameters->get("initial").section();
+ assert(!initialSection.isNull());
+
+ for (int iPoint=0; iPoint < numPoints; ++iPoint) {
+ const int p_bc = _points[iPoint]; // Get point label
+
+ initialSection->restrictPoint(p_bc, &values[0], values.size());
+ valueSection->updateAddPoint(p_bc, &values[0]);
+ } // for
+ } // if
+
+
+ // Contribution from rate of change of value
+ if (0 != _dbRate) {
+ double_array values(numBCDOF);
+ double tRate = 0.0;
+ const ALE::Obj<RealSection>& rateSection =
+ _parameters->get("rate").section();
+ assert(!rateSection.isNull());
+ const ALE::Obj<RealSection>& rateTimeSection =
+ _parameters->get("rate time").section();
+ assert(!rateTimeSection.isNull());
+
+ for (int iPoint=0; iPoint < numPoints; ++iPoint) {
+ const int p_bc = _points[iPoint]; // Get point label
+
+ rateSection->restrictPoint(p_bc, &values[0], values.size());
+ rateTimeSection->restrictPoint(p_bc, &tRate, 1);
+ if (t > tRate) { // rate of change integrated over time
+ values *= (t - tRate);
+ valueSection->updateAddPoint(p_bc, &values[0]);
} // if
- } // if
-
- for (int iDOF=0; iDOF < numForceDOF; ++iDOF)
- forcesVertex[_forceDOF[iDOF]] = amplitudeVertex[iDOF]*scale;
- } // for
-} // interateResidualAssembled
+ } // for
+ } // if
+
+ // Contribution from change of value
+ if (0 != _dbChange) {
+ double_array values(numBCDOF);
+ double tChange = 0.0;
+ const ALE::Obj<RealSection>& changeSection =
+ _parameters->get("change").section();
+ assert(!changeSection.isNull());
+ const ALE::Obj<RealSection>& changeTimeSection =
+ _parameters->get("change time").section();
+ assert(!changeTimeSection.isNull());
+ for (int iPoint=0; iPoint < numPoints; ++iPoint) {
+ const int p_bc = _points[iPoint]; // Get point label
+
+ changeSection->restrictPoint(p_bc, &values[0], values.size());
+ changeTimeSection->restrictPoint(p_bc, &tChange, 1);
+ if (t > tChange) { // change in value over time
+ double scale = 1.0;
+ if (0 != _dbTimeHistory) {
+ double tDim = t - tChange;
+ _normalizer->dimensionalize(&tDim, 1, timeScale);
+ const int err = _dbTimeHistory->query(&scale, tDim);
+ if (0 != err) {
+ std::ostringstream msg;
+ msg << "ADD SOMETHING HERE";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // if
+ values *= scale;
+ valueSection->updateAddPoint(p_bc, &values[0]);
+ } // if
+ } // for
+ } // if
+} // _calculateValue
+
// ----------------------------------------------------------------------
// Verify configuration is acceptable.
void
-verifyConfiguration(const topology::Mesh& mesh) const
+pylith::bc::PointForce::verifyConfiguration(const topology::Mesh& mesh) const
{ // verifyConfiguration
+ if (0 != _dbRate && 0 == _dbRateTime) {
+ std::ostringstream msg;
+ msg << "Point force boundary condition '" << label() << "',\n has a rate "
+ << "of change spatial database but no rate of change start time "
+ << "spatial database.";
+ throw std::runtime_error(msg.str());
+ } // if
+ if (0 == _dbRate && 0 != _dbRateTime) {
+ std::ostringstream msg;
+ msg << "Point force boundary condition '" << label() << "',\n has a rate "
+ << "of change start time spatial database but no rate of change "
+ << "spatial database.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ if (0 != _dbChange && 0 == _dbChangeTime) {
+ std::ostringstream msg;
+ msg << "Point force boundary condition '" << label() << "',\n has a "
+ << "change in value spatial database but change in value start time "
+ << "spatial database.";
+ throw std::runtime_error(msg.str());
+ } // if
+ if (0 == _dbChange && 0 != _dbChangeTime) {
+ std::ostringstream msg;
+ msg << "Point force boundary condition '" << label() << "',\n has a "
+ << "change in value start time spatial database but change in value "
+ << "spatial database.";
+ throw std::runtime_error(msg.str());
+ } // if
+ if (0 == _dbChange && 0 != _dbTimeHistory) {
+ std::ostringstream msg;
+ msg << "Point force boundary condition '" << label() << "',\n has a "
+ << "time history database but not change in value spatial database.";
+ throw std::runtime_error(msg.str());
+} // if
} // verifyConfiguration
// ----------------------------------------------------------------------
@@ -155,90 +303,104 @@
} // _getPoints
// ----------------------------------------------------------------------
-// Setup initial and rate of change databases for querying.
+// Query databases for values.
void
-pylith::bc::PointForce::_setupQueryDatabases(void)
-{ // _setupQueryDatabases
- assert(0 != _db);
+pylith::bc::PointForce::_queryDatabases(const topology::Mesh& mesh,
+ const double valueScale,
+ const char* fieldName)
+{ // _queryDatabases
+ assert(0 != _normalizer);
+ const double timeScale = _normalizer->timeScale();
+ const double rateScale = valueScale / timeScale;
- const int numForceDOF = _forceDOF.size();
- char** valueNames = (numForceDOF > 0) ? new char*[numForceDOF] : 0;
- for (int i=0; i < numForceDOF; ++i) {
+ const int numPoints = _points.size();
+ const int numBCDOF = _bcDOF.size();
+ char** valueNames = (numBCDOF > 0) ? new char*[numBCDOF] : 0;
+ for (int i=0; i < numBCDOF; ++i) {
std::ostringstream name;
- name << "dof-" << _forceDOF[i];
+ name << "dof-" << _bcDOF[i];
const int size = 1 + name.str().length();
valueNames[i] = new char[size];
strcpy(valueNames[i], name.str().c_str());
} // for
- // Setup initial database.
- _db->open();
- _db->queryVals(const_cast<const char**>(valueNames), numForceDOF);
+ delete _parameters;
+ _parameters = new topology::Fields<topology::Field<topology::Mesh> >(mesh);
- // Setup rate database, if provided.
- if (0 != _dbRate) {
- _dbRate->open();
- _dbRate->queryVals((const char**) valueNames, numForceDOF);
+ if (0 != _dbInitial) { // Setup initial values, if provided.
+ std::string fieldLabel = std::string("initial_") + std::string(fieldName);
+ _parameters->add("initial", fieldLabel.c_str());
+ topology::Field<topology::Mesh>& initial =
+ _parameters->get("initial");
+ initial.newSection(_points, numBCDOF);
+ initial.allocate();
+ initial.scale(valueScale);
+ initial.vectorFieldType(topology::FieldBase::VECTOR);
+
+ _dbInitial->queryVals(valueNames, numBCDOF);
+ _queryDB(&initial, _dbInitial, numBCDOF, valueScale);
} // if
- for (int i=0; i < numForceDOF; ++i) {
- delete[] valueNames[i]; valueNames[i] = 0;
- } // for
- delete[] valueNames; valueNames = 0;
-} // _setupQueryDatabases
+} // _queryDatabases
+
+
// ----------------------------------------------------------------------
-// Query initial and rate of change databases for values.
+// Query database for values.
void
-pylith::bc::PointForce::_queryDatabases(const topology::Mesh& mesh)
-{ // _queryDatabases
- assert(0 != _db);
+pylith::bc::PointForce::_queryDB(topology::Field<topology::Mesh>* field,
+ spatialdata::spatialdb::SpatialDB* const db,
+ const int querySize,
+ const double scale)
+{ // _queryDB
+ assert(0 != field);
+ assert(0 != db);
+ assert(0 != _normalizer);
+ const topology::Mesh& mesh = field->mesh();
+ const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
+ assert(0 != cs);
+ const int spaceDim = cs->spaceDim();
+
+ assert(0 != _normalizer);
+ const double lengthScale = _normalizer->lengthScale();
+
+ double_array coordsVertex(spaceDim);
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
assert(!sieveMesh.isNull());
-
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
- assert(0 != cs);
- const int spaceDim = cs->spaceDim();
+ const ALE::Obj<RealSection>& section = field->section();
+ assert(!section.isNull());
- assert(0 != _normalizer);
- const double lengthScale = _normalizer->lengthScale();
- const double velocityScale =
- _normalizer->lengthScale() / _normalizer->timeScale();
+ double_array valuesVertex(querySize);
+ db->open();
const int numPoints = _points.size();
- const int numForceDOF = _forceDOF.size();
- _valuesInitial.resize(numPoints*numForceDOF);
- if (0 != _dbRate)
- _valuesRate.resize(numPoints*numForceDOF);
-
- double_array queryValues(numForceDOF);
- double_array vCoordsGlobal(spaceDim);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
- // Get coordinates of vertex
+ // Get dimensionalized coordinates of vertex
coordinates->restrictPoint(_points[iPoint],
- &vCoordsGlobal[0], vCoordsGlobal.size());
- _normalizer->dimensionalize(&vCoordsGlobal[0], vCoordsGlobal.size(),
+ &coordsVertex[0], coordsVertex.size());
+ _normalizer->dimensionalize(&coordsVertex[0], coordsVertex.size(),
lengthScale);
- int err = _db->query(&queryValues[0], numForceDOF,
- &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
+ int err = db->query(&valuesVertex[0], valuesVertex.size(),
+ &coordsVertex[0], coordsVertex.size(), cs);
if (err) {
std::ostringstream msg;
- msg << "Could not find initial values at (";
+ msg << "Error querying for '" << field->label() << "' at (";
for (int i=0; i < spaceDim; ++i)
- msg << " " << vCoordsGlobal[i];
- msg << ") using spatial database " << _db->label() << ".";
+ msg << " " << coordsVertex[i];
+ msg << ") using spatial database " << db->label() << ".";
throw std::runtime_error(msg.str());
} // if
- for (int iDOF=0; iDOF < numForceDOF; ++iDOF)
- _valuesInitial[numForceDOF*iPoint+iDOF] =
- _normalizer->nondimensionalize(queryValues[iDOF], lengthScale);
+ _normalizer->nondimensionalize(&valuesVertex[0], valuesVertex.size(),
+ scale);
+ section->updatePoint(_points[iPoint], &valuesVertex[0]);
} // for
- _db->close();
-} // _queryDatabases
+ db->close();
+} // _queryDB
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/bc/PointForce.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/PointForce.hh 2009-06-05 19:26:11 UTC (rev 15124)
+++ short/3D/PyLith/trunk/libsrc/bc/PointForce.hh 2009-06-06 00:19:51 UTC (rev 15125)
@@ -20,6 +20,8 @@
// Include directives ---------------------------------------------------
#include "BoundaryCondition.hh" // ISA BoundaryCondition
+#include "pylith/topology/Mesh.hh" // ISA Integrator<Quadrature<Mesh> >
+#include "pylith/feassemble/Quadrature.hh" // ISA Integrator<Quadrature<Mesh >
#include "pylith/feassemble/Integrator.hh" // ISA Integrator
#include "pylith/topology/topologyfwd.hh" // USES Fields<Mesh>
@@ -27,7 +29,7 @@
// PointForce ------------------------------------------------------
class pylith::bc::PointForce : public BoundaryCondition,
- public feassemble::Integrator
+ public feassemble::Integrator<feassemble::Quadrature<topology::Mesh> >
{ // class PointForce
friend class TestPointForce; // unit testing
@@ -53,6 +55,43 @@
void forceDOF(const int* flags,
const int size);
+ /** Set database for initial values.
+ *
+ * @param db Spatial database
+ */
+ void dbInitial(spatialdata::spatialdb::SpatialDB* const db);
+
+ /** Set database for rate of change of values.
+ *
+ * @param db Spatial database
+ */
+ void dbRate(spatialdata::spatialdb::SpatialDB* const db);
+
+ /** Set database for start time of rate change.
+ *
+ * @param db Spatial database
+ */
+ void dbRateTime(spatialdata::spatialdb::SpatialDB* const db);
+
+ /** Set database for change in values.
+ *
+ * @param db Spatial database
+ */
+ void dbChange(spatialdata::spatialdb::SpatialDB* const db);
+
+ /** Set database for start time of change in values.
+ *
+ * @param db Spatial database
+ */
+ void dbChangeTime(spatialdata::spatialdb::SpatialDB* const db);
+
+ /** Set database for temporal evolution of change in value.
+ *
+ * @param db Time history database.
+ */
+ void dbTimeHistory(spatialdata::spatialdb::TimeHistory* const db);
+
+
/** Initialize boundary condition.
*
* @param mesh PETSc mesh
@@ -67,9 +106,9 @@
* @param t Current time
* @param fields Solution fields
*/
- void integrateResidual(const topology::Field<topology::Mesh>& residual,
- const double t,
- topology::SolutionFields* const fields);
+ void integrateResidualAssembled(topology::Field<topology::Mesh>* residual,
+ const double t,
+ topology::SolutionFields* const fields);
/** Verify configuration is acceptable.
*
@@ -86,40 +125,68 @@
*/
void _getPoints(const topology::Mesh& mesh);
- /// Setup databases for querying for point forces.
- void _setupQueryDatabases(void);
-
- /** Query databases for point forces.
+ /** Query databases for time dependent parameters.
*
* @param mesh Finite-element mesh.
+ * @param valueScale Dimension scale for value.
+ * @param fieldName Name of field associated with value.
*/
- void _queryDatabases(const topology::Mesh& mesh);
+ void _queryDatabases(const topology::Mesh& mesh,
+ const double valueScale,
+ const char* fieldName);
- // NOT IMPLEMENTED ////////////////////////////////////////////////////
-private :
+ /** Wuery 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::Mesh>* field,
+ spatialdata::spatialdb::SpatialDB* const db,
+ const int querySize,
+ const double scale);
- /// Not implemented
- PointForce(const PointForce& m);
+ /** Calculate spatial and temporal variation of value over the list
+ * of points.
+ *
+ * @param t Current time.
+ */
+ void _calculateValue(const double t);
- /// Not implemented
- const PointForce& operator=(const PointForce& m);
-
// PROTECTED MEMBERS //////////////////////////////////////////////////
protected :
/// Parameters for point forces.
topology::Fields<topology::Field<topology::Mesh> >* _parameters;
- /// Start time for point forces.
- spatialdata::spatialdb::SpatialDB* _dbStartTime;
+ /// Spatial database for initial values.
+ spatialdata::spatialdb::SpatialDB* _dbInitial;
- /// Temporal evolution for amplitude of point forces.
- spatialdata::spatialdb::SpatialDB* _dbTimeAmp;
+ /// Spatial database for rate of change of values.
+ spatialdata::spatialdb::SpatialDB* _dbRate;
+
+ /// Spatial database for start time of rate change.
+ spatialdata::spatialdb::SpatialDB* _dbRateTime;
+
+ /// Spatial database for change in value.
+ spatialdata::spatialdb::SpatialDB* _dbChange;
+
+ /// Spatial database for start time of change in value.
+ spatialdata::spatialdb::SpatialDB* _dbChangeTime;
+
+ /// Temporal evolution of amplitude for change in value;
+ spatialdata::spatialdb::TimeHistory* _dbTimeHistory;
-
int_array _points; ///< Points for forces.
- int_array _forceDOF; ///< Indices of degrees of freedom with forces.
+ int_array _bcDOF; ///< Indices of degrees of freedom with boundary condition.
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ PointForce(const PointForce&); ///< Not implemented.
+ const PointForce& operator=(const PointForce&); ///< Not implemented.
+
}; // class PointForce
#include "PointForce.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/bc/PointForce.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/PointForce.icc 2009-06-05 19:26:11 UTC (rev 15124)
+++ short/3D/PyLith/trunk/libsrc/bc/PointForce.icc 2009-06-06 00:19:51 UTC (rev 15125)
@@ -14,5 +14,47 @@
#error "PointForce.icc can only be included from PointForce.hh"
#endif
+// Set database for initial values.
+inline
+void
+pylith::bc::PointForce::dbInitial(spatialdata::spatialdb::SpatialDB* const db) {
+ _dbInitial = db;
+}
+// Set database for rate of change of values.
+inline
+void
+pylith::bc::PointForce::dbRate(spatialdata::spatialdb::SpatialDB* const db) {
+ _dbRate = db;
+}
+
+// Set database for start time of rate change.
+inline
+void
+pylith::bc::PointForce::dbRateTime(spatialdata::spatialdb::SpatialDB* const db) {
+ _dbRateTime = db;
+}
+
+// Set database for change in values.
+inline
+void
+pylith::bc::PointForce::dbChange(spatialdata::spatialdb::SpatialDB* const db) {
+ _dbChange = db;
+}
+
+// Set database for start time of change in values.
+inline
+void
+pylith::bc::PointForce::dbChangeTime(spatialdata::spatialdb::SpatialDB* const db) {
+ _dbChangeTime = db;
+}
+
+// Set database for temporal evolution of change in value.
+inline
+void
+pylith::bc::PointForce::dbTimeHistory(spatialdata::spatialdb::TimeHistory* const db) {
+ _dbTimeHistory = db;
+}
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Field.cc 2009-06-05 19:26:11 UTC (rev 15124)
+++ short/3D/PyLith/trunk/libsrc/topology/Field.cc 2009-06-06 00:19:51 UTC (rev 15125)
@@ -106,7 +106,8 @@
} // newSection
// ----------------------------------------------------------------------
-// Create sieve section and set chart and fiber dimesion.
+// Create sieve section and set chart and fiber dimesion for a
+// sequence of points.
template<typename mesh_type>
void
pylith::topology::Field<mesh_type>::newSection(
@@ -120,9 +121,8 @@
logger.stagePush("Field");
if (fiberDim < 0) {
std::ostringstream msg;
- msg
- << "Fiber dimension (" << fiberDim << ") for field '" << _label
- << "' must be nonnegative.";
+ msg << "Fiber dimension (" << fiberDim << ") for field '" << _label
+ << "' must be nonnegative.";
throw std::runtime_error(msg.str());
} // if
@@ -135,15 +135,48 @@
*std::max_element(points->begin(), points->end());
_section->setChart(chart_type(pointMin, pointMax+1));
_section->setFiberDimension(points, fiberDim);
- } else {
- // Create empty chart
+ } else // Create empty chart
_section->setChart(chart_type(0, 0));
- _section->setFiberDimension(points, fiberDim);
- } // if/else
+
logger.stagePop();
} // newSection
// ----------------------------------------------------------------------
+// Create sieve section and set chart and fiber dimesion for a list of
+// points.
+template<typename mesh_type>
+void
+pylith::topology::Field<mesh_type>::newSection(const int_array& points,
+ const int fiberDim)
+{ // newSection
+ typedef typename mesh_type::SieveMesh::point_type point_type;
+
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ std::cout << "Making Field " << _label << " section type 1b" << std::endl;
+ logger.stagePush("Field");
+ if (fiberDim < 0) {
+ std::ostringstream msg;
+ msg << "Fiber dimension (" << fiberDim << ") for field '" << _label
+ << "' must be nonnegative.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ _section = new RealSection(_mesh.comm(), _mesh.debug());
+
+ const int npts = points.size();
+ if (npts > 0) {
+ const point_type pointMin = points.min();
+ const point_type pointMax = points.max();
+ _section->setChart(chart_type(pointMin, pointMax+1));
+ for (int i=0; i < npts; ++i)
+ _section->setFiberDimension(points[i], fiberDim);
+ } else // create empty chart
+ _section->setChart(chart_type(0, 0));
+
+ logger.stagePop();
+} // newSection
+
+// ----------------------------------------------------------------------
// Create sieve section and set chart and fiber dimesion.
template<typename mesh_type>
void
Modified: short/3D/PyLith/trunk/libsrc/topology/Field.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Field.hh 2009-06-05 19:26:11 UTC (rev 15124)
+++ short/3D/PyLith/trunk/libsrc/topology/Field.hh 2009-06-06 00:19:51 UTC (rev 15125)
@@ -25,6 +25,7 @@
// Include directives ---------------------------------------------------
#include "FieldBase.hh" // ISA FieldBase
+#include "pylith/utils/arrayfwd.hh" // HASA int_array
#include "pylith/utils/petscfwd.h" // HASA PetscVec
#include <petscmesh.hh>
@@ -147,7 +148,8 @@
/// Create sieve section.
void newSection(void);
- /** Create sieve section and set chart and fiber dimesion.
+ /** Create sieve section and set chart and fiber dimesion for
+ * sequence of points.
*
* @param points Points over which to define section.
* @param dim Fiber dimension for section.
@@ -155,6 +157,15 @@
void newSection(const ALE::Obj<label_sequence>& points,
const int fiberDim);
+ /** Create sieve section and set chart and fiber dimesion for a list
+ * of points.
+ *
+ * @param points Points over which to define section.
+ * @param dim Fiber dimension for section.
+ */
+ void newSection(const int_array& points,
+ const int fiberDim);
+
/** Create sieve section and set chart and fiber dimesion.
*
* @param domain Type of points over which to define section.
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.cc 2009-06-05 19:26:11 UTC (rev 15124)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.cc 2009-06-06 00:19:51 UTC (rev 15125)
@@ -21,6 +21,8 @@
#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
+#define FIELD_SPLIT
+
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::topology::TestFieldMesh );
@@ -168,6 +170,54 @@
} // testNewSectionPoints
// ----------------------------------------------------------------------
+// Test newSection(int_array).
+void
+pylith::topology::TestFieldMesh::testNewSectionPointsArray(void)
+{ // testNewSectionPointsArray
+ const int fiberDim = 2;
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+
+ Field<Mesh> field(mesh);
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const int npts = vertices->size() / 2;
+ int_array pointsIn(npts);
+ int_array pointsOut(vertices->size() - npts);
+ int count = 0;
+ size_t iIn = 0;
+ size_t iOut = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ if (count % 2 == 0)
+ pointsIn[iIn++] = *v_iter;
+ else
+ pointsOut[iOut++] = *v_iter;
+ ++count;
+ } // for
+ CPPUNIT_ASSERT_EQUAL(iIn, pointsIn.size());
+ CPPUNIT_ASSERT_EQUAL(iOut, pointsOut.size());
+
+ field.newSection(pointsIn, fiberDim);
+ field.allocate();
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ // Points in array should have a fiber dimension of fiberDim.
+ for (int i=0; i < pointsIn.size(); ++i)
+ CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(pointsIn[i]));
+
+ // Points not int array should have a fiber dimension of zero.
+ for (int i=0; i < pointsOut.size(); ++i)
+ CPPUNIT_ASSERT_EQUAL(0, section->getFiberDimension(pointsOut[i]));
+} // testNewSectionPointsArray
+
+// ----------------------------------------------------------------------
// Test newSection(domain).
void
pylith::topology::TestFieldMesh::testNewSectionDomain(void)
@@ -935,7 +985,7 @@
v_iter != vertices->end();
++v_iter, ++iV) {
section->addConstraintDimension(*v_iter, nconstraints[iV]);
-#if 1 // TODO: FIELD SPLIT
+#if defined(FIELD_SPLIT)
section->addConstraintDimension(*v_iter, nconstraints[iV], fibration);
#endif
} // for
@@ -947,29 +997,27 @@
v_iter != vertices->end();
++v_iter, index += nconstraints[i++]) {
section->setConstraintDof(*v_iter, &constraints[index]);
-#if 1 // TODO: FIELD SPLIT
- std::cout << "nconstraints["<<i<<"]: " << nconstraints[i]
- << ", constraintDimension: " << section->getConstraintDimension(*v_iter) << std::endl;
+#if defined(FIELD_SPLIT)
section->setConstraintDof(*v_iter, &constraints[index], fibration);
#endif
} // for
} // Setup source field
-#if 0 // TODO: FIELD SPLIT
+#if defined(FIELD_SPLIT)
const ALE::Obj<Mesh::RealSection>& section = fieldSrc.section();
CPPUNIT_ASSERT(!section.isNull());
CPPUNIT_ASSERT_EQUAL(numFibrations, section->getNumSpaces());
const ALE::Obj<Mesh::RealSection>& sectionSplit = section->getFibration(0);
CPPUNIT_ASSERT(!sectionSplit.isNull());
- section->view("FULL FIELD");
- sectionSplit->view("FIBRATION 0");
+ section->view("FULL FIELD"); // TEMPORARY
+ sectionSplit->view("FIBRATION 0"); // TEMPORARY
CPPUNIT_ASSERT(!vertices.isNull());
int iV = 0;
for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
v_iter != vertices->end();
++v_iter) {
- CPPUNIT_ASSERT_EQUAL(fiberDim, sectionSplit->getFiberDimension(*v_iter));
+ CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(*v_iter, fibration));
CPPUNIT_ASSERT_EQUAL(nconstraints[iV++],
section->getConstraintDimension(*v_iter, fibration));
} // for
@@ -1013,9 +1061,8 @@
v_iter != vertices->end();
++v_iter) {
section->addConstraintDimension(*v_iter, nconstraints[iV++]);
-#if 0 // TODO: FIELD SPLIT
- section->addConstraintDimension(*v_iter, nconstraints[iV++],
- fibration);
+#if defined(FIELD_SPLIT)
+ section->addConstraintDimension(*v_iter, nconstraints[iV++], fibration);
#endif
} // for
fieldSrc.allocate();
@@ -1026,13 +1073,13 @@
v_iter != vertices->end();
++v_iter, index += nconstraints[i++]) {
section->setConstraintDof(*v_iter, &constraints[index]);
-#if 0 // TODO: FIELD SPLIT
+#if defined(FIELD_SPLIT)
section->setConstraintDof(*v_iter, &constraints[index], fibration);
#endif
} // for
} // Setup source field
-#if 0 // TODO: FIELD SPLIT
+#if defined(FIELD_SPLIT)
Field<Mesh> field(mesh);
field.cloneSection(fieldSrc);
const ALE::Obj<Mesh::RealSection>& section = field.section();
@@ -1040,17 +1087,17 @@
CPPUNIT_ASSERT_EQUAL(numFibrations, section->getNumSpaces());
const ALE::Obj<Mesh::RealSection>& sectionSplit = section->getFibration(0);
CPPUNIT_ASSERT(!sectionSplit.isNull());
- section->view("FULL FIELD");
- sectionSplit->view("FIBRATION 0");
+ section->view("FULL FIELD"); // TEMPORARY
+ sectionSplit->view("FIBRATION 0"); // TEMPORARY
+
int iV = 0;
for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
v_iter != vertices->end();
++v_iter) {
- CPPUNIT_ASSERT_EQUAL(fiberDim,
+ CPPUNIT_ASSERT_EQUAL(fiberDim,
section->getFiberDimension(*v_iter, fibration));
CPPUNIT_ASSERT_EQUAL(nconstraints[iV++],
- section->getConstraintDimension(*v_iter,
- fibration));
+ section->getConstraintDimension(*v_iter, fibration));
} // for
#endif
} // testCloneSectionSplit
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.hh 2009-06-05 19:26:11 UTC (rev 15124)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.hh 2009-06-06 00:19:51 UTC (rev 15125)
@@ -49,6 +49,7 @@
CPPUNIT_TEST( testSpaceDim );
CPPUNIT_TEST( testNewSection );
CPPUNIT_TEST( testNewSectionPoints );
+ CPPUNIT_TEST( testNewSectionPointsArray );
CPPUNIT_TEST( testNewSectionDomain );
CPPUNIT_TEST( testNewSectionChart );
CPPUNIT_TEST( testCloneSection );
@@ -103,6 +104,9 @@
/// Test newSection(points).
void testNewSectionPoints(void);
+ /// Test newSection(int_array).
+ void testNewSectionPointsArray(void);
+
/// Test newSection(domain).
void testNewSectionDomain(void);
More information about the CIG-COMMITS
mailing list