[cig-commits] r21469 - short/3D/PyLith/trunk/libsrc/pylith/problems
brad at geodynamics.org
brad at geodynamics.org
Thu Mar 7 16:02:57 PST 2013
Author: brad
Date: 2013-03-07 16:02:57 -0800 (Thu, 07 Mar 2013)
New Revision: 21469
Modified:
short/3D/PyLith/trunk/libsrc/pylith/problems/Explicit.cc
short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.cc
short/3D/PyLith/trunk/libsrc/pylith/problems/Implicit.cc
Log:
Code cleanup.
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/Explicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/Explicit.cc 2013-03-07 23:39:15 UTC (rev 21468)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/Explicit.cc 2013-03-08 00:02:57 UTC (rev 21469)
@@ -23,10 +23,6 @@
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
// Constructor
pylith::problems::Explicit::Explicit(void)
{ // constructor
@@ -43,7 +39,7 @@
void
pylith::problems::Explicit::calcRateFields(void)
{ // calcRateFields
- assert(0 != _fields);
+ assert(_fields);
// vel(t) = (disp(t+dt) - disp(t-dt)) / (2*dt)
// = (dispIncr(t+dt) + disp(t) - disp(t-dt)) / (2*dt)
@@ -57,71 +53,64 @@
topology::Field<topology::Mesh>& dispIncr = _fields->get("dispIncr(t->t+dt)");
const spatialdata::geocoords::CoordSys* cs = dispIncr.mesh().coordsys();
- assert(0 != cs);
+ assert(cs);
const int spaceDim = cs->spaceDim();
+
+ PetscErrorCode err = 0;
// Get sections.
- PetscSection dispIncrSection = dispIncr.petscSection();
- Vec dispIncrVec = dispIncr.localVector();
- PetscScalar *dispIncrArray;
- assert(dispIncrSection);assert(dispIncrVec);
+ PetscSection dispIncrSection = dispIncr.petscSection();assert(dispIncrSection);
+ PetscVec dispIncrVec = dispIncr.localVector();assert(dispIncrVec);
+ PetscScalar *dispIncrArray = NULL;
+ err = VecGetArray(dispIncrVec, &dispIncrArray);CHECK_PETSC_ERROR(err);
- PetscSection dispTSection = _fields->get("disp(t)").petscSection();
- Vec dispTVec = _fields->get("disp(t)").localVector();
- PetscScalar *dispTArray;
- assert(dispTSection);assert(dispTVec);
+ PetscSection dispTSection = _fields->get("disp(t)").petscSection();assert(dispTSection);
+ PetscVec dispTVec = _fields->get("disp(t)").localVector();assert(dispTVec);
+ PetscScalar *dispTArray = NULL;
+ err = VecGetArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
- PetscSection dispTmdtSection = _fields->get("disp(t-dt)").petscSection();
- Vec dispTmdtVec = _fields->get("disp(t-dt)").localVector();
- PetscScalar *dispTmdtArray;
- assert(dispTmdtSection);assert(dispTmdtVec);
+ PetscSection dispTmdtSection = _fields->get("disp(t-dt)").petscSection();assert(dispTmdtSection);
+ PetscVec dispTmdtVec = _fields->get("disp(t-dt)").localVector();assert(dispTmdtVec);
+ PetscScalar *dispTmdtArray = NULL;
+ err = VecGetArray(dispTmdtVec, &dispTmdtArray);CHECK_PETSC_ERROR(err);
- PetscSection velSection = _fields->get("velocity(t)").petscSection();
- Vec velVec = _fields->get("velocity(t)").localVector();
- PetscScalar *velArray;
- assert(velSection);assert(velVec);
+ PetscSection velSection = _fields->get("velocity(t)").petscSection();assert(velSection);
+ PetscVec velVec = _fields->get("velocity(t)").localVector();assert(velVec);
+ PetscScalar *velArray = NULL;
+ err = VecGetArray(velVec, &velArray);CHECK_PETSC_ERROR(err);
- PetscSection accSection = _fields->get("acceleration(t)").petscSection();
- Vec accVec = _fields->get("acceleration(t)").localVector();
- PetscScalar *accArray;
- assert(accSection);assert(accVec);
+ PetscSection accSection = _fields->get("acceleration(t)").petscSection();assert(accSection);
+ PetscVec accVec = _fields->get("acceleration(t)").localVector();assert(accVec);
+ PetscScalar *accArray = NULL;
+ err = VecGetArray(accVec, &accArray);CHECK_PETSC_ERROR(err);
// Get mesh vertices.
- DM dmMesh = dispIncr.mesh().dmMesh();
- PetscInt vStart, vEnd;
- PetscErrorCode err;
-
- assert(dmMesh);
+ PetscDM dmMesh = dispIncr.mesh().dmMesh();assert(dmMesh);
+ PetscInt vStart, vEnd;
err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
- err = VecGetArray(dispIncrVec, &dispIncrArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(dispTmdtVec, &dispTmdtArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(velVec, &velArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(accVec, &accArray);CHECK_PETSC_ERROR(err);
for(PetscInt v = vStart; v < vEnd; ++v) {
PetscInt didof, dioff;
-
err = PetscSectionGetDof(dispIncrSection, v, &didof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(dispIncrSection, v, &dioff);CHECK_PETSC_ERROR(err);
assert(spaceDim == didof);
- PetscInt dtdof, dtoff;
+ PetscInt dtdof, dtoff;
err = PetscSectionGetDof(dispTSection, v, &dtdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(dispTSection, v, &dtoff);CHECK_PETSC_ERROR(err);
assert(spaceDim == dtdof);
- PetscInt dmdof, dmoff;
+ PetscInt dmdof, dmoff;
err = PetscSectionGetDof(dispTmdtSection, v, &dmdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(dispTmdtSection, v, &dmoff);CHECK_PETSC_ERROR(err);
assert(spaceDim == dmdof);
- PetscInt vdof, voff;
+ PetscInt vdof, voff;
err = PetscSectionGetDof(velSection, v, &vdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(velSection, v, &voff);CHECK_PETSC_ERROR(err);
assert(spaceDim == vdof);
- PetscInt adof, aoff;
+ PetscInt adof, aoff;
err = PetscSectionGetDof(accSection, v, &adof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(accSection, v, &aoff);CHECK_PETSC_ERROR(err);
assert(spaceDim == adof);
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.cc 2013-03-07 23:39:15 UTC (rev 21468)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.cc 2013-03-08 00:02:57 UTC (rev 21469)
@@ -142,8 +142,8 @@
pylith::problems::Formulation::meshIntegrators(IntegratorMesh** integrators,
const int numIntegrators)
{ // meshIntegrators
- assert( (0 == integrators && 0 == numIntegrators) ||
- (0 != integrators && 0 < numIntegrators) );
+ assert( (!integrators && 0 == numIntegrators) ||
+ (integrators && 0 < numIntegrators) );
_meshIntegrators.resize(numIntegrators);
for (int i=0; i < numIntegrators; ++i)
_meshIntegrators[i] = integrators[i];
@@ -155,8 +155,8 @@
pylith::problems::Formulation::submeshIntegrators(IntegratorSubMesh** integrators,
const int numIntegrators)
{ // submeshIntegrators
- assert( (0 == integrators && 0 == numIntegrators) ||
- (0 != integrators && 0 < numIntegrators) );
+ assert( (!integrators && 0 == numIntegrators) ||
+ (integrators && 0 < numIntegrators) );
_submeshIntegrators.resize(numIntegrators);
for (int i=0; i < numIntegrators; ++i)
_submeshIntegrators[i] = integrators[i];
@@ -184,8 +184,8 @@
const PylithScalar t,
const PylithScalar dt)
{ // updateSettings
- assert(0 != jacobian);
- assert(0 != fields);
+ assert(jacobian);
+ assert(fields);
assert(dt > 0.0);
_jacobian = jacobian;
@@ -203,8 +203,8 @@
const PylithScalar t,
const PylithScalar dt)
{ // updateSettings
- assert(0 != jacobian);
- assert(0 != fields);
+ assert(jacobian);
+ assert(fields);
assert(dt > 0.0);
_jacobianLumped = jacobian;
@@ -219,10 +219,10 @@
pylith::problems::Formulation::reformResidual(const PetscVec* tmpResidualVec,
const PetscVec* tmpSolutionVec)
{ // reformResidual
- assert(0 != _fields);
+ assert(_fields);
// Update section view of field.
- if (0 != tmpSolutionVec) {
+ if (tmpSolutionVec) {
topology::Field<topology::Mesh>& solution = _fields->solution();
solution.scatterVectorToSection(*tmpSolutionVec);
} // if
@@ -251,13 +251,13 @@
residual.complete();
// Update PETSc view of residual
- if (0 != tmpResidualVec)
+ if (tmpResidualVec)
residual.scatterSectionToVector(*tmpResidualVec);
else
residual.scatterSectionToVector();
// TODO: Move this to SolverLinear
- if (0 != tmpResidualVec)
+ if (tmpResidualVec)
VecScale(*tmpResidualVec, -1.0);
} // reformResidual
@@ -270,7 +270,7 @@
assert(0 != _fields);
// Update section view of field.
- if (0 != tmpSolutionVec) {
+ if (tmpSolutionVec) {
topology::Field<topology::Mesh>& solution = _fields->solution();
solution.scatterVectorToSection(*tmpSolutionVec);
} // if
@@ -315,8 +315,8 @@
void
pylith::problems::Formulation::reformJacobianLumped(void)
{ // reformJacobianLumped
- assert(0 != _jacobianLumped);
- assert(0 != _fields);
+ assert(_jacobianLumped);
+ assert(_fields);
// Set jacobian to zero.
_jacobianLumped->zero();
@@ -426,7 +426,6 @@
const int numTimeSteps = 1;
writer.open(mesh, numTimeSteps);
-
topology::Field<topology::Mesh>& solution = _fields->solution();
solution.scatterVectorToSection(*solutionVec);
writer.writeVertexField(0.0, solution, mesh);
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/Implicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/Implicit.cc 2013-03-07 23:39:15 UTC (rev 21468)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/Implicit.cc 2013-03-08 00:02:57 UTC (rev 21469)
@@ -51,43 +51,42 @@
topology::Field<topology::Mesh>& dispIncr = _fields->get("dispIncr(t->t+dt)");
const spatialdata::geocoords::CoordSys* cs = dispIncr.mesh().coordsys();
- assert(0 != cs);
+ assert(cs);
const int spaceDim = cs->spaceDim();
+
+ PetscErrorCode err = 0;
// Get sections.
- PetscSection dispIncrSection = dispIncr.petscSection();
- Vec dispIncrVec = dispIncr.localVector();
- PetscScalar *dispIncrArray;
- assert(dispIncrSection);assert(dispIncrVec);
-
- PetscSection velSection = _fields->get("velocity(t)").petscSection();
- Vec velVec = _fields->get("velocity(t)").localVector();
- PetscScalar *velArray;
- assert(velSection);assert(velVec);
+ PetscSection dispIncrSection = dispIncr.petscSection();assert(dispIncrSection);
+ PetscVec dispIncrVec = dispIncr.localVector();assert(dispIncrVec);
+ PetscScalar *dispIncrArray = NULL;
+ err = VecGetArray(dispIncrVec, &dispIncrArray);CHECK_PETSC_ERROR(err);
+
+ PetscSection velSection = _fields->get("velocity(t)").petscSection();assert(velSection);
+ PetscVec velVec = _fields->get("velocity(t)").localVector();assert(velVec);
+ PetscScalar *velArray = NULL;
+ err = VecGetArray(velVec, &velArray);CHECK_PETSC_ERROR(err);
// Get mesh vertices.
- DM dmMesh = dispIncr.dmMesh();
- PetscInt vStart, vEnd;
- PetscErrorCode err;
-
- assert(dmMesh);
+ PetscDM dmMesh = dispIncr.dmMesh();assert(dmMesh);
+ PetscInt vStart, vEnd;
err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
- err = VecGetArray(dispIncrVec, &dispIncrArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(velVec, &velArray);CHECK_PETSC_ERROR(err);
+
for(PetscInt v = vStart; v < vEnd; ++v) {
PetscInt dof, off, vdof, voff;
-
err = PetscSectionGetDof(dispIncrSection, v, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(dispIncrSection, v, &off);CHECK_PETSC_ERROR(err);
err = PetscSectionGetDof(velSection, v, &vdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(velSection, v, &voff);CHECK_PETSC_ERROR(err);
assert(dof == spaceDim);assert(vdof == spaceDim);
+
for(PetscInt d = 0; d < dof; ++d) {
velArray[voff+d] = dispIncrArray[off+d] / dt;
}
} // for
err = VecRestoreArray(dispIncrVec, &dispIncrArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(velVec, &velArray);CHECK_PETSC_ERROR(err);
+
PetscLogFlops((vEnd - vStart) * spaceDim);
} // calcRateFields
More information about the CIG-COMMITS
mailing list