[cig-commits] r15898 - in short/3D/PyLith/branches/pylith-friction: . libsrc/faults libsrc/problems
brad at geodynamics.org
brad at geodynamics.org
Thu Oct 29 16:28:13 PDT 2009
Author: brad
Date: 2009-10-29 16:28:13 -0700 (Thu, 29 Oct 2009)
New Revision: 15898
Modified:
short/3D/PyLith/branches/pylith-friction/TODO
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.hh
short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc
Log:
Initial implementation of extracting Jacobian diagonal.
Modified: short/3D/PyLith/branches/pylith-friction/TODO
===================================================================
--- short/3D/PyLith/branches/pylith-friction/TODO 2009-10-29 15:30:22 UTC (rev 15897)
+++ short/3D/PyLith/branches/pylith-friction/TODO 2009-10-29 23:28:13 UTC (rev 15898)
@@ -4,15 +4,13 @@
Data structure to hold label of fault Lagrange vertices and
conventional vertices
- array = numVertics x 3 (i, j, k)
+ array of structure (lagrange, positive, negative, fault)
Status?
Field split
uniform global refinement for tets with faults
customized SNES line search
- restrictOperator() for sparse matrix
- updateAllVisitor for Mesh
----------------------------------------------------------------------
FRICTION
@@ -48,9 +46,6 @@
integrateJacobian (lumped)
adjustSolnLumped
- ElasticityExplicit
- integrateJacobian (lumped)
-
AbsorbingDampers
integrateJacobian (lumped)
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc 2009-10-29 15:30:22 UTC (rev 15897)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc 2009-10-29 23:28:13 UTC (rev 15898)
@@ -678,7 +678,7 @@
const ALE::Obj<RealSection>& jacobianSection =
_fields->get("Jacobian diagonal").section();
assert(!jacobianSection.isNull());
- _updateJacobianDiagonal(*fields, jacobian);
+ _updateJacobianDiagonal(*fields);
slipSection->view("SLIP");
areaSection->view("AREA");
@@ -1653,8 +1653,7 @@
// associated with Lagrange vertex k.
void
pylith::faults::FaultCohesiveDynL::_updateJacobianDiagonal(
- const topology::SolutionFields& fields,
- const topology::Jacobian& jacobian)
+ const topology::SolutionFields& fields)
{ // _updateJacobianDiagonal
assert(0 != _fields);
@@ -1670,54 +1669,39 @@
cellsCohesive->end();
const int cellsCohesiveSize = cellsCohesive->size();
- const int spaceDim = _quadrature->spaceDim();
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
+ const int spaceDim = _quadrature->spaceDim();
const int numConstraintVert = _quadrature->numBasis();
const int numCorners = 3*numConstraintVert; // cohesive cell
- double_array matrixCell(numCorners*spaceDim * numCorners*spaceDim);
// Get section information
- const ALE::Obj<RealSection>& solutionSection = fields.solution().section();
- assert(!solutionSection.isNull());
-
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
-
- double_array jacobianDiagCell(numConstraintVert*2*spaceDim);
+ double_array jacobianDiagCell(numCorners*spaceDim);
const ALE::Obj<RealSection>& jacobianDiagSection =
+ fields.get("Jacobian diagonal").section();
+ assert(!jacobianDiagSection.isNull());
+ topology::Mesh::RestrictVisitor jacobianDiagVisitor(*jacobianDiagSection,
+ jacobianDiagCell.size(),
+ &jacobianDiagCell[0]);
+
+ double_array jacobianDiagFaultCell(numConstraintVert*2*spaceDim);
+ const ALE::Obj<RealSection>& jacobianDiagFaultSection =
_fields->get("Jacobian diagonal").section();
- topology::Mesh::UpdateAllVisitor jacobianDiagVisitor(*jacobianDiagSection,
- &jacobianDiagCell[0]);
+ topology::Mesh::UpdateAllVisitor jacobianDiagFaultVisitor(*jacobianDiagFaultSection,
+ &jacobianDiagFaultCell[0]);
- const PetscMat jacobianMatrix = jacobian.matrix();
- assert(0 != jacobianMatrix);
- const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
- solutionSection);
- assert(!globalOrder.isNull());
- // We would need to request unique points here if we had an interpolated mesh
- topology::Mesh::IndicesVisitor jacobianVisitor(*solutionSection,
- *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
-
for (SieveMesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
c_iter != cellsCohesiveEnd;
++c_iter) {
const SieveMesh::point_type c_fault = _cohesiveToFault[*c_iter];
- // Insert cell contribution into PETSc Matrix
- jacobianVisitor.clear();
-#if 0 // MISSING
- PetscErrorCode err = restrictOperator(jacobianMatrix,
- *sieveMesh->getSieve(),
- jacobianVisitor, *c_iter,
- &matrixCell[0]);
-#endif
+ jacobianDiagVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, jacobianDiagVisitor);
- for (int iConstraint=0, iJacobian=0;
+ for (int iConstraint=0, indexF=0;
iConstraint < numConstraintVert;
- ++iConstraint, iJacobian += 2*spaceDim) {
+ ++iConstraint) {
// Blocks in cell matrix associated with normal cohesive
// vertices i and j and constraint vertex k
const int indexI = iConstraint;
@@ -1725,25 +1709,25 @@
const int indexK = iConstraint + 2*numConstraintVert;
// Diagonal for vertex i
- int row = indexI*spaceDim;
- int col = row;
- for (int iDim=0; iDim < spaceDim; ++iDim)
- jacobianDiagCell[iJacobian+iDim] =
- matrixCell[(row+iDim)*numCorners*spaceDim+col+iDim];
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ jacobianDiagFaultCell[indexF+iDim] =
+ jacobianDiagCell[indexI*spaceDim+iDim];
+ assert(jacobianDiagFaultCell[indexF+iDim] > 0.0);
+ } // for
+ indexF += spaceDim;
// Diagonal for vertex j
- row = indexJ*spaceDim;
- col = row;
- for (int iDim=0; iDim < spaceDim; ++iDim)
- jacobianDiagCell[iJacobian+spaceDim+iDim] =
- matrixCell[(row+iDim)*numCorners*spaceDim+col+iDim];
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ jacobianDiagFaultCell[indexF+iDim] =
+ jacobianDiagCell[indexJ*spaceDim+iDim];
+ assert(jacobianDiagFaultCell[indexF+iDim] > 0.0);
+ } // for
+ indexF += spaceDim;
} // for
// Insert cell contribution into
jacobianDiagVisitor.clear();
-#if 0 // MISSING
- sieveMesh->updateAll(*c_iter, jacobianDiagVisitor);
-#endif
+ sieveMesh->updateClosure(c_fault, jacobianDiagFaultVisitor);
} // for
} // _updateJacobianDiagonal
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.hh 2009-10-29 15:30:22 UTC (rev 15897)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.hh 2009-10-29 23:28:13 UTC (rev 15898)
@@ -257,10 +257,8 @@
* associated with Lagrange vertex k.
*
* @param fields Solution fields.
- * @param jacobian System jacobian.
*/
- void _updateJacobianDiagonal(const topology::SolutionFields& fields,
- const topology::Jacobian& jacobian);
+ void _updateJacobianDiagonal(const topology::SolutionFields& fields);
/** Compute change in tractions on fault surface using solution.
*
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc 2009-10-29 15:30:22 UTC (rev 15897)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc 2009-10-29 23:28:13 UTC (rev 15898)
@@ -234,6 +234,23 @@
// Assemble jacobian.
_jacobian->assemble("final_assembly");
+
+ // Extract diagonal :KLUDGE: We extract the diagonal even if we
+ // don't need to.
+ if (!_fields->hasField("Jacobian diagonal")) {
+ _fields->add("Jacobian diagonal", "jacobian_diagonal");
+ topology::Field<topology::Mesh>& jacobianDiag =
+ _fields->get("Jacobian diagonal");
+ const topology::Field<topology::Mesh>& disp =
+ _fields->get("disp(t)");
+ jacobianDiag.cloneSection(disp);
+ jacobianDiag.createVector();
+ jacobianDiag.createScatter();
+ } // if
+ topology::Field<topology::Mesh>& jacobianDiag =
+ _fields->get("Jacobian diagonal");
+ MatGetDiagonal(_jacobian->matrix(), jacobianDiag.vector());
+ jacobianDiag.scatterVectorToSection();
} // reformJacobian
// ----------------------------------------------------------------------
More information about the CIG-COMMITS
mailing list