[cig-commits] r14941 - in short/3D/PyLith/trunk: libsrc/bc unittests/libtests/bc unittests/libtests/bc/data

brad at geodynamics.org brad at geodynamics.org
Sat May 9 11:17:18 PDT 2009


Author: brad
Date: 2009-05-09 11:17:16 -0700 (Sat, 09 May 2009)
New Revision: 14941

Modified:
   short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.hh
   short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.hh
   short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.hh
   short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.hh
   short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.hh
   short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc
   short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.hh
   short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py
Log:
Updated absorbing dampers for increment solution.

Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc	2009-05-09 18:17:16 UTC (rev 14941)
@@ -283,12 +283,12 @@
   topology::SubMesh::UpdateAddVisitor residualVisitor(*residualSection,
 						      &_cellVector[0]);
 
-  const topology::Field<topology::Mesh>& dispTpdt = 
-    fields->get("disp(t+dt)");
-  const ALE::Obj<RealSection>& dispTpdtSection = dispTpdt.section();
-  assert(!dispTpdtSection.isNull());
-  topology::Mesh::RestrictVisitor dispTpdtVisitor(*dispTpdtSection, 
-						  numBasis*spaceDim);
+  const topology::Field<topology::Mesh>& dispT = 
+    fields->get("disp(t)");
+  const ALE::Obj<RealSection>& dispTSection = dispT.section();
+  assert(!dispTSection.isNull());
+  topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection, 
+					       numBasis*spaceDim);
   const topology::Field<topology::Mesh>& dispTmdt = 
     fields->get("disp(t-dt)");
   const ALE::Obj<RealSection>& dispTmdtSection = dispTmdt.section();
@@ -310,16 +310,15 @@
     _resetCellVector();
 
     // Restrict input fields to cell
-    dispTpdtVisitor.clear();
-    submesh->restrictClosure(*c_iter, dispTpdtVisitor);
-    const double* dispTpdtCell = dispTpdtVisitor.getValues();
+    dispTVisitor.clear();
+    submesh->restrictClosure(*c_iter, dispTVisitor);
+    const double* dispTCell = dispTVisitor.getValues();
 
     dispTmdtVisitor.clear();
     submesh->restrictClosure(*c_iter, dispTmdtVisitor);
     const double* dispTmdtCell = dispTmdtVisitor.getValues();
 
-    dampersSection->restrictPoint(*c_iter, 
-				  &dampersCell[0], dampersCell.size());
+    dampersSection->restrictPoint(*c_iter, &dampersCell[0], dampersCell.size());
 
     // Get cell geometry information that depends on cell
     const double_array& basis = _quadrature->basis();
@@ -337,7 +336,8 @@
           for (int iDim=0; iDim < spaceDim; ++iDim)
             _cellVector[iBasis*spaceDim+iDim] += 
 	      dampersCell[iQuad*spaceDim+iDim] *
-	      valIJ * (dispTmdtCell[jBasis*spaceDim+iDim]);
+	      valIJ * (dispTmdtCell[jBasis*spaceDim+iDim] - \
+		       dispTCell[jBasis*spaceDim+iDim]);
         } // for
       } // for
 

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc	2009-05-09 18:17:16 UTC (rev 14941)
@@ -312,13 +312,10 @@
     // Setup fields
     CPPUNIT_ASSERT(0 != fields);
     fields->add("residual", "residual");
-    fields->add("disp(t+dt)", "displacement");
+    fields->add("dispIncr(t->t+dt)", "displacement_increment");
     fields->add("disp(t)", "displacement");
     fields->add("disp(t-dt)", "displacement");
-    fields->solutionName("disp(t+dt)");
-    const char* history[] = { "disp(t+dt)", "disp(t)", "disp(t-dt)" };
-    const int historySize = 3;
-    fields->createHistory(history, historySize);
+    fields->solutionName("dispIncr(t->t+dt)");
   
     topology::Field<topology::Mesh>& residual = fields->get("residual");
     const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
@@ -333,20 +330,20 @@
 
     const int totalNumVertices = sieveMesh->depthStratum(0)->size();
     const int numMeshCells = sieveMesh->heightStratum(0)->size();
-    const int fieldSize = _data->spaceDim * totalNumVertices;
-    const ALE::Obj<RealSection>& dispTpdtSection = 
-      fields->get("disp(t+dt)").section();
+ int fieldSize = _data->spaceDim * totalNumVertices;
+    const ALE::Obj<RealSection>& dispTIncrSection = 
+      fields->get("dispIncr(t->t+dt)").section();
     const ALE::Obj<RealSection>& dispTSection = 
       fields->get("disp(t)").section();
     const ALE::Obj<RealSection>& dispTmdtSection = 
       fields->get("disp(t-dt)").section();
-    CPPUNIT_ASSERT(!dispTpdtSection.isNull());
+    CPPUNIT_ASSERT(!dispTIncrSection.isNull());
     CPPUNIT_ASSERT(!dispTSection.isNull());
     CPPUNIT_ASSERT(!dispTmdtSection.isNull());
     const int offset = numMeshCells;
     for (int iVertex=0; iVertex < totalNumVertices; ++iVertex) {
-      dispTpdtSection->updatePoint(iVertex+offset, 
-				   &_data->fieldTpdt[iVertex*_data->spaceDim]);
+      dispTIncrSection->updatePoint(iVertex+offset, 
+				   &_data->fieldTIncr[iVertex*_data->spaceDim]);
       dispTSection->updatePoint(iVertex+offset, 
 				&_data->fieldT[iVertex*_data->spaceDim]);
       dispTmdtSection->updatePoint(iVertex+offset, 

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.cc	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.cc	2009-05-09 18:17:16 UTC (rev 14941)
@@ -26,7 +26,7 @@
   id(0),
   label(""),
   dt(0),
-  fieldTpdt(0),
+  fieldTIncr(0),
   fieldT(0),
   fieldTmdt(0),
   spaceDim(0),

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.hh	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.hh	2009-05-09 18:17:16 UTC (rev 14941)
@@ -56,7 +56,7 @@
   /// @name Input fields
   //@{
   double dt; ///< Time step
-  double* fieldTpdt; ///< Input field at time t+dt.
+  double* fieldTIncr; ///< Input increment field for time to to t+dt.
   double* fieldT; ///< Input field at time t.
   double* fieldTmdt; ///< Input field at time t-dt.
   //@}

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.cc	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.cc	2009-05-09 18:17:16 UTC (rev 14941)
@@ -85,7 +85,7 @@
   3.1,  0.3,  7.2,
   3.3,  0.1,  7.6,
 };
-const double pylith::bc::AbsorbingDampersDataHex8::_fieldTpdt[] = {
+const double pylith::bc::AbsorbingDampersDataHex8::_fieldTIncr[] = {
   1.2,  1.1,  3.4,
   1.5,  1.0,  4.0,
   1.8,  0.9,  4.6,
@@ -122,18 +122,18 @@
   7.5e+06,  1.25e+07,  7.5e+06,
 };
 const double pylith::bc::AbsorbingDampersDataHex8::_valsResidual[] = {
-  9.50000021e+06,  2.33333337e+07,  2.65000006e+07,
-  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-  2.10000004e+07,  4.00000007e+07,  5.70000012e+07,
-  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-  1.15000002e+07,  1.66666670e+07,  3.05000006e+07,
-  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-  1.10000002e+07,  1.83333337e+07,  2.95000006e+07,
-  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-  2.40000005e+07,  3.00000007e+07,  6.30000012e+07,
-  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
-  1.30000002e+07,  1.16666670e+07,  3.35000006e+07,
-  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
+   -2.75000007e+06,    1.25000003e+06,   -5.50000014e+06,
+    0.00000000e+00,    0.00000000e+00,    0.00000000e+00,
+   -7.50000018e+06,    2.50000005e+06,   -1.50000004e+07,
+    0.00000000e+00,    0.00000000e+00,    0.00000000e+00,
+   -4.75000010e+06,    1.25000003e+06,   -9.50000021e+06,
+    0.00000000e+00,    0.00000000e+00,    0.00000000e+00,
+   -4.25000008e+06,    1.25000002e+06,   -8.50000015e+06,
+    0.00000000e+00,    0.00000000e+00,    0.00000000e+00,
+   -1.05000002e+07,    2.50000005e+06,   -2.10000004e+07,
+    0.00000000e+00,    0.00000000e+00,    0.00000000e+00,
+   -6.25000011e+06,    1.25000003e+06,   -1.25000002e+07,
+    0.00000000e+00,    0.00000000e+00,    0.00000000e+00,
 };
 const double pylith::bc::AbsorbingDampersDataHex8::_valsJacobian[] = {
   3.33333333e+06, 0.0, 0.0, // 0x
@@ -586,7 +586,7 @@
   label = const_cast<char*>(_label);
 
   dt = _dt;
-  fieldTpdt = const_cast<double*>(_fieldTpdt);
+  fieldTIncr = const_cast<double*>(_fieldTIncr);
   fieldT = const_cast<double*>(_fieldT);
   fieldTmdt = const_cast<double*>(_fieldTmdt);
 

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.hh	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.hh	2009-05-09 18:17:16 UTC (rev 14941)
@@ -48,7 +48,7 @@
   static const char* _label;
 
   static const double _dt;
-  static const double _fieldTpdt[];
+  static const double _fieldTIncr[];
   static const double _fieldT[];
   static const double _fieldTmdt[];
 

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.cc	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.cc	2009-05-09 18:17:16 UTC (rev 14941)
@@ -45,7 +45,7 @@
   1.3,
   1.5,
 };
-const double pylith::bc::AbsorbingDampersDataLine2::_fieldTpdt[] = {
+const double pylith::bc::AbsorbingDampersDataLine2::_fieldTIncr[] = {
   1.2,
   1.5,
   1.8,
@@ -67,9 +67,9 @@
   17.5e+6,
 };
 const double pylith::bc::AbsorbingDampersDataLine2::_valsResidual[] = {
-  12.5e+6*1.0/0.5,
+  12.5e+6*(1.0-1.1)/0.5,
   0.0,
-  17.5e+6*1.2/0.5,
+  17.5e+6*(1.2-1.5)/0.5,
 };
 const double pylith::bc::AbsorbingDampersDataLine2::_valsJacobian[] = {
   12.5e+6/0.5, 0.0, 0.0,
@@ -93,7 +93,7 @@
   label = const_cast<char*>(_label);
 
   dt = _dt;
-  fieldTpdt = const_cast<double*>(_fieldTpdt);
+  fieldTIncr = const_cast<double*>(_fieldTIncr);
   fieldT = const_cast<double*>(_fieldT);
   fieldTmdt = const_cast<double*>(_fieldTmdt);
 

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.hh	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.hh	2009-05-09 18:17:16 UTC (rev 14941)
@@ -48,7 +48,7 @@
   static const char* _label;
 
   static const double _dt;
-  static const double _fieldTpdt[];
+  static const double _fieldTIncr[];
   static const double _fieldT[];
   static const double _fieldTmdt[];
 

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.cc	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.cc	2009-05-09 18:17:16 UTC (rev 14941)
@@ -53,7 +53,7 @@
   1.9,  2.4,
   2.1,  2.5,
 };
-const double pylith::bc::AbsorbingDampersDataQuad4::_fieldTpdt[] = {
+const double pylith::bc::AbsorbingDampersDataQuad4::_fieldTIncr[] = {
   1.2,  1.6,
   1.5,  2.4,
   1.8,  2.0,
@@ -78,12 +78,12 @@
   1.25e+07, 7.5e+06,
 };
 const double pylith::bc::AbsorbingDampersDataQuad4::_valsResidual[] = {
-  2.625e+07, 3.15e+07,
-  2.625e+07, 3.15e+07,
+  -3.75000000e+06,    7.50000000e+05,
+  -3.75000000e+06,    7.50000000e+05,
   0.0, 0.0, 
   0.0, 0.0,
-  3.625e+07, 3.0e+07,
-  3.625e+07, 3.0e+07,
+  -1.37500000e+07,   -6.75000000e+06,
+  -1.37500000e+07,   -6.75000000e+06,
 };
 const double pylith::bc::AbsorbingDampersDataQuad4::_valsJacobian[] = {
   1.25e+07, 0.0, // 0x
@@ -176,7 +176,7 @@
   label = const_cast<char*>(_label);
 
   dt = _dt;
-  fieldTpdt = const_cast<double*>(_fieldTpdt);
+  fieldTIncr = const_cast<double*>(_fieldTIncr);
   fieldT = const_cast<double*>(_fieldT);
   fieldTmdt = const_cast<double*>(_fieldTmdt);
 

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.hh	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.hh	2009-05-09 18:17:16 UTC (rev 14941)
@@ -48,7 +48,7 @@
   static const char* _label;
 
   static const double _dt;
-  static const double _fieldTpdt[];
+  static const double _fieldTIncr[];
   static const double _fieldT[];
   static const double _fieldTmdt[];
 

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.cc	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.cc	2009-05-09 18:17:16 UTC (rev 14941)
@@ -52,7 +52,7 @@
   1.5,  2.2,  4.0,
   1.7,  2.3,  4.4,
 };
-const double pylith::bc::AbsorbingDampersDataTet4::_fieldTpdt[] = {
+const double pylith::bc::AbsorbingDampersDataTet4::_fieldTIncr[] = {
   1.2,  1.6,  3.4,
   1.5,  2.4,  4.0,
   1.8,  2.0,  4.6,
@@ -73,10 +73,10 @@
   1.25e+07,  7.5e+06,  7.5e+06
 };
 const double pylith::bc::AbsorbingDampersDataTet4::_valsResidual[] = {
-  4.861111111e+06,  5.83333333e+06,  8.33333333e+06,
+  -1.11111111e+06,    4.16666667e+05,   -1.33333333e+06,
   0.0,              0.0,             0.0,
-  4.861111111e+06,  5.83333333e+06,  8.33333333e+06,
-  4.861111111e+06,  5.83333333e+06,  8.33333333e+06,
+  -1.11111111e+06,    4.16666667e+05,   -1.33333333e+06,
+  -1.11111111e+06,    4.16666667e+05,   -1.33333333e+06,
   0.0,              0.0,             0.0,  
 };
 const double pylith::bc::AbsorbingDampersDataTet4::_valsJacobian[] = {
@@ -173,7 +173,7 @@
   label = const_cast<char*>(_label);
 
   dt = _dt;
-  fieldTpdt = const_cast<double*>(_fieldTpdt);
+  fieldTIncr = const_cast<double*>(_fieldTIncr);
   fieldT = const_cast<double*>(_fieldT);
   fieldTmdt = const_cast<double*>(_fieldTmdt);
 

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.hh	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.hh	2009-05-09 18:17:16 UTC (rev 14941)
@@ -48,7 +48,7 @@
   static const char* _label;
 
   static const double _dt;
-  static const double _fieldTpdt[];
+  static const double _fieldTIncr[];
   static const double _fieldT[];
   static const double _fieldTmdt[];
 

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc	2009-05-09 18:17:16 UTC (rev 14941)
@@ -50,7 +50,7 @@
   1.5,  2.2,
   1.7,  2.3,
 };
-const double pylith::bc::AbsorbingDampersDataTri3::_fieldTpdt[] = {
+const double pylith::bc::AbsorbingDampersDataTri3::_fieldTIncr[] = {
   1.2,  1.6,
   1.5,  2.4,
   1.8,  2.0,
@@ -72,9 +72,9 @@
 };
 const double pylith::bc::AbsorbingDampersDataTri3::_valsResidual[] = {
   0.0, 0.0,
-  2.4e+07,  1.0e+07,
+  -6.00000000e+06,   -1.00000000e+06,
   0.0, 0.0,
-  2.4e+07,  1.0e+07,
+  -6.00000000e+06,   -1.00000000e+06,
 };
 const double pylith::bc::AbsorbingDampersDataTri3::_valsJacobian[] = {
   0.0, 0.0, // 0x
@@ -127,7 +127,7 @@
   label = const_cast<char*>(_label);
 
   dt = _dt;
-  fieldTpdt = const_cast<double*>(_fieldTpdt);
+  fieldTIncr = const_cast<double*>(_fieldTIncr);
   fieldT = const_cast<double*>(_fieldT);
   fieldTmdt = const_cast<double*>(_fieldTmdt);
 

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.hh	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.hh	2009-05-09 18:17:16 UTC (rev 14941)
@@ -48,7 +48,7 @@
   static const char* _label;
 
   static const double _dt;
-  static const double _fieldTpdt[];
+  static const double _fieldTIncr[];
   static const double _fieldT[];
   static const double _fieldTmdt[];
 

Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py	2009-05-09 02:37:34 UTC (rev 14940)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py	2009-05-09 18:17:16 UTC (rev 14941)
@@ -31,8 +31,8 @@
   edgeLen = 2.0**0.5
   N0 = 0.5
   N1 = 0.5
-  dispX = 0.5*(N0+N1)*(N0*(1.1)+N1*(1.3)) / (2.0*dt)
-  dispY = 0.5*(N0+N1)*(N0*(1.8)+N1*(2.2)) / (2.0*dt)
+  dispX = 0.5*(N0+N1)*(N0*(1.1-1.3)+N1*(1.3-1.7)) / (2.0*dt)
+  dispY = 0.5*(N0+N1)*(N0*(1.8-2.1)+N1*(2.2-2.3)) / (2.0*dt)
   normal = [0.5**0.5, -0.5**0.5]
 
   constNormal = density*vp
@@ -78,12 +78,12 @@
   dt = 0.25
   N0 = 0.5
   N1 = 0.5
-  disp1X = 0.5*(N0+N1)*(N0*(1.1) + N1*(1.0)) / (2.0*dt)
-  disp1Y = 0.5*(N0+N1)*(N0*(1.8) + N1*(2.4)) / (2.0*dt)
+  disp1X = 0.5*(N0+N1)*(N0*(1.1-1.3) + N1*(1.0-1.1)) / (2.0*dt)
+  disp1Y = 0.5*(N0+N1)*(N0*(1.8-2.1) + N1*(2.4-2.0)) / (2.0*dt)
   normal1 = [-1.0, 0.0]
 
-  disp2X = 0.5*(N0+N1)*(N0*(1.4) + N1*(1.5)) / (2.0*dt)
-  disp2Y = 0.5*(N0+N1)*(N0*(2.4) + N1*(1.6)) / (2.0*dt)
+  disp2X = 0.5*(N0+N1)*(N0*(1.4-1.9) + N1*(1.5-2.1)) / (2.0*dt)
+  disp2Y = 0.5*(N0+N1)*(N0*(2.4-2.4) + N1*(1.6-2.5)) / (2.0*dt)
   normal2 = [1.0, 0.0]
 
   edgeLen = 2.0
@@ -145,9 +145,9 @@
   N0 = 1.0/3.0
   N1 = 1.0/3.0
   N2 = 1.0/3.0
-  dispX = (N0+N1+N2)/3.0 * (N0*(1.3)+N1*(1.2)+N2*(1.0)) / (2.0*dt)
-  dispY = (N0+N1+N2)/3.0 * (N0*(2.2)+N1*(2.4)+N2*(2.4)) / (2.0*dt)
-  dispZ = (N0+N1+N2)/3.0 * (N0*(3.6)+N1*(3.4)+N2*(3.0)) / (2.0*dt)
+  dispX = (N0+N1+N2)/3.0 * (N0*(1.3-1.7)+N1*(1.2-1.5)+N2*(1.0-1.1)) / (2.0*dt)
+  dispY = (N0+N1+N2)/3.0 * (N0*(2.2-2.3)+N1*(2.4-2.2)+N2*(2.4-2.0)) / (2.0*dt)
+  dispZ = (N0+N1+N2)/3.0 * (N0*(3.6-4.4)+N1*(3.4-4.0)+N2*(3.0-3.2)) / (2.0*dt)
   normal = [-1.0, 0.0, 0.0]
   tangent1 = [0.0, -1.0, 0.0]
   tangent2 = [0.0, 0.0, 1.0]
@@ -233,6 +233,18 @@
               [1.9,  0.6,  4.8],
               [2.0,  0.4,  5.0],
               [2.1,  0.2,  5.2]]
+  dispT = [[1.1,  2.3,  3.2],
+           [1.3,  2.1,  3.6],
+           [1.5,  1.9,  4.0],
+           [1.7,  1.7,  4.4],
+           [1.9,  1.5,  4.8],
+           [2.1,  1.3,  5.2],
+           [2.3,  1.1,  5.6],
+           [2.5,  0.9,  6.0],
+           [2.7,  0.7,  6.4],
+           [2.9,  0.5,  6.8],
+           [3.1,  0.3,  7.2],
+           [3.3,  0.1,  7.6]]
   normal = [0.0, 1.0, 0.0]
   tangent1 = [-1.0, 0.0, 0.0]
   tangent2 = [0.0, 0.0, 1.0]
@@ -257,18 +269,18 @@
       N1 = b[1]
       N2 = b[2]
       N3 = b[3]
-      d0x = dispTmdt[cell[0]][0]/(2.0*dt)
-      d0y = dispTmdt[cell[0]][1]/(2.0*dt)
-      d0z = dispTmdt[cell[0]][2]/(2.0*dt)
-      d1x = dispTmdt[cell[1]][0]/(2.0*dt)
-      d1y = dispTmdt[cell[1]][1]/(2.0*dt)
-      d1z = dispTmdt[cell[1]][2]/(2.0*dt)
-      d2x = dispTmdt[cell[2]][0]/(2.0*dt)
-      d2y = dispTmdt[cell[2]][1]/(2.0*dt)
-      d2z = dispTmdt[cell[2]][2]/(2.0*dt)
-      d3x = dispTmdt[cell[3]][0]/(2.0*dt)
-      d3y = dispTmdt[cell[3]][1]/(2.0*dt)
-      d3z = dispTmdt[cell[3]][2]/(2.0*dt)
+      d0x = (dispTmdt[cell[0]][0] - dispT[cell[0]][0])/(2.0*dt)
+      d0y = (dispTmdt[cell[0]][1] - dispT[cell[0]][1])/(2.0*dt)
+      d0z = (dispTmdt[cell[0]][2] - dispT[cell[0]][2])/(2.0*dt)
+      d1x = (dispTmdt[cell[1]][0] - dispT[cell[1]][0])/(2.0*dt)
+      d1y = (dispTmdt[cell[1]][1] - dispT[cell[1]][1])/(2.0*dt)
+      d1z = (dispTmdt[cell[1]][2] - dispT[cell[1]][2])/(2.0*dt)
+      d2x = (dispTmdt[cell[2]][0] - dispT[cell[2]][0])/(2.0*dt)
+      d2y = (dispTmdt[cell[2]][1] - dispT[cell[2]][1])/(2.0*dt)
+      d2z = (dispTmdt[cell[2]][2] - dispT[cell[2]][2])/(2.0*dt)
+      d3x = (dispTmdt[cell[3]][0] - dispT[cell[3]][0])/(2.0*dt)
+      d3y = (dispTmdt[cell[3]][1] - dispT[cell[3]][1])/(2.0*dt)
+      d3z = (dispTmdt[cell[3]][2] - dispT[cell[3]][2])/(2.0*dt)
       dispX = N0*d0x + N1*d1x + N2*d2x + N3*d3x
       dispY = N0*d0y + N1*d1y + N2*d2y + N3*d3y
       dispZ = N0*d0z + N1*d1z + N2*d2z + N3*d3z
@@ -334,7 +346,7 @@
       print "  %16.8e" % v
   print "values for residual:"
   for v in numpy.ravel(residual):
-      print "  %16.8e" % v
+      print "  %16.8e," % v
   print "values for jacobian:"
   for j in numpy.ravel(jacobian):
       print "  %16.8e" % j
@@ -343,8 +355,8 @@
 # ----------------------------------------------------------------------
 #calcTri3()
 #calcQuad4()
-#calcTet4()
-calcHex8()
+calcTet4()
+#calcHex8()
 
   
 # End of file 



More information about the CIG-COMMITS mailing list