[cig-commits] r22566 - in short/3D/PyLith/trunk/libsrc/pylith: feassemble materials topology
brad at geodynamics.org
brad at geodynamics.org
Wed Jul 10 15:33:04 PDT 2013
Author: brad
Date: 2013-07-10 15:33:03 -0700 (Wed, 10 Jul 2013)
New Revision: 22566
Modified:
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/CoordsVisitor.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/CoordsVisitor.icc
short/3D/PyLith/trunk/libsrc/pylith/topology/VisitorMesh.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/VisitorMesh.icc
Log:
Additional use of optimized getClosure().
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc 2013-07-10 22:33:03 UTC (rev 22566)
@@ -160,6 +160,7 @@
const scalar_array& quadWts = _quadrature->quadWts();
assert(quadWts.size() == numQuadPts);
const int numBasis = _quadrature->numBasis();
+ const int numCorners = _quadrature->refGeometry().numCorners();
const int spaceDim = _quadrature->spaceDim();
const int cellDim = _quadrature->cellDim();
const int tensorSize = _material->tensorSize();
@@ -211,15 +212,19 @@
const PetscInt numCells = _materialIS->size();
// Setup field visitors.
+ scalar_array accCell(numBasis*spaceDim);
topology::VecVisitorMesh accVisitor(fields->get("acceleration(t)"));
+ scalar_array velCell(numBasis*spaceDim);
topology::VecVisitorMesh velVisitor(fields->get("velocity(t)"));
+ scalar_array dispCell(numBasis*spaceDim);
scalar_array dispAdjCell(numBasis*spaceDim);
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
topology::VecVisitorMesh residualVisitor(residual);
+ scalar_array coordsCell(numBasis*spaceDim); // :KULDGE: Update numBasis to numCorners after implementing higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
assert(_normalizer);
@@ -237,11 +242,6 @@
_logger->eventBegin(computeEvent);
#endif
- PetscScalar *coordsCell = new PetscScalar[64];
- PetscScalar *accCell = new PetscScalar[64];
- PetscScalar *velCell = new PetscScalar[64];
- PetscScalar *dispCell = new PetscScalar[64];
-
// Loop over cells
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
@@ -249,10 +249,8 @@
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
- PetscInt coordsSize = 64;
- coordsVisitor.getClosure((PetscScalar **) &coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- _quadrature->computeGeometry(coordsCell, coordsSize, cell);
- //coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ _quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -271,15 +269,10 @@
_resetCellVector();
// Restrict input fields to cell
- PetscInt accSize = 64;
- accVisitor.getClosure((PetscScalar **) &accCell, &accSize, cell);assert(accCell);assert(numBasis*spaceDim == accSize);
+ accVisitor.getClosure(&accCell, cell);
+ velVisitor.getClosure(&velCell, cell);
+ dispVisitor.getClosure(&dispCell, cell);
- PetscInt velSize = 64;
- velVisitor.getClosure((PetscScalar **) &velCell, &velSize, cell);assert(velCell);assert(numBasis*spaceDim == velSize);
-
- PetscInt dispSize = 64;
- dispVisitor.getClosure((PetscScalar **) &dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -343,12 +336,9 @@
// Numerical damping. Compute displacements adjusted by velocity
// times normalized viscosity.
- for(PetscInt i = 0; i < dispSize; ++i) {
+ for(PetscInt i = 0, dispSize = dispCell.size(); i < dispSize; ++i) {
dispAdjCell[i] = dispCell[i] + viscosity * velCell[i];
} // for
- //accVisitor.restoreClosure(&accCell, &accSize, cell);
- //velVisitor.restoreClosure(&velCell, &velSize, cell);
- //dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numQuadPts*(4+numBasis*3));
@@ -380,11 +370,6 @@
#endif
} // for
- delete [] coordsCell;
- delete [] accCell;
- delete [] velCell;
- delete [] dispCell;
-
#if !defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numCells*numQuadPts*(4+numBasis*3));
_logger->eventEnd(computeEvent);
@@ -435,6 +420,7 @@
const scalar_array& quadWts = _quadrature->quadWts();
assert(quadWts.size() == numQuadPts);
const int numBasis = _quadrature->numBasis();
+ const int numCorners = _quadrature->refGeometry().numCorners();
const int spaceDim = _quadrature->spaceDim();
const int cellDim = _quadrature->cellDim();
if (cellDim != spaceDim)
@@ -460,6 +446,7 @@
PetscScalar* jacobianCell = NULL;
PetscInt jacobianSize = 0;
+ scalar_array coordsCell(numBasis*spaceDim); // :KLUDGE: numBasis to numCorners after switching to higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
_logger->eventEnd(setupEvent);
@@ -474,11 +461,8 @@
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
- PetscScalar* coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- _quadrature->computeGeometry(coordsCell, coordsSize, cell);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ _quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc 2013-07-10 22:33:03 UTC (rev 22566)
@@ -155,6 +155,7 @@
const scalar_array& quadWts = _quadrature->quadWts();
assert(quadWts.size() == numQuadPts);
const int numBasis = _quadrature->numBasis();
+ const int numCorners = _quadrature->refGeometry().numCorners();
const int spaceDim = _quadrature->spaceDim();
const int cellDim = _quadrature->cellDim();
const int tensorSize = _material->tensorSize();
@@ -204,11 +205,19 @@
const PetscInt numCells = _materialIS->size();
// Setup field visitors.
- scalar_array dispAdjCell(numBasis*spaceDim);
+ scalar_array accCell(numBasis*spaceDim);
topology::VecVisitorMesh accVisitor(fields->get("acceleration(t)"));
+
+ scalar_array velCell(numBasis*spaceDim);
topology::VecVisitorMesh velVisitor(fields->get("velocity(t)"));
+
+ scalar_array dispCell(numBasis*spaceDim);
+ scalar_array dispAdjCell(numBasis*spaceDim);
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
+
topology::VecVisitorMesh residualVisitor(residual);
+
+ scalar_array coordsCell(numBasis*spaceDim); // :KULDGE: Update numBasis to numCorners after implementing higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
assert(_normalizer);
@@ -229,10 +238,8 @@
const PetscInt cell = cells[c];
// Compute geometry information for current cell
- PetscScalar *coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- _quadrature->computeGeometry(coordsCell, coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ _quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
// Get state variables for cell.
_material->retrievePropsAndVars(cell);
@@ -241,18 +248,10 @@
_resetCellVector();
// Restrict input fields to cell
- PetscScalar* accCell = NULL;
- PetscInt accSize = 0;
- accVisitor.getClosure(&accCell, &accSize, cell);assert(accCell);assert(numBasis*spaceDim == accSize);
+ accVisitor.getClosure(&accCell, cell);
+ velVisitor.getClosure(&velCell, cell);
+ dispVisitor.getClosure(&dispCell, cell);
- PetscScalar* velCell = NULL;
- PetscInt velSize = 0;
- velVisitor.getClosure(&velCell, &velSize, cell);assert(velCell);assert(numBasis*spaceDim == velSize);
-
- PetscScalar* dispCell = NULL;
- PetscInt dispSize = 0;
- dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
-
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
const scalar_array& basisDeriv = _quadrature->basisDeriv();
@@ -312,15 +311,12 @@
// Numerical damping. Compute displacements adjusted by velocity
// times normalized viscosity.
- for (PetscInt i = 0; i < dispSize; ++i) {
+ for (PetscInt i = 0, dispSize = dispCell.size(); i < dispSize; ++i) {
dispAdjCell[i] = dispCell[i] + viscosity * velCell[i];
} // for
- accVisitor.restoreClosure(&accCell, &accSize, cell);
- velVisitor.restoreClosure(&velCell, &velSize, cell);
- dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
// Compute B(transpose) * sigma, first computing strains
- _calcDeformation(&deformCell, basisDeriv, coordsCell, &dispAdjCell[0], numBasis, numQuadPts, spaceDim);
+ _calcDeformation(&deformCell, basisDeriv, &coordsCell[0], &dispAdjCell[0], numBasis, numQuadPts, spaceDim);
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
const scalar_array& stressCell = _material->calcStress(strainCell, true);
@@ -328,8 +324,6 @@
// Assemble cell contribution into field
residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), cell, ADD_VALUES);
-
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
} // for
_logger->eventEnd(computeEvent);
@@ -379,6 +373,7 @@
const scalar_array& quadWts = _quadrature->quadWts();
assert(quadWts.size() == numQuadPts);
const int numBasis = _quadrature->numBasis();
+ const int numCorners = _quadrature->refGeometry().numCorners();
const int spaceDim = _quadrature->spaceDim();
const int cellDim = _quadrature->cellDim();
const int tensorSize = _material->tensorSize();
@@ -402,6 +397,7 @@
scalar_array valuesIJ(numBasis);
topology::VecVisitorMesh jacobianVisitor(*jacobian);
+ scalar_array coordsCell(numBasis*spaceDim); // :KULDGE: Update numBasis to numCorners after implementing higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
_logger->eventEnd(setupEvent);
@@ -411,11 +407,8 @@
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
// Compute geometry information for current cell
- PetscScalar* coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- _quadrature->computeGeometry(coordsCell, coordsSize, cell);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ _quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
// Get state variables for cell.
_material->retrievePropsAndVars(cell);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc 2013-07-10 22:33:03 UTC (rev 22566)
@@ -54,6 +54,7 @@
const int pylith::feassemble::ElasticityExplicitTet4::_cellDim = 3;
const int pylith::feassemble::ElasticityExplicitTet4::_tensorSize = 6;
const int pylith::feassemble::ElasticityExplicitTet4::_numBasis = 4;
+const int pylith::feassemble::ElasticityExplicitTet4::_numCorners = 4;
const int pylith::feassemble::ElasticityExplicitTet4::_numQuadPts = 1;
// ----------------------------------------------------------------------
@@ -170,6 +171,7 @@
const int cellDim = _cellDim;
const int tensorSize = _tensorSize;
const int numBasis = _numBasis;
+ const int numCorners = _numCorners;
const int numQuadPts = _numQuadPts;
const int cellVectorSize = numBasis*spaceDim;
/** :TODO:
@@ -197,15 +199,19 @@
const PetscInt numCells = _materialIS->size();
// Setup field visitors.
+ scalar_array accCell(numBasis*spaceDim);
topology::VecVisitorMesh accVisitor(fields->get("acceleration(t)"));
+ scalar_array velCell(numBasis*spaceDim);
topology::VecVisitorMesh velVisitor(fields->get("velocity(t)"));
+ scalar_array dispCell(numBasis*spaceDim);
scalar_array dispAdjCell(numBasis*spaceDim);
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
topology::VecVisitorMesh residualVisitor(residual);
+ scalar_array coordsCell(numCorners*spaceDim);
topology::CoordsVisitor coordsVisitor(dmMesh);
assert(_normalizer);
@@ -231,17 +237,9 @@
#endif
// Restrict input fields to cell
- PetscScalar* accCell = NULL;
- PetscInt accSize = 0;
- accVisitor.getClosure(&accCell, &accSize, cell);assert(accCell);assert(numBasis*spaceDim == accSize);
-
- PetscScalar* velCell = NULL;
- PetscInt velSize = 0;
- velVisitor.getClosure(&velCell, &velSize, cell);assert(velCell);assert(numBasis*spaceDim == velSize);
-
- PetscScalar* dispCell = NULL;
- PetscInt dispSize = 0;
- dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
+ accVisitor.getClosure(&accCell, cell);
+ velVisitor.getClosure(&velCell, cell);
+ dispVisitor.getClosure(&dispCell, cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
@@ -249,10 +247,8 @@
#endif
// Compute geometry information for current cell
- PetscScalar *coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- const PylithScalar volume = _volume(coordsCell, coordsSize);assert(volume > 0.0);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ const PylithScalar volume = _volume(coordsCell);assert(volume > 0.0);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -303,7 +299,7 @@
// Compute action for inertial terms
const PylithScalar wtVertex = density[0] * volume / 4.0;
- assert(cellVectorSize == dispSize);
+ assert(cellVectorSize == dispCell.size());
for(PetscInt i = 0; i < cellVectorSize; ++i) {
_cellVector[i] -= wtVertex * accCell[i];
} // for
@@ -319,9 +315,6 @@
for(PetscInt i = 0; i < cellVectorSize; ++i) {
dispAdjCell[i] = dispCell[i] + viscosity * velCell[i];
} // for
- accVisitor.restoreClosure(&accCell, &accSize, cell);
- velVisitor.restoreClosure(&velCell, &velSize, cell);
- dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
// Compute B(transpose) * sigma, first computing strains
const PylithScalar x0 = coordsCell[0];
@@ -412,8 +405,6 @@
_logger->eventBegin(updateEvent);
#endif
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
-
// Assemble cell contribution into field
residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), cell, ADD_VALUES);
@@ -475,6 +466,7 @@
const int spaceDim = _spaceDim;
const int cellDim = _cellDim;
const int numBasis = _numBasis;
+ const int numCorners = _numCorners;
if (cellDim != spaceDim)
throw std::logic_error("Don't know how to integrate elasticity " \
"contribution to Jacobian matrix for cells with " \
@@ -496,6 +488,7 @@
PetscScalar* jacobianCell = NULL;
PetscInt jacobianSize = 0;
+ scalar_array coordsCell(numCorners*spaceDim);
topology::CoordsVisitor coordsVisitor(dmMesh);
_logger->eventEnd(setupEvent);
@@ -509,11 +502,8 @@
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
- PetscScalar* coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- const PylithScalar volume = _volume(coordsCell, coordsSize);assert(volume > 0.0);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ const PylithScalar volume = _volume(coordsCell);assert(volume > 0.0);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -584,10 +574,9 @@
// ----------------------------------------------------------------------
// Compute volume of tetrahedral cell.
PylithScalar
-pylith::feassemble::ElasticityExplicitTet4::_volume(const PylithScalar coordinatesCell[],
- const int coordinatesSize) const
+pylith::feassemble::ElasticityExplicitTet4::_volume(const scalar_array& coordinatesCell) const
{ // __volume
- assert(12 == coordinatesSize);
+ assert(12 == coordinatesCell.size());
const PylithScalar x0 = coordinatesCell[0];
const PylithScalar y0 = coordinatesCell[1];
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.hh 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.hh 2013-07-10 22:33:03 UTC (rev 22566)
@@ -138,11 +138,9 @@
/** Compute volume of tetrahedral cell.
*
* @param coordinatesCell Coordinates of vertices of cell.
- * @param coordinatesSize Size of array.
* @returns Volume of cell.
*/
- PylithScalar _volume(const PylithScalar coordinatesCell[],
- const int coordinatesSize) const;
+ PylithScalar _volume(const scalar_array& coordinatesCell) const;
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
@@ -154,6 +152,7 @@
static const int _cellDim;
static const int _tensorSize;
static const int _numBasis;
+ static const int _numCorners;
static const int _numQuadPts;
// NOT IMPLEMENTED //////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc 2013-07-10 22:33:03 UTC (rev 22566)
@@ -54,6 +54,7 @@
const int pylith::feassemble::ElasticityExplicitTri3::_cellDim = 2;
const int pylith::feassemble::ElasticityExplicitTri3::_tensorSize = 3;
const int pylith::feassemble::ElasticityExplicitTri3::_numBasis = 3;
+const int pylith::feassemble::ElasticityExplicitTri3::_numCorners = 3;
const int pylith::feassemble::ElasticityExplicitTri3::_numQuadPts = 1;
// ----------------------------------------------------------------------
@@ -96,7 +97,7 @@
_dtm1 = dt;
_dt = dt;
assert(_dt == _dtm1); // For now, don't allow variable time step
- if (0 != _material)
+ if (_material)
_material->timeStep(_dt);
PYLITH_METHOD_END;
@@ -169,6 +170,7 @@
const int cellDim = _cellDim;
const int tensorSize = _tensorSize;
const int numBasis = _numBasis;
+ const int numCorners = _numCorners;
const int numQuadPts = _numQuadPts;
const int cellVectorSize = numBasis*spaceDim;
/** :TODO:
@@ -196,15 +198,19 @@
const PetscInt numCells = _materialIS->size();
// Setup field visitors.
+ scalar_array accCell(numBasis*spaceDim);
topology::VecVisitorMesh accVisitor(fields->get("acceleration(t)"));
+ scalar_array velCell(numBasis*spaceDim);
topology::VecVisitorMesh velVisitor(fields->get("velocity(t)"));
+ scalar_array dispCell(numBasis*spaceDim);
scalar_array dispAdjCell(numBasis*spaceDim);
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
topology::VecVisitorMesh residualVisitor(residual);
+ scalar_array coordsCell(numCorners*spaceDim);
topology::CoordsVisitor coordsVisitor(dmMesh);
assert(_normalizer);
@@ -222,11 +228,6 @@
_logger->eventBegin(computeEvent);
#endif
- PetscScalar *coordsCell = new PetscScalar[64];
- PetscScalar *accCell = new PetscScalar[64];
- PetscScalar *velCell = new PetscScalar[64];
- PetscScalar *dispCell = new PetscScalar[64];
-
// Loop over cells
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
@@ -235,24 +236,18 @@
#endif
// Restrict input fields to cell
- PetscInt accSize = 64;
- accVisitor.getClosure(&accCell, &accSize, cell);assert(accCell);assert(numBasis*spaceDim == accSize);
+ accVisitor.getClosure(&accCell, cell);
+ velVisitor.getClosure(&velCell, cell);
+ dispVisitor.getClosure(&dispCell, cell);
- PetscInt velSize = 64;
- velVisitor.getClosure(&velCell, &velSize, cell);assert(velCell);assert(numBasis*spaceDim == velSize);
-
- PetscInt dispSize = 64;
- dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(geometryEvent);
#endif
// Compute geometry information for current cell
- PetscInt coordsSize = 64;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- const PylithScalar area = _area(coordsCell, coordsSize);assert(area > 0.0);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ const PylithScalar area = _area(coordsCell);assert(area > 0.0);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -303,7 +298,7 @@
// Compute action for inertial terms
const PylithScalar wtVertex = density[0] * area / 3.0;
- assert(cellVectorSize == dispSize);
+ assert(cellVectorSize == dispCell.size());
for(PetscInt i = 0; i < cellVectorSize; ++i) {
_cellVector[i] -= wtVertex * accCell[i];
} // for
@@ -319,9 +314,6 @@
for(PetscInt i = 0; i < cellVectorSize; ++i) {
dispAdjCell[i] = dispCell[i] + viscosity * velCell[i];
} // for
- //accVisitor.restoreClosure(&accCell, &accSize, cell);
- //velVisitor.restoreClosure(&velCell, &velSize, cell);
- //dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
// Compute B(transpose) * sigma, first computing strains
const PylithScalar x0 = coordsCell[0];
@@ -374,8 +366,6 @@
_logger->eventBegin(updateEvent);
#endif
- //coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
-
// Assemble cell contribution into field
residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), cell, ADD_VALUES);
@@ -384,11 +374,6 @@
#endif
} // for
- delete [] coordsCell;
- delete [] accCell;
- delete [] velCell;
- delete [] dispCell;
-
#if !defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numCells*(2 + numBasis*spaceDim*2 + 34+30));
_logger->eventEnd(computeEvent);
@@ -443,6 +428,7 @@
const int spaceDim = _spaceDim;
const int cellDim = _cellDim;
const int numBasis = _numBasis;
+ const int numCorners = _numCorners;
if (cellDim != spaceDim)
throw std::logic_error("Don't know how to integrate elasticity " \
"contribution to Jacobian matrix for cells with " \
@@ -464,6 +450,7 @@
PetscScalar* jacobianCell = NULL;
PetscInt jacobianSize = 0;
+ scalar_array coordsCell(numCorners*spaceDim);
topology::CoordsVisitor coordsVisitor(dmMesh);
_logger->eventEnd(setupEvent);
@@ -477,11 +464,8 @@
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
- PetscScalar* coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- const PylithScalar area = _area(coordsCell, coordsSize);assert(area > 0.0);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ const PylithScalar area = _area(coordsCell);assert(area > 0.0);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -550,10 +534,9 @@
// ----------------------------------------------------------------------
// Compute area of triangular cell.
PylithScalar
-pylith::feassemble::ElasticityExplicitTri3::_area(const PylithScalar coordinatesCell[],
- const int coordinatesSize) const
+pylith::feassemble::ElasticityExplicitTri3::_area(const scalar_array& coordinatesCell) const
{ // __area
- assert(6 == coordinatesSize);
+ assert(6 == coordinatesCell.size());
const PylithScalar x0 = coordinatesCell[0];
const PylithScalar y0 = coordinatesCell[1];
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.hh 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.hh 2013-07-10 22:33:03 UTC (rev 22566)
@@ -138,11 +138,9 @@
/** Compute area of triangular cell.
*
* @param coordinatesCell Coordinates of vertices of cell.
- * @param coordinatesSize Size of array.
- * @returns Volume of cell.
+ * @returns Area of cell.
*/
- PylithScalar _area(const PylithScalar coordinatesCell[],
- const int coordinatesSize) const;
+ PylithScalar _area(const scalar_array& coordinatesCell) const;
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
@@ -154,6 +152,7 @@
static const int _cellDim;
static const int _tensorSize;
static const int _numBasis;
+ static const int _numCorners;
static const int _numQuadPts;
// NOT IMPLEMENTED //////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc 2013-07-10 22:33:03 UTC (rev 22566)
@@ -48,8 +48,6 @@
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
-//#define PRECOMPUTE_GEOMETRY
-
// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::ElasticityImplicit::ElasticityImplicit(void) :
@@ -176,9 +174,15 @@
const PetscInt numCells = _materialIS->size();
// Setup field visitors.
+ scalar_array dispCell(numBasis*spaceDim);
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
+
+ scalar_array dispIncrCell(numBasis*spaceDim);
topology::VecVisitorMesh dispIncrVisitor(fields->get("dispIncr(t->t+dt)"));
+
topology::VecVisitorMesh residualVisitor(residual);
+
+ scalar_array coordsCell(numBasis*spaceDim); // :KLUDGE: numBasis to numCorners after switching to higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
assert(_normalizer);
@@ -192,15 +196,8 @@
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
// Compute geometry information for current cell
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(cell);
-#else
- PetscScalar *coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- _quadrature->computeGeometry(coordsCell, coordsSize, cell);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
-#endif
+ coordsVisitor.getClosure(&coordsCell, cell);
+ _quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
// Get state variables for cell.
_material->retrievePropsAndVars(cell);
@@ -209,14 +206,9 @@
_resetCellVector();
// Restrict input fields to cell
- PetscScalar* dispCell = NULL;
- PetscInt dispSize = 0;
- dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
+ dispVisitor.getClosure(&dispCell, cell);
+ dispIncrVisitor.getClosure(&dispIncrCell, cell);
- PetscScalar* dispIncrCell = NULL;
- PetscInt dispIncrSize = 0;
- dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);assert(dispIncrCell);assert(numBasis*spaceDim == dispIncrSize);
-
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
const scalar_array& basisDeriv = _quadrature->basisDeriv();
@@ -224,11 +216,9 @@
const scalar_array& quadPtsNondim = _quadrature->quadPts();
// Compute current estimate of displacement at time t+dt using solution increment.
- for(PetscInt i = 0; i < dispSize; ++i) {
+ for(PetscInt i = 0, dispSize = dispCell.size(); i < dispSize; ++i) {
dispTpdtCell[i] = dispCell[i] + dispIncrCell[i];
} // for
- dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
- dispIncrVisitor.restoreClosure(&dispIncrCell, &dispIncrSize, cell);
// Compute body force vector if gravity is being used.
if (_gravityField) {
@@ -314,6 +304,7 @@
const scalar_array& quadWts = _quadrature->quadWts();
assert(quadWts.size() == numQuadPts);
const int numBasis = _quadrature->numBasis();
+ const int numCorners = _quadrature->refGeometry().numCorners();
const int spaceDim = _quadrature->spaceDim();
const int cellDim = _quadrature->cellDim();
const int tensorSize = _material->tensorSize();
@@ -357,8 +348,13 @@
const PetscInt numCells = _materialIS->size();
// Setup field visitors.
+ scalar_array dispCell(numBasis*spaceDim);
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
+
+ scalar_array dispIncrCell(numBasis*spaceDim);
topology::VecVisitorMesh dispIncrVisitor(fields->get("dispIncr(t->t+dt)"));
+
+ scalar_array coordsCell(numBasis*spaceDim); // :KLUDGE: numBasis to numCorners after switching to higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
// Get sparse matrix
@@ -375,16 +371,10 @@
// Loop over cells
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
+
// Compute geometry information for current cell
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(cell);
-#else
- PetscScalar* coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- _quadrature->computeGeometry(coordsCell, coordsSize, cell);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
-#endif
+ coordsVisitor.getClosure(&coordsCell, cell);
+ _quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
// Get physical properties and state variables for cell.
_material->retrievePropsAndVars(cell);
@@ -393,25 +383,18 @@
_resetCellMatrix();
// Restrict input fields to cell
- PetscScalar* dispCell = NULL;
- PetscInt dispSize = 0;
- dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
+ dispVisitor.getClosure(&dispCell, cell);
+ dispIncrVisitor.getClosure(&dispIncrCell, cell);
- PetscScalar* dispIncrCell = NULL;
- PetscInt dispIncrSize = 0;
- dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);assert(dispIncrCell);assert(numBasis*spaceDim == dispIncrSize);
-
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
const scalar_array& basisDeriv = _quadrature->basisDeriv();
const scalar_array& jacobianDet = _quadrature->jacobianDet();
// Compute current estimate of displacement at time t+dt using solution increment.
- for(PetscInt i = 0; i < dispSize; ++i) {
+ for(PetscInt i = 0, dispSize = dispCell.size(); i < dispSize; ++i) {
dispTpdtCell[i] = dispCell[i] + dispIncrCell[i];
} // for
- dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
- dispIncrVisitor.restoreClosure(&dispIncrCell, &dispIncrSize, cell);
// Compute strains
calcTotalStrainFn(&strainCell, basisDeriv, &dispTpdtCell[0], numBasis, spaceDim, numQuadPts);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc 2013-07-10 22:33:03 UTC (rev 22566)
@@ -134,6 +134,7 @@
const scalar_array& quadWts = _quadrature->quadWts();
assert(quadWts.size() == numQuadPts);
const int numBasis = _quadrature->numBasis();
+ const int numCorners = _quadrature->refGeometry().numCorners();
const int spaceDim = _quadrature->spaceDim();
const int cellDim = _quadrature->cellDim();
const int tensorSize = _material->tensorSize();
@@ -179,9 +180,15 @@
const PetscInt numCells = _materialIS->size();
// Setup field visitors.
+ scalar_array dispCell(numBasis*spaceDim);
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
+
+ scalar_array dispIncrCell(numBasis*spaceDim);
topology::VecVisitorMesh dispIncrVisitor(fields->get("dispIncr(t->t+dt)"));
+
topology::VecVisitorMesh residualVisitor(residual);
+
+ scalar_array coordsCell(numBasis*spaceDim); // :KULDGE: Update numBasis to numCorners after implementing higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
assert(_normalizer);
@@ -198,10 +205,8 @@
const PetscInt cell = cells[c];
// Compute geometry information for current cell
- PetscScalar *coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- _quadrature->computeGeometry(coordsCell, coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ _quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
// Get state variables for cell.
_material->retrievePropsAndVars(cell);
@@ -210,14 +215,9 @@
_resetCellVector();
// Restrict input fields to cell
- PetscScalar* dispCell = NULL;
- PetscInt dispSize = 0;
- dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
+ dispVisitor.getClosure(&dispCell, cell);
+ dispIncrVisitor.getClosure(&dispIncrCell, cell);
- PetscScalar* dispIncrCell = NULL;
- PetscInt dispIncrSize = 0;
- dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);assert(dispIncrCell);assert(numBasis*spaceDim == dispIncrSize);
-
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
const scalar_array& basisDeriv = _quadrature->basisDeriv();
@@ -226,11 +226,9 @@
// Compute current estimate of displacement at time t+dt using
// solution increment.
- for(PetscInt i = 0; i < dispSize; ++i) {
+ for(PetscInt i = 0, dispSize = dispCell.size(); i < dispSize; ++i) {
dispTpdtCell[i] = dispCell[i] + dispIncrCell[i];
} // for
- dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
- dispIncrVisitor.restoreClosure(&dispIncrCell, &dispIncrSize, cell);
// Compute body force vector if gravity is being used.
if (_gravityField) {
@@ -263,7 +261,7 @@
// Compute B(transpose) * sigma, first computing deformation
// tensor and strains
- _calcDeformation(&deformCell, basisDeriv, coordsCell, &dispTpdtCell[0], numBasis, numQuadPts, spaceDim);
+ _calcDeformation(&deformCell, basisDeriv, &coordsCell[0], &dispTpdtCell[0], numBasis, numQuadPts, spaceDim);
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
const scalar_array& stressCell = _material->calcStress(strainCell, true);
@@ -277,8 +275,6 @@
#endif
// Assemble cell contribution into field
residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), cell, ADD_VALUES);
-
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
} // for
_logger->eventEnd(computeEvent);
@@ -363,8 +359,13 @@
const PetscInt numCells = _materialIS->size();
// Setup field visitors.
+ scalar_array dispCell(numBasis*spaceDim);
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
+
+ scalar_array dispIncrCell(numBasis*spaceDim);
topology::VecVisitorMesh dispIncrVisitor(fields->get("dispIncr(t->t+dt)"));
+
+ scalar_array coordsCell(numBasis*spaceDim); // :KULDGE: Update numBasis to numCorners after implementing higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
// Get sparse matrix
@@ -381,11 +382,10 @@
// Loop over cells
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
+
// Compute geometry information for current cell
- PetscScalar* coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- _quadrature->computeGeometry(coordsCell, coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ _quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
// Get state variables for cell.
_material->retrievePropsAndVars(cell);
@@ -394,14 +394,9 @@
_resetCellMatrix();
// Restrict input fields to cell
- PetscScalar* dispCell = NULL;
- PetscInt dispSize = 0;
- dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
+ dispVisitor.getClosure(&dispCell, cell);
+ dispIncrVisitor.getClosure(&dispIncrCell, cell);
- PetscScalar* dispIncrCell = NULL;
- PetscInt dispIncrSize = 0;
- dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);assert(dispIncrCell);assert(numBasis*spaceDim == dispIncrSize);
-
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
const scalar_array& basisDeriv = _quadrature->basisDeriv();
@@ -409,14 +404,12 @@
// Compute current estimate of displacement at time t+dt using
// solution increment.
- for(PetscInt i = 0; i < dispSize; ++i) {
+ for(PetscInt i = 0, dispSize = dispCell.size(); i < dispSize; ++i) {
dispTpdtCell[i] = dispCell[i] + dispIncrCell[i];
} // for
- dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
- dispIncrVisitor.restoreClosure(&dispIncrCell, &dispIncrSize, cell);
// Compute deformation tensor, strains, and stresses
- _calcDeformation(&deformCell, basisDeriv, coordsCell, &dispTpdtCell[0], numBasis, numQuadPts, spaceDim);
+ _calcDeformation(&deformCell, basisDeriv, &coordsCell[0], &dispTpdtCell[0], numBasis, numQuadPts, spaceDim);
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
// Get "elasticity" matrix at quadrature points for this cell
@@ -463,8 +456,6 @@
// Assemble cell contribution into PETSc matrix.
jacobianVisitor.setClosure(&_cellMatrix[0], _cellMatrix.size(), cell, ADD_VALUES);
-
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
} // for
_logger->eventEnd(computeEvent);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc 2013-07-10 22:33:03 UTC (rev 22566)
@@ -34,7 +34,7 @@
pylith::feassemble::GeometryQuad2D::GeometryQuad2D(void) :
CellGeometry(QUADRILATERAL, 2)
{ // constructor
- const PylithScalar vertices[] = {
+ const PylithScalar vertices[4*2] = {
-1.0, -1.0,
+1.0, -1.0,
+1.0, +1.0,
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc 2013-07-10 22:33:03 UTC (rev 22566)
@@ -172,6 +172,7 @@
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
+ const int numCorners = _quadrature->refGeometry().numCorners();
const int tensorSize = _material->tensorSize();
totalStrain_fn_type calcTotalStrainFn;
if (1 == cellDim) {
@@ -197,35 +198,31 @@
const PetscInt numCells = _materialIS->size();
// Setup visitors.
+ scalar_array dispCell(numBasis*spaceDim);
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
+
+ scalar_array coordsCell(numCorners*spaceDim);
topology::CoordsVisitor coordsVisitor(dmMesh);
// Loop over cells
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
// Retrieve geometry information for current cell
- PetscScalar* coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- _quadrature->computeGeometry(coordsCell, coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ _quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
const scalar_array& basisDeriv = _quadrature->basisDeriv();
// Get physical properties and state variables for cell.
_material->retrievePropsAndVars(cell);
// Restrict input fields to cell
- PetscScalar* dispCell = NULL;
- PetscInt dispSize = 0;
- dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
+ dispVisitor.getClosure(&dispCell, cell);
// Compute strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispCell, numBasis, spaceDim, numQuadPts);
+ calcTotalStrainFn(&strainCell, basisDeriv, &dispCell[0], numBasis, spaceDim, numQuadPts);
// Update material state
_material->updateStateVars(strainCell, cell);
-
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
- dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
} // for
PYLITH_METHOD_END;
@@ -503,6 +500,7 @@
const int cellDim = _quadrature->cellDim();
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
+ const int numCorners = _quadrature->refGeometry().numCorners();
const int spaceDim = _quadrature->spaceDim();
const int tensorSize = _material->tensorSize();
totalStrain_fn_type calcTotalStrainFn;
@@ -533,32 +531,29 @@
const PetscInt numCells = _materialIS->size();
// Setup field visitors.
+ scalar_array dispCell(numBasis*spaceDim);
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
topology::VecVisitorMesh fieldVisitor(*field);
PetscScalar* fieldArray = fieldVisitor.localArray();
+ scalar_array coordsCell(numBasis*spaceDim); // :KULDGE: Update numBasis to numCorners after implementing higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
// Loop over cells
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
+
// Retrieve geometry information for current cell
- PetscScalar *coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- _quadrature->computeGeometry(coordsCell, coordsSize, cell);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ _quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
// Get cell geometry information that depends on cell
- PetscScalar* dispCell = NULL;
- PetscInt dispSize = 0;
- dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
+ dispVisitor.getClosure(&dispCell, cell);
const scalar_array& basisDeriv = _quadrature->basisDeriv();
// Compute strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispCell, numBasis, spaceDim, numQuadPts);
- dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
+ calcTotalStrainFn(&strainCell, basisDeriv, &dispCell[0], numBasis, spaceDim, numQuadPts);
const PetscInt off = fieldVisitor.sectionOffset(cell);
assert(tensorCellSize == fieldVisitor.sectionDof(cell));
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc 2013-07-10 22:33:03 UTC (rev 22566)
@@ -81,6 +81,7 @@
const int cellDim = _quadrature->cellDim();
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
+ const int numCorners = _quadrature->refGeometry().numCorners();
const int spaceDim = _quadrature->spaceDim();
const int tensorSize = _material->tensorSize();
totalStrain_fn_type calcTotalStrainFn;
@@ -109,36 +110,31 @@
const PetscInt* cells = _materialIS->points();
const PetscInt numCells = _materialIS->size();
- scalar_array dispTCell(numBasis*spaceDim);
+ scalar_array dispCell(numBasis*spaceDim);
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
+ scalar_array coordsCell(numBasis*spaceDim); // :KULDGE: Update numBasis to numCorners after implementing higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
// Loop over cells
for (PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
+
// Retrieve geometry information for current cell
- PetscScalar* coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- _quadrature->computeGeometry(coordsCell, coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ _quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
const scalar_array& basisDeriv = _quadrature->basisDeriv();
- PetscScalar* dispCell = NULL;
- PetscInt dispSize = 0;
- dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
+ dispVisitor.getClosure(&dispCell, cell);
// Compute deformation tensor.
- _calcDeformation(&deformCell, basisDeriv, coordsCell, dispCell, numBasis, numQuadPts, spaceDim);
+ _calcDeformation(&deformCell, basisDeriv, &coordsCell[0], &dispCell[0], numBasis, numQuadPts, spaceDim);
// Compute strains
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
// Update material state
_material->updateStateVars(strainCell, cell);
-
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
- dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
} // for
PYLITH_METHOD_END;
@@ -162,6 +158,7 @@
const int cellDim = _quadrature->cellDim();
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
+ const int numCorners = _quadrature->refGeometry().numCorners();
const int spaceDim = _quadrature->spaceDim();
const int tensorSize = _material->tensorSize();
totalStrain_fn_type calcTotalStrainFn;
@@ -193,34 +190,30 @@
const PetscInt numCells = _materialIS->size();
// Setup field visitors.
+ scalar_array dispCell(numBasis*spaceDim);
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
topology::VecVisitorMesh fieldVisitor(*field);
PetscScalar* fieldArray = fieldVisitor.localArray();
+ scalar_array coordsCell(numBasis*spaceDim); // :KULDGE: Update numBasis to numCorners after implementing higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
// Loop over cells
for (PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
+
// Retrieve geometry information for current cell
- PetscScalar *coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- _quadrature->computeGeometry(coordsCell, coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ _quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
const scalar_array& basisDeriv = _quadrature->basisDeriv();
// Restrict input fields to cell
- PetscScalar* dispCell = NULL;
- PetscInt dispSize = 0;
- dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
+ dispVisitor.getClosure(&dispCell, cell);
// Compute deformation tensor.
- _calcDeformation(&deformCell, basisDeriv, coordsCell, dispCell, numBasis, numQuadPts, spaceDim);
+ _calcDeformation(&deformCell, basisDeriv, &coordsCell[0], &dispCell[0], numBasis, numQuadPts, spaceDim);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
- dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
-
// Compute strains
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc 2013-07-10 22:33:03 UTC (rev 22566)
@@ -411,7 +411,9 @@
const int spaceDim = quadrature->spaceDim();
const int numBasis = quadrature->numBasis();
+ const int numCorners = quadrature->refGeometry().numCorners();
+ scalar_array coordsCell(numBasis*spaceDim); // :KULDGE: Update numBasis to numCorners after implementing higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
PylithScalar dtStable = pylith::PYLITH_MAXSCALAR;
@@ -421,11 +423,8 @@
retrievePropsAndVars(cell);
- PetscScalar *coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- const PylithScalar minCellWidth = quadrature->minCellWidth(coordsCell, numBasis, spaceDim);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ const PylithScalar minCellWidth = quadrature->minCellWidth(&coordsCell[0], numBasis, spaceDim);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const PylithScalar dt =
@@ -560,6 +559,7 @@
const int numQuadPts = quadrature->numQuadPts();
const int spaceDim = quadrature->spaceDim();
const int numBasis = quadrature->numBasis();
+ const int numCorners = quadrature->refGeometry().numCorners();
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();assert(cs);
@@ -569,6 +569,7 @@
const PetscInt* cells = _materialIS->points();
const PetscInt numCells = _materialIS->size();
+ scalar_array coordsCell(numBasis*spaceDim); // :KULDGE: Update numBasis to numCorners after implementing higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
// Create arrays for querying
@@ -629,11 +630,8 @@
const PetscInt cell = cells[c];
// Compute geometry information for current cell
- PetscScalar *coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- quadrature->computeGeometry(coordsCell, coordsSize, cell);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
// Dimensionalize coordinates for querying
const scalar_array& quadPtsNonDim = quadrature->quadPts();
@@ -693,6 +691,7 @@
const int numQuadPts = quadrature->numQuadPts();
const int spaceDim = quadrature->spaceDim();
const int numBasis = quadrature->numBasis();
+ const int numCorners = quadrature->refGeometry().numCorners();
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();assert(cs);
@@ -702,6 +701,7 @@
const PetscInt* cells = _materialIS->points();
const PetscInt numCells = _materialIS->size();
+ scalar_array coordsCell(numBasis*spaceDim); // :KULDGE: Update numBasis to numCorners after implementing higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
// Create arrays for querying
@@ -762,11 +762,8 @@
const PetscInt cell = cells[c];
// Compute geometry information for current cell
- PetscScalar* coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- quadrature->computeGeometry(coordsCell, coordsSize, cell);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
// Dimensionalize coordinates for querying
const scalar_array& quadPtsNonDim = quadrature->quadPts();
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc 2013-07-10 22:33:03 UTC (rev 22566)
@@ -123,6 +123,7 @@
// Get quadrature information
const int numQuadPts = quadrature->numQuadPts();
const int numBasis = quadrature->numBasis();
+ const int numCorners = quadrature->refGeometry().numCorners();
const int spaceDim = quadrature->spaceDim();
// Get cells associated with material
@@ -145,6 +146,7 @@
topology::VecVisitorMesh propertiesVisitor(*_properties);
PetscScalar* propertiesArray = propertiesVisitor.localArray();
+ scalar_array coordsCell(numBasis*spaceDim); // :KULDGE: Update numBasis to numCorners after implementing higher order
topology::CoordsVisitor coordsVisitor(dmMesh);
// Create arrays for querying.
@@ -199,12 +201,10 @@
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
+
// Compute geometry information for current cell
- PetscScalar* coordsCell = NULL;
- PetscInt coordsSize = 0;
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
- quadrature->computeGeometry(coordsCell, coordsSize, cell);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ coordsVisitor.getClosure(&coordsCell, cell);
+ quadrature->computeGeometry(&coordsCell[0], coordsCell.size(), cell);
const scalar_array& quadPtsNonDim = quadrature->quadPts();
quadPtsGlobal = quadPtsNonDim;
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/CoordsVisitor.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/CoordsVisitor.hh 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/CoordsVisitor.hh 2013-07-10 22:33:03 UTC (rev 22566)
@@ -29,6 +29,7 @@
#include "topologyfwd.hh" // forward declarations
#include "pylith/utils/petscfwd.h" // HASA PetscVec, PetscSection
+#include "pylith/utils/arrayfwd.hh" // USES scalar_array
// CoordsVisitor ----------------------------------------------------------
/** @brief Helper class for accessing coordinates in a finite-element mesh.
@@ -77,6 +78,8 @@
/** Get coordinates array associated with closure.
*
+ * @pre Must be followed by call to restoreClosure().
+ *
* @param coordsCell Array of coordinates for cell.
* @param coordsSize Size of coordinates array.
* @param cell Finite-element cell.
@@ -85,8 +88,18 @@
PetscInt* coordsSize,
const PetscInt cell) const;
+ /** Get coordinates array associated with closure.
+ *
+ * @param coords Array of coordinates for cell.
+ * @param cell Finite-element cell.
+ */
+ void getClosure(scalar_array* coordsCell,
+ const PetscInt cell) const;
+
/** Restore coordinates array associated with closure.
*
+ * @pre Must be preceded by call to getClosure().
+ *
* @param coordsCell Array of coordinates for cell.
* @param coordsSize Size of coordinates array.
* @param cell Finite-element cell.
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/CoordsVisitor.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/CoordsVisitor.icc 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/CoordsVisitor.icc 2013-07-10 22:33:03 UTC (rev 22566)
@@ -121,6 +121,22 @@
} // getClosure
// ----------------------------------------------------------------------
+// Get coordinates array associated with closure.
+inline
+void
+pylith::topology::CoordsVisitor::getClosure(scalar_array* coords,
+ const PetscInt cell) const
+{ // getClosure
+ assert(_dm);
+ assert(_section);
+ assert(_localVec);
+ assert(coords);
+ PetscScalar* coordsCell = &(*coords)[0];
+ PetscInt coordsSize = coords->size();
+ PetscErrorCode err = DMPlexVecGetClosure(_dm, _section, _localVec, cell, &coordsSize, &coordsCell);PYLITH_CHECK_ERROR(err);
+} // getClosure
+
+// ----------------------------------------------------------------------
// Restore coordinates array associated with closure.
inline
void
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/VisitorMesh.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/VisitorMesh.hh 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/VisitorMesh.hh 2013-07-10 22:33:03 UTC (rev 22566)
@@ -30,6 +30,7 @@
#include "topologyfwd.hh" // forward declarations
#include "pylith/utils/petscfwd.h" // HASA PetscVec, PetscSection
+#include "pylith/utils/arrayfwd.hh" // USES scalar_array
// VecVisitorMesh ----------------------------------------------------------
/** @brief Helper class for accessing field values at points in a
@@ -101,7 +102,7 @@
/** Get array of values associated with closure.
*
- * @pre Must be followed by call to getClosure().
+ * @pre Must be followed by call to restoreClosure().
*
* @param valuesCell Array of values for cell.
* @param valuesSize Size of values array.
@@ -111,6 +112,14 @@
PetscInt* valuesSize,
const PetscInt cell) const;
+ /** Get array of values associated with closure.
+ *
+ * @param values Array of values for cell.
+ * @param cell Finite-element cell.
+ */
+ void getClosure(scalar_array* values,
+ const PetscInt cell) const;
+
/** Restore array of values associated with closure.
*
* @pre Must be preceded by call to getClosure().
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/VisitorMesh.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/VisitorMesh.icc 2013-07-10 16:32:56 UTC (rev 22565)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/VisitorMesh.icc 2013-07-10 22:33:03 UTC (rev 22566)
@@ -153,6 +153,22 @@
} // getClosure
// ----------------------------------------------------------------------
+// Get coordinates array associated with closure.
+inline
+void
+pylith::topology::VecVisitorMesh::getClosure(scalar_array* values,
+ const PetscInt cell) const
+{ // getClosure
+ assert(_dm);
+ assert(_section);
+ assert(_localVec);
+ assert(values);
+ PetscScalar* valuesCell = &(*values)[0];
+ PetscInt valuesSize = values->size();
+ PetscErrorCode err = DMPlexVecGetClosure(_dm, _section, _localVec, cell, &valuesSize, &valuesCell);PYLITH_CHECK_ERROR(err);
+} // getClosure
+
+// ----------------------------------------------------------------------
// Restore coordinates array associated with closure.
inline
void
More information about the CIG-COMMITS
mailing list