[cig-commits] r15091 - in short/3D/PyLith/trunk: . libsrc/bc libsrc/faults libsrc/feassemble libsrc/topology modulesrc/feassemble pylith/problems unittests/pytests/topology
brad at geodynamics.org
brad at geodynamics.org
Sat May 30 20:27:41 PDT 2009
Author: brad
Date: 2009-05-30 20:27:40 -0700 (Sat, 30 May 2009)
New Revision: 15091
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/bc/DirichletBC.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc
short/3D/PyLith/trunk/libsrc/topology/Field.cc
short/3D/PyLith/trunk/libsrc/topology/Field.hh
short/3D/PyLith/trunk/modulesrc/feassemble/Integrator.i
short/3D/PyLith/trunk/pylith/problems/Formulation.py
short/3D/PyLith/trunk/unittests/pytests/topology/TestMeshField.py
Log:
Worked on setting up field split.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2009-05-30 21:42:31 UTC (rev 15090)
+++ short/3D/PyLith/trunk/TODO 2009-05-31 03:27:40 UTC (rev 15091)
@@ -7,11 +7,10 @@
field split
Brad
- symmetric matrix
+ field split (setup fields)
+ symmetric matrix (use in examples, testing)
full-scale testing
cleanup
- Questions for Matt:
- How to use symmetric matrix (updateOperator)?
test uniform refinement
PointForce
Modified: short/3D/PyLith/trunk/libsrc/bc/DirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/DirichletBC.cc 2009-05-30 21:42:31 UTC (rev 15090)
+++ short/3D/PyLith/trunk/libsrc/bc/DirichletBC.cc 2009-05-31 03:27:40 UTC (rev 15091)
@@ -85,6 +85,7 @@
const ALE::Obj<RealSection>& section = field.section();
assert(!section.isNull());
+ // Set constraints in field
const int numPoints = _points.size();
_offsetLocal.resize(numPoints);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
@@ -104,6 +105,20 @@
_offsetLocal[iPoint] = curNumConstraints;
section->addConstraintDimension(_points[iPoint], numFixedDOF);
} // for
+
+ // Set constraints in fibration 0 (split field) if it exists.
+ if (section->getNumSpaces() > 0) {
+ const int fibration = 0;
+ const ALE::Obj<RealSection>& splitSection =
+ section->getFibration(fibration);
+ for (int iPoint=0; iPoint < numPoints; ++iPoint) {
+ const int fiberDim = splitSection->getFiberDimension(_points[iPoint]);
+ const int curNumConstraints =
+ splitSection->getConstraintDimension(_points[iPoint]);
+ assert(curNumConstraints + numFixedDOF <= fiberDim);
+ splitSection->addConstraintDimension(_points[iPoint], numFixedDOF);
+ } // for
+ } // if
} // setConstraintSizes
// ----------------------------------------------------------------------
@@ -118,6 +133,10 @@
const ALE::Obj<RealSection>& section = field.section();
assert(!section.isNull());
+ const int fibration = (section->getNumSpaces() > 0) ? 0 : -1;
+ const ALE::Obj<RealSection>& splitSection = (section->getNumSpaces() > 0) ?
+ section->getFibration(fibration) : 0;
+
const int numPoints = _points.size();
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const SieveMesh::point_type point = _points[iPoint];
@@ -164,6 +183,9 @@
// Update list of constrained DOF
section->setConstraintDof(point, &allFixedDOF[0]);
+
+ if (fibration >= 0)
+ splitSection->setConstraintDof(point, &allFixedDOF[0]);
} // for
} // setConstraints
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2009-05-30 21:42:31 UTC (rev 15090)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2009-05-31 03:27:40 UTC (rev 15091)
@@ -149,6 +149,56 @@
} // initialize
// ----------------------------------------------------------------------
+void
+pylith::faults::FaultCohesiveKin::splitFields(topology::Field<topology::Mesh>* field)
+{ // splitFields
+ assert(0 != field);
+
+ const ALE::Obj<RealSection>& section = field->section();
+ assert(!section.isNull());
+ if (0 == section->getNumSpaces())
+ return; // Return if there are no fibrations
+
+ const int fibrationDisp = 0;
+ const int fibrationLagrange = 1;
+ const ALE::Obj<RealSection>& splitSection =
+ section->getFibration(fibrationLagrange);
+ assert(!splitSection.isNull());
+
+ // Get domain Sieve mesh
+ const ALE::Obj<SieveMesh>& sieveMesh = field->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+
+ // Get fault Sieve mesh
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
+
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ assert(!vertices.isNull());
+ 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) != renumberingEnd) {
+ const int vertexFault = renumbering[*v_iter];
+ const int vertexMesh = *v_iter;
+ const int fiberDim = section->getFiberDimension(vertexMesh);
+ assert(fiberDim > 0);
+ // Reset displacement fibration fiber dimension to zero.
+ section->setFiberDimension(vertexMesh, 0, fibrationDisp);
+ // Set Langrange fibration fiber dimension.
+ splitSection->setFiberDimension(vertexMesh, fiberDim, fibrationLagrange);
+ } // if
+} // splitFields
+
+// ----------------------------------------------------------------------
// Integrate contribution of cohesive cells to residual term that
// require assembly across processors.
void
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2009-05-30 21:42:31 UTC (rev 15090)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2009-05-31 03:27:40 UTC (rev 15091)
@@ -90,6 +90,12 @@
const double upDir[3],
const double normalDir[3]);
+ /** Split solution fields for separate preconditioning.
+ *
+ * @param field Solution field.
+ */
+ void splitFields(topology::Field<topology::Mesh>* field);
+
/** Integrate contributions to residual term (r) for operator that
* require assembly across processors.
*
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2009-05-30 21:42:31 UTC (rev 15090)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2009-05-31 03:27:40 UTC (rev 15091)
@@ -116,6 +116,13 @@
virtual
void initialize(const topology::Mesh& mesh);
+ /** Split solution fields for separate preconditioning.
+ *
+ * @param field Solution field.
+ */
+ virtual
+ void splitFields(topology::Field<topology::Mesh>* field);
+
/** Integrate contributions to residual term (r) for operator.
*
* @param residual Field containing values for residual
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc 2009-05-30 21:42:31 UTC (rev 15090)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc 2009-05-31 03:27:40 UTC (rev 15091)
@@ -44,6 +44,13 @@
void
pylith::feassemble::Integrator<quadrature_type>::initialize(const topology::Mesh& mesh) {
} // initialize
+
+// Split solution fields for separate preconditioning.
+template<typename quadrature_type>
+inline
+void
+pylith::feassemble::Integrator<quadrature_type>::splitFields(topology::Field<topology::Mesh>* field) {
+} // splitFields
// Integrate contributions to residual term (r) for operator.
template<typename quadrature_type>
Modified: short/3D/PyLith/trunk/libsrc/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Field.cc 2009-05-30 21:42:31 UTC (rev 15090)
+++ short/3D/PyLith/trunk/libsrc/topology/Field.cc 2009-05-31 03:27:40 UTC (rev 15091)
@@ -219,6 +219,7 @@
_section->setAtlas(srcSection->getAtlas());
_section->allocateStorage();
_section->setBC(srcSection->getBC());
+ _section->copyFibration(srcSection);
if (0 != src._scatter) {
_scatter = src._scatter;
@@ -640,5 +641,28 @@
err = VecDestroy(localVec); CHECK_PETSC_ERROR(err);
} // scatterVectorToSection
+// ----------------------------------------------------------------------
+// Setup split field with all entries set to a default space of 0.
+template<typename mesh_type>
+void
+pylith::topology::Field<mesh_type>::splitDefault(void)
+{ // splitDefault
+ assert(!_section.isNull());
+ _section->addSpace(); // displacements
+ _section->addSpace(); // Lagrange multipliers
+ // Assume fiber dimension is uniform
+ const chart_type& chart = _section->getChart();
+ const int fiberDim = (chart.size() > 0) ?
+ _section->getFiberDimension(*chart.begin()) : 0;
+
+ const int fibration = 0;
+ const typename chart_type::const_iterator chartBegin = chart.begin();
+ const typename chart_type::const_iterator chartEnd = chart.end();
+ for (typename chart_type::const_iterator c_iter = chart.begin();
+ c_iter != chartEnd;
+ ++c_iter)
+ _section->setFiberDimension(*c_iter, fiberDim, fibration);
+} // splitDefault
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/topology/Field.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Field.hh 2009-05-30 21:42:31 UTC (rev 15090)
+++ short/3D/PyLith/trunk/libsrc/topology/Field.hh 2009-05-31 03:27:40 UTC (rev 15091)
@@ -266,6 +266,9 @@
*/
void scatterVectorToSection(const PetscVec vector) const;
+ /// Setup split field with all entries set to a default space of 0.
+ void splitDefault(void);
+
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/Integrator.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/Integrator.i 2009-05-30 21:42:31 UTC (rev 15090)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/Integrator.i 2009-05-31 03:27:40 UTC (rev 15091)
@@ -97,6 +97,13 @@
virtual
void initialize(const pylith::topology::Mesh& mesh);
+ /** Split solution fields for separate preconditioning.
+ *
+ * @param field Solution field.
+ */
+ virtual
+ void splitFields(pylith::topology::Field<pylith::topology::Mesh>* field);
+
/** Integrate contributions to residual term (r) for operator.
*
* @param residual Field containing values for residual
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2009-05-30 21:42:31 UTC (rev 15090)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2009-05-31 03:27:40 UTC (rev 15091)
@@ -57,7 +57,10 @@
##
## \b Properties
## @li \b matrix_type Type of PETSc sparse matrix.
- ## @li \b view_jacobian Flag to output Jacobian matrix when it is reformed.
+ ## @li \b split_fields Split solution fields into displacements
+ ## and Lagrange multipliers for separate preconditioning.
+ ## @li \b view_jacobian Flag to output Jacobian matrix when it is
+ ## reformed.
##
## \b Facilities
## @li \b time_step Time step size manager.
@@ -70,6 +73,10 @@
matrixType = pyre.inventory.str("matrix_type", default="unknown")
matrixType.meta['tip'] = "Type of PETSc sparse matrix."
+ splitFields = pyre.inventory.bool("split_fields", default=False)
+ splitFields.meta['tip'] = "Split solution fields into displacements "\
+ "and Lagrange multipliers for separate preconditioning."
+
viewJacobian = pyre.inventory.bool("view_jacobian", default=False)
viewJacobian.meta['tip'] = "Write Jacobian matrix to binary file."
@@ -219,12 +226,16 @@
solution.vectorFieldType(solution.VECTOR)
solution.scale(lengthScale.value)
solution.newSection(solution.VERTICES_FIELD, dimension)
+ if self.splitFields:
+ solution.splitDefault()
for constraint in self.constraints:
constraint.setConstraintSizes(solution)
+ if self.splitFields:
+ for integrator in self.integratorsMesh + self.integratorsSubMesh:
+ integrator.splitFields(solution)
solution.allocate()
for constraint in self.constraints:
constraint.setConstraints(solution)
- solution = self.fields.get("dispIncr(t->t+dt)")
solution.createVector()
solution.createScatter()
@@ -350,6 +361,7 @@
"""
PetscComponent._configure(self)
self.matrixType = self.inventory.matrixType
+ self.splitFields = self.inventory.splitFields
self.timeStep = self.inventory.timeStep
self.solver = self.inventory.solver
self.output = self.inventory.output
Modified: short/3D/PyLith/trunk/unittests/pytests/topology/TestMeshField.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/topology/TestMeshField.py 2009-05-30 21:42:31 UTC (rev 15090)
+++ short/3D/PyLith/trunk/unittests/pytests/topology/TestMeshField.py 2009-05-31 03:27:40 UTC (rev 15091)
@@ -137,12 +137,12 @@
return
- def test_newSectionField(self):
+ def test_cloneSectionField(self):
"""
Test newSection(field).
"""
fieldB = MeshField(self.mesh)
- fieldB.newSection(self.field)
+ fieldB.cloneSection(self.field)
# No test of result
return
More information about the CIG-COMMITS
mailing list