[cig-commits] r7724 - in short/3D/PyLith/trunk: . libsrc/bc
brad at geodynamics.org
brad at geodynamics.org
Fri Jul 20 12:07:36 PDT 2007
Author: brad
Date: 2007-07-20 12:07:36 -0700 (Fri, 20 Jul 2007)
New Revision: 7724
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh
Log:
Finished initial implementation of AbsorbingDampers. Need unit tests and bindings.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-07-20 18:35:22 UTC (rev 7723)
+++ short/3D/PyLith/trunk/TODO 2007-07-20 19:07:36 UTC (rev 7724)
@@ -35,32 +35,6 @@
relative to orientation and if orientation coincides with
preexisting constraint, then constrain Lagrange multiplier DOF.
-----------------------------------------------------------------------
-List of missing features (EqSim/PyLith 0.8) and new features.
-
- Missing features
- EqSim
- Absorbing boundaries
- Dynamic fault interface conditions
- Fault constitutive models
- Output of surface/fault information
- PyLith 0.8
- Traction boundary conditions
- Viscoelastic material models (several)
- Output of stress/strain information
-
- New features
- EqSim
- Multiple cell types
- 1-D and 2-D simulations
- Importing of meshes from CUBIT
- Exporting to VTK files
- Availability as open-source with documentation
- PyLith 0.8
- Dislocation-based fault implementation
- Importing of meshes from CUBIT and LaGriT
- User-friendly specification of boundary conditions and parameters
-
======================================================================
MAIN PRIORITIES (Brad)
======================================================================
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2007-07-20 18:35:22 UTC (rev 7723)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2007-07-20 19:07:36 UTC (rev 7724)
@@ -16,6 +16,7 @@
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
#include "pylith/feassemble/CellGeometry.hh" // USES CellGeometry
+#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
@@ -56,9 +57,9 @@
"a vector with 3 components.");
// Extract submesh associated with boundary
- const ALE::Obj<ALE::Mesh> boundaryMesh =
+ _boundaryMesh =
ALE::Selection<ALE::Mesh>::submesh(mesh, mesh->getIntSection(_label));
- if (boundaryMesh.isNull()) {
+ if (_boundaryMesh.isNull()) {
std::ostringstream msg;
msg << "Could not construct boundary mesh for absorbing boundary "
<< "condition '" << _label << "'.";
@@ -66,11 +67,11 @@
} // if
// check compatibility of quadrature and boundary mesh
- if (_quadrature->cellDim() != boundaryMesh->getDimension()) {
+ if (_quadrature->cellDim() != _boundaryMesh->getDimension()) {
std::ostringstream msg;
msg << "Quadrature is incompatible with cells for absorbing boundary "
<< "condition '" << _label << "'.\n"
- << "Dimension of boundary mesh: " << boundaryMesh->getDimension()
+ << "Dimension of boundary mesh: " << _boundaryMesh->getDimension()
<< ", dimension of quadrature: " << _quadrature->cellDim()
<< ".";
throw std::runtime_error(msg.str());
@@ -79,7 +80,7 @@
// Get 'surface' cells (1 dimension lower than top-level cells)
const ALE::Obj<ALE::Mesh::label_sequence>& cells =
- boundaryMesh->heightStratum(1);
+ _boundaryMesh->heightStratum(1);
assert(!cells.isNull());
const Mesh::label_sequence::iterator cellsBegin = cells->begin();
const Mesh::label_sequence::iterator cellsEnd = cells->end();
@@ -148,7 +149,7 @@
for(Mesh::label_sequence::iterator c_iter = cells->begin();
c_iter != cells->end();
++c_iter) {
- _quadrature->computeGeometry(mesh, coordinates, *c_iter);
+ _quadrature->computeGeometry(_boundaryMesh, coordinates, *c_iter);
const double_array& quadPts = _quadrature->quadPts();
const double_array& quadPtsRef = _quadrature->quadPtsRef();
@@ -193,22 +194,194 @@
// 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)
+ const ALE::Obj<real_section_type>& residual,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh)
{ // integrateResidual
- throw std::logic_error("Not implemented.");
+ assert(0 != _quadrature);
+ assert(!_boundaryMesh.isNull());
+ assert(!residual.isNull());
+ assert(0 != fields);
+ assert(!mesh.isNull());
+
+ PetscErrorCode err = 0;
+
+ // Get cell information
+ const ALE::Obj<ALE::Mesh::label_sequence>& cells =
+ _boundaryMesh->heightStratum(1);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ const ALE::Obj<real_section_type>& dispTpdt = fields->getHistoryItem(0);
+ const ALE::Obj<real_section_type>& dispTmdt = fields->getHistoryItem(2);
+ assert(!dispTpdt.isNull());
+ assert(!dispTmdt.isNull());
+
+ // Get parameters used in integration.
+ const double dt = _dt;
+ assert(dt > 0);
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+
+ // Allocate vectors for cell values.
+ _initCellVector();
+ const int cellVecSize = numBasis*spaceDim;
+ double_array dispTpdtCell(cellVecSize);
+ double_array dispTmdtCell(cellVecSize);
+ double_array dampingConstsCell(numQuadPts*spaceDim);
+
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(_boundaryMesh, coordinates, *c_iter);
+
+ // Reset element vector to zero
+ _resetCellVector();
+
+ // Restrict input fields to cell
+ _boundaryMesh->restrict(dispTpdt, *c_iter, &dispTpdtCell[0], cellVecSize);
+ _boundaryMesh->restrict(dispTmdt, *c_iter, &dispTmdtCell[0], cellVecSize);
+ _boundaryMesh->restrict(_dampingConsts, *c_iter,
+ &dampingConstsCell[0], dampingConstsCell.size());
+
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
+ // Compute action for absorbing bc terms
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt =
+ quadWts[iQuad] * jacobianDet[iQuad] / (2.0 * dt);
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const double valI = wt*basis[iQuad*numBasis+iBasis];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const double valIJ = valI * basis[iQuad*numBasis+jBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _cellVector[iBasis*spaceDim+iDim] +=
+ dampingConstsCell[iQuad*spaceDim+iDim] *
+ valIJ * (- dispTpdtCell[jBasis*spaceDim+iDim]
+ + dispTmdtCell[jBasis*spaceDim+iDim]);
+ } // for
+ } // for
+ } // for
+ err = PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(5*spaceDim))));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ // Assemble cell contribution into field
+ mesh->updateAdd(residual, *c_iter, _cellVector);
+ } // for
} // integrateResidual
// ----------------------------------------------------------------------
// Integrate contributions to Jacobian matrix (A) associated with
void
pylith::bc::AbsorbingDampers::integrateJacobian(
- PetscMat* mat,
+ PetscMat* jacobian,
const double t,
topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh)
{ // integrateJacobian
- throw std::logic_error("Not implemented.");
+ assert(0 != _quadrature);
+ assert(!_boundaryMesh.isNull());
+ assert(0 != jacobian);
+ assert(0 != fields);
+ assert(!mesh.isNull());
+
+ PetscErrorCode err = 0;
+
+ // Get cell information
+ const ALE::Obj<ALE::Mesh::label_sequence>& cells =
+ _boundaryMesh->heightStratum(1);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ const ALE::Obj<real_section_type>& dispT = fields->getHistoryItem(1);
+ assert(!dispT.isNull());
+
+ // Get parameters used in integration.
+ const double dt = _dt;
+ assert(dt > 0);
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+
+ // Allocate vectors for cell values.
+ _initCellVector();
+ const int cellVecSize = numBasis*spaceDim;
+ double_array dampingConstsCell(numQuadPts*spaceDim);
+
+ // Allocate vector for cell values (if necessary)
+ _initCellMatrix();
+
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(_boundaryMesh, coordinates, *c_iter);
+
+ // Reset element vector to zero
+ _resetCellMatrix();
+
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
+ // Compute Jacobian for absorbing bc terms
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt =
+ quadWts[iQuad] * jacobianDet[iQuad] / (2.0 * dt);
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const double valI = wt*basis[iQ+iBasis];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const double valIJ = valI * basis[iQ+jBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ const int iBlock = (iBasis*spaceDim + iDim) * (numBasis*spaceDim);
+ const int jBlock = (jBasis*spaceDim + iDim);
+ _cellMatrix[iBlock+jBlock] +=
+ valIJ * dampingConstsCell[iQuad*spaceDim+iDim];
+ } // for
+ } // for
+ } // for
+ } // for
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(1+2*spaceDim))));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ // Assemble cell contribution into PETSc Matrix
+ const ALE::Obj<Mesh::order_type>& globalOrder =
+ mesh->getFactory()->getGlobalOrder(mesh, "default", dispT);
+ assert(!globalOrder.isNull());
+
+ err = updateOperator(*jacobian, mesh, dispT, globalOrder,
+ *c_iter, _cellMatrix, ADD_VALUES);
+ if (err)
+ throw std::runtime_error("Update to PETSc Mat failed.");
+ } // for
+
+ _needNewJacobian = false;
} // integrateJacobian
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh 2007-07-20 18:35:22 UTC (rev 7723)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh 2007-07-20 19:07:36 UTC (rev 7724)
@@ -29,21 +29,21 @@
* 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]
+ * A = 1/(2*dt) [C]
*
- * r = (1/(dt*dt) [M])(- {u(t+dt)} + 2/(dt*dt){u(t)} - {u(t-dt)}) - [K]{u(t)}
+ * r = (1/(2*dt) [C])(- {u(t+dt)} + {u(t-dt)})
*
- * Translational inertia.
+ * Damping matrix.
* - Residual action over cell
* \f[
- * \int_{V^e} \rho N^p \sum_q N^q u_i^q \, dV
+ * \int_{S^e} \rho c_i N^p \sum_q N^q u_i^q \, dS
* \f]
* - Jacobian action over cell
* \f[
- * \int_{V^e} (\rho N^q N^q)_i \, dV
+ * \int_{S^e} (\rho c_i N^q N^q)_i \, dS
* \f]
*
- * See governing equations section of user manual for more
+ * See boundary conditions section of user manual for more
* information.
*/
@@ -94,9 +94,13 @@
/** Integrate contributions to residual term (r) for operator.
*
* @param residual Field containing values for residual
+ * @param t Current time
+ * @param fields Solution fields
* @param mesh Finite-element mesh
*/
void integrateResidual(const ALE::Obj<real_section_type>& residual,
+ const double t,
+ topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh);
/** Integrate contributions to Jacobian matrix (A) associated with
@@ -130,6 +134,9 @@
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
+ /// Mesh of absorbing boundary
+ ALE::Obj<Mesh> _boundaryMesh;
+
/// Damping constants in global coordinates at integration points.
ALE::Obj<real_section_type> _dampingConsts;
More information about the cig-commits
mailing list