[cig-commits] r20675 - in short/3D/PyLith/trunk: libsrc/pylith/feassemble libsrc/pylith/topology unittests/libtests/feassemble
knepley at geodynamics.org
knepley at geodynamics.org
Mon Sep 3 07:57:48 PDT 2012
Author: knepley
Date: 2012-09-03 07:57:47 -0700 (Mon, 03 Sep 2012)
New Revision: 20675
Modified:
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc
Log:
Initial test of new IntegrateResidual work, and put more operations in Field
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc 2012-09-03 03:54:59 UTC (rev 20674)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc 2012-09-03 14:57:47 UTC (rev 20675)
@@ -111,7 +111,7 @@
// ----------------------------------------------------------------------
// Integrate constributions to residual term (r) for operator.
void
-pylith::feassemble::ElasticityImplicit::integrateResidual(
+pylith::feassemble::ElasticityImplicit::integrateResidualOld(
const topology::Field<topology::Mesh>& residual,
const PylithScalar t,
topology::SolutionFields* const fields)
@@ -308,6 +308,203 @@
_logger->eventEnd(computeEvent);
} // integrateResidual
+void
+pylith::feassemble::ElasticityImplicit::integrateResidual(
+ const topology::Field<topology::Mesh>& residual,
+ const PylithScalar t,
+ topology::SolutionFields* const fields)
+{ // integrateResidual
+ /// Member prototype for _elasticityResidualXD()
+ typedef void (pylith::feassemble::ElasticityImplicit::*elasticityResidual_fn_type)
+ (const scalar_array&);
+ PetscErrorCode err;
+
+ assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(0 != _logger);
+ assert(0 != fields);
+
+ const int setupEvent = _logger->eventId("ElIR setup");
+ const int geometryEvent = _logger->eventId("ElIR geometry");
+ const int computeEvent = _logger->eventId("ElIR compute");
+ const int restrictEvent = _logger->eventId("ElIR restrict");
+ const int stateVarsEvent = _logger->eventId("ElIR stateVars");
+ const int stressEvent = _logger->eventId("ElIR stress");
+ const int updateEvent = _logger->eventId("ElIR update");
+
+ _logger->eventBegin(setupEvent);
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const scalar_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ const int tensorSize = _material->tensorSize();
+ if (cellDim != spaceDim)
+ throw std::logic_error("Integration for cells with spatial dimensions "
+ "different than the spatial dimension of the "
+ "domain not implemented yet.");
+
+ // Set variables dependent on dimension of cell
+ totalStrain_fn_type calcTotalStrainFn;
+ elasticityResidual_fn_type elasticityResidualFn;
+ if (1 == cellDim) {
+ elasticityResidualFn =
+ &pylith::feassemble::ElasticityImplicit::_elasticityResidual1D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D;
+ } else if (2 == cellDim) {
+ elasticityResidualFn =
+ &pylith::feassemble::ElasticityImplicit::_elasticityResidual2D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D;
+ } else if (3 == cellDim) {
+ elasticityResidualFn =
+ &pylith::feassemble::ElasticityImplicit::_elasticityResidual3D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain3D;
+ } else
+ assert(0);
+
+ // Allocate vectors for cell values.
+ scalar_array dispTpdtCell(numBasis*spaceDim);
+ scalar_array strainCell(numQuadPts*tensorSize);
+ strainCell = 0.0;
+ scalar_array gravVec(spaceDim);
+ scalar_array quadPtsGlobal(numQuadPts*spaceDim);
+
+ // Get cell information
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ assert(dmMesh);
+ const int materialId = _material->id();
+ err = DMComplexGetStratumIS(dmMesh, "material-id", materialId, &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
+ // Get sections
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
+
+ scalar_array dispTIncrCell(numBasis*spaceDim);
+ topology::Field<topology::Mesh>& dispTIncr = fields->get("dispIncr(t->t+dt)");
+ PetscSection dispTIncrSection = dispTIncr.petscSection();
+ Vec dispTIncrVec = dispTIncr.localVector();
+ assert(dispTIncrSection);assert(dispTIncrVec);
+
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+
+ scalar_array coordinatesCell(numBasis*spaceDim);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
+
+ assert(0 != _normalizer);
+ const PylithScalar lengthScale = _normalizer->lengthScale();
+ const PylithScalar gravityScale =
+ _normalizer->pressureScale() / (_normalizer->lengthScale() *
+ _normalizer->densityScale());
+
+ _logger->eventEnd(setupEvent);
+ _logger->eventBegin(computeEvent);
+
+ // 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
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, cell);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+#endif
+
+ // Get state variables for cell.
+ _material->retrievePropsAndVars(cell);
+
+ // Reset element vector to zero
+ _resetCellVector();
+
+ // Restrict input fields to cell
+ const PetscScalar *dispTArray, *dispTIncrArray;
+ PetscInt dispTSize, dispTIncrSize;
+ err = DMComplexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, dispTIncrSection, dispTIncrVec, cell, &dispTIncrSize, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ assert(dispTSize == dispTIncrSize);
+
+ // 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();
+ const scalar_array& quadPtsNondim = _quadrature->quadPts();
+
+ // Compute current estimate of displacement at time t+dt using solution increment.
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispTpdtCell[i] = dispTArray[i] + dispTIncrArray[i];}
+ err = DMComplexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, dispTIncrSection, dispTIncrVec, cell, &dispTIncrSize, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+
+ // Compute body force vector if gravity is being used.
+ if (0 != _gravityField) {
+ const spatialdata::geocoords::CoordSys* cs = fields->mesh().coordsys();
+ assert(0 != cs);
+
+ // Get density at quadrature points for this cell
+ const scalar_array& density = _material->calcDensity();
+
+ quadPtsGlobal = quadPtsNondim;
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(), lengthScale);
+
+ // Compute action for element body forces
+ spatialdata::spatialdb::SpatialDB* db = _gravityField;
+ for (int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
+ const int err = db->query(&gravVec[0], gravVec.size(), &quadPtsGlobal[0], spaceDim, cs);
+ if (err) throw std::runtime_error("Unable to get gravity vector for point.");
+ _normalizer->nondimensionalize(&gravVec[0], gravVec.size(), gravityScale);
+ const PylithScalar wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad];
+ for (int iBasis = 0, iQ = iQuad * numBasis; iBasis < numBasis; ++iBasis) {
+ const PylithScalar valI = wt * basis[iQ + iBasis];
+ for (int iDim = 0; iDim < spaceDim; ++iDim) {
+ _cellVector[iBasis * spaceDim + iDim] += valI * gravVec[iDim];
+ } // for
+ } // for
+ } // for
+ PetscLogFlops(numQuadPts * (2 + numBasis * (1 + 2 * spaceDim)));
+ }
+
+ // residualSection->view("After gravity contribution");
+ // Compute B(transpose) * sigma, first computing strains
+ calcTotalStrainFn(&strainCell, basisDeriv, dispTpdtCell, numBasis, numQuadPts);
+ const scalar_array& stressCell = _material->calcStress(strainCell, true);
+
+ CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell);
+
+#if 0 // DEBUGGING
+ std::cout << "Updating residual for cell " << cell << std::endl;
+ for(PetscInt i = 0; i < spaceDim*numBasis; ++i) {
+ std::cout << " v["<<i<<"]: " << _cellVector[i] << std::endl;
+ }
+#endif
+ // Assemble cell contribution into field
+ err = DMComplexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
+ }
+
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ _logger->eventEnd(computeEvent);
+} // integrateResidual
+
// ----------------------------------------------------------------------
// Compute stiffness matrix.
void
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.hh 2012-09-03 03:54:59 UTC (rev 20674)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.hh 2012-09-03 14:57:47 UTC (rev 20675)
@@ -112,6 +112,9 @@
void integrateResidual(const topology::Field<topology::Mesh>& residual,
const PylithScalar t,
topology::SolutionFields* const fields);
+ void integrateResidualOld(const topology::Field<topology::Mesh>& residual,
+ const PylithScalar t,
+ topology::SolutionFields* const fields);
/** Integrate contributions to Jacobian matrix (A) associated with
* operator.
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-09-03 03:54:59 UTC (rev 20674)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-09-03 14:57:47 UTC (rev 20675)
@@ -418,9 +418,9 @@
for (typename chart_type::const_iterator c_iter = chartBegin;
c_iter != chartEnd;
++c_iter) {
- const int fiberDim = srcSection->getFiberDimension(*c_iter);
- if (fiberDim > 0)
- _section->setFiberDimension(*c_iter, fiberDim);
+ const int fiberDim = srcSection->getFiberDimension(*c_iter);
+ if (fiberDim > 0)
+ _section->setFiberDimension(*c_iter, fiberDim);
} // for
const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
assert(!sieveMesh.isNull());
@@ -428,9 +428,19 @@
_section->setBC(srcSection->getBC());
_section->copySpaces(srcSection); // :BUG: NEED TO REBUILD SPACES
} // if/else
+ PetscSection section = src.petscSection();
+ PetscSection newSection;
+ PetscErrorCode err;
+
+ if (_dm) {
+ err = PetscSectionClone(section, &newSection);CHECK_PETSC_ERROR(err);
+ err = DMSetDefaultSection(_dm, newSection);CHECK_PETSC_ERROR(err);
+ err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
+ err = DMCreateLocalVector(_dm, &_localVec);CHECK_PETSC_ERROR(err);
+ err = PetscPrintf(PETSC_COMM_SELF, "Created new section %p, glboal vector %p, and local vector %p\n", newSection, _globalVec, _localVec);CHECK_PETSC_ERROR(err);
+ }
// Reuse scatters in clone
- PetscErrorCode err = 0;
if (src._scatters.size() > 0) {
const typename scatter_map_type::const_iterator scattersEnd = src._scatters.end();
for (typename scatter_map_type::const_iterator s_iter=src._scatters.begin();
@@ -663,6 +673,29 @@
{ // zero
if (!_section.isNull())
_section->zero(); // Does not zero BC.
+ if (_localVec) {
+ PetscSection section;
+ PetscInt pStart, pEnd, dof = 0;
+ PetscErrorCode err;
+
+ err = DMGetDefaultSection(_dm, §ion);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetChart(section, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+
+ // Assume fiber dimension is uniform
+ if (pEnd > pStart) {err = PetscSectionGetDof(section, pStart, &dof);CHECK_PETSC_ERROR(err);}
+ scalar_array values(dof);
+ values *= 0.0;
+
+ for(PetscInt p = pStart; p < pEnd; ++p) {
+ PetscInt pdof;
+
+ err = PetscSectionGetDof(section, p, &pdof);CHECK_PETSC_ERROR(err);
+ if (pdof > 0) {
+ assert(dof == pdof);
+ err = DMComplexVecSetClosure(_dm, section, _localVec, p, &values[0], INSERT_VALUES);CHECK_PETSC_ERROR(err);
+ } // if
+ } // for
+ }
} // zero
// ----------------------------------------------------------------------
@@ -691,6 +724,9 @@
} // if
} // for
} // if
+ if (_localVec) {
+ PetscErrorCode err = VecSet(_localVec, 0.0);CHECK_PETSC_ERROR(err);
+ }
} // zeroAll
// ----------------------------------------------------------------------
@@ -710,6 +746,15 @@
sieveMesh->getRecvOverlap(),
_section, _section);
+ if (_dm) {
+ // Not sure if DMLocalToLocal() would work
+ PetscErrorCode err;
+
+ err = DMLocalToGlobalBegin(_dm, _localVec, ADD_VALUES, _globalVec);CHECK_PETSC_ERROR(err);
+ err = DMLocalToGlobalEnd(_dm, _localVec, ADD_VALUES, _globalVec);CHECK_PETSC_ERROR(err);
+ err = DMGlobalToLocalBegin(_dm, _globalVec, ADD_VALUES, _localVec);CHECK_PETSC_ERROR(err);
+ err = DMGlobalToLocalEnd(_dm, _globalVec, ADD_VALUES, _localVec);CHECK_PETSC_ERROR(err);
+ }
logger.stagePop();
} // complete
@@ -1318,6 +1363,11 @@
INSERT_VALUES, SCATTER_FORWARD); CHECK_PETSC_ERROR(err);
err = VecScatterEnd(sinfo.scatter, sinfo.scatterVec, vector,
INSERT_VALUES, SCATTER_FORWARD); CHECK_PETSC_ERROR(err);
+
+ if (_dm) {
+ err = DMLocalToGlobalBegin(_dm, _localVec, ADD_VALUES, _globalVec);CHECK_PETSC_ERROR(err);
+ err = DMLocalToGlobalEnd(_dm, _localVec, ADD_VALUES, _globalVec);CHECK_PETSC_ERROR(err);
+ }
} // scatterSectionToVector
// ----------------------------------------------------------------------
@@ -1351,6 +1401,11 @@
INSERT_VALUES, SCATTER_REVERSE); CHECK_PETSC_ERROR(err);
err = VecScatterEnd(sinfo.scatter, vector, sinfo.scatterVec,
INSERT_VALUES, SCATTER_REVERSE); CHECK_PETSC_ERROR(err);
+
+ if (_dm) {
+ err = DMGlobalToLocalBegin(_dm, _globalVec, INSERT_VALUES, _localVec);CHECK_PETSC_ERROR(err);
+ err = DMGlobalToLocalEnd(_dm, _globalVec, INSERT_VALUES, _localVec);CHECK_PETSC_ERROR(err);
+ }
} // scatterVectorToSection
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc 2012-09-03 03:54:59 UTC (rev 20674)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc 2012-09-03 14:57:47 UTC (rev 20675)
@@ -173,10 +173,15 @@
const PylithScalar* valsE = _data->valsResidual;
const int sizeE = _data->spaceDim * _data->numVertices;
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
- const PylithScalar* vals = residualSection->restrictSpace();
- const int size = residualSection->sizeWithBC();
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(residualSection);
+ err = PetscSectionGetStorageSize(residualSection, &size);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
const PylithScalar tolerance = (sizeof(double) == sizeof(PylithScalar)) ? 1.0e-06 : 4.0e-05;
@@ -185,6 +190,7 @@
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
+ err = VecRestoreArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
} // testIntegrateResidual
// ----------------------------------------------------------------------
@@ -275,6 +281,9 @@
CPPUNIT_ASSERT_EQUAL(pylith::PYLITH_MAXSCALAR, stableTimeStep);
} // testStableTimeStep
+extern PetscErrorCode DMComplexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]);
+extern PetscErrorCode DMComplexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]);
+
// ----------------------------------------------------------------------
// Initialize elasticity integrator.
void
@@ -301,6 +310,10 @@
new SieveMesh::sieve_type(mesh->comm());
CPPUNIT_ASSERT(!sieve.isNull());
+ mesh->createDMMesh(_data->cellDim);
+ DM dmMesh = mesh->dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+
// Cells and vertices
const bool interpolate = false;
ALE::Obj<SieveFlexMesh::sieve_type> s =
@@ -317,7 +330,11 @@
sieveMesh->stratify();
ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, _data->spaceDim,
_data->vertices);
+ PetscErrorCode err;
+ err = DMComplexBuildFromCellList_Private(dmMesh, _data->numCells, _data->numVertices, _data->numBasis, _data->cells);CHECK_PETSC_ERROR(err);
+ err = DMComplexBuildCoordinates_Private(dmMesh, _data->spaceDim, _data->numCells, _data->numVertices, _data->vertices);CHECK_PETSC_ERROR(err);
+
// Material ids
const ALE::Obj<SieveMesh::label_sequence>& cells =
sieveMesh->heightStratum(0);
@@ -330,7 +347,13 @@
e_iter != cells->end();
++e_iter)
sieveMesh->setValue(labelMaterials, *e_iter, _data->matId);
+ PetscInt cStart, cEnd;
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ err = DMComplexSetLabelValue(dmMesh, "material-id", c, _data->matId);CHECK_PETSC_ERROR(err);
+ }
+
// Setup quadrature
_quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
_data->basisDerivRef, _data->numQuadPts,
@@ -369,6 +392,7 @@
residual.newSection(topology::FieldBase::VERTICES_FIELD, _data->spaceDim);
residual.allocate();
residual.zero();
+ std::cout << residual.label() << " section " << residual.petscSection() << std::endl;
fields->copyLayout("residual");
const int fieldSize = _data->spaceDim * _data->numVertices;
@@ -385,6 +409,18 @@
dispTIncrSection->updatePoint(iVertex+offset,
&_data->fieldTIncr[iVertex*_data->spaceDim]);
} // for
+
+ PetscSection dispTSectionP = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ CPPUNIT_ASSERT(dispTSectionP);CPPUNIT_ASSERT(dispTVec);
+ PetscSection dispTIncrSectionP = dispTIncr.petscSection();
+ Vec dispTIncrVec = dispTIncr.localVector();
+ CPPUNIT_ASSERT(dispTIncrSectionP);CPPUNIT_ASSERT(dispTIncrVec);
+ for(int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
+ err = DMComplexVecSetClosure(dmMesh, dispTSectionP, dispTVec, iVertex+offset, &_data->fieldT[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, dispTIncrSectionP, dispTIncrVec, iVertex+offset, &_data->fieldTIncr[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ } // for
+
} // _initialize
More information about the CIG-COMMITS
mailing list