[cig-commits] r8220 - in short/3D/PyLith/trunk: libsrc/bc
unittests/libtests/bc unittests/libtests/bc/data
willic3 at geodynamics.org
willic3 at geodynamics.org
Tue Nov 6 18:34:02 PST 2007
Author: willic3
Date: 2007-11-06 18:34:01 -0800 (Tue, 06 Nov 2007)
New Revision: 8220
Modified:
short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.hh
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.hh
Log:
Fixes to Neumann BC and corresponding changes to unit tests.
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2007-11-07 00:33:42 UTC (rev 8219)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2007-11-07 02:34:01 UTC (rev 8220)
@@ -17,12 +17,9 @@
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
#include "pylith/feassemble/CellGeometry.hh" // USES CellGeometry
#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
-#include "pylith/utils/array.hh" // USES double_array
-#include <petscmat.h> // USES PETSc Mat
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
-#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include <Selection.hh> // USES submesh algorithms
@@ -69,6 +66,8 @@
throw std::runtime_error(msg.str());
} // if
+ //_boundaryMesh->view("TRACTION BOUNDARY MESH");
+
// check compatibility of quadrature and boundary mesh
if (_quadrature->cellDim() != _boundaryMesh->getDimension()) {
std::ostringstream msg;
@@ -82,19 +81,23 @@
const int numCorners = _quadrature->numBasis();
// Get 'surface' cells (1 dimension lower than top-level cells)
- const ALE::Obj<sieve_type>& sieve = _boundaryMesh->getSieve();
const ALE::Obj<Mesh::label_sequence>& cells = _boundaryMesh->heightStratum(1);
+ assert(!cells.isNull());
+
const Mesh::label_sequence::iterator cellsBegin = cells->begin();
const Mesh::label_sequence::iterator cellsEnd = cells->end();
- std::cout << "cellsBegin: " << *cellsBegin << std::endl;
- std::cout << "cellsEnd: " << *cellsEnd << std::endl;
+ // std::cout << "cellsBegin: " << *cellsBegin << std::endl;
+ // std::cout << "cellsEnd: " << *cellsEnd << std::endl;
+ const ALE::Obj<sieve_type>& sieve = _boundaryMesh->getSieve();
+ const int boundaryDepth = _boundaryMesh->depth()-1; //depth of boundary cells
+ assert(!sieve.isNull());
// Make sure surface cells are compatible with quadrature.
- 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()-1)->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 Neumann traction "
@@ -112,10 +115,11 @@
const int numQuadPts = _quadrature->numQuadPts();
const int spaceDim = cs->spaceDim();
const int fiberDim = spaceDim * numQuadPts;
- _tractionGlobal = new real_section_type(mesh->comm(), mesh->debug());
- assert(!_tractionGlobal.isNull());
- _tractionGlobal->setFiberDimension(cells, fiberDim);
- mesh->allocate(_tractionGlobal);
+ _tractionsGlobal = new real_section_type(_boundaryMesh->comm(),
+ _boundaryMesh->debug());
+ assert(!_tractionsGlobal.isNull());
+ _tractionsGlobal->setFiberDimension(cells, fiberDim);
+ _boundaryMesh->allocate(_tractionsGlobal);
// Containers for orientation information
const int orientationSize = spaceDim * spaceDim;
@@ -126,11 +130,6 @@
double_array orientation(orientationSize);
double_array cellVertices(numBasis*spaceDim);
- // Get mesh coordinates.
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- // coordinates->view("Mesh coordinates from Neumann::initialize");
-
// open database with traction information
// NEED TO SET NAMES BASED ON DIMENSION OF BOUNDARY
// 1-D problem = {'normal-traction'}
@@ -163,12 +162,18 @@
// Containers for database query results and quadrature coordinates in
// reference geometry.
double_array tractionDataLocal(spaceDim);
+ double_array quadPtRef(spaceDim);
const double_array& quadPtsRef = _quadrature->quadPtsRef();
- // double_array quadPtRef(spaceDim);
// Container for cell tractions rotated to global coordinates.
double_array cellTractionsGlobal(fiberDim);
+ // Get mesh coordinates.
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ // coordinates->view("Mesh coordinates from Neumann::initialize");
+
// Loop over cells in boundary mesh, compute orientations, and then
// compute corresponding traction vector in global coordinates
// (store values in _tractionGlobal).
@@ -207,20 +212,11 @@
// Compute Jacobian and determinant at quadrature point, then get
// orientation.
- double_array quadPtRef(&quadPtsRef[iRef], cellDim);
+ memcpy(&quadPtRef[0], &quadPtsRef[iRef], cellDim*sizeof(double));
cellGeometry.jacobian(&jacobian, &jacobianDet, cellVertices, quadPtRef);
cellGeometry.orientation(&orientation, jacobian, jacobianDet, upDir);
+ orientation /= jacobianDet;
- // Normalize orientation.
- for(int iDim=0; iDim < spaceDim; ++iDim) {
- double mag = 0.0;
- for(int jDim=0; jDim < spaceDim; ++jDim)
- mag += orientation[iDim*spaceDim+jDim]*orientation[iDim*spaceDim+jDim];
- mag = sqrt(mag);
- for(int jDim=0; jDim < spaceDim; ++jDim)
- orientation[iDim*spaceDim+jDim] /= mag;
- } // for
-
// Rotate traction vector from local coordinate system to global
// coordinate system
for(int iDim = 0; iDim < spaceDim; ++iDim) {
@@ -230,10 +226,10 @@
} // for
} // for
- // Update tractionGlobal
- _boundaryMesh->update(_tractionGlobal, *c_iter, &cellTractionsGlobal[0]);
+ // Update tractionsGlobal
+ _tractionsGlobal->updatePoint(*c_iter, &cellTractionsGlobal[0]);
} // for
- _tractionGlobal->view("Global tractions from Neumann::initialize");
+ // _tractionsGlobal->view("Global tractions from Neumann::initialize");
_db->close();
} // initialize
@@ -278,7 +274,7 @@
// Allocate vectors for cell values.
_initCellVector();
const int cellVecSize = numBasis*spaceDim;
- double_array tractionCell(numQuadPts*spaceDim);
+ double_array tractionsCell(numQuadPts*spaceDim);
// Loop over faces and integrate contribution from each face
for (Mesh::label_sequence::iterator c_iter=cellsBegin;
@@ -291,8 +287,8 @@
_resetCellVector();
// Restrict tractions to cell
- _boundaryMesh->restrict(_tractionGlobal, *c_iter,
- &tractionCell[0], tractionCell.size());
+ _boundaryMesh->restrict(_tractionsGlobal, *c_iter,
+ &tractionsCell[0], tractionsCell.size());
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -307,14 +303,14 @@
const double valIJ = valI * basis[iQuad*numBasis+jBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] +=
- tractionCell[iQuad*spaceDim+iDim] * valIJ;
+ tractionsCell[iQuad*spaceDim+iDim] * valIJ;
} // for
} // for
} // for
PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(1+numBasis*(1+2*spaceDim))));
// Assemble cell contribution into field
- mesh->updateAdd(residual, *c_iter, _cellVector);
+ _boundaryMesh->updateAdd(residual, *c_iter, _cellVector);
} // for
} // integrateResidual
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.hh 2007-11-07 00:33:42 UTC (rev 8219)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.hh 2007-11-07 02:34:01 UTC (rev 8220)
@@ -22,8 +22,8 @@
#include "BoundaryCondition.hh" // ISA BoundaryCondition
#include "pylith/feassemble/Integrator.hh" // ISA Integrator
-// #include "pylith/utils/array.hh" // USES std::vector, double_array, int_array
-// #include "pylith/utils/sievetypes.hh" // USES real_section_type
+#include "pylith/utils/array.hh" // USES std::vector, double_array, int_array
+#include "pylith/utils/sievetypes.hh" // USES real_section_type
/// Namespace for pylith package
namespace pylith {
@@ -107,7 +107,7 @@
ALE::Obj<Mesh> _boundaryMesh;
/// Traction vector in global coordinates at integration points.
- ALE::Obj<real_section_type> _tractionGlobal;
+ ALE::Obj<real_section_type> _tractionsGlobal;
}; // class Neumann
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2007-11-07 00:33:42 UTC (rev 8219)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2007-11-07 02:34:01 UTC (rev 8220)
@@ -101,7 +101,7 @@
// Check traction values
const int numQuadPts = _data->numQuadPts;
const int fiberDim = numQuadPts * spaceDim;
- double_array tractionCell(fiberDim);
+ double_array tractionsCell(fiberDim);
int index = 0;
const double tolerance = 1.0e-06;
@@ -109,13 +109,13 @@
c_iter != cells->end();
++c_iter) {
- bc._boundaryMesh->restrict(bc._tractionGlobal, *c_iter,
- &tractionCell[0], tractionCell.size());
+ bc._boundaryMesh->restrict(bc._tractionsGlobal, *c_iter,
+ &tractionsCell[0], tractionsCell.size());
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
for (int iDim =0; iDim < spaceDim; ++iDim) {
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->tractionCell[index],
- tractionCell[iQuad*spaceDim+iDim],
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->tractionsCell[index],
+ tractionsCell[iQuad*spaceDim+iDim],
tolerance);
++index;
} // for
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.cc 2007-11-07 00:33:42 UTC (rev 8219)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.cc 2007-11-07 02:34:01 UTC (rev 8220)
@@ -32,7 +32,7 @@
numVertices(0),
numCorners(0),
cells(0),
- tractionCell(0),
+ tractionsCell(0),
valsResidual(0)
{ // constructor
} // constructor
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.hh 2007-11-07 00:33:42 UTC (rev 8219)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.hh 2007-11-07 02:34:01 UTC (rev 8220)
@@ -62,7 +62,7 @@
int numVertices; ///< Expected number of vertices in the mesh.
int* numCorners; ///< Expected number of vertices for each boundary cell.
int* cells; ///< Expected array of vertices defining each boundary cell.
- double* tractionCell; ///< Expected traction values at quadrature points.
+ double* tractionsCell; ///< Expected traction values at quadrature points.
double* valsResidual; ///< Expected residual at each vertex.
//@}
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.cc 2007-11-07 00:33:42 UTC (rev 8219)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.cc 2007-11-07 02:34:01 UTC (rev 8220)
@@ -82,14 +82,14 @@
const int pylith::bc::NeumannDataHex8::_numCorners[] = { 4, 4 };
const int pylith::bc::NeumannDataHex8::_cells[] = { 0, 4, 6, 2,
4, 8, 10, 6 };
-const double pylith::bc::NeumannDataHex8::_tractionCell[] = { 1.0, 0.0, 0.0,
- 1.0, 0.0, 0.0,
- 1.0, 0.0, 0.0,
- 1.0, 0.0, 0.0,
- 1.0, 0.0, 0.0,
- 1.0, 0.0, 0.0,
- 1.0, 0.0, 0.0,
- 1.0, 0.0, 0.0};
+const double pylith::bc::NeumannDataHex8::_tractionsCell[] = { 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0};
const double pylith::bc::NeumannDataHex8::_valsResidual[] = { 2.0, 0.0, 0.0,
0.0, 0.0, 0.0,
2.0, 0.0, 0.0,
@@ -123,7 +123,7 @@
numVertices = _numVertices;
numCorners = const_cast<int*>(_numCorners);
cells = const_cast<int*>(_cells);
- tractionCell = const_cast<double*>(_tractionCell);
+ tractionsCell = const_cast<double*>(_tractionsCell);
valsResidual = const_cast<double*>(_valsResidual);
} // constructor
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.hh 2007-11-07 00:33:42 UTC (rev 8219)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.hh 2007-11-07 02:34:01 UTC (rev 8220)
@@ -58,7 +58,7 @@
static const int _numVertices; ///< Expected number of vertices in the mesh.
static const int _numCorners[]; ///< Expected number of vertices for each boundary cell.
static const int _cells[]; ///< Expected array of vertices defining each boundary cell.
- static const double _tractionCell[]; ///< Expected traction values at quadrature points.
+ static const double _tractionsCell[]; ///< Expected traction values at quadrature points.
static const double _valsResidual[]; ///< Expected residual at each vertex.
};
More information about the cig-commits
mailing list