[cig-commits] r16232 - short/3D/PyLith/trunk/libsrc/feassemble
brad at geodynamics.org
brad at geodynamics.org
Thu Feb 4 16:17:29 PST 2010
Author: brad
Date: 2010-02-04 16:17:29 -0800 (Thu, 04 Feb 2010)
New Revision: 16232
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
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
Log:
Optimized routines using constants for known cellDim and spaceDim.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2010-02-04 19:29:24 UTC (rev 16231)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2010-02-05 00:17:29 UTC (rev 16232)
@@ -673,17 +673,20 @@
const int stressSize = 3;
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const int iQs = iQuad*stressSize;
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad*stressSize+0];
- const double s22 = stress[iQuad*stressSize+1];
- const double s12 = stress[iQuad*stressSize+2];
+ const double s11 = stress[iQs ];
+ const double s22 = stress[iQs+1];
+ const double s12 = stress[iQs+2];
for (int iBasis=0, iQ=iQuad*numBasis*spaceDim;
iBasis < numBasis;
++iBasis) {
- const double N1 = wt*basisDeriv[iQ+iBasis*spaceDim ];
- const double N2 = wt*basisDeriv[iQ+iBasis*spaceDim+1];
- _cellVector[iBasis*spaceDim ] -= N1*s11 + N2*s12;
- _cellVector[iBasis*spaceDim+1] -= N1*s12 + N2*s22;
+ const int iBlock = iBasis*spaceDim;
+ const double N1 = wt*basisDeriv[iQ+iBlock ];
+ const double N2 = wt*basisDeriv[iQ+iBlock+1];
+
+ _cellVector[iBlock ] -= N1*s11 + N2*s12;
+ _cellVector[iBlock+1] -= N1*s12 + N2*s22;
} // for
} // for
PetscLogFlops(numQuadPts*(1+numBasis*(8+2+9)));
@@ -695,34 +698,37 @@
pylith::feassemble::IntegratorElasticity::_elasticityResidual3D(
const double_array& stress)
{ // _elasticityResidual3D
+ const int spaceDim = 3;
+ const int cellDim = 3;
+ const int stressSize = 6;
+
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
- const int cellDim = _quadrature->cellDim();
const double_array& quadWts = _quadrature->quadWts();
const double_array& jacobianDet = _quadrature->jacobianDet();
const double_array& basisDeriv = _quadrature->basisDeriv();
- assert(3 == cellDim);
+ assert(_quadrature->spaceDim() == spaceDim);
+ assert(_quadrature->cellDim() == cellDim);
assert(quadWts.size() == numQuadPts);
- const int stressSize = 6;
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const int iQs = iQuad * stressSize;
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad*stressSize+0];
- const double s22 = stress[iQuad*stressSize+1];
- const double s33 = stress[iQuad*stressSize+2];
- const double s12 = stress[iQuad*stressSize+3];
- const double s23 = stress[iQuad*stressSize+4];
- const double s13 = stress[iQuad*stressSize+5];
+ const double s11 = stress[iQs ];
+ const double s22 = stress[iQs+1];
+ const double s33 = stress[iQs+2];
+ const double s12 = stress[iQs+3];
+ const double s23 = stress[iQs+4];
+ const double s13 = stress[iQs+5];
for (int iBasis=0, iQ=iQuad*numBasis*spaceDim;
- iBasis < numBasis;
- ++iBasis) {
+ iBasis < numBasis;
+ ++iBasis) {
const int iBlock = iBasis*spaceDim;
- const double N1 = wt*basisDeriv[iQ+iBasis*spaceDim+0];
- const double N2 = wt*basisDeriv[iQ+iBasis*spaceDim+1];
- const double N3 = wt*basisDeriv[iQ+iBasis*spaceDim+2];
+ const double N1 = wt*basisDeriv[iQ+iBlock+0];
+ const double N2 = wt*basisDeriv[iQ+iBlock+1];
+ const double N3 = wt*basisDeriv[iQ+iBlock+2];
_cellVector[iBlock ] -= N1*s11 + N2*s12 + N3*s13;
_cellVector[iBlock+1] -= N1*s12 + N2*s22 + N3*s23;
@@ -910,17 +916,19 @@
pylith::feassemble::IntegratorElasticity::_elasticityJacobian3D(
const double_array& elasticConsts)
{ // _elasticityJacobian3D
+ const int spaceDim = 3;
+ const int cellDim = 3;
+ const int numConsts = 21;
+
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
- const int cellDim = _quadrature->cellDim();
const double_array& quadWts = _quadrature->quadWts();
const double_array& jacobianDet = _quadrature->jacobianDet();
const double_array& basisDeriv = _quadrature->basisDeriv();
- assert(3 == cellDim);
+ assert(_quadrature->spaceDim() == spaceDim);
+ assert(_quadrature->cellDim() == cellDim);
assert(quadWts.size() == numQuadPts);
- const int numConsts = 21;
// Compute Jacobian for consistent tangent matrix
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
@@ -1098,12 +1106,12 @@
assert(0 != strain);
const int dim = 2;
+ const int strainSize = 3;
assert(basisDeriv.size() == numQuadPts*numBasis*dim);
assert(disp.size() == numBasis*dim);
(*strain) = 0.0;
- const int strainSize = 3;
for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
for (int iBasis=0, iQ=iQuad*numBasis*dim; iBasis < numBasis; ++iBasis) {
(*strain)[iQuad*strainSize+0] +=
@@ -1128,29 +1136,29 @@
assert(0 != strain);
const int dim = 3;
+ const int strainSize = 6;
assert(basisDeriv.size() == numQuadPts*numBasis*dim);
assert(disp.size() == numBasis*dim);
(*strain) = 0.0;
- const int strainSize = 6;
for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
for (int iBasis=0, iQ=iQuad*numBasis*dim; iBasis < numBasis; ++iBasis) {
- (*strain)[iQuad*strainSize+0] +=
- basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim ];
- (*strain)[iQuad*strainSize+1] +=
- basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+1];
- (*strain)[iQuad*strainSize+2] +=
- basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim+2];
- (*strain)[iQuad*strainSize+3] +=
- 0.5 * (basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim ] +
- basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+1]);
- (*strain)[iQuad*strainSize+4] +=
- 0.5 * (basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim+1] +
- basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+2]);
- (*strain)[iQuad*strainSize+5] +=
- 0.5 * (basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim ] +
- basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+2]);
+ (*strain)[iQuad*strainSize ] += basisDeriv[iQ+iBasis*dim]
+ * disp[iBasis*dim];
+ (*strain)[iQuad*strainSize+1] += basisDeriv[iQ+iBasis*dim+1]
+ * disp[iBasis*dim+1];
+ (*strain)[iQuad*strainSize+2] += basisDeriv[iQ+iBasis*dim+2]
+ * disp[iBasis*dim+2];
+ (*strain)[iQuad*strainSize+3] +=
+ 0.5 * (basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim ] +
+ basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+1]);
+ (*strain)[iQuad*strainSize+4] +=
+ 0.5 * (basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim+1] +
+ basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+2]);
+ (*strain)[iQuad*strainSize+5] +=
+ 0.5 * (basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim ] +
+ basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+2]);
} // for
} // calcTotalStrain3D
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc 2010-02-04 19:29:24 UTC (rev 16231)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc 2010-02-05 00:17:29 UTC (rev 16232)
@@ -47,20 +47,20 @@
pylith::feassemble::Quadrature0D::computeGeometry(const double_array& coordinatesCell,
const int cell)
{ // computeGeometry
- const int cellDim = _quadRefCell.cellDim();
- const int spaceDim = _quadRefCell.spaceDim();
- const int numQuadPts = _quadRefCell.numQuadPts();
- const int numBasis = _quadRefCell.numBasis();
+ 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(coordinatesCell.size() == numBasis*spaceDim);
const double_array& basisDerivRef = _quadRefCell.basisDerivRef();
- assert(0 == cellDim);
- assert(1 == numQuadPts);
- assert(1 == numBasis);
-
zero();
- assert(numBasis*1 == coordinatesCell.size());
_quadPts = coordinatesCell;
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc 2010-02-04 19:29:24 UTC (rev 16231)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc 2010-02-05 00:17:29 UTC (rev 16232)
@@ -49,8 +49,9 @@
pylith::feassemble::Quadrature1D::computeGeometry(const double_array& coordinatesCell,
const int cell)
{ // computeGeometry
- const int cellDim = _quadRefCell.cellDim();
- const int spaceDim = _quadRefCell.spaceDim();
+ const int cellDim = 1;
+ const int spaceDim = 1;
+
const int numQuadPts = _quadRefCell.numQuadPts();
const int numBasis = _quadRefCell.numBasis();
@@ -59,21 +60,21 @@
const double_array& basisDerivRef = _quadRefCell.basisDerivRef();
const CellGeometry& geometry = _quadRefCell.refGeometry();
+ assert(_quadRefCell.cellDim() == cellDim);
+ assert(_quadRefCell.spaceDim() == spaceDim);
assert(numBasis*spaceDim == coordinatesCell.size());
- assert(1 == cellDim);
- assert(1 == spaceDim);
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[iQuadPt*numBasis+iBasis]*coordinatesCell[iBasis];
+ _quadPts[iQuadPt] += basis[iQ+iBasis]*coordinatesCell[iBasis];
#else
geometry.coordsRefToGlobal(&_quadPts[iQuadPt], &quadPtsRef[iQuadPt],
&coordinatesCell[0], spaceDim);
@@ -83,8 +84,7 @@
// 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[iQuadPt*numBasis+iBasis] * coordinatesCell[iBasis];
+ _jacobian[iQuadPt] += basisDerivRef[iQ+iBasis] * coordinatesCell[iBasis];
// Compute determinant of Jacobian at quadrature point
// |J| = j00
@@ -111,9 +111,8 @@
// coordinates
// dNi/dx = dNi/dp dp/dx + dNi/dq dq/dx + dNi/dr dr/dx
for (int iBasis=0; iBasis < numBasis; ++iBasis)
- _basisDeriv[iQuadPt*numBasis+iBasis] +=
- basisDerivRef[iQuadPt*numBasis+iBasis] *
- _jacobianInv[iQuadPt];
+ _basisDeriv[iQ+iBasis] +=
+ basisDerivRef[iQ+iBasis] * _jacobianInv[iQuadPt];
} // for
PetscLogFlops(numQuadPts * (1 + numBasis * 4));
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc 2010-02-04 19:29:24 UTC (rev 16231)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc 2010-02-05 00:17:29 UTC (rev 16232)
@@ -49,8 +49,9 @@
pylith::feassemble::Quadrature1Din2D::computeGeometry(const double_array& coordinatesCell,
const int cell)
{ // computeGeometry
- const int cellDim = _quadRefCell.cellDim();
- const int spaceDim = _quadRefCell.spaceDim();
+ const int cellDim = 1;
+ const int spaceDim = 2;
+
const int numQuadPts = _quadRefCell.numQuadPts();
const int numBasis = _quadRefCell.numBasis();
@@ -59,21 +60,22 @@
const double_array& basisDerivRef = _quadRefCell.basisDerivRef();
const CellGeometry& geometry = _quadRefCell.refGeometry();
+ assert(_quadRefCell.cellDim() == cellDim);
+ assert(_quadRefCell.spaceDim() == spaceDim);
assert(numBasis*spaceDim == coordinatesCell.size());
- assert(1 == cellDim);
- assert(2 == spaceDim);
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 double valueBasis = basis[iQuadPt*numBasis+iBasis];
+ const double valueBasis = basis[iQ+iBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
_quadPts[iQuadPt*spaceDim+iDim] +=
valueBasis * coordinatesCell[iBasis*spaceDim+iDim];
@@ -91,7 +93,7 @@
// 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 double deriv = basisDerivRef[iQuadPt*numBasis+iBasis];
+ const double deriv = basisDerivRef[iQ+iBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
_jacobian[iQuadPt*spaceDim+iDim] +=
deriv * coordinatesCell[iBasis*spaceDim+iDim];
@@ -123,12 +125,15 @@
// 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 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[iQuadPt*numBasis*spaceDim+iBasis*spaceDim+iDim] +=
- basisDerivRef[iQuadPt*numBasis*cellDim+iBasis*cellDim+jDim] *
- _jacobianInv[iQuadPt*cellDim*spaceDim+jDim*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 +
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc 2010-02-04 19:29:24 UTC (rev 16231)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc 2010-02-05 00:17:29 UTC (rev 16232)
@@ -49,9 +49,9 @@
pylith::feassemble::Quadrature1Din3D::computeGeometry(const double_array& coordinatesCell,
const int cell)
{ // computeGeometry
+ const int cellDim = 1;
+ const int spaceDim = 3;
- const int cellDim = _quadRefCell.cellDim();
- const int spaceDim = _quadRefCell.spaceDim();
const int numQuadPts = _quadRefCell.numQuadPts();
const int numBasis = _quadRefCell.numBasis();
@@ -60,14 +60,15 @@
const double_array& basisDerivRef = _quadRefCell.basisDerivRef();
const CellGeometry& geometry = _quadRefCell.refGeometry();
+ assert(_quadRefCell.cellDim() == cellDim);
+ assert(_quadRefCell.spaceDim() == spaceDim);
assert(numBasis*spaceDim == coordinatesCell.size());
- assert(1 == cellDim);
- assert(3 == spaceDim);
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)
@@ -75,7 +76,7 @@
// y = sum[i=0,n-1] (Ni * yi)
// z = sum[i=0,n-1] (Ni * zi)
for (int iBasis=0; iBasis < numBasis; ++iBasis) {
- const double valueBasis = basis[iQuadPt*numBasis+iBasis];
+ const double valueBasis = basis[iQ+iBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
_quadPts[iQuadPt*spaceDim+iDim] +=
valueBasis * coordinatesCell[iBasis*spaceDim+iDim];
@@ -95,7 +96,7 @@
// 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 double deriv = basisDerivRef[iQuadPt*numBasis+iBasis];
+ const double deriv = basisDerivRef[iQ+iBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
_jacobian[iQuadPt*spaceDim+iDim] +=
deriv * coordinatesCell[iBasis*spaceDim+iDim];
@@ -127,12 +128,15 @@
// 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 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[iQuadPt*numBasis*spaceDim+iBasis*spaceDim+iDim] +=
- basisDerivRef[iQuadPt*numBasis*cellDim+iBasis*cellDim+jDim] *
- _jacobianInv[iQuadPt*cellDim*spaceDim+jDim*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 +
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc 2010-02-04 19:29:24 UTC (rev 16231)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc 2010-02-05 00:17:29 UTC (rev 16232)
@@ -49,8 +49,9 @@
pylith::feassemble::Quadrature2D::computeGeometry(const double_array& coordinatesCell,
const int cell)
{ // computeGeometry
- const int cellDim = _quadRefCell.cellDim();
- const int spaceDim = _quadRefCell.spaceDim();
+ const int spaceDim = 2;
+ const int cellDim = 2;
+
const int numQuadPts = _quadRefCell.numQuadPts();
const int numBasis = _quadRefCell.numBasis();
@@ -59,9 +60,9 @@
const double_array& basisDerivRef = _quadRefCell.basisDerivRef();
const CellGeometry& geometry = _quadRefCell.refGeometry();
+ assert(_quadRefCell.cellDim() == cellDim);
+ assert(_quadRefCell.spaceDim() == spaceDim);
assert(numBasis*spaceDim == coordinatesCell.size());
- assert(2 == cellDim);
- assert(2 == spaceDim);
zero();
// Loop over quadrature points
@@ -118,13 +119,6 @@
&_jacobianDet[iQuadPt],
&coordinatesCell[0], &quadPtsRef[iQuadPt*cellDim], spaceDim);
_checkJacobianDet(_jacobianDet[iQuadPt], cell);
-
- 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 double det = _jacobianDet[iQuadPt];
#endif
// Compute inverse of Jacobian at quadrature point
@@ -138,12 +132,14 @@
// 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 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[iQuadPt*numBasis*spaceDim+iBasis*spaceDim+iDim] +=
- basisDerivRef[iQuadPt*numBasis*cellDim+iBasis*cellDim+jDim] *
- _jacobianInv[iQuadPt*cellDim*spaceDim+jDim*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 +
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc 2010-02-04 19:29:24 UTC (rev 16231)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc 2010-02-05 00:17:29 UTC (rev 16232)
@@ -51,9 +51,9 @@
pylith::feassemble::Quadrature2Din3D::computeGeometry(const double_array& coordinatesCell,
const int cell)
{ // computeGeometry
+ const int cellDim = 2;
+ const int spaceDim = 3;
- const int cellDim = _quadRefCell.cellDim();
- const int spaceDim = _quadRefCell.spaceDim();
const int numQuadPts = _quadRefCell.numQuadPts();
const int numBasis = _quadRefCell.numBasis();
@@ -63,9 +63,9 @@
const CellGeometry& geometry = _quadRefCell.refGeometry();
const double minJacobian = _quadRefCell.minJacobian();
+ assert(_quadRefCell.cellDim() == cellDim);
+ assert(_quadRefCell.spaceDim() == spaceDim);
assert(numBasis*spaceDim == coordinatesCell.size());
- assert(2 == cellDim);
- assert(3 == spaceDim);
zero();
// Loop over quadrature points
@@ -140,14 +140,6 @@
&_jacobianDet[iQuadPt],
&coordinatesCell[0], &quadPtsRef[iQuadPt*cellDim], spaceDim);
_checkJacobianDet(_jacobianDet[iQuadPt], cell);
-
- 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;
#endif
// Compute inverse of Jacobian at quadrature point
@@ -219,12 +211,14 @@
// 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 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[iQuadPt*numBasis*spaceDim+iBasis*spaceDim+iDim] +=
- basisDerivRef[iQuadPt*numBasis*cellDim + iBasis*cellDim+jDim] *
- _jacobianInv[iQuadPt*cellDim*spaceDim+jDim*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 +
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc 2010-02-04 19:29:24 UTC (rev 16231)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc 2010-02-05 00:17:29 UTC (rev 16232)
@@ -49,9 +49,11 @@
pylith::feassemble::Quadrature3D::computeGeometry(const double_array& coordinatesCell,
const int cell)
{ // computeGeometry
+ const int spaceDim = 3;
+ const int cellDim = 3;
- const int cellDim = _quadRefCell.cellDim();
- const int spaceDim = _quadRefCell.spaceDim();
+ assert(cellDim == _quadRefCell.cellDim());
+ assert(spaceDim == _quadRefCell.spaceDim());
const int numQuadPts = _quadRefCell.numQuadPts();
const int numBasis = _quadRefCell.numBasis();
@@ -61,8 +63,6 @@
const CellGeometry& geometry = _quadRefCell.refGeometry();
assert(numBasis*spaceDim == coordinatesCell.size());
- assert(3 == cellDim);
- assert(3 == spaceDim);
zero();
// Loop over quadrature points
@@ -137,18 +137,6 @@
&_jacobianDet[iQuadPt],
&coordinatesCell[0], &quadPtsRef[iQuadPt*cellDim], spaceDim);
_checkJacobianDet(_jacobianDet[iQuadPt], cell);
-
- 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 double det = _jacobianDet[iQuadPt];
#endif
// Compute inverse of Jacobian at quadrature point
@@ -174,16 +162,17 @@
// 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 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[iQuadPt*numBasis*spaceDim+iBasis*spaceDim+iDim] +=
- basisDerivRef[iQuadPt*numBasis*cellDim+iBasis*cellDim+jDim] *
- _jacobianInv[iQuadPt*cellDim*spaceDim+jDim*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));
+ PetscLogFlops(numQuadPts*(2+36 + numBasis*spaceDim*cellDim*4));
} // computeGeometry
More information about the CIG-COMMITS
mailing list