[cig-commits] r21536 - in short/3D/PyLith/trunk: libsrc/pylith/feassemble unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Thu Mar 14 11:47:29 PDT 2013
Author: brad
Date: 2013-03-14 11:47:28 -0700 (Thu, 14 Mar 2013)
New Revision: 21536
Modified:
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.icc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature0D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature0D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din2D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din2D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din3D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din3D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2Din3D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2Din3D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature3D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature3D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureEngine.hh
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
Log:
Added compute geometry with C array to remove copy.
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.hh 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.hh 2013-03-14 18:47:28 UTC (rev 21536)
@@ -123,7 +123,7 @@
* @param cells Finite-element cells for geometry.
*/
void computeGeometry(const mesh_type& mesh,
- const ALE::Obj<typename mesh_type::SieveMesh::label_sequence>& cells);
+ const ALE::Obj<typename mesh_type::SieveMesh::label_sequence>& cells);
/** Compute geometric quantities for each cell.
*
@@ -150,6 +150,16 @@
void computeGeometry(const scalar_array& coordinatesCell,
const int cell);
+ /** Compute geometric quantities for a cell at quadrature points.
+ *
+ * @param coordinatesCell Array of coordinates of cell's vertices.
+ * @param coordinatesSize Size of coordinates array.
+ * @param cell Finite-element cell
+ */
+ void computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell);
+
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.icc 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.icc 2013-03-14 18:47:28 UTC (rev 21536)
@@ -89,10 +89,23 @@
void
pylith::feassemble::Quadrature<mesh_type>::computeGeometry(const scalar_array& coordinatesCell,
const int cell) {
- assert(0 != _engine);
+ assert(_engine);
_engine->computeGeometry(coordinatesCell, cell);
}
+// Compute geometric quantities for a cell at quadrature points.
+template<typename mesh_type>
+void
+pylith::feassemble::Quadrature<mesh_type>::computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell)
+{ // computeGeometry
+ assert(_engine);
+ _engine->computeGeometry(coordinatesCell, coordinatesSize, cell);
+} // computeGeometry
+
+
+
#endif
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature0D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature0D.cc 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature0D.cc 2013-03-14 18:47:28 UTC (rev 21536)
@@ -79,4 +79,39 @@
} // computeGeometry
+// ----------------------------------------------------------------------
+// Compute geometric quantities for a cell at quadrature points.
+void
+pylith::feassemble::Quadrature0D::computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell)
+{ // computeGeometry
+ assert(coordinatesCell);
+
+ const int cellDim = 0;
+ const int spaceDim = 1;
+ const int numQuadPts = 1;
+ const int numBasis = 1;
+
+ assert(_quadRefCell.cellDim() == cellDim);
+ assert(_quadRefCell.spaceDim() == spaceDim);
+ assert(_quadRefCell.numQuadPts() == numQuadPts);
+ assert(_quadRefCell.numBasis() == numBasis);
+ assert(coordinatesSize == numBasis*spaceDim);
+
+ const scalar_array& basisDerivRef = _quadRefCell.basisDerivRef();
+
+ zero();
+
+ _quadPts = scalar_array(coordinatesCell, coordinatesSize);
+
+ _jacobian[0] = 1.0;
+ _jacobianDet[0] = 1.0;
+ _jacobianInv[0] = 1.0;
+ _basisDeriv[0] = basisDerivRef[0];
+
+ PetscLogFlops(0);
+} // computeGeometry
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature0D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature0D.hh 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature0D.hh 2013-03-14 18:47:28 UTC (rev 21536)
@@ -61,6 +61,16 @@
void computeGeometry(const scalar_array& coordinatesCell,
const int cell);
+ /** Compute geometric quantities for a cell at quadrature points.
+ *
+ * @param coordinatesCell Array of coordinates of cell's vertices.
+ * @param coordinatesSize Size of coordinates array.
+ * @param cell Finite-element cell
+ */
+ void computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell);
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1D.cc 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1D.cc 2013-03-14 18:47:28 UTC (rev 21536)
@@ -124,4 +124,82 @@
} // computeGeometry
+// ----------------------------------------------------------------------
+// Compute geometric quantities for a cell at quadrature points.
+void
+pylith::feassemble::Quadrature1D::computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell)
+{ // computeGeometry
+ assert(coordinatesCell);
+
+ const int cellDim = 1;
+ const int spaceDim = 1;
+
+ const int numQuadPts = _quadRefCell.numQuadPts();
+ const int numBasis = _quadRefCell.numBasis();
+
+ const scalar_array& basis = _quadRefCell.basis();
+ const scalar_array& quadPtsRef = _quadRefCell.quadPtsRef();
+ const scalar_array& basisDerivRef = _quadRefCell.basisDerivRef();
+ const CellGeometry& geometry = _quadRefCell.refGeometry();
+
+ assert(_quadRefCell.cellDim() == cellDim);
+ assert(_quadRefCell.spaceDim() == spaceDim);
+ assert(numBasis*spaceDim == coordinatesSize);
+
+ zero();
+
+ // Loop over quadrature points
+ for (int iQuadPt=0; iQuadPt < numQuadPts; ++iQuadPt) {
+ const int iQ = iQuadPt*numBasis;
+
+ // Compute coordinates of quadrature point in cell
+#if defined(ISOPARAMETRIC)
+ // x = sum[i=0,n-1] (Ni * xi)
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ _quadPts[iQuadPt] += basis[iQ+iBasis]*coordinatesCell[iBasis];
+#else
+ geometry.ptsRefToGlobal(&_quadPts[iQuadPt], &quadPtsRef[iQuadPt],
+ coordinatesCell, spaceDim, 1);
+#endif
+
+#if defined(ISOPARAMETRIC)
+ // Compute Jacobian at quadrature point
+ // J = dx/dp = sum[i=0,n-1] (dNi/dp * xi)
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ _jacobian[iQuadPt] += basisDerivRef[iQ+iBasis] * coordinatesCell[iBasis];
+
+ // Compute determinant of Jacobian at quadrature point
+ // |J| = j00
+ const PylithScalar det = _jacobian[iQuadPt];
+ _checkJacobianDet(det, cell);
+ _jacobianDet[iQuadPt] = _jacobian[iQuadPt];
+#else
+ // Compute Jacobian and determinant of Jacobian at quadrature point
+ geometry.jacobian(&_jacobian[iQuadPt], &_jacobianDet[iQuadPt],
+ coordinatesCell, &quadPtsRef[iQuadPt], spaceDim, 1);
+ _checkJacobianDet(_jacobianDet[iQuadPt], cell);
+#endif
+
+ // Compute inverse of Jacobian at quadrature point
+ // Jinv = 1/j00
+ _jacobianInv[iQuadPt] = 1.0 / _jacobianDet[iQuadPt];
+
+ 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)
+ _basisDeriv[iQ+iBasis] +=
+ basisDerivRef[iQ+iBasis] * _jacobianInv[iQuadPt];
+ } // for
+
+ PetscLogFlops(numQuadPts * (1 + numBasis * 4));
+} // computeGeometry
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1D.hh 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1D.hh 2013-03-14 18:47:28 UTC (rev 21536)
@@ -57,6 +57,16 @@
void computeGeometry(const scalar_array& coordinatesCell,
const int cell);
+ /** Compute geometric quantities for a cell at quadrature points.
+ *
+ * @param coordinatesCell Array of coordinates of cell's vertices.
+ * @param coordinatesSize Size of coordinates array.
+ * @param cell Finite-element cell
+ */
+ void computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell);
+
// PRIVATE METHODS //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din2D.cc 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din2D.cc 2013-03-14 18:47:28 UTC (rev 21536)
@@ -149,4 +149,107 @@
} // computeGeometry
+// ----------------------------------------------------------------------
+// Compute geometric quantities for a cell at quadrature points.
+void
+pylith::feassemble::Quadrature1Din2D::computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell)
+{ // computeGeometry
+ assert(coordinatesCell);
+
+ const int cellDim = 1;
+ const int spaceDim = 2;
+
+ const int numQuadPts = _quadRefCell.numQuadPts();
+ const int numBasis = _quadRefCell.numBasis();
+
+ const scalar_array& basis = _quadRefCell.basis();
+ const scalar_array& quadPtsRef = _quadRefCell.quadPtsRef();
+ const scalar_array& basisDerivRef = _quadRefCell.basisDerivRef();
+ const CellGeometry& geometry = _quadRefCell.refGeometry();
+
+ assert(_quadRefCell.cellDim() == cellDim);
+ assert(_quadRefCell.spaceDim() == spaceDim);
+ assert(numBasis*spaceDim == coordinatesSize);
+
+ zero();
+
+ // Loop over quadrature points
+ for (int iQuadPt=0; iQuadPt < numQuadPts; ++iQuadPt) {
+ const int iQ = iQuadPt*numBasis;
+
+ // Compute coordinates of quadrature point in cell
+#if defined(ISOPARAMETRIC)
+ // x = sum[i=0,n-1] (Ni * xi)
+ // y = sum[i=0,n-1] (Ni * yi)
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const PylithScalar valueBasis = basis[iQ+iBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _quadPts[iQuadPt*spaceDim+iDim] +=
+ valueBasis * coordinatesCell[iBasis*spaceDim+iDim];
+ } // for
+#else
+ geometry.ptsRefToGlobal(&_quadPts[iQuadPt*spaceDim],
+ &quadPtsRef[iQuadPt*cellDim],
+ coordinatesCell, spaceDim, 1);
+#endif
+
+#if defined(ISOPARAMETRIC)
+ // Compute Jacobian at quadrature point
+ // J = [dx/dp
+ // dy/dp]
+ // dx/dp = sum[i=0,n-1] (dNi/dp * xi)
+ // dy/dp = sum[i=0,n-1] (dNi/dp * yi)
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const PylithScalar deriv = basisDerivRef[iQ+iBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _jacobian[iQuadPt*spaceDim+iDim] +=
+ deriv * coordinatesCell[iBasis*spaceDim+iDim];
+ } // for
+
+ // Compute determinant of Jacobian at quadrature point
+ // |J| = sqrt(transpose(J) J)
+ PylithScalar det = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ det += _jacobian[iQuadPt*spaceDim+iDim] *
+ _jacobian[iQuadPt*spaceDim+iDim];
+ det = sqrt(det);
+ _checkJacobianDet(det, cell);
+ _jacobianDet[iQuadPt] = det;
+#else
+ // Compute Jacobian and determinant of Jacobian at quadrature point
+ geometry.jacobian(&_jacobian[iQuadPt*cellDim*spaceDim],
+ &_jacobianDet[iQuadPt],
+ coordinatesCell, &quadPtsRef[iQuadPt*cellDim],
+ spaceDim, 1);
+ _checkJacobianDet(_jacobianDet[iQuadPt], cell);
+#endif
+
+ // Compute inverse of Jacobian at quadrature point
+ // Jinv = 1.0/[J]
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _jacobianInv[iQuadPt*spaceDim+iDim] =
+ 1.0 / _jacobian[iQuadPt*spaceDim+iDim];
+
+ // Compute derivatives of basis functions with respect to global
+ // coordinates
+ // dNi/dx = dNi/dp dp/dx + dNi/dq dq/dx + dNi/dr dr/dx
+ const int iJ = iQuadPt*cellDim*spaceDim;
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const int iD = iQuadPt*numBasis*spaceDim + iBasis*spaceDim;
+ const int iDR = iQuadPt*numBasis*cellDim + iBasis*cellDim;
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int jDim=0; jDim < cellDim; ++jDim)
+ _basisDeriv[iD+iDim] +=
+ basisDerivRef[iDR+jDim] * _jacobianInv[iJ+jDim*spaceDim+iDim];
+ } // for
+ } // for
+
+ PetscLogFlops(numQuadPts * (1 + numBasis*spaceDim*2 +
+ spaceDim*1 +
+ numBasis*spaceDim*cellDim*2));
+} // computeGeometry
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din2D.hh 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din2D.hh 2013-03-14 18:47:28 UTC (rev 21536)
@@ -57,6 +57,16 @@
void computeGeometry(const scalar_array& coordinatesCell,
const int cell);
+ /** Compute geometric quantities for a cell at quadrature points.
+ *
+ * @param coordinatesCell Array of coordinates of cell's vertices.
+ * @param coordinatesSize Size of coordinates array.
+ * @param cell Finite-element cell
+ */
+ void computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell);
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din3D.cc 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din3D.cc 2013-03-14 18:47:28 UTC (rev 21536)
@@ -153,4 +153,111 @@
} // computeGeometry
+// ----------------------------------------------------------------------
+// Compute geometric quantities for a cell at quadrature points.
+void
+pylith::feassemble::Quadrature1Din3D::computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell)
+{ // computeGeometry
+ assert(coordinatesCell);
+
+ const int cellDim = 1;
+ const int spaceDim = 3;
+
+ const int numQuadPts = _quadRefCell.numQuadPts();
+ const int numBasis = _quadRefCell.numBasis();
+
+ const scalar_array& basis = _quadRefCell.basis();
+ const scalar_array& quadPtsRef = _quadRefCell.quadPtsRef();
+ const scalar_array& basisDerivRef = _quadRefCell.basisDerivRef();
+ const CellGeometry& geometry = _quadRefCell.refGeometry();
+
+ assert(_quadRefCell.cellDim() == cellDim);
+ assert(_quadRefCell.spaceDim() == spaceDim);
+ assert(numBasis*spaceDim == coordinatesSize);
+
+ zero();
+
+ // Loop over quadrature points
+ for (int iQuadPt=0; iQuadPt < numQuadPts; ++iQuadPt) {
+ const int iQ = iQuadPt*numBasis;
+
+ // Compute coordinates of quadrature point in cell
+#if defined(ISOPARAMETRIC)
+ // x = sum[i=0,n-1] (Ni * xi)
+ // y = sum[i=0,n-1] (Ni * yi)
+ // z = sum[i=0,n-1] (Ni * zi)
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const PylithScalar valueBasis = basis[iQ+iBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _quadPts[iQuadPt*spaceDim+iDim] +=
+ valueBasis * coordinatesCell[iBasis*spaceDim+iDim];
+ } // for
+#else
+ geometry.ptsRefToGlobal(&_quadPts[iQuadPt*spaceDim],
+ &quadPtsRef[iQuadPt*cellDim],
+ coordinatesCell, spaceDim, 1);
+#endif
+
+#if defined(ISOPARAMETRIC)
+ // Compute Jacobian at quadrature point
+ // J = [dx/dp
+ // dy/dp
+ // dz/dp]
+ // dx/dp = sum[i=0,n-1] (dNi/dp * xi)
+ // dy/dp = sum[i=0,n-1] (dNi/dp * yi)
+ // dz/dp = sum[i=0,n-1] (dNi/dp * zi)
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const PylithScalar deriv = basisDerivRef[iQ+iBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _jacobian[iQuadPt*spaceDim+iDim] +=
+ deriv * coordinatesCell[iBasis*spaceDim+iDim];
+ } // for
+
+ // Compute determinant of Jacobian at quadrature point
+ // |J| = sqrt(transpose(J) J)
+ PylithScalar det = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ det += _jacobian[iQuadPt*spaceDim+iDim] *
+ _jacobian[iQuadPt*spaceDim+iDim];
+ det = sqrt(det);
+ _checkJacobianDet(det, cell);
+ _jacobianDet[iQuadPt] = det;
+#else
+ // Compute Jacobian and determinant of Jacobian at quadrature point
+ geometry.jacobian(&_jacobian[iQuadPt*cellDim*spaceDim],
+ &_jacobianDet[iQuadPt],
+ coordinatesCell, &quadPtsRef[iQuadPt*cellDim],
+ spaceDim, 1);
+ _checkJacobianDet(_jacobianDet[iQuadPt], cell);
+#endif
+
+ // Compute inverse of Jacobian at quadrature point
+ // Jinv = 1.0/[J]
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _jacobianInv[iQuadPt*spaceDim+iDim] =
+ 1.0 / _jacobian[iQuadPt*spaceDim+iDim];
+
+ // Compute derivatives of basis functions with respect to global
+ // coordinates
+ // dNi/dx = dNi/dp dp/dx + dNi/dq dq/dx + dNi/dr dr/dx
+ const int iJ = iQuadPt*cellDim*spaceDim;
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const int iD = iQuadPt*numBasis*spaceDim + iBasis*spaceDim;
+ const int iDR = iQuadPt*numBasis*cellDim + iBasis*cellDim;
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int jDim=0; jDim < cellDim; ++jDim)
+ _basisDeriv[iD+iDim] +=
+ basisDerivRef[iDR+jDim] * _jacobianInv[iJ+jDim*spaceDim+iDim];
+ } // for
+ } // for
+
+ PetscLogFlops(numQuadPts * (1 + numBasis*spaceDim*2 +
+ spaceDim*1 +
+ numBasis*spaceDim*cellDim*2));
+
+} // computeGeometry
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din3D.hh 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din3D.hh 2013-03-14 18:47:28 UTC (rev 21536)
@@ -58,6 +58,16 @@
void computeGeometry(const scalar_array& coordinatesCell,
const int cell);
+ /** Compute geometric quantities for a cell at quadrature points.
+ *
+ * @param coordinatesCell Array of coordinates of cell's vertices.
+ * @param coordinatesSize Size of coordinates array.
+ * @param cell Finite-element cell
+ */
+ void computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell);
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2D.cc 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2D.cc 2013-03-14 18:47:28 UTC (rev 21536)
@@ -161,4 +161,120 @@
} // computeGeometry
+// ----------------------------------------------------------------------
+// Compute geometric quantities for a cell at quadrature points.
+void
+pylith::feassemble::Quadrature2D::computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell)
+{ // computeGeometry
+ assert(coordinatesCell);
+
+ const int spaceDim = 2;
+ const int cellDim = 2;
+
+ const int numQuadPts = _quadRefCell.numQuadPts();
+ const int numBasis = _quadRefCell.numBasis();
+
+ const scalar_array& basis = _quadRefCell.basis();
+ const scalar_array& quadPtsRef = _quadRefCell.quadPtsRef();
+ const scalar_array& basisDerivRef = _quadRefCell.basisDerivRef();
+ const CellGeometry& geometry = _quadRefCell.refGeometry();
+
+ assert(_quadRefCell.cellDim() == cellDim);
+ assert(_quadRefCell.spaceDim() == spaceDim);
+ assert(numBasis*spaceDim == coordinatesSize);
+
+ zero();
+
+ // Loop over quadrature points
+ for (int iQuadPt=0; iQuadPt < numQuadPts; ++iQuadPt) {
+
+ // Compute coordinates of quadrature point in cell
+#if defined(ISOPARAMETRIC)
+ // x = sum[i=0,n-1] (Ni * xi)
+ // y = sum[i=0,n-1] (Ni * yi)
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const PylithScalar valueBasis = basis[iQuadPt*numBasis+iBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _quadPts[iQuadPt*spaceDim+iDim] +=
+ valueBasis * coordinatesCell[iBasis*spaceDim+iDim];
+ } // for
+#else
+ geometry.ptsRefToGlobal(&_quadPts[iQuadPt*spaceDim],
+ &quadPtsRef[iQuadPt*cellDim],
+ coordinatesCell, spaceDim, 1);
+#endif
+
+#if defined(ISOPARAMETRIC)
+ // Compute Jacobian at quadrature point
+ // J = [dx/dp dx/dq]
+ // [dy/dp dy/dq]
+ // dx/dp = sum[i=0,n-1] (dNi/dp * xi)
+ // dx/dq = sum[i=0,n-1] (dNi/dq * xi)
+ // dy/dp = sum[i=0,n-1] (dNi/dp * yi)
+ // dy/dq = sum[i=0,n-1] (dNi/dq * yi)
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ for (int iCol=0; iCol < cellDim; ++iCol) {
+ const PylithScalar deriv =
+ basisDerivRef[iQuadPt*numBasis*spaceDim+iBasis*cellDim+iCol];
+ for (int iRow=0; iRow < spaceDim; ++iRow)
+ _jacobian[iQuadPt*cellDim*spaceDim+iRow*cellDim+iCol] +=
+ deriv * coordinatesCell[iBasis*spaceDim+iRow];
+ } // for
+
+ // Compute determinant of Jacobian at quadrature point
+ // |J| = j00*j11-j01*j10
+ const int iJ = iQuadPt*cellDim*spaceDim;
+ const int i00 = iJ + 0*spaceDim + 0;
+ const int i01 = iJ + 0*spaceDim + 1;
+ const int i10 = iJ + 1*spaceDim + 0;
+ const int i11 = iJ + 1*spaceDim + 1;
+ const PylithScalar det =
+ _jacobian[i00]*_jacobian[i11] -
+ _jacobian[i01]*_jacobian[i10];
+ _checkJacobianDet(det, cell);
+ _jacobianDet[iQuadPt] = det;
+#else
+ // Compute Jacobian and determinant of Jacobian at quadrature point
+ const int iJ = iQuadPt*cellDim*spaceDim;
+ const int i00 = iJ + 0*spaceDim + 0;
+ const int i01 = iJ + 0*spaceDim + 1;
+ const int i10 = iJ + 1*spaceDim + 0;
+ const int i11 = iJ + 1*spaceDim + 1;
+ geometry.jacobian(&_jacobian[iQuadPt*cellDim*spaceDim],
+ &_jacobianDet[iQuadPt],
+ coordinatesCell, &quadPtsRef[iQuadPt*cellDim],
+ spaceDim, 1);
+ _checkJacobianDet(_jacobianDet[iQuadPt], cell);
+ const PylithScalar det = _jacobianDet[iQuadPt];
+#endif
+
+ // Compute inverse of Jacobian at quadrature point
+ // Jinv = 1/det*[ j11 -j01]
+ // [-j10 j00]
+ _jacobianInv[i00] = _jacobian[i11] / det;
+ _jacobianInv[i01] = -_jacobian[i01] / det;
+ _jacobianInv[i10] = -_jacobian[i10] / det;
+ _jacobianInv[i11] = _jacobian[i00] / det;
+
+ // 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) {
+ const int iD = iQuadPt*numBasis*spaceDim + iBasis*spaceDim;
+ const int iDR = iQuadPt*numBasis*cellDim + iBasis*cellDim;
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int jDim=0; jDim < cellDim; ++jDim)
+ _basisDeriv[iD+iDim] += basisDerivRef[iDR+jDim] *
+ _jacobianInv[iJ+jDim*spaceDim+iDim];
+ } // for
+ } // for
+
+ PetscLogFlops(numQuadPts*(4 +
+ numBasis*spaceDim*2 +
+ numBasis*spaceDim*cellDim*2));
+} // computeGeometry
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2D.hh 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2D.hh 2013-03-14 18:47:28 UTC (rev 21536)
@@ -57,6 +57,16 @@
void computeGeometry(const scalar_array& coordinatesCell,
const int cell);
+ /** Compute geometric quantities for a cell at quadrature points.
+ *
+ * @param coordinatesCell Array of coordinates of cell's vertices.
+ * @param coordinatesSize Size of coordinates array.
+ * @param cell Finite-element cell
+ */
+ void computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell);
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2Din3D.cc 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2Din3D.cc 2013-03-14 18:47:28 UTC (rev 21536)
@@ -241,4 +241,198 @@
} // computeGeometry
+// ----------------------------------------------------------------------
+// Compute geometric quantities for a cell at quadrature points.
+void
+pylith::feassemble::Quadrature2Din3D::computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell)
+{ // computeGeometry
+ assert(coordinatesCell);
+
+ const int cellDim = 2;
+ const int spaceDim = 3;
+
+ const int numQuadPts = _quadRefCell.numQuadPts();
+ const int numBasis = _quadRefCell.numBasis();
+
+ const scalar_array& basis = _quadRefCell.basis();
+ const scalar_array& quadPtsRef = _quadRefCell.quadPtsRef();
+ const scalar_array& basisDerivRef = _quadRefCell.basisDerivRef();
+ const CellGeometry& geometry = _quadRefCell.refGeometry();
+ const PylithScalar minJacobian = _quadRefCell.minJacobian();
+
+ assert(_quadRefCell.cellDim() == cellDim);
+ assert(_quadRefCell.spaceDim() == spaceDim);
+ assert(numBasis*spaceDim == coordinatesSize);
+
+ zero();
+
+ // Loop over quadrature points
+ for (int iQuadPt=0; iQuadPt < numQuadPts; ++iQuadPt) {
+
+ // Compute coordinates of quadrature point in cell
+#if defined(ISOPARAMETRIC)
+ // x = sum[i=0,n-1] (Ni * xi)
+ // y = sum[i=0,n-1] (Ni * yi)
+ // z = sum[i=0,n-1] (Ni * zi)
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const PylithScalar valueBasis = basis[iQuadPt*numBasis+iBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _quadPts[iQuadPt*spaceDim+iDim] +=
+ valueBasis * coordinatesCell[iBasis*spaceDim+iDim];
+ } // for
+#else
+ geometry.ptsRefToGlobal(&_quadPts[iQuadPt*spaceDim],
+ &quadPtsRef[iQuadPt*cellDim],
+ coordinatesCell, spaceDim, 1);
+#endif
+
+#if defined(ISOPARAMETRIC)
+ // Compute Jacobian at quadrature point
+ // J = [dx/dp dx/dq]
+ // [dy/dp dy/dq]
+ // [dz/dp dz/dq]
+ // dx/dp = sum[i=0,n-1] (dNi/dp * xi)
+ // dx/dq = sum[i=0,n-1] (dNi/dq * xi)
+ // dy/dp = sum[i=0,n-1] (dNi/dp * yi)
+ // dy/dq = sum[i=0,n-1] (dNi/dq * yi)
+ // dz/dp = sum[i=0,n-1] (dNi/dp * zi)
+ // dz/dq = sum[i=0,n-1] (dNi/dq * zi)
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ for (int iCol=0; iCol < cellDim; ++iCol) {
+ const PylithScalar deriv =
+ basisDerivRef[iQuadPt*numBasis*cellDim+iBasis*cellDim+iCol];
+ for (int iRow=0; iRow < spaceDim; ++iRow)
+ _jacobian[iQuadPt*cellDim*spaceDim+iRow*cellDim+iCol] +=
+ deriv * coordinatesCell[iBasis*+spaceDim+iRow];
+ } // for
+
+ // Compute determinant of Jacobian at quadrature point
+ // |J| = sqrt(transpose(J) J)
+ const int iJ = iQuadPt*cellDim*spaceDim;
+ const int i00 = iJ + 0*cellDim + 0;
+ const int i01 = iJ + 0*cellDim + 1;
+ const int i10 = iJ + 1*cellDim + 0;
+ const int i11 = iJ + 1*cellDim + 1;
+ const int i20 = iJ + 2*cellDim + 0;
+ const int i21 = iJ + 2*cellDim + 1;
+ // JJ = transpose(J) J
+ const PylithScalar jj00 =
+ _jacobian[i00]*_jacobian[i00] +
+ _jacobian[i10]*_jacobian[i10] +
+ _jacobian[i20]*_jacobian[i20];
+ const PylithScalar jj10 =
+ _jacobian[i00]*_jacobian[i01] +
+ _jacobian[i10]*_jacobian[i11] +
+ _jacobian[i20]*_jacobian[i21];
+ const PylithScalar jj01 = jj10;
+ const PylithScalar jj11 =
+ _jacobian[i01]*_jacobian[i01] +
+ _jacobian[i11]*_jacobian[i11] +
+ _jacobian[i21]*_jacobian[i21];
+ const PylithScalar det = sqrt(jj00*jj11 - jj01*jj10);
+ _checkJacobianDet(det, cell);
+ _jacobianDet[iQuadPt] = det;
+#else
+ // Compute Jacobian and determinant of Jacobian at quadrature point
+ const int iJ = iQuadPt*cellDim*spaceDim;
+ const int i00 = iJ + 0*cellDim + 0;
+ const int i01 = iJ + 0*cellDim + 1;
+ const int i10 = iJ + 1*cellDim + 0;
+ const int i11 = iJ + 1*cellDim + 1;
+ const int i20 = iJ + 2*cellDim + 0;
+ const int i21 = iJ + 2*cellDim + 1;
+ geometry.jacobian(&_jacobian[iQuadPt*cellDim*spaceDim],
+ &_jacobianDet[iQuadPt],
+ coordinatesCell, &quadPtsRef[iQuadPt*cellDim],
+ spaceDim, 1);
+ _checkJacobianDet(_jacobianDet[iQuadPt], cell);
+#endif
+
+ // Compute inverse of Jacobian at quadrature point
+ const PylithScalar d01 =
+ _jacobian[i00]*_jacobian[i11] -
+ _jacobian[i10]*_jacobian[i01];
+ const PylithScalar d12 =
+ _jacobian[i10]*_jacobian[i21] -
+ _jacobian[i20]*_jacobian[i11];
+ const PylithScalar d02 =
+ _jacobian[i00]*_jacobian[i21] -
+ _jacobian[i20]*_jacobian[i01];
+ if (fabs(d01) > minJacobian) {
+ // Jinv00 = 1/d01 * J11
+ // Jinv01 = 1/d01 * -J01
+ // Jinv10 = 1/d01 * -J10
+ // Jinv11 = 1/d01 * J00
+ _jacobianInv[iJ+0] = _jacobian[i11] / d01; // Jinv00
+ _jacobianInv[iJ+1] = -_jacobian[i01] / d01; // Jinv01
+ _jacobianInv[iJ+3] = -_jacobian[i10] / d01; // Jinv10
+ _jacobianInv[iJ+4] = _jacobian[i00] / d01; // Jinv11
+ if (fabs(d12) > minJacobian) {
+ // Jinv02 = 1/d12 -J11
+ // Jinv12 = 1/d12 J10
+ _jacobianInv[iJ+2] = -_jacobian[i11] / d12; // Jinv02
+ _jacobianInv[iJ+5] = _jacobian[i10] / d12; // Jinv12
+
+ } else if (fabs(d02) > minJacobian) {
+ // Jinv02 = 1/d02 -J01
+ // Jinv12 = 1/d02 J00
+ _jacobianInv[iJ+2] = -_jacobian[i01] / d02; // Jinv02
+ _jacobianInv[iJ+5] = _jacobian[i00] / d02; // Jinv12
+ } else {
+ _jacobianInv[iJ+2] = 0.0; // Jinv02
+ _jacobianInv[iJ+5] = 0.0; // Jinv12
+ } // if/else
+ } else if (fabs(d02) > minJacobian) {
+ // Jinv00 = 1/d02 * J21
+ // Jinv02 = 1/d02 * -J01
+ // Jinv10 = 1/d02 * -J20
+ // Jinv12 = 1/d02 * J00
+ _jacobianInv[iJ+0] = _jacobian[i21] / d02; // Jinv00
+ _jacobianInv[iJ+2] = -_jacobian[i01] / d02; // Jinv02
+ _jacobianInv[iJ+3] = -_jacobian[i20] / d02; // Jinv10
+ _jacobianInv[iJ+5] = _jacobian[i00] / d02; // Jinv12
+ if (fabs(d12) > minJacobian) {
+ // Jinv01 = 1/d12 J21
+ // Jinv11 = 1/d12 -J20
+ _jacobianInv[iJ+1] = -_jacobian[i21] / d12; // Jinv01
+ _jacobianInv[iJ+4] = _jacobian[i20] / d12; // Jinv11
+ } else {
+ _jacobianInv[iJ+1] = 0.0; // Jinv01
+ _jacobianInv[iJ+4] = 0.0; // Jinv11
+ } // if/else
+ } else if (fabs(d12) > minJacobian) {
+ _jacobianInv[iJ+0] = 0.0; // Jinv00
+ _jacobianInv[iJ+3] = 0.0; // Jinv10
+ // Jinv01 = 1/d12 J21
+ // Jinv02 = 1/d12 -J11
+ // Jinv11 = 1/d12 -J20
+ // Jinv12 = 1/d12 J10
+ _jacobianInv[iJ+1] = _jacobian[i21] / d12; // Jinv01
+ _jacobianInv[iJ+2] = -_jacobian[i11] / d12; // Jinv02
+ _jacobianInv[iJ+4] = -_jacobian[i20] / d12; // Jinv11
+ _jacobianInv[iJ+5] = _jacobian[i10] / d12; // Jinv12
+ } else
+ throw std::runtime_error("Could not invert Jacobian.");
+
+ // 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) {
+ const int iD = iQuadPt*numBasis*spaceDim + iBasis*spaceDim;
+ const int iDR = iQuadPt*numBasis*cellDim + iBasis*cellDim;
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int jDim=0; jDim < cellDim; ++jDim)
+ _basisDeriv[iD+iDim] += basisDerivRef[iDR+jDim] *
+ _jacobianInv[iJ+jDim*spaceDim+iDim];
+ } // for
+ } // for
+
+ PetscLogFlops(numQuadPts*(15 +
+ numBasis*spaceDim*2 +
+ numBasis*spaceDim*cellDim*2));
+} // computeGeometry
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2Din3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2Din3D.hh 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2Din3D.hh 2013-03-14 18:47:28 UTC (rev 21536)
@@ -57,6 +57,16 @@
void computeGeometry(const scalar_array& coordinatesCell,
const int cell);
+ /** Compute geometric quantities for a cell at quadrature points.
+ *
+ * @param coordinatesCell Array of coordinates of cell's vertices.
+ * @param coordinatesSize Size of coordinates array.
+ * @param cell Finite-element cell
+ */
+ void computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell);
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature3D.cc 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature3D.cc 2013-03-14 18:47:28 UTC (rev 21536)
@@ -194,4 +194,153 @@
} // computeGeometry
+// ----------------------------------------------------------------------
+// Compute geometric quantities for a cell at quadrature points.
+void
+pylith::feassemble::Quadrature3D::computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell)
+{ // computeGeometry
+ assert(coordinatesCell);
+
+ const int spaceDim = 3;
+ const int cellDim = 3;
+
+ assert(cellDim == _quadRefCell.cellDim());
+ assert(spaceDim == _quadRefCell.spaceDim());
+ const int numQuadPts = _quadRefCell.numQuadPts();
+ const int numBasis = _quadRefCell.numBasis();
+
+ const scalar_array& basis = _quadRefCell.basis();
+ const scalar_array& quadPtsRef = _quadRefCell.quadPtsRef();
+ const scalar_array& basisDerivRef = _quadRefCell.basisDerivRef();
+ const CellGeometry& geometry = _quadRefCell.refGeometry();
+
+ assert(numBasis*spaceDim == coordinatesSize);
+
+ zero();
+
+ // Loop over quadrature points
+ for (int iQuadPt=0; iQuadPt < numQuadPts; ++iQuadPt) {
+
+ // Compute coordinates of quadrature point in cell
+#if defined(ISOPARAMETRIC)
+ // x = sum[i=0,n-1] (Ni * xi)
+ // y = sum[i=0,n-1] (Ni * yi)
+ // z = sum[i=0,n-1] (Ni * zi)
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const PylithScalar valueBasis = basis[iQuadPt*numBasis+iBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _quadPts[iQuadPt*spaceDim+iDim] +=
+ valueBasis * coordinatesCell[iBasis*spaceDim+iDim];
+ } // for
+#else
+ geometry.ptsRefToGlobal(&_quadPts[iQuadPt*spaceDim],
+ &quadPtsRef[iQuadPt*cellDim],
+ coordinatesCell, spaceDim, 1);
+#endif
+
+#if defined(ISOPARAMETRIC)
+ // Compute Jacobian at quadrature point
+ // J = [dx/dp dx/dq dx/dr]
+ // [dy/dp dy/dq dy/dr]
+ // [dz/dp dz/dq dz/dr]
+ // dx/dp = sum[i=0,n-1] (dNi/dp * xi)
+ // dx/dq = sum[i=0,n-1] (dNi/dq * xi)
+ // dx/dr = sum[i=0,n-1] (dNi/dr * xi)
+ // dy/dp = sum[i=0,n-1] (dNi/dp * yi)
+ // dy/dq = sum[i=0,n-1] (dNi/dq * yi)
+ // dy/dr = sum[i=0,n-1] (dNi/dr * yi)
+ // dz/dp = sum[i=0,n-1] (dNi/dp * zi)
+ // dz/dq = sum[i=0,n-1] (dNi/dq * zi)
+ // dz/dr = sum[i=0,n-1] (dNi/dr * zi)
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ for (int iCol=0; iCol < cellDim; ++iCol) {
+ const PylithScalar deriv =
+ basisDerivRef[iQuadPt*numBasis*spaceDim+iBasis*cellDim+iCol];
+ for (int iRow=0; iRow < spaceDim; ++iRow)
+ _jacobian[iQuadPt*cellDim*spaceDim+iRow*cellDim+iCol] +=
+ deriv * coordinatesCell[iBasis*spaceDim+iRow];
+ } // for
+
+ // Compute determinant of Jacobian at quadrature point
+ // |J| = j00*(j11*j22-j12*j21) +
+ // -j01*(j10*j22-j12*j20) +
+ // j02*(j10*j21-j11*j20)
+ const int iJ = iQuadPt*cellDim*spaceDim;
+ const int i00 = iJ + 0*spaceDim + 0;
+ const int i01 = iJ + 0*spaceDim + 1;
+ const int i02 = iJ + 0*spaceDim + 2;
+ const int i10 = iJ + 1*spaceDim + 0;
+ const int i11 = iJ + 1*spaceDim + 1;
+ const int i12 = iJ + 1*spaceDim + 2;
+ const int i20 = iJ + 2*spaceDim + 0;
+ const int i21 = iJ + 2*spaceDim + 1;
+ const int i22 = iJ + 2*spaceDim + 2;
+ const PylithScalar det =
+ _jacobian[i00]*(_jacobian[i11]*_jacobian[i22] -
+ _jacobian[i12]*_jacobian[i21]) -
+ _jacobian[i01]*(_jacobian[i10]*_jacobian[i22] -
+ _jacobian[i12]*_jacobian[i20]) +
+ _jacobian[i02]*(_jacobian[i10]*_jacobian[i21] -
+ _jacobian[i11]*_jacobian[i20]);
+ _checkJacobianDet(det, cell);
+ _jacobianDet[iQuadPt] = det;
+#else
+ // Compute Jacobian and determinant of Jacobian at quadrature point
+ const int iJ = iQuadPt*cellDim*spaceDim;
+ const int i00 = iJ + 0*spaceDim + 0;
+ const int i01 = iJ + 0*spaceDim + 1;
+ const int i02 = iJ + 0*spaceDim + 2;
+ const int i10 = iJ + 1*spaceDim + 0;
+ const int i11 = iJ + 1*spaceDim + 1;
+ const int i12 = iJ + 1*spaceDim + 2;
+ const int i20 = iJ + 2*spaceDim + 0;
+ const int i21 = iJ + 2*spaceDim + 1;
+ const int i22 = iJ + 2*spaceDim + 2;
+ geometry.jacobian(&_jacobian[iQuadPt*cellDim*spaceDim],
+ &_jacobianDet[iQuadPt],
+ coordinatesCell, &quadPtsRef[iQuadPt*cellDim],
+ spaceDim, 1);
+ _checkJacobianDet(_jacobianDet[iQuadPt], cell);
+ const PylithScalar det = _jacobianDet[iQuadPt];
+#endif
+
+ // Compute inverse of Jacobian at quadrature point
+ _jacobianInv[i00] = (_jacobian[i11]*_jacobian[i22] -
+ _jacobian[i12]*_jacobian[i21]) / det;
+ _jacobianInv[i01] = (_jacobian[i02]*_jacobian[i21] -
+ _jacobian[i01]*_jacobian[i22]) / det;
+ _jacobianInv[i02] = (_jacobian[i01]*_jacobian[i12] -
+ _jacobian[i02]*_jacobian[i11]) / det;
+ _jacobianInv[i10] = (_jacobian[i12]*_jacobian[i20] -
+ _jacobian[i10]*_jacobian[i22]) / det;
+ _jacobianInv[i11] = (_jacobian[i00]*_jacobian[i22] -
+ _jacobian[i02]*_jacobian[i20]) / det;
+ _jacobianInv[i12] = (_jacobian[i02]*_jacobian[i10] -
+ _jacobian[i00]*_jacobian[i12]) / det;
+ _jacobianInv[i20] = (_jacobian[i10]*_jacobian[i21] -
+ _jacobian[i11]*_jacobian[i20]) / det;
+ _jacobianInv[i21] = (_jacobian[i01]*_jacobian[i20] -
+ _jacobian[i00]*_jacobian[i21]) / det;
+ _jacobianInv[i22] = (_jacobian[i00]*_jacobian[i11] -
+ _jacobian[i01]*_jacobian[i10]) / det;
+
+ // 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) {
+ const int iD = iQuadPt*numBasis*spaceDim + iBasis*spaceDim;
+ const int iDR = iQuadPt*numBasis*cellDim + iBasis*cellDim;
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int jDim=0; jDim < cellDim; ++jDim)
+ _basisDeriv[iD+iDim] += basisDerivRef[iDR+jDim] *
+ _jacobianInv[iJ+jDim*spaceDim+iDim];
+ } // for
+ } // for
+
+ PetscLogFlops(numQuadPts*(2+36 + numBasis*spaceDim*cellDim*4));
+} // computeGeometry
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature3D.hh 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature3D.hh 2013-03-14 18:47:28 UTC (rev 21536)
@@ -57,6 +57,16 @@
void computeGeometry(const scalar_array& coordinatesCell,
const int cell);
+ /** Compute geometric quantities for a cell at quadrature points.
+ *
+ * @param coordinatesCell Array of coordinates of cell's vertices.
+ * @param coordinatesSize Size of coordinates array.
+ * @param cell Finite-element cell
+ */
+ void computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell);
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureEngine.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureEngine.hh 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureEngine.hh 2013-03-14 18:47:28 UTC (rev 21536)
@@ -99,6 +99,17 @@
void computeGeometry(const scalar_array& coordinatesCell,
const int cell) = 0;
+ /** Compute geometric quantities for a cell at quadrature points.
+ *
+ * @param coordinatesCell Array of coordinates of cell's vertices.
+ * @param coordinatesSize Size of coordinates array.
+ * @param cell Finite-element cell
+ */
+ virtual
+ void computeGeometry(const PylithScalar* coordinatesCell,
+ const int coordinatesSize,
+ const int cell) = 0;
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc 2013-03-14 14:57:10 UTC (rev 21535)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc 2013-03-14 18:47:28 UTC (rev 21536)
@@ -238,13 +238,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
quadrature.retrieveGeometry(c);
#else
- PetscScalar *coords = PETSC_NULL;
- PetscInt coordsSize;
+ PetscScalar *coordsArray = NULL;
+ PetscInt coordsSize = 0;
- err = DMPlexVecGetClosure(dmMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
- err = DMPlexVecRestoreClosure(dmMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- quadrature.computeGeometry(coordinatesCell, c);
+ err = DMPlexVecGetClosure(dmMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
+ quadrature.computeGeometry(coordsArray, coordsSize, c);
+ err = DMPlexVecRestoreClosure(dmMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
#endif
const scalar_array& quadPts = quadrature.quadPts();
@@ -296,7 +295,8 @@
const int spaceDim = data.spaceDim;
const int numCells = data.numCells;
- scalar_array vertCoords(data.vertices, numBasis*spaceDim);
+ const PylithScalar* vertCoords = data.vertices;
+ const int vertCoordsSize = numBasis*spaceDim;
const PylithScalar* quadPtsE = data.quadPts;
const PylithScalar* jacobianE = data.jacobian;
const PylithScalar* jacobianDetE = data.jacobianDet;
@@ -347,7 +347,7 @@
#endif
quadrature.initializeGeometry();
- quadrature.computeGeometry(vertCoords, 0);
+ quadrature.computeGeometry(vertCoords, vertCoordsSize, 0);
size_t size = 0;
More information about the CIG-COMMITS
mailing list