[cig-commits] r6909 - in short/3D/PyLith/trunk: . libsrc/bc
pylith/bc
brad at geodynamics.org
brad at geodynamics.org
Wed May 16 17:26:36 PDT 2007
Author: brad
Date: 2007-05-16 17:26:35 -0700 (Wed, 16 May 2007)
New Revision: 6909
Modified:
short/3D/PyLith/trunk/DEPENDENCIES
short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.cc
short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.hh
short/3D/PyLith/trunk/libsrc/bc/Dirichlet.cc
short/3D/PyLith/trunk/libsrc/bc/Dirichlet.hh
short/3D/PyLith/trunk/pylith/bc/BoundaryCondition.py
short/3D/PyLith/trunk/pylith/bc/Dirichlet.py
Log:
Initial C++ implementation of Dirichlet boundary conditions. Will need update if/when Matt fixes sections to allow overlapping Dirichlet BC domains. Updated Python implementation acccordingly. Still need bindings and unit tests.
Modified: short/3D/PyLith/trunk/DEPENDENCIES
===================================================================
--- short/3D/PyLith/trunk/DEPENDENCIES 2007-05-16 22:09:21 UTC (rev 6908)
+++ short/3D/PyLith/trunk/DEPENDENCIES 2007-05-17 00:26:35 UTC (rev 6909)
@@ -1,3 +1,10 @@
+TEMPORARY KLUDGES
+
+#define NEW_SECTION in $PETSC_DIR/include/petscmesh.h
+
+ Matt will change the section type to this new section type after his
+ meeting in Turkey.
+
REQUIRED DEPENDENCIES
MPI
@@ -10,5 +17,7 @@
OPTIONAL DEPENDENCIES
cppunit
+CUBIT
+PETSc w/TetGen
doxygen
Modified: short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.cc 2007-05-16 22:09:21 UTC (rev 6908)
+++ short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.cc 2007-05-17 00:26:35 UTC (rev 6909)
@@ -31,11 +31,20 @@
} // destructor
// ----------------------------------------------------------------------
-// Set constrained degrees of freedom in field.
+// Set number of degrees of freedom that are constrained at points in field.
void
-pylith::bc::BoundaryCondition::setConstraints(
+pylith::bc::BoundaryCondition::setConstraintSizes(
const ALE::Obj<real_section_type>& field,
const ALE::Obj<ALE::Mesh>& mesh)
+{ // setConstraintSizes
+} // setConstraintSizes
+
+// ----------------------------------------------------------------------
+// Set which degrees of freedom are constrained at points in field.
+void
+pylith::bc::BoundaryCondition::setConstraints(
+ const ALE::Obj<real_section_type>& field,
+ const ALE::Obj<ALE::Mesh>& mesh)
{ // setConstraints
} // setConstraints
Modified: short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.hh 2007-05-16 22:09:21 UTC (rev 6908)
+++ short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.hh 2007-05-17 00:26:35 UTC (rev 6909)
@@ -98,12 +98,21 @@
void initialize(const ALE::Obj<ALE::Mesh>& mesh,
const spatialdata::geocoords::CoordSys* cs) = 0;
- /** Set constrained degrees of freedom in field.
+ /** Set number of degrees of freedom that are constrained at points in field.
*
* @param field Solution field
* @param mesh PETSc mesh
*/
virtual
+ void setConstraintSizes(const ALE::Obj<real_section_type>& field,
+ const ALE::Obj<ALE::Mesh>& mesh);
+
+ /** Set which degrees of freedom are constrained at points in field.
+ *
+ * @param field Solution field
+ * @param mesh PETSc mesh
+ */
+ virtual
void setConstraints(const ALE::Obj<real_section_type>& field,
const ALE::Obj<ALE::Mesh>& mesh);
Modified: short/3D/PyLith/trunk/libsrc/bc/Dirichlet.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Dirichlet.cc 2007-05-16 22:09:21 UTC (rev 6908)
+++ short/3D/PyLith/trunk/libsrc/bc/Dirichlet.cc 2007-05-17 00:26:35 UTC (rev 6909)
@@ -14,6 +14,13 @@
#include "Dirichlet.hh" // implementation of object methods
+#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+
+#include <assert.h> // USES assert()
+#include <stdexcept> // USES std::runtime_error
+#include <sstream> // USES std::ostringstream
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::bc::Dirichlet::Dirichlet(void)
@@ -32,14 +39,98 @@
pylith::bc::Dirichlet::initialize(const ALE::Obj<ALE::Mesh>& mesh,
const spatialdata::geocoords::CoordSys* cs)
{ // initialize
+ assert(0 != _db);
+ assert(!mesh.isNull());
+ assert(0 != cs);
+
+ // Get points associated with boundary condition
+ const ALE::Obj<int_section_type>& groupField = mesh->getIntSection(_label);
+ if (groupField.isNull()) {
+ std::ostringstream msg;
+ msg << "Could not find group of points '" << _label << "' in mesh.";
+ throw std::runtime_error(msg.str());
+ } // if
+ assert(!groupField.isNull());
+ const int_section_type::chart_type& chart = groupField->getChart();
+ const int numPoints = chart.size();
+ _points.resize(numPoints);
+ int i = 0;
+ for(int_section_type::chart_type::iterator c_iter = chart.begin();
+ c_iter != chart.end();
+ ++c_iter)
+ _points[i] = *c_iter;
+
+ // Get values for degrees of freedom
+ const int numFixedDOF = _fixedDOF.size();
+ char** valueNames = (numFixedDOF > 0) ? new char*[numFixedDOF] : 0;
+ for (int i=0; i < numFixedDOF; ++i) {
+ std::ostringstream name;
+ name << "dof-" << _fixedDOF[i];
+ const int size = 1 + name.str().length();
+ valueNames[i] = new char[size];
+ strcpy(valueNames[i], name.str().c_str());
+ } // for
+ _db->open();
+ _db->queryVals((const char**) valueNames, numFixedDOF);
+ for (int i=0; i < numFixedDOF; ++i) {
+ delete[] valueNames[i]; valueNames[i] = 0;
+ } // for
+ delete[] valueNames; valueNames = 0;
+
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ const int spaceDim = cs->spaceDim();
+
+ _values.resize(numPoints*numFixedDOF);
+ double_array queryValues(numFixedDOF);
+ for (int iPoint=0, i=0; iPoint < numPoints; ++iPoint) {
+ // Get coordinates of vertex
+ const real_section_type::value_type* vCoords =
+ coordinates->restrictPoint(_points[iPoint]);
+ int err = _db->query(&queryValues[0], numFixedDOF, vCoords, spaceDim, cs);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Could not find values at (";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << vCoords[i];
+ msg << ") using spatial database " << _db->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
+ _values[i++] = queryValues[iDOF];
+ } // for
+ _db->close();
} // initialize
// ----------------------------------------------------------------------
-// Set constrained degrees of freedom in field.
+// Set number of degrees of freedom that are constrained at points in field.
void
+pylith::bc::Dirichlet::setConstraintSizes(const ALE::Obj<real_section_type>& field,
+ const ALE::Obj<ALE::Mesh>& mesh)
+{ // setConstraintSizes
+ assert(!field.isNull());
+ assert(!mesh.isNull());
+
+ const int numPoints = _points.size();
+ const int numFixedDOF = _fixedDOF.size();
+ for (int iPoint=0; iPoint < numPoints; ++iPoint)
+ field->setConstraintDimension(_points[iPoint], numFixedDOF);
+} // setConstraintSizes
+
+// ----------------------------------------------------------------------
+// Set which degrees of freedom are constrained at points in field.
+void
pylith::bc::Dirichlet::setConstraints(const ALE::Obj<real_section_type>& field,
const ALE::Obj<ALE::Mesh>& mesh)
{ // setConstraints
+ assert(!field.isNull());
+ assert(!mesh.isNull());
+
+ const int numPoints = _points.size();
+ const int numFixedDOF = _fixedDOF.size();
+ for (int iPoint=0; iPoint < numPoints; ++iPoint)
+ field->setConstraintDof(_points[iPoint], &_fixedDOF[0]);
} // setConstraints
// ----------------------------------------------------------------------
@@ -49,6 +140,13 @@
const ALE::Obj<real_section_type>& field,
const ALE::Obj<ALE::Mesh>& mesh)
{ // setField
+ assert(!field.isNull());
+ assert(!mesh.isNull());
+
+ const int numPoints = _points.size();
+ const int numFixedDOF = _fixedDOF.size();
+ for (int iPoint=0, i=0; iPoint < numPoints; ++iPoint, i+=numFixedDOF)
+ field->updatePointBC(_points[iPoint], &_values[i]);
} // setField
Modified: short/3D/PyLith/trunk/libsrc/bc/Dirichlet.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Dirichlet.hh 2007-05-16 22:09:21 UTC (rev 6908)
+++ short/3D/PyLith/trunk/libsrc/bc/Dirichlet.hh 2007-05-17 00:26:35 UTC (rev 6909)
@@ -65,11 +65,21 @@
void initialize(const ALE::Obj<ALE::Mesh>& mesh,
const spatialdata::geocoords::CoordSys* cs);
- /** Set constrained degrees of freedom in field.
+ /** Set number of degrees of freedom that are constrained at points in field.
*
* @param field Solution field
* @param mesh PETSc mesh
*/
+ virtual
+ void setConstraintSizes(const ALE::Obj<real_section_type>& field,
+ const ALE::Obj<ALE::Mesh>& mesh);
+
+ /** Set which degrees of freedom are constrained at points in field.
+ *
+ * @param field Solution field
+ * @param mesh PETSc mesh
+ */
+ virtual
void setConstraints(const ALE::Obj<real_section_type>& field,
const ALE::Obj<ALE::Mesh>& mesh);
Modified: short/3D/PyLith/trunk/pylith/bc/BoundaryCondition.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/BoundaryCondition.py 2007-05-16 22:09:21 UTC (rev 6908)
+++ short/3D/PyLith/trunk/pylith/bc/BoundaryCondition.py 2007-05-17 00:26:35 UTC (rev 6909)
@@ -88,9 +88,16 @@
return
+ def setConstraintSizes(self, field, mesh):
+ """
+ Set number of constraints at points in field.
+ """
+ return
+
+
def setConstraints(self, field, mesh):
"""
- Set constraints in field.
+ Set which degrees of freedom are constrained at points in field.
"""
return
Modified: short/3D/PyLith/trunk/pylith/bc/Dirichlet.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/Dirichlet.py 2007-05-16 22:09:21 UTC (rev 6908)
+++ short/3D/PyLith/trunk/pylith/bc/Dirichlet.py 2007-05-17 00:26:35 UTC (rev 6909)
@@ -91,9 +91,18 @@
return
+ def setConstraintSizes(self, field, mesh):
+ """
+ Set number of constraints at points in field.
+ """
+ assert(None != self.cppHandle)
+ self.cppHandle.setConstraintSizes(field)
+ return
+
+
def setConstraints(self, field, mesh):
"""
- Set constraints in field.
+ Set which degrees of freedom are constrained at points in field.
"""
assert(None != self.cppHandle)
self.cppHandle.setConstraints(field)
More information about the cig-commits
mailing list