[cig-commits] r13143 - short/3D/PyLith/trunk/libsrc/feassemble
knepley at geodynamics.org
knepley at geodynamics.org
Tue Oct 28 20:54:13 PDT 2008
Author: knepley
Date: 2008-10-28 20:54:13 -0700 (Tue, 28 Oct 2008)
New Revision: 13143
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
Log:
I changed some restrictClosure()s and get some speedup
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2008-10-28 21:47:05 UTC (rev 13142)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2008-10-29 03:54:13 UTC (rev 13143)
@@ -195,6 +195,12 @@
totalStrain = 0.0;
PetscLogEventEnd(setupEvent,0,0,0,0);
+ ALE::ISieveVisitor::RestrictVisitor<real_section_type> rV(*dispTBctpdt, cellVecSize, &dispTBctpdtCell[0]);
+ if (mesh->depth() > 1) {
+ //ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::RestrictVisitor<Section> > pV((int) pow((double) mesh->getSieve()->getMaxConeSize(), this->depth())+1, rV, true);
+ throw ALE::Exception("Need to reorganize to use a different visitor class");
+ }
+
// Loop over cells
int c_index = 0;
for (Mesh::label_sequence::iterator c_iter=cells->begin();
@@ -215,7 +221,7 @@
// Restrict input fields to cell
PetscLogEventBegin(restrictEvent,0,0,0,0);
- mesh->restrictClosure(dispTBctpdt, *c_iter, &dispTBctpdtCell[0], cellVecSize);
+ mesh->restrictClosure(*c_iter, rV);
PetscLogEventEnd(restrictEvent,0,0,0,0);
// Get cell geometry information that depends on cell
@@ -278,6 +284,7 @@
PetscLogEventBegin(updateEvent,0,0,0,0);
mesh->updateAdd(residual, *c_iter, _cellVector);
PetscLogEventEnd(updateEvent,0,0,0,0);
+ rV.clear();
} // for
} // integrateResidual
@@ -373,6 +380,7 @@
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
visitor_type iV(*dispTBctpdt, *globalOrder, (int) pow(mesh->getSieve()->getMaxConeSize(), mesh->depth())*spaceDim);
+ ALE::ISieveVisitor::RestrictVisitor<real_section_type> rV(*dispTBctpdt, cellVecSize, &dispTBctpdtCell[0]);
// Loop over cells
int c_index = 0;
@@ -389,7 +397,7 @@
_resetCellMatrix();
// Restrict input fields to cell
- mesh->restrictClosure(dispTBctpdt, *c_iter, &dispTBctpdtCell[0], cellVecSize);
+ mesh->restrictClosure(*c_iter, rV);
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -426,7 +434,7 @@
&lwork, &lierr);
if (lierr)
throw std::runtime_error("Lapack SVD failed");
- minSV = svalues[n-1];
+ minSV = svalues[n-7];
maxSV = svalues[0];
std::cout << "Element " << c_index << std::endl;
for(int i = 0; i < n; ++i)
@@ -443,6 +451,7 @@
if (err)
throw std::runtime_error("Update to PETSc Mat failed.");
iV.clear();
+ rV.clear();
} // for
_needNewJacobian = false;
_material->resetNeedNewJacobian();
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2008-10-28 21:47:05 UTC (rev 13142)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2008-10-29 03:54:13 UTC (rev 13143)
@@ -285,22 +285,27 @@
*std::max_element(cells->begin(), cells->end())+1));
_quadPtsPre->setFiberDimension(cells, _numQuadPts*_spaceDim);
_quadPtsPre->allocatePoint();
+ _quadPtsPreV = new ALE::ISieveVisitor::RestrictVisitor<real_section_type>(*_quadPtsPre, _numQuadPts*_spaceDim, &_quadPts[0]);
_jacobianPre->getAtlas()->setAtlas(_quadPtsPre->getAtlas()->getAtlas());
_jacobianPre->getAtlas()->allocatePoint();
_jacobianPre->setFiberDimension(cells, _numQuadPts*_cellDim*_spaceDim);
_jacobianPre->allocatePoint();
+ _jacobianPreV = new ALE::ISieveVisitor::RestrictVisitor<real_section_type>(*_jacobianPre, _numQuadPts*_cellDim*_spaceDim, &_jacobian[0]);
_jacobianDetPre->getAtlas()->setAtlas(_quadPtsPre->getAtlas()->getAtlas());
_jacobianDetPre->getAtlas()->allocatePoint();
_jacobianDetPre->setFiberDimension(cells, _numQuadPts);
_jacobianDetPre->allocatePoint();
+ _jacobianDetPreV = new ALE::ISieveVisitor::RestrictVisitor<real_section_type>(*_jacobianDetPre, _numQuadPts, &_jacobianDet[0]);
_jacobianInvPre->setAtlas(_jacobianPre->getAtlas());
_jacobianInvPre->setFiberDimension(cells, _numQuadPts*_cellDim*_spaceDim);
_jacobianInvPre->allocatePoint();
+ _jacobianInvPreV = new ALE::ISieveVisitor::RestrictVisitor<real_section_type>(*_jacobianInvPre, _numQuadPts*_cellDim*_spaceDim, &_jacobianInv[0]);
//_jITag = _jacobianInvPre->copyCustomAtlas(_jacobianPre, _jTag);
_basisDerivPre->getAtlas()->setAtlas(_quadPtsPre->getAtlas()->getAtlas());
_basisDerivPre->getAtlas()->allocatePoint();
_basisDerivPre->setFiberDimension(cells, _numQuadPts*_numBasis*_spaceDim);
_basisDerivPre->allocatePoint();
+ _basisDerivPreV = new ALE::ISieveVisitor::RestrictVisitor<real_section_type>(*_basisDerivPre, _numQuadPts*_numBasis*_spaceDim, &_basisDeriv[0]);
const int ncells = cells->size();
const int nbytes =
@@ -349,31 +354,20 @@
const Mesh::point_type& cell,
const int c)
{ // retrieveGeometry
- const real_section_type::value_type* values =
- mesh->restrictClosure(_quadPtsPre, cell);
- int size = _numQuadPts * _spaceDim;
- assert(size == _quadPtsPre->getFiberDimension(cell));
- memcpy(&_quadPts[0], &values[0], size*sizeof(double));
+ _quadPtsPreV->clear();
+ mesh->restrictClosure(cell, *_quadPtsPreV);
- values = mesh->restrictClosure(_jacobianPre, cell);
- size = _numQuadPts * _cellDim * _spaceDim;
- assert(size == _jacobianPre->getFiberDimension(cell));
- memcpy(&_jacobian[0], &values[0], size*sizeof(double));
+ _jacobianPreV->clear();
+ mesh->restrictClosure(cell, *_jacobianPreV);
- values = mesh->restrictClosure(_jacobianDetPre, cell);
- size = _numQuadPts;
- assert(size == _jacobianDetPre->getFiberDimension(cell));
- memcpy(&_jacobianDet[0], &values[0], size*sizeof(double));
+ _jacobianDetPreV->clear();
+ mesh->restrictClosure(cell, *_jacobianDetPreV);
- values = mesh->restrictClosure(_jacobianInvPre, cell);
- size = _numQuadPts * _cellDim * _spaceDim;
- assert(size == _jacobianInvPre->getFiberDimension(cell));
- memcpy(&_jacobianInv[0], &values[0], size*sizeof(double));
+ _jacobianInvPreV->clear();
+ mesh->restrictClosure(cell, *_jacobianInvPreV);
- values = mesh->restrictClosure(_basisDerivPre, cell);
- size = _numQuadPts * _numBasis * _spaceDim;
- assert(size == _basisDerivPre->getFiberDimension(cell));
- memcpy(&_basisDeriv[0], &values[0], size*sizeof(double));
+ _basisDerivPreV->clear();
+ mesh->restrictClosure(cell, *_basisDerivPreV);
} // retrieveGeometry
// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2008-10-28 21:47:05 UTC (rev 13142)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2008-10-29 03:54:13 UTC (rev 13143)
@@ -382,10 +382,15 @@
/* Precomputation sections */
int _qTag, _jTag, _jDTag, _jITag, _bTag;
Obj<real_section_type> _quadPtsPre;
+ Obj<ALE::ISieveVisitor::RestrictVisitor<real_section_type> > _quadPtsPreV;
Obj<real_section_type> _jacobianPre;
+ Obj<ALE::ISieveVisitor::RestrictVisitor<real_section_type> > _jacobianPreV;
Obj<real_section_type> _jacobianDetPre;
+ Obj<ALE::ISieveVisitor::RestrictVisitor<real_section_type> > _jacobianDetPreV;
Obj<real_section_type> _jacobianInvPre;
+ Obj<ALE::ISieveVisitor::RestrictVisitor<real_section_type> > _jacobianInvPreV;
Obj<real_section_type> _basisDerivPre;
+ Obj<ALE::ISieveVisitor::RestrictVisitor<real_section_type> > _basisDerivPreV;
bool _precomputed; ///< True if we have computed geometry info
bool _checkConditioning; ///< True if checking for ill-conditioning
More information about the CIG-COMMITS
mailing list