[cig-commits] r8180 - short/3D/PyLith/trunk/libsrc/feassemble
knepley at geodynamics.org
knepley at geodynamics.org
Fri Oct 26 15:06:29 PDT 2007
Author: knepley
Date: 2007-10-26 15:06:28 -0700 (Fri, 26 Oct 2007)
New Revision: 8180
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
Log:
Added optimization for restrict/update
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-10-26 21:16:08 UTC (rev 8179)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-10-26 22:06:28 UTC (rev 8180)
@@ -176,14 +176,16 @@
const int spaceDim = _quadrature->spaceDim();
#ifdef FASTER
- static std::map<int, int> tags;
- int c = 0;
+ int c = 0;
- if (tags.find(_material->id()) == tags.end()) {
- tags[_material->id()] = mesh->calculateCustomAtlas(dispTBctpdt, cells);
- residual->copyCustomAtlas(dispTBctpdt, tags[_material->id()]);
+ if (_dTags.find(_material->id()) == _dTags.end()) {
+ _dTags[_material->id()] = mesh->calculateCustomAtlas(dispTBctpdt, cells);
}
- const int tag = tags[_material->id()];
+ if (_rTags.find(_material->id()) == _rTags.end()) {
+ _rTags[_material->id()] = residual->copyCustomAtlas(dispTBctpdt, _dTags[_material->id()]);
+ }
+ const int dTag = _dTags[_material->id()];
+ const int rTag = _rTags[_material->id()];
#endif
// Precompute the geometric and function space information
_quadrature->precomputeGeometry(mesh, coordinates, cells);
@@ -208,7 +210,7 @@
++c_iter) {
// Compute geometry information for current cell
PetscLogEventBegin(cellGeomEvent,0,0,0,0);
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter);
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c);
PetscLogEventEnd(cellGeomEvent,0,0,0,0);
// Get state variables for cell.
@@ -222,7 +224,7 @@
// Restrict input fields to cell
PetscLogEventBegin(restrictEvent,0,0,0,0);
#ifdef FASTER
- mesh->restrict(dispTBctpdt, tag, c, &dispTBctpdtCell[0], cellVecSize);
+ mesh->restrict(dispTBctpdt, dTag, c, &dispTBctpdtCell[0], cellVecSize);
#else
mesh->restrict(dispTBctpdt, *c_iter, &dispTBctpdtCell[0], cellVecSize);
#endif
@@ -280,7 +282,7 @@
// Assemble cell contribution into field
PetscLogEventBegin(updateEvent,0,0,0,0);
#ifdef FASTER
- mesh->updateAdd(residual, tag, c++, _cellVector);
+ mesh->updateAdd(residual, rTag, c++, _cellVector);
#else
mesh->updateAdd(residual, *c_iter, _cellVector);
#endif
@@ -375,12 +377,20 @@
totalStrain[iQuad] = 0.0;
} // for
+#ifdef FASTER
+ if (_dTags.find(_material->id()) == _dTags.end()) {
+ _dTags[_material->id()] = mesh->calculateCustomAtlas(dispTBctpdt, cells);
+ }
+ const int dTag = _dTags[_material->id()];
+#endif
+
// Loop over cells
+ int c = 0;
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
- ++c_iter) {
+ ++c_iter, ++c) {
// Compute geometry information for current cell
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter);
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c);
// Get state variables for cell.
_material->getStateVarsCell(*c_iter, numQuadPts);
@@ -389,7 +399,11 @@
_resetCellMatrix();
// Restrict input fields to cell
+#ifdef FASTER
+ mesh->restrict(dispTBctpdt, dTag, c, &dispTBctpdtCell[0], cellVecSize);
+#else
mesh->restrict(dispTBctpdt, *c_iter, &dispTBctpdtCell[0], cellVecSize);
+#endif
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -479,14 +493,23 @@
totalStrain[iQuad] = 0.0;
} // for
+#ifdef FASTER
+ const int dTag = _dTags[_material->id()];
+#endif
+
+ int c = 0;
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
- ++c_iter) {
+ ++c_iter, ++c) {
// Compute geometry information for current cell
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter);
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c);
// Restrict input fields to cell
+#ifdef FASTER
+ mesh->restrict(disp, dTag, c, &dispCell[0], cellVecSize);
+#else
mesh->restrict(disp, *c_iter, &dispCell[0], cellVecSize);
+#endif
// Get cell geometry information that depends on cell
const double_array& basisDeriv = _quadrature->basisDeriv();
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh 2007-10-26 21:16:08 UTC (rev 8179)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh 2007-10-26 22:06:28 UTC (rev 8180)
@@ -212,6 +212,10 @@
/// Elastic material associated with integrator
materials::ElasticMaterial* _material;
+ /// Optimization
+ std::map<int, int> _dTags; // dispTBctpdt tags
+ std::map<int, int> _rTags; // residual tags
+
}; // ElasticityImplicit
#endif // pylith_feassemble_elasticityimplicit_hh
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2007-10-26 21:16:08 UTC (rev 8179)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2007-10-26 22:06:28 UTC (rev 8180)
@@ -281,18 +281,23 @@
_quadPtsPre->setFiberDimension(cells, _numQuadPts*_spaceDim);
_quadPtsPre->allocatePoint();
+ _qTag = mesh->calculateCustomAtlas(_quadPtsPre, cells);
_jacobianPre->getAtlas()->setAtlas(_quadPtsPre->getAtlas()->getAtlas());
_jacobianPre->setFiberDimension(cells, _numQuadPts*_cellDim*_spaceDim);
_jacobianPre->allocatePoint();
+ _jTag = mesh->calculateCustomAtlas(_jacobianPre, cells);
_jacobianDetPre->getAtlas()->setAtlas(_quadPtsPre->getAtlas()->getAtlas());
_jacobianDetPre->setFiberDimension(cells, _numQuadPts);
_jacobianDetPre->allocatePoint();
+ _jDTag = mesh->calculateCustomAtlas(_jacobianDetPre, cells);
_jacobianInvPre->setAtlas(_jacobianPre->getAtlas());
_jacobianInvPre->setFiberDimension(cells, _numQuadPts*_cellDim*_spaceDim);
_jacobianInvPre->allocatePoint();
+ _jITag = _jacobianInvPre->copyCustomAtlas(_jacobianPre, _jTag);
_basisDerivPre->getAtlas()->setAtlas(_quadPtsPre->getAtlas()->getAtlas());
_basisDerivPre->setFiberDimension(cells, _numQuadPts*_numBasis*_spaceDim);
_basisDerivPre->allocatePoint();
+ _bTag = mesh->calculateCustomAtlas(_basisDerivPre, cells);
for(Mesh::label_sequence::iterator c_iter = cells->begin();
c_iter != end;
@@ -321,36 +326,58 @@
void
pylith::feassemble::Quadrature::retrieveGeometry(
const ALE::Obj<Mesh>& mesh,
- const ALE::Obj<real_section_type>& coordinates,
- const Mesh::point_type& cell)
+ const ALE::Obj<real_section_type>& coordinates,
+ const Mesh::point_type& cell,
+ const int c)
{ // retrieveGeometry
- // Do they have a fast copy?
+#define FASTER
+#ifdef FASTER
const real_section_type::value_type* values =
+ mesh->restrict(_quadPtsPre, _qTag, c);
+#else
+ const real_section_type::value_type* values =
_quadPtsPre->restrictPoint(cell);
+#endif
int size = _numQuadPts * _spaceDim;
assert(size == _quadPtsPre->getFiberDimension(cell));
for(int i=0; i < size; ++i)
_quadPts[i] = values[i];
-
+
+#ifdef FASTER
+ values = mesh->restrict(_jacobianPre, _jTag, c);
+#else
values = _jacobianPre->restrictPoint(cell);
+#endif
size = _numQuadPts * _cellDim * _spaceDim;
assert(size == _jacobianPre->getFiberDimension(cell));
for(int i=0; i < size; ++i)
_jacobian[i] = values[i];
+#ifdef FASTER
+ values = mesh->restrict(_jacobianDetPre, _jDTag, c);
+#else
values = _jacobianDetPre->restrictPoint(cell);
+#endif
size = _numQuadPts;
assert(size == _jacobianDetPre->getFiberDimension(cell));
for(int i=0; i < size; ++i)
_jacobianDet[i] = values[i];
+#ifdef FASTER
+ values = mesh->restrict(_jacobianInvPre, _jITag, c);
+#else
values = _jacobianInvPre->restrictPoint(cell);
+#endif
size = _numQuadPts * _cellDim * _spaceDim;
assert(size == _jacobianInvPre->getFiberDimension(cell));
for(int i=0; i < size; ++i)
_jacobianInv[i] = values[i];
+#ifdef FASTER
+ values = mesh->restrict(_basisDerivPre, _bTag, c);
+#else
values = _basisDerivPre->restrictPoint(cell);
+#endif
size = _numQuadPts * _numBasis * _spaceDim;
assert(size == _basisDerivPre->getFiberDimension(cell));
for(int i=0; i < size; ++i)
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2007-10-26 21:16:08 UTC (rev 8179)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2007-10-26 22:06:28 UTC (rev 8180)
@@ -233,7 +233,8 @@
*/
void retrieveGeometry(const ALE::Obj<Mesh>& mesh,
const ALE::Obj<real_section_type>& coordinates,
- const Mesh::point_type& cell);
+ const Mesh::point_type& cell,
+ const int c);
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
@@ -368,6 +369,7 @@
bool _precomputed;
/* Precomputation sections */
+ int _qTag, _jTag, _jDTag, _jITag, _bTag;
Obj<real_section_type> _quadPtsPre;
Obj<real_section_type> _jacobianPre;
Obj<real_section_type> _jacobianDetPre;
More information about the cig-commits
mailing list