[cig-commits] r7235 - in short/3D/PyLith/trunk: . libsrc/feassemble
libsrc/materials
brad at geodynamics.org
brad at geodynamics.org
Thu Jun 14 10:12:12 PDT 2007
Author: brad
Date: 2007-06-14 10:12:12 -0700 (Thu, 14 Jun 2007)
New Revision: 7235
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh
Log:
Fixed bug associated with derivatives of basis functions not being with respect to global coordinates but with respect to local cell coordinates. Derivatives with respect to the global coordinates are computed as part of Quadrature::computeGeometry(). Need to update unit tests to account for fixing bug (unit tests had bug too). Need to add unit test exposing problem.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-06-14 17:10:08 UTC (rev 7234)
+++ short/3D/PyLith/trunk/TODO 2007-06-14 17:12:12 UTC (rev 7235)
@@ -2,6 +2,11 @@
MAIN PRIORITIES (Brad)
======================================================================
+Create unit tests for Quadrature::basisDeriv().
+
+Add unit test for verifying cell orientation doesn't affect
+contribution to Jacobian operator.
+
FaultCohesiveKin::integrateJacobian()
Need to prevent multiple contributions of orientation information for
each vertex, because we want to use ADD_VALUES in updateOperator().
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2007-06-14 17:10:08 UTC (rev 7234)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2007-06-14 17:12:12 UTC (rev 7235)
@@ -50,9 +50,7 @@
_quadWts(q._quadWts),
_basis(q._basis),
_basisDerivRef(q._basisDerivRef),
-#if 0
_basisDeriv(q._basisDeriv),
-#endif
_jacobian(q._jacobian),
_jacobianInv(q._jacobianInv),
_jacobianDet(q._jacobianDet),
@@ -107,27 +105,27 @@
} // if
if (cellDim > 0) {
- int size = numBasis * cellDim;
+ int size = numBasis * cellDim; assert(size > 0);
_vertices.resize(size);
for (int i=0; i < size; ++i)
_vertices[i] = vertices[i];
- size = numBasis * numQuadPts;
+ size = numBasis * numQuadPts; assert(size > 0);
_basis.resize(size);
for (int i=0; i < size; ++i)
_basis[i] = basis[i];
- size = numBasis * numQuadPts * cellDim;
+ size = numQuadPts * numBasis * cellDim; assert(size > 0);
_basisDerivRef.resize(size);
for (int i=0; i < size; ++i)
_basisDerivRef[i] = basisDerivRef[i];
- size = numQuadPts * cellDim;
+ size = numQuadPts * cellDim; assert(size > 0);
_quadPtsRef.resize(size);
for (int i=0; i < size; ++i)
_quadPtsRef[i] = quadPtsRef[i];
- size = numQuadPts;
+ size = numQuadPts; assert(size > 0);
_quadWts.resize(size);
for (int i=0; i < size; ++i)
_quadWts[i] = quadWts[i];
@@ -138,22 +136,20 @@
_spaceDim = spaceDim;
// Allocate for Jacobian and its inverse
- size = numQuadPts * cellDim * spaceDim;
+ size = numQuadPts * cellDim * spaceDim; assert(size > 0);
_jacobian.resize(size);
_jacobianInv.resize(size);
// Allocate for Jacobian determinant
- size = numQuadPts;
+ size = numQuadPts; assert(size > 0);
_jacobianDet.resize(size);
-#if 0
// Allocate for basis derivatives (in global coordinates)
- size = numBasis * numQuadPts * spaceDim;
+ size = numQuadPts * numBasis * spaceDim; assert(size > 0);
_basisDeriv.resize(size);
-#endif
// Allocate for quad pts
- size = numQuadPts*spaceDim;
+ size = numQuadPts*spaceDim; assert(size > 0);
_quadPts.resize(size);
} else {
if (1 != numBasis ||
@@ -209,14 +205,12 @@
size = 1;
_jacobianDet.resize(size);
-#if 0
// Allocate for basis derivatives (in global coordinates)
- size = numBasis * numQuadPts * spaceDim;
+ size = numQuadPts * numBasis * spaceDim; assert(size > 0);
_basisDeriv.resize(size);
-#endif
// Allocate for quad pts
- size = spaceDim;
+ size = spaceDim; assert(size > 0);
_quadPts.resize(size);
} // else
} // initialize
@@ -246,9 +240,7 @@
_jacobian = 0.0;
_jacobianDet = 0.0;
_jacobianInv = 0.0;
-#if 0
_basisDeriv = 0.0;
-#endif
_quadPts = 0.0;
} // _resetGeometry
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2007-06-14 17:10:08 UTC (rev 7234)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2007-06-14 17:12:12 UTC (rev 7235)
@@ -320,7 +320,7 @@
* size = numQuadPts * numBasis * cellDim
* index = iQuadPt*numBasis*cellDim + iBasis*cellDim + iDim
*/
- //double_array _basisDeriv;
+ double_array _basisDeriv;
/** Array of Jacobians evaluated at quadrature points.
*
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc 2007-06-14 17:10:08 UTC (rev 7234)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc 2007-06-14 17:12:12 UTC (rev 7235)
@@ -60,7 +60,7 @@
inline
const pylith::double_array&
pylith::feassemble::Quadrature::basisDeriv(void) const {
- return _basisDerivRef;
+ return _basisDeriv;
}
// Get Jacobians evaluated at quadrature points.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc 2007-06-14 17:10:08 UTC (rev 7234)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc 2007-06-14 17:12:12 UTC (rev 7235)
@@ -62,9 +62,7 @@
_jacobian[0] = 1.0;
_jacobianDet[0] = 1.0;
_jacobianInv[0] = 1.0;
-#if 0
_basisDeriv[0] = _basisDerivRef[0];
-#endif
} // computeGeometry
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc 2007-06-14 17:10:08 UTC (rev 7234)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc 2007-06-14 17:12:12 UTC (rev 7235)
@@ -80,16 +80,17 @@
// Jinv = 1/j00
_jacobianInv[iQuadPt] = 1.0/_jacobianDet[iQuadPt];
-#if 0
+ assert(_numQuadPts*_numBasis*_spaceDim == _basisDeriv.size());
+ assert(_numQuadPts*_numBasis*_cellDim == _basisDerivRef.size());
+ assert(_numQuadPts*_cellDim*_spaceDim == _jacobianInv.size());
+
// Compute derivatives of basis functions with respect to global
// coordinates
// dNi/dx = dNi/dp dp/dx + dNi/dq dq/dx + dNi/dr dr/dx
for (int iBasis=0; iBasis < _numBasis; ++iBasis)
- for (int iDim=0; iDim < _spaceDim; ++iDim)
- _basisDeriv[iQuadPt*_numBasis*_spaceDim+iBasis*_spaceDim+iDim] +=
- _basisDerivRef[iQuadPt*_numBasis*+iBasis] *
- _jacobianInv[iQuadPt*_spaceDim+iDim];
-#endif
+ _basisDeriv[iQuadPt*_numBasis+iBasis] +=
+ _basisDerivRef[iQuadPt*_numBasis+iBasis] *
+ _jacobianInv[iQuadPt];
} // for
} // computeGeometry
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc 2007-06-14 17:10:08 UTC (rev 7234)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc 2007-06-14 17:12:12 UTC (rev 7235)
@@ -94,7 +94,6 @@
_jacobianInv[iQuadPt*_spaceDim+iDim] =
1.0 / _jacobian[iQuadPt*_spaceDim+iDim];
-#if 0
// Compute derivatives of basis functions with respect to global
// coordinates
// dNi/dx = dNi/dp dp/dx + dNi/dq dq/dx + dNi/dr dr/dx
@@ -104,7 +103,6 @@
_basisDeriv[iQuadPt*_numBasis*_spaceDim+iBasis*_spaceDim+iDim] +=
_basisDerivRef[iQuadPt*_numBasis*_cellDim+iBasis*_cellDim+jDim] *
_jacobianInv[iQuadPt*_cellDim*_spaceDim+jDim*_spaceDim+iDim];
-#endif
} // for
} // computeGeometry
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc 2007-06-14 17:10:08 UTC (rev 7234)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc 2007-06-14 17:12:12 UTC (rev 7235)
@@ -97,7 +97,6 @@
_jacobianInv[iQuadPt*_spaceDim+iDim] =
1.0 / _jacobian[iQuadPt*_spaceDim+iDim];
-#if 0
// Compute derivatives of basis functions with respect to global
// coordinates
// dNi/dx = dNi/dp dp/dx + dNi/dq dq/dx + dNi/dr dr/dx
@@ -107,7 +106,6 @@
_basisDeriv[iQuadPt*_numBasis*_spaceDim+iBasis*_spaceDim+iDim] +=
_basisDerivRef[iQuadPt*_numBasis*_cellDim+iBasis*_cellDim+jDim] *
_jacobianInv[iQuadPt*_cellDim*_spaceDim+jDim*_spaceDim+iDim];
-#endif
} // for
} // computeGeometry
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc 2007-06-14 17:10:08 UTC (rev 7234)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc 2007-06-14 17:12:12 UTC (rev 7235)
@@ -104,7 +104,6 @@
_jacobianInv[i11] = _jacobian[i00] / det;
PetscLogFlops(_numBasis*2 + _numBasis*_cellDim*_spaceDim*2 + 21);
-#if 0
// Compute derivatives of basis functions with respect to global
// coordinates
// dNi/dx = dNi/dp dp/dx + dNi/dq dq/dx + dNi/dr dr/dx
@@ -114,7 +113,6 @@
_basisDeriv[iQuadPt*_numBasis*_spaceDim+iBasis*_spaceDim+iDim] +=
_basisDerivRef[iQuadPt*_numBasis*_cellDim+iBasis*_cellDim+jDim] *
_jacobianInv[iQuadPt*_cellDim*_spaceDim+jDim*_spaceDim+iDim];
-#endif
} // for
} // computeGeometry
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc 2007-06-14 17:10:08 UTC (rev 7234)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc 2007-06-14 17:12:12 UTC (rev 7235)
@@ -182,7 +182,6 @@
} else
throw std::runtime_error("Could not invert Jacobian.");
-#if 0
// Compute derivatives of basis functions with respect to global
// coordinates
// dNi/dx = dNi/dp dp/dx + dNi/dq dq/dx + dNi/dr dr/dx
@@ -192,7 +191,6 @@
_basisDeriv[iQuadPt*_numBasis*_spaceDim+iBasis*_spaceDim+iDim] +=
_basisDerivRef[iQuadPt*_numBasis*_cellDim + iBasis*_cellDim+jDim] *
_jacobianInv[iQuadPt*_cellDim*_spaceDim+jDim*_spaceDim+iDim];
-#endif
} // for
} // computeGeometry
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc 2007-06-14 17:10:08 UTC (rev 7234)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc 2007-06-14 17:12:12 UTC (rev 7235)
@@ -133,7 +133,6 @@
_jacobianInv[i22] = (_jacobian[i00]*_jacobian[i11] -
_jacobian[i01]*_jacobian[i10]) / det;
-#if 0
// Compute derivatives of basis functions with respect to global
// coordinates
// dNi/dx = dNi/dp dp/dx + dNi/dq dq/dx + dNi/dr dr/dx
@@ -143,7 +142,6 @@
_basisDeriv[iQuadPt*_numBasis*_spaceDim+iBasis*_spaceDim+iDim] +=
_basisDerivRef[iQuadPt*_numBasis*_cellDim+iBasis*_cellDim+jDim] *
_jacobianInv[iQuadPt*_cellDim*_spaceDim+jDim*_spaceDim+iDim];
-#endif
} // for
} // computeGeometry
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh 2007-06-14 17:10:08 UTC (rev 7234)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh 2007-06-14 17:12:12 UTC (rev 7235)
@@ -22,9 +22,7 @@
#include "Material.hh" // ISA Material
-#include "pylith/utils/arrayfwd.hh" // USES double_array
#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
-#include <vector> // USES std::vector
/// Namespace for pylith package
namespace pylith {
More information about the cig-commits
mailing list