[cig-commits] r22004 - in short/3D/PyLith/trunk: libsrc/pylith/bc libsrc/pylith/faults libsrc/pylith/feassemble libsrc/pylith/friction libsrc/pylith/materials libsrc/pylith/topology modulesrc/feassemble unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Wed May 8 15:30:44 PDT 2013
Author: brad
Date: 2013-05-08 15:30:43 -0700 (Wed, 08 May 2013)
New Revision: 22004
Modified:
short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/CellGeometry.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.cc
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/libsrc/pylith/feassemble/QuadratureRefCell.cc
short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellQpQsIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
short/3D/PyLith/trunk/modulesrc/feassemble/CellGeometry.i
short/3D/PyLith/trunk/modulesrc/feassemble/GeometryHex3D.i
short/3D/PyLith/trunk/modulesrc/feassemble/GeometryLine1D.i
short/3D/PyLith/trunk/modulesrc/feassemble/GeometryLine2D.i
short/3D/PyLith/trunk/modulesrc/feassemble/GeometryLine3D.i
short/3D/PyLith/trunk/modulesrc/feassemble/GeometryPoint1D.i
short/3D/PyLith/trunk/modulesrc/feassemble/GeometryPoint2D.i
short/3D/PyLith/trunk/modulesrc/feassemble/GeometryPoint3D.i
short/3D/PyLith/trunk/modulesrc/feassemble/GeometryQuad2D.i
short/3D/PyLith/trunk/modulesrc/feassemble/GeometryQuad3D.i
short/3D/PyLith/trunk/modulesrc/feassemble/GeometryTet3D.i
short/3D/PyLith/trunk/modulesrc/feassemble/GeometryTri2D.i
short/3D/PyLith/trunk/modulesrc/feassemble/GeometryTri3D.i
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestCellGeometry.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestGeometryPoint1D.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestGeometryPoint2D.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestGeometryPoint3D.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticity.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticityLgDeform.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature0D.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadratureEngine.cc
Log:
Code cleanup. Improve quadrature interface to eliminate copies.
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -34,12 +34,10 @@
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
-#include <cstring> // USES memcpy()
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
-//#define PRECOMPUTE_GEOMETRY
//#define DETAILED_EVENT_LOGGING
// ----------------------------------------------------------------------
@@ -134,7 +132,6 @@
// Container for data returned in query of database
scalar_array queryData(numValues);
- scalar_array quadPtRef(cellDim);
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Container for damping constants for current cell
@@ -142,7 +139,6 @@
topology::Field<topology::SubMesh>& dampingConsts = _parameters->get("damping constants");
topology::VecVisitorMesh dampingConstsVisitor(dampingConsts);
- scalar_array coordinatesCell(numBasis*spaceDim);
PetscScalar *coordsCell = NULL;
PetscInt coordsSize = 0;
topology::CoordsVisitor coordsVisitor(dmSubMesh);
@@ -153,7 +149,7 @@
assert(_normalizer->timeScale() > 0);
const PylithScalar velocityScale = _normalizer->lengthScale() / _normalizer->timeScale();
- const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
+ const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();assert(cs);
// Compute quadrature information
_quadrature->initializeGeometry();
@@ -163,11 +159,7 @@
for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry information for current cell
coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
- for (int i=0; i < coordsSize; ++i) { // :TODO: Remove copy
- coordinatesCell[i] = coordsCell[i];
- } // for
_quadrature->computeGeometry(coordsCell, coordsSize, c);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
const PetscInt doff = dampingConstsVisitor.sectionOffset(c);
assert(fiberDim == dampingConstsVisitor.sectionDof(c));
@@ -204,8 +196,7 @@
dampingConstsLocal[spaceDim-1] = constNormal;
// Compute normal/tangential orientation
- memcpy(&quadPtRef[0], &quadPtsRef[iQuad*cellDim], cellDim*sizeof(PylithScalar));
- cellGeometry.jacobian(&jacobian, &jacobianDet, coordinatesCell, quadPtRef);
+ cellGeometry.jacobian(&jacobian, &jacobianDet, coordsCell, numBasis, spaceDim, &quadPtsRef[iQuad*cellDim], cellDim);
cellGeometry.orientation(&orientation, jacobian, jacobianDet, up);
assert(jacobianDet > 0.0);
orientation /= jacobianDet;
@@ -219,6 +210,7 @@
dampingConstsArray[doff+iQuad*spaceDim+iDim] = fabs(dampingConstsArray[doff+iQuad*spaceDim+iDim]);
} // for
} // for
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
} // for
_db->close();
@@ -277,11 +269,9 @@
submeshIS.deallocate();
-#if !defined(PRECOMPUTE_GEOMETRY)
PetscScalar* coordsCell = NULL;
PetscInt coordsSize = 0;
topology::CoordsVisitor coordsVisitor(dmSubMesh);
-#endif
// Get 'surface' cells (1 dimension lower than top-level cells)
topology::Stratum cellsStratum(dmSubMesh, topology::Stratum::HEIGHT, 1);
@@ -298,13 +288,9 @@
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
-#if defined(PRECOMPUTE_GEOMETRY)
-#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
-#else
coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
_quadrature->computeGeometry(coordsCell, coordsSize, c);
coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
-#endif
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -423,11 +409,9 @@
submeshIS.deallocate();
-#if !defined(PRECOMPUTE_GEOMETRY)
PetscScalar *coordsCell = NULL;
PetscInt coordsSize = 0;
topology::CoordsVisitor coordsVisitor(dmSubMesh);
-#endif
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -439,13 +423,9 @@
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
-#if defined(PRECOMPUTE_GEOMETRY)
-#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
-#else
coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
_quadrature->computeGeometry(coordsCell, coordsSize, c);
coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
-#endif
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -565,11 +545,9 @@
// Allocate matrix for cell values.
_initCellMatrix();
-#if !defined(PRECOMPUTE_GEOMETRY)
PetscScalar *coordsCell = NULL;
PetscInt coordsSize = 0;
topology::CoordsVisitor coordsVisitor(dmSubMesh);
-#endif
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -581,13 +559,10 @@
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
-#if defined(PRECOMPUTE_GEOMETRY)
-#error("Code for PRECOMPUTE_GEOMETRY not implemented")
-#else
coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
_quadrature->computeGeometry(coordsCell, coordsSize, c);
coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
-#endif
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
_logger->eventBegin(restrictEvent);
@@ -703,11 +678,9 @@
submeshIS.deallocate();
-#if !defined(PRECOMPUTE_GEOMETRY)
PetscScalar* coordsCell = NULL;
PetscInt coordsSize = 0;
topology::CoordsVisitor coordsVisitor(dmSubMesh);
-#endif
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -719,13 +692,10 @@
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
-#if defined(PRECOMPUTE_GEOMETRY)
-#error("Code for PRECOMPUTE_GEOMETRY not implemented");
-#else
coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
_quadrature->computeGeometry(coordsCell, coordsSize, c);
coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
-#endif
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
_logger->eventBegin(restrictEvent);
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -39,8 +39,6 @@
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
-//#define PRECOMPUTE_GEOMETRY
-
// ----------------------------------------------------------------------
// Default constructor.
pylith::bc::Neumann::Neumann(void)
@@ -123,21 +121,15 @@
topology::VecVisitorSubMesh residualVisitor(residual, submeshIS);
submeshIS.deallocate();
-#if !defined(PRECOMPUTE_GEOMETRY)
PetscScalar* coordsCell = NULL;
PetscInt coordsSize = 0;
topology::CoordsVisitor coordsVisitor(dmSubMesh);
-#endif
// Loop over faces and integrate contribution from each face
for(PetscInt c = cStart; c < cEnd; ++c) {
-#if defined(PRECOMPUTE_GEOMETRY)
-#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
-#else
coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
_quadrature->computeGeometry(coordsCell, coordsSize, c);
coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
-#endif
// Reset element vector to zero
_resetCellVector();
@@ -437,20 +429,14 @@
// Compute quadrature information
_quadrature->initializeGeometry();
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->computeGeometry(*_boundaryMesh, cells);
-#endif
// Loop over cells in boundary mesh and perform queries.
for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry information for current cell
-#if defined(PRECOMPUTE_GEOMETRY)
-#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
-#else
coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
_quadrature->computeGeometry(coordsCell, coordsSize, c);
coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
-#endif
+
const scalar_array& quadPtsNondim = _quadrature->quadPts();
quadPtsGlobal = quadPtsNondim;
_normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(), lengthScale);
@@ -511,7 +497,6 @@
const int numBasis = _quadrature->numBasis();
const int numQuadPts = _quadrature->numQuadPts();
const int spaceDim = cellGeometry.spaceDim();
- scalar_array quadPtRef(cellDim);
const scalar_array& quadPtsRef = _quadrature->quadPtsRef();
// Containers for orientation information
@@ -522,7 +507,6 @@
scalar_array orientation(orientationSize);
// Get coordinates.
- scalar_array coordinatesCell(numBasis*spaceDim);
PetscScalar* coordsCell = NULL;
PetscInt coordsSize = 0;
topology::CoordsVisitor coordsVisitor(dmSubMesh);
@@ -546,23 +530,12 @@
// rotate corresponding traction vector from local to global coordinates.
for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry information for current cell
-#if defined(PRECOMPUTE_GEOMETRY)
-#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
-#else
coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
- for (PetscInt i=0; i < coordsSize; ++i) { // :TODO: Remove copy.
- coordinatesCell[i] = coordsCell[i];
- } // for
- _quadrature->computeGeometry(coordinatesCell, c);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
-#endif
+ _quadrature->computeGeometry(coordsCell, coordsSize, c);
+
for(int iQuad=0, iRef=0, iSpace=0; iQuad < numQuadPts; ++iQuad, iRef+=cellDim, iSpace+=spaceDim) {
// Compute Jacobian and determinant at quadrature point, then get orientation.
- memcpy(&quadPtRef[0], &quadPtsRef[iRef], cellDim*sizeof(PylithScalar));
-#if defined(PRECOMPUTE_GEOMETRY)
-#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
-#endif
- cellGeometry.jacobian(&jacobian, &jacobianDet, coordinatesCell, quadPtRef);
+ cellGeometry.jacobian(&jacobian, &jacobianDet, coordsCell, numBasis, spaceDim, &quadPtsRef[iRef], cellDim);
cellGeometry.orientation(&orientation, jacobian, jacobianDet, up);
assert(jacobianDet > 0.0);
orientation /= jacobianDet;
@@ -615,6 +588,7 @@
} // for
} // if
} // for
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
} // for
delete initialVisitor; initialVisitor = 0;
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -1273,7 +1273,6 @@
const scalar_array& quadWts = _quadrature->quadWts();
scalar_array jacobian(jacobianSize);
PylithScalar jacobianDet = 0;
- scalar_array refCoordsVertex(cohesiveDim);
assert(cohesiveDim == _quadrature->cellDim());
assert(spaceDim == _quadrature->spaceDim());
@@ -1310,7 +1309,6 @@
PetscScalar* orientationArray = orientationVisitor.localArray();
// Get section containing coordinates of vertices
- scalar_array coordinatesCell(numBasis * spaceDim); // :TODO: Remove copy.
topology::CoordsVisitor coordsVisitor(faultDMMesh);
PetscScalar *coordsCell = NULL;
PetscInt coordsSize = 0;
@@ -1325,9 +1323,6 @@
// Get orientations at fault cell's vertices.
coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
- for(PetscInt i = 0; i < coordsSize; ++i) { // :TODO: Remove copy.
- coordinatesCell[i] = coordsCell[i];
- } // for
PetscErrorCode err = DMPlexGetTransitiveClosure(faultDMMesh, c, PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
@@ -1342,8 +1337,7 @@
closureSize = q;
for(PetscInt v = 0; v < closureSize; ++v) {
// Compute Jacobian and determinant of Jacobian at vertex
- memcpy(&refCoordsVertex[0], &verticesRef[v * cohesiveDim], cohesiveDim * sizeof(PylithScalar));
- cellGeometry.jacobian(&jacobian, &jacobianDet, coordinatesCell, refCoordsVertex);
+ cellGeometry.jacobian(&jacobian, &jacobianDet, coordsCell, numBasis, spaceDim, &verticesRef[v*cohesiveDim], cohesiveDim);
// Compute orientation
cellGeometry.orientation(&orientationVertex, jacobian, jacobianDet, up);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/CellGeometry.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/CellGeometry.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/CellGeometry.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -126,13 +126,19 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
virtual
void jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const = 0;
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const = 0;
/** Compute Jacobian at location in cell.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -43,7 +43,6 @@
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
-//#define PRECOMPUTE_GEOMETRY
//#define DETAILED_EVENT_LOGGING
// ----------------------------------------------------------------------
@@ -248,13 +247,10 @@
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(cell);
-#else
coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
_quadrature->computeGeometry(coordsCell, coordsSize, cell);
coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
-#endif
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
_logger->eventBegin(stateVarsEvent);
@@ -355,7 +351,7 @@
#endif
// Compute B(transpose) * sigma, first computing strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispAdjCell, numBasis, numQuadPts);
+ calcTotalStrainFn(&strainCell, basisDeriv, &dispAdjCell[0], numBasis, spaceDim, numQuadPts);
const scalar_array& stressCell = _material->calcStress(strainCell, false);
#if defined(DETAILED_EVENT_LOGGING)
@@ -469,13 +465,9 @@
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(cell);
-#else
coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
_quadrature->computeGeometry(coordsCell, coordsSize, cell);
coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
-#endif
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -42,8 +42,6 @@
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
-//#define PRECOMPUTE_GEOMETRY
-
// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::ElasticityExplicitLgDeform::ElasticityExplicitLgDeform(void) :
@@ -238,14 +236,8 @@
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
// Compute geometry information for current cell
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(cell);
-#else
coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
- for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coordsCell[i];} // :TODO: Remove copy
_quadrature->computeGeometry(coordsCell, coordsSize, cell);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
-#endif
// Get state variables for cell.
_material->retrievePropsAndVars(cell);
@@ -327,7 +319,7 @@
dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
// Compute B(transpose) * sigma, first computing strains
- _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispAdjCell, numBasis, numQuadPts, spaceDim);
+ _calcDeformation(&deformCell, basisDeriv, coordsCell, &dispAdjCell[0], numBasis, numQuadPts, spaceDim);
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
const scalar_array& stressCell = _material->calcStress(strainCell, true);
@@ -335,6 +327,8 @@
// Assemble cell contribution into field
residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), cell, ADD_VALUES);
+
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
} // for
_logger->eventEnd(computeEvent);
@@ -418,13 +412,9 @@
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
// Compute geometry information for current cell
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(cell);
-#else
coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
_quadrature->computeGeometry(coordsCell, coordsSize, cell);
coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
-#endif
// Get state variables for cell.
_material->retrievePropsAndVars(cell);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -261,7 +261,7 @@
// residualSection->view("After gravity contribution");
// Compute B(transpose) * sigma, first computing strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispTpdtCell, numBasis, numQuadPts);
+ calcTotalStrainFn(&strainCell, basisDeriv, &dispTpdtCell[0], numBasis, spaceDim, numQuadPts);
const scalar_array& stressCell = _material->calcStress(strainCell, true);
CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell);
@@ -416,7 +416,7 @@
dispIncrVisitor.restoreClosure(&dispIncrCell, &dispIncrSize, cell);
// Compute strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispTpdtCell, numBasis, numQuadPts);
+ calcTotalStrainFn(&strainCell, basisDeriv, &dispTpdtCell[0], numBasis, spaceDim, numQuadPts);
// Get "elasticity" matrix at quadrature points for this cell
const scalar_array& elasticConsts = _material->calcDerivElastic(strainCell);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -45,8 +45,6 @@
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
-//#define PRECOMPUTE_GEOMETRY
-
// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::ElasticityImplicitLgDeform::ElasticityImplicitLgDeform(void) :
@@ -180,19 +178,12 @@
// Setup field visitors.
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
- PetscScalar* dispCell = NULL;
- PetscInt dispSize = 0;
topology::VecVisitorMesh dispIncrVisitor(fields->get("dispIncr(t->t+dt)"));
- PetscScalar* dispIncrCell = NULL;
- PetscInt dispIncrSize = 0;
topology::VecVisitorMesh residualVisitor(residual);
- scalar_array coordinatesCell(numBasis*spaceDim);
topology::CoordsVisitor coordsVisitor(dmMesh);
- PetscScalar *coordsCell = NULL;
- PetscInt coordsSize = 0;
assert(_normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -206,15 +197,12 @@
// Loop over cells
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
+
// Compute geometry information for current cell
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(cell);
-#else
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+ PetscScalar *coordsCell = NULL;
+ PetscInt coordsSize = 0;
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
_quadrature->computeGeometry(coordsCell, coordsSize, cell);
- for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coordsCell[i];} // :TODO: Remove copy
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
-#endif
// Get state variables for cell.
_material->retrievePropsAndVars(cell);
@@ -223,8 +211,13 @@
_resetCellVector();
// Restrict input fields to cell
- dispVisitor.getClosure(&dispCell, &dispSize, cell);
- dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);
+ PetscScalar* dispCell = NULL;
+ PetscInt dispSize = 0;
+ dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);
+
+ PetscScalar* dispIncrCell = NULL;
+ PetscInt dispIncrSize = 0;
+ dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);assert(dispIncrCell);
assert(dispSize == dispIncrSize);
// Get cell geometry information that depends on cell
@@ -272,7 +265,7 @@
// Compute B(transpose) * sigma, first computing deformation
// tensor and strains
- _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispTpdtCell, numBasis, numQuadPts, spaceDim);
+ _calcDeformation(&deformCell, basisDeriv, coordsCell, &dispTpdtCell[0], numBasis, numQuadPts, spaceDim);
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
const scalar_array& stressCell = _material->calcStress(strainCell, true);
@@ -286,6 +279,8 @@
#endif
// Assemble cell contribution into field
residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), cell, ADD_VALUES);
+
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
} // for
_logger->eventEnd(computeEvent);
@@ -398,15 +393,10 @@
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
// Compute geometry information for current cell
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(cell);
-#else
coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
_quadrature->computeGeometry(coordsCell, coordsSize, cell);
- for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coordsCell[i];} // :TODO: Remove copy.
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
-#endif
+
// Get state variables for cell.
_material->retrievePropsAndVars(cell);
@@ -432,7 +422,7 @@
dispIncrVisitor.restoreClosure(&dispIncrCell, &dispIncrSize, cell);
// Compute deformation tensor, strains, and stresses
- _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispTpdtCell, numBasis, numQuadPts, spaceDim);
+ _calcDeformation(&deformCell, basisDeriv, coordsCell, &dispTpdtCell[0], numBasis, numQuadPts, spaceDim);
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
// Get "elasticity" matrix at quadrature points for this cell
@@ -477,6 +467,8 @@
// Assemble cell contribution into PETSc matrix.
jacobianVisitor.setClosure(&_cellMatrix[0], _cellMatrix.size(), cell, ADD_VALUES);
+
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
} // for
_logger->eventEnd(computeEvent);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -169,16 +169,22 @@
void
pylith::feassemble::GeometryHex3D::jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const
{ // jacobian
assert(jacobian);
assert(det);
+ assert(vertices);
+ assert(location);
- assert( (numCorners()*spaceDim() == vertices.size()) || // linear hex
- ((numCorners()+19)*spaceDim() == vertices.size()) ); // quadratic hex
- assert(cellDim() == location.size());
- assert(spaceDim()*cellDim() == jacobian->size());
+ assert(this->numCorners() == numVertices || // linear hex
+ this->numCorners()+1 == numVertices); // quadratic hex
+ assert(this->spaceDim() == spaceDim);
+ assert(this->cellDim() == cellDim);
+ assert(spaceDim*cellDim == jacobian->size());
const PylithScalar x0 = vertices[0];
const PylithScalar y0 = vertices[1];
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryHex3D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -91,12 +91,18 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const;
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
/** Compute Jacobian at location in cell.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -72,9 +72,9 @@
const int dim,
const int npts) const
{ // ptsRefToGlobal
- assert(0 != ptsGlobal);
- assert(0 != ptsRef);
- assert(0 != vertices);
+ assert(ptsGlobal);
+ assert(ptsRef);
+ assert(vertices);
assert(1 == dim);
assert(spaceDim() == dim);
@@ -94,15 +94,22 @@
void
pylith::feassemble::GeometryLine1D::jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
- assert( (numCorners()*spaceDim() == vertices.size()) || // linear edge
- ((numCorners()+1)*spaceDim() == vertices.size()) ); // quadratic edge
- assert(spaceDim()*cellDim() == jacobian->size());
+ assert(this->numCorners() == numVertices || // linear
+ this->numCorners()+1 == numVertices); // quadratic
+ assert(this->spaceDim() == spaceDim);
+ assert(this->cellDim() == cellDim);
+ assert(spaceDim*cellDim == jacobian->size());
const PylithScalar x0 = vertices[0];
const PylithScalar x1 = vertices[1];
@@ -121,10 +128,10 @@
const int dim,
const int npts) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
- assert(0 != vertices);
- assert(0 != location);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
assert(1 == dim);
assert(spaceDim() == dim);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine1D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -75,12 +75,18 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const;
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
/** Compute Jacobian at location in cell.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -72,9 +72,9 @@
const int dim,
const int npts) const
{ // ptsRefToGlobal
- assert(0 != ptsGlobal);
- assert(0 != ptsRef);
- assert(0 != vertices);
+ assert(ptsGlobal);
+ assert(ptsRef);
+ assert(vertices);
assert(2 == dim);
assert(spaceDim() == dim);
@@ -101,15 +101,22 @@
void
pylith::feassemble::GeometryLine2D::jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
- assert( (numCorners()*spaceDim() == vertices.size()) || // linear
- ((numCorners()+1)*spaceDim() == vertices.size()) ); // quadratic
- assert(spaceDim()*cellDim() == jacobian->size());
+ assert(this->numCorners() == numVertices || // linear
+ this->numCorners()+1 == numVertices); // quadratic
+ assert(this->spaceDim() == spaceDim);
+ assert(this->cellDim() == cellDim);
+ assert(spaceDim*cellDim == jacobian->size());
const PylithScalar x0 = vertices[0];
const PylithScalar y0 = vertices[1];
@@ -135,10 +142,10 @@
const int dim,
const int npts) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
- assert(0 != vertices);
- assert(0 != location);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
assert(2 == dim);
assert(spaceDim() == dim);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine2D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -75,12 +75,18 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const;
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
/** Compute Jacobian at location in cell.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -72,9 +72,9 @@
const int dim,
const int npts) const
{ // ptsRefToGlobal
- assert(0 != ptsGlobal);
- assert(0 != ptsRef);
- assert(0 != vertices);
+ assert(ptsGlobal);
+ assert(ptsRef);
+ assert(vertices);
assert(3 == dim);
assert(spaceDim() == dim);
@@ -104,16 +104,23 @@
// Compute Jacobian at location in cell.
void
pylith::feassemble::GeometryLine3D::jacobian(scalar_array* jacobian,
- PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const
+ PylithScalar* det,
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
- assert( (numCorners()*spaceDim() == vertices.size()) || // linear edge
- ((numCorners()+1)*spaceDim() == vertices.size()) ); // quadratic edge
- assert(spaceDim()*cellDim() == jacobian->size());
+ assert(this->numCorners() == numVertices || // linear
+ this->numCorners()+1 == numVertices); // quadratic
+ assert(this->spaceDim() == spaceDim);
+ assert(this->cellDim() == cellDim);
+ assert(spaceDim*cellDim == jacobian->size());
const PylithScalar x0 = vertices[0];
const PylithScalar y0 = vertices[1];
@@ -142,10 +149,10 @@
const int dim,
const int npts) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
- assert(0 != vertices);
- assert(0 != location);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
assert(3 == dim);
assert(spaceDim() == dim);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryLine3D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -75,12 +75,18 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const;
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
/** Compute Jacobian at location in cell.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -65,9 +65,9 @@
const int dim,
const int npts) const
{ // ptsRefToGlobal
- assert(0 != ptsGlobal);
- assert(0 != ptsRef);
- assert(0 != vertices);
+ assert(ptsGlobal);
+ assert(ptsRef);
+ assert(vertices);
assert(1 == dim);
assert(spaceDim() == dim);
@@ -79,12 +79,15 @@
// Compute Jacobian at location in cell.
void
pylith::feassemble::GeometryPoint1D::jacobian(scalar_array* jacobian,
- PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const
+ PylithScalar* det,
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
+ assert(jacobian);
+ assert(det);
assert(1 == jacobian->size());
@@ -102,10 +105,10 @@
const int dim,
const int npts) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
- assert(0 != vertices);
- assert(0 != ptsRef);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(ptsRef);
assert(1 == dim);
assert(spaceDim() == dim);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint1D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -75,12 +75,18 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const;
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
/** Compute Jacobian at location in cell.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -65,9 +65,9 @@
const int dim,
const int npts) const
{ // ptsRefToGlobal
- assert(0 != ptsGlobal);
- assert(0 != ptsRef);
- assert(0 != vertices);
+ assert(ptsGlobal);
+ assert(ptsRef);
+ assert(vertices);
assert(2 == dim);
assert(spaceDim() == dim);
@@ -80,12 +80,15 @@
// Compute Jacobian at location in cell.
void
pylith::feassemble::GeometryPoint2D::jacobian(scalar_array* jacobian,
- PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const
+ PylithScalar* det,
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
+ assert(jacobian);
+ assert(det);
assert(1 == jacobian->size());
@@ -103,10 +106,10 @@
const int dim,
const int npts) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
- assert(0 != vertices);
- assert(0 != location);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
assert(2 == dim);
assert(spaceDim() == dim);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint2D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -75,12 +75,18 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const;
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
/** Compute Jacobian at location in cell.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -65,9 +65,9 @@
const int dim,
const int npts) const
{ // ptsRefToGlobal
- assert(0 != ptsGlobal);
- assert(0 != ptsRef);
- assert(0 != vertices);
+ assert(ptsGlobal);
+ assert(ptsRef);
+ assert(vertices);
assert(3 == dim);
assert(spaceDim() == dim);
@@ -80,12 +80,15 @@
// Compute Jacobian at location in cell.
void
pylith::feassemble::GeometryPoint3D::jacobian(scalar_array* jacobian,
- PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const
+ PylithScalar* det,
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
+ assert(jacobian);
+ assert(det);
assert(1 == jacobian->size());
@@ -103,10 +106,10 @@
const int dim,
const int npts) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
- assert(0 != vertices);
- assert(0 != location);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
assert(3 == dim);
assert(spaceDim() == dim);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryPoint3D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -75,12 +75,18 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const;
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
/** Compute Jacobian at location in cell.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -74,9 +74,9 @@
const int dim,
const int npts) const
{ // ptsRefToGlobal
- assert(0 != ptsGlobal);
- assert(0 != ptsRef);
- assert(0 != vertices);
+ assert(ptsGlobal);
+ assert(ptsRef);
+ assert(vertices);
assert(2 == dim);
assert(spaceDim() == dim);
@@ -115,17 +115,23 @@
// Compute Jacobian at location in cell.
void
pylith::feassemble::GeometryQuad2D::jacobian(scalar_array* jacobian,
- PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const
+ PylithScalar* det,
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
- assert( (numCorners()*spaceDim() == vertices.size()) || // linear quad
- ((numCorners()+5)*spaceDim() == vertices.size()) ); // quadratic quad
- assert(cellDim() == location.size());
- assert(spaceDim()*cellDim() == jacobian->size());
+ assert(this->numCorners() == numVertices || // linear
+ this->numCorners()+1 == numVertices); // quadratic
+ assert(this->spaceDim() == spaceDim);
+ assert(this->cellDim() == cellDim);
+ assert(spaceDim*cellDim == jacobian->size());
const PylithScalar x0 = vertices[0];
const PylithScalar y0 = vertices[1];
@@ -169,10 +175,10 @@
const int dim,
const int npts) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
- assert(0 != vertices);
- assert(0 != ptsRef);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(ptsRef);
assert(2 == dim);
assert(spaceDim() == dim);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad2D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -93,12 +93,18 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const;
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
/** Compute Jacobian at location in cell.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -74,9 +74,9 @@
const int dim,
const int npts) const
{ // ptsRefToGlobal
- assert(0 != ptsGlobal);
- assert(0 != ptsRef);
- assert(0 != vertices);
+ assert(ptsGlobal);
+ assert(ptsRef);
+ assert(vertices);
assert(3 == dim);
assert(spaceDim() == dim);
@@ -125,17 +125,23 @@
// Compute Jacobian at location in cell.
void
pylith::feassemble::GeometryQuad3D::jacobian(scalar_array* jacobian,
- PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const
+ PylithScalar* det,
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
- assert( (numCorners()*spaceDim() == vertices.size()) || // linear quad
- ((numCorners()+5)*spaceDim() == vertices.size()) ); // quadratic quad
- assert(cellDim() == location.size());
- assert(spaceDim()*cellDim() == jacobian->size());
+ assert(this->numCorners() == numVertices || // linear
+ this->numCorners()+1 == numVertices); // quadratic
+ assert(this->spaceDim() == spaceDim);
+ assert(this->cellDim() == cellDim);
+ assert(spaceDim*cellDim == jacobian->size());
const PylithScalar x0 = vertices[0];
const PylithScalar y0 = vertices[1];
@@ -199,10 +205,10 @@
const int dim,
const int npts) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
- assert(0 != vertices);
- assert(0 != ptsRef);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(ptsRef);
assert(3 == dim);
assert(spaceDim() == dim);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryQuad3D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -91,12 +91,18 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const;
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
/** Compute Jacobian at location in cell.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -74,9 +74,9 @@
const int dim,
const int npts) const
{ // ptsRefToGlobal
- assert(0 != ptsGlobal);
- assert(0 != ptsRef);
- assert(0 != vertices);
+ assert(ptsGlobal);
+ assert(ptsRef);
+ assert(vertices);
assert(3 == dim);
assert(spaceDim() == dim);
@@ -129,14 +129,22 @@
void
pylith::feassemble::GeometryTet3D::jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
- assert(numCorners()*spaceDim() == vertices.size());
- assert(spaceDim()*cellDim() == jacobian->size());
+ assert(this->numCorners() == numVertices || // linear
+ this->numCorners()+1 == numVertices); // quadratic
+ assert(this->spaceDim() == spaceDim);
+ assert(this->cellDim() == cellDim);
+ assert(spaceDim*cellDim == jacobian->size());
const PylithScalar x0 = vertices[0];
const PylithScalar y0 = vertices[1];
@@ -185,10 +193,10 @@
const int dim,
const int npts) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
- assert(0 != vertices);
- assert(0 != ptsRef);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(ptsRef);
assert(3 == dim);
assert(spaceDim() == dim);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTet3D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -86,12 +86,18 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const;
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
/** Compute Jacobian at location in cell.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -73,9 +73,9 @@
const int dim,
const int npts) const
{ // ptsRefToGlobal
- assert(0 != ptsGlobal);
- assert(0 != ptsRef);
- assert(0 != vertices);
+ assert(ptsGlobal);
+ assert(ptsRef);
+ assert(vertices);
assert(2 == dim);
assert(spaceDim() == dim);
@@ -109,13 +109,22 @@
void
pylith::feassemble::GeometryTri2D::jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const
{ // jacobian
- assert(0 != jacobian);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
- assert(numCorners()*spaceDim() == vertices.size());
- assert(spaceDim()*cellDim() == jacobian->size());
+ assert(this->numCorners() == numVertices || // linear
+ this->numCorners()+1 == numVertices); // quadratic
+ assert(this->spaceDim() == spaceDim);
+ assert(this->cellDim() == cellDim);
+ assert(spaceDim*cellDim == jacobian->size());
const PylithScalar x0 = vertices[0];
const PylithScalar y0 = vertices[1];
@@ -148,10 +157,10 @@
const int dim,
const int npts) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
- assert(0 != vertices);
- assert(0 != location);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
assert(2 == dim);
assert(spaceDim() == dim);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri2D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -92,12 +92,18 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const;
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
/** Compute Jacobian at location in cell.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -73,9 +73,9 @@
const int dim,
const int npts) const
{ // ptsRefToGlobal
- assert(0 != ptsGlobal);
- assert(0 != ptsRef);
- assert(0 != vertices);
+ assert(ptsGlobal);
+ assert(ptsRef);
+ assert(vertices);
assert(3 == dim);
assert(spaceDim() == dim);
@@ -115,13 +115,22 @@
void
pylith::feassemble::GeometryTri3D::jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const
{ // jacobian
- assert(0 != jacobian);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
- assert(numCorners()*spaceDim() == vertices.size());
- assert(spaceDim()*cellDim() == jacobian->size());
+ assert(this->numCorners() == numVertices || // linear
+ this->numCorners()+1 == numVertices); // quadratic
+ assert(this->spaceDim() == spaceDim);
+ assert(this->cellDim() == cellDim);
+ assert(spaceDim*cellDim == jacobian->size());
const PylithScalar x0 = vertices[0];
const PylithScalar y0 = vertices[1];
@@ -171,10 +180,10 @@
const int dim,
const int npts) const
{ // jacobian
- assert(0 != jacobian);
- assert(0 != det);
- assert(0 != vertices);
- assert(0 != location);
+ assert(jacobian);
+ assert(det);
+ assert(vertices);
+ assert(location);
assert(3 == dim);
assert(spaceDim() == dim);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/GeometryTri3D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -92,12 +92,18 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(scalar_array* jacobian,
PylithScalar* det,
- const scalar_array& vertices,
- const scalar_array& location) const;
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
/** Compute Jacobian at location in cell.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -36,7 +36,6 @@
#include "pylith/utils/array.hh" // USES scalar_array
#include "pylith/utils/EventLogger.hh" // USES EventLogger
-#include <cstring> // USES memcpy()
#include <strings.h> // USES strcasecmp()
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
@@ -186,40 +185,35 @@
const PetscInt numCells = materialIS.size();
// Setup visitors.
- scalar_array dispTCell(numBasis*spaceDim);
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
- PetscScalar* dispCell = NULL;
- PetscInt dispSize = 0;
-
topology::CoordsVisitor coordsVisitor(dmMesh);
- PetscScalar* coordsCell = NULL;
- PetscInt coordsSize = 0;
// Loop over cells
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
// Retrieve geometry information for current cell
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+ PetscScalar* coordsCell = NULL;
+ PetscInt coordsSize = 0;
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
_quadrature->computeGeometry(coordsCell, coordsSize, cell);
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ const scalar_array& basisDeriv = _quadrature->basisDeriv();
// Get physical properties and state variables for cell.
_material->retrievePropsAndVars(cell);
// Restrict input fields to cell
- dispVisitor.getClosure(&dispCell, &dispSize, cell);
- assert(numBasis*spaceDim == dispSize);
- for(PetscInt i = 0; i < dispSize; ++i) {dispTCell[i] = dispCell[i];} // :TODO: Remove copy
- dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
+ PetscScalar* dispCell = NULL;
+ PetscInt dispSize = 0;
+ dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
- // Get cell geometry information that depends on cell
- const scalar_array& basisDeriv = _quadrature->basisDeriv();
-
// Compute strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispTCell, numBasis, numQuadPts);
+ calcTotalStrainFn(&strainCell, basisDeriv, dispCell, numBasis, spaceDim, numQuadPts);
// Update material state
_material->updateStateVars(strainCell, cell);
+
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
} // for
PYLITH_METHOD_END;
@@ -539,17 +533,14 @@
_quadrature->computeGeometry(coordsCell, coordsSize, cell);
coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
- // Restrict input fields to cell
+ // Get cell geometry information that depends on cell
dispVisitor.getClosure(&dispCell, &dispSize, cell);
assert(numBasis*spaceDim == dispSize);
- for(PetscInt i = 0; i < dispSize; ++i) {dispCellTmp[i] = dispCell[i];} // :TODO: Remove copy
- dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
-
- // Get cell geometry information that depends on cell
const scalar_array& basisDeriv = _quadrature->basisDeriv();
// Compute strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispCellTmp, numBasis, numQuadPts);
+ calcTotalStrainFn(&strainCell, basisDeriv, dispCell, numBasis, spaceDim, numQuadPts);
+ dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
const PetscInt off = fieldVisitor.sectionOffset(cell);
assert(tensorCellSize == fieldVisitor.sectionDof(cell));
@@ -973,8 +964,9 @@
void
pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D(scalar_array* strain,
const scalar_array& basisDeriv,
- const scalar_array& disp,
+ const PylithScalar* disp,
const int numBasis,
+ const int spaceDim,
const int numQuadPts)
{ // calcTotalStrain1D
assert(strain);
@@ -982,7 +974,7 @@
const int dim = 1;
assert(basisDeriv.size() == numQuadPts*numBasis*dim);
- assert(disp.size() == numBasis*dim);
+ assert(dim == spaceDim);
(*strain) = 0.0;
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
@@ -995,8 +987,9 @@
void
pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D(scalar_array* strain,
const scalar_array& basisDeriv,
- const scalar_array& disp,
+ const PylithScalar* disp,
const int numBasis,
+ const int spaceDim,
const int numQuadPts)
{ // calcTotalStrain2D
assert(strain);
@@ -1005,18 +998,15 @@
const int strainSize = 3;
assert(basisDeriv.size() == numQuadPts*numBasis*dim);
- assert(disp.size() == numBasis*dim);
+ assert(dim == spaceDim);
(*strain) = 0.0;
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] +=
- 0.5 * (basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim ] +
- basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+1]);
+ (*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] += 0.5 * (basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim ] +
+ basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+1]);
} // for
} // calcTotalStrain2D
@@ -1024,8 +1014,9 @@
void
pylith::feassemble::IntegratorElasticity::_calcTotalStrain3D(scalar_array* strain,
const scalar_array& basisDeriv,
- const scalar_array& disp,
+ const PylithScalar* disp,
const int numBasis,
+ const int spaceDim,
const int numQuadPts)
{ // calcTotalStrain3D
assert(strain);
@@ -1034,26 +1025,20 @@
const int strainSize = 6;
assert(basisDeriv.size() == numQuadPts*numBasis*dim);
- assert(disp.size() == numBasis*dim);
+ assert(dim == spaceDim);
(*strain) = 0.0;
for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
for (int iBasis=0, iQ=iQuad*numBasis*dim; iBasis < numBasis; ++iBasis) {
- (*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]);
+ (*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/pylith/feassemble/IntegratorElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -51,8 +51,9 @@
typedef void (*totalStrain_fn_type)(scalar_array*,
const scalar_array&,
- const scalar_array&,
+ const PylithScalar*,
const int,
+ const int,
const int);
@@ -196,15 +197,16 @@
* @param strain Strain tensor at quadrature points.
* @param basisDeriv Derivatives of basis functions at quadrature points.
* @param disp Displacement at vertices of cell.
- * @param dimension Dimension of cell.
* @param numBasis Number of basis functions for cell.
+ * @param spaceDim Spatial dimension.
* @param numQuadPts Number of quadrature points.
*/
static
void _calcTotalStrain1D(scalar_array* strain,
const scalar_array& basisDeriv,
- const scalar_array& disp,
+ const PylithScalar* disp,
const int numBasis,
+ const int spaceDim,
const int numQuadPts);
/** Compute total strain in at quadrature points of a cell.
@@ -213,13 +215,15 @@
* @param basisDeriv Derivatives of basis functions at quadrature points.
* @param disp Displacement at vertices of cell.
* @param numBasis Number of basis functions for cell.
+ * @param spaceDim Spatial dimension.
* @param numQuadPts Number of quadrature points.
*/
static
void _calcTotalStrain2D(scalar_array* strain,
const scalar_array& basisDeriv,
- const scalar_array& disp,
+ const PylithScalar* disp,
const int numBasis,
+ const int spaceDim,
const int numQuadPts);
/** Compute total strain in at quadrature points of a cell.
@@ -228,20 +232,17 @@
* @param basisDeriv Derivatives of basis functions at quadrature points.
* @param disp Displacement at vertices of cell.
* @param numBasis Number of basis functions for cell.
+ * @param spaceDim Spatial dimension.
* @param numQuadPts Number of quadrature points.
*/
static
void _calcTotalStrain3D(scalar_array* strain,
const scalar_array& basisDeriv,
- const scalar_array& disp,
+ const PylithScalar* disp,
const int numBasis,
+ const int spaceDim,
const int numQuadPts);
-// PROTECTED TYPEDEFS ///////////////////////////////////////////////////
-protected :
-
- typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -36,7 +36,6 @@
#include "pylith/utils/array.hh" // USES scalar_array
#include "pylith/utils/EventLogger.hh" // USES EventLogger
-#include <cstring> // USES memcpy()
#include <strings.h> // USES strcasecmp()
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
@@ -126,26 +125,22 @@
// Retrieve geometry information for current cell
coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
_quadrature->computeGeometry(coordsCell, coordsSize, cell);
- for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coordsCell[i];} // :TODO: Remove copy
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ const scalar_array& basisDeriv = _quadrature->basisDeriv();
- // Restrict input fields to cell
dispVisitor.getClosure(&dispCell, &dispSize, cell);
assert(numBasis*spaceDim == dispSize);
- for(PetscInt i = 0; i < dispSize; ++i) {dispCellTmp[i] = dispCell[i];} // :TODO: Remove copy
- dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
-
- // Get cell geometry information that depends on cell
- const scalar_array& basisDeriv = _quadrature->basisDeriv();
// Compute deformation tensor.
- _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispCellTmp, numBasis, numQuadPts, spaceDim);
+ _calcDeformation(&deformCell, basisDeriv, coordsCell, dispCell, numBasis, numQuadPts, spaceDim);
// Compute strains
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
// Update material state
_material->updateStateVars(strainCell, cell);
+
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
} // for
PYLITH_METHOD_END;
@@ -185,7 +180,6 @@
} // else
// Allocate arrays for cell data.
- scalar_array dispCellTmp(numBasis*spaceDim);
scalar_array deformCell(numQuadPts*spaceDim*spaceDim);
const int tensorCellSize = numQuadPts*tensorSize;
scalar_array strainCell(tensorCellSize);
@@ -202,38 +196,34 @@
// Setup field visitors.
topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
- PetscScalar* dispCell = NULL;
- PetscInt dispSize = 0;
topology::VecVisitorMesh fieldVisitor(*field);
PetscScalar* fieldArray = fieldVisitor.localArray();
scalar_array coordinatesCell(numBasis*spaceDim);
topology::CoordsVisitor coordsVisitor(dmMesh);
- PetscScalar *coordsCell = NULL;
- PetscInt coordsSize = 0;
// Loop over cells
for (PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
// Retrieve geometry information for current cell
- coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+ PetscScalar *coordsCell = NULL;
+ PetscInt coordsSize = 0;
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
_quadrature->computeGeometry(coordsCell, coordsSize, cell);
- for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coordsCell[i];} // :TODO: Remove copy
- coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ const scalar_array& basisDeriv = _quadrature->basisDeriv();
// Restrict input fields to cell
- dispVisitor.getClosure(&dispCell, &dispSize, cell);
- assert(numBasis*spaceDim == dispSize);
- for(PetscInt i = 0; i < dispSize; ++i) {dispCellTmp[i] = dispCell[i];} // :TODO: Remove copy
- dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
+ PetscScalar* dispCell = NULL;
+ PetscInt dispSize = 0;
+ dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
- // Get cell geometry information that depends on cell
- const scalar_array& basisDeriv = _quadrature->basisDeriv();
-
// Compute deformation tensor.
- _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispCellTmp, numBasis, numQuadPts, spaceDim);
+ _calcDeformation(&deformCell, basisDeriv, coordsCell, dispCell, numBasis, numQuadPts, spaceDim);
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+ dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
+
// Compute strains
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
@@ -1197,8 +1187,8 @@
void
pylith::feassemble::IntegratorElasticityLgDeform::_calcDeformation(scalar_array* deform,
const scalar_array& basisDeriv,
- const scalar_array& vertices,
- const scalar_array& disp,
+ const PylithScalar* vertices,
+ const PylithScalar* disp,
const int numBasis,
const int numQuadPts,
const int dim)
@@ -1206,7 +1196,6 @@
assert(deform);
assert(basisDeriv.size() == numQuadPts*numBasis*dim);
- assert(disp.size() == numBasis*dim);
const int deformSize = dim*dim;
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -188,8 +188,8 @@
static
void _calcDeformation(scalar_array* deform,
const scalar_array& basisDeriv,
- const scalar_array& vertices,
- const scalar_array& disp,
+ const PylithScalar* vertices,
+ const PylithScalar* disp,
const int numBasis,
const int numQuadPts,
const int dim);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -29,7 +29,6 @@
#include "Quadrature2Din3D.hh"
#include "Quadrature3D.hh"
-#include <cstring> // USES memcpy()
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -113,14 +113,6 @@
/// Deallocate temporary storage.
void clear(void);
- /** Compute geometric quantities for cell.
- *
- * @param coordinatesCell Coordinates of vertices in cell.
- * @param cell Finite-element cell
- */
- 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.
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.icc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.icc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -75,16 +75,6 @@
return _engine->jacobianDet();
}
-// Compute geometric quantities for each cell.
-template<typename mesh_type>
-inline
-void
-pylith::feassemble::Quadrature<mesh_type>::computeGeometry(const scalar_array& coordinatesCell,
- const int cell) {
- assert(_engine);
- _engine->computeGeometry(coordinatesCell, cell);
-}
-
// Compute geometric quantities for a cell at quadrature points.
template<typename mesh_type>
void
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature0D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature0D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature0D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -50,38 +50,6 @@
// ----------------------------------------------------------------------
// Compute geometric quantities for a cell at quadrature points.
void
-pylith::feassemble::Quadrature0D::computeGeometry(const scalar_array& coordinatesCell,
- const int cell)
-{ // computeGeometry
- 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 scalar_array& basisDerivRef = _quadRefCell.basisDerivRef();
-
- zero();
-
- _quadPts = coordinatesCell;
-
- _jacobian[0] = 1.0;
- _jacobianDet[0] = 1.0;
- _jacobianInv[0] = 1.0;
- _basisDeriv[0] = basisDerivRef[0];
-
- PetscLogFlops(0);
-} // computeGeometry
-
-
-// ----------------------------------------------------------------------
-// Compute geometric quantities for a cell at quadrature points.
-void
pylith::feassemble::Quadrature0D::computeGeometry(const PylithScalar* coordinatesCell,
const int coordinatesSize,
const int cell)
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature0D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature0D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature0D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -55,14 +55,6 @@
/** Compute geometric quantities for a cell at quadrature points.
*
- * @param coordinatesCell Coordinates of cell's vertices.
- * @param cell Finite-element cell
- */
- 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
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -52,81 +52,6 @@
// ----------------------------------------------------------------------
// Compute geometric quantities for a cell at quadrature points.
void
-pylith::feassemble::Quadrature1D::computeGeometry(const scalar_array& coordinatesCell,
- const int cell)
-{ // computeGeometry
- 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 == coordinatesCell.size());
-
- 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[0], 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[0], &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
-
-
-// ----------------------------------------------------------------------
-// Compute geometric quantities for a cell at quadrature points.
-void
pylith::feassemble::Quadrature1D::computeGeometry(const PylithScalar* coordinatesCell,
const int coordinatesSize,
const int cell)
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -51,14 +51,6 @@
/** Compute geometric quantities for a cell at quadrature points.
*
- * @param coordinatesCell Coordinates of cell's vertices.
- * @param cell Finite-element cell
- */
- 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
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din2D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din2D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -52,106 +52,6 @@
// ----------------------------------------------------------------------
// Compute geometric quantities for a cell at quadrature points.
void
-pylith::feassemble::Quadrature1Din2D::computeGeometry(const scalar_array& coordinatesCell,
- const int cell)
-{ // computeGeometry
- 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 == coordinatesCell.size());
-
- 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[0], 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[0], &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
-
-
-// ----------------------------------------------------------------------
-// Compute geometric quantities for a cell at quadrature points.
-void
pylith::feassemble::Quadrature1Din2D::computeGeometry(const PylithScalar* coordinatesCell,
const int coordinatesSize,
const int cell)
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din2D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din2D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -51,14 +51,6 @@
/** Compute geometric quantities for a cell at quadrature points.
*
- * @param coordinatesCell Coordinates of cell's vertices.
- * @param cell Finite-element cell
- */
- 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
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din3D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din3D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -52,110 +52,6 @@
// ----------------------------------------------------------------------
// Compute geometric quantities for a cell at quadrature points.
void
-pylith::feassemble::Quadrature1Din3D::computeGeometry(const scalar_array& coordinatesCell,
- const int cell)
-{ // computeGeometry
- 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 == coordinatesCell.size());
-
- 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[0], 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[0], &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
-
-
-// ----------------------------------------------------------------------
-// Compute geometric quantities for a cell at quadrature points.
-void
pylith::feassemble::Quadrature1Din3D::computeGeometry(const PylithScalar* coordinatesCell,
const int coordinatesSize,
const int cell)
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din3D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature1Din3D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -52,14 +52,6 @@
/** Compute geometric quantities for a cell at quadrature points.
*
- * @param coordinatesCell Coordinates of cell's vertices.
- * @param cell Finite-element cell
- */
- 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
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -52,118 +52,6 @@
// ----------------------------------------------------------------------
// Compute geometric quantities for a cell at quadrature points.
void
-pylith::feassemble::Quadrature2D::computeGeometry(const scalar_array& coordinatesCell,
- const int cell)
-{ // computeGeometry
- 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 == coordinatesCell.size());
- 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[0], 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[0], &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
-
-
-// ----------------------------------------------------------------------
-// Compute geometric quantities for a cell at quadrature points.
-void
pylith::feassemble::Quadrature2D::computeGeometry(const PylithScalar* coordinatesCell,
const int coordinatesSize,
const int cell)
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -51,14 +51,6 @@
/** Compute geometric quantities for a cell at quadrature points.
*
- * @param coordinatesCell Coordinates of cell's vertices.
- * @param cell Finite-element cell
- */
- 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
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2Din3D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2Din3D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -54,196 +54,6 @@
// ----------------------------------------------------------------------
// Compute geometric quantities for a cell at quadrature points.
void
-pylith::feassemble::Quadrature2Din3D::computeGeometry(const scalar_array& coordinatesCell,
- const int cell)
-{ // computeGeometry
- 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 == coordinatesCell.size());
- 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[0], 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[0], &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
-
-
-// ----------------------------------------------------------------------
-// Compute geometric quantities for a cell at quadrature points.
-void
pylith::feassemble::Quadrature2Din3D::computeGeometry(const PylithScalar* coordinatesCell,
const int coordinatesSize,
const int cell)
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2Din3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2Din3D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature2Din3D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -51,14 +51,6 @@
/** Compute geometric quantities for a cell at quadrature points.
*
- * @param coordinatesCell Coordinates of cell's vertices.
- * @param cell Finite-element cell
- */
- 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
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature3D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature3D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -52,151 +52,6 @@
// ----------------------------------------------------------------------
// Compute geometric quantities for a cell at quadrature points.
void
-pylith::feassemble::Quadrature3D::computeGeometry(const scalar_array& coordinatesCell,
- const int cell)
-{ // computeGeometry
- 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 == coordinatesCell.size());
- 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[0], 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[0], &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
-
-
-// ----------------------------------------------------------------------
-// Compute geometric quantities for a cell at quadrature points.
-void
pylith::feassemble::Quadrature3D::computeGeometry(const PylithScalar* coordinatesCell,
const int coordinatesSize,
const int cell)
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature3D.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature3D.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -51,14 +51,6 @@
/** Compute geometric quantities for a cell at quadrature points.
*
- * @param coordinatesCell Coordinates of cell's vertices.
- * @param cell Finite-element cell
- */
- 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
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureEngine.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureEngine.hh 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureEngine.hh 2013-05-08 22:30:43 UTC (rev 22004)
@@ -92,15 +92,6 @@
/** Compute geometric quantities for a cell at quadrature points.
*
- * @param coordinatesCell Coordinates of cell's vertices.
- * @param cell Finite-element cell
- */
- virtual
- 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
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureRefCell.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureRefCell.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/QuadratureRefCell.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -24,7 +24,6 @@
#include "pylith/utils/error.h" // USES PYLITH_METHOD_BEGIN/END
-#include <cstring> // USES memcpy()
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -32,7 +32,6 @@
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
-#include <cstring> // USES memcpy()
#include <strings.h> // USES strcasecmp()
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -34,7 +34,6 @@
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
-#include <cstring> // USES memcpy()
#include <strings.h> // USES strcasecmp()
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellIsotropic3D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -31,7 +31,6 @@
#include "petsc.h" // USES PetscLogFlops
#include <cassert> // USES assert()
-#include <cstring> // USES memcpy()
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellPlaneStrain.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -31,7 +31,6 @@
#include "petsc.h" // USES PetscLogFlops
#include <cassert> // USES assert()
-#include <cstring> // USES memcpy()
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellQpQsIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellQpQsIsotropic3D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/GenMaxwellQpQsIsotropic3D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -31,7 +31,6 @@
#include "petsc.h" // USES PetscLogFlops
#include <cassert> // USES assert()
-#include <cstring> // USES memcpy()
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -31,14 +31,11 @@
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
-#include <cstring> // USES memcpy()
#include <strings.h> // USES strcasecmp()
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
-//#define PRECOMPUTE_GEOMETRY
-
// ----------------------------------------------------------------------
// Default constructor.
pylith::materials::Material::Material(const int dimension,
@@ -146,11 +143,9 @@
topology::VecVisitorMesh propertiesVisitor(*_properties);
PetscScalar* propertiesArray = propertiesVisitor.localArray();
-#if !defined(PRECOMPUTE_GEOMETRY)
topology::CoordsVisitor coordsVisitor(dmMesh);
PetscScalar* coordsArray = NULL;
PetscInt coordsSize = 0;
-#endif
// Create arrays for querying.
const int numDBProperties = _metadata.numDBProperties();
@@ -205,13 +200,9 @@
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
// Compute geometry information for current cell
-#if defined(PRECOMPUTE_GEOMETRY)
- quadrature->retrieveGeometry(cell);
-#else
coordsVisitor.getClosure(&coordsArray, &coordsSize, cell);
quadrature->computeGeometry(coordsArray, coordsSize, cell);
coordsVisitor.restoreClosure(&coordsArray, &coordsSize, cell);
-#endif
const scalar_array& quadPtsNonDim = quadrature->quadPts();
quadPtsGlobal = quadPtsNonDim;
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -1336,17 +1336,11 @@
// remove old scatter
ScatterInfo& sinfo = _scatters[context];
PetscErrorCode err = 0;
- if (sinfo.vector) {
- err = VecDestroy(&sinfo.vector);PYLITH_CHECK_ERROR(err);
- } // if
- if (sinfo.scatter) {
- err = VecScatterDestroy(&sinfo.scatter);PYLITH_CHECK_ERROR(err);
- } // if
+ err = DMDestroy(&sinfo.dm);PYLITH_CHECK_ERROR(err);
+ err = VecDestroy(&sinfo.vector);PYLITH_CHECK_ERROR(err);
+ err = VecScatterDestroy(&sinfo.scatter);PYLITH_CHECK_ERROR(err);
+ err = VecDestroy(&sinfo.scatterVec);PYLITH_CHECK_ERROR(err);
- if (sinfo.scatterVec) {
- err = VecDestroy(&sinfo.scatterVec);PYLITH_CHECK_ERROR(err);
- } // if
-
_scatters.erase(context);
isNewScatter = true;
} // if
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/CellGeometry.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/CellGeometry.i 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/CellGeometry.i 2013-05-08 22:30:43 UTC (rev 22004)
@@ -117,14 +117,20 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
virtual
void jacobian(pylith::scalar_array* jacobian,
PylithScalar* det,
- const pylith::scalar_array& vertices,
- const pylith::scalar_array& location) const = 0;
-
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const = 0;
+
/** Compute Jacobian at location in cell.
*
* @param jacobian Jacobian at location.
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/GeometryHex3D.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/GeometryHex3D.i 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/GeometryHex3D.i 2013-05-08 22:30:43 UTC (rev 22004)
@@ -68,13 +68,19 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(pylith::scalar_array* jacobian,
PylithScalar* det,
- const pylith::scalar_array& vertices,
- const pylith::scalar_array& location) const;
-
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
+
/** Compute Jacobian at location in cell.
*
* @param jacobian Jacobian at location.
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/GeometryLine1D.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/GeometryLine1D.i 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/GeometryLine1D.i 2013-05-08 22:30:43 UTC (rev 22004)
@@ -68,13 +68,19 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(pylith::scalar_array* jacobian,
PylithScalar* det,
- const pylith::scalar_array& vertices,
- const pylith::scalar_array& location) const;
-
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
+
/** Compute Jacobian at location in cell.
*
* @param jacobian Jacobian at location.
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/GeometryLine2D.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/GeometryLine2D.i 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/GeometryLine2D.i 2013-05-08 22:30:43 UTC (rev 22004)
@@ -68,13 +68,19 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(pylith::scalar_array* jacobian,
PylithScalar* det,
- const pylith::scalar_array& vertices,
- const pylith::scalar_array& location) const;
-
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
+
/** Compute Jacobian at location in cell.
*
* @param jacobian Jacobian at location.
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/GeometryLine3D.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/GeometryLine3D.i 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/GeometryLine3D.i 2013-05-08 22:30:43 UTC (rev 22004)
@@ -68,13 +68,19 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(pylith::scalar_array* jacobian,
PylithScalar* det,
- const pylith::scalar_array& vertices,
- const pylith::scalar_array& location) const;
-
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
+
/** Compute Jacobian at location in cell.
*
* @param jacobian Jacobian at location.
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/GeometryPoint1D.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/GeometryPoint1D.i 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/GeometryPoint1D.i 2013-05-08 22:30:43 UTC (rev 22004)
@@ -68,13 +68,19 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(pylith::scalar_array* jacobian,
PylithScalar* det,
- const pylith::scalar_array& vertices,
- const pylith::scalar_array& location) const;
-
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
+
/** Compute Jacobian at location in cell.
*
* @param jacobian Jacobian at location.
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/GeometryPoint2D.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/GeometryPoint2D.i 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/GeometryPoint2D.i 2013-05-08 22:30:43 UTC (rev 22004)
@@ -68,13 +68,19 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(pylith::scalar_array* jacobian,
PylithScalar* det,
- const pylith::scalar_array& vertices,
- const pylith::scalar_array& location) const;
-
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
+
/** Compute Jacobian at location in cell.
*
* @param jacobian Jacobian at location.
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/GeometryPoint3D.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/GeometryPoint3D.i 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/GeometryPoint3D.i 2013-05-08 22:30:43 UTC (rev 22004)
@@ -68,13 +68,19 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(pylith::scalar_array* jacobian,
PylithScalar* det,
- const pylith::scalar_array& vertices,
- const pylith::scalar_array& location) const;
-
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
+
/** Compute Jacobian at location in cell.
*
* @param jacobian Jacobian at location.
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/GeometryQuad2D.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/GeometryQuad2D.i 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/GeometryQuad2D.i 2013-05-08 22:30:43 UTC (rev 22004)
@@ -68,13 +68,19 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(pylith::scalar_array* jacobian,
PylithScalar* det,
- const pylith::scalar_array& vertices,
- const pylith::scalar_array& location) const;
-
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
+
/** Compute Jacobian at location in cell.
*
* @param jacobian Jacobian at location.
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/GeometryQuad3D.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/GeometryQuad3D.i 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/GeometryQuad3D.i 2013-05-08 22:30:43 UTC (rev 22004)
@@ -68,13 +68,19 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(pylith::scalar_array* jacobian,
PylithScalar* det,
- const pylith::scalar_array& vertices,
- const pylith::scalar_array& location) const;
-
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
+
/** Compute Jacobian at location in cell.
*
* @param jacobian Jacobian at location.
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/GeometryTet3D.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/GeometryTet3D.i 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/GeometryTet3D.i 2013-05-08 22:30:43 UTC (rev 22004)
@@ -68,13 +68,19 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(pylith::scalar_array* jacobian,
PylithScalar* det,
- const pylith::scalar_array& vertices,
- const pylith::scalar_array& location) const;
-
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
+
/** Compute Jacobian at location in cell.
*
* @param jacobian Jacobian at location.
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/GeometryTri2D.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/GeometryTri2D.i 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/GeometryTri2D.i 2013-05-08 22:30:43 UTC (rev 22004)
@@ -68,13 +68,19 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(pylith::scalar_array* jacobian,
PylithScalar* det,
- const pylith::scalar_array& vertices,
- const pylith::scalar_array& location) const;
-
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
+
/** Compute Jacobian at location in cell.
*
* @param jacobian Jacobian at location.
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/GeometryTri3D.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/GeometryTri3D.i 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/GeometryTri3D.i 2013-05-08 22:30:43 UTC (rev 22004)
@@ -68,13 +68,19 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
+ * @param numVertices Number of vertices in cell.
+ * @param spaceDim Spatial dimension of coordinates.
* @param location Location in reference cell at which to compute Jacobian.
+ * @param cellDim Dimension of reference cell.
*/
void jacobian(pylith::scalar_array* jacobian,
PylithScalar* det,
- const pylith::scalar_array& vertices,
- const pylith::scalar_array& location) const;
-
+ const PylithScalar* vertices,
+ const int numVertices,
+ const int spaceDim,
+ const PylithScalar* location,
+ const int cellDim) const;
+
/** Compute Jacobian at location in cell.
*
* @param jacobian Jacobian at location.
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestCellGeometry.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestCellGeometry.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestCellGeometry.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -283,11 +283,8 @@
scalar_array jacobian(cellDim*spaceDim);
PylithScalar det = 0;
for (int iLoc=0; iLoc < numLocs; ++iLoc) {
- scalar_array vertices(_data->vertices, numCorners*spaceDim);
- scalar_array location(&_data->locations[iLoc*cellDim], cellDim);
+ _object->jacobian(&jacobian, &det, _data->vertices, numCorners, spaceDim, &_data->locations[iLoc*cellDim], cellDim);
- _object->jacobian(&jacobian, &det, vertices, location);
-
const int size = jacobian.size();
const int index = iLoc*cellDim*spaceDim;
const PylithScalar tolerance = 1.0e-06;
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestGeometryPoint1D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestGeometryPoint1D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestGeometryPoint1D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -94,10 +94,7 @@
scalar_array jacobian(1);
PylithScalar det = 0.0;
for (int iLoc=0; iLoc < numLocs; ++iLoc) {
- scalar_array vertices(data.vertices, numCorners*spaceDim);
- scalar_array location(&data.locations[iLoc], 1);
-
- geometry.jacobian(&jacobian, &det, vertices, location);
+ geometry.jacobian(&jacobian, &det, data.vertices, numCorners, spaceDim, &data.locations[iLoc], 1);
CPPUNIT_ASSERT_EQUAL(PylithScalar(1.0), jacobian[0]);
CPPUNIT_ASSERT_EQUAL(PylithScalar(1.0), det);
} //for
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestGeometryPoint2D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestGeometryPoint2D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestGeometryPoint2D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -94,10 +94,7 @@
scalar_array jacobian(1);
PylithScalar det = 0.0;
for (int iLoc=0; iLoc < numLocs; ++iLoc) {
- scalar_array vertices(data.vertices, numCorners*spaceDim);
- scalar_array location(&data.locations[iLoc], 1);
-
- geometry.jacobian(&jacobian, &det, vertices, location);
+ geometry.jacobian(&jacobian, &det, data.vertices, numCorners, spaceDim, &data.locations[iLoc], 1);
CPPUNIT_ASSERT_EQUAL(PylithScalar(1.0), jacobian[0]);
CPPUNIT_ASSERT_EQUAL(PylithScalar(1.0), det);
} //for
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestGeometryPoint3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestGeometryPoint3D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestGeometryPoint3D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -97,7 +97,7 @@
scalar_array vertices(data.vertices, numCorners*spaceDim);
scalar_array location(&data.locations[iLoc], 1);
- geometry.jacobian(&jacobian, &det, vertices, location);
+ geometry.jacobian(&jacobian, &det, data.vertices, numCorners, spaceDim, &data.locations[iLoc], 1);
CPPUNIT_ASSERT_EQUAL(PylithScalar(1.0), jacobian[0]);
CPPUNIT_ASSERT_EQUAL(PylithScalar(1.0), det);
} //for
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticity.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticity.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -43,22 +43,21 @@
const int dim = 1;
const int numBasis = 2;
const int numQuadPts = 2;
- const PylithScalar basisDerivVals[] = {
+ const PylithScalar basisDerivVals[numQuadPts*numBasis*dim] = {
-0.50, 0.50,
-0.25, 0.25 };
const int tensorSize = 1;
// Let u(x) = 1 + 0.5 * x
- const PylithScalar dispVals[] = { 0.5, 1.5 };
- const PylithScalar strainE[] = { 0.5, 0.25 };
+ const PylithScalar disp[numBasis*dim] = { 0.5, 1.5 };
+ const PylithScalar strainE[numQuadPts*tensorSize] = { 0.5, 0.25 };
const int size = numQuadPts * tensorSize;
scalar_array strain(size);
scalar_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
- scalar_array disp(dispVals, numBasis*dim);
- IntegratorElasticity::_calcTotalStrain1D(&strain, basisDeriv, disp, numBasis, numQuadPts);
+ IntegratorElasticity::_calcTotalStrain1D(&strain, basisDeriv, disp, numBasis, dim, numQuadPts);
const PylithScalar tolerance = 1.0e-06;
CPPUNIT_ASSERT_EQUAL(size, int(strain.size()));
@@ -85,7 +84,7 @@
const int dim = 2;
const int numBasis = 3;
const int numQuadPts = 2;
- const PylithScalar basisDerivVals[] = {
+ const PylithScalar basisDerivVals[numQuadPts*numBasis*dim] = {
+1.0, 0.0, 0.0, +1.0, -1.0, -1.0,
+2.0, 0.0, 0.0, +2.0, -2.0, -2.0
};
@@ -93,12 +92,12 @@
// Let ux(x,y) = +0.4 + 0.3*x + 0.8*y
// Ley uy(x,y) = -2.0 + 0.5*x - 0.2*y
- const PylithScalar dispVals[] = {
+ const PylithScalar disp[numBasis*dim] = {
0.7, -1.5,
1.2, -2.2,
0.4, -2.0
};
- const PylithScalar strainE[] = {
+ const PylithScalar strainE[numQuadPts*tensorSize] = {
0.3, -0.2, 0.65,
0.6, -0.4, 1.3
};
@@ -107,9 +106,8 @@
scalar_array strain(size);
scalar_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
- scalar_array disp(dispVals, numBasis*dim);
- IntegratorElasticity::_calcTotalStrain2D(&strain, basisDeriv, disp, numBasis, numQuadPts);
+ IntegratorElasticity::_calcTotalStrain2D(&strain, basisDeriv, disp, numBasis, dim, numQuadPts);
const PylithScalar tolerance = 1.0e-06;
CPPUNIT_ASSERT_EQUAL(size, int(strain.size()));
@@ -138,7 +136,7 @@
const int dim = 3;
const int numBasis = 4;
const int numQuadPts = 2;
- const PylithScalar basisDerivVals[] = {
+ const PylithScalar basisDerivVals[numQuadPts*numBasis*dim] = {
+1.0, 0.0, 0.0, // Quad pt 0
0.0, +1.0, 0.0,
0.0, 0.0, +1.0,
@@ -153,13 +151,13 @@
// Let ux(x,y,z) = +0.4 + 0.3*x + 0.8*y + 0.4*z
// Ley uy(x,y,z) = -2.0 + 0.5*x - 0.2*y + 1.2*z
// Ley uz(x,y,z) = -1.0 + 0.2*x - 0.7*y - 0.3*z
- const PylithScalar dispVals[] = {
+ const PylithScalar disp[numBasis*dim] = {
0.7, -1.5, -0.8,
1.2, -2.2, -1.7,
0.8, -0.8, -1.3,
0.4, -2.0, -1.0
};
- const PylithScalar strainE[] = {
+ const PylithScalar strainE[numQuadPts*tensorSize] = {
0.3, -0.2, -0.3, 0.65, 0.25, 0.3,
0.9, -0.6, -0.9, 1.95, 0.75, 0.9
};
@@ -168,9 +166,8 @@
scalar_array strain(size);
scalar_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
- scalar_array disp(dispVals, numBasis*dim);
- IntegratorElasticity::_calcTotalStrain3D(&strain, basisDeriv, disp, numBasis, numQuadPts);
+ IntegratorElasticity::_calcTotalStrain3D(&strain, basisDeriv, disp, numBasis, dim, numQuadPts);
const PylithScalar tolerance = 1.0e-06;
CPPUNIT_ASSERT_EQUAL(size, int(strain.size()));
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticityLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticityLgDeform.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticityLgDeform.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -43,8 +43,8 @@
const int dim = 1;
const int numBasis = 2;
const int numQuadPts = 2;
- const PylithScalar verticesVals[] = { -1.0, 1.0 };
- const PylithScalar basisDerivVals[] = {
+ const PylithScalar vertices[numBasis*dim] = { -1.0, 1.0 };
+ const PylithScalar basisDerivVals[numQuadPts*numBasis*dim] = {
-0.50, 0.50,
-0.50, 0.50 };
const int tensorSize = 1;
@@ -52,17 +52,14 @@
const int size = numQuadPts * dim*dim;
scalar_array deform(size);
scalar_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
- scalar_array vertices(verticesVals, numBasis*dim);
const PylithScalar tolerance = 1.0e-06;
{ // Rigid body translation
// u(x) = 1.0
- const PylithScalar dispVals[] = { 1.0, 1.0 };
- const PylithScalar deformE[] = { 1.0, 1.0 };
+ const PylithScalar disp[numBasis*dim] = { 1.0, 1.0 };
+ const PylithScalar deformE[size] = { 1.0, 1.0 };
- scalar_array disp(dispVals, numBasis*dim);
-
IntegratorElasticityLgDeform::_calcDeformation(&deform, basisDeriv, vertices, disp, numBasis, numQuadPts, dim);
CPPUNIT_ASSERT_EQUAL(size, int(deform.size()));
@@ -72,11 +69,9 @@
{ // Uniform strain
// u(x) = 0.1*x
- const PylithScalar dispVals[] = { -0.1, 0.1 };
- const PylithScalar deformE[] = { 1.1, 1.1 };
+ const PylithScalar disp[numBasis*dim] = { -0.1, 0.1 };
+ const PylithScalar deformE[size] = { 1.1, 1.1 };
- scalar_array disp(dispVals, numBasis*dim);
-
IntegratorElasticityLgDeform::_calcDeformation(&deform, basisDeriv, vertices, disp, numBasis, numQuadPts, dim);
CPPUNIT_ASSERT_EQUAL(size, int(deform.size()));
@@ -103,12 +98,12 @@
const int dim = 2;
const int numBasis = 3;
const int numQuadPts = 2;
- const PylithScalar verticesVals[] = {
+ const PylithScalar vertices[numBasis*dim] = {
1.0, 0.0,
0.0, 1.0,
0.0, 0.0,
};
- const PylithScalar basisDerivVals[] = {
+ const PylithScalar basisDerivVals[numQuadPts*numBasis*dim] = {
+1.0, 0.0, 0.0, +1.0, -1.0, -1.0,
+1.0, 0.0, 0.0, +1.0, -1.0, -1.0
};
@@ -117,25 +112,22 @@
const int size = numQuadPts * dim*dim;
scalar_array deform(size);
scalar_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
- scalar_array vertices(verticesVals, numBasis*dim);
const PylithScalar tolerance = 1.0e-06;
{ // Rigid body translation
// ux(x,y) = 0.5
// uy(x,y) = 0.2
- const PylithScalar dispVals[] = {
+ const PylithScalar disp[numBasis*dim] = {
0.5, 0.2,
0.5, 0.2,
0.5, 0.2,
};
- const PylithScalar deformE[] = {
+ const PylithScalar deformE[size] = {
1.0, 0.0, 0.0, 1.0,
1.0, 0.0, 0.0, 1.0,
};
- scalar_array disp(dispVals, numBasis*dim);
-
IntegratorElasticityLgDeform::_calcDeformation(&deform, basisDeriv, vertices, disp, numBasis, numQuadPts, dim);
CPPUNIT_ASSERT_EQUAL(size, int(deform.size()));
@@ -149,7 +141,7 @@
// theta = pi/6
const PylithScalar pi = 4.0*atan(1.0);
const PylithScalar theta = pi / 6.0;
- const PylithScalar dispVals[] = {
+ const PylithScalar disp[numBasis*dim] = {
0.5+cos(theta)*1.0+sin(theta)*0.0-1.0,
0.2-sin(theta)*1.0+cos(theta)*0.0-0.0,
@@ -159,13 +151,11 @@
0.5+cos(theta)*0.0+sin(theta)*0.0-0.0,
0.2-sin(theta)*0.0+cos(theta)*0.0-0.0,
};
- const PylithScalar deformE[] = {
+ const PylithScalar deformE[size] = {
cos(theta), sin(theta), -sin(theta), cos(theta),
cos(theta), sin(theta), -sin(theta), cos(theta),
};
- scalar_array disp(dispVals, numBasis*dim);
-
IntegratorElasticityLgDeform::_calcDeformation(&deform, basisDeriv, vertices, disp, numBasis, numQuadPts, dim);
CPPUNIT_ASSERT_EQUAL(size, int(deform.size()));
@@ -176,18 +166,16 @@
{ // Uniform strain
// Let ux(x,y) = +0.4 + 0.3*x + 0.8*y
// Ley uy(x,y) = -2.0 + 0.5*x - 0.2*y
- const PylithScalar dispVals[] = {
+ const PylithScalar disp[numBasis*dim] = {
0.7, -1.5,
1.2, -2.2,
0.4, -2.0
};
- const PylithScalar deformE[] = {
+ const PylithScalar deformE[size] = {
1.3, 0.8, 0.5, 0.8,
1.3, 0.8, 0.5, 0.8,
};
- scalar_array disp(dispVals, numBasis*dim);
-
IntegratorElasticityLgDeform::_calcDeformation(&deform, basisDeriv, vertices, disp, numBasis, numQuadPts, dim);
CPPUNIT_ASSERT_EQUAL(size, int(deform.size()));
@@ -224,7 +212,7 @@
const int dim = 3;
const int numBasis = 8;
const int numQuadPts = 1;
- const PylithScalar verticesVals[] = {
+ const PylithScalar vertices[numBasis*dim] = {
-1.0, -1.0, -1.0,
+1.0, -1.0, -1.0,
+1.0, +1.0, -1.0,
@@ -234,7 +222,7 @@
+1.0, +1.0, +1.0,
-1.0, +1.0, +1.0,
};
- const PylithScalar basisDerivVals[] = {
+ const PylithScalar basisDerivVals[numQuadPts*numBasis*dim] = {
-0.125, -0.125, -0.125,
+0.125, -0.125, -0.125,
+0.125, +0.125, -0.125,
@@ -249,7 +237,6 @@
const int size = numQuadPts * dim*dim;
scalar_array deform(size);
scalar_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
- scalar_array vertices(verticesVals, numBasis*dim);
const PylithScalar tolerance = 1.0e-06;
@@ -257,7 +244,7 @@
// ux(x,y,z) = 0.5
// uy(x,y,z) = 0.2
// uz(x,y,z) = 0.3
- const PylithScalar dispVals[] = {
+ const PylithScalar disp[numBasis*dim] = {
0.5, 0.2, 0.3,
0.5, 0.2, 0.3,
0.5, 0.2, 0.3,
@@ -267,13 +254,10 @@
0.5, 0.2, 0.3,
0.5, 0.2, 0.3,
};
- const PylithScalar deformE[] = {
+ const PylithScalar deformE[size] = {
1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
- 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
};
- scalar_array disp(dispVals, numBasis*dim);
-
IntegratorElasticityLgDeform::_calcDeformation(&deform, basisDeriv, vertices, disp, numBasis, numQuadPts, dim);
CPPUNIT_ASSERT_EQUAL(size, int(deform.size()));
@@ -285,7 +269,7 @@
// Let ux(x,y,z) = +0.4 + 0.3*x + 0.8*y - 0.2*z
// Ley uy(x,y,z) = -2.0 + 0.5*x - 0.2*y + 0.6*z
// Ley uz(x,y,z) = +0.7 + 0.8*x - 0.9*y - 0.1*z
- const PylithScalar dispVals[] = {
+ const PylithScalar disp[numBasis*dim] = {
+0.4+0.3*-1+0.8*-1-0.2*-1,
-2.0+0.5*-1-0.2*-1+0.6*-1,
+0.7+0.8*-1-0.9*-1-0.1*-1,
@@ -319,14 +303,12 @@
+0.7+0.8*-1-0.9*+1-0.1*+1,
};
- const PylithScalar deformE[] = {
+ const PylithScalar deformE[size] = {
1.3, 0.8, -0.2,
0.5, 0.8, 0.6,
0.8, -0.9, 0.9,
};
- scalar_array disp(dispVals, numBasis*dim);
-
IntegratorElasticityLgDeform::_calcDeformation(&deform, basisDeriv, vertices, disp, numBasis, numQuadPts, dim);
CPPUNIT_ASSERT_EQUAL(size, int(deform.size()));
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature0D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature0D.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature0D.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -81,7 +81,7 @@
Quadrature0D engine(refCell);
engine.initialize();
- engine.computeGeometry(vertCoords, 0);
+ engine.computeGeometry(&vertCoords[0], vertCoords.size(), 0);
const PylithScalar tolerance = 1.0e-06;
CPPUNIT_ASSERT_DOUBLES_EQUAL(quadPts[0], engine._quadPts[0], tolerance);
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadratureEngine.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadratureEngine.cc 2013-05-08 22:27:32 UTC (rev 22003)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadratureEngine.cc 2013-05-08 22:30:43 UTC (rev 22004)
@@ -198,7 +198,7 @@
// Check values from computeGeometry()
engine->initialize();
- engine->computeGeometry(vertCoords, 0);
+ engine->computeGeometry(&vertCoords[0], vertCoords.size(), 0);
const PylithScalar tolerance = 1.0e-06;
int size = numQuadPts * spaceDim;
More information about the CIG-COMMITS
mailing list