[cig-commits] r4571 - in short/3D/PyLith/trunk: libsrc/feassemble
unittests/libtests/feassemble unittests/libtests/feassemble/data
baagaard at geodynamics.org
baagaard at geodynamics.org
Tue Sep 19 20:54:08 PDT 2006
Author: baagaard
Date: 2006-09-19 20:54:07 -0700 (Tue, 19 Sep 2006)
New Revision: 4571
Added:
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/quadrature1d.py
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature1D.cc
Log:
Reordered basis and basisDeriv in Quadrature. Implemented test for quadratic 1-D quadrature.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2006-09-19 19:50:58 UTC (rev 4570)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2006-09-20 03:54:07 UTC (rev 4571)
@@ -61,26 +61,26 @@
* weights of the quadrature points.
*
* @param basis Array of basis functions evaluated at quadrature pts
- * N0Qp0, N0Qp1, ...
- * N1Qp0, N1Qp1, ...
+ * N0Qp0, N1Qp0, ...
+ * N0Qp1, N1Qp1, ...
* ...
- * size = numCorners * numQuadPts
- * index = iBasis*numQuadPts + iQuadPt
+ * size = numQuadPts * numCorners
+ * index = iQuadPt*numCorners + iBasis
*
* @param basisDeriv Array of basis function derivaties evaluated
* at quadrature pts
- * N0xQp0, N0yQp0, N0zQp0, N0xQp1, N0yQp1, N0zQp1, ...
- * N1xQp0, N1yQp0, N1zQp0, N1xQp1, N1yQp1, N1zQp1, ...
+ * N0xQp0, N0yQp0, N0zQp0, N1xQp0, N1yQp0, N1zQp0, ...
+ * N0xQp1, N0yQp1, N0zQp1, N1xQp1, N1yQp1, N1zQp1, ...
* ...
- * size = numCorners * numQuadPts * spaceDim
- * index = iVertex*numQuadPts*spaceDim + iQuadPt*spaceDim + iDim
+ * size = numCorners * numQuadPts * spaceDim
+ * index = iQuadPt*numCorners*spaceDim + iBasis*spaceDim + iDim
*
* @param quadPts Array of coordinates of quadrature points in
* reference cell
* Qp0x, Qp0y, Qp0z
* Qp1x, Qp1y, Qp1z
* size = numQuadPts * numDims
- * index = iQuadPts*numDims + iDim
+ * index = iQuadPt*numDims + iDim
*
* @param quadWts Array of weights of quadrature points
* WtQp0, WtQp1, ...
@@ -142,21 +142,21 @@
/** Array of basis functions evaluated at the quadrature points.
*
- * N0Qp0, N0Qp1, ...
- * N1Qp0, N1Qp1, ...
+ * N0Qp0, N1Qp0, ...
+ * N0Qp1, N1Qp1, ...
*
- * size = numCorners * numQuadPts
- * index = iBasis*numQuadPts + iQuadPt
+ * size = numQuadPts * numCorners
+ * index = iQuadPt*numCorners + iBasis
*/
double* _basis; ///< Array of basis fns evaluated at quad pts
/** Array of basis functions evaluated at the quadrature points.
*
- * N0xQp0, N0yQp0, N0zQp0, N0xQp1, N0yQp1, N0zQp1, ...
- * N1xQp0, N1yQp0, N1zQp0, N1xQp1, N1yQp1, N1zQp1, ...
+ * N0xQp0, N0yQp0, N0zQp0, N1xQp0, N1yQp0, N1zQp0, ...
+ * N0xQp1, N0yQp1, N0zQp1, N1xQp1, N1yQp1, N1zQp1, ...
*
- * size = numCorners * numQuadPts * spaceDim
- * index = iVertex*numQuadPts*spaceDim + iQuadPt*spaceDim + iDim
+ * size = numQuadPts * numCorners * spaceDim
+ * index = iQuadPt*numCorners*spaceDim + iBasis*spaceDim + iDim
*/
double* _basisDeriv;
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc 2006-09-19 19:50:58 UTC (rev 4570)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc 2006-09-20 03:54:07 UTC (rev 4571)
@@ -67,21 +67,22 @@
_quadPts[iQuadPt] = 0.0;
for (int iVertex=0; iVertex < _numCorners; ++iVertex)
_quadPts[iQuadPt] +=
- _basis[iVertex*_numQuadPts+iQuadPt]*vertCoords[iVertex];
+ _basis[iQuadPt*_numCorners+iVertex]*vertCoords[iVertex];
// Compute Jacobian at quadrature point
// J = dx/dp = sum[i=1,n] (dNi/dp * xi)
_jacobian[iQuadPt] = 0.0;
for (int iVertex=0; iVertex < _numCorners; ++iVertex)
_jacobian[iQuadPt] +=
- _basisDeriv[iVertex*_numQuadPts+iQuadPt] * vertCoords[iVertex];
+ _basisDeriv[iQuadPt*_numCorners+iVertex] * vertCoords[iVertex];
// Compute determinant of Jacobian at quadrature point
// |J| = j00
const double det = _jacobian[iQuadPt];
if (det < _jacobianTol) {
std::ostringstream msg;
- msg << "Determinant of Jacobian is below minimum permissible value!\n";
+ msg << "Determinant of Jacobian (" << det << ") is below minimum\n"
+ << "permissible value (" << _jacobianTol << ")!\n";
throw std::runtime_error(msg.str());
} // for
_jacobianDet[iQuadPt] = _jacobian[iQuadPt];
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature1D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature1D.cc 2006-09-19 19:50:58 UTC (rev 4570)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature1D.cc 2006-09-20 03:54:07 UTC (rev 4571)
@@ -34,20 +34,21 @@
void
pylith::feassemble::TestQuadrature1D::testLinear(void)
{ // testLinear
- // Basis fns: N0 = 1-p, N1 = p
- // Quadrature points: (p=0.5, wt=1.0)
+ // Basis fns: N0 = 0.5*(1-p), N1 = 0.5*(1+p)
+ // Quadrature points: (p=0.0, wt=2.0)
const int cellDim = 1;
const int numCorners = 2;
const int numQuadPts = 1;
const int spaceDim = 1;
const double basis[] = { 0.5, 0.5 };
- const double basisDeriv[] = { -1.0, 1.0 };
- const double quadPtsRef[] = { 0.5 };
- const double quadWts[] = { 1.0 };
- const double jacobianTol = 1.0;
+ const double basisDeriv[] = { -0.5, 0.5 };
+ const double quadPtsRef[] = { 0.0 };
+ const double quadWts[] = { 2.0 };
+ const double jacobianTol = 1.0e-06;
Quadrature1D q;
+ q.jacobianTolerance(jacobianTol);
q.initialize(basis, basisDeriv, quadPtsRef, quadWts,
cellDim, numCorners, numQuadPts, spaceDim);
@@ -57,9 +58,9 @@
const double vertCoords[] = { -0.25, 2.0 };
const int cells[] = { 0, 1 };
const double quadPts[] = { 0.875 };
- const double jacobian[] = { 2.25 };
- const double jacobianInv[] = { 1.0/2.25 };
- const double jacobianDet[] = { 2.25 };
+ const double jacobian[] = { 1.125 };
+ const double jacobianInv[] = { 1.0/1.125 };
+ const double jacobianDet[] = { 1.125 };
// Create mesh with test cell
typedef ALE::Sieve<int, int, int> sieve_type;
@@ -114,7 +115,93 @@
void
pylith::feassemble::TestQuadrature1D::testQuadratic(void)
{ // testQuadratic
- CPPUNIT_ASSERT(false);
+ // Basis fns:
+ // N0 = 0.5*p*(1-p)
+ // N1 = 1-p^2
+ // N2 = 0.5*p*(1+p)
+ // Quadrature points
+ // p=+1/sqrt(3), wt=1
+ // p=-1/sqrt(3), wt=1
+
+
+ const int cellDim = 1;
+ const int numCorners = 3;
+ const int numQuadPts = 2;
+ const int spaceDim = 1;
+ const double basis[] = {
+ 4.55341801e-01, 6.66666667e-01, -1.22008468e-01,
+ -1.22008468e-01, 6.66666667e-01, 4.55341801e-01,
+ };
+ const double basisDeriv[] = {
+ -1.07735027e+00, 1.15470054e+00, -7.73502692e-02,
+ 7.73502692e-02, -1.15470054e+00, 1.07735027e+00,
+ };
+ const double quadPtsRef[] = { -1.0/sqrt(3.0), +1.0/sqrt(3.0) };
+ const double quadWts[] = { 1.0 };
+ const double jacobianTol = -10.0;
+
+ Quadrature1D q;
+ q.jacobianTolerance(jacobianTol);
+ q.initialize(basis, basisDeriv, quadPtsRef, quadWts,
+ cellDim, numCorners, numQuadPts, spaceDim);
+
+ // Test cell: x0=-0.25, x1=2.0, x2=2.0
+ const int numVertices = 3;
+ const int numCells = 1;
+ const double vertCoords[] = { -0.25, 0.875, 2.0 };
+ const int cells[] = { 0, 1, 2 };
+ const double quadPts[] = { 2.25480947e-01, 1.52451905e+00 };
+ const double jacobian[] = { 1.125, 1.125 };
+ const double jacobianInv[] = { 1.0/1.125, 1.0/1.125 };
+ const double jacobianDet[] = { 1.125, 1.125 };
+
+ // Create mesh with test cell
+ typedef ALE::Sieve<int, int, int> sieve_type;
+ typedef ALE::New::Topology<int, sieve_type> topology_type;
+ typedef ALE::New::Section<topology_type, double> section_type;
+ ALE::Obj<ALE::Mesh> mesh = ALE::Mesh(PETSC_COMM_WORLD, cellDim);
+ ALE::Obj<sieve_type> sieve = new sieve_type(mesh->comm());
+ ALE::Obj<topology_type> topology = new topology_type(mesh->comm());
+
+ const bool interpolate = false;
+ ALE::New::SieveBuilder<sieve_type>::buildTopology(sieve, cellDim, numCells,
+ (int*) cells, numVertices, interpolate, numCorners);
+ sieve->stratify();
+ topology->setPatch(0, sieve);
+ topology->stratify();
+ mesh->setTopologyNew(topology);
+ ALE::New::SieveBuilder<sieve_type>::buildCoordinates(
+ mesh->getSection("coordinates"), spaceDim, vertCoords);
+
+ // Check values from _computeGeometry()
+ const ALE::Mesh::topology_type::patch_type patch = 0;
+ const ALE::Obj<topology_type::label_sequence>& elements =
+ topology->heightStratum(patch, 0);
+ const topology_type::label_sequence::iterator e_iter = elements->begin();
+ const ALE::Obj<ALE::Mesh::section_type>& coordinates =
+ mesh->getSection("coordinates");
+ q._computeGeometry(coordinates, *e_iter);
+
+ const double tolerance = 1.0e-06;
+ int size = 1;
+ CPPUNIT_ASSERT(0 != q._quadPts);
+ for (int i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(quadPts[i], q._quadPts[i], tolerance);
+
+ size = 1;
+ CPPUNIT_ASSERT(0 != q._jacobian);
+ for (int i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(jacobian[i], q._jacobian[i], tolerance);
+
+ size = 1;
+ CPPUNIT_ASSERT(0 != q._jacobianInv);
+ for (int i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(jacobianInv[i], q._jacobianInv[i], tolerance);
+
+ size = 1;
+ CPPUNIT_ASSERT(0 != q._jacobianDet);
+ for (int i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(jacobianDet[i], q._jacobianDet[i], tolerance);
} // testQuadratic
// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/quadrature1d.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/quadrature1d.py 2006-09-19 19:50:58 UTC (rev 4570)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/quadrature1d.py 2006-09-20 03:54:07 UTC (rev 4571)
@@ -0,0 +1,51 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+def N0(p):
+ return -0.5*p*(1-p)
+
+
+def N1(p):
+ return (1-p**2)
+
+
+def N2(p):
+ return 0.5*p*(1+p)
+
+
+def N0p(p):
+ return +0.5*p - 0.5*(1-p)
+
+
+def N1p(p):
+ return -2.0*p
+
+
+def N2p(p):
+ return 0.5*p + 0.5*(1+p)
+
+
+
+q = [-1.0/3.0**0.5, 1.0/3.0**0.5]
+#q = [-1.0, 0.0, 1.0]
+print "basis = {"
+for qp in q:
+ print "%16.8e, %16.8e, %16.8e," % (N0(qp), N1(qp), N2(qp))
+print "}"
+
+print "basisDeriv = {"
+for qp in q:
+ print "%16.8e, %16.8e, %16.8e," % (N0p(qp), N1p(qp), N2p(qp))
+print "}"
+
+
+# End of file
More information about the cig-commits
mailing list