[cig-commits] r16475 - in short/3D/PyLith/trunk: . libsrc libsrc/bc libsrc/feassemble libsrc/meshio libsrc/problems libsrc/topology modulesrc/topology pylith pylith/topology
brad at geodynamics.org
brad at geodynamics.org
Wed Mar 31 12:08:40 PDT 2010
Author: brad
Date: 2010-03-31 12:08:39 -0700 (Wed, 31 Mar 2010)
New Revision: 16475
Added:
short/3D/PyLith/trunk/libsrc/topology/ReverseCuthillMcKee.cc
short/3D/PyLith/trunk/libsrc/topology/ReverseCuthillMcKee.hh
short/3D/PyLith/trunk/modulesrc/topology/ReverseCuthillMcKee.i
short/3D/PyLith/trunk/pylith/topology/ReverseCuthillMcKee.py
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc
short/3D/PyLith/trunk/libsrc/meshio/MeshBuilder.cc
short/3D/PyLith/trunk/libsrc/problems/Explicit.cc
short/3D/PyLith/trunk/libsrc/topology/Makefile.am
short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh
short/3D/PyLith/trunk/modulesrc/topology/Makefile.am
short/3D/PyLith/trunk/modulesrc/topology/topology.i
short/3D/PyLith/trunk/pylith/Makefile.am
short/3D/PyLith/trunk/pylith/topology/MeshImporter.py
Log:
Move mesh reordering to ReverseCuthillMcKee object. Integrated mesh ordering into MeshImporter.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2010-03-31 16:52:14 UTC (rev 16474)
+++ short/3D/PyLith/trunk/TODO 2010-03-31 19:08:39 UTC (rev 16475)
@@ -12,7 +12,7 @@
+ Need settings for Schur complement
* Drucker-Prager elastoplastic
- Add isJacobianSymmetric() to Integrator (default is true, need to query material model).
+ Need to link isJacobianSymmetric() in Integrator with materials.
Need to pass hint when creating Jacobian.
* Friction
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2010-03-31 16:52:14 UTC (rev 16474)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2010-03-31 19:08:39 UTC (rev 16475)
@@ -126,6 +126,7 @@
topology/SolutionFields.cc \
topology/Distributor.cc \
topology/RefineUniform.cc \
+ topology/ReverseCuthillMcKee.cc \
utils/EventLogger.cc \
utils/TestArray.cc
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2010-03-31 16:52:14 UTC (rev 16474)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2010-03-31 19:08:39 UTC (rev 16475)
@@ -273,7 +273,6 @@
// 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();
@@ -294,25 +293,12 @@
assert(!residualSection.isNull());
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 velCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& velSection =
+ fields->get("velocity(t)").section();
+ assert(!velSection.isNull());
+ RestrictVisitor velVisitor(*velSection, velCell.size(), &velCell[0]);
- double_array dispTCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
- 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());
- RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
- dispTmdtCell.size(), &dispTmdtCell[0]);
-
#if !defined(PRECOMPUTE_GEOMETRY)
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
@@ -326,10 +312,6 @@
_logger->eventBegin(computeEvent);
#endif
- // Get parameters used in integration.
- const double dt = _dt;
- assert(dt > 0);
-
for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
c_iter != cellsEnd;
++c_iter) {
@@ -354,15 +336,9 @@
_resetCellVector();
// Restrict input fields to cell
- dispTIncrVisitor.clear();
- sieveSubMesh->restrictClosure(*c_iter, dispTIncrVisitor);
+ velVisitor.clear();
+ sieveSubMesh->restrictClosure(*c_iter, velVisitor);
- dispTVisitor.clear();
- sieveSubMesh->restrictClosure(*c_iter, dispTVisitor);
-
- dispTmdtVisitor.clear();
- sieveSubMesh->restrictClosure(*c_iter, dispTmdtVisitor);
-
dampersSection->restrictPoint(*c_iter, &dampersCell[0], dampersCell.size());
#if defined(DETAILED_EVENT_LOGGING)
@@ -375,10 +351,8 @@
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];
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
for (int iBasis=0; iBasis < numBasis; ++iBasis) {
const double valI = wt*basis[iQuad*numBasis+iBasis];
@@ -393,8 +367,7 @@
} // for
#if defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(numBasis*spaceDim*3 +
- numQuadPts*(1+numBasis*(1+numBasis*(3*spaceDim))));
+ PetscLogFlops(numQuadPts*(1+numBasis*(1+numBasis*(3*spaceDim))));
_logger->eventEnd(computeEvent);
_logger->eventBegin(updateEvent);
#endif
@@ -409,8 +382,7 @@
} // for
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*(numBasis*spaceDim*3 +
- numQuadPts*(1+numBasis*(1+numBasis*(3*spaceDim)))));
+ PetscLogFlops(cells->size()*numQuadPts*(1+numBasis*(1+numBasis*(3*spaceDim))));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidual
@@ -447,7 +419,6 @@
// 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();
@@ -468,25 +439,12 @@
assert(!residualSection.isNull());
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 velCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& velSection =
+ fields->get("velocity(t)").section();
+ assert(!velSection.isNull());
+ RestrictVisitor velVisitor(*velSection, velCell.size(), &velCell[0]);
- double_array dispTCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
- 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());
- RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
- dispTmdtCell.size(), &dispTmdtCell[0]);
-
#if !defined(PRECOMPUTE_GEOMETRY)
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
@@ -500,10 +458,6 @@
_logger->eventBegin(computeEvent);
#endif
- // Get parameters used in integration.
- const double dt = _dt;
- assert(dt > 0);
-
for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
c_iter != cellsEnd;
++c_iter) {
@@ -528,12 +482,9 @@
_resetCellVector();
// Restrict input fields to cell
- dispTVisitor.clear();
- sieveSubMesh->restrictClosure(*c_iter, dispTVisitor);
+ velVisitor.clear();
+ sieveSubMesh->restrictClosure(*c_iter, velVisitor);
- dispTmdtVisitor.clear();
- sieveSubMesh->restrictClosure(*c_iter, dispTmdtVisitor);
-
dampersSection->restrictPoint(*c_iter, &dampersCell[0], dampersCell.size());
#if defined(DETAILED_EVENT_LOGGING)
@@ -546,7 +497,6 @@
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];
const int iQ = iQuad * numBasis;
@@ -562,8 +512,7 @@
} // for
} // for
#if defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(numBasis*spaceDim*3 +
- numQuadPts*(1+numBasis+numBasis*(1+spaceDim*3)));
+ PetscLogFlops(numQuadPts*(1+numBasis+numBasis*(1+spaceDim*3)));
_logger->eventEnd(computeEvent);
_logger->eventBegin(updateEvent);
#endif
@@ -577,8 +526,7 @@
#endif
} // for
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*(numBasis*spaceDim*3 +
- numQuadPts*(1+numBasis+numBasis*(1+spaceDim*3))));
+ PetscLogFlops(cells->size()*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-31 16:52:14 UTC (rev 16474)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2010-03-31 19:08:39 UTC (rev 16475)
@@ -172,7 +172,6 @@
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();
@@ -185,26 +184,18 @@
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 accCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& accSection =
+ fields->get("acceleration(t)").section();
+ assert(!accSection.isNull());
+ RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
- double_array dispTCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
- 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());
- RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
- dispTmdtCell.size(), &dispTmdtCell[0]);
-
+ double_array dispCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispSection =
+ fields->get("disp(t)").section();
+ assert(!dispSection.isNull());
+ RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
+
const ALE::Obj<RealSection>& residualSection = residual.section();
UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
@@ -221,11 +212,6 @@
_normalizer->pressureScale() / (_normalizer->lengthScale() *
_normalizer->densityScale());
- // Get parameters used in integration.
- const double dt = _dt;
- const double dt2 = dt*dt;
- assert(dt > 0);
-
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(computeEvent);
@@ -263,15 +249,12 @@
_resetCellVector();
// Restrict input fields to cell
- dispTIncrVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTIncrVisitor);
+ accVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, accVisitor);
- dispTVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTVisitor);
+ dispVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispVisitor);
- dispTmdtVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTmdtVisitor);
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -315,7 +298,6 @@
} // 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];
@@ -330,14 +312,13 @@
} // for
} // for
#if defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(numBasis*spaceDim*2 +
- numQuadPts*(3+numBasis*(1+numBasis*(2*spaceDim))));
+ PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(2*spaceDim))));
_logger->eventEnd(computeEvent);
_logger->eventBegin(stressEvent);
#endif
// Compute B(transpose) * sigma, first computing strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispTCell,
+ calcTotalStrainFn(&strainCell, basisDeriv, dispCell,
numBasis, numQuadPts);
const double_array& stressCell = _material->calcStress(strainCell, true);
@@ -362,8 +343,7 @@
} // for
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*(numBasis*spaceDim*2 +
- numQuadPts*(3+numBasis*(1+numBasis*(2*spaceDim)))));
+ PetscLogFlops(cells->size()*numQuadPts*(3+numBasis*(1+numBasis*(2*spaceDim))));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidual
@@ -443,7 +423,6 @@
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();
@@ -456,26 +435,18 @@
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 accCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& accSection =
+ fields->get("acceleration(t)").section();
+ assert(!accSection.isNull());
+ RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
- double_array dispTCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
- RestrictVisitor dispTVisitor(*dispTSection,
- dispTCell.size(), &dispTCell[0]);
+ double_array dispCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispSection =
+ fields->get("disp(t)").section();
+ assert(!dispSection.isNull());
+ RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
- double_array dispTmdtCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispTmdtSection =
- fields->get("disp(t-dt)").section();
- assert(!dispTmdtSection.isNull());
- RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
- dispTmdtCell.size(), &dispTmdtCell[0]);
-
const ALE::Obj<RealSection>& residualSection = residual.section();
UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
@@ -493,9 +464,6 @@
_normalizer->densityScale());
// Get parameters used in integration.
- const double dt = _dt;
- const double dt2 = dt*dt;
- assert(dt > 0);
double_array valuesIJ(numBasis);
_logger->eventEnd(setupEvent);
@@ -535,15 +503,12 @@
_resetCellVector();
// Restrict input fields to cell
- dispTIncrVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTIncrVisitor);
+ accVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, accVisitor);
- dispTVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTVisitor);
+ dispVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispVisitor);
- dispTmdtVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTmdtVisitor);
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -587,7 +552,6 @@
} // 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) {
@@ -606,13 +570,13 @@
accCell[iBasis*spaceDim+iDim];
#if defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(numQuadPts*(4+numBasis*3) + numBasis*spaceDim*3);
+ PetscLogFlops(numQuadPts*(4+numBasis*3));
_logger->eventEnd(computeEvent);
_logger->eventBegin(stressEvent);
#endif
// Compute B(transpose) * sigma, first computing strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispTCell,
+ calcTotalStrainFn(&strainCell, basisDeriv, dispCell,
numBasis, numQuadPts);
const double_array& stressCell = _material->calcStress(strainCell, true);
@@ -638,7 +602,7 @@
} // for
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*(numQuadPts*(4+numBasis*3)+numBasis*spaceDim*3));
+ PetscLogFlops(cells->size()*numQuadPts*(4+numBasis*3));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidualLumped
@@ -678,9 +642,6 @@
"contribution to Jacobian matrix for cells with " \
"different dimensions than the spatial dimension.");
- // Allocate vectors for cell data.
- double_array dispTCell(numBasis*spaceDim);
-
// Get cell information
const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
assert(!sieveMesh.isNull());
@@ -692,9 +653,10 @@
const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
// Get sections
- const ALE::Obj<RealSection>& dispTSection =
+ double_array dispCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispSection =
fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
+ assert(!dispSection.isNull());
// Get sparse matrix
const PetscMat jacobianMat = jacobian->matrix();
@@ -706,10 +668,10 @@
assert(dt > 0);
const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", dispTSection);
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", dispSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
- topology::Mesh::IndicesVisitor jacobianVisitor(*dispTSection, *globalOrder,
+ topology::Mesh::IndicesVisitor jacobianVisitor(*dispSection, *globalOrder,
(int) pow(sieveMesh->getSieve()->getMaxConeSize(),
sieveMesh->depth())*spaceDim);
@@ -860,10 +822,10 @@
double_array valuesIJ(numBasis);
// Get sections
- double_array dispTCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispTSection =
+ double_array dispCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispSection =
fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
+ assert(!dispSection.isNull());
const ALE::Obj<RealSection>& jacobianSection = jacobian->section();
assert(!jacobianSection.isNull());
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc 2010-03-31 16:52:14 UTC (rev 16474)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc 2010-03-31 19:08:39 UTC (rev 16475)
@@ -151,7 +151,6 @@
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();
@@ -164,24 +163,18 @@
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 accCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& accSection =
+ fields->get("acceleration(t)").section();
+ assert(!accSection.isNull());
+ RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
- double_array dispTCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
- RestrictVisitor dispTVisitor(*dispTSection, dispTCell.size(), &dispTCell[0]);
+ double_array dispCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispSection =
+ fields->get("disp(t)").section();
+ assert(!dispSection.isNull());
+ RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
- double_array dispTmdtCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispTmdtSection =
- fields->get("disp(t-dt)").section();
- assert(!dispTmdtSection.isNull());
- RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
- dispTmdtCell.size(), &dispTmdtCell[0]);
const ALE::Obj<RealSection>& residualSection = residual.section();
UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
@@ -198,11 +191,6 @@
_normalizer->pressureScale() / (_normalizer->lengthScale() *
_normalizer->densityScale());
- // Get parameters used in integration.
- const double dt = _dt;
- const double dt2 = dt*dt;
- assert(dt > 0);
-
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(computeEvent);
@@ -239,15 +227,12 @@
_resetCellVector();
// Restrict input fields to cell
- dispTIncrVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTIncrVisitor);
+ accVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, accVisitor);
- dispTVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTVisitor);
-
- dispTmdtVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTmdtVisitor);
-
+ dispVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispVisitor);
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -283,7 +268,6 @@
} // if
// Compute action for inertial terms
- 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)
@@ -292,7 +276,7 @@
wtVertex * accCell[jBasis*spaceDim+iDim];
#if defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(3 + numBasis*spaceDim*3 + numBasis*numBasis*spaceDim*2);
+ PetscLogFlops(3 + numBasis*numBasis*spaceDim*2);
_logger->eventEnd(computeEvent);
_logger->eventBegin(stressEvent);
#endif
@@ -332,21 +316,21 @@
const double d4 = (x1*y2+x0*(y1-y2)-x2*y1-(x1-x2)*y0) / scaleB;
assert(strainCell.size() == 6);
- strainCell[0] = b1 * dispTCell[0] + b2 * dispTCell[3] + b3 * dispTCell[6]
- + b4 * dispTCell[9];
- strainCell[1] = c3 * dispTCell[7] + c2 * dispTCell[4] + c4 * dispTCell[10]
- + c1 * dispTCell[1];
- strainCell[2] = d3 * dispTCell[8] + d2 * dispTCell[5] + d1 * dispTCell[2]
- + d4 * dispTCell[11];
- strainCell[3] = (c4 * dispTCell[9] + b3 * dispTCell[7] + c3 * dispTCell[6]
- + b2 * dispTCell[4] + c2 * dispTCell[3] + b4 * dispTCell[10] + b1
- * dispTCell[1] + c1 * dispTCell[0]) / 2.0;
- strainCell[4] = (c3 * dispTCell[8] + d3 * dispTCell[7] + c2 * dispTCell[5]
- + d2 * dispTCell[4] + c1 * dispTCell[2] + c4 * dispTCell[11] + d4
- * dispTCell[10] + d1 * dispTCell[1]) / 2.0;
- strainCell[5] = (d4 * dispTCell[9] + b3 * dispTCell[8] + d3 * dispTCell[6]
- + b2 * dispTCell[5] + d2 * dispTCell[3] + b1 * dispTCell[2] + b4
- * dispTCell[11] + d1 * dispTCell[0]) / 2.0;
+ strainCell[0] = b1 * dispCell[0] + b2 * dispCell[3] + b3 * dispCell[6]
+ + b4 * dispCell[9];
+ strainCell[1] = c3 * dispCell[7] + c2 * dispCell[4] + c4 * dispCell[10]
+ + c1 * dispCell[1];
+ strainCell[2] = d3 * dispCell[8] + d2 * dispCell[5] + d1 * dispCell[2]
+ + d4 * dispCell[11];
+ strainCell[3] = (c4 * dispCell[9] + b3 * dispCell[7] + c3 * dispCell[6]
+ + b2 * dispCell[4] + c2 * dispCell[3] + b4 * dispCell[10] + b1
+ * dispCell[1] + c1 * dispCell[0]) / 2.0;
+ strainCell[4] = (c3 * dispCell[8] + d3 * dispCell[7] + c2 * dispCell[5]
+ + d2 * dispCell[4] + c1 * dispCell[2] + c4 * dispCell[11] + d4
+ * dispCell[10] + d1 * dispCell[1]) / 2.0;
+ strainCell[5] = (d4 * dispCell[9] + b3 * dispCell[8] + d3 * dispCell[6]
+ + b2 * dispCell[5] + d2 * dispCell[3] + b1 * dispCell[2] + b4
+ * dispCell[11] + d1 * dispCell[0]) / 2.0;
const double_array& stressCell = _material->calcStress(strainCell, true);
@@ -386,10 +370,7 @@
} // for
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*(3 +
- numBasis*spaceDim*3 +
- numBasis*numBasis*spaceDim*2 +
- 196+84));
+ PetscLogFlops(cells->size()*(3 + numBasis*numBasis*spaceDim*2 + 196+84));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidual
@@ -445,8 +426,6 @@
"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);
@@ -466,28 +445,27 @@
assert(!sieve.isNull());
// Get sections
- const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
- topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection,
- numBasis*spaceDim,
- &dispTCell[0]);
- const ALE::Obj<RealSection>& dispTmdtSection =
- fields->get("disp(t-dt)").section();
- assert(!dispTmdtSection.isNull());
- topology::Mesh::RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
- numBasis*spaceDim,
- &dispTmdtCell[0]);
+ double_array accCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& accSection =
+ fields->get("acceleration(t)").section();
+ assert(!accSection.isNull());
+ RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
+
+ double_array dispCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispSection =
+ fields->get("disp(t)").section();
+ assert(!dispSection.isNull());
+ RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[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();
@@ -496,9 +474,6 @@
_normalizer->densityScale());
// Get parameters used in integration.
- const double dt = _dt;
- const double dt2 = dt*dt;
- assert(dt > 0);
double_array valuesIJ(numBasis);
_logger->eventEnd(setupEvent);
@@ -519,20 +494,20 @@
coordsVisitor.clear();
sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- dispTVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTVisitor);
+ accVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, accVisitor);
- dispTmdtVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTmdtVisitor);
+ dispVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispVisitor);
#else
coordsVisitor.clear();
sieve->orientedConeOpt(*c_iter, coordsVisitor, numBasis, spaceDim);
- dispTVisitor.clear();
- sieve->orientedConeOpt(*c_iter, dispTVisitor, numBasis, spaceDim);
+ accVisitor.clear();
+ sieve->orientedConeOpt(*c_iter, accVisitor, numBasis, spaceDim);
- dispTmdtVisitor.clear();
- sieve->orientedConeOpt(*c_iter, dispTmdtVisitor, numBasis, spaceDim);
+ dispVisitor.clear();
+ sieve->orientedConeOpt(*c_iter, dispVisitor, numBasis, spaceDim);
#endif
#if defined(DETAILED_EVENT_LOGGING)
@@ -589,11 +564,11 @@
} // if
// Compute action for inertial terms
- const double wtVertex = density[0] * volume / (4.0 * dt2);
- _cellVector += wtVertex * (dispTCell - dispTmdtCell);
+ const double wtVertex = density[0] * volume / 4.0;
+ _cellVector -= wtVertex * accCell;
#if defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(2 + numBasis*spaceDim*3);
+ PetscLogFlops(2 + numBasis*spaceDim*2);
_logger->eventEnd(computeEvent);
_logger->eventBegin(stressEvent);
#endif
@@ -633,21 +608,21 @@
const double d4 = (x1*y2+x0*(y1-y2)-x2*y1-(x1-x2)*y0) / scaleB;
assert(strainCell.size() == 6);
- strainCell[0] = b1 * dispTCell[0] + b2 * dispTCell[3] + b3 * dispTCell[6]
- + b4 * dispTCell[9];
- strainCell[1] = c3 * dispTCell[7] + c2 * dispTCell[4] + c4 * dispTCell[10]
- + c1 * dispTCell[1];
- strainCell[2] = d3 * dispTCell[8] + d2 * dispTCell[5] + d1 * dispTCell[2]
- + d4 * dispTCell[11];
- strainCell[3] = (c4 * dispTCell[9] + b3 * dispTCell[7] + c3 * dispTCell[6]
- + b2 * dispTCell[4] + c2 * dispTCell[3] + b4 * dispTCell[10] + b1
- * dispTCell[1] + c1 * dispTCell[0]) / 2.0;
- strainCell[4] = (c3 * dispTCell[8] + d3 * dispTCell[7] + c2 * dispTCell[5]
- + d2 * dispTCell[4] + c1 * dispTCell[2] + c4 * dispTCell[11] + d4
- * dispTCell[10] + d1 * dispTCell[1]) / 2.0;
- strainCell[5] = (d4 * dispTCell[9] + b3 * dispTCell[8] + d3 * dispTCell[6]
- + b2 * dispTCell[5] + d2 * dispTCell[3] + b1 * dispTCell[2] + b4
- * dispTCell[11] + d1 * dispTCell[0]) / 2.0;
+ strainCell[0] = b1 * dispCell[0] + b2 * dispCell[3] + b3 * dispCell[6]
+ + b4 * dispCell[9];
+ strainCell[1] = c3 * dispCell[7] + c2 * dispCell[4] + c4 * dispCell[10]
+ + c1 * dispCell[1];
+ strainCell[2] = d3 * dispCell[8] + d2 * dispCell[5] + d1 * dispCell[2]
+ + d4 * dispCell[11];
+ strainCell[3] = (c4 * dispCell[9] + b3 * dispCell[7] + c3 * dispCell[6]
+ + b2 * dispCell[4] + c2 * dispCell[3] + b4 * dispCell[10] + b1
+ * dispCell[1] + c1 * dispCell[0]) / 2.0;
+ strainCell[4] = (c3 * dispCell[8] + d3 * dispCell[7] + c2 * dispCell[5]
+ + d2 * dispCell[4] + c1 * dispCell[2] + c4 * dispCell[11] + d4
+ * dispCell[10] + d1 * dispCell[1]) / 2.0;
+ strainCell[5] = (d4 * dispCell[9] + b3 * dispCell[8] + d3 * dispCell[6]
+ + b2 * dispCell[5] + d2 * dispCell[3] + b1 * dispCell[2] + b4
+ * dispCell[11] + d1 * dispCell[0]) / 2.0;
const double_array& stressCell = _material->calcStress(strainCell, true);
@@ -688,7 +663,7 @@
} // for
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*(2 + numBasis*spaceDim*3));
+ PetscLogFlops(cells->size()*(2 + numBasis*spaceDim*2 + 196+84));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidualLumped
@@ -728,9 +703,6 @@
"contribution to Jacobian matrix for cells with " \
"different dimensions than the spatial dimension.");
- // Allocate vectors for cell data.
- double_array dispTCell(numBasis*spaceDim);
-
// Get cell information
const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
assert(!sieveMesh.isNull());
@@ -742,9 +714,10 @@
const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
// Get sections
- const ALE::Obj<RealSection>& dispTSection =
+ double_array dispCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispSection =
fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
+ assert(!dispSection.isNull());
// Get sparse matrix
const PetscMat jacobianMat = jacobian->matrix();
@@ -756,10 +729,10 @@
assert(dt > 0);
const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", dispTSection);
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", dispSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
- topology::Mesh::IndicesVisitor jacobianVisitor(*dispTSection, *globalOrder,
+ topology::Mesh::IndicesVisitor jacobianVisitor(*dispSection, *globalOrder,
(int) pow(sieveMesh->getSieve()->getMaxConeSize(),
sieveMesh->depth())*spaceDim);
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshBuilder.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshBuilder.cc 2010-03-31 16:52:14 UTC (rev 16474)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshBuilder.cc 2010-03-31 19:08:39 UTC (rev 16475)
@@ -177,16 +177,10 @@
logger.stagePush("MeshCoordinates");
ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, spaceDim,
&(*coordinates)[0]);
- logger.stagePop();
+ logger.stagePop(); // MeshCoordinates
- logger.stagePush("MeshReordering");
- ALE::Obj<ALE::Ordering<>::perm_type> reordering = new ALE::Ordering<>::perm_type(sieveMesh->comm(), sieveMesh->debug());
+ logger.stagePop(); // Mesh
- ALE::Ordering<>::calculateMeshReordering(sieveMesh, reordering);
- sieveMesh->relabel(*reordering);
- logger.stagePop();
- logger.stagePop();
-
sieveMesh->getFactory()->clear();
} // buildMesh
Modified: short/3D/PyLith/trunk/libsrc/problems/Explicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Explicit.cc 2010-03-31 16:52:14 UTC (rev 16474)
+++ short/3D/PyLith/trunk/libsrc/problems/Explicit.cc 2010-03-31 19:08:39 UTC (rev 16475)
@@ -47,6 +47,8 @@
// = (dispIncr(t+dt) - disp(t) + disp(t-dt)) / (dt*dt)
const double dt = _dt;
+ const double dt2 = dt*dt;
+ const double twodt = 2.0*dt;
topology::Field<topology::Mesh>& dispIncr = _fields->get("dispIncr(t->t+dt)");
const spatialdata::geocoords::CoordSys* cs = dispIncr.mesh().coordsys();
@@ -107,8 +109,8 @@
dispTmdtSection->restrictPoint(*v_iter, &dispTmdtVertex[0],
dispTmdtVertex.size());
- velVertex = (dispIncrVertex + dispTVertex - dispTmdtVertex) / (2.0 * dt);
- accVertex = (dispIncrVertex - dispTVertex + dispTmdtVertex) / (dt * dt);
+ velVertex = (dispIncrVertex + dispTVertex - dispTmdtVertex) / twodt;
+ accVertex = (dispIncrVertex - dispTVertex + dispTmdtVertex) / dt2;
assert(velSection->getFiberDimension(*v_iter) == spaceDim);
velSection->updatePoint(*v_iter, &velVertex[0]);
Modified: short/3D/PyLith/trunk/libsrc/topology/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Makefile.am 2010-03-31 16:52:14 UTC (rev 16474)
+++ short/3D/PyLith/trunk/libsrc/topology/Makefile.am 2010-03-31 19:08:39 UTC (rev 16475)
@@ -26,6 +26,7 @@
Mesh.icc \
MeshOps.hh \
RefineUniform.hh \
+ ReverseCuthillMcKee.hh \
SolutionFields.hh \
SubMesh.hh \
SubMesh.icc \
Added: short/3D/PyLith/trunk/libsrc/topology/ReverseCuthillMcKee.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/ReverseCuthillMcKee.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/ReverseCuthillMcKee.cc 2010-03-31 19:08:39 UTC (rev 16475)
@@ -0,0 +1,45 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "ReverseCuthillMcKee.hh" // implementation of class methods
+
+#include "pylith/topology/Mesh.hh" // USES Mesh
+
+#include <cassert> // USES assert()
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+
+// ----------------------------------------------------------------------
+// Set vertices and cells in mesh.
+void
+pylith::topology::ReverseCuthillMcKee::reorder(topology::Mesh* mesh)
+{ // reorder
+ assert(0 != mesh);
+
+ //logger.stagePush("MeshReordering");
+
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
+ assert(!sieveMesh.isNull());
+ ALE::Obj<ALE::Ordering<>::perm_type> reordering =
+ new ALE::Ordering<>::perm_type(sieveMesh->comm(), sieveMesh->debug());
+
+ ALE::Ordering<>::calculateMeshReordering(sieveMesh, reordering);
+ sieveMesh->relabel(*reordering);
+
+ //logger.stagePop();
+} // reorder
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/topology/ReverseCuthillMcKee.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/ReverseCuthillMcKee.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/ReverseCuthillMcKee.hh 2010-03-31 19:08:39 UTC (rev 16475)
@@ -0,0 +1,48 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/topology/ReverseCuthillMcKee.hh
+ *
+ * @brief Interface to PETSc reverse Cuthill-McKee reordering.
+ */
+
+#if !defined(pylith_topology_reversecuthillmckee_hh)
+#define pylith_topology_reversecuthillmckee_hh
+
+// Include directives ---------------------------------------------------
+#include "topologyfwd.hh" // forward declarations
+
+// ReverseCuthillMcKee --------------------------------------------------
+/// Interface to PETSc reverse Cuthill-McKee reordering.
+class pylith::topology::ReverseCuthillMcKee
+{ // ReverseCuthillMcKee
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /** Reorder vertices and cells of mesh using PETSc routines
+ * implementing reverse Cuthill-McKee algorithm.
+ *
+ * @param mesh PyLith finite-element mesh.
+ */
+ static
+ void reorder(topology::Mesh* mesh);
+
+}; // ReverseCuthillMcKee
+
+#endif // pylith_topology_reversecuthillmckee_hh
+
+
+// End of file
+
+
Modified: short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh 2010-03-31 16:52:14 UTC (rev 16474)
+++ short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh 2010-03-31 19:08:39 UTC (rev 16475)
@@ -42,6 +42,8 @@
class MeshRefiner;
class RefineUniform;
+ class ReverseCuthillMcKee;
+
} // topology
} // pylith
Modified: short/3D/PyLith/trunk/modulesrc/topology/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/Makefile.am 2010-03-31 16:52:14 UTC (rev 16474)
+++ short/3D/PyLith/trunk/modulesrc/topology/Makefile.am 2010-03-31 19:08:39 UTC (rev 16475)
@@ -28,7 +28,8 @@
SolutionFields.i \
Jacobian.i \
Distributor.i \
- RefineUniform.i
+ RefineUniform.i \
+ ReverseCuthillMcKee.i
swig_generated = \
topology_wrap.cxx \
Added: short/3D/PyLith/trunk/modulesrc/topology/ReverseCuthillMcKee.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/ReverseCuthillMcKee.i (rev 0)
+++ short/3D/PyLith/trunk/modulesrc/topology/ReverseCuthillMcKee.i 2010-03-31 19:08:39 UTC (rev 16475)
@@ -0,0 +1,43 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file modulesrc/topology/ReverseCuthillMcKee.hh
+ *
+ * @brief Python interface to C++ PyLith ReverseCuthillMcKee object.
+ */
+
+namespace pylith {
+ namespace topology {
+
+ // ReverseCuthillMcKee ----------------------------------------------
+ class ReverseCuthillMcKee
+ { // ReverseCuthillMcKee
+
+ // PUBLIC METHODS /////////////////////////////////////////////////
+ public :
+
+ /** Reorder vertices and cells of mesh using PETSc routines
+ * implementing reverse Cuthill-McKee algorithm.
+ *
+ * @param mesh PyLith finite-element mesh.
+ */
+ static
+ void reorder(topology::Mesh* mesh);
+
+ }; // ReverseCuthillMcKee
+
+ } // topology
+} // pylith
+
+
+// End of file
Modified: short/3D/PyLith/trunk/modulesrc/topology/topology.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/topology.i 2010-03-31 16:52:14 UTC (rev 16474)
+++ short/3D/PyLith/trunk/modulesrc/topology/topology.i 2010-03-31 19:08:39 UTC (rev 16475)
@@ -25,6 +25,7 @@
#include "pylith/topology/Jacobian.hh"
#include "pylith/topology/Distributor.hh"
#include "pylith/topology/RefineUniform.hh"
+#include "pylith/topology/ReverseCuthillMcKee.hh"
%}
%include "exception.i"
@@ -60,6 +61,7 @@
%include "Jacobian.i"
%include "Distributor.i"
%include "RefineUniform.i"
+%include "ReverseCuthillMcKee.i"
// Template instatiation
%template(MeshField) pylith::topology::Field<pylith::topology::Mesh>;
Modified: short/3D/PyLith/trunk/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/pylith/Makefile.am 2010-03-31 16:52:14 UTC (rev 16474)
+++ short/3D/PyLith/trunk/pylith/Makefile.am 2010-03-31 19:08:39 UTC (rev 16475)
@@ -146,6 +146,7 @@
topology/MeshImporterDist.py \
topology/MeshRefiner.py \
topology/RefineUniform.py \
+ topology/ReverseCuthillMcKee.py \
utils/__init__.py \
utils/CheckpointTimer.py \
utils/CppData.py \
Modified: short/3D/PyLith/trunk/pylith/topology/MeshImporter.py
===================================================================
--- short/3D/PyLith/trunk/pylith/topology/MeshImporter.py 2010-03-31 16:52:14 UTC (rev 16474)
+++ short/3D/PyLith/trunk/pylith/topology/MeshImporter.py 2010-03-31 19:08:39 UTC (rev 16475)
@@ -46,7 +46,7 @@
import pyre.inventory
- reorderMesh = pyre.inventory.bool("reorder_mesh", default=True)
+ reorderMesh = pyre.inventory.bool("reorder_mesh", default=False)
reorderMesh.meta['tip'] = "Reorder mesh using reverse Cuthill-McKee."
from pylith.meshio.MeshIOAscii import MeshIOAscii
@@ -93,7 +93,11 @@
if self.debug:
mesh.view("Finite-element mesh.")
- # :TODO: Reorder mesh
+ # Reorder mesh
+ if self.reorderMesh:
+ from pylith.topology.ReverseCuthillMcKee import ReverseCuthillMcKee
+ ordering = ReverseCuthillMcKee()
+ ordering.reorder(mesh)
# Adjust topology
self._debug.log(resourceUsageString())
@@ -128,6 +132,7 @@
self.reader = self.inventory.reader
self.distributor = self.inventory.distributor
self.refiner = self.inventory.refiner
+ self.reorderMesh = self.inventory.reorderMesh
return
Added: short/3D/PyLith/trunk/pylith/topology/ReverseCuthillMcKee.py
===================================================================
--- short/3D/PyLith/trunk/pylith/topology/ReverseCuthillMcKee.py (rev 0)
+++ short/3D/PyLith/trunk/pylith/topology/ReverseCuthillMcKee.py 2010-03-31 19:08:39 UTC (rev 16475)
@@ -0,0 +1,44 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/topology/ReverseCuthillMcKee.py
+##
+## @brief Python interface to PETSc reverse Cuthill-McKee reordering
+## of mesh cells and vertices.
+
+from topology import ReverseCuthillMcKee as ModuleReverseCuthillMcKee
+
+# ReverseCuthillMcKee class
+class ReverseCuthillMcKee(ModuleReverseCuthillMcKee):
+ """
+ Python interface to PETSc reverse Cuthill-McKee reordering of mesh
+ cells and vertices.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self):
+ """
+ Constructor.
+ """
+ return
+
+
+ def reorder(self, mesh):
+ """
+ Set communicator.
+ """
+ ModuleReverseCuthillMcKee.reorder(mesh)
+ return
+
+
+# End of file
More information about the CIG-COMMITS
mailing list