[cig-commits] r8177 - in short/3D/PyLith/trunk: libsrc/bc
unittests/libtests/bc unittests/libtests/bc/data
brad at geodynamics.org
brad at geodynamics.org
Thu Oct 25 14:56:25 PDT 2007
Author: brad
Date: 2007-10-25 14:56:24 -0700 (Thu, 25 Oct 2007)
New Revision: 8177
Added:
short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampersQuad4.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampersQuad4.hh
short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.hh
Modified:
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py
short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4.mesh
Log:
More work on absorbing boundary unit tests. Implemented tests for mesh with quad4 cells.
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2007-10-23 00:22:46 UTC (rev 8176)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2007-10-25 21:56:24 UTC (rev 8177)
@@ -173,9 +173,12 @@
<< "using spatial database " << _db->label() << ".";
throw std::runtime_error(msg.str());
} // if
- dampingConstsLocal[index ] = queryData[0]*queryData[1];
- for (int iDim=1; iDim < spaceDim; ++iDim)
- dampingConstsLocal[index+iDim] = queryData[0]*queryData[2];
+ const double constTangential = queryData[0]*queryData[2];
+ const double constNormal = queryData[0]*queryData[1];
+ const int numTangential = spaceDim-1;
+ for (int iDim=0; iDim < numTangential; ++iDim)
+ dampingConstsLocal[iDim] = constTangential;
+ dampingConstsLocal[spaceDim-1] = constNormal;
// Compute normal/tangential orientation
_boundaryMesh->restrict(coordinates, *c_iter,
@@ -185,6 +188,7 @@
cellGeometry.orientation(&orientation, jacobian, jacobianDet,
upDir);
orientation /= jacobianDet;
+
dampingConstsGlobal = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim)
@@ -274,6 +278,7 @@
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt =
quadWts[iQuad] * jacobianDet[iQuad] / (2.0 * dt);
+
for (int iBasis=0; iBasis < numBasis; ++iBasis) {
const double valI = wt*basis[iQuad*numBasis+iBasis];
for (int jBasis=0; jBasis < numBasis; ++jBasis) {
@@ -281,10 +286,11 @@
for (int iDim=0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] +=
dampingConstsCell[iQuad*spaceDim+iDim] *
- valIJ * (- dispTpdtCell[jBasis*spaceDim+iDim]
+ valIJ * (-dispTpdtCell[jBasis*spaceDim+iDim]
+ dispTmdtCell[jBasis*spaceDim+iDim]);
} // for
} // for
+
} // for
PetscLogFlopsNoCheck(numQuadPts*(3+numBasis*(1+numBasis*(5*spaceDim))));
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am 2007-10-23 00:22:46 UTC (rev 8176)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am 2007-10-25 21:56:24 UTC (rev 8177)
@@ -24,6 +24,7 @@
TestAbsorbingDampers.cc \
TestAbsorbingDampersLine2.cc \
TestAbsorbingDampersTri3.cc \
+ TestAbsorbingDampersQuad4.cc \
TestBoundaryCondition.cc \
TestDirichlet.cc \
TestDirichletLine2.cc \
@@ -42,6 +43,7 @@
TestAbsorbingDampers.hh \
TestAbsorbingDampersLine2.hh \
TestAbsorbingDampersTri3.hh \
+ TestAbsorbingDampersQuad4.hh \
TestBoundaryCondition.hh \
TestDirichlet.hh \
TestDirichletLine2.hh \
@@ -60,6 +62,7 @@
data/AbsorbingDampersData.cc \
data/AbsorbingDampersDataLine2.cc \
data/AbsorbingDampersDataTri3.cc \
+ data/AbsorbingDampersDataQuad4.cc \
data/DirichletData.cc \
data/DirichletDataLine2.cc \
data/DirichletDataLine2b.cc \
@@ -76,6 +79,7 @@
data/AbsorbingDampersData.hh \
data/AbsorbingDampersDataLine2.hh \
data/AbsorbingDampersDataTri3.hh \
+ data/AbsorbingDampersDataQuad4.hh \
data/DirichletData.hh \
data/DirichletDataLine2.hh \
data/DirichletDataLine2b.hh \
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2007-10-23 00:22:46 UTC (rev 8176)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2007-10-25 21:56:24 UTC (rev 8177)
@@ -143,6 +143,8 @@
const int size = residual->sizeWithBC();
CPPUNIT_ASSERT_EQUAL(sizeE, size);
+ //residual->view("RESIDUAL");
+
const double tolerance = 1.0e-06;
for (int i=0; i < size; ++i)
if (fabs(valsE[i]) > 1.0)
@@ -205,6 +207,14 @@
for (int iCol=0; iCol < ncols; ++iCol)
cols[iCol] = iCol;
MatGetValues(jDense, nrows, &rows[0], ncols, &cols[0], &vals[0]);
+
+#if 0
+ std::cout << "JACOBIAN\n";
+ for (int iRow=0, i=0; iRow < nrows; ++iRow)
+ for (int iCol=0; iCol < ncols; ++iCol, ++i)
+ std::cout << " iRow: " << iRow << ", iCol: " << iCol << ", value: " << vals[i] << std::endl;
+#endif
+
const double tolerance = 1.0e-06;
for (int iRow=0; iRow < nrows; ++iRow)
for (int iCol=0; iCol < ncols; ++iCol) {
@@ -283,6 +293,7 @@
fields->copyLayout("residual");
const int totalNumVertices = (*mesh)->depthStratum(0)->size();
+ const int numMeshCells = (*mesh)->heightStratum(0)->size();
const int fieldSize = _data->spaceDim * totalNumVertices;
const ALE::Obj<real_section_type>& dispTpdt = fields->getReal("dispTpdt");
const ALE::Obj<real_section_type>& dispT = fields->getReal("dispT");
@@ -290,7 +301,7 @@
CPPUNIT_ASSERT(!dispTpdt.isNull());
CPPUNIT_ASSERT(!dispT.isNull());
CPPUNIT_ASSERT(!dispTmdt.isNull());
- const int offset = _data->numCells;
+ const int offset = numMeshCells;
for (int iVertex=0; iVertex < totalNumVertices; ++iVertex) {
dispTpdt->updatePoint(iVertex+offset,
&_data->fieldTpdt[iVertex*_data->spaceDim]);
Added: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampersQuad4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampersQuad4.cc 2007-10-23 00:22:46 UTC (rev 8176)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampersQuad4.cc 2007-10-25 21:56:24 UTC (rev 8177)
@@ -0,0 +1,38 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestAbsorbingDampersQuad4.hh" // Implementation of class methods
+
+#include "data/AbsorbingDampersDataQuad4.hh" // USES AbsorbingDampersDataQuad4
+
+#include "pylith/feassemble/Quadrature1Din2D.hh" // USES Quadrature1Din2D
+#include "pylith/feassemble/GeometryLine2D.hh" // USES GeometryLine2D
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::bc::TestAbsorbingDampersQuad4 );
+
+// ----------------------------------------------------------------------
+// Setup testing data.
+void
+pylith::bc::TestAbsorbingDampersQuad4::setUp(void)
+{ // setUp
+ _data = new AbsorbingDampersDataQuad4();
+ _quadrature = new feassemble::Quadrature1Din2D();
+ CPPUNIT_ASSERT(0 != _quadrature);
+ feassemble::GeometryLine2D geometry;
+ _quadrature->refGeometry(&geometry);
+} // setUp
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampersQuad4.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampersQuad4.hh 2007-10-23 00:22:46 UTC (rev 8176)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampersQuad4.hh 2007-10-25 21:56:24 UTC (rev 8177)
@@ -0,0 +1,55 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/bc/TestAbsorbingDampersQuad4.hh
+ *
+ * @brief C++ TestAbsorbingDampers object.
+ *
+ * C++ unit testing for AbsorbingDampers for mesh with 2-D quad cells.
+ */
+
+#if !defined(pylith_bc_testabsorbingdampersquad4_hh)
+#define pylith_bc_testabsorbingdampersquad4_hh
+
+#include "TestAbsorbingDampers.hh" // ISA TestAbsorbingDampers
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace bc {
+ class TestAbsorbingDampersQuad4;
+ } // bc
+} // pylith
+
+/// C++ unit testing for AbsorbingDampers for mesh with 2-D quad cells.
+class pylith::bc::TestAbsorbingDampersQuad4 : public TestAbsorbingDampers
+{ // class TestAbsorbingDampers
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestAbsorbingDampersQuad4 );
+ CPPUNIT_TEST( testInitialize );
+ CPPUNIT_TEST( testIntegrateResidual );
+ CPPUNIT_TEST( testIntegrateJacobian );
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Setup testing data.
+ void setUp(void);
+
+}; // class TestAbsorbingDampersQuad4
+
+#endif // pylith_bc_absorbingdampersquad4_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.cc 2007-10-23 00:22:46 UTC (rev 8176)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.cc 2007-10-25 21:56:24 UTC (rev 8177)
@@ -0,0 +1,199 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include "AbsorbingDampersDataQuad4.hh"
+
+const char* pylith::bc::AbsorbingDampersDataQuad4::_meshFilename =
+ "data/quad4.mesh";
+
+const int pylith::bc::AbsorbingDampersDataQuad4::_numBasis = 2;
+const int pylith::bc::AbsorbingDampersDataQuad4::_numQuadPts = 1;
+const double pylith::bc::AbsorbingDampersDataQuad4::_quadPts[] = {
+ 0.0,
+};
+const double pylith::bc::AbsorbingDampersDataQuad4::_quadWts[] = {
+ 2.0,
+};
+const double pylith::bc::AbsorbingDampersDataQuad4::_basis[] = {
+ 0.5,
+ 0.5,
+};
+const double pylith::bc::AbsorbingDampersDataQuad4::_basisDerivRef[] = {
+ -0.5,
+ 0.5,
+};
+
+const char* pylith::bc::AbsorbingDampersDataQuad4::_spatialDBFilename = "data/elasticplanestrain.spatialdb";
+const int pylith::bc::AbsorbingDampersDataQuad4::_id = 2;
+const char* pylith::bc::AbsorbingDampersDataQuad4::_label = "bc2";
+
+const double pylith::bc::AbsorbingDampersDataQuad4::_dt = 0.25;
+const double pylith::bc::AbsorbingDampersDataQuad4::_fieldTmdt[] = {
+ 1.0, 2.4,
+ 1.1, 1.8,
+ 1.2, 2.4,
+ 1.3, 2.2,
+ 1.4, 2.4,
+ 1.5, 1.6,
+};
+const double pylith::bc::AbsorbingDampersDataQuad4::_fieldT[] = {
+ 1.1, 2.0,
+ 1.3, 2.1,
+ 1.5, 2.2,
+ 1.7, 2.3,
+ 1.9, 2.4,
+ 2.1, 2.5,
+};
+const double pylith::bc::AbsorbingDampersDataQuad4::_fieldTpdt[] = {
+ 1.2, 1.6,
+ 1.5, 2.4,
+ 1.8, 2.0,
+ 2.1, 2.4,
+ 2.4, 2.4,
+ 2.7, 3.4,
+};
+
+const int pylith::bc::AbsorbingDampersDataQuad4::_spaceDim = 2;
+const int pylith::bc::AbsorbingDampersDataQuad4::_cellDim = 1;
+const int pylith::bc::AbsorbingDampersDataQuad4::_numVertices = 4;
+const int pylith::bc::AbsorbingDampersDataQuad4::_numCells = 2;
+const int pylith::bc::AbsorbingDampersDataQuad4::_numCorners = 2;
+const int pylith::bc::AbsorbingDampersDataQuad4::_cells[] = {
+ 3, 2,
+ 6, 7,
+};
+
+
+const double pylith::bc::AbsorbingDampersDataQuad4::_dampingConsts[] = {
+ 1.25e+07, 7.5e+06,
+ 1.25e+07, 7.5e+06,
+};
+const double pylith::bc::AbsorbingDampersDataQuad4::_valsResidual[] = {
+ -7.5e+06, 1.5e+06,
+ -7.5e+06, 1.5e+06,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ -2.75e+07, -1.35e+07,
+ -2.75e+07, -1.35e+07,
+};
+const double pylith::bc::AbsorbingDampersDataQuad4::_valsJacobian[] = {
+ 1.25e+07, 0.0, // 0x
+ 1.25e+07, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 7.5e+06, // 0y
+ 0.0, 7.5e+06,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 1.25e+07, 0.0, // 1x
+ 1.25e+07, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 7.5e+06, // 1y
+ 0.0, 7.5e+06,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0, // 2x
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0, // 2y
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0, // 3x
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0, // 3y
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0, // 4x
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 1.25e+07, 0.0,
+ 1.25e+07, 0.0,
+ 0.0, 0.0, // 4y
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 7.5e+06,
+ 0.0, 7.5e+06,
+ 0.0, 0.0, // 5x
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 1.25e+07, 0.0,
+ 1.25e+07, 0.0,
+ 0.0, 0.0, // 5y
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 7.5e+06,
+ 0.0, 7.5e+06,
+};
+
+pylith::bc::AbsorbingDampersDataQuad4::AbsorbingDampersDataQuad4(void)
+{ // constructor
+ meshFilename = const_cast<char*>(_meshFilename);
+
+ numBasis = _numBasis;
+ numQuadPts = _numQuadPts;
+ quadPts = const_cast<double*>(_quadPts);
+ quadWts = const_cast<double*>(_quadWts);
+ basis = const_cast<double*>(_basis);
+ basisDerivRef = const_cast<double*>(_basisDerivRef);
+
+ spatialDBFilename = const_cast<char*>(_spatialDBFilename);
+ id = _id;
+ label = const_cast<char*>(_label);
+
+ dt = _dt;
+ fieldTpdt = const_cast<double*>(_fieldTpdt);
+ fieldT = const_cast<double*>(_fieldT);
+ fieldTmdt = const_cast<double*>(_fieldTmdt);
+
+ spaceDim = _spaceDim;
+ cellDim = _cellDim;
+ numVertices = _numVertices;
+ numCells = _numCells;
+ numCorners = _numCorners;
+ cells = const_cast<int*>(_cells);
+
+ dampingConsts = const_cast<double*>(_dampingConsts);
+ valsResidual = const_cast<double*>(_valsResidual);
+ valsJacobian = const_cast<double*>(_valsJacobian);
+} // constructor
+
+pylith::bc::AbsorbingDampersDataQuad4::~AbsorbingDampersDataQuad4(void)
+{}
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.hh 2007-10-23 00:22:46 UTC (rev 8176)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataQuad4.hh 2007-10-25 21:56:24 UTC (rev 8177)
@@ -0,0 +1,70 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_bc_absorbingdampersdataquad4_hh)
+#define pylith_bc_absorbingdampersdataquad4_hh
+
+#include "AbsorbingDampersData.hh"
+
+namespace pylith {
+ namespace bc {
+ class AbsorbingDampersDataQuad4;
+ } // pylith
+} // bc
+
+class pylith::bc::AbsorbingDampersDataQuad4 : public AbsorbingDampersData
+{
+
+public:
+
+ /// Constructor
+ AbsorbingDampersDataQuad4(void);
+
+ /// Destructor
+ ~AbsorbingDampersDataQuad4(void);
+
+private:
+
+ static const char* _meshFilename;
+
+ static const int _numBasis;
+ static const int _numQuadPts;
+ static const double _quadPts[];
+ static const double _quadWts[];
+ static const double _basis[];
+ static const double _basisDerivRef[];
+
+ static const char* _spatialDBFilename;
+ static const int _id;
+ static const char* _label;
+
+ static const double _dt;
+ static const double _fieldTpdt[];
+ static const double _fieldT[];
+ static const double _fieldTmdt[];
+
+ static const int _spaceDim;
+ static const int _cellDim;
+ static const int _numVertices;
+ static const int _numCells;
+ static const int _numCorners;
+ static const int _cells[];
+
+ static const double _dampingConsts[];
+ static const double _valsResidual[];
+ static const double _valsJacobian[];
+
+};
+
+#endif // pylith_bc_absorbingdampersdataquad4_hh
+
+// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc 2007-10-23 00:22:46 UTC (rev 8176)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc 2007-10-25 21:56:24 UTC (rev 8177)
@@ -32,7 +32,8 @@
0.5,
};
-const char* pylith::bc::AbsorbingDampersDataTri3::_spatialDBFilename = "data/elasticplanestrain.spatialdb";
+const char* pylith::bc::AbsorbingDampersDataTri3::_spatialDBFilename =
+ "data/elasticplanestrain.spatialdb";
const int pylith::bc::AbsorbingDampersDataTri3::_id = 2;
const char* pylith::bc::AbsorbingDampersDataTri3::_label = "bc";
@@ -71,9 +72,9 @@
};
const double pylith::bc::AbsorbingDampersDataTri3::_valsResidual[] = {
0.0, 0.0,
- -6.0e+06, 1.0e+06,
+ -1.2e+07, -2.0e+06,
0.0, 0.0,
- -6.0e+06, 1.0e+06,
+ -1.2e+07, -2.0e+06,
};
const double pylith::bc::AbsorbingDampersDataTri3::_valsJacobian[] = {
0.0, 0.0, // 0x
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py 2007-10-23 00:22:46 UTC (rev 8176)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py 2007-10-25 21:56:24 UTC (rev 8177)
@@ -24,30 +24,30 @@
absorbing dampers on mesh with tri3 cells.
"""
+ dt = 0.25
density = 2500.0
vp = 5000.0
vs = 3000.0
- velocityX = 0.3
- velocityY = 0.2
- edgeLen = 0.5**0.5
- dt = 0.25
- normal = [0.5**0.5, -0.5**0.5]
+ edgeLen = 2.0**0.5
N0 = 0.5
N1 = 0.5
+ velocityX = 0.5*(N0+N1)*(N0*(1.5-1.1)+N1*(2.1-1.3)) / (2.0*dt)
+ velocityY = 0.5*(N0+N1)*(N0*(2.4-1.8)+N1*(2.4-2.2)) / (2.0*dt)
+ normal = [0.5**0.5, -0.5**0.5]
constNormal = density*vp
constTangential = density*vs
dampingConsts = [abs(constNormal*normal[0] - constTangential*normal[1]),
abs(constNormal*normal[1] + constTangential*normal[0])]
- residualX = -dampingConsts[0]*velocityX*edgeLen/(2.0*dt)
- residualY = -dampingConsts[1]*velocityY*edgeLen/(2.0*dt)
+ residualX = -dampingConsts[0]*velocityX*edgeLen
+ residualY = -dampingConsts[1]*velocityY*edgeLen
residual = [residualX, residualY,
residualX, residualY]
- j00 = 2*edgeLen*N0**2 / (2.0*dt)
- j01 = 2*edgeLen*N0*N1 / (2.0*dt)
+ j00 = edgeLen*N0**2 / (2.0*dt)
+ j01 = edgeLen*N0*N1 / (2.0*dt)
j10 = j01
- j11 = 2*edgeLen*N1**2 / (2.0*dt)
+ j11 = edgeLen*N1**2 / (2.0*dt)
jacobian = [dampingConsts[0]*j00, dampingConsts[1]*j00,
dampingConsts[0]*j01, dampingConsts[1]*j01,
dampingConsts[0]*j10, dampingConsts[1]*j10,
@@ -66,7 +66,73 @@
# ----------------------------------------------------------------------
+def calcQuad4():
+ """
+ Calculate damping constants, residual, and Jacobian values for
+ absorbing dampers on mesh with quad4 cells.
+ """
+
+ density = 2500.0
+ vp = 5000.0
+ vs = 3000.0
+ dt = 0.25
+ N0 = 0.5
+ N1 = 0.5
+ velocity1X = 0.5*(N0+N1)*(N0*(1.5-1.1) + N1*(1.2-1.0)) / (2.0*dt)
+ velocity1Y = 0.5*(N0+N1)*(N0*(2.4-1.8) + N1*(1.6-2.4)) / (2.0*dt)
+ normal1 = [-1.0, 0.0]
+
+ velocity2X = (N0*(2.4-1.4) + N1*(2.7-1.5)) / (2.0*dt)
+ velocity2Y = (N0*(2.4-2.4) + N1*(3.4-1.6)) / (2.0*dt)
+ normal2 = [1.0, 0.0]
+
+ edgeLen = 2.0
+
+ constNormal = density*vp
+ constTangential = density*vs
+ dampingConsts = [abs(constNormal*normal1[0] - constTangential*normal1[1]),
+ abs(constNormal*normal1[1] + constTangential*normal1[0]),
+ abs(constNormal*normal2[0] - constTangential*normal2[1]),
+ abs(constNormal*normal2[1] + constTangential*normal2[0])]
+ residual1X = -dampingConsts[0]*velocity1X*edgeLen
+ residual1Y = -dampingConsts[1]*velocity1Y*edgeLen
+ residual2X = -dampingConsts[2]*velocity2X*edgeLen
+ residual2Y = -dampingConsts[3]*velocity2Y*edgeLen
+ residual = [residual1X, residual1Y,
+ residual1X, residual1Y,
+ residual2X, residual2Y,
+ residual2X, residual2Y]
+
+ j00 = edgeLen*N0**2 / (2.0*dt)
+ j01 = edgeLen*N0*N1 / (2.0*dt)
+ j10 = j01
+ j11 = edgeLen*N1**2 / (2.0*dt)
+ jacobian = [dampingConsts[0]*j00, dampingConsts[1]*j00,
+ dampingConsts[0]*j01, dampingConsts[1]*j01,
+ dampingConsts[0]*j10, dampingConsts[1]*j10,
+ dampingConsts[0]*j11, dampingConsts[1]*j11,
+ dampingConsts[2]*j00, dampingConsts[3]*j00,
+ dampingConsts[2]*j01, dampingConsts[3]*j01,
+ dampingConsts[2]*j10, dampingConsts[3]*j10,
+ dampingConsts[2]*j11, dampingConsts[3]*j11]
+ print "Absorbing boundary for quad4mesh"
+ print "damping constants:"
+ for v in dampingConsts:
+ print " %16.8e" % v
+ print "velocity:"
+ print " velocity1: ",velocity1X," ",velocity1Y
+ print " velocity2: ",velocity2X," ",velocity2Y
+ print "values for residual:"
+ for v in residual:
+ print " %16.8e" % v
+ print "values for jacobian:"
+ for j in jacobian:
+ print " %16.8e" % j
+
+
+# ----------------------------------------------------------------------
calcTri3()
+calcQuad4()
# End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4.mesh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4.mesh 2007-10-23 00:22:46 UTC (rev 8176)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/quad4.mesh 2007-10-25 21:56:24 UTC (rev 8177)
@@ -1,3 +1,19 @@
+// Original mesh
+//
+// 1 ----- 3 ----- 5
+// | | |
+// | 0 | 1 |
+// | | |
+// 0 ----- 2 ----- 4
+//
+// Sieve mesh
+//
+// 3 ----- 5 ----- 7
+// | | |
+// | 0 | 1 |
+// | | |
+// 2 ----- 4 ----- 6
+//
mesh = {
dimension = 2
use-index-zero = true
@@ -33,4 +49,12 @@
0 1 4
}
}
+ group = {
+ name = bc2
+ type = vertices
+ count = 4
+ indices = {
+ 0 1 4 5
+ }
+ }
}
More information about the cig-commits
mailing list