[cig-commits] r13018 - in short/3D/PyLith/trunk: libsrc/faults libsrc/feassemble unittests/libtests/faults

[knepley at geodynamics.org]

Author: knepley
Date: 2008-10-10 18:21:34 -0700 (Fri, 10 Oct 2008)
New Revision: 13018

Modified:
   short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
   short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
Log:
Corrected mistake in FaultCohesiveKin, fixed fault tests


Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc	2008-10-09 20:39:51 UTC (rev 13017)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc	2008-10-11 01:21:34 UTC (rev 13018)
@@ -864,13 +864,19 @@
   const Mesh::label_sequence::iterator  verticesEnd = vertices->end();
   const int                             numVertices = vertices->size();
   Mesh::renumbering_type&               renumbering = _faultMesh->getRenumbering();
+  Mesh::point_type                      firstFaultVertex = -1;
 
   const int fiberDim = solution->getFiberDimension(*vertices->begin());
   double_array tractionValues(fiberDim);
 
   // Allocate buffer for tractions field (if nec.).
+  for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != verticesEnd; ++v_iter) {
+    if (renumbering.find(*v_iter) != renumbering.end()) {
+      firstFaultVertex = renumbering[*v_iter];
+    }
+  }
   if (tractions->isNull() ||
-      fiberDim != (*tractions)->getFiberDimension(renumbering[*vertices->begin()])) {
+      fiberDim != (*tractions)->getFiberDimension(firstFaultVertex)) {
     *tractions = new real_section_type(_faultMesh->comm(), _faultMesh->debug());
     int minE = _faultMesh->getSieve()->getChart().min();
     int maxE = _faultMesh->getSieve()->getChart().max();

Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc	2008-10-09 20:39:51 UTC (rev 13017)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc	2008-10-11 01:21:34 UTC (rev 13018)
@@ -406,6 +406,30 @@
 
     CALL_MEMBER_FN(*this, elasticityJacobianFn)(elasticConsts);
 
+#if 0
+    for(int i = 0; i < numBasis*spaceDim; ++i) {
+      for(int j = 0; j < numBasis*spaceDim; ++j) {
+        if (_cellMatrix[i*numBasis*spaceDim+j] < 1.0e10) {
+          
+        }
+      }
+    }
+    int     n = numBasis*spaceDim, lwork = 5*n, idummy, lierr;
+    double *elemMat = new double[n*n];
+    double *svalues = new double[n];
+    double *work    = new double[lwork];
+    double  sdummy;
+
+    for(int i = 0; i < n; ++i) {elemMat[i] = _cellMatrix[i];}
+    LAPACKgesvd_("N", "N", &n, &n, elemMat, &n, svalues, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, &lierr);
+    if (lierr) {throw std::runtime_error("Lapack SVD failed");}
+    minSV = svalues[n-1];
+    maxSV = svalues[0];
+    delete [] elemMat;
+    delete [] svalues;
+    delete [] work;
+#endif
+
     // Assemble cell contribution into field.  Not sure if this is correct for
     // global stiffness matrix.
     PetscErrorCode err = updateOperator(*mat, *mesh->getSieve(), iV, *c_iter, _cellMatrix, ADD_VALUES);

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc	2008-10-09 20:39:51 UTC (rev 13017)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc	2008-10-11 01:21:34 UTC (rev 13018)
@@ -131,6 +131,7 @@
   FaultCohesiveKin fault;
   _initialize(&mesh, &fault);
   
+  Mesh::renumbering_type& renumbering = fault._faultMesh->getRenumbering();
   const ALE::Obj<Mesh::label_sequence>& vertices = 
     fault._faultMesh->depthStratum(0);
   const Mesh::label_sequence::iterator verticesEnd = vertices->end();
@@ -138,7 +139,8 @@
   for (Mesh::label_sequence::iterator v_iter=vertices->begin();
        v_iter != verticesEnd;
        ++v_iter, ++iVertex) {
-    CPPUNIT_ASSERT_EQUAL(_data->constraintVertices[iVertex],
+    CPPUNIT_ASSERT(renumbering.find(_data->constraintVertices[iVertex]) != renumbering.end());
+    CPPUNIT_ASSERT_EQUAL(renumbering[_data->constraintVertices[iVertex]],
 			 *v_iter);
   } // for
   CPPUNIT_ASSERT_EQUAL(_data->numConstraintVert, iVertex);
@@ -646,6 +648,7 @@
   ALE::Obj<real_section_type> tractions =
     new real_section_type(fault._faultMesh->comm(), fault._faultMesh->debug());
   CPPUNIT_ASSERT(!tractions.isNull());
+  Mesh::renumbering_type& renumbering = fault._faultMesh->getRenumbering();
   const ALE::Obj<Mesh::label_sequence>& vertices = 
     fault._faultMesh->depthStratum(0);
   const Mesh::label_sequence::iterator verticesEnd = vertices->end();
@@ -660,16 +663,27 @@
   for (Mesh::label_sequence::iterator v_iter=vertices->begin();
        v_iter != verticesEnd;
        ++v_iter, ++iVertex) {
+    Mesh::point_type meshVertex = -1;
+    bool found = false;
+
+    for(Mesh::renumbering_type::const_iterator r_iter = renumbering.begin(); r_iter != renumbering.end(); ++r_iter) {
+      if (r_iter->second == *v_iter) {
+        meshVertex = r_iter->first;
+        found      = true;
+        break;
+      }
+    }
+    CPPUNIT_ASSERT(found);
     int fiberDim = tractions->getFiberDimension(*v_iter);
     CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
     const real_section_type::value_type* vertexTractions = 
       tractions->restrictPoint(*v_iter);
     CPPUNIT_ASSERT(0 != vertexTractions);
 
-    fiberDim = solution->getFiberDimension(*v_iter);
+    fiberDim = solution->getFiberDimension(meshVertex);
     CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
     const real_section_type::value_type* vertexSolution = 
-      solution->restrictPoint(*v_iter);
+      solution->restrictPoint(meshVertex);
     CPPUNIT_ASSERT(0 != vertexSolution);
 
     const double scale = _data->pseudoStiffness / _data->area[iVertex];