[cig-commits] r5104 - short/3D/PyLith/trunk/libsrc/feassemble
brad at geodynamics.org
brad at geodynamics.org
Thu Oct 26 16:37:16 PDT 2006
Author: brad
Date: 2006-10-26 16:37:16 -0700 (Thu, 26 Oct 2006)
New Revision: 5104
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia3D.cc
Log:
Continued implementing integration of interial action.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc 2006-10-26 23:36:47 UTC (rev 5103)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc 2006-10-26 23:37:16 UTC (rev 5104)
@@ -47,4 +47,54 @@
_quadrature = (0 != q) ? q->clone() : 0;
} // quadrature
+// ----------------------------------------------------------------------
+// Initialize vector containing result of integration action for cell.
+void
+pylith::feassemble::Integrator::_initCellVector(void)
+{ // _initCellVector
+ assert(0 != _quadrature);
+ const int size = _quadrature->spaceDim() * _quadrature->numCorners();
+ _cellVector = (size > 0) ? new real_section_type::value_type[size] : 0;
+ for (int i=0; i < size; ++i)
+ _cellVector[i] = 0.0;
+} // _initCellVector
+
+// ----------------------------------------------------------------------
+// Zero out vector containing result of integration actions for cell.
+void
+pylith::feassemble::Integrator::_resetCellVector(void)
+{ // _resetCellVector
+ assert(0 != _quadrature);
+ const int size = _quadrature->spaceDim() * _quadrature->numCorners();
+ for (int i=0; i < size; ++i)
+ _cellVector[i] = 0.0;
+} // _resetCellVector
+
+// ----------------------------------------------------------------------
+// Initialize matrix containing result of integration for cell.
+void
+pylith::feassemble::Integrator::_initCellMatrix(void)
+{ // _initCellMatrix
+ assert(0 != _quadrature);
+ const int size =
+ _quadrature->spaceDim() * _quadrature->numCorners() *
+ _quadrature->spaceDim() * _quadrature->numCorners();
+ _cellMatrix = (size > 0) ? new real_section_type::value_type[size] : 0;
+ for (int i=0; i < size; ++i)
+ _cellMatrix[i] = 0.0;
+} // _initCellMatrix
+
+// ----------------------------------------------------------------------
+// Zero out matrix containing result of integration for cell.
+void
+pylith::feassemble::Integrator::_resetCellMatrix(void)
+{ // _resetCellMatrix
+ assert(0 != _quadrature);
+ const int size =
+ _quadrature->spaceDim() * _quadrature->numCorners() *
+ _quadrature->spaceDim() * _quadrature->numCorners();
+ for (int i=0; i < size; ++i)
+ _cellMatrix[i] = 0.0;
+} // _resetCellMatrix
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2006-10-26 23:36:47 UTC (rev 5103)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2006-10-26 23:37:16 UTC (rev 5104)
@@ -35,6 +35,14 @@
{ // Integrator
friend class TestIntegrator; // unit testing
+// PUBLIC TYPEDEFS //////////////////////////////////////////////////////
+public :
+
+ typedef ALE::Mesh Mesh;
+ typedef Mesh::topology_type topology_type;
+ typedef topology_type::point_type point_type;
+ typedef Mesh::real_section_type real_section_type;
+
// PUBLIC MEMBERS ///////////////////////////////////////////////////////
public :
@@ -56,9 +64,9 @@
* @param coordinates Field of cell vertex coordinates
*/
virtual
- void integrateAction(const ALE::Obj<ALE::Mesh::real_section_type>& fieldOut,
- const ALE::Obj<ALE::Mesh::real_section_type>& fieldIn,
- const ALE::Obj<ALE::Mesh::real_section_type>& coordinates) = 0;
+ void integrateAction(const ALE::Obj<real_section_type>& fieldOut,
+ const ALE::Obj<real_section_type>& fieldIn,
+ const ALE::Obj<real_section_type>& coordinates) = 0;
/** Compute matrix associated with operator.
*
@@ -67,7 +75,7 @@
*/
virtual
void integrate(PetscMat* mat,
- const ALE::Obj<ALE::Mesh::real_section_type>& coordinates) = 0;
+ const ALE::Obj<real_section_type>& coordinates) = 0;
/** Set quadrature for integrating finite-element quantities.
*
@@ -84,17 +92,35 @@
*/
Integrator(const Integrator& i);
+ /// Initialize vector containing result of integration action for cell.
+ void _initCellVector(void);
+
+ /// Zero out vector containing result of integration actions for cell.
+ void _resetCellVector(void);
+
+ /// Initialize matrix containing result of integration for cell.
+ void _initCellMatrix(void);
+
+ /// Zero out matrix containing result of integration for cell.
+ void _resetCellMatrix(void);
+
// PRIVATE METHODS //////////////////////////////////////////////////////
private :
/// Not implemented
const Integrator& operator=(const Integrator&);
-// PRIVATE MEMBERS //////////////////////////////////////////////////////
-private :
+// PROTECTED MEMBERS ////////////////////////////////////////////////////
+protected :
Quadrature* _quadrature; ///< Quadrature for integrating finite-element
+ /// Vector local to cell containing result of integration action
+ real_section_type::value_type* _cellVector;
+
+ /// Matrix local to cell containing result of integration
+ real_section_type::value_type* _cellMatrix;
+
}; // Integrator
#endif // pylith_feassemble_integrator_hh
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia3D.cc 2006-10-26 23:36:47 UTC (rev 5103)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia3D.cc 2006-10-26 23:37:16 UTC (rev 5104)
@@ -14,6 +14,8 @@
#include "IntegratorInertia3D.hh" // implementation of class methods
+#include "Quadrature.hh" // USES Quadrature
+
#include "petscmat.h" // USES PetscMat
#include <assert.h> // USES assert()
@@ -40,17 +42,65 @@
// ----------------------------------------------------------------------
// Integrate inertial term for 3-D finite elements.
void
-pylith::feassemble::IntegratorInertia3D::integrateAction(const ALE::Obj<ALE::Mesh::real_section_type>& fieldOut,
- const ALE::Obj<ALE::Mesh::real_section_type>& fieldIn,
- const ALE::Obj<ALE::Mesh::real_section_type>& coordinates)
+pylith::feassemble::IntegratorInertia3D::integrateAction(
+ const ALE::Obj<real_section_type>& fieldOut,
+ const ALE::Obj<real_section_type>& fieldIn,
+ const ALE::Obj<real_section_type>& coordinates)
{ // integrateAction
+ assert(0 != _quadrature);
+
+ const topology_type::patch_type patch = 0;
+ const ALE::Obj<topology_type>& topology = fieldIn->getTopology();
+ const ALE::Obj<topology_type::label_sequence>& cells =
+ topology->heightStratum(patch, 0);
+ const topology_type::label_sequence::iterator cellsEnd = cells->end();
+
+ _initCellVector();
+
+ for (topology_type::label_sequence::iterator cellIter=cells->begin();
+ cellIter != cellsEnd;
+ ++cellIter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(coordinates, *cellIter);
+
+ // Reset element vector to zero
+ _resetCellVector();
+
+ // Restrict input field to cell
+ const real_section_type::value_type* fieldInCell =
+ fieldIn->restrict(patch, *cellIter);
+
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double* basis = _quadrature->basis();
+ const double* quadPts = _quadrature->quadPts();
+ const double* quadWts = _quadrature->quadWts();
+ const double* jacobianDet = _quadrature->jacobianDet();
+ const int numBasis = _quadrature->numCorners();
+ const int spaceDim = _quadrature->spaceDim();
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad]*jacobianDet[iQuad];
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iC = iBasis*spaceDim;
+ double val = wt*basis[iQ+iBasis];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ val *= basis[iQ+jBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _cellVector[iC+iDim] += val*fieldInCell[iC+iDim];
+ } // for
+ } // for
+ } // for
+
+ // Assemble cell contribution into field
+ fieldOut->updateAdd(patch, *cellIter, _cellVector);
+ } // for
} // integrateAction
// ----------------------------------------------------------------------
// Compute matrix associated with operator.
void
pylith::feassemble::IntegratorInertia3D::integrate(PetscMat* mat,
- const ALE::Obj<ALE::Mesh::real_section_type>& coordinates)
+ const ALE::Obj<real_section_type>& coordinates)
{ // integrate
} // integrate
More information about the cig-commits
mailing list