[cig-commits] r15309 - short/3D/PyLith/trunk/libsrc/faults

[brad at geodynamics.org]

Author: brad
Date: 2009-06-16 16:22:38 -0700 (Tue, 16 Jun 2009)
New Revision: 15309

Modified:
   short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
   short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
Log:
In progress bug fix for faults and nonlinear solve.

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc	2009-06-16 23:21:56 UTC (rev 15308)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc	2009-06-16 23:22:38 UTC (rev 15309)
@@ -245,6 +245,8 @@
   // Allocate vectors for cell values
   double_array orientationCell(numConstraintVert*orientationSize);
   double_array dispTCell(numCorners*spaceDim);
+  double_array dispTIncrCell(numCorners*spaceDim);
+  double_array dispTpdtCell(numBasis*spaceDim);
   double_array residualCell(numCorners*spaceDim);
 
   // Tributary area for the current for each vertex.
@@ -290,6 +292,15 @@
 					       dispTCell.size(), 
 					       &dispTCell[0]);
 
+#if 0 // :TODO: Need to use solution
+  topology::Field<topology::Mesh>& dispTIncr = fields->get("dispIncr(t->t+dt)");
+  const ALE::Obj<RealSection>& dispTIncrSection = dispTIncr.section();
+  assert(!dispTIncrSection.isNull());  
+  topology::Mesh::RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
+					       dispTIncrCell.size(), 
+					       &dispTIncrCell[0]);
+#endif
+
   const ALE::Obj<RealSection>& residualSection = residual.section();
   topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
 						   &residualCell[0]);
@@ -342,6 +353,16 @@
     dispTVisitor.clear();
     sieveMesh->restrictClosure(*c_iter, dispTVisitor);
     
+#if 0 // :TODO: need to use solution
+    // Get dispIncr(t) at cohesive cell's vertices.
+    dispTIncrVisitor.clear();
+    sieveMesh->restrictClosure(*c_iter, dispTIncrVisitor);
+#endif
+    
+    // Compute current estimate of displacement at time t+dt using
+    // solution increment.
+    dispTpdtCell = dispTCell + dispTIncrCell;
+
     for (int iConstraint=0; iConstraint < numConstraintVert; ++iConstraint) {
       // Blocks in cell matrix associated with normal cohesive
       // vertices i and j and constraint vertex k
@@ -365,13 +386,25 @@
 	    -orientationVertex[kDim*spaceDim+iDim] * wt;
       } // for
       
-	// Entries associated with constraint forces applied at node j
+      // Entries associated with constraint forces applied at node j
       for (int jDim=0; jDim < spaceDim; ++jDim) {
 	for (int kDim=0; kDim < spaceDim; ++kDim)
 	  residualCell[indexJ*spaceDim+jDim] -=
 	    dispTCell[indexK*spaceDim+kDim] * 
 	    orientationVertex[kDim*spaceDim+jDim] * wt;
       } // for
+      
+#if 0 // :TODO: Is this a missing -C u term??
+      // Entries associated with relative displacements between node i
+      // and node j for constraint node k
+      for (int kDim=0; kDim < spaceDim; ++kDim) {
+	for (int iDim=0; iDim < spaceDim; ++iDim)
+	  residualCell[indexK*spaceDim+kDim] -=
+	    (dispTpdtCell[indexI*spaceDim+iDim] - 
+	     dispTpdtCell[indexJ*spaceDim+iDim]) *
+	    orientationVertex[kDim*spaceDim+iDim] * wt;
+      } // for
+#endif
     } // for
 
 #if 0 // DEBUGGING
@@ -383,7 +416,7 @@
       std::cout << "  slip["<<i<<"]: " << cellSlip[i] << std::endl;
     }
     for(int i = 0; i < numCorners*spaceDim; ++i) {
-      std::cout << "  soln["<<i<<"]: " << dispTCell[i] << std::endl;
+      std::cout << "  soln["<<i<<"]: " << dispTpdtCell[i] << std::endl;
     }
     for(int i = 0; i < numCorners*spaceDim; ++i) {
       std::cout << "  v["<<i<<"]: " << residualCell[i] << std::endl;

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh	2009-06-16 23:21:56 UTC (rev 15308)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh	2009-06-16 23:22:38 UTC (rev 15309)
@@ -24,14 +24,13 @@
  * one side of the fault, the corresponding entries on the other side
  * of the fault, and then the corresponding constraint vertices.
  *
- * [ K  aC^T ] [ U    ] = [ Fe ]
- * [ C   0   ] [ Fi/a ] = [ D  ]
+ * [ K   C^T ] [ U  ] = [ Fe ]
+ * [ C   0   ] [ Fi ] = [ D  ]
  *
  * where K is the stiffness matrix, C is the matrix of Lagrange
  * constraints, U is the displacement field, Fe is the vector of
  * external forces, Fi is the vector of Lagrange multipers (forces), D
- * is the fault slip, and "a" is the conditioning value (taken to be
- * the shear modulus).
+ * is the fault slip.
  */
 
 #if !defined(pylith_faults_faultcohesivekin_hh)