[cig-commits] r14805 - in short/3D/PyLith/branches/pylith-swig: libsrc/faults libsrc/utils unittests/libtests/faults
brad at geodynamics.org
brad at geodynamics.org
Mon Apr 27 14:28:51 PDT 2009
Author: brad
Date: 2009-04-27 14:28:50 -0700 (Mon, 27 Apr 2009)
New Revision: 14805
Modified:
short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/branches/pylith-swig/libsrc/utils/sievetypes.hh
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/Makefile.am
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKin.hh
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinHex8.hh
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinLine2.cc
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinLine2.hh
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinQuad4.hh
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinQuad4e.hh
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinTri3.cc
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinTri3.hh
Log:
Worked on updating fault unit tests.
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.cc 2009-04-27 21:01:09 UTC (rev 14804)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.cc 2009-04-27 21:28:50 UTC (rev 14805)
@@ -121,7 +121,7 @@
const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
assert(!faultSieveMesh.isNull());
- const ALE::Obj<sieve_type>& sieve = sieveMesh->getSieve();
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
assert(!sieve.isNull());
const ALE::Obj<SieveSubMesh::sieve_type> ifaultSieve =
faultSieveMesh->getSieve();
@@ -229,8 +229,8 @@
TopologyOps::PointSet replaceCells;
TopologyOps::PointSet noReplaceCells;
TopologyOps::PointSet replaceVertices;
- ALE::ISieveVisitor::PointRetriever<sieve_type> sV2(std::max(1, ifaultSieve->getMaxSupportSize()));
- ALE::ISieveVisitor::NConeRetriever<sieve_type> cV2(*ifaultSieve, (size_t) pow(std::max(1, ifaultSieve->getMaxConeSize()), faultSieveMesh->depth()));
+ ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> sV2(std::max(1, ifaultSieve->getMaxSupportSize()));
+ ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type> cV2(*ifaultSieve, (size_t) pow(std::max(1, ifaultSieve->getMaxConeSize()), faultSieveMesh->depth()));
std::set<Mesh::point_type> faceSet;
const SieveSubMesh::label_sequence::const_iterator facesEnd = faces->end();
@@ -280,7 +280,7 @@
for(int c = 0; c < coneSize2; ++c)
std::cout << " " << cellCone[c] << std::endl;
std::cout << " fault cell support:" << std::endl;
- ALE::ISieveVisitor::PointRetriever<sieve_type> sV(std::max(1, ifaultSieve->getMaxSupportSize()));
+ ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> sV(std::max(1, ifaultSieve->getMaxSupportSize()));
ifaultSieve->support(face, sV);
const int supportSize2 = sV.getSize();
const point_type *cellSupport = sV.getPoints();
@@ -530,7 +530,7 @@
} // if/else
} // if/else
} // for
- ReplaceVisitor<sieve_type,std::map<Mesh::point_type,Mesh::point_type> > rVc(vertexRenumber, std::max(1, sieve->getMaxConeSize()), debug);
+ ReplaceVisitor<SieveMesh::sieve_type,std::map<Mesh::point_type,Mesh::point_type> > rVc(vertexRenumber, std::max(1, sieve->getMaxConeSize()), debug);
rCellsEnd = replaceCells.end();
for (TopologyOps::PointSet::const_iterator c_iter = replaceCells.begin();
@@ -544,7 +544,7 @@
} // if
rVc.clear();
} // for
- ReplaceVisitor<sieve_type,std::map<Mesh::point_type,Mesh::point_type> > rVs(cellRenumber, std::max(1, sieve->getMaxSupportSize()), debug);
+ ReplaceVisitor<SieveMesh::sieve_type,std::map<Mesh::point_type,Mesh::point_type> > rVs(cellRenumber, std::max(1, sieve->getMaxSupportSize()), debug);
rVerticesEnd = replaceVertices.end();
for (TopologyOps::PointSet::const_iterator v_iter = replaceVertices.begin();
@@ -642,9 +642,9 @@
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
assert(!sieveMesh.isNull());
ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
+ faultSieveMesh.destroy();
- const ALE::Obj<sieve_type>& sieve = sieveMesh->getSieve();
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
assert(!sieve.isNull());
faultSieveMesh =
new SieveSubMesh(mesh.comm(), mesh.dimension()-1, mesh.debug());
@@ -675,7 +675,7 @@
MPI_INT, MPI_SUM, sieve->comm());
int face = globalSieveEnd + globalFaceOffset - numFaces;
- ALE::ISieveVisitor::PointRetriever<sieve_type> cV(sieve->getMaxConeSize());
+ ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> cV(sieve->getMaxConeSize());
for(SieveMesh::label_sequence::iterator c_iter = cBegin;
c_iter != cEnd;
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc 2009-04-27 21:01:09 UTC (rev 14804)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc 2009-04-27 21:28:50 UTC (rev 14805)
@@ -90,6 +90,7 @@
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
assert(0 != cs);
+ delete _faultMesh; _faultMesh = new topology::SubMesh();
CohesiveTopology::createFaultParallel(_faultMesh, &_cohesiveToFault,
mesh, id(), _useLagrangeConstraints());
//_faultMesh->getLabel("height")->view("Fault mesh height");
@@ -379,10 +380,12 @@
const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
SieveSubMesh::renumbering_type& renumbering = faultSieveMesh->getRenumbering();
+ const SieveSubMesh::renumbering_type::const_iterator renumberingEnd =
+ renumbering.end();
for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
++v_iter)
- if (renumbering.find(*v_iter) != renumbering.end()) {
+ if (renumbering.find(*v_iter) != renumberingEnd) {
const int vertexFault = renumbering[*v_iter];
const int vertexMesh = *v_iter;
const double* slipVertex = slipSection->restrictPoint(vertexFault);
@@ -701,7 +704,7 @@
} else if (0 == strcasecmp("traction_change", name)) {
*fieldType = VECTOR_FIELD;
- const ALE::Obj<real_section_type>& solution = fields->getSolution();
+ const ALE::Obj<RealSection>& solution = fields->getSolution();
_calcTractionsChange(&_bufferVertexVector, mesh, solution);
_bufferTmp = _bufferVertexVector;
scale = _normalizer->pressureScale();
@@ -831,7 +834,7 @@
assert(!sieve.isNull());
typedef ALE::SieveAlg<SieveSubMesh> SieveAlg;
- ALE::ISieveVisitor::NConeRetriever<sieve_type> ncV(*sieve, (size_t) pow(sieve->getMaxConeSize(), std::max(0, faultSieveMesh->depth())));
+ ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type> ncV(*sieve, (size_t) pow(sieve->getMaxConeSize(), std::max(0, faultSieveMesh->depth())));
for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
c_iter != cellsEnd;
@@ -840,17 +843,34 @@
coordinatesVisitor.clear();
faultSieveMesh->restrictClosure(*c_iter, coordinatesVisitor);
+ ncV.clear();
ALE::ISieveTraversal<SieveSubMesh::sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
const int coneSize = ncV.getSize();
const Mesh::point_type *cone = ncV.getPoints();
- for(int v = 0; v < coneSize; ++v) {
+ for (int v=0; v < coneSize; ++v) {
// Compute Jacobian and determinant of Jacobian at vertex
memcpy(&refCoordsVertex[0], &verticesRef[v*cohesiveDim],
cohesiveDim*sizeof(double));
cellGeometry.jacobian(&jacobian, &jacobianDet, coordinatesCell,
refCoordsVertex);
+ for (int ii=0; ii < numBasis; ++ii) {
+ std::cout << " vertex " << ii << ": ";
+ for (int jj=0; jj < spaceDim; ++jj)
+ std::cout << " " << coordinatesCell[ii*spaceDim+jj];
+ std::cout << std::endl;
+ } // for
+ std::cout << " location vertex: ";
+ for (int jj=0; jj < cohesiveDim; ++jj)
+ std::cout << " " << refCoordsVertex[jj];
+ std::cout << std::endl;
+ std::cout << " jacobian: ";
+ for (int jj=0; jj < jacobianSize; ++jj)
+ std::cout << " " << jacobian[jj];
+ std::cout << std::endl;
+
+
// Compute orientation
cellGeometry.orientation(&orientationVertex, jacobian, jacobianDet,
upDirArray);
@@ -858,9 +878,10 @@
// Update orientation
orientationSection->updateAddPoint(cone[v], &orientationVertex[0]);
} // for
- ncV.clear();
} // for
+ orientation.view("ORIENTATION BEFORE COMPLETE");
+
// Assemble orientation information
orientation.complete();
@@ -922,7 +943,7 @@
PetscLogFlops(5 + count * 3);
} // if
- //_orientation->view("ORIENTATION");
+ orientation.view("ORIENTATION");
} // _calcOrientation
// ----------------------------------------------------------------------
@@ -947,8 +968,8 @@
const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
// Allocate stiffness field.
- _fields->add("stiffness", "stiffness");
- topology::Field<topology::SubMesh>& stiffness = _fields->get("stiffness");
+ _fields->add("pseudostiffness", "pseudostiffness");
+ topology::Field<topology::SubMesh>& stiffness = _fields->get("pseudostiffness");
stiffness.newSection(topology::FieldBase::VERTICES_FIELD, 1);
stiffness.allocate();
stiffness.zero();
@@ -1001,6 +1022,8 @@
PetscLogFlops(count * 2);
matDB->close();
+
+ stiffness.view("PSEUDO STIFFNESS");
} // _calcConditioning
// ----------------------------------------------------------------------
@@ -1123,14 +1146,16 @@
const int numFaultVertices = fvertices->size();
Mesh::renumbering_type& renumbering = faultSieveMesh->getRenumbering();
+ const SieveSubMesh::renumbering_type::const_iterator renumberingEnd =
+ renumbering.end();
#if 0 // MOVE TO SEPARATE CHECK METHOD
// Check fault mesh and volume mesh coordinates
- const ALE::Obj<real_section_type>& coordinates = mesh->getRealSection("coordinates");
- const ALE::Obj<real_section_type>& fCoordinates = _faultMesh->getRealSection("coordinates");
+ const ALE::Obj<RealSection>& coordinates = mesh->getRealSection("coordinates");
+ const ALE::Obj<RealSection>& fCoordinates = _faultMesh->getRealSection("coordinates");
for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != verticesEnd; ++v_iter) {
- if (renumbering.find(*v_iter) != renumbering.end()) {
+ if (renumbering.find(*v_iter) != renumberingEnd) {
const int v = *v_iter;
const int dim = coordinates->getFiberDimension(*v_iter);
const double *a = coordinates->restrictPoint(*v_iter);
@@ -1162,7 +1187,7 @@
for (SieveMesh::label_sequence::iterator v_iter = verticesBegin;
v_iter != verticesEnd;
++v_iter)
- if (renumbering.find(*v_iter) != renumbering.end()) {
+ if (renumbering.find(*v_iter) != renumberingEnd) {
const int vertexMesh = *v_iter;
const int vertexFault = renumbering[*v_iter];
assert(fiberDim == solutionSection->getFiberDimension(vertexMesh));
@@ -1170,14 +1195,11 @@
assert(1 == stiffnessSection->getFiberDimension(vertexFault));
assert(1 == areaSection->getFiberDimension(vertexFault));
- const real_section_type::value_type* solutionVertex =
- solutionSection->restrictPoint(vertexMesh);
+ const double* solutionVertex = solutionSection->restrictPoint(vertexMesh);
assert(0 != solutionVertex);
- const real_section_type::value_type* stiffnessVertex =
- stiffnessSection->restrictPoint(vertexFault);
+ const double* stiffnessVertex = stiffnessSection->restrictPoint(vertexFault);
assert(0 != stiffnessVertex);
- const real_section_type::value_type* areaVertex =
- areaSection->restrictPoint(vertexFault);
+ const double* areaVertex = areaSection->restrictPoint(vertexFault);
assert(0 != areaVertex);
const double scale = stiffnessVertex[0] / areaVertex[0];
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/utils/sievetypes.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/utils/sievetypes.hh 2009-04-27 21:01:09 UTC (rev 14804)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/utils/sievetypes.hh 2009-04-27 21:28:50 UTC (rev 14805)
@@ -25,9 +25,6 @@
typedef ALE::IMesh<> Mesh;
typedef ALE::IMesh<ALE::LabelSifter<int, Mesh::point_type> > SubMesh;
- typedef Mesh::sieve_type sieve_type;
- typedef Mesh::real_section_type real_section_type;
- typedef Mesh::int_section_type int_section_type;
} // pylith
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/Makefile.am 2009-04-27 21:01:09 UTC (rev 14804)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/Makefile.am 2009-04-27 21:28:50 UTC (rev 14805)
@@ -29,10 +29,10 @@
TestLiuCosSlipFn.cc \
TestEqKinSrc.cc \
TestFaultCohesiveKin.cc \
+ TestFaultCohesiveKinLine2.cc \
+ TestFaultCohesiveKinTri3.cc \
test_faults.cc
-# TestFaultCohesiveKinLine2.cc \
-# TestFaultCohesiveKinTri3.cc \
# TestFaultCohesiveKinTri3d.cc \
# TestFaultCohesiveKinQuad4.cc \
# TestFaultCohesiveKinQuad4e.cc \
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKin.cc 2009-04-27 21:01:09 UTC (rev 14804)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKin.cc 2009-04-27 21:28:50 UTC (rev 14805)
@@ -24,6 +24,7 @@
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+#include "pylith/topology/Jacobian.hh" // USES Jacobian
#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
@@ -46,7 +47,8 @@
pylith::faults::TestFaultCohesiveKin::setUp(void)
{ // setUp
_data = 0;
- _quadrature = 0;
+ _quadrature = new feassemble::Quadrature<topology::SubMesh>();
+ CPPUNIT_ASSERT(0 != _quadrature);
const int nsrcs = 1;
_eqsrcs.resize(nsrcs);
_eqsrcs[0] = new EqKinSrc();
@@ -137,9 +139,11 @@
{ // testInitialize
topology::Mesh mesh;
FaultCohesiveKin fault;
- _initialize(&mesh, &fault);
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &fault, &fields);
- //mesh.view(_data->meshFilename);
+ mesh.view(_data->meshFilename);
+ fault._faultMesh->view("FAULT MESH");
const ALE::Obj<SieveSubMesh>& faultSieveMesh = fault._faultMesh->sieveMesh();
CPPUNIT_ASSERT(!faultSieveMesh.isNull());
@@ -205,7 +209,7 @@
// Check pseudoStiffness
const ALE::Obj<RealSection>& stiffnessSection =
- fault._fields->get("stiffness").section();
+ fault._fields->get("pseudostiffness").section();
CPPUNIT_ASSERT(!stiffnessSection.isNull());
iVertex = 0;
for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
@@ -228,24 +232,13 @@
{ // testIntegrateResidual
topology::Mesh mesh;
FaultCohesiveKin fault;
- _initialize(&mesh, &fault);
-
- spatialdata::geocoords::CSCart cs;
- cs.setSpaceDim(mesh.dimension());
- cs.initialize();
-
- // Setup fields
topology::SolutionFields fields(mesh);
- fields.add("residual", "residual");
- fields.add("solution", "displacement");
- fields.solutionName("solution");
+ _initialize(&mesh, &fault, &fields);
const int spaceDim = _data->spaceDim;
topology::Field<topology::Mesh>& residual = fields.get("residual");
- residual.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
- residual.allocate();
- fields.copyLayout("residual");
const ALE::Obj<RealSection>& residualSection = residual.section();
+ CPPUNIT_ASSERT(!residualSection.isNull());
const ALE::Obj<RealSection>& solutionSection =
fields.get("solution").section();
@@ -350,68 +343,52 @@
void
pylith::faults::TestFaultCohesiveKin::testIntegrateJacobian(void)
{ // testIntegrateJacobian
-#if 0
- ALE::Obj<Mesh> mesh;
+ topology::Mesh mesh;
FaultCohesiveKin fault;
- _initialize(&mesh, &fault);
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &fault, &fields);
- // Setup fields
- topology::FieldsManager fields(mesh);
- fields.addReal("residual");
- fields.addReal("solution");
- fields.solutionField("solution");
-
- const ALE::Obj<real_section_type>& residual = fields.getReal("residual");
- CPPUNIT_ASSERT(!residual.isNull());
- const int spaceDim = _data->spaceDim;
- residual->setChart(mesh->getSieve()->getChart());
- residual->setFiberDimension(mesh->depthStratum(0), spaceDim);
- mesh->allocate(residual);
- residual->zero();
- fields.copyLayout("residual");
+ const ALE::Obj<RealSection>& solutionSection = fields.get("solution").section();
+ CPPUNIT_ASSERT(!solutionSection.isNull());
- const ALE::Obj<real_section_type>& solution = fields.getReal("solution");
- CPPUNIT_ASSERT(!solution.isNull());
-
- const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
+ const int spaceDim = _data->spaceDim;
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices = sieveMesh->depthStratum(0);
CPPUNIT_ASSERT(!vertices.isNull());
- const Mesh::label_sequence::iterator verticesBegin = vertices->begin();
- const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+ const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
int iVertex = 0;
- for (Mesh::label_sequence::iterator v_iter=verticesBegin;
+ for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
++v_iter, ++iVertex) {
- solution->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+ solutionSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
} // for
- PetscMat jacobian;
- PetscErrorCode err = MeshCreateMatrix(mesh, solution, MATMPIBAIJ, &jacobian);
- CPPUNIT_ASSERT(0 == err);
+ topology::Jacobian jacobian(fields);
const double t = 2.134;
- fault.integrateJacobian(&jacobian, t, &fields, mesh);
+ fault.integrateJacobian(&jacobian, t, &fields);
CPPUNIT_ASSERT_EQUAL(false, fault.needNewJacobian());
- err = MatAssemblyBegin(jacobian, MAT_FINAL_ASSEMBLY);
- CPPUNIT_ASSERT(0 == err);
- err = MatAssemblyEnd(jacobian, MAT_FINAL_ASSEMBLY);
- CPPUNIT_ASSERT(0 == err);
+ jacobian.assemble("final_assembly");
//MatView(jacobian, PETSC_VIEWER_STDOUT_WORLD); // DEBUGGING
const double* valsE = _data->valsJacobian;
- const int nrowsE = solution->sizeWithBC();
+ const int nrowsE = solutionSection->sizeWithBC();
const int ncolsE = nrowsE;
int nrows = 0;
int ncols = 0;
- MatGetSize(jacobian, &nrows, &ncols);
+ PetscMat jacobianMat = jacobian.matrix();
+ MatGetSize(jacobianMat, &nrows, &ncols);
CPPUNIT_ASSERT_EQUAL(nrowsE, nrows);
CPPUNIT_ASSERT_EQUAL(ncolsE, ncols);
PetscMat jDense;
PetscMat jSparseAIJ;
- MatConvert(jacobian, MATSEQAIJ, MAT_INITIAL_MATRIX, &jSparseAIJ);
+ MatConvert(jacobianMat, MATSEQAIJ, MAT_INITIAL_MATRIX, &jSparseAIJ);
MatConvert(jSparseAIJ, MATSEQDENSE, MAT_INITIAL_MATRIX, &jDense);
double_array vals(nrows*ncols);
@@ -439,7 +416,6 @@
MatDestroy(jDense);
MatDestroy(jSparseAIJ);
CPPUNIT_ASSERT_EQUAL(false, fault.needNewJacobian());
-#endif
} // testIntegrateJacobian
// ----------------------------------------------------------------------
@@ -447,50 +423,37 @@
void
pylith::faults::TestFaultCohesiveKin::testIntegrateResidualAssembled(void)
{ // testIntegrateResidualAssembled
-#if 0
- ALE::Obj<Mesh> mesh;
+ topology::Mesh mesh;
FaultCohesiveKin fault;
- _initialize(&mesh, &fault);
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &fault, &fields);
- spatialdata::geocoords::CSCart cs;
- cs.setSpaceDim((mesh)->getDimension());
- cs.initialize();
-
- // Setup fields
- topology::FieldsManager fields(mesh);
- fields.addReal("residual");
- fields.addReal("solution");
- fields.solutionField("solution");
-
- const ALE::Obj<real_section_type>& residual = fields.getReal("residual");
- CPPUNIT_ASSERT(!residual.isNull());
const int spaceDim = _data->spaceDim;
- residual->setChart(mesh->getSieve()->getChart());
- residual->setFiberDimension(mesh->depthStratum(0), spaceDim);
- mesh->allocate(residual);
- residual->zero();
- fields.copyLayout("residual");
+ topology::Field<topology::Mesh>& residual = fields.get("residual");
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ CPPUNIT_ASSERT(!residualSection.isNull());
- const ALE::Obj<real_section_type>& solution = fields.getReal("solution");
- CPPUNIT_ASSERT(!solution.isNull());
+ const ALE::Obj<RealSection>& solutionSection = fields.get("solution").section();
+ CPPUNIT_ASSERT(!solutionSection.isNull());
- const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices = sieveMesh->depthStratum(0);
CPPUNIT_ASSERT(!vertices.isNull());
- const Mesh::label_sequence::iterator verticesBegin = vertices->begin();
- const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+ const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
int iVertex = 0;
- for (Mesh::label_sequence::iterator v_iter=verticesBegin;
+ for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
- ++v_iter, ++iVertex) {
- solution->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
- } // for
+ ++v_iter, ++iVertex)
+ solutionSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
const double t = 2.134;
const double dt = 0.01;
fault.timeStep(dt);
{ // Integrate residual with solution (as opposed to solution increment).
fault.useSolnIncr(false);
- fault.integrateResidualAssembled(residual, t, &fields, mesh, &cs);
+ fault.integrateResidualAssembled(residual, t, &fields);
//residual->view("RESIDUAL"); // DEBUGGING
@@ -499,13 +462,12 @@
iVertex = 0;
const int fiberDimE = spaceDim;
const double tolerance = 1.0e-06;
- for (Mesh::label_sequence::iterator v_iter=verticesBegin;
+ for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
++v_iter, ++iVertex) {
- const int fiberDim = residual->getFiberDimension(*v_iter);
+ const int fiberDim = residualSection->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
- const real_section_type::value_type* vals =
- residual->restrictPoint(*v_iter);
+ const double* vals = residualSection->restrictPoint(*v_iter);
CPPUNIT_ASSERT(0 != vals);
const bool isConstraint = _isConstraintVertex(*v_iter);
@@ -527,10 +489,10 @@
} // for
} // Integrate residual with solution (as opposed to solution increment).
- residual->zero();
+ residual.zero();
{ // Integrate residual with solution increment.
fault.useSolnIncr(true);
- fault.integrateResidualAssembled(residual, t, &fields, mesh, &cs);
+ fault.integrateResidualAssembled(residual, t, &fields);
//residual->view("RESIDUAL"); // DEBUGGING
@@ -539,13 +501,12 @@
iVertex = 0;
const int fiberDimE = spaceDim;
const double tolerance = 1.0e-06;
- for (Mesh::label_sequence::iterator v_iter=verticesBegin;
+ for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
++v_iter, ++iVertex) {
- const int fiberDim = residual->getFiberDimension(*v_iter);
+ const int fiberDim = residualSection->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
- const real_section_type::value_type* vals =
- residual->restrictPoint(*v_iter);
+ const double* vals = residualSection->restrictPoint(*v_iter);
CPPUNIT_ASSERT(0 != vals);
const bool isConstraint = _isConstraintVertex(*v_iter);
@@ -566,7 +527,6 @@
} // if/else
} // for
} // Integrate residual with solution increment.
-#endif
} // testIntegrateResidualAssembled
// ----------------------------------------------------------------------
@@ -574,68 +534,52 @@
void
pylith::faults::TestFaultCohesiveKin::testIntegrateJacobianAssembled(void)
{ // testIntegrateJacobianAssembled
-#if 0
- ALE::Obj<Mesh> mesh;
+ topology::Mesh mesh;
FaultCohesiveKin fault;
- _initialize(&mesh, &fault);
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &fault, &fields);
- // Setup fields
- topology::FieldsManager fields(mesh);
- fields.addReal("residual");
- fields.addReal("solution");
- fields.solutionField("solution");
-
- const ALE::Obj<real_section_type>& residual = fields.getReal("residual");
- CPPUNIT_ASSERT(!residual.isNull());
const int spaceDim = _data->spaceDim;
- residual->setChart(mesh->getSieve()->getChart());
- residual->setFiberDimension(mesh->depthStratum(0), spaceDim);
- mesh->allocate(residual);
- residual->zero();
- fields.copyLayout("residual");
+ const ALE::Obj<RealSection>& solutionSection = fields.get("solution").section();
+ CPPUNIT_ASSERT(!solutionSection.isNull());
- const ALE::Obj<real_section_type>& solution = fields.getReal("solution");
- CPPUNIT_ASSERT(!solution.isNull());
-
- const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices = sieveMesh->depthStratum(0);
CPPUNIT_ASSERT(!vertices.isNull());
- const Mesh::label_sequence::iterator verticesBegin = vertices->begin();
- const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+ const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
int iVertex = 0;
- for (Mesh::label_sequence::iterator v_iter=verticesBegin;
+ for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
- ++v_iter, ++iVertex) {
- solution->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
- } // for
-
- PetscMat jacobian;
- PetscErrorCode err = MeshCreateMatrix(mesh, solution, MATMPIBAIJ, &jacobian);
- CPPUNIT_ASSERT(0 == err);
+ ++v_iter, ++iVertex)
+ solutionSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+ topology::Jacobian jacobian(fields);
+
const double t = 2.134;
- fault.integrateJacobian(&jacobian, t, &fields, mesh);
+ fault.integrateJacobian(&jacobian, t, &fields);
CPPUNIT_ASSERT_EQUAL(false, fault.needNewJacobian());
- err = MatAssemblyBegin(jacobian, MAT_FINAL_ASSEMBLY);
- CPPUNIT_ASSERT(0 == err);
- err = MatAssemblyEnd(jacobian, MAT_FINAL_ASSEMBLY);
- CPPUNIT_ASSERT(0 == err);
+ jacobian.assemble("final_assembly");
//MatView(jacobian, PETSC_VIEWER_STDOUT_WORLD); // DEBUGGING
const double* valsE = _data->valsJacobian;
- const int nrowsE = solution->sizeWithBC();
+ const int nrowsE = solutionSection->sizeWithBC();
const int ncolsE = nrowsE;
+ PetscMat jacobianMat = jacobian.matrix();
+
int nrows = 0;
int ncols = 0;
- MatGetSize(jacobian, &nrows, &ncols);
+ MatGetSize(jacobianMat, &nrows, &ncols);
CPPUNIT_ASSERT_EQUAL(nrowsE, nrows);
CPPUNIT_ASSERT_EQUAL(ncolsE, ncols);
PetscMat jDense;
PetscMat jSparseAIJ;
- MatConvert(jacobian, MATSEQAIJ, MAT_INITIAL_MATRIX, &jSparseAIJ);
+ MatConvert(jacobianMat, MATSEQAIJ, MAT_INITIAL_MATRIX, &jSparseAIJ);
MatConvert(jSparseAIJ, MATSEQDENSE, MAT_INITIAL_MATRIX, &jDense);
double_array vals(nrows*ncols);
@@ -668,74 +612,65 @@
MatDestroy(jDense);
MatDestroy(jSparseAIJ);
CPPUNIT_ASSERT_EQUAL(false, fault.needNewJacobian());
-#endif
} // testIntegrateJacobianAssembled
// ----------------------------------------------------------------------
-// Test updateState().
+// Test updateStateVars().
void
-pylith::faults::TestFaultCohesiveKin::testUpdateState(void)
-{ // testUpdateState
-#if 0
- ALE::Obj<Mesh> mesh;
+pylith::faults::TestFaultCohesiveKin::testUpdateStateVars(void)
+{ // testUpdateStateVars
+ topology::Mesh mesh;
FaultCohesiveKin fault;
- _initialize(&mesh, &fault);
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &fault, &fields);
- spatialdata::geocoords::CSCart cs;
- cs.setSpaceDim((mesh)->getDimension());
- cs.initialize();
-
- // Setup fields
- topology::FieldsManager fields(mesh);
- fields.addReal("residual");
- fields.addReal("solution");
- fields.solutionField("solution");
-
- const ALE::Obj<real_section_type>& residual = fields.getReal("residual");
- CPPUNIT_ASSERT(!residual.isNull());
const int spaceDim = _data->spaceDim;
- residual->setChart(mesh->getSieve()->getChart());
- residual->setFiberDimension(mesh->depthStratum(0), spaceDim);
- mesh->allocate(residual);
- residual->zero();
- fields.copyLayout("residual");
-
- const ALE::Obj<real_section_type>& solution = fields.getReal("solution");
+ const ALE::Obj<RealSection>& solutionSection = fields.get("solution").section();
+ CPPUNIT_ASSERT(!solutionSection.isNull());
{ // setup solution
- CPPUNIT_ASSERT(!solution.isNull());
- solution->zero();
+ solutionSection->zero();
- const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices = sieveMesh->depthStratum(0);
CPPUNIT_ASSERT(!vertices.isNull());
- const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+ const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
int iVertex = 0;
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
- ++v_iter, ++iVertex) {
- solution->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
- } // for
+ ++v_iter, ++iVertex)
+ solutionSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
} // setup solution
+ topology::Field<topology::Mesh>& residual = fields.get("residual");
const double t = 2.134;
const double dt = 0.01;
fault.useSolnIncr(false);
fault.timeStep(dt);
- fault.integrateResidualAssembled(residual, t, &fields, mesh, &cs);
- fault.updateState(t, &fields, mesh);
+ fault.integrateResidualAssembled(residual, t, &fields);
+ fault.updateStateVars(t, &fields);
- const ALE::Obj<Mesh::label_sequence>& vertices =
- fault._faultMesh->depthStratum(0);
- const Mesh::label_sequence::iterator verticesEnd = vertices->end();
- Mesh::renumbering_type& renumbering = fault._faultMesh->getRenumbering();
+ CPPUNIT_ASSERT(0 != fault._faultMesh);
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = fault._faultMesh->sieveMesh();
+ CPPUNIT_ASSERT(!faultSieveMesh.isNull());
+ const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
+ faultSieveMesh->depthStratum(0);
+ const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
+ SieveSubMesh::renumbering_type& renumbering = faultSieveMesh->getRenumbering();
// Compute expected cumulative slip using eqsrcs
- ALE::Obj<real_section_type> cumSlipE =
- new real_section_type(fault._faultMesh->comm(), fault._faultMesh->debug());
- CPPUNIT_ASSERT(!cumSlipE.isNull());
- cumSlipE->setChart(fault._faultMesh->getSieve()->getChart());
- cumSlipE->setFiberDimension(vertices, spaceDim);
- fault._faultMesh->allocate(cumSlipE);
+ topology::Field<topology::SubMesh> cumSlipE(*fault._faultMesh);
+ cumSlipE.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
+ cumSlipE.allocate();
+ const ALE::Obj<RealSection> cumSlipESection = cumSlipE.section();
+ CPPUNIT_ASSERT(!cumSlipESection.isNull());
+ const ALE::Obj<RealSection> cumSlipSection =
+ fault._fields->get("cumulative slip").section();
+ CPPUNIT_ASSERT(!cumSlipSection.isNull());
+
const FaultCohesiveKin::srcs_type::const_iterator srcsEnd = fault._eqSrcs.end();
for (FaultCohesiveKin::srcs_type::iterator s_iter=fault._eqSrcs.begin();
s_iter != srcsEnd;
@@ -743,17 +678,17 @@
EqKinSrc* src = s_iter->second;
assert(0 != src);
if (t >= src->originTime())
- src->slip(cumSlipE, t, fault._faultMesh);
+ src->slip(&cumSlipE, t);
} // for
int iVertex = 0;
const double tolerance = 1.0e-06;
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
++v_iter, ++iVertex) {
- const Mesh::point_type meshVertex = _data->constraintVertices[iVertex];
+ const SieveSubMesh::point_type meshVertex = _data->constraintVertices[iVertex];
bool found = false;
- for(Mesh::renumbering_type::const_iterator r_iter = renumbering.begin();
+ for(SieveSubMesh::renumbering_type::const_iterator r_iter = renumbering.begin();
r_iter != renumbering.end();
++r_iter) {
if (r_iter->second == *v_iter) {
@@ -764,14 +699,12 @@
CPPUNIT_ASSERT(found);
// Check _cumSlip
- int fiberDim = fault._cumSlip->getFiberDimension(*v_iter);
+ int fiberDim = cumSlipSection->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const real_section_type::value_type* slipV =
- fault._cumSlip->restrictPoint(*v_iter);
+ const double* slipV = cumSlipSection->restrictPoint(*v_iter);
CPPUNIT_ASSERT(0 != slipV);
- const real_section_type::value_type* slipE =
- cumSlipE->restrictPoint(*v_iter);
+ const double* slipE = cumSlipESection->restrictPoint(*v_iter);
CPPUNIT_ASSERT(0 != slipE);
for (int iDim=0; iDim < spaceDim; ++iDim) {
@@ -781,187 +714,188 @@
CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE[iDim], slipV[iDim], tolerance);
} // for
} // for
-#endif
-} // testUpdateState
+} // testUpdateStateVars
// ----------------------------------------------------------------------
// Test calcTractionsChange().
void
pylith::faults::TestFaultCohesiveKin::testCalcTractionsChange(void)
{ // testCalcTractionsChange
-#if 0
- ALE::Obj<Mesh> mesh;
+ topology::Mesh mesh;
FaultCohesiveKin fault;
- _initialize(&mesh, &fault);
-
- // Setup fields
- topology::FieldsManager fields(mesh);
- fields.addReal("solution");
- fields.solutionField("solution");
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &fault, &fields);
const int spaceDim = _data->spaceDim;
- const ALE::Obj<real_section_type>& solution = fields.getReal("solution");
- { // setup solution
- CPPUNIT_ASSERT(!solution.isNull());
- solution->setChart(mesh->getSieve()->getChart());
- solution->setFiberDimension(mesh->depthStratum(0), spaceDim);
- mesh->allocate(solution);
- solution->zero();
- fields.copyLayout("solution");
-
- const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const Mesh::label_sequence::iterator verticesEnd = vertices->end();
- int iVertex = 0;
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) {
- solution->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
- } // for
- } // setup solution
+ const ALE::Obj<RealSection>& solutionSection = fields.get("solution").section();
+ CPPUNIT_ASSERT(!solutionSection.isNull());
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ int iVertex = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
+ v_iter != verticesEnd;
+ ++v_iter, ++iVertex) {
+ solutionSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+ } // for
- ALE::Obj<real_section_type> tractions =
- new real_section_type(fault._faultMesh->comm(), fault._faultMesh->debug());
- CPPUNIT_ASSERT(!tractions.isNull());
- Mesh::renumbering_type& renumbering = fault._faultMesh->getRenumbering();
- const ALE::Obj<Mesh::label_sequence>& vertices =
- fault._faultMesh->depthStratum(0);
- const Mesh::label_sequence::iterator verticesEnd = vertices->end();
- tractions->setChart(fault._faultMesh->getSieve()->getChart());
- tractions->setFiberDimension(vertices, spaceDim);
- fault._faultMesh->allocate(tractions);
+ CPPUNIT_ASSERT(0 != fault._faultMesh);
+ topology::Field<topology::SubMesh> tractions(*fault._faultMesh);
+ tractions.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
+ tractions.allocate();
+ tractions.zero();
+ const ALE::Obj<RealSection>& tractionsSection = tractions.section();
+ CPPUNIT_ASSERT(!tractionsSection.isNull());
+
const double t = 0;
- fault.updateState(t, &fields, mesh);
- fault._calcTractionsChange(&tractions, mesh, solution);
+ fault.updateStateVars(t, &fields);
+ fault._calcTractionsChange(&tractions, fields.get("solution"));
- int iVertex = 0;
+ iVertex = 0;
const double tolerance = 1.0e-06;
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = fault._faultMesh->sieveMesh();
+ CPPUNIT_ASSERT(!faultSieveMesh.isNull());
+ SieveSubMesh::renumbering_type& renumbering = faultSieveMesh->getRenumbering();
+ const SieveMesh::renumbering_type::const_iterator rEnd = renumbering.end();
+ for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
++v_iter, ++iVertex) {
- Mesh::point_type meshVertex = -1;
+ SieveMesh::point_type meshVertex = -1;
bool found = false;
- for(Mesh::renumbering_type::const_iterator r_iter = renumbering.begin(); r_iter != renumbering.end(); ++r_iter) {
+ for (SieveMesh::renumbering_type::const_iterator r_iter = renumbering.begin();
+ r_iter != rEnd;
+ ++r_iter) {
if (r_iter->second == *v_iter) {
meshVertex = r_iter->first;
- found = true;
+ found = true;
break;
- }
- }
+ } // if
+ } // for
CPPUNIT_ASSERT(found);
- int fiberDim = tractions->getFiberDimension(*v_iter);
+ int fiberDim = tractionsSection->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const real_section_type::value_type* vertexTractions =
- tractions->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vertexTractions);
+ const double* tractionsVertex = tractionsSection->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != tractionsVertex);
- fiberDim = solution->getFiberDimension(meshVertex);
+ fiberDim = solutionSection->getFiberDimension(meshVertex);
CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const real_section_type::value_type* vertexSolution =
- solution->restrictPoint(meshVertex);
- CPPUNIT_ASSERT(0 != vertexSolution);
+ const double* solutionVertex = solutionSection->restrictPoint(meshVertex);
+ CPPUNIT_ASSERT(0 != solutionVertex);
const double scale = _data->pseudoStiffness / _data->area[iVertex];
for (int iDim=0; iDim < spaceDim; ++iDim) {
- const double tractionE = vertexSolution[iDim] * scale;
+ const double tractionE = solutionVertex[iDim] * scale;
if (tractionE > 1.0)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vertexTractions[iDim]/tractionE,
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, tractionsVertex[iDim]/tractionE,
tolerance);
else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(tractionE, vertexTractions[iDim],
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(tractionE, tractionsVertex[iDim],
tolerance);
} // for
} // for
-#endif
} // testCalcTractionsChange
// ----------------------------------------------------------------------
// Initialize FaultCohesiveKin interface condition.
void
pylith::faults::TestFaultCohesiveKin::_initialize(
- topology::Mesh* mesh,
- FaultCohesiveKin* const fault) const
+ topology::Mesh* const mesh,
+ FaultCohesiveKin* const fault,
+ topology::SolutionFields* const fields) const
{ // _initialize
CPPUNIT_ASSERT(0 != mesh);
CPPUNIT_ASSERT(0 != fault);
+ CPPUNIT_ASSERT(0 != fields);
CPPUNIT_ASSERT(0 != _data);
CPPUNIT_ASSERT(0 != _quadrature);
- try {
- meshio::MeshIOAscii iohandler;
- iohandler.filename(_data->meshFilename);
- iohandler.read(mesh);
-
- //(*mesh)->setDebug(true); // DEBUGGING
-
- spatialdata::geocoords::CSCart cs;
- cs.setSpaceDim(mesh->dimension());
- cs.initialize();
-
- _quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
- _data->basisDeriv,
- _data->numQuadPts, _data->numBasis, _data->cellDim,
- _data->quadPts, _data->numQuadPts, _data->cellDim,
- _data->quadWts, _data->numQuadPts,
- _data->spaceDim);
-
- // Setup earthquake source
- spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
- spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
- ioFinalSlip.filename(_data->finalSlipFilename);
- dbFinalSlip.ioHandler(&ioFinalSlip);
+ meshio::MeshIOAscii iohandler;
+ iohandler.filename(_data->meshFilename);
+ iohandler.read(mesh);
- spatialdata::spatialdb::SimpleDB dbSlipTime("slip time");
- spatialdata::spatialdb::SimpleIOAscii ioSlipTime;
- ioSlipTime.filename(_data->slipTimeFilename);
- dbSlipTime.ioHandler(&ioSlipTime);
+ //(*mesh)->setDebug(true); // DEBUGGING
- spatialdata::spatialdb::SimpleDB dbRiseTime("rise time");
- spatialdata::spatialdb::SimpleIOAscii ioRiseTime;
- ioRiseTime.filename(_data->riseTimeFilename);
- dbRiseTime.ioHandler(&ioRiseTime);
-
- const int nsrcs = _eqsrcs.size();
- CPPUNIT_ASSERT(nsrcs == _slipfns.size());
- EqKinSrc** sources = new EqKinSrc*[nsrcs];
- char** names = new char*[nsrcs];
- for (int i=0; i < nsrcs; ++i) {
- _slipfns[i]->dbFinalSlip(&dbFinalSlip);
- _slipfns[i]->dbSlipTime(&dbSlipTime);
- _slipfns[i]->dbRiseTime(&dbRiseTime);
-
- _eqsrcs[i]->slipfn(_slipfns[i]);
- sources[i] = _eqsrcs[i];
- names[i] = new char[2];
- names[i][0] = 'a' + i;
- names[i][1] = '\0';
- } // for
+ spatialdata::geocoords::CSCart cs;
+ cs.setSpaceDim(mesh->dimension());
+ cs.initialize();
+ mesh->coordsys(&cs);
- fault->id(_data->id);
- fault->label(_data->label);
- fault->quadrature(_quadrature);
+ _quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
+ _data->basisDeriv,
+ _data->numQuadPts, _data->numBasis, _data->cellDim,
+ _data->quadPts, _data->numQuadPts, _data->cellDim,
+ _data->quadWts, _data->numQuadPts,
+ _data->spaceDim);
+
+ // Setup earthquake source
+ spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
+ spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
+ ioFinalSlip.filename(_data->finalSlipFilename);
+ dbFinalSlip.ioHandler(&ioFinalSlip);
+
+ spatialdata::spatialdb::SimpleDB dbSlipTime("slip time");
+ spatialdata::spatialdb::SimpleIOAscii ioSlipTime;
+ ioSlipTime.filename(_data->slipTimeFilename);
+ dbSlipTime.ioHandler(&ioSlipTime);
+
+ spatialdata::spatialdb::SimpleDB dbRiseTime("rise time");
+ spatialdata::spatialdb::SimpleIOAscii ioRiseTime;
+ ioRiseTime.filename(_data->riseTimeFilename);
+ dbRiseTime.ioHandler(&ioRiseTime);
+
+ const int nsrcs = _eqsrcs.size();
+ CPPUNIT_ASSERT(nsrcs == _slipfns.size());
+ EqKinSrc** sources = new EqKinSrc*[nsrcs];
+ char** names = new char*[nsrcs];
+ for (int i=0; i < nsrcs; ++i) {
+ _slipfns[i]->dbFinalSlip(&dbFinalSlip);
+ _slipfns[i]->dbSlipTime(&dbSlipTime);
+ _slipfns[i]->dbRiseTime(&dbRiseTime);
- fault->eqsrcs(const_cast<const char**>(names), sources, nsrcs);
- fault->adjustTopology(mesh, _flipFault);
-
- const double upDir[] = { 0.0, 0.0, 1.0 };
- const double normalDir[] = { 1.0, 0.0, 0.0 };
-
- spatialdata::spatialdb::SimpleDB dbMatProp("material properties");
- spatialdata::spatialdb::SimpleIOAscii ioMatProp;
- ioMatProp.filename(_data->matPropsFilename);
- dbMatProp.ioHandler(&ioMatProp);
-
- fault->initialize(*mesh, upDir, normalDir, &dbMatProp);
-
- delete[] sources; sources = 0;
- for (int i=0; i < nsrcs; ++i)
- delete[] names[i];
- delete[] names; names = 0;
- } catch (const ALE::Exception& err) {
- throw std::runtime_error(err.msg());
- } // catch
+ _eqsrcs[i]->slipfn(_slipfns[i]);
+ sources[i] = _eqsrcs[i];
+ names[i] = new char[2];
+ names[i][0] = 'a' + i;
+ names[i][1] = '\0';
+ } // for
+
+ fault->id(_data->id);
+ fault->label(_data->label);
+ fault->quadrature(_quadrature);
+
+ fault->eqsrcs(const_cast<const char**>(names), sources, nsrcs);
+ fault->adjustTopology(mesh, _flipFault);
+
+ const double upDir[] = { 0.0, 0.0, 1.0 };
+ const double normalDir[] = { 1.0, 0.0, 0.0 };
+
+ spatialdata::spatialdb::SimpleDB dbMatProp("material properties");
+ spatialdata::spatialdb::SimpleIOAscii ioMatProp;
+ ioMatProp.filename(_data->matPropsFilename);
+ dbMatProp.ioHandler(&ioMatProp);
+
+ fault->initialize(*mesh, upDir, normalDir, &dbMatProp);
+
+ delete[] sources; sources = 0;
+ for (int i=0; i < nsrcs; ++i)
+ delete[] names[i];
+ delete[] names; names = 0;
+
+ // Setup fields
+ fields->add("residual", "residual");
+ fields->add("solution", "displacement");
+ fields->solutionName("solution");
+
+ const int spaceDim = _data->spaceDim;
+ topology::Field<topology::Mesh>& residual = fields->get("residual");
+ residual.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
+ residual.allocate();
+ fields->copyLayout("residual");
} // _initialize
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKin.hh 2009-04-27 21:01:09 UTC (rev 14804)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKin.hh 2009-04-27 21:28:50 UTC (rev 14805)
@@ -101,8 +101,8 @@
/// Test integrateJacobianAssembled().
void testIntegrateJacobianAssembled(void);
- /// Test updateState().
- void testUpdateState(void);
+ /// Test updateStateVars().
+ void testUpdateStateVars(void);
/// Test _calcTractionsChange().
void testCalcTractionsChange(void);
@@ -114,9 +114,11 @@
*
* @param mesh PETSc mesh to initialize
* @param fault Cohesive fault interface condition to initialize.
+ * @param fields Solution fields.
*/
- void _initialize(topology::Mesh* mesh,
- FaultCohesiveKin* const fault) const;
+ void _initialize(topology::Mesh* const mesh,
+ FaultCohesiveKin* const fault,
+ topology::SolutionFields* const fields) const;
/** Determine if vertex is a Lagrange multiplier constraint vertex.
*
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinHex8.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinHex8.hh 2009-04-27 21:01:09 UTC (rev 14804)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinHex8.hh 2009-04-27 21:28:50 UTC (rev 14805)
@@ -40,7 +40,7 @@
CPPUNIT_TEST( testInitialize );
CPPUNIT_TEST( testIntegrateResidual );
CPPUNIT_TEST( testIntegrateJacobian );
- CPPUNIT_TEST( testUpdateState );
+ CPPUNIT_TEST( testUpdateStateVars );
CPPUNIT_TEST( testCalcTractionsChange );
CPPUNIT_TEST_SUITE_END();
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinLine2.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinLine2.cc 2009-04-27 21:01:09 UTC (rev 14804)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinLine2.cc 2009-04-27 21:28:50 UTC (rev 14805)
@@ -16,7 +16,8 @@
#include "data/CohesiveKinDataLine2.hh" // USES CohesiveKinDataLine2
-#include "pylith/feassemble/Quadrature0D.hh" // USES Quadrature0D
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/feassemble/Quadrature.hh" // USES Quadrature<SubMesh>
#include "pylith/feassemble/GeometryPoint1D.hh" // USES GeometryPoint1D
// ----------------------------------------------------------------------
@@ -29,8 +30,8 @@
{ // setUp
TestFaultCohesiveKin::setUp();
_data = new CohesiveKinDataLine2();
- _quadrature = new feassemble::Quadrature0D();
- CPPUNIT_ASSERT(0 != _quadrature);
+
+ assert(0 != _quadrature);
feassemble::GeometryPoint1D geometry;
_quadrature->refGeometry(&geometry);
} // setUp
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinLine2.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinLine2.hh 2009-04-27 21:01:09 UTC (rev 14804)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinLine2.hh 2009-04-27 21:28:50 UTC (rev 14805)
@@ -40,7 +40,7 @@
CPPUNIT_TEST( testInitialize );
CPPUNIT_TEST( testIntegrateResidual );
CPPUNIT_TEST( testIntegrateJacobian );
- CPPUNIT_TEST( testUpdateState );
+ CPPUNIT_TEST( testUpdateStateVars );
CPPUNIT_TEST( testCalcTractionsChange );
CPPUNIT_TEST_SUITE_END();
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinQuad4.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinQuad4.hh 2009-04-27 21:01:09 UTC (rev 14804)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinQuad4.hh 2009-04-27 21:28:50 UTC (rev 14805)
@@ -40,7 +40,7 @@
CPPUNIT_TEST( testInitialize );
CPPUNIT_TEST( testIntegrateResidual );
CPPUNIT_TEST( testIntegrateJacobian );
- CPPUNIT_TEST( testUpdateState );
+ CPPUNIT_TEST( testUpdateStateVars );
CPPUNIT_TEST( testCalcTractionsChange );
CPPUNIT_TEST_SUITE_END();
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinQuad4e.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinQuad4e.hh 2009-04-27 21:01:09 UTC (rev 14804)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinQuad4e.hh 2009-04-27 21:28:50 UTC (rev 14805)
@@ -40,7 +40,7 @@
CPPUNIT_TEST( testInitialize );
CPPUNIT_TEST( testIntegrateResidual );
CPPUNIT_TEST( testIntegrateJacobian );
- CPPUNIT_TEST( testUpdateState );
+ CPPUNIT_TEST( testUpdateStateVars );
CPPUNIT_TEST( testCalcTractionsChange );
CPPUNIT_TEST_SUITE_END();
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinTri3.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinTri3.cc 2009-04-27 21:01:09 UTC (rev 14804)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinTri3.cc 2009-04-27 21:28:50 UTC (rev 14805)
@@ -16,7 +16,8 @@
#include "data/CohesiveKinDataTri3.hh" // USES CohesiveKinDataTri3
-#include "pylith/feassemble/Quadrature1Din2D.hh" // USES Quadrature1Din2D
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/feassemble/Quadrature.hh" // USES Quadrature<SubMesh>
#include "pylith/feassemble/GeometryLine2D.hh" // USES GeometryLine2D
// ----------------------------------------------------------------------
@@ -29,7 +30,7 @@
{ // setUp
TestFaultCohesiveKin::setUp();
_data = new CohesiveKinDataTri3();
- _quadrature = new feassemble::Quadrature1Din2D();
+
CPPUNIT_ASSERT(0 != _quadrature);
feassemble::GeometryLine2D geometry;
_quadrature->refGeometry(&geometry);
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinTri3.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinTri3.hh 2009-04-27 21:01:09 UTC (rev 14804)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesiveKinTri3.hh 2009-04-27 21:28:50 UTC (rev 14805)
@@ -40,7 +40,7 @@
CPPUNIT_TEST( testInitialize );
CPPUNIT_TEST( testIntegrateResidual );
CPPUNIT_TEST( testIntegrateJacobian );
- CPPUNIT_TEST( testUpdateState );
+ CPPUNIT_TEST( testUpdateStateVars );
CPPUNIT_TEST( testCalcTractionsChange );
CPPUNIT_TEST_SUITE_END();
More information about the CIG-COMMITS
mailing list