[cig-commits] r7781 - in short/3D/PyLith/trunk: libsrc/bc
unittests/libtests/bc
willic3 at geodynamics.org
willic3 at geodynamics.org
Tue Aug 7 08:48:25 PDT 2007
Author: willic3
Date: 2007-08-07 08:48:24 -0700 (Tue, 07 Aug 2007)
New Revision: 7781
Modified:
short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannHex8.hh
Log:
More work and debugging on traction BC.
Stuff checked in is broken right now.
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2007-08-07 00:47:49 UTC (rev 7780)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2007-08-07 15:48:24 UTC (rev 7781)
@@ -19,6 +19,7 @@
#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
@@ -81,12 +82,14 @@
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);
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;
// Make sure surface cells are compatible with quadrature.
- const ALE::Obj<sieve_type>& sieve = _boundaryMesh->getSieve();
assert(!sieve.isNull());
for (Mesh::label_sequence::iterator c_iter=cellsBegin;
c_iter != cellsEnd;
@@ -126,6 +129,7 @@
// 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
@@ -161,7 +165,6 @@
double_array tractionDataLocal(spaceDim);
const double_array& quadPtsRef = _quadrature->quadPtsRef();
// double_array quadPtRef(spaceDim);
- const double_array& quadPts = _quadrature->quadPts();
// Container for cell tractions rotated to global coordinates.
double_array cellTractionsGlobal(fiberDim);
@@ -172,8 +175,20 @@
for(Mesh::label_sequence::iterator c_iter = cellsBegin;
c_iter != cellsEnd;
++c_iter) {
+ std::cout << "c_iter: " << *c_iter << std::endl;
_quadrature->computeGeometry(_boundaryMesh, coordinates, *c_iter);
- mesh->restrict(coordinates, *c_iter, &cellVertices[0], cellVertices.size());
+ const double_array& quadPts = _quadrature->quadPts();
+ const real_section_type::value_type* cellVert = mesh->restrict(coordinates, *c_iter);
+ // For some reason this method doesn' work. I need to find out why.
+ // mesh->restrict(coordinates, *c_iter, &cellVertices[0], cellVertices.size());
+ std::cout << "cellVertices: " << std::endl;
+ for(int iTest = 0; iTest < numBasis; ++iTest) {
+ for(int iDim = 0; iDim < spaceDim; ++iDim) {
+ cellVertices[iDim+spaceDim*iTest] = cellVert[iDim+spaceDim*iTest];
+ std::cout << " " << cellVertices[iDim+spaceDim*iTest];
+ } // for
+ std::cout << std::endl;
+ } // for
for(int iQuad = 0, iRef=0, iSpace=0; iQuad < numQuadPts; ++iQuad, iRef+=cellDim, iSpace+=spaceDim) {
// Get traction vector in local coordinate system at quadrature point
@@ -196,6 +211,16 @@
cellGeometry.jacobian(&jacobian, &jacobianDet, cellVertices, quadPtRef);
cellGeometry.orientation(&orientation, jacobian, jacobianDet, upDir);
+ // 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
cellTractionsGlobal = 0.0;
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2007-08-07 00:47:49 UTC (rev 7780)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2007-08-07 15:48:24 UTC (rev 7781)
@@ -21,6 +21,7 @@
#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+#include <Selection.hh> // USES submesh algorithms
#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
@@ -72,25 +73,33 @@
const ALE::Obj<Mesh::label_sequence>& cells = bc._boundaryMesh->heightStratum(1);
const int numBoundaryCells = cells->size();
CPPUNIT_ASSERT_EQUAL(_data->numBoundaryCells, numBoundaryCells);
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ coordinates->view("Mesh coordinates from TestNeumann::testInitialize");
+
+ const int spaceDim = _data->spaceDim;
+ const int numBasis = bc._quadrature->numBasis();
+ double_array cellVertices(numBasis*spaceDim);
int iCell = 0;
- int i = 0;
for(Mesh::label_sequence::iterator c_iter = cells->begin();
c_iter != cells->end();
++c_iter) {
const int numCorners = sieve->nCone(*c_iter, mesh->depth())->size();
CPPUNIT_ASSERT_EQUAL(_data->numCorners[iCell++], numCorners);
- const ALE::Obj<sieve_type::traits::coneSequence>& cone =
- sieve->cone(*c_iter);
- for(sieve_type::traits::coneSequence::iterator v_iter = cone->begin();
- v_iter != cone->end();
- ++v_iter)
- CPPUNIT_ASSERT_EQUAL(_data->cells[i++], *v_iter);
+ const real_section_type::value_type* cellVert = mesh->restrict(coordinates, *c_iter);
+ std::cout << "c_iter " << *c_iter << " vertex coords:" << std::endl;
+ for(int iVert = 0; iVert < numCorners; ++iVert) {
+ for(int iDim = 0; iDim < spaceDim; ++iDim) {
+ cellVertices[iDim+spaceDim*iVert] = cellVert[iDim+spaceDim*iVert];
+ std::cout << " " << cellVertices[iDim+spaceDim*iVert];
+ } // for
+ std::cout << std::endl;
+ } // for
} // for
// Check traction values
- int numQuadPts = _data->numQuadPts;
- int spaceDim = _data->spaceDim;
- int fiberDim = numQuadPts * spaceDim;
+ const int numQuadPts = _data->numQuadPts;
+ const int fiberDim = numQuadPts * spaceDim;
double_array tractionCell(fiberDim);
int index = 0;
const double tolerance = 1.0e-06;
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannHex8.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannHex8.hh 2007-08-07 00:47:49 UTC (rev 7780)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumannHex8.hh 2007-08-07 15:48:24 UTC (rev 7781)
@@ -36,7 +36,7 @@
// CPPUNIT TEST SUITE /////////////////////////////////////////////////
CPPUNIT_TEST_SUITE( TestNeumannHex8 );
- // CPPUNIT_TEST( testInitialize );
+ CPPUNIT_TEST( testInitialize );
// CPPUNIT_TEST( testIntegrateResidual );
CPPUNIT_TEST_SUITE_END();
More information about the cig-commits
mailing list