[cig-commits] r16453 - in short/3D/PyLith/trunk/libsrc: bc feassemble
brad at geodynamics.org
brad at geodynamics.org
Wed Mar 24 17:14:12 PDT 2010
Author: brad
Date: 2010-03-24 17:14:12 -0700 (Wed, 24 Mar 2010)
New Revision: 16453
Modified:
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc
Log:
Fixed bugs in computing residual for explicit time stepping. Missing dispIncr term (only showed up when using nonlinear solver).
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2010-03-24 19:22:49 UTC (rev 16452)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2010-03-25 00:14:12 UTC (rev 16453)
@@ -38,6 +38,7 @@
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Mesh::UpdateAddVisitor UpdateAddVisitor;
// ----------------------------------------------------------------------
// Default constructor.
@@ -272,6 +273,7 @@
// Allocate vectors for cell values.
_initCellVector();
double_array dampersCell(numQuadPts*spaceDim);
+ double_array velCell(numBasis*spaceDim);
// Get cell information
const ALE::Obj<SieveSubMesh>& sieveSubMesh = _boundaryMesh->sieveMesh();
@@ -287,24 +289,30 @@
_parameters->get("damping constants").section();
assert(!dampersSection.isNull());
+ // Use _cellVector for cell residual.
const ALE::Obj<RealSection>& residualSection = residual.section();
assert(!residualSection.isNull());
- topology::SubMesh::UpdateAddVisitor residualVisitor(*residualSection,
- &_cellVector[0]);
-
- const topology::Field<topology::Mesh>& dispT =
- fields->get("disp(t)");
- const ALE::Obj<RealSection>& dispTSection = dispT.section();
+ UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
+
+ double_array dispTIncrCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispTIncrSection =
+ fields->get("dispIncr(t->t+dt)").section();
+ assert(!dispTIncrSection.isNull());
+ RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
+ dispTIncrCell.size(), &dispTIncrCell[0]);
+
+ double_array dispTCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
assert(!dispTSection.isNull());
- topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection,
- numBasis*spaceDim);
- const topology::Field<topology::Mesh>& dispTmdt =
- fields->get("disp(t-dt)");
- const ALE::Obj<RealSection>& dispTmdtSection = dispTmdt.section();
+ RestrictVisitor dispTVisitor(*dispTSection, dispTCell.size(), &dispTCell[0]);
+
+ double_array dispTmdtCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispTmdtSection =
+ fields->get("disp(t-dt)").section();
assert(!dispTmdtSection.isNull());
- topology::Mesh::RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
- numBasis*spaceDim);
-
+ RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
+ dispTmdtCell.size(), &dispTmdtCell[0]);
+
#if !defined(PRECOMPUTE_GEOMETRY)
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
@@ -346,13 +354,14 @@
_resetCellVector();
// Restrict input fields to cell
+ dispTIncrVisitor.clear();
+ sieveSubMesh->restrictClosure(*c_iter, dispTIncrVisitor);
+
dispTVisitor.clear();
sieveSubMesh->restrictClosure(*c_iter, dispTVisitor);
- const double* dispTCell = dispTVisitor.getValues();
dispTmdtVisitor.clear();
sieveSubMesh->restrictClosure(*c_iter, dispTmdtVisitor);
- const double* dispTmdtCell = dispTmdtVisitor.getValues();
dampersSection->restrictPoint(*c_iter, &dampersCell[0], dampersCell.size());
@@ -366,25 +375,26 @@
const double_array& jacobianDet = _quadrature->jacobianDet();
// Compute action for absorbing bc terms
+ velCell = (dispTCell + dispTIncrCell - dispTmdtCell) / (2.0*dt);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt =
- quadWts[iQuad] * jacobianDet[iQuad] / (2.0 * dt);
+ quadWts[iQuad] * jacobianDet[iQuad];
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] +=
+ _cellVector[iBasis*spaceDim+iDim] -=
dampersCell[iQuad*spaceDim+iDim] *
- valIJ * (dispTmdtCell[jBasis*spaceDim+iDim] - \
- dispTCell[jBasis*spaceDim+iDim]);
+ valIJ * velCell[jBasis*spaceDim+iDim];
} // for
} // for
} // for
#if defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(5*spaceDim))));
+ PetscLogFlops(numBasis*spaceDim*3 +
+ numQuadPts*(1+numBasis*(1+numBasis*(3*spaceDim))));
_logger->eventEnd(computeEvent);
_logger->eventBegin(updateEvent);
#endif
@@ -399,7 +409,8 @@
} // for
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*numQuadPts*(3+numBasis*(1+numBasis*(5*spaceDim))));
+ PetscLogFlops(cells->size()*(numBasis*spaceDim*3 +
+ numQuadPts*(1+numBasis*(1+numBasis*(3*spaceDim)))));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidual
@@ -436,6 +447,7 @@
// Allocate vectors for cell values.
_initCellVector();
double_array dampersCell(numQuadPts*spaceDim);
+ double_array velCell(numBasis*spaceDim);
// Get cell information
const ALE::Obj<SieveSubMesh>& sieveSubMesh = _boundaryMesh->sieveMesh();
@@ -451,30 +463,36 @@
_parameters->get("damping constants").section();
assert(!dampersSection.isNull());
+ // Use _cellVector for cell values.
const ALE::Obj<RealSection>& residualSection = residual.section();
assert(!residualSection.isNull());
- topology::SubMesh::UpdateAddVisitor residualVisitor(*residualSection,
- &_cellVector[0]);
+ UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
- const topology::Field<topology::Mesh>& dispT =
- fields->get("disp(t)");
- const ALE::Obj<RealSection>& dispTSection = dispT.section();
+ double_array dispTIncrCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispTIncrSection =
+ fields->get("dispIncr(t->t+dt)").section();
+ assert(!dispTIncrSection.isNull());
+ RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
+ dispTIncrCell.size(), &dispTIncrCell[0]);
+
+ double_array dispTCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
assert(!dispTSection.isNull());
- topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection,
- numBasis*spaceDim);
- const topology::Field<topology::Mesh>& dispTmdt =
- fields->get("disp(t-dt)");
- const ALE::Obj<RealSection>& dispTmdtSection = dispTmdt.section();
+ RestrictVisitor dispTVisitor(*dispTSection, dispTCell.size(), &dispTCell[0]);
+
+ double_array dispTmdtCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispTmdtSection =
+ fields->get("disp(t-dt)").section();
assert(!dispTmdtSection.isNull());
- topology::Mesh::RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
- numBasis*spaceDim);
+ RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
+ dispTmdtCell.size(), &dispTmdtCell[0]);
#if !defined(PRECOMPUTE_GEOMETRY)
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveSubMesh->getRealSection("coordinates");
RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ coordinatesCell.size(), &coordinatesCell[0]);
#endif
_logger->eventEnd(setupEvent);
@@ -512,11 +530,9 @@
// Restrict input fields to cell
dispTVisitor.clear();
sieveSubMesh->restrictClosure(*c_iter, dispTVisitor);
- const double* dispTCell = dispTVisitor.getValues();
dispTmdtVisitor.clear();
sieveSubMesh->restrictClosure(*c_iter, dispTmdtVisitor);
- const double* dispTmdtCell = dispTmdtVisitor.getValues();
dampersSection->restrictPoint(*c_iter, &dampersCell[0], dampersCell.size());
@@ -530,8 +546,9 @@
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);
+ velCell = (dispTCell + dispTIncrCell - dispTmdtCell) / (2.0*dt);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
const int iQ = iQuad * numBasis;
double valJ = 0.0;
for (int jBasis = 0; jBasis < numBasis; ++jBasis)
@@ -540,13 +557,13 @@
for (int iBasis = 0; iBasis < numBasis; ++iBasis) {
const double valIJ = basis[iQ + iBasis] * valJ;
for (int iDim = 0; iDim < spaceDim; ++iDim)
- _cellVector[iBasis * spaceDim + iDim] += valIJ * dampersCell[iQuad
- * spaceDim + iDim] * (dispTmdtCell[iBasis * spaceDim + iDim]
- - dispTCell[iBasis * spaceDim + iDim]);
+ _cellVector[iBasis*spaceDim+iDim] -= valIJ *
+ dampersCell[iQuad*spaceDim+iDim] * velCell[iBasis*spaceDim+iDim];
} // for
} // for
#if defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(numQuadPts*(4+numBasis+numBasis*(1+spaceDim*4)));
+ PetscLogFlops(numBasis*spaceDim*3 +
+ numQuadPts*(1+numBasis+numBasis*(1+spaceDim*3)));
_logger->eventEnd(computeEvent);
_logger->eventBegin(updateEvent);
#endif
@@ -560,7 +577,8 @@
#endif
} // for
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*numQuadPts*(4+numBasis+numBasis*(1+spaceDim*4)));
+ PetscLogFlops(cells->size()*(numBasis*spaceDim*3 +
+ numQuadPts*(1+numBasis+numBasis*(1+spaceDim*3))));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidualLumped
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2010-03-24 19:22:49 UTC (rev 16452)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2010-03-25 00:14:12 UTC (rev 16453)
@@ -40,6 +40,8 @@
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
+typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Mesh::UpdateAddVisitor UpdateAddVisitor;
// ----------------------------------------------------------------------
// Constructor
@@ -166,12 +168,11 @@
} // switch
// Allocate vectors for cell values.
- double_array dispTCell(numBasis*spaceDim);
- double_array dispTmdtCell(numBasis*spaceDim);
double_array strainCell(numQuadPts*tensorSize);
strainCell = 0.0;
double_array gravVec(spaceDim);
double_array quadPtsGlobal(numQuadPts*spaceDim);
+ double_array accCell(numBasis*spaceDim);
// Get cell information
const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
@@ -184,28 +185,35 @@
const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
// Get sections
+ double_array dispTIncrCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispTIncrSection =
+ fields->get("dispIncr(t->t+dt)").section();
+ assert(!dispTIncrSection.isNull());
+ RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
+ dispTIncrCell.size(), &dispTIncrCell[0]);
+
+ double_array dispTCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
assert(!dispTSection.isNull());
- topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection,
- numBasis*spaceDim,
- &dispTCell[0]);
+ RestrictVisitor dispTVisitor(*dispTSection,
+ dispTCell.size(), &dispTCell[0]);
+
+ double_array dispTmdtCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& dispTmdtSection =
fields->get("disp(t-dt)").section();
assert(!dispTmdtSection.isNull());
- topology::Mesh::RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
- numBasis*spaceDim,
- &dispTmdtCell[0]);
+ RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
+ dispTmdtCell.size(), &dispTmdtCell[0]);
+
const ALE::Obj<RealSection>& residualSection = residual.section();
- topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
- &_cellVector[0]);
+ UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
assert(0 != _normalizer);
const double lengthScale = _normalizer->lengthScale();
@@ -245,6 +253,7 @@
// Get state variables for cell.
_material->retrievePropsAndVars(*c_iter);
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(stateVarsEvent);
_logger->eventBegin(restrictEvent);
@@ -254,8 +263,12 @@
_resetCellVector();
// Restrict input fields to cell
+ dispTIncrVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTIncrVisitor);
+
dispTVisitor.clear();
sieveMesh->restrictClosure(*c_iter, dispTVisitor);
+
dispTmdtVisitor.clear();
sieveMesh->restrictClosure(*c_iter, dispTmdtVisitor);
@@ -302,22 +315,23 @@
} // if
// Compute action for inertial terms
+ accCell = (dispTIncrCell - dispTCell + dispTmdtCell) / dt2;
const double_array& density = _material->calcDensity();
for (int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad]
- / dt2;
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad];
for (int iBasis = 0; iBasis < numBasis; ++iBasis) {
- const double valI = wt * basis[iQuad * 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];
+ const double valIJ = valI * basis[iQuad*numBasis+jBasis];
for (int iDim = 0; iDim < spaceDim; ++iDim)
- _cellVector[iBasis * spaceDim + iDim] += valIJ * (dispTCell[jBasis
- * spaceDim + iDim] - dispTmdtCell[jBasis * spaceDim + iDim]);
+ _cellVector[iBasis*spaceDim+iDim] -= valIJ *
+ accCell[jBasis*spaceDim+iDim];
} // for
} // for
} // for
#if defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(6*spaceDim))));
+ PetscLogFlops(numBasis*spaceDim*2 +
+ numQuadPts*(3+numBasis*(1+numBasis*(2*spaceDim))));
_logger->eventEnd(computeEvent);
_logger->eventBegin(stressEvent);
#endif
@@ -348,7 +362,8 @@
} // for
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*numQuadPts*(3+numBasis*(1+numBasis*(6*spaceDim))));
+ PetscLogFlops(cells->size()*(numBasis*spaceDim*2 +
+ numQuadPts*(3+numBasis*(1+numBasis*(2*spaceDim)))));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidual
@@ -424,12 +439,11 @@
} // if/else
// Allocate vectors for cell values.
- double_array dispTCell(numBasis*spaceDim);
- double_array dispTmdtCell(numBasis*spaceDim);
double_array strainCell(numQuadPts*tensorSize);
strainCell = 0.0;
double_array gravVec(spaceDim);
double_array quadPtsGlobal(numQuadPts*spaceDim);
+ double_array accCell(numBasis*spaceDim);
// Get cell information
const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
@@ -442,28 +456,35 @@
const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
// Get sections
+ double_array dispTIncrCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispTIncrSection =
+ fields->get("dispIncr(t->t+dt)").section();
+ assert(!dispTIncrSection.isNull());
+ RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
+ dispTIncrCell.size(), &dispTIncrCell[0]);
+
+ double_array dispTCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
assert(!dispTSection.isNull());
- topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection,
- numBasis*spaceDim,
- &dispTCell[0]);
+ RestrictVisitor dispTVisitor(*dispTSection,
+ dispTCell.size(), &dispTCell[0]);
+
+ double_array dispTmdtCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& dispTmdtSection =
fields->get("disp(t-dt)").section();
assert(!dispTmdtSection.isNull());
- topology::Mesh::RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
- numBasis*spaceDim,
- &dispTmdtCell[0]);
+ RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
+ dispTmdtCell.size(), &dispTmdtCell[0]);
+
const ALE::Obj<RealSection>& residualSection = residual.section();
- topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
- &_cellVector[0]);
-
+ UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
+
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
assert(0 != _normalizer);
const double lengthScale = _normalizer->lengthScale();
@@ -514,10 +535,15 @@
_resetCellVector();
// Restrict input fields to cell
+ dispTIncrVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTIncrVisitor);
+
dispTVisitor.clear();
sieveMesh->restrictClosure(*c_iter, dispTVisitor);
+
dispTmdtVisitor.clear();
sieveMesh->restrictClosure(*c_iter, dispTmdtVisitor);
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -553,7 +579,7 @@
for (int iBasis = 0, iQ = iQuad * numBasis; iBasis < numBasis; ++iBasis) {
const double valI = wt * basis[iQ + iBasis];
for (int iDim = 0; iDim < spaceDim; ++iDim) {
- _cellVector[iBasis * spaceDim + iDim] += valI * gravVec[iDim];
+ _cellVector[iBasis*spaceDim+iDim] += valI * gravVec[iDim];
} // for
} // for
} // for
@@ -561,11 +587,11 @@
} // if
// Compute action for inertial terms
+ accCell = (dispTIncrCell - dispTCell + dispTmdtCell) / dt2;
const double_array& density = _material->calcDensity();
valuesIJ = 0.0;
for (int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad]
- / dt2;
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad];
const int iQ = iQuad * numBasis;
double valJ = 0.0;
for (int jBasis = 0; jBasis < numBasis; ++jBasis)
@@ -576,8 +602,8 @@
} // for
for (int iBasis = 0; iBasis < numBasis; ++iBasis)
for (int iDim = 0; iDim < spaceDim; ++iDim)
- _cellVector[iBasis*spaceDim+iDim] += valuesIJ[iBasis] *
- (dispTCell[iBasis*spaceDim+iDim] - dispTmdtCell[iBasis*spaceDim+iDim]);
+ _cellVector[iBasis*spaceDim+iDim] -= valuesIJ[iBasis] *
+ accCell[iBasis*spaceDim+iDim];
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numQuadPts*(4+numBasis*3) + numBasis*spaceDim*3);
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc 2010-03-24 19:22:49 UTC (rev 16452)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc 2010-03-25 00:14:12 UTC (rev 16453)
@@ -38,6 +38,8 @@
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
+typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Mesh::UpdateAddVisitor UpdateAddVisitor;
// ----------------------------------------------------------------------
const int pylith::feassemble::ElasticityExplicitTet4::_spaceDim = 3;
@@ -145,12 +147,11 @@
"domain not implemented yet.");
// Allocate vectors for cell values.
- double_array dispTCell(numBasis*spaceDim);
- double_array dispTmdtCell(numBasis*spaceDim);
double_array strainCell(numQuadPts*tensorSize);
strainCell = 0.0;
double_array gravVec(spaceDim);
double_array quadPtsGlobal(numQuadPts*spaceDim);
+ double_array accCell(numBasis*spaceDim);
// Get cell information
const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
@@ -163,28 +164,33 @@
const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
// Get sections
+ double_array dispTIncrCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispTIncrSection =
+ fields->get("dispIncr(t->t+dt)").section();
+ assert(!dispTIncrSection.isNull());
+ RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
+ dispTIncrCell.size(), &dispTIncrCell[0]);
+
+ double_array dispTCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
assert(!dispTSection.isNull());
- topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection,
- numBasis*spaceDim,
- &dispTCell[0]);
+ RestrictVisitor dispTVisitor(*dispTSection, dispTCell.size(), &dispTCell[0]);
+
+ double_array dispTmdtCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& dispTmdtSection =
fields->get("disp(t-dt)").section();
assert(!dispTmdtSection.isNull());
- topology::Mesh::RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
- numBasis*spaceDim,
- &dispTmdtCell[0]);
+ RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
+ dispTmdtCell.size(), &dispTmdtCell[0]);
const ALE::Obj<RealSection>& residualSection = residual.section();
- topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
- &_cellVector[0]);
-
+ UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
+
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
assert(0 != _normalizer);
const double lengthScale = _normalizer->lengthScale();
@@ -233,8 +239,12 @@
_resetCellVector();
// Restrict input fields to cell
+ dispTIncrVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTIncrVisitor);
+
dispTVisitor.clear();
sieveMesh->restrictClosure(*c_iter, dispTVisitor);
+
dispTmdtVisitor.clear();
sieveMesh->restrictClosure(*c_iter, dispTmdtVisitor);
@@ -273,15 +283,16 @@
} // if
// Compute action for inertial terms
- const double wtVertex = density[0] * volume / (16.0 * dt2);
+ accCell = (dispTIncrCell - dispTCell + dispTmdtCell) / dt2;
+ const double wtVertex = density[0] * volume / 16.0;
for (int iBasis = 0; iBasis < numBasis; ++iBasis)
for (int jBasis = 0; jBasis < numBasis; ++jBasis)
for (int iDim = 0; iDim < spaceDim; ++iDim)
- _cellVector[iBasis * spaceDim + iDim] += wtVertex * (dispTCell[jBasis
- * spaceDim + iDim] - dispTmdtCell[jBasis * spaceDim + iDim]);
+ _cellVector[iBasis*spaceDim+iDim] -=
+ wtVertex * accCell[jBasis*spaceDim+iDim];
#if defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(3 + numBasis*numBasis*spaceDim*3);
+ PetscLogFlops(3 + numBasis*spaceDim*3 + numBasis*numBasis*spaceDim*2);
_logger->eventEnd(computeEvent);
_logger->eventBegin(stressEvent);
#endif
@@ -375,7 +386,10 @@
} // for
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*(196+84));
+ PetscLogFlops(cells->size()*(3 +
+ numBasis*spaceDim*3 +
+ numBasis*numBasis*spaceDim*2 +
+ 196+84));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidual
More information about the CIG-COMMITS
mailing list