[cig-commits] r7708 - in short/3D/PyLith/trunk/libsrc: bc feassemble
brad at geodynamics.org
brad at geodynamics.org
Wed Jul 18 18:20:36 PDT 2007
Author: brad
Date: 2007-07-18 18:20:35 -0700 (Wed, 18 Jul 2007)
New Revision: 7708
Added:
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh
Modified:
short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.cc
short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.hh
short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.cc
short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.hh
short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.icc
Log:
Worked on factoring out boundary orientation information from FaultCohesive into CellGeometry (should allow automatic setting of which orientation function to call). Removed orientation stuff from BoundaryCondition. Did some cleanup and added notes in Neumann. Started work on C++ implementation of AbsorbingDampers absorbing boundary condition.
Added: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2007-07-19 01:17:32 UTC (rev 7707)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2007-07-19 01:20:35 UTC (rev 7708)
@@ -0,0 +1,141 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "AbsorbingDampers.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::AbsorbingDampers::AbsorbingDampers(void)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::bc::AbsorbingDampers::~AbsorbingDampers(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Initialize boundary condition. Determine orienation and compute traction
+// vector at integration points.
+void
+pylith::bc::AbsorbingDampers::initialize(const ALE::Obj<ALE::Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs,
+ const double_array& upDir);
+{ // initialize
+ assert(0 != _quadrature);
+ assert(0 != _db);
+ assert(!mesh.isNull());
+ assert(0 != cs);
+
+ if (3 != upDir.size())
+ throw std::runtime_error("Up direction for boundary orientation must be "
+ "a vector with 3 components.");
+
+ // Extract submesh associated with boundary
+ const ALE::Obj<ALE::Mesh> boundaryMesh =
+ ALE::Selection<ALE::Mesh>::submesh(mesh, mesh->getIntSection(_label));
+ if (boundaryMesh.isNull()) {
+ std::ostringstream msg;
+ msg << "Could not construct boundary mesh for absorbing boundary "
+ << "condition '" << _label << "'.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ // Get 'surface' cells (1 dimension lower than top-level cells)
+ const ALE::Obj<sieve_type>& sieve = boundaryMesh->getSieve();
+ const ALE::Obj<Mesh::label_sequence>& cells = boundaryMesh->heightStratum(1);
+
+ // Get damping constants at each quadrature point
+ const int spaceDim = cs->spaceDim();
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int fiberDim = spaceDim * numQuadPts;
+ _dampingConsts = new real_section_type(mesh->comm(), mesh->debug());
+ assert(!_dampingConsts.isNull());
+ _dampingConsts->setFiberDimension(cells, fiberDim);
+ _dampingConsts->allocate();
+
+ const int numBasis = _quadrature->numBasis();
+ const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
+ const double_array& verticesRef = _quadrature->vertices();
+ const int jacobianSize = (cellDim > 0) ? spaceDim * cellDim : 1;
+ double_array jacobian(jacobianSize);
+ double jacobianDet = 0;
+ double_array cellVertices(numBasis*spaceDim);
+
+ // Set up orientation information
+ const int orientationSize = spaceDim * spaceDim;
+ double orientation[orientationSize];
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+
+ // open database with material property information
+ _db->open();
+ const char* valuesNames =
+ (cellDim > 0) ? { "vp", "vs", "density" } : { "vp", "density" };
+ const int numValues = (cellDim > 0) ? 3 : 2;
+ _db->queryVals((const char**) valueNames, numValues);
+
+ const Mesh::label_sequence::iterator cellsBegin = cells->begin();
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+ for(Mesh::label_sequence::iterator c_iter = cells->begin();
+ c_iter != cells->end();
+ ++c_iter) {
+ _quadrature->computeGeometry(mesh, coordinates, *c_iter);
+ for(int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
+ } // for
+ } // for
+
+ _db->close();
+} // initialize
+
+// ----------------------------------------------------------------------
+// Integrate contributions to residual term (r) for operator.
+void
+pylith::bc::AbsorbingDampers::integrateResidual(
+ const ALE::Obj<real_section_type>& residual,
+ const ALE::Obj<Mesh>& mesh)
+{ // integrateResidual
+ throw std::logic_error("Not implemented.");
+} // integrateResidual
+
+// ----------------------------------------------------------------------
+// Integrate contributions to Jacobian matrix (A) associated with
+void
+pylith::bc::AbsorbingDampers::integrateJacobian(
+ PetscMat* mat,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh)
+{ // integrateJacobian
+ throw std::logic_error("Not implemented.");
+} // integrateJacobian
+
+// ----------------------------------------------------------------------
+// Verify configuration is acceptable.
+void
+pylith::bc::AbsorbingDampers::verifyConfiguration(const ALE::Obj<Mesh>& mesh)
+{ // verifyConfiguration
+ throw std::logic_error("Not implemented.");
+} // verifyConfiguration
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh 2007-07-19 01:17:32 UTC (rev 7707)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh 2007-07-19 01:20:35 UTC (rev 7708)
@@ -0,0 +1,137 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/bc/AbsorbingDampers.hh
+ *
+ * @brief C++ implementation of an absorbing boundary with simple
+ * dampers.
+ *
+ * Computes contributions to terms A and r in
+ *
+ * A(t) u(t+dt) = b(u(t), u(t-dt)),
+ *
+ * r = b - a u0(t+dt)
+ *
+ * where A(t) is a sparse matrix or vector, u(t+dt) is the field we
+ * want to compute at time t+dt, b is a vector that depends on the
+ * field at time t and t-dt, and u0 is zero at unknown DOF and set to
+ * the known values at the constrained DOF.
+ *
+ * Contributions from elasticity include the intertial and stiffness
+ * terms, so this object computes the following portions of A and r:
+ *
+ * A = 1/(dt*dt) [M]
+ *
+ * r = (1/(dt*dt) [M])(- {u(t+dt)} + 2/(dt*dt){u(t)} - {u(t-dt)}) - [K]{u(t)}
+ *
+ * Translational inertia.
+ * - Residual action over cell
+ * \f[
+ * \int_{V^e} \rho N^p \sum_q N^q u_i^q \, dV
+ * \f]
+ * - Jacobian action over cell
+ * \f[
+ * \int_{V^e} (\rho N^q N^q)_i \, dV
+ * \f]
+ *
+ * See governing equations section of user manual for more
+ * information.
+ */
+
+#if !defined(pylith_bc_absorbingdampers_hh)
+#define pylith_bc_absorbingdampers_hh
+
+#include "BoundaryCondition.hh" // ISA BoundaryCondition
+#include "pylith/feassemble/Integrator.hh" // ISA Integrator
+
+#include "pylith/utils/array.hh" // USES std::vector, double_array, int_array
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace bc {
+ class AbsorbingDampers;
+ class TestAbsorbingDampers; // unit testing
+ } // bc
+} // pylith
+
+
+/// C++ implementation of AbsorbingDampers boundary conditions.
+class pylith::bc::AbsorbingDampers : public BoundaryCondition,
+ public feassemble::Integrator
+{ // class AbsorbingDampers
+ friend class TestAbsorbingDampers; // unit testing
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Default constructor.
+ AbsorbingDampers(void);
+
+ /// Destructor.
+ ~AbsorbingDampers(void);
+
+ /** Initialize boundary condition.
+ *
+ * @param mesh PETSc mesh
+ * @param cs Coordinate system for mesh
+ * @param upDir Direction perpendicular to horizontal surface tangent
+ * direction that is not collinear with surface normal.
+ */
+ void initialize(const ALE::Obj<ALE::Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs,
+ const double_array& upDir);
+
+ /** Integrate contributions to residual term (r) for operator.
+ *
+ * @param residual Field containing values for residual
+ * @param mesh Finite-element mesh
+ */
+ void integrateResidual(const ALE::Obj<real_section_type>& residual,
+ const ALE::Obj<Mesh>& mesh);
+
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
+ *
+ * @param mat Sparse matrix
+ * @param t Current time
+ * @param fields Solution fields
+ * @param mesh Finite-element mesh
+ */
+ void integrateJacobian(PetscMat* mat,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh);
+
+ /** Verify configuration is acceptable.
+ *
+ * @param mesh Finite-element mesh
+ */
+ void verifyConfiguration(const ALE::Obj<Mesh>& mesh);
+
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ /// Not implemented
+ AbsorbingDampers(const AbsorbingDampers& m);
+
+ /// Not implemented
+ const AbsorbingDampers& operator=(const AbsorbingDampers& m);
+
+ // PRIVATE MEMBERS ////////////////////////////////////////////////////
+private :
+
+}; // class AbsorbingDampers
+
+#endif // pylith_bc_absorbingdampers_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.cc 2007-07-19 01:17:32 UTC (rev 7707)
+++ short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.cc 2007-07-19 01:20:35 UTC (rev 7708)
@@ -30,109 +30,5 @@
_db = 0;
} // destructor
-// ----------------------------------------------------------------------
-// Compute weighted orientation of cell face for 1-D elements.
-void
-pylith::bc::BoundaryCondition::_orient1D(double_array* orientation,
- const double_array& jacobian,
- const double jacobianDet,
- const double_array& upDir)
-{ // _orient1D
- assert(0 != orientation);
- assert(1 == orientation->size());
- (*orientation) = 1.0;
-} // _orient1D
-
-// ----------------------------------------------------------------------
-// Compute weighted orientation of cell face for 2-D elements.
-void
-pylith::bc::BoundaryCondition::_orient2D(double_array* orientation,
- const double_array& jacobian,
- const double jacobianDet,
- const double_array& upDir)
-{ // _orient2D
- const int orientSize = 4;
- assert(0 != orientation);
- assert(orientSize == orientation->size());
- const int jacobianSize = 2;
- assert(jacobianSize == jacobian.size());
- // cellDim is 1
- const int spaceDim = 2;
-
- const double j1 = jacobian[0];
- const double j2 = jacobian[1];
- (*orientation)[0] = j1;
- (*orientation)[1] = j2;
- (*orientation)[2] = j2;
- (*orientation)[3] = -j1;
-} // _orient2D
-
-// ----------------------------------------------------------------------
-// Compute weighted orientation of cell face for 3-D elements.
-void
-pylith::bc::BoundaryCondition::_orient3D(double_array* orientation,
- const double_array& jacobian,
- const double jacobianDet,
- const double_array& upDir)
-{ // _orient3D
- const int orientSize = 9;
- assert(0 != orientation);
- assert(orientSize == orientation->size());
- const int jacobianSize = 6;
- assert(jacobianSize == jacobian.size());
- assert(3 == upDir.size());
-
- const int cellDim = 2;
- const int spaceDim = 3;
-
- const double j00 = jacobian[0];
- const double j01 = jacobian[1];
- const double j10 = jacobian[2];
- const double j11 = jacobian[3];
- const double j20 = jacobian[4];
- const double j21 = jacobian[5];
-
- // Compute normal using Jacobian
- double r0 = j10*j21 - j20*j11;
- double r1 = -j00*j21 + j20*j01;
- double r2 = j00*j11 - j10*j01;
- // Make unit vector
- double mag = sqrt(r0*r0 + r1*r1 + r2*r2);
- assert(mag > 0.0);
- r0 /= mag;
- r1 /= mag;
- r2 /= mag;
-
- // Compute along-strike direction; cross product of "up" and normal
- double p0 = upDir[1]*r2 - upDir[2]*r1;
- double p1 = -upDir[0]*r2 + upDir[2]*r0;
- double p2 = upDir[0]*r1 - upDir[1]*r0;
- // Make unit vector
- mag = sqrt(p0*p0 + p1*p1 + p2*p2);
- assert(mag > 0.0);
- p0 /= mag;
- p1 /= mag;
- p2 /= mag;
-
- // Compute up-dip direction; cross product of normal and along-strike
- const double q0 = r1*p2 - r2*p1;
- const double q1 = -r0*p2 + r2*p0;
- const double q2 = r0*p1 - r1*p0;
- mag = sqrt(q0*q0 + q1*q1 + q2*q2);
- assert(mag > 0.0);
-
- const double wt = jacobianDet;
- (*orientation)[0] = p0*wt;
- (*orientation)[1] = p1*wt;
- (*orientation)[2] = p2*wt;
- (*orientation)[3] = q0*wt;
- (*orientation)[4] = q1*wt;
- (*orientation)[5] = q2*wt;
- (*orientation)[6] = r0*wt;
- (*orientation)[7] = r1*wt;
- (*orientation)[8] = r2*wt;
-} // _orient3D
-
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.hh 2007-07-19 01:17:32 UTC (rev 7707)
+++ short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.hh 2007-07-19 01:20:35 UTC (rev 7708)
@@ -99,84 +99,6 @@
void initialize(const ALE::Obj<ALE::Mesh>& mesh,
const spatialdata::geocoords::CoordSys* cs) = 0;
- // PROTECTED TYPEDEFS /////////////////////////////////////////////////
-protected :
-
- /// Function type for orientation methods.
- typedef void (*orient_fn_type)(double_array*,
- const double_array&,
- const double,
- const double_array&);
-
- // PROTECTED METHODS //////////////////////////////////////////////////
-protected :
-
- /** Compute weighted orientation of cell face for
- * 1-D elements. Orientation is either at vertices or quadrature
- * points, depending on whether the arguments have been evaluated at
- * the vertices or quadrature points.
- *
- * The orientation is returned as an array of direction cosines.
- *
- * size = spaceDim*spaceDim
- * index = iDir*spaceDim + iComponent
- *
- * @param orientation Array of direction cosines.
- * @param jacobian Jacobian matrix at point.
- * @param jacobianDet Determinant of Jacobian matrix at point.
- * @param upDir Direction perpendicular to horizontal tangent that
- * is not collinear with face normal (usually "up" direction).
- */
- static
- void _orient1D(double_array* orientation,
- const double_array& jacobian,
- const double jacobianDet,
- const double_array& upDir);
-
- /** Compute weighted orientation of cell face for
- * 2-D elements. Orientation is either at vertices or quadrature
- * points, depending on whether the arguments have been evaluated at
- * the vertices or quadrature points.
- *
- * The orientation is returned as an array of direction cosines.
- *
- * size = spaceDim*spaceDim
- * index = iDir*spaceDim + iComponent
- *
- * @param orientation Array of direction cosines.
- * @param jacobian Jacobian matrix at point.
- * @param jacobianDet Determinant of Jacobian matrix at point.
- * @param upDir Direction perpendicular to horizontal tangent that is
- * not collinear with face normal (usually "up" direction).
- */
- static
- void _orient2D(double_array* orientation,
- const double_array& jacobian,
- const double jacobianDet,
- const double_array& upDir);
-
- /** Compute weighted orientation of cell face for
- * 3-D elements. Orientation is either at vertices or quadrature
- * points, depending on whether the arguments have been evaluated at
- * the vertices or quadrature points.
- *
- * The orientation is returned as an array of direction cosines.
- *
- * size = spaceDim*spaceDim
- * index = iDir*spaceDim + iComponent
- *
- * @param orientation Array of direction cosines.
- * @param jacobian Jacobian matrix at point.
- * @param jacobianDet Determinant of Jacobian matrix at point.
- * @param upDir Direction perpendicular to horizontal tangent that is
- * not collinear with face normal (usually "up" direction).
- */
- static
- void _orient3D(double_array* orientation,
- const double_array& jacobian,
- const double jacobianDet,
- const double_array& upDir);
-
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2007-07-19 01:17:32 UTC (rev 7707)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2007-07-19 01:20:35 UTC (rev 7708)
@@ -34,13 +34,12 @@
} // destructor
// ----------------------------------------------------------------------
-// Initialize boundary condition. Determine orienation and compute stress
+// Initialize boundary condition. Determine orienation and compute traction
// vector at integration points.
void
pylith::bc::Neumann::initialize(const ALE::Obj<ALE::Mesh>& mesh,
const spatialdata::geocoords::CoordSys* cs,
- const double_array& upDir,
- const double_array& normalDir);
+ const double_array& upDir);
{ // initialize
assert(0 != _quadrature);
assert(0 != _db);
@@ -50,226 +49,108 @@
if (3 != upDir.size())
throw std::runtime_error("Up direction for surface orientation must be "
"a vector with 3 components.");
- if (3 != normalDir.size())
- throw std::runtime_error("Normal direction for surface orientation must be "
- "a vector with 3 components.");
// Extract submesh associated with surface
- const ALE::Obj<ALE::Mesh> submesh =
+ const ALE::Obj<ALE::Mesh> boundaryMesh =
ALE::Selection<ALE::Mesh>::submesh(mesh, mesh->getIntSection(_label));
- if (submesh.isNull()) {
+ if (boundaryMesh.isNull()) {
std::ostringstream msg;
- msg << "Could not find group of points '" << _label << "' in mesh.";
+ msg << "Could not construct boundary mesh for Neumann traction "
+ << "boundary condition '" << _label << "'.";
throw std::runtime_error(msg.str());
} // if
- assert(!submesh.isNull());
// Get 'surface' cells (1 dimension lower than top-level cells)
- const ALE::Obj<sieve_type>& sieve = submesh->getSieve();
- const ALE::Obj<Mesh::label_sequence>& cells = submesh->heightStratum(1);
+ const ALE::Obj<sieve_type>& sieve = boundaryMesh->getSieve();
+ const ALE::Obj<Mesh::label_sequence>& cells = boundaryMesh->heightStratum(1);
const int numCells = cells->size();
- // Create section for stress vector in global coordinates
+ // Create section for traction vector in global coordinates
const int spaceDim = cs->spaceDim();
const int numQuadPts = _quadrature->numQuadPts();
- const int stressVecSize = spaceDim * numQuadPts;
- // Not sure this part is correct
- _stressVecGlobal = new real_section_type(mesh->comm(), mesh->debug());
- assert(!_stressVecGlobal.isNull());
- const Mesh::label_sequence::iterator cellsBegin = cells->begin();
- const Mesh::label_sequence::iterator cellsEnd = cells->end();
- const int numCorners = sieve->nCone(*c_iter, mesh->depth())->size();
+ const int fiberDim = spaceDim * numQuadPts;
+ _tractionGlobal = new real_section_type(mesh->comm(), mesh->debug());
+ assert(!_tractionGlobal.isNull());
+ _tractionGlobal->setFiberDimension(cells, fiberDim);
+ _tractionGlobal->allocate();
- for(Mesh::label_sequence::iterator c_iter = cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- _stressVecGlobal->setFiberDimension(*c_iter, stressVecSize);
- }
- mesh->allocate(_stressVecGlobal);
+ const int numBasis = _quadrature->numBasis();
+ const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
+ const double_array& verticesRef = _quadrature->vertices();
+ const int jacobianSize = (cellDim > 0) ? spaceDim * cellDim : 1;
+ double_array jacobian(jacobianSize);
+ double jacobianDet = 0;
+ double_array cellVertices(numBasis*spaceDim);
// Set up orientation information
const int orientationSize = spaceDim * spaceDim;
double orientation[orientationSize];
- const int surfaceDim = mesh->getDimension()-1;
- orient_fn_type orientFn;
const ALE::Obj<real_section_type>& coordinates =
mesh->getRealSection("coordinates");
- switch (surfaceDim)
- { // switch
- case 0 :
- orientFn = _orient1D;
- break;
- case 1 :
- orientFn = _orient2D;
- break;
- case 2 :
- orientFn = _orient3D;
- break;
- default :
- assert(0);
- } // switch
- // Loop over cell faces, compute orientations, and then compute
- // corresponding stress vector in global coordinates.
- // Store values in _stressVecGlobal.
+ // setup database with traction information
+ // NEED TO SET NAMES BASED ON DIMENSION OF BOUNDARY
+ // 1-D problem = {'normal-traction'}
+ // 2-D problem = {'shear-traction', 'normal-traction'}
+ // 3-D problem = {'horiz-shear-traction', 'vert-shear-traction', 'normal-traction'}
+ _db->open();
+ _db->queryVals((const char**) valueNames, numFixedDOF);
+
+ // Loop over cells in boundary mesh, compute orientations, and then
+ // compute corresponding traction vector in global coordinates
+ // (store values in _tractionGlobal).
+ const Mesh::label_sequence::iterator cellsBegin = cells->begin();
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
for(Mesh::label_sequence::iterator c_iter = cells->begin();
c_iter != cells->end();
++c_iter) {
_quadrature->computeGeometry(mesh, coordinates, *c_iter);
for(int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
- // FINISH FIXING FROM HERE
+ // Get traction vector in local coordinate system at quadrature point
- 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 orientation of boundary at quadrature point
+ // Call restrict with cell to get face vertices
+ // Compute Jacobian at quadrature point
+ // Compute orientation using Jacobian
- // Get values for degrees of freedom
- 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;
+ // Rotate traction vector from local coordinate system to global
+ // coordinate system
- 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];
+ // Store traction vector in global coordinate system
+ } // for
} // for
+
_db->close();
} // initialize
// ----------------------------------------------------------------------
-// Set number of degrees of freedom that are constrained at points in field.
+// Integrate contributions to residual term (r) for operator.
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());
+pylith::bc::Neumann::integrateResidual(
+ const ALE::Obj<real_section_type>& residual,
+ const ALE::Obj<Mesh>& mesh)
+{ // integrateResidual
+ throw std::logic_error("Not implemented.");
+} // integrateResidual
- const int numFixedDOF = _fixedDOF.size();
- if (0 == numFixedDOF)
- return;
-
- const int numPoints = _points.size();
- _offsetLocal.resize(numPoints);
- for (int iPoint=0; iPoint < numPoints; ++iPoint) {
- const int fiberDim = field->getFiberDimension(_points[iPoint]);
- const int curNumConstraints = field->getConstraintDimension(_points[iPoint]);
- if (curNumConstraints + numFixedDOF > fiberDim) {
- std::ostringstream msg;
- msg << "Found overly constrained point while setting up constraints for Dirichlet "
- << "boundary condition '" << _label << "'.\n" << "Number of DOF at point "
- << _points[iPoint] << " is " << fiberDim << " and number of attempted constraints is "
- << curNumConstraints+numFixedDOF << ".";
- throw std::runtime_error(msg.str());
- } // if
- _offsetLocal[iPoint] = curNumConstraints;
- field->addConstraintDimension(_points[iPoint], numFixedDOF);
- } // for
-} // setConstraintSizes
-
// ----------------------------------------------------------------------
-// Set which degrees of freedom are constrained at points in field.
+// Integrate contributions to Jacobian matrix (A) associated with
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());
+pylith::bc::Neumann::integrateJacobian(PetscMat* mat,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh)
+{ // integrateJacobian
+ throw std::logic_error("Not implemented.");
+} // integrateJacobian
- const int numFixedDOF = _fixedDOF.size();
- if (0 == numFixedDOF)
- return;
-
- const int numPoints = _points.size();
- for (int iPoint=0; iPoint < numPoints; ++iPoint) {
- const Mesh::point_type point = _points[iPoint];
-
- // Get list of currently constrained DOF
- const int* curFixedDOF = field->getConstraintDof(point);
- const int numTotalConstrained = field->getConstraintDimension(point);
-
- // Create array holding all constrained DOF
- int_array allFixedDOF(curFixedDOF, numTotalConstrained);
-
- // Add in the ones for this Dirichlet BC
- for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
- allFixedDOF[_offsetLocal[iPoint]+iDOF] = _fixedDOF[iDOF];
-
- // Fill in rest of values not yet set (will be set by another Dirichlet BC)
- for (int iDOF=_offsetLocal[iPoint]+numFixedDOF;
- iDOF < numTotalConstrained;
- ++iDOF)
- allFixedDOF[_offsetLocal[iPoint]+iDOF] = 999;
-
- // Sort list of constrained DOF
- // I need these sorted for my update algorithms to work properly
- std::sort(&allFixedDOF[0], &allFixedDOF[numTotalConstrained]);
-
- // Update list of constrained DOF
- field->setConstraintDof(point, &allFixedDOF[0]);
- } // for
-} // setConstraints
-
// ----------------------------------------------------------------------
-// Set values in field.
+// Verify configuration is acceptable.
void
-pylith::bc::Dirichlet::setField(const double t,
- const ALE::Obj<real_section_type>& field,
- const ALE::Obj<ALE::Mesh>& mesh)
-{ // setField
- assert(!field.isNull());
- assert(!mesh.isNull());
+pylith::bc::Neumann::verifyConfiguration(const ALE::Obj<Mesh>& mesh)
+{ // verifyConfiguration
+ throw std::logic_error("Not implemented.");
+} // verifyConfiguration
- const int numFixedDOF = _fixedDOF.size();
- if (0 == numFixedDOF)
- return;
- const int numPoints = _points.size();
- for (int iPoint=0, i=0; iPoint < numPoints; ++iPoint) {
- const Mesh::point_type point = _points[iPoint];
- const int fiberDimension = field->getFiberDimension(point);
- double_array allValues(fiberDimension);
- mesh->restrict(field, point, &allValues[0], fiberDimension);
- for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
- allValues[_fixedDOF[iDOF]] = _values[i++];
- field->updatePointAll(_points[iPoint], &allValues[0]);
- } // for
-} // setField
-
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.hh 2007-07-19 01:17:32 UTC (rev 7707)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.hh 2007-07-19 01:20:35 UTC (rev 7708)
@@ -52,16 +52,12 @@
*
* @param mesh PETSc mesh
* @param cs Coordinate system for mesh
- * @param upDir Direction perpendicular to surface tangent direction that
- * is not collinear with surface normal. Surface tangent direction is
- * generally meant to be the horizontal direction, but this will not be
- * true for a horizontal surface. In that case it lies along the x-axis.
- * @param normalDir Outward-directed surface normal.
+ * @param upDir Direction perpendicular to horizontal surface tangent
+ * direction that is not collinear with surface normal.
*/
void initialize(const ALE::Obj<ALE::Mesh>& mesh,
const spatialdata::geocoords::CoordSys* cs,
- const double_array& upDir,
- const double_array& normalDir);
+ const double_array& upDir);
/** Integrate contributions to residual term (r) for operator.
*
@@ -71,20 +67,39 @@
void integrateResidual(const ALE::Obj<real_section_type>& residual,
const ALE::Obj<Mesh>& mesh);
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
+ *
+ * @param mat Sparse matrix
+ * @param t Current time
+ * @param fields Solution fields
+ * @param mesh Finite-element mesh
+ */
+ void integrateJacobian(PetscMat* mat,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh);
+
+ /** Verify configuration is acceptable.
+ *
+ * @param mesh Finite-element mesh
+ */
+ void verifyConfiguration(const ALE::Obj<Mesh>& mesh);
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
/// Not implemented
- Neumann(const Neumann& m);
+ Neumann(const Neumann&);
/// Not implemented
- const Neumann& operator=(const Neumann& m);
+ const Neumann& operator=(const Neumann&);
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
- /// Stress vector in global coordinates at integration points
- ALE::Obj<real_section_type> _stressVecGlobal;
+ /// Traction vector in global coordinates at integration points.
+ ALE::Obj<real_section_type> _tractionGlobal;
}; // class Neumann
Modified: short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.cc 2007-07-19 01:17:32 UTC (rev 7707)
+++ short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.cc 2007-07-19 01:20:35 UTC (rev 7708)
@@ -14,15 +14,39 @@
#include "CellGeometry.hh" // implementation of class methods
+#include "pylith/utils/array.hh" // USES double_array
+
+#include <iostream> // USES std::cerr
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::feassemble::CellGeometry::CellGeometry(const int cellDim,
- const int spaceDim,
- const int numCorners) :
+ const int spaceDim,
+ const int numCorners) :
+ _orientFn(0),
_cellDim(cellDim),
_spaceDim(spaceDim),
_numCorners(numCorners)
{ // constructor
+ switch (cellDim)
+ { // switch
+ case 0 :
+ _orientFn = _orient1D;
+ break;
+ case 1 :
+ _orientFn = _orient2D;
+ break;
+ case 2 :
+ _orientFn = _orient3D;
+ break;
+ case 3 :
+ break;
+ default:
+ std::cerr
+ << "Could not find orientation function for cell with dimension "
+ << cellDim << ".";
+ assert(0);
+ } // switch
} // constructor
// ----------------------------------------------------------------------
@@ -34,11 +58,116 @@
// ----------------------------------------------------------------------
// Copy constructor.
pylith::feassemble::CellGeometry::CellGeometry(const CellGeometry& g) :
+ _orientFn(g._orientFn),
_cellDim(g._cellDim),
_spaceDim(g._spaceDim),
_numCorners(g._numCorners)
{ // copy constructor
} // copy constructor
+// ----------------------------------------------------------------------
+// Compute weighted orientation of boundary for 1-D cell.
+void
+pylith::feassemble::CellGeometry::_orient1D(double_array* orientation,
+ const double_array& jacobian,
+ const double jacobianDet,
+ const double_array& upDir)
+{ // _orient1D
+ assert(0 != orientation);
+ assert(1 == orientation->size());
+ (*orientation) = 1.0;
+} // _orient1D
+
+// ----------------------------------------------------------------------
+// Compute weighted orientation of boundary for 2-D cell.
+void
+pylith::feassemble::CellGeometry::_orient2D(double_array* orientation,
+ const double_array& jacobian,
+ const double jacobianDet,
+ const double_array& upDir)
+{ // _orient2D
+ const int orientSize = 4;
+ assert(0 != orientation);
+ assert(orientSize == orientation->size());
+ const int jacobianSize = 2;
+ assert(jacobianSize == jacobian.size());
+ // cellDim is 1
+ const int spaceDim = 2;
+
+ const double j1 = jacobian[0];
+ const double j2 = jacobian[1];
+ (*orientation)[0] = j1;
+ (*orientation)[1] = j2;
+ (*orientation)[2] = j2;
+ (*orientation)[3] = -j1;
+} // _orient2D
+
+// ----------------------------------------------------------------------
+// Compute weighted orientation of boundary for 3-D cell.
+void
+pylith::feassemble::CellGeometry::_orient3D(double_array* orientation,
+ const double_array& jacobian,
+ const double jacobianDet,
+ const double_array& upDir)
+{ // _orient3D
+ const int orientSize = 9;
+ assert(0 != orientation);
+ assert(orientSize == orientation->size());
+ const int jacobianSize = 6;
+ assert(jacobianSize == jacobian.size());
+ assert(3 == upDir.size());
+
+ const int cellDim = 2;
+ const int spaceDim = 3;
+
+ const double j00 = jacobian[0];
+ const double j01 = jacobian[1];
+ const double j10 = jacobian[2];
+ const double j11 = jacobian[3];
+ const double j20 = jacobian[4];
+ const double j21 = jacobian[5];
+
+ // Compute normal using Jacobian
+ double r0 = j10*j21 - j20*j11;
+ double r1 = -j00*j21 + j20*j01;
+ double r2 = j00*j11 - j10*j01;
+ // Make unit vector
+ double mag = sqrt(r0*r0 + r1*r1 + r2*r2);
+ assert(mag > 0.0);
+ r0 /= mag;
+ r1 /= mag;
+ r2 /= mag;
+
+ // Compute along-strike direction; cross product of "up" and normal
+ double p0 = upDir[1]*r2 - upDir[2]*r1;
+ double p1 = -upDir[0]*r2 + upDir[2]*r0;
+ double p2 = upDir[0]*r1 - upDir[1]*r0;
+ // Make unit vector
+ mag = sqrt(p0*p0 + p1*p1 + p2*p2);
+ assert(mag > 0.0);
+ p0 /= mag;
+ p1 /= mag;
+ p2 /= mag;
+
+ // Compute up-dip direction; cross product of normal and along-strike
+ const double q0 = r1*p2 - r2*p1;
+ const double q1 = -r0*p2 + r2*p0;
+ const double q2 = r0*p1 - r1*p0;
+ mag = sqrt(q0*q0 + q1*q1 + q2*q2);
+ assert(mag > 0.0);
+
+ const double wt = jacobianDet;
+ (*orientation)[0] = p0*wt;
+ (*orientation)[1] = p1*wt;
+ (*orientation)[2] = p2*wt;
+ (*orientation)[3] = q0*wt;
+ (*orientation)[4] = q1*wt;
+ (*orientation)[5] = q2*wt;
+ (*orientation)[6] = r0*wt;
+ (*orientation)[7] = r1*wt;
+ (*orientation)[8] = r2*wt;
+} // _orient3D
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.hh 2007-07-19 01:17:32 UTC (rev 7707)
+++ short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.hh 2007-07-19 01:20:35 UTC (rev 7708)
@@ -93,6 +93,24 @@
const double_array& vertices,
const double_array& location) const = 0;
+ /** Compute orientation of cell at location.
+ *
+ * The orientation is returned as an array of direction cosines.
+ *
+ * size = spaceDim*spaceDim
+ * index = iDir*spaceDim + iComponent
+ *
+ * @param orientation Array of direction cosines.
+ * @param jacobian Jacobian matrix at point.
+ * @param jacobianDet Determinant of Jacobian matrix at point.
+ * @param upDir Direction perpendicular to horizontal direction that is
+ * not collinear with cell normal (usually "up" direction).
+ */
+ void orientation(double_array* orientation,
+ const double_array& jacobian,
+ const double jacobianDet,
+ const double_array& upDir);
+
// PROTECTED ////////////////////////////////////////////////////////////
protected :
@@ -108,13 +126,95 @@
/// Not implemented
const CellGeometry& operator=(const CellGeometry& );
+// PRIVATE TYPEDEFS /////////////////////////////////////////////////////
+private :
+
+ /// Function type for orientation methods.
+ typedef void (*orient_fn_type)(double_array*,
+ const double_array&,
+ const double,
+ const double_array&);
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+ /** Compute weighted orientation of fault for cohesive cell between
+ * 1-D elements. Orientation is either at vertices or quadrature
+ * points, depending on whether the arguments have been evaluated at
+ * the vertices or quadrature points.
+ *
+ * The orientation is returned as an array of direction cosines.
+ *
+ * size = spaceDim*spaceDim
+ * index = iDir*spaceDim + iComponent
+ *
+ * @param orientation Array of direction cosines.
+ * @param jacobian Jacobian matrix at point.
+ * @param jacobianDet Determinant of Jacobian matrix at point.
+ * @param upDir Direction perpendicular to along-strike direction that is
+ * not collinear with fault normal (usually "up" direction but could
+ * be up-dip direction).
+ */
+ static
+ void _orient1D(double_array* orientation,
+ const double_array& jacobian,
+ const double jacobianDet,
+ const double_array& upDir);
+
+ /** Compute weighted orientation of fault for cohesive cell between
+ * 2-D elements. Orientation is either at vertices or quadrature
+ * points, depending on whether the arguments have been evaluated at
+ * the vertices or quadrature points.
+ *
+ * The orientation is returned as an array of direction cosines.
+ *
+ * size = spaceDim*spaceDim
+ * index = iDir*spaceDim + iComponent
+ *
+ * @param orientation Array of direction cosines.
+ * @param jacobian Jacobian matrix at point.
+ * @param jacobianDet Determinant of Jacobian matrix at point.
+ * @param upDir Direction perpendicular to along-strike direction that is
+ * not collinear with fault normal (usually "up" direction but could
+ * be up-dip direction).
+ */
+ static
+ void _orient2D(double_array* orientation,
+ const double_array& jacobian,
+ const double jacobianDet,
+ const double_array& upDir);
+
+ /** Compute weighted orientation of fault for cohesive cell between
+ * 3-D elements. Orientation is either at vertices or quadrature
+ * points, depending on whether the arguments have been evaluated at
+ * the vertices or quadrature points.
+ *
+ * The orientation is returned as an array of direction cosines.
+ *
+ * size = spaceDim*spaceDim
+ * index = iDir*spaceDim + iComponent
+ *
+ * @param orientation Array of direction cosines.
+ * @param jacobian Jacobian matrix at point.
+ * @param jacobianDet Determinant of Jacobian matrix at point.
+ * @param upDir Direction perpendicular to along-strike direction that is
+ * not collinear with fault normal (usually "up" direction but could
+ * be up-dip direction).
+ */
+ static
+ void _orient3D(double_array* orientation,
+ const double_array& jacobian,
+ const double jacobianDet,
+ const double_array& upDir);
+
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
+ orient_fn_type _orientFn; ///< Function for computing orientation of cell
int _cellDim; ///< Dimension of cell.
int _spaceDim; ///< Dimension of coordinate space.
int _numCorners; ///< Number of corners in cell.
-
+
}; // CellGeometry
#include "CellGeometry.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.icc 2007-07-19 01:17:32 UTC (rev 7707)
+++ short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.icc 2007-07-19 01:20:35 UTC (rev 7708)
@@ -14,6 +14,8 @@
#error "CellGeometry.icc must be included only from CellGeometry.hh"
#else
+#include <assert.h> // USES assert()
+
// Get dimension of cell.
inline
int
@@ -35,7 +37,18 @@
return _numCorners;
} // numCorners
+// Compute orientation of cell at location.
+inline
+void
+pylith::feassemble::CellGeometry::orientation(double_array* orientation,
+ const double_array& jacobian,
+ const double jacobianDet,
+ const double_array& upDir) {
+ if (0 != _orientFn)
+ _orientFn(orientation, jacobian, jacobianDet, upDir);
+} // orientation
+
#endif
// End of file
More information about the cig-commits
mailing list