[cig-commits] r7248 - in short/3D/PyLith/trunk: . libsrc/faults
pylith/problems tests/1d/line2 tests/1d/line3 tests/2d/quad4
tests/2d/tri3 tests/3d/tet4
brad at geodynamics.org
brad at geodynamics.org
Thu Jun 14 15:52:19 PDT 2007
Author: brad
Date: 2007-06-14 15:52:19 -0700 (Thu, 14 Jun 2007)
New Revision: 7248
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
short/3D/PyLith/trunk/pylith/problems/Formulation.py
short/3D/PyLith/trunk/pylith/problems/Implicit.py
short/3D/PyLith/trunk/tests/1d/line2/Makefile.am
short/3D/PyLith/trunk/tests/1d/line3/Makefile.am
short/3D/PyLith/trunk/tests/2d/quad4/dislocation.cfg
short/3D/PyLith/trunk/tests/2d/tri3/dislocation.cfg
short/3D/PyLith/trunk/tests/3d/tet4/dislocation.cfg
Log:
Fixed bug that prevented faults with multiple cohesive cells. Still need to create unit tests to verify implementation works. Faults seem to be working for 1-D, 2-D, and 3-D meshes.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-06-14 21:40:29 UTC (rev 7247)
+++ short/3D/PyLith/trunk/TODO 2007-06-14 22:52:19 UTC (rev 7248)
@@ -2,19 +2,12 @@
MAIN PRIORITIES (Brad)
======================================================================
-FaultCohesiveKin::integrateJacobian()
- Need to prevent multiple contributions of orientation information for
- each vertex, because we want to use ADD_VALUES in updateOperator().
+Unit tests with multiple cohesive cells.
+ tri3
+ quad4
+ tet4
+ hex8
- Make slip directions consistent (left-lateral, fault-opening) across
- 0-D, 1-D, and 2-D. Need to adjust slip for 1-D (orient2D) to be
- left-lateral instead of right lateral (p is lateral, q is out of
- plane, r is normal in order to be consistent with 2-D faults).
-
-Dirichlet
- Add check for overly constrained points.
- (number of previously constrained DOF + number added <= fiberDim)
-
5. Additional unit tests
b. ElasticityExplicit and ElasticityImplicit
i. multiple materials
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-06-14 21:40:29 UTC (rev 7247)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-06-14 22:52:19 UTC (rev 7248)
@@ -226,6 +226,50 @@
} // for
_eqsrc->initialize(mesh, *_faultMesh, _constraintVert, cs);
+
+ // Establish pairing between constraint vertices and first cell they
+ // appear in to prevent overlap in integrating Jacobian
+ const int noCell = -1;
+ _constraintCell = new int_section_type(mesh->comm(), mesh->debug());
+ assert(!_constraintCell.isNull());
+ for (std::set<Mesh::point_type>::const_iterator v_iter=vertCohesiveBegin;
+ v_iter != vertCohesiveEnd;
+ ++v_iter)
+ _constraintCell->setFiberDimension(*v_iter, 1);
+ mesh->allocate(_constraintCell);
+ // Set values to noCell
+ for (std::set<Mesh::point_type>::const_iterator v_iter=vertCohesiveBegin;
+ v_iter != vertCohesiveEnd;
+ ++v_iter)
+ _constraintCell->updatePoint(*v_iter, &noCell);
+
+ for (Mesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
+ c_iter != cellsCohesiveEnd;
+ ++c_iter) {
+ const ALE::Obj<sieve_type::traits::coneSequence>& cone =
+ sieve->cone(*c_iter);
+ assert(!cone.isNull());
+ const sieve_type::traits::coneSequence::iterator vBegin = cone->begin();
+ const sieve_type::traits::coneSequence::iterator vEnd = cone->end();
+ const int coneSize = cone->size();
+ assert(coneSize % 3 == 0);
+ sieve_type::traits::coneSequence::iterator v_iter = vBegin;
+ // Skip over non-constraint vertices
+ for (int i=0, numSkip=2*coneSize/3; i < numSkip; ++i)
+ ++v_iter;
+ // If haven't set cell-constraint pair, then set it for current
+ // cell, otherwise move on.
+ for(int i=0, numConstraintVert=coneSize/3;
+ i < numConstraintVert;
+ ++i, ++v_iter) {
+ const int_section_type::value_type* curCell =
+ _constraintCell->restrictPoint(*v_iter);
+ if (noCell == *curCell) {
+ int point = *c_iter;
+ _constraintCell->updatePoint(*v_iter, &point);
+ } // if
+ } // for
+ } // for
} // initialize
// ----------------------------------------------------------------------
@@ -298,16 +342,26 @@
const int numConstraintVert = _quadrature->numBasis();
const int numCorners = 3*numConstraintVert; // cohesive cell
double_array cellMatrix(numCorners*spaceDim * numCorners*spaceDim);
+ double_array cellOrientation(numConstraintVert*orientationSize);
+ int_array cellConstraintCell(numConstraintVert);
for (Mesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
c_iter != cellsCohesiveEnd;
++c_iter) {
cellMatrix = 0.0;
// Get orientations at cells vertices (only valid at constraint vertices)
- const real_section_type::value_type* cellOrientation =
- mesh->restrict(_orientation, *c_iter);
+ mesh->restrict(_orientation, *c_iter, &cellOrientation[0],
+ cellOrientation.size());
+ // Get constraint/cell pairings (only valid at constraint vertices)
+ mesh->restrict(_constraintCell, *c_iter, &cellConstraintCell[0],
+ cellConstraintCell.size());
+
for (int iConstraint=0; iConstraint < numConstraintVert; ++iConstraint) {
+ // Skip setting values if they are set by another cell
+ if (cellConstraintCell[iConstraint] != *c_iter)
+ continue;
+
// Blocks in cell matrix associated with normal cohesive
// vertices i and j and constraint vertex k
const int indexI = iConstraint;
@@ -347,13 +401,9 @@
// Assemble cell contribution into PETSc Matrix
const ALE::Obj<Mesh::order_type>& globalOrder =
mesh->getFactory()->getGlobalOrder(mesh, "default", disp);
- // Update values (do not add)
-
- // Most integrators will call the PETSc updateOperator() routine
- // with ADD_VALUES, but with Lagrange multipler constraints, we
- // call updateOperator() with INSERT_VALUES.
-
- // :BUG: NEED TO USE INSERT_VALUES HERE
+ // Note: We are not really adding values because we prevent
+ // overlap across cells. We use ADD_VALUES for compatibility with
+ // the other integrators.
err = updateOperator(*mat, mesh, disp, globalOrder,
*c_iter, &cellMatrix[0], ADD_VALUES);
if (err)
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2007-06-14 21:40:29 UTC (rev 7247)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2007-06-14 22:52:19 UTC (rev 7248)
@@ -134,10 +134,10 @@
/// Fault vertices associated with constraints
std::set<Mesh::point_type> _constraintVert;
- /// Label of cell used to compute Jacobian for each vertex (must
- /// prevent overlap so only allow 1 cell to contribute for each
- /// vertex).
- ALE::Obj<int_section_type> _vertexCells;
+ /// Label of cell used to compute Jacobian for each constraint vertex (must
+ /// prevent overlap so that only allow 1 cell will contribute for
+ /// each vertex).
+ ALE::Obj<int_section_type> _constraintCell;
}; // class FaultCohesiveKin
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2007-06-14 21:40:29 UTC (rev 7247)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2007-06-14 22:52:19 UTC (rev 7248)
@@ -157,6 +157,11 @@
self.fields.allocate(solnName)
for constraint in self.constraints:
constraint.setConstraints(self.fields.getReal(solnName))
+
+ # BEGIN TEMPORARY
+ import pylith.topology.topology as bindings
+ bindings.sectionView(self.fields.getReal(solnName), "solution")
+ # END TEMPORARY
return
Modified: short/3D/PyLith/trunk/pylith/problems/Implicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Implicit.py 2007-06-14 21:40:29 UTC (rev 7247)
+++ short/3D/PyLith/trunk/pylith/problems/Implicit.py 2007-06-14 22:52:19 UTC (rev 7248)
@@ -95,7 +95,6 @@
self.fields.copyLayout(self.solnField['name'])
dispTBctpdt = self.fields.getReal("dispTBctpdt")
- import pylith.topology.topology as bindings
self.jacobian = mesh.createMatrix(self.fields.getReal("dispTBctpdt"))
Modified: short/3D/PyLith/trunk/tests/1d/line2/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/tests/1d/line2/Makefile.am 2007-06-14 21:40:29 UTC (rev 7247)
+++ short/3D/PyLith/trunk/tests/1d/line2/Makefile.am 2007-06-14 22:52:19 UTC (rev 7248)
@@ -30,8 +30,10 @@
matprops.spatialdb
noinst_TMP = \
- axialextension.vtk \
- dislocation.vtk
+ axialextension_t0.vtk \
+ axialextension_t1.vtk \
+ dislocation_t0.vtk \
+ dislocation_t1.vtk
TESTS_ENVIRONMENT = $(PYTHON)
Modified: short/3D/PyLith/trunk/tests/1d/line3/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/tests/1d/line3/Makefile.am 2007-06-14 21:40:29 UTC (rev 7247)
+++ short/3D/PyLith/trunk/tests/1d/line3/Makefile.am 2007-06-14 22:52:19 UTC (rev 7248)
@@ -30,8 +30,10 @@
matprops.spatialdb
noinst_TMP = \
- axialextension.vtk \
- dislocation.vtk
+ axialextension_t0.vtk \
+ axialextension_t1.vtk \
+ dislocation_t0.vtk \
+ dislocation_t1.vtk
TESTS_ENVIRONMENT = $(PYTHON)
Modified: short/3D/PyLith/trunk/tests/2d/quad4/dislocation.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/dislocation.cfg 2007-06-14 21:40:29 UTC (rev 7247)
+++ short/3D/PyLith/trunk/tests/2d/quad4/dislocation.cfg 2007-06-14 22:52:19 UTC (rev 7248)
@@ -56,7 +56,7 @@
# boundary conditions
# ----------------------------------------------------------------------
[pylithapp.timedependent.bc.neg]
-fixed_dof = [0]
+fixed_dof = [0, 1]
id = 10
label = x_neg
db.label = Dirichlet BC -x edge
@@ -64,7 +64,7 @@
db.query_type = linear
[pylithapp.timedependent.bc.pos]
-fixed_dof = [0]
+fixed_dof = [0, 1]
id = 11
label = x_pos
db.label = Dirichlet BC +x edge
Modified: short/3D/PyLith/trunk/tests/2d/tri3/dislocation.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/dislocation.cfg 2007-06-14 21:40:29 UTC (rev 7247)
+++ short/3D/PyLith/trunk/tests/2d/tri3/dislocation.cfg 2007-06-14 22:52:19 UTC (rev 7248)
@@ -55,7 +55,7 @@
# boundary conditions
# ----------------------------------------------------------------------
[pylithapp.timedependent.bc.bc]
-fixed_dof = [0]
+fixed_dof = [0, 1]
id = 10
label = end points
db.label = Dirichlet BC
Modified: short/3D/PyLith/trunk/tests/3d/tet4/dislocation.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/3d/tet4/dislocation.cfg 2007-06-14 21:40:29 UTC (rev 7247)
+++ short/3D/PyLith/trunk/tests/3d/tet4/dislocation.cfg 2007-06-14 22:52:19 UTC (rev 7248)
@@ -55,7 +55,7 @@
# boundary conditions
# ----------------------------------------------------------------------
[pylithapp.timedependent.bc.bc]
-fixed_dof = [0]
+fixed_dof = [0, 1, 2]
id = 10
label = end points
db.label = Dirichlet BC
More information about the cig-commits
mailing list