[cig-commits] r20683 - in short/3D/PyLith/trunk: libsrc/pylith/bc libsrc/pylith/problems libsrc/pylith/topology unittests/libtests/faults
knepley at geodynamics.org
knepley at geodynamics.org
Tue Sep 4 15:51:09 PDT 2012
Author: knepley
Date: 2012-09-04 15:51:09 -0700 (Tue, 04 Sep 2012)
New Revision: 20683
Modified:
short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc
short/3D/PyLith/trunk/libsrc/pylith/problems/Implicit.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc
short/3D/PyLith/trunk/libsrc/pylith/topology/Jacobian.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
Log:
Fixed bug in Field, added dmMesh() to Field, converted sections in Implicit.cc and Jacobian.cc
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc 2012-09-04 20:41:19 UTC (rev 20682)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc 2012-09-04 22:51:09 UTC (rev 20683)
@@ -264,6 +264,7 @@
PetscErrorCode err;
assert(fieldSectionP);
+ err = VecGetArray(localVec, &array);CHECK_PETSC_ERROR(err);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const SieveMesh::point_type p_bc = _points[iPoint];
@@ -282,11 +283,10 @@
fieldSection->updatePointAll(p_bc, &fieldVertex[0]);
PetscScalar *array;
- err = VecGetArray(localVec, &array);CHECK_PETSC_ERROR(err);
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
- array[_bcDOF[iDOF]] = parametersVertex[valueIndex+iDOF];
- err = VecRestoreArray(localVec, &array);CHECK_PETSC_ERROR(err);
+ array[_bcDOF[iDOF]+off] = parametersVertex[valueIndex+iDOF];
} // for
+ err = VecRestoreArray(localVec, &array);CHECK_PETSC_ERROR(err);
} // setField
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/Implicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/Implicit.cc 2012-09-04 20:41:19 UTC (rev 20682)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/Implicit.cc 2012-09-04 22:51:09 UTC (rev 20683)
@@ -56,36 +56,40 @@
const int spaceDim = cs->spaceDim();
// Get sections.
- scalar_array dispIncrVertex(spaceDim);
- const ALE::Obj<RealSection>& dispIncrSection = dispIncr.section();
- assert(!dispIncrSection.isNull());
+ PetscSection dispIncrSection = dispIncr.petscSection();
+ Vec dispIncrVec = dispIncr.localVector();
+ PetscScalar *dispIncrArray;
+ assert(dispIncrSection);assert(dispIncrVec);
- scalar_array velVertex(spaceDim);
- const ALE::Obj<RealSection>& velSection =
- _fields->get("velocity(t)").section();
- assert(!velSection.isNull());
+ PetscSection velSection = _fields->get("velocity(t)").petscSection();
+ Vec velVec = _fields->get("velocity(t)").localVector();
+ PetscScalar *velArray;
+ assert(velSection);assert(velVec);
// Get mesh vertices.
- const ALE::Obj<SieveMesh>& sieveMesh = dispIncr.mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- assert(!vertices.isNull());
- const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
-
- for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter) {
- dispIncrSection->restrictPoint(*v_iter, &dispIncrVertex[0],
- dispIncrVertex.size());
- velVertex = dispIncrVertex / dt;
-
- assert(velSection->getFiberDimension(*v_iter) == spaceDim);
- velSection->updatePointAll(*v_iter, &velVertex[0]);
+ DM dmMesh = dispIncr.dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+
+ assert(dmMesh);
+ err = DMComplexGetDepthStratum(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
- PetscLogFlops(vertices->size() * spaceDim);
-
+ err = VecRestoreArray(dispIncrVec, &dispIncrArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(velVec, &velArray);CHECK_PETSC_ERROR(err);
+ PetscLogFlops((vEnd - vStart) * spaceDim);
} // calcRateFields
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-09-04 20:41:19 UTC (rev 20682)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-09-04 22:51:09 UTC (rev 20683)
@@ -751,8 +751,8 @@
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);
+ err = DMGlobalToLocalBegin(_dm, _globalVec, INSERT_VALUES, _localVec);CHECK_PETSC_ERROR(err);
+ err = DMGlobalToLocalEnd(_dm, _globalVec, INSERT_VALUES, _localVec);CHECK_PETSC_ERROR(err);
}
logger.stagePop();
} // complete
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2012-09-04 20:41:19 UTC (rev 20682)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2012-09-04 22:51:09 UTC (rev 20683)
@@ -130,6 +130,12 @@
*/
const mesh_type& mesh(void) const;
+ /** Get DM associated with field.
+ *
+ * @returns DM
+ */
+ DM dmMesh(void) const;
+
/** Set label for field.
*
* @param value Label for field.
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc 2012-09-04 20:41:19 UTC (rev 20682)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc 2012-09-04 22:51:09 UTC (rev 20683)
@@ -67,6 +67,13 @@
return _mesh;
}
+template<typename mesh_type, typename section_type>
+inline
+DM
+pylith::topology::Field<mesh_type, section_type>::dmMesh(void) const {
+ return _dm;
+}
+
// Get label for field.
template<typename mesh_type, typename section_type>
inline
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Jacobian.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Jacobian.cc 2012-09-04 20:41:19 UTC (rev 20682)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Jacobian.cc 2012-09-04 22:51:09 UTC (rev 20683)
@@ -60,8 +60,7 @@
_matrix(0),
_valuesChanged(true)
{ // constructor
- const ALE::Obj<SubMesh::SieveMesh>& sieveMesh = field.mesh().sieveMesh();
- const ALE::Obj<SubMesh::RealSection>& fieldSection = field.section();
+ DM dmMesh = field.dmMesh();
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Jacobian");
@@ -69,8 +68,7 @@
// dimension, otherwise use a block size of 1.
const int blockFlag = (blockOkay) ? -1 : 1;
- PetscErrorCode err = DMMeshCreateMatrix(sieveMesh, fieldSection,
- matrixType, &_matrix, blockFlag);
+ PetscErrorCode err = DMCreateMatrix(dmMesh, matrixType, &_matrix);
CHECK_PETSC_ERROR_MSG(err, "Could not create PETSc sparse matrix "
"associated with subsystem Jacobian.");
logger.stagePop();
@@ -209,15 +207,16 @@
pylith::topology::Jacobian::verifySymmetry(void) const
{ // verifySymmetry
const PetscMat matSparse = _matrix;
+ PetscErrorCode err;
int nrows = 0;
int ncols = 0;
- MatGetSize(matSparse, &nrows, &ncols);
+ err = MatGetSize(matSparse, &nrows, &ncols);CHECK_PETSC_ERROR(err);
PetscMat matDense;
PetscMat matSparseAIJ;
- MatConvert(matSparse, MATSEQAIJ, MAT_INITIAL_MATRIX, &matSparseAIJ);
- MatConvert(matSparseAIJ, MATSEQDENSE, MAT_INITIAL_MATRIX, &matDense);
+ err = MatConvert(matSparse, MATSEQAIJ, MAT_INITIAL_MATRIX, &matSparseAIJ);CHECK_PETSC_ERROR(err);
+ err = MatConvert(matSparseAIJ, MATSEQDENSE, MAT_INITIAL_MATRIX, &matDense);CHECK_PETSC_ERROR(err);
scalar_array vals(nrows*ncols);
int_array rows(nrows);
@@ -226,7 +225,7 @@
rows[iRow] = iRow;
for (int iCol=0; iCol < ncols; ++iCol)
cols[iCol] = iCol;
- MatGetValues(matDense, nrows, &rows[0], ncols, &cols[0], &vals[0]);
+ err = MatGetValues(matDense, nrows, &rows[0], ncols, &cols[0], &vals[0]);CHECK_PETSC_ERROR(err);
const PylithScalar tolerance = 1.0e-06;
bool isSymmetric = true;
for (int iRow=0; iRow < nrows; ++iRow)
@@ -236,24 +235,24 @@
const PylithScalar valIJ = vals[indexIJ];
const PylithScalar valJI = vals[indexJI];
if (fabs(valIJ) > 1.0)
- if (fabs(1.0 - valJI/valIJ) > tolerance) {
- std::cerr << "Mismatch: "
- << "(" << iRow << ", " << iCol << ") = " << valIJ
- << ", (" << iCol << ", " << iRow << ") = " << valJI
- << std::endl;
- isSymmetric = false;
- } // if
+ if (fabs(1.0 - valJI/valIJ) > tolerance) {
+ std::cerr << "Mismatch: "
+ << "(" << iRow << ", " << iCol << ") = " << valIJ
+ << ", (" << iCol << ", " << iRow << ") = " << valJI
+ << std::endl;
+ isSymmetric = false;
+ } // if
else
- if (fabs(valJI - valIJ) > tolerance) {
- std::cerr << "Mismatch: "
- << "(" << iRow << ", " << iCol << ") = " << valIJ
- << ", (" << iCol << ", " << iRow << ") = " << valJI
- << std::endl;
- isSymmetric = false;
- } // if
+ if (fabs(valJI - valIJ) > tolerance) {
+ std::cerr << "Mismatch: "
+ << "(" << iRow << ", " << iCol << ") = " << valIJ
+ << ", (" << iCol << ", " << iRow << ") = " << valJI
+ << std::endl;
+ isSymmetric = false;
+ } // if
} // for
- MatDestroy(&matDense);
- MatDestroy(&matSparseAIJ);
+ err = MatDestroy(&matDense);CHECK_PETSC_ERROR(err);
+ err = MatDestroy(&matSparseAIJ);CHECK_PETSC_ERROR(err);
if (!isSymmetric)
throw std::runtime_error("Jacobian matrix is not symmetric.");
} // verifySymmetry
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2012-09-04 20:41:19 UTC (rev 20682)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2012-09-04 22:51:09 UTC (rev 20683)
@@ -728,7 +728,7 @@
meshio::MeshIOAscii iohandler;
iohandler.filename(_data->meshFilename);
iohandler.read(mesh);
-
+
//mesh->debug(true); // DEBUGGING
_quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
More information about the CIG-COMMITS
mailing list