[cig-commits] r8175 - in short/3D/PyLith/trunk: libsrc/bc
unittests/libtests/bc unittests/libtests/bc/data
brad at geodynamics.org
brad at geodynamics.org
Mon Oct 22 16:50:29 PDT 2007
Author: brad
Date: 2007-10-22 16:50:29 -0700 (Mon, 22 Oct 2007)
New Revision: 8175
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/AbsorbingDampersDataLine2.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py
Log:
Fixed bugs associated with getting cells/vertices for boundary mesh. Fixed incorrect values for testing absorbing dampers with tri3 cells.
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2007-10-22 23:49:47 UTC (rev 8174)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2007-10-22 23:50:29 UTC (rev 8175)
@@ -66,7 +66,7 @@
throw std::runtime_error(msg.str());
} // if
- _boundaryMesh->view("ABSORBING BOUNDARY MESH");
+ //_boundaryMesh->view("ABSORBING BOUNDARY MESH");
// check compatibility of quadrature and boundary mesh
if (_quadrature->cellDim() != _boundaryMesh->getDimension()) {
@@ -81,17 +81,19 @@
const int numCorners = _quadrature->numBasis();
const ALE::Obj<ALE::Mesh::label_sequence>& cells =
- _boundaryMesh->heightStratum(0);
+ _boundaryMesh->heightStratum(1);
+
assert(!cells.isNull());
const Mesh::label_sequence::iterator cellsBegin = cells->begin();
const Mesh::label_sequence::iterator cellsEnd = cells->end();
const ALE::Obj<sieve_type>& sieve = _boundaryMesh->getSieve();
+ const int boundaryDepth = _boundaryMesh->depth()-1; // depth of bndry cells
assert(!sieve.isNull());
for (Mesh::label_sequence::iterator c_iter=cellsBegin;
c_iter != cellsEnd;
++c_iter) {
- const int cellNumCorners = sieve->nCone(*c_iter,
- _boundaryMesh->depth())->size();
+ const int cellNumCorners = (_boundaryMesh->getDimension() > 0) ?
+ sieve->nCone(*c_iter, boundaryDepth)->size() : 1;
if (numCorners != cellNumCorners) {
std::ostringstream msg;
msg << "Quadrature is incompatible with cell for absorbing boundary "
@@ -147,7 +149,7 @@
double_array dampingConstsGlobal(fiberDim);
const ALE::Obj<real_section_type>& coordinates =
- _boundaryMesh->getRealSection("coordinates");
+ mesh->getRealSection("coordinates");
assert(!coordinates.isNull());
for(Mesh::label_sequence::iterator c_iter = cells->begin();
@@ -182,11 +184,15 @@
cellGeometry.jacobian(&jacobian, &jacobianDet, cellVertices, quadPtRef);
cellGeometry.orientation(&orientation, jacobian, jacobianDet,
upDir);
+ orientation /= jacobianDet;
dampingConstsGlobal = 0.0;
- for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim)
dampingConstsGlobal[iDim] +=
dampingConstsLocal[jDim]*orientation[iDim*spaceDim+jDim];
+ // Ensure damping constants are positive
+ dampingConstsGlobal[iDim] = fabs(dampingConstsGlobal[iDim]);
+ } // for
} // for
_dampingConsts->updatePoint(*c_iter, &dampingConstsGlobal[0]);
} // for
@@ -213,13 +219,13 @@
// Get cell information
const ALE::Obj<ALE::Mesh::label_sequence>& cells =
- _boundaryMesh->heightStratum(0);
+ _boundaryMesh->heightStratum(1);
assert(!cells.isNull());
const Mesh::label_sequence::iterator cellsEnd = cells->end();
// Get sections
const ALE::Obj<real_section_type>& coordinates =
- _boundaryMesh->getRealSection("coordinates");
+ mesh->getRealSection("coordinates");
assert(!coordinates.isNull());
const ALE::Obj<real_section_type>& dispTpdt = fields->getHistoryItem(0);
const ALE::Obj<real_section_type>& dispTmdt = fields->getHistoryItem(2);
@@ -306,13 +312,13 @@
// Get cell information
const ALE::Obj<ALE::Mesh::label_sequence>& cells =
- _boundaryMesh->heightStratum(0);
+ _boundaryMesh->heightStratum(1);
assert(!cells.isNull());
const Mesh::label_sequence::iterator cellsEnd = cells->end();
// Get sections
const ALE::Obj<real_section_type>& coordinates =
- _boundaryMesh->getRealSection("coordinates");
+ mesh->getRealSection("coordinates");
assert(!coordinates.isNull());
const ALE::Obj<real_section_type>& dispT = fields->getHistoryItem(1);
assert(!dispT.isNull());
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2007-10-22 23:49:47 UTC (rev 8174)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2007-10-22 23:50:29 UTC (rev 8175)
@@ -77,7 +77,7 @@
const int cellDim = boundaryMesh->getDimension();
const ALE::Obj<sieve_type>& sieve = boundaryMesh->getSieve();
- const ALE::Obj<Mesh::label_sequence>& cells = boundaryMesh->heightStratum(0);
+ const ALE::Obj<Mesh::label_sequence>& cells = boundaryMesh->heightStratum(1);
const int numVertices = boundaryMesh->depthStratum(0)->size();
const int numCells = cells->size();
@@ -85,13 +85,15 @@
CPPUNIT_ASSERT_EQUAL(_data->numVertices, numVertices);
CPPUNIT_ASSERT_EQUAL(_data->numCells, numCells);
- const int numCorners = sieve->nCone(*cells->begin(),
- boundaryMesh->depth())->size();
- CPPUNIT_ASSERT_EQUAL(_data->numCorners, numCorners);
+ const int boundaryDepth = boundaryMesh->depth()-1; // depth of bndry cells
int iCell = 0;
for(Mesh::label_sequence::iterator c_iter = cells->begin();
c_iter != cells->end();
++c_iter) {
+ const int numCorners = (boundaryMesh->getDimension() > 0) ?
+ sieve->nCone(*c_iter, boundaryDepth)->size() : 1;
+ CPPUNIT_ASSERT_EQUAL(_data->numCorners, numCorners);
+
const ALE::Obj<sieve_type::traits::coneSequence>& cone =
sieve->cone(*c_iter);
for(sieve_type::traits::coneSequence::iterator v_iter = cone->begin();
@@ -233,7 +235,7 @@
// Setup mesh
meshio::MeshIOAscii iohandler;
iohandler.filename(_data->meshFilename);
- iohandler.debug(true);
+ //iohandler.debug(true);
iohandler.read(mesh);
CPPUNIT_ASSERT(!mesh->isNull());
(*mesh)->getFactory()->clear();
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.cc 2007-10-22 23:49:47 UTC (rev 8174)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataLine2.cc 2007-10-22 23:50:29 UTC (rev 8175)
@@ -35,7 +35,7 @@
const char* pylith::bc::AbsorbingDampersDataLine2::_label = "bc0";
const double pylith::bc::AbsorbingDampersDataLine2::_dt = 0.25;
-const double pylith::bc::AbsorbingDampersDataLine2::_fieldTpdt[] = {
+const double pylith::bc::AbsorbingDampersDataLine2::_fieldTmdt[] = {
1.0,
1.1,
1.2,
@@ -45,7 +45,7 @@
1.3,
1.5,
};
-const double pylith::bc::AbsorbingDampersDataLine2::_fieldTmdt[] = {
+const double pylith::bc::AbsorbingDampersDataLine2::_fieldTpdt[] = {
1.2,
1.5,
1.8,
@@ -67,9 +67,9 @@
17.5e+6,
};
const double pylith::bc::AbsorbingDampersDataLine2::_valsResidual[] = {
- 12.5e+6*0.2/0.5,
+ -12.5e+6*0.2/0.5,
0.0,
- 17.5e+6*0.6/0.5,
+ -17.5e+6*0.6/0.5,
};
const double pylith::bc::AbsorbingDampersDataLine2::_valsJacobian[] = {
12.5e+6/0.5, 0.0, 0.0,
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc 2007-10-22 23:49:47 UTC (rev 8174)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTri3.cc 2007-10-22 23:50:29 UTC (rev 8175)
@@ -37,23 +37,23 @@
const char* pylith::bc::AbsorbingDampersDataTri3::_label = "bc";
const double pylith::bc::AbsorbingDampersDataTri3::_dt = 0.25;
-const double pylith::bc::AbsorbingDampersDataTri3::_fieldTpdt[] = {
- 1.0,
- 1.1,
- 1.2,
- 1.3,
+const double pylith::bc::AbsorbingDampersDataTri3::_fieldTmdt[] = {
+ 1.0, 2.4,
+ 1.1, 1.8,
+ 1.2, 2.4,
+ 1.3, 2.2,
};
const double pylith::bc::AbsorbingDampersDataTri3::_fieldT[] = {
- 1.1,
- 1.3,
- 1.5,
- 1.7,
+ 1.1, 2.0,
+ 1.3, 2.1,
+ 1.5, 2.2,
+ 1.7, 2.3,
};
-const double pylith::bc::AbsorbingDampersDataTri3::_fieldTmdt[] = {
- 1.2,
- 1.5,
- 1.8,
- 2.1,
+const double pylith::bc::AbsorbingDampersDataTri3::_fieldTpdt[] = {
+ 1.2, 1.6,
+ 1.5, 2.4,
+ 1.8, 2.0,
+ 2.1, 2.4,
};
const int pylith::bc::AbsorbingDampersDataTri3::_spaceDim = 2;
@@ -67,31 +67,47 @@
const double pylith::bc::AbsorbingDampersDataTri3::_dampingConsts[] = {
- 12.5e+6, 7.5e+6,
+ 1.41421356e+07, 3.53553391e+06
};
const double pylith::bc::AbsorbingDampersDataTri3::_valsResidual[] = {
0.0, 0.0,
- 2.4e+07, 6.0e+06,
+ -6.0e+06, 1.0e+06,
0.0, 0.0,
- 2.4e+07, 6.0e+06
+ -6.0e+06, 1.0e+06,
};
const double pylith::bc::AbsorbingDampersDataTri3::_valsJacobian[] = {
- 0.0, 0.0, // 0 // FIX THESE
+ 0.0, 0.0, // 0x
0.0, 0.0,
0.0, 0.0,
0.0, 0.0,
- 0.0, 0.0, // 1
+ 0.0, 0.0, // 0y
0.0, 0.0,
0.0, 0.0,
0.0, 0.0,
- 0.0, 0.0, // 2
+ 0.0, 0.0, // 1x
+ 1.0e+07, 0.0,
0.0, 0.0,
+ 1.0e+07, 0.0,
+ 0.0, 0.0, // 1y
+ 0.0, 2.5e+06,
0.0, 0.0,
+ 0.0, 2.5e+06,
+ 0.0, 0.0, // 2x
0.0, 0.0,
- 0.0, 0.0, // 3
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, // 3x
+ 1.0e+07, 0.0,
+ 0.0, 0.0,
+ 1.0e+07, 0.0,
+ 0.0, 0.0, // 3y
+ 0.0, 2.5e+06,
+ 0.0, 0.0,
+ 0.0, 2.5e+06,
};
pylith::bc::AbsorbingDampersDataTri3::AbsorbingDampersDataTri3(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py 2007-10-22 23:49:47 UTC (rev 8174)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/absorbingdampers.py 2007-10-22 23:50:29 UTC (rev 8175)
@@ -27,28 +27,42 @@
density = 2500.0
vp = 5000.0
vs = 3000.0
- velocity = 1.2
+ velocityX = 0.3
+ velocityY = 0.2
edgeLen = 0.5**0.5
dt = 0.25
normal = [0.5**0.5, -0.5**0.5]
+ N0 = 0.5
+ N1 = 0.5
- dampingConsts = [density*vp, density*vs]
- residualNormal = dampingConsts[0]*velocity*edgeLen/(2.0*dt)
- residualTangential = dampingConsts[1]*velocity*edgeLen/(2.0*dt)
- residual = [residualNormal*normal[0] - residualTangential*normal[1],
- residualNormal*normal[1] + residualTangential*normal[0],
- residualNormal*normal[0] - residualTangential*normal[1],
- residualNormal*normal[1] + residualTangential*normal[0]]
+ 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)
+ residual = [residualX, residualY,
+ residualX, residualY]
- jacobian = []
+ j00 = 2*edgeLen*N0**2 / (2.0*dt)
+ j01 = 2*edgeLen*N0*N1 / (2.0*dt)
+ j10 = j01
+ j11 = 2*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]
print "Absorbing boundary for tri3 mesh"
print "damping constants:"
for v in dampingConsts:
print " %16.8e" % v
- print "residual:"
+ print "values for residual:"
for v in residual:
print " %16.8e" % v
+ print "values for jacobian:"
+ for j in jacobian:
+ print " %16.8e" % j
# ----------------------------------------------------------------------
More information about the cig-commits
mailing list