[cig-commits] r19042 - in short/3D/PyLith/branches/v1.6-revisedfault: libsrc/pylith/faults unittests/libtests/faults unittests/libtests/faults/data

brad at geodynamics.org brad at geodynamics.org
Fri Oct 7 10:38:43 PDT 2011


Author: brad
Date: 2011-10-07 10:38:42 -0700 (Fri, 07 Oct 2011)
New Revision: 19042

Modified:
   short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveDyn.cc
   short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc
   short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveDyn.cc
   short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
   short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc
   short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
   short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
   short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc
   short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/cohesivedyn.py
Log:
Updated spontaneous rupture implementation for collocated quadrature points and cells vertices. Updated corresponding test data.

Modified: short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveDyn.cc	2011-10-07 16:25:42 UTC (rev 19041)
+++ short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveDyn.cc	2011-10-07 17:38:42 UTC (rev 19042)
@@ -177,215 +177,126 @@
   _logger->eventBegin(setupEvent);
 
   // Get cell geometry information that doesn't depend on cell
-  const int numQuadPts = _quadrature->numQuadPts();
-  const double_array& quadWts = _quadrature->quadWts();
-  assert(quadWts.size() == numQuadPts);
-  const int numBasis = _quadrature->numBasis();
   const int spaceDim = _quadrature->spaceDim();
-  const int cellDim = _quadrature->cellDim();
-  assert(cellDim == spaceDim-1);
 
-  // Get cohesive cell information
-  const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
-  assert(!sieveMesh.isNull());
-  const ALE::Obj<SieveMesh::label_sequence>& cells = 
-    sieveMesh->getLabelStratum("material-id", id());
-  assert(!cells.isNull());
-  const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
-  const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
-  const int numCells = cells->size();
-
   // Get sections associated with cohesive cells
-  double_array residualCell(3*numBasis*spaceDim);
+  double_array residualVertexN(spaceDim);
+  double_array residualVertexP(spaceDim);
   const ALE::Obj<RealSection>& residualSection = residual.section();
   assert(!residualSection.isNull());
-  UpdateAddVisitor residualVisitor(*residualSection, &residualCell[0]);
 
-  // Get fault cell information
-  const ALE::Obj<SieveMesh>& faultSieveMesh = _faultMesh->sieveMesh();
-  assert(!faultSieveMesh.isNull());
-  const ALE::Obj<SieveSubMesh::label_sequence>& faultCells =
-    faultSieveMesh->heightStratum(0);
-  assert(!faultCells.isNull());
-  assert(faultCells->size() == cells->size());
+  const ALE::Obj<RealSection>& areaSection = _fields->get("area").section();
+  assert(!areaSection.isNull());
 
-  // Get sections associated with fault cells
-  double_array coordinatesCell(numBasis*spaceDim);
-  const ALE::Obj<RealSection>& coordinates = 
-    faultSieveMesh->getRealSection("coordinates");
-  assert(!coordinates.isNull());
-  RestrictVisitor coordsVisitor(*coordinates, 
-				coordinatesCell.size(), &coordinatesCell[0]);
+  const ALE::Obj<RealSection>& initialTractionsSection = 
+    _fields->get("initial traction").section();
+  assert(!initialTractionsSection.isNull());
 
-  double_array dispRelCell(numBasis*spaceDim);
   const ALE::Obj<RealSection>& dispRelSection = 
     _fields->get("relative disp").section();
   assert(!dispRelSection.isNull());
-  RestrictVisitor dispRelVisitor(*dispRelSection, 
-				 dispRelCell.size(), &dispRelCell[0]);
 
-  double_array initialTractionsCell(numBasis*spaceDim);
-  const ALE::Obj<RealSection>& initialTractionsSection = 
-    _fields->get("initial traction").section();
-  assert(!initialTractionsSection.isNull());
-  RestrictVisitor initialTractionsVisitor(*initialTractionsSection, 
-					  initialTractionsCell.size(),
-					  &initialTractionsCell[0]);
-
-  double_array orientationCell(numBasis*spaceDim*spaceDim);
   const ALE::Obj<RealSection>& orientationSection = 
     _fields->get("orientation").section();
   assert(!orientationSection.isNull());
-  RestrictVisitor orientationVisitor(*orientationSection, 
-				     orientationCell.size(),
-				     &orientationCell[0]);
 
+  // Get fault information
+  const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
+  assert(!sieveMesh.isNull());
+  const ALE::Obj<SieveMesh::order_type>& globalOrder =
+      sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
+					      residualSection);
+  assert(!globalOrder.isNull());
+
   _logger->eventEnd(setupEvent);
 #if !defined(DETAILED_EVENT_LOGGING)
   _logger->eventBegin(computeEvent);
 #endif
 
-  // Loop over cells
-  for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
-	 c_iter != cellsEnd;
-       ++c_iter) {
-    topology::SubMesh::SieveMesh::point_type c_fault = 
-      _cohesiveToFault[*c_iter];
+  // Loop over fault vertices
+  const int numVertices = _cohesiveVertices.size();
+  for (int iVertex=0; iVertex < numVertices; ++iVertex) {
+    const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
+    const int v_fault = _cohesiveVertices[iVertex].fault;
+    const int v_negative = _cohesiveVertices[iVertex].negative;
+    const int v_positive = _cohesiveVertices[iVertex].positive;
 
-    // Compute geometry information for current cell
-#if defined(DETAILED_EVENT_LOGGING)
-    _logger->eventBegin(geometryEvent);
-#endif
-#if defined(PRECOMPUTE_GEOMETRY)
-    _quadrature->retrieveGeometry(c_fault);
-#else
-    coordsVisitor.clear();
-    faultSieveMesh->restrictClosure(c_fault, coordsVisitor);
-    _quadrature->computeGeometry(coordinatesCell, c_fault);
-#endif
+    // Compute contribution only if Lagrange constraint is local.
+    if (!globalOrder->isLocal(v_lagrange))
+      continue;
 
 #if defined(DETAILED_EVENT_LOGGING)
-    _logger->eventEnd(geometryEvent);
     _logger->eventBegin(restrictEvent);
 #endif
 
-    // Restrict input fields to cell
-    initialTractionsVisitor.clear();
-    faultSieveMesh->restrictClosure(c_fault, initialTractionsVisitor);
+    // Get initial tractions at fault vertex.
+    assert(spaceDim == initialTractionsSection->getFiberDimension(v_fault));
+    const double* initialTractionsVertex = 
+      initialTractionsSection->restrictPoint(v_fault);
+    assert(initialTractionsVertex);
 
-    dispRelVisitor.clear();
-    faultSieveMesh->restrictClosure(c_fault, dispRelVisitor);
+    // Get relative dislplacement at fault vertex.
+    assert(spaceDim == dispRelSection->getFiberDimension(v_fault));
+    const double* dispRelVertex = dispRelSection->restrictPoint(v_fault);
+    assert(dispRelVertex);
 
-    orientationVisitor.clear();
-    faultSieveMesh->restrictClosure(c_fault, orientationVisitor);
+    // Get area associated with fault vertex.
+    assert(1 == areaSection->getFiberDimension(v_fault));
+    assert(areaSection->restrictPoint(v_fault));
+    const double areaVertex = *areaSection->restrictPoint(v_fault);
 
-    // Get cell geometry information that depends on cell
-    const double_array& basis = _quadrature->basis();
-    const double_array& basisDeriv = _quadrature->basisDeriv();
-    const double_array& jacobianDet = _quadrature->jacobianDet();
+    // Get orientation associated with fault vertex.
+    assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
+    const double* orientationVertex = 
+      orientationSection->restrictPoint(v_fault);
+    assert(orientationVertex);
 
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventEnd(restrictEvent);
     _logger->eventBegin(computeEvent);
 #endif
 
-    residualCell = 0.0;
-
-    // Compute action for positive side of fault and Lagrange constraints
-    for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
-      const double wt = quadWts[iQuad] * jacobianDet[iQuad];
-      const int iQ = iQuad * numBasis;
-      for (int iBasis=0; iBasis < numBasis; ++iBasis) {
-        const double valI = wt*basis[iQ+iBasis];
-
-	// Add entries to residual at
-	// iN = DOF at negative node
-	// iP = DOF at positive node
-	// iL = DOF at constraint node
-
-	// Indices for negative vertex
-	const int iBN = 0*numBasis*spaceDim + iBasis*spaceDim;
-	
-	// Indices for positive vertex
-	const int iBP = 1*numBasis*spaceDim + iBasis*spaceDim;
-	
-        for (int jBasis=0; jBasis < numBasis; ++jBasis) {
-          const double valIJ = valI * basis[iQ+jBasis];
-
-	  // Indices for fault vertex
-	  const int jB = jBasis*spaceDim;
-
-          for (int iDim=0; iDim < spaceDim; ++iDim) {
-	    // Initial (external) tractions oppose (internal)
-	    // tractions associated with Lagrange multiplier, so these
-	    // terms have the opposite sign as the integration of the
-	    // Lagrange multipliers in FaultCohesiveLagrange.
-
-	    // negative side of the fault
-            residualCell[iBN + iDim] -= valIJ * initialTractionsCell[jB + iDim];
-	    
-	    // positive side of the fault
-            residualCell[iBP + iDim] += valIJ * initialTractionsCell[jB + iDim];
-	    
-#if 0
-	    std::cout << "iBasis: " << iBasis
-		      << ", jBasis: " << jBasis
-		      << ", iDim: " << iDim
-		      << ", valIJ: " << valIJ
-		      << ", jacobianDet: " << jacobianDet[iQuad]
-		      << ", initialTractions: " << initialTractionsCell[jB + iDim]
-		      << ", residualN: " << residualCell[iBN + iDim]
-		      << ", residualP: " << residualCell[iBP + iDim]
-		      << std::endl;
-#endif
-
-	  } // for
-        } // for
-      } // for
+    // Initial (external) tractions oppose (internal) tractions
+    // associated with Lagrange multiplier, so these terms have the
+    // opposite sign as the integration of the Lagrange multipliers in
+    // FaultCohesiveLagrange.
+    for (int iDim=0; iDim < spaceDim; ++iDim) {
+      residualVertexP[iDim] = areaVertex * initialTractionsVertex[iDim];
     } // for
+    residualVertexN = -residualVertexP;
 
-
     // Only apply initial tractions if there is no opening.
     // If there is opening, zero out initial tractions
-    for (int iBasis=0; iBasis < numBasis; ++iBasis) {
-      const int iB = iBasis*spaceDim;
-      const int iO = iBasis*spaceDim*spaceDim;
-      
-      double slipNormal = 0.0;
-      const int indexN = spaceDim - 1;
-      for (int jDim=0; jDim < spaceDim; ++jDim) {
-	slipNormal += 
-	  orientationCell[iO + indexN*spaceDim+jDim]*dispRelCell[iB+jDim];
-      } // for
-
-      if (slipNormal > _zeroTolerance) {
-	// Indices for negative vertex
-	const int iBN = 0*numBasis*spaceDim + iBasis*spaceDim;
-	
-	// Indices for positive vertex
-	const int iBP = 1*numBasis*spaceDim + iBasis*spaceDim;
-	
-	for (int iDim=0; iDim < spaceDim; ++iDim) {
-	  residualCell[iBN+iDim] = 0.0;
-	  residualCell[iBP+iDim] = 0.0;
-	} // for
-      } // if
+    double slipNormal = 0.0;
+    const int indexN = spaceDim - 1;
+    for (int jDim=0; jDim < spaceDim; ++jDim) {
+      slipNormal += orientationVertex[indexN*spaceDim+jDim]*dispRelVertex[jDim];
     } // for
 
+    if (slipNormal > _zeroTolerance) {
+      residualVertexN = 0.0;
+      residualVertexP = 0.0;
+    } // if
+
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventEnd(computeEvent);
     _logger->eventBegin(updateEvent);
 #endif
 
-    // Assemble cell contribution into field
-    residualVisitor.clear();
-    sieveMesh->updateClosure(*c_iter, residualVisitor);
+    // Assemble contributions into field
+    assert(residualVertexN.size() == 
+	   residualSection->getFiberDimension(v_negative));
+    residualSection->updateAddPoint(v_negative, &residualVertexN[0]);
 
+    assert(residualVertexP.size() == 
+	   residualSection->getFiberDimension(v_positive));
+    residualSection->updateAddPoint(v_positive, &residualVertexP[0]);
+
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventEnd(updateEvent);
 #endif
   } // for
-  PetscLogFlops(numCells*spaceDim*spaceDim*8);
+  PetscLogFlops(numVertices*spaceDim*(2+spaceDim*2));
 
 #if !defined(DETAILED_EVENT_LOGGING)
   _logger->eventEnd(computeEvent);

Modified: short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2011-10-07 16:25:42 UTC (rev 19041)
+++ short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2011-10-07 17:38:42 UTC (rev 19042)
@@ -213,7 +213,6 @@
   _logger->eventBegin(setupEvent);
 
   // Get cell geometry information that doesn't depend on cell
-  const int numBasis = _quadrature->numBasis();
   const int spaceDim = _quadrature->spaceDim();
 
   // Get sections associated with cohesive cells
@@ -500,9 +499,9 @@
 #endif
 
   } // for
+  PetscLogFlops(numVertices*spaceDim*2);
 
 #if !defined(DETAILED_EVENT_LOGGING)
-  PetscLogFlops(numVertices*spaceDim*2);
   _logger->eventEnd(computeEvent);
 #endif
 
@@ -722,7 +721,6 @@
     
 
 #if defined(DETAILED_EVENT_LOGGING)
-    PetscLogFlops(spaceDim*6);
     _logger->eventEnd(computeEvent);
     _logger->eventBegin(updateEvent);
 #endif
@@ -752,9 +750,9 @@
 #endif
   } // for
   err = MatDestroy(&jacobianNP); CHECK_PETSC_ERROR(err);
+  PetscLogFlops(numVertices*spaceDim*6);
 
 #if !defined(DETAILED_EVENT_LOGGING)
-    PetscLogFlops(numVertices*spaceDim*6);
     _logger->eventEnd(computeEvent);
 #endif
 } // calcPreconditioner

Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveDyn.cc	2011-10-07 16:25:42 UTC (rev 19041)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveDyn.cc	2011-10-07 17:38:42 UTC (rev 19042)
@@ -737,6 +737,8 @@
   disp.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
   disp.allocate();
   fields->copyLayout("disp(t)");
+
+  fault->verifyConfiguration(*mesh);
 } // _initialize
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataHex8.cc	2011-10-07 16:25:42 UTC (rev 19041)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataHex8.cc	2011-10-07 17:38:42 UTC (rev 19042)
@@ -49,10 +49,10 @@
 const int pylith::faults::CohesiveDynDataHex8::_numQuadPts = 4;
 
 const double pylith::faults::CohesiveDynDataHex8::_quadPts[] = {
-  -0.57735027, -0.57735027,
-  +0.57735027, -0.57735027,
-  +0.57735027, +0.57735027,
-  -0.57735027, +0.57735027,
+  -1.0, -1.0,
+  +1.0, -1.0,
+  +1.0, +1.0,
+  -1.0, +1.0
 };
 
 const double pylith::faults::CohesiveDynDataHex8::_quadWts[] = {
@@ -60,10 +60,10 @@
 };
 
 const double pylith::faults::CohesiveDynDataHex8::_basis[] = {
-  0.62200847,  0.16666667,  0.0446582,   0.16666667,
-  0.16666667,  0.62200847,  0.16666667,   0.0446582,
-  0.0446582,   0.16666667,  0.62200847,  0.16666667,
-  0.16666667,   0.0446582,  0.16666667,  0.62200847,
+  1.0, 0.0, 0.0, 0.0,
+  0.0, 1.0, 0.0, 0.0,
+  0.0, 0.0, 1.0, 0.0,
+  0.0, 0.0, 0.0, 1.0,
 };
 
 const double pylith::faults::CohesiveDynDataHex8::_basisDeriv[] = {
@@ -1174,18 +1174,18 @@
    1.200000000000,   2.200000000000,   0.200000000000,
    1.300000000000,   2.300000000000,   0.300000000000,
    1.400000000000,   2.400000000000,   0.400000000000,
-   1.500000000000,   3.738231538218,   1.301744825995,
-   1.600000000000,   3.315253904282,   1.405182104855,
-   1.700000000000,   3.631013149063,   1.616032753440,
-   1.800000000000,   3.652003834523,   1.629677922638,
+   1.500000000000,   3.745293503266,   1.312227205500,
+   1.600000000000,   3.294980394575,   1.386664200491,
+   1.700000000000,   3.591763265955,   1.582543631245,
+   1.800000000000,   3.703107773803,   1.671356663600,
    1.900000000000,   2.900000000000,   0.900000000000,
    1.000000000000,   2.000000000000,   0.000000000000,
    1.100000000000,   2.100000000000,   0.100000000000,
    1.200000000000,   2.200000000000,   0.200000000000,
-   1.300000000000,   1.061768461782,  -0.501744825995,
-   1.500000000000,   1.784746095718,  -0.305182104855,
-   1.700000000000,   1.768986850937,  -0.216032753440,
-   1.900000000000,   2.047996165477,   0.070322077362,
+   1.300000000000,   1.054706496734,  -0.512227205500,
+   1.500000000000,   1.805019605425,  -0.286664200491,
+   1.700000000000,   1.808236734045,  -0.182543631245,
+   1.900000000000,   1.996892226197,   0.028643336400,
    1.400000000000,   0.328479469127,  -1.239953753608,
    1.600000000000,   0.293941190358,  -1.262585961634,
    1.800000000000,   0.259995967920,  -1.286431883494,
@@ -1193,10 +1193,10 @@
 };
 
 const double pylith::faults::CohesiveDynDataHex8::_slipSlipE[] = {
-   2.476463076435,  -1.603489651991,   0.000000000000,
-   1.430507808565,  -1.610364209711,   0.000000000000,
-   1.862026298126,  -1.832065506881,   0.000000000000,
-   1.704007669046,  -1.659355845275,   0.000000000000,
+   2.490587006532,  -1.624454410999,   0.000000000000,
+   1.389960789150,  -1.573328400983,   0.000000000000,
+   1.783526531910,  -1.765087262491,   0.000000000000,
+   1.806215547605,  -1.742713327201,   0.000000000000,
 };
 
 // ----------------------------------------------------------------------
@@ -1232,18 +1232,18 @@
    1.200000000000,   2.200000000000,   0.200000000000,
    1.300000000000,   2.300000000000,   0.300000000000,
    1.400000000000,   2.400000000000,   0.400000000000,
-   3.486516791199,   3.439968915124,   1.042649773060,
-   3.314817860545,   3.003266771802,   1.131804624352,
-   3.634863362999,   3.321964812539,   1.325307625637,
-   3.486832213768,   3.146857179684,   1.155004889953,
+   3.473386524898,   3.450458773757,   1.061077907335,
+   3.302801176428,   2.966499592247,   1.097759351340,
+   3.627558189277,   3.252968860584,   1.264881344260,
+   3.475698668581,   3.241605816971,   1.233144951440,
    1.900000000000,   2.900000000000,   0.900000000000,
    1.000000000000,   2.000000000000,   0.000000000000,
    1.100000000000,   2.100000000000,   0.100000000000,
    1.200000000000,   2.200000000000,   0.200000000000,
-  -0.686516791199,   1.360031084876,  -0.242649773060,
-  -0.214817860545,   2.096733228198,  -0.031804624352,
-  -0.234863362999,   2.078035187461,   0.074692374363,
-   0.213167786232,   2.553142820316,   0.544995110047,
+  -0.673386524898,   1.349541226243,  -0.261077907335,
+  -0.202801176428,   2.133500407753,   0.002240648660,
+  -0.227558189277,   2.147031139416,   0.135118655740,
+   0.224301331419,   2.458394183029,   0.466855048560,
   -4.400000000000,  -2.400000000000,  -3.400000000000,
   -4.600000000000,  -2.600000000000,  -3.600000000000,
   -4.800000000000,  -2.800000000000,  -3.800000000000,
@@ -1251,10 +1251,10 @@
 };
 
 const double pylith::faults::CohesiveDynDataHex8::_slipOpenE[] = {
-   1.879937830249,  -1.085299546121,   3.973033582398,
-   0.806533543603,  -1.063609248704,   3.429635721090,
-   1.243929625077,  -1.250615251274,   3.869726725998,
-   0.693714359369,  -0.710009779906,   3.373664427535,
+   1.900917547514,  -1.122155814671,   3.946773049797,
+   0.732999184495,  -0.995518702680,   3.405602352856,
+   1.105937721169,  -1.129762688519,   3.855116378553,
+   0.883211633942,  -0.866289902880,   3.351397337162,
 };
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc	2011-10-07 16:25:42 UTC (rev 19041)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc	2011-10-07 17:38:42 UTC (rev 19042)
@@ -59,24 +59,24 @@
 
 const int pylith::faults::CohesiveDynDataQuad4::_numBasis = 2;
 
-const int pylith::faults::CohesiveDynDataQuad4::_numQuadPts = 1;
+const int pylith::faults::CohesiveDynDataQuad4::_numQuadPts = 2;
 
 const double pylith::faults::CohesiveDynDataQuad4::_quadPts[] = {
-  0.0,
+  -1.0, 1.0,
 };
 
 const double pylith::faults::CohesiveDynDataQuad4::_quadWts[] = {
-  2.0,
+  1.0, 1.0
 };
 
 const double pylith::faults::CohesiveDynDataQuad4::_basis[] = {
-  0.5,
-  0.5
+  1.0, 0.0,
+  0.0, 1.0,
 };
 
 const double pylith::faults::CohesiveDynDataQuad4::_basisDeriv[] = {
-  -0.5,
-   0.5
+  -0.5, 0.5,
+  -0.5, 0.5,
 };
 
 const double pylith::faults::CohesiveDynDataQuad4::_verticesRef[] = {
@@ -332,20 +332,19 @@
 const double pylith::faults::CohesiveDynDataQuad4::_fieldIncrSlipE[] = {
    1.100000000000,   2.100000000000,
    1.200000000000,   2.200000000000,
-   1.300000000000,   2.468370545405,
-   1.400000000000,   3.350253402750,
+   1.300000000000,   2.406970297611,
+   1.400000000000,   3.369389712309,
    1.500000000000,   2.500000000000,
    1.600000000000,   2.600000000000,
-   1.700000000000,   2.531629454595,
-   1.900000000000,   1.949746597250,
+   1.700000000000,   2.593029702389,
+   1.900000000000,   1.930610287691,
    1.800000000000,  -3.440000000000,
    1.000000000000,  -3.600000000000,
 };
 
-// Update slip values based on changes in Lagrange multiplier values
 const double pylith::faults::CohesiveDynDataQuad4::_slipSlipE[] = {
-   0.336741090809,   0.000000000000,
-   1.900506805500,   0.000000000000,
+   0.213940595221,   0.000000000000,
+   1.938779424618,   0.000000000000,
 };
 
 // ----------------------------------------------------------------------
@@ -369,19 +368,19 @@
 const double pylith::faults::CohesiveDynDataQuad4::_fieldIncrOpenE[] = {
    1.100000000000,   2.100000000000,
    1.200000000000,   2.200000000000,
-   3.837558167495,   2.192904776328,
-   4.297022233334,   3.773022694556,
+   3.815834136232,   2.039404156843,
+   4.287400342924,   3.820863468454,
    1.500000000000,   2.500000000000,
    1.600000000000,   2.600000000000,
-  -0.837558167495,   2.807095223672,
-  -0.997022233334,   1.526977305444,
+  -0.815834136232,   2.960595843157,
+  -0.987400342924,   1.479136531546,
   -8.800000000000,  -9.800000000000,
   -8.000000000000,  -9.000000000000,
 };
 
 const double pylith::faults::CohesiveDynDataQuad4::_slipOpenE[] = {
-  -0.214190447343,   5.075116334990,
-   2.746045389111,   5.794044466667,
+  -0.521191686313,   5.031668272463,
+   2.841726936907,   5.774800685848,
 };
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTet4.cc	2011-10-07 16:25:42 UTC (rev 19041)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTet4.cc	2011-10-07 17:38:42 UTC (rev 19042)
@@ -45,24 +45,34 @@
 
 const int pylith::faults::CohesiveDynDataTet4::_numBasis = 3;
 
-const int pylith::faults::CohesiveDynDataTet4::_numQuadPts = 1;
+const int pylith::faults::CohesiveDynDataTet4::_numQuadPts = 3;
 
 const double pylith::faults::CohesiveDynDataTet4::_quadPts[] = {
-  -3.33333333e-01,  -3.33333333e-01,
+ -1.00000000e+00, -1.00000000e+00,
+  1.00000000e+00, -1.00000000e+00,
+ -1.00000000e+00,  1.00000000e+00,
 };
 
 const double pylith::faults::CohesiveDynDataTet4::_quadWts[] = {
-  2.0,
+  2.0/3.0, 2.0/3.0, 2.0/3.0,
 };
 
 const double pylith::faults::CohesiveDynDataTet4::_basis[] = {
-  3.33333333e-01,  3.33333333e-01,
-  3.33333333e-01,};
+  1.0, 0.0, 0.0,
+  0.0, 1.0, 0.0,
+  0.0, 0.0, 1.0,
+};
 
 const double pylith::faults::CohesiveDynDataTet4::_basisDeriv[] = {
  -0.50000000e+00, -0.50000000e+00,
   0.50000000e+00,  0.00000000e+00,
   0.00000000e+00,  0.50000000e+00,
+ -0.50000000e+00, -0.50000000e+00,
+  0.50000000e+00,  0.00000000e+00,
+  0.00000000e+00,  0.50000000e+00,
+ -0.50000000e+00, -0.50000000e+00,
+  0.50000000e+00,  0.00000000e+00,
+  0.00000000e+00,  0.50000000e+00,
 };
 
 const double pylith::faults::CohesiveDynDataTet4::_verticesRef[] = {
@@ -528,22 +538,22 @@
 // Output
 const double pylith::faults::CohesiveDynDataTet4::_fieldIncrSlipE[] = {
    1.100000000000,   2.100000000000,   3.100000000000,
-   1.200000000000,   2.433903699854,   3.235410526853,
-   1.300000000000,   2.308933719863,   3.314292166303,
-   1.400000000000,   2.517316497418,   3.471580524964,
+   1.200000000000,   2.427020189783,   3.226804829918,
+   1.300000000000,   2.281534813168,   3.292288303722,
+   1.400000000000,   2.543753403100,   3.499170314087,
    1.500000000000,   2.500000000000,   3.500000000000,
-   1.600000000000,   2.366096300146,   3.564589473147,
-   1.800000000000,   2.791066280137,   3.785707833697,
-   1.000000000000,   1.882683502582,   2.928419475036,
+   1.600000000000,   2.372979810217,   3.573195170082,
+   1.800000000000,   2.818465186832,   3.807711696278,
+   1.000000000000,   1.856246596900,   2.900829685913,
    9.700000000000,  -1.935106067360,  -1.748282570405,
    9.900000000000,  -1.959241090539,  -1.782841275376,
    9.100000000000,  -1.865391385620,  -1.642919108291,
 };
 
 const double pylith::faults::CohesiveDynDataTet4::_slipSlipE[] = {
-   0.467807399707,  -0.070821053706,   0.000000000000,
-   0.017867439727,  -0.028584332606,   0.000000000000,
-   0.234632994837,  -0.143161049929,   0.000000000000,
+   0.454040379566,  -0.053609659837,   0.000000000000,
+  -0.036930373664,   0.015423392555,   0.000000000000,
+   0.287506806199,  -0.198340628174,   0.000000000000,
 };
 
 // ----------------------------------------------------------------------
@@ -567,22 +577,22 @@
 // Output
 const double pylith::faults::CohesiveDynDataTet4::_fieldIncrOpenE[] = {
    1.100000000000,   2.100000000000,   3.100000000000,
-   2.475088478561,   2.213019936374,   2.881149654501,
-   2.516111596990,   1.955359213432,   2.906382357903,
-   2.454638633592,   2.326073874055,   3.211721621635,
+   2.478935676274,   2.200668834288,   2.866497129446,
+   2.514833997358,   1.906576121577,   2.868854721277,
+   2.460386652668,   2.373508489140,   3.258697327233,
    1.500000000000,   2.500000000000,   3.500000000000,
-   0.324911521439,   2.586980063626,   3.918850345499,
-   0.583888403010,   3.144640786568,   4.193617642097,
-  -0.054638633592,   2.073926125945,   3.188278378365,
+   0.321064323726,   2.599331165712,   3.933502870554,
+   0.585166002642,   3.193423878423,   4.231145278723,
+  -0.060386652668,   2.026491510860,   3.141302672767,
   -7.700000000000,  -8.700000000000,  -9.700000000000,
   -7.900000000000,  -8.900000000000,  -9.900000000000,
   -7.100000000000,  -8.100000000000,  -9.100000000000,
 };
 
 const double pylith::faults::CohesiveDynDataTet4::_slipOpenE[] = {
-   0.026039872749,   0.637700690999,   2.550176957122,
-  -0.689281573136,   0.787235284195,   2.432223193980,
-  -0.147852251890,   0.376556756729,   2.109277267185,
+   0.001337668577,   0.667005741108,   2.557871352547,
+  -0.786847756846,   0.862290557446,   2.429667994716,
+  -0.052983021720,   0.282605345535,   2.120773305336,
 };
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3.cc	2011-10-07 16:25:42 UTC (rev 19041)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3.cc	2011-10-07 17:38:42 UTC (rev 19042)
@@ -61,24 +61,24 @@
 
 const int pylith::faults::CohesiveDynDataTri3::_numBasis = 2;
 
-const int pylith::faults::CohesiveDynDataTri3::_numQuadPts = 1;
+const int pylith::faults::CohesiveDynDataTri3::_numQuadPts = 2;
 
 const double pylith::faults::CohesiveDynDataTri3::_quadPts[] = {
-  0.0,
+  -1.0, 1.0,
 };
 
 const double pylith::faults::CohesiveDynDataTri3::_quadWts[] = {
-  2.0,
+  1.0, 1.0
 };
 
 const double pylith::faults::CohesiveDynDataTri3::_basis[] = {
-  0.5,
-  0.5
+  1.0, 0.0,
+  0.0, 1.0,
 };
 
 const double pylith::faults::CohesiveDynDataTri3::_basisDeriv[] = {
-  -0.5,
-   0.5
+  -0.5, 0.5,
+  -0.5, 0.5,
 };
 
 const double pylith::faults::CohesiveDynDataTri3::_verticesRef[] = {
@@ -297,18 +297,18 @@
 // Output
 const double pylith::faults::CohesiveDynDataTri3::_fieldIncrSlipE[] = {
    9.100000000000,   7.100000000000,
-   9.200000000000,   7.375196378327,
-   9.300000000000,   8.288777189348,
+   9.200000000000,   7.390546440275,
+   9.300000000000,   8.283993111958,
    9.400000000000,   7.400000000000,
-   9.500000000000,   7.324803621673,
-   9.700000000000,   6.711222810652,
+   9.500000000000,   7.309453559725,
+   9.700000000000,   6.716006888042,
    1.600000000000,  -3.480000000000,
    1.800000000000,  -3.440000000000,
 };
 
 const double pylith::faults::CohesiveDynDataTri3::_slipSlipE[] = {
-   0.350392756653,   0.000000000000,
-   1.977554378695,   0.000000000000,
+   0.381092880550,   0.000000000000,
+   1.967986223916,   0.000000000000,
 };
 
 // ----------------------------------------------------------------------
@@ -329,18 +329,18 @@
 // Output
 const double pylith::faults::CohesiveDynDataTri3::_fieldIncrOpenE[] = {
    9.100000000000,   7.100000000000,
-  11.868906669946,   7.109969358633,
-  12.354802377535,   8.769332161050,
+  11.874337677762,   7.148344513504,
+  12.357207850137,   8.757371967576,
    9.400000000000,   7.400000000000,
-   6.831093330054,   7.590030641367,
-   6.645197622465,   6.230667838950,
+   6.825662322238,   7.551655486496,
+   6.642792149863,   6.242628032424,
   -8.600000000000,  -9.600000000000,
   -8.800000000000,  -9.800000000000,
 };
 
 const double pylith::faults::CohesiveDynDataTri3::_slipOpenE[] = {
-  -0.180061282734,   5.337813339891,
-   2.938664322101,   6.109604755069,
+  -0.103310972991,   5.348675355523,
+   2.914743935152,   6.114415700274,
 };
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc	2011-10-07 16:25:42 UTC (rev 19041)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc	2011-10-07 17:38:42 UTC (rev 19042)
@@ -76,24 +76,24 @@
 
 const int pylith::faults::CohesiveDynDataTri3d::_numBasis = 2;
 
-const int pylith::faults::CohesiveDynDataTri3d::_numQuadPts = 1;
+const int pylith::faults::CohesiveDynDataTri3d::_numQuadPts = 2;
 
 const double pylith::faults::CohesiveDynDataTri3d::_quadPts[] = {
-  0.0,
+  -1.0, 1.0,
 };
 
 const double pylith::faults::CohesiveDynDataTri3d::_quadWts[] = {
-  2.0,
+  1.0, 1.0
 };
 
 const double pylith::faults::CohesiveDynDataTri3d::_basis[] = {
-  0.5,
-  0.5
+  1.0, 0.0,
+  0.0, 1.0,
 };
 
 const double pylith::faults::CohesiveDynDataTri3d::_basisDeriv[] = {
-  -0.5,
-   0.5
+  -0.5, 0.5,
+  -0.5, 0.5,
 };
 
 const double pylith::faults::CohesiveDynDataTri3d::_verticesRef[] = {
@@ -487,23 +487,23 @@
 // Output
 const double pylith::faults::CohesiveDynDataTri3d::_fieldIncrSlipE[] = {
    1.100000000000,   2.100000000000,
-   1.563917511528,   1.836082488472,
-   1.300000000000,   1.726466384513,
+   2.005761009464,   1.394238990536,
+   1.300000000000,   1.928356338235,
    1.400000000000,   2.400000000000,
-   1.389619276805,   2.500000000000,
+   1.216364402612,   2.500000000000,
    1.600000000000,   2.600000000000,
-   1.336082488472,   3.063917511528,
-   1.900000000000,   3.473533615487,
-   1.210380723195,   2.100000000000,
+   0.894238990536,   3.505761009464,
+   1.900000000000,   3.271643661765,
+   1.383635597388,   2.100000000000,
    4.520000000000,  -1.920000000000,
    1.000000000000,  -1.600000000000,
   -0.560000000000,   0.200000000000,
 };
 
 const double pylith::faults::CohesiveDynDataTri3d::_slipSlipE[] = {
-  -1.029314160777,   0.000000000000,
-  -1.147067230974,   0.000000000000,
-   0.220761446391,   0.000000000000,
+  -2.279036295232,   0.000000000000,
+  -0.743287323529,   0.000000000000,
+   0.567271194775,   0.000000000000,
 };
 
 // ----------------------------------------------------------------------
@@ -528,23 +528,23 @@
 // Output
 const double pylith::faults::CohesiveDynDataTri3d::_fieldIncrOpenE[] = {
    1.100000000000,   2.100000000000,
-   7.276939399405,   4.204545518954,
-   5.678082695935,   4.684667319277,
+   9.940107961522,   3.662120348918,
+   5.688525058023,   5.314075601410,
    1.400000000000,   2.400000000000,
-   4.440285622520,   4.762452035962,
+   3.622067204316,   6.167076609904,
    1.600000000000,   2.600000000000,
-  -4.376939399405,   0.695454481046,
-  -2.478082695935,   0.515332680723,
-  -1.840285622520,  -0.162452035962,
+  -7.040107961522,   1.237879651082,
+  -2.488525058023,  -0.114075601410,
+  -1.022067204316,  -1.567076609904,
    3.800000000000,  -4.800000000000,
   -3.000000000000,  -4.000000000000,
   -3.200000000000,  -4.200000000000,
 };
 
 const double pylith::faults::CohesiveDynDataTri3d::_slipOpenE[] = {
-  -5.759234657059,  11.428945575657,
-   4.769334638554,   8.756165391870,
-  -5.880571245040,   4.524904071924,
+ -10.292628788528,  14.428129643052,
+   6.028151202820,   8.777050116045,
+  -4.244134408633,   7.334153219809,
 };
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/cohesivedyn.py
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/cohesivedyn.py	2011-10-07 16:25:42 UTC (rev 19041)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/cohesivedyn.py	2011-10-07 17:38:42 UTC (rev 19042)
@@ -58,10 +58,10 @@
                               [8.8, 9.8]])
         fieldIncr = numpy.array([[1.6, 2.6],
                                  [1.8, 2.8]])
-        L = numpy.array([[0.5, 0.0, 0.5, 0.0,],
-                         [0.0, 0.5, 0.0, 0.5,],
-                         [0.5, 0.0, 0.5, 0.0,],
-                         [0.0, 0.5, 0.0, 0.5,],]);
+        L = numpy.array([[1.0, 0.0, 0.0, 0.0,],
+                         [0.0, 1.0, 0.0, 0.0,],
+                         [0.0, 0.0, 1.0, 0.0,],
+                         [0.0, 0.0, 0.0, 1.0,],]);
         C = numpy.array([[0.0, -1.0, 0.0, 0.0,],
                          [-1.0, 0.0, 0.0, 0.0,],
                          [0.0, 0.0, 0.0, -1.0,],
@@ -124,12 +124,12 @@
                                  [1.0, 0.1],
                                  [1.2, 0.2]])
 
-        L = numpy.array([[1.0, 0.0, 0.5, 0.0, 0.5, 0.0],
-                         [0.0, 1.0, 0.0, 0.5, 0.0, 0.5],
-                         [0.5, 0.0, 0.5, 0.0, 0.0, 0.0],
-                         [0.0, 0.5, 0.0, 0.5, 0.0, 0.0],
-                         [0.5, 0.0, 0.0, 0.0, 0.5, 0.0],
-                         [0.0, 0.5, 0.0, 0.0, 0.0, 0.5]])
+        L = numpy.array([[2.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+                         [0.0, 2.0, 0.0, 0.0, 0.0, 0.0],
+                         [0.0, 0.0, 1.0, 0.0, 0.0, 0.0],
+                         [0.0, 0.0, 0.0, 1.0, 0.0, 0.0],
+                         [0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
+                         [0.0, 0.0, 0.0, 0.0, 0.0, 1.0]])
 
         C = numpy.array([[+0.70710678118654757, -0.70710678118654757, 0.0, 0.0, 0.0, 0.0,],
                          [-0.70710678118654757, -0.70710678118654757, 0.0, 0.0, 0.0, 0.0,],
@@ -208,10 +208,10 @@
                               [8.0, 9.0]])
         fieldIncr = numpy.array([[1.8, 2.8],
                                  [1.0, 2.0]])
-        L = numpy.array([[0.5, 0.0, 0.5, 0.0,],
-                         [0.0, 0.5, 0.0, 0.5,],
-                         [0.5, 0.0, 0.5, 0.0,],
-                         [0.0, 0.5, 0.0, 0.5,],]);
+        L = numpy.array([[1.0, 0.0, 0.0, 0.0,],
+                         [0.0, 1.0, 0.0, 0.0,],
+                         [0.0, 0.0, 1.0, 0.0,],
+                         [0.0, 0.0, 0.0, 1.0,],]);
         C = numpy.array([[0.0, -1.0, 0.0, 0.0,],
                          [-1.0, 0.0, 0.0, 0.0,],
                          [0.0, 0.0, 0.0, -1.0,],
@@ -347,15 +347,15 @@
                                  [9.9, 2.9, 3.9],
                                  [9.1, 2.1, 3.1]])
         
-        L = numpy.array([[1.0/9.0,0,0, 1.0/9.0,0,0, 1.0/9.0,0,0,],
-                         [0,1.0/9.0,0, 0,1.0/9.0,0, 0,1.0/9.0,0,],
-                         [0,0,1.0/9.0, 0,0,1.0/9.0, 0,0,1.0/9.0,],
-                         [1.0/9.0,0,0, 1.0/9.0,0,0, 1.0/9.0,0,0,],
-                         [0,1.0/9.0,0, 0,1.0/9.0,0, 0,1.0/9.0,0,],
-                         [0,0,1.0/9.0, 0,0,1.0/9.0, 0,0,1.0/9.0,],
-                         [1.0/9.0,0,0, 1.0/9.0,0,0, 1.0/9.0,0,0,],
-                         [0,1.0/9.0,0, 0,1.0/9.0,0, 0,1.0/9.0,0,],
-                         [0,0,1.0/9.0, 0,0,1.0/9.0, 0,0,1.0/9.0,]])
+        L = numpy.array([[1.0/3.0,0,0, 0.0,0,0, 0.0,0,0,],
+                         [0,1.0/3.0,0, 0,0.0,0, 0,0.0,0,],
+                         [0,0,1.0/3.0, 0,0,0.0, 0,0,0.0,],
+                         [0.0,0,0, 1.0/3.0,0,0, 0.0,0,0,],
+                         [0,0.0,0, 0,1.0/3.0,0, 0,0.0,0,],
+                         [0,0,0.0, 0,0,1.0/3.0, 0,0,0.0,],
+                         [0.0,0,0, 0.0,0,0, 1.0/3.0,0,0,],
+                         [0,0.0,0, 0,0.0,0, 0,1.0/3.0,0,],
+                         [0,0,0.0, 0,0,0.0, 0,0,1.0/3.0,]])
 
         Cv = numpy.array([[ 0, -1, 0,],
                           [ 0, 0, +1,],
@@ -434,9 +434,9 @@
         m = 12
         DOF = 3
 
-        a0 = 4.0/9.0
-        a1 = 2.0/9.0
-        a2 = 1.0/9.0
+        a0 = 1.0
+        a1 = 0.0
+        a2 = 0.0
         L = numpy.array([[a0, 0, 0, a1, 0, 0, a1, 0, 0, a2, 0, 0],
                          [0, a0, 0, 0, a1, 0, 0, a1, 0, 0, a2, 0],
                          [0, 0, a0, 0, 0, a1, 0, 0, a1, 0, 0, a2],



More information about the CIG-COMMITS mailing list