[cig-commits] r6625 - short/3D/PyLith/trunk/libsrc/feassemble
willic3 at geodynamics.org
willic3 at geodynamics.org
Fri Apr 20 13:49:20 PDT 2007
Author: willic3
Date: 2007-04-20 13:49:20 -0700 (Fri, 20 Apr 2007)
New Revision: 6625
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc
Log:
Started fixing this routine for new PETSc usage and to use arrays/vectors.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc 2007-04-20 20:33:46 UTC (rev 6624)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc 2007-04-20 20:49:20 UTC (rev 6625)
@@ -16,6 +16,7 @@
#include "Quadrature.hh" // USES Quadrature
#include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
+#include "pylith/utils/array.hh" // USES double_array
#include "petscmat.h" // USES PetscMat
#include "spatialdata/spatialdb/SpatialDB.hh"
@@ -67,54 +68,64 @@
pylith::feassemble::ImplicitElasticity::integrateConstant(
const ALE::Obj<real_section_type>& b,
const ALE::Obj<real_section_type>& dispT,
- const ALE::Obj<real_section_type>& grav,
- const ALE::Obj<real_section_type>& coordinates)
+ const ALE::Obj<Mesh>& mesh)
{ // integrateConstant
assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(!b.isNull());
+ assert(!dispT.isNull());
+ assert(!mesh.isNull());
+ /* Comment out the guts of this code until we figure out where grav
+ section is coming from. Code should be essentially correct, though.
+ It may also be necessary to put in a switch to turn gravity off.
+
// Get information about section
- const topology_type::patch_type patch = 0;
- const ALE::Obj<topology_type>& topology = dispT->getTopology();
- const ALE::Obj<topology_type::label_sequence>& cells =
- topology->heightStratum(patch, 0);
- const topology_type::label_sequence::iterator cellsEnd = cells->end();
+ const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
- // Allocate vector for cell values (if necessary)
- _initCellVector();
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
// Get cell geometry information that doesn't depend on cell
const int numQuadPts = _quadrature->numQuadPts();
- const double* quadWts = _quadrature->quadWts();
- const int numBasis = _quadrature->numCorners();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
const int cellDim = _quadrature->cellDim();
- for (topology_type::label_sequence::iterator cellIter=cells->begin();
+ // Allocate vector for cell values (if necessary)
+ _initCellVector();
+
+ // Loop over cells.
+ for (Mesh::label_sequence::iterator cellIter=cells->begin();
cellIter != cellsEnd;
++cellIter) {
// Compute geometry information for current cell
- _quadrature->computeGeometry(coordinates, *cellIter);
+ _quadrature->computeGeometry(mesh, coordinates, *cellIter);
+ // Set cell data in material
+ _material->initCellData(*cellIter, numQuadPts);
+
// Reset element vector to zero
_resetCellVector();
- // Restrict input fields to cell
- const real_section_type::value_type* dispTCell =
- dispT->restrict(patch, *cellIter);
-
// Get cell geometry information that depends on cell
- const double* basis = _quadrature->basis();
- const double* jacobianDet = _quadrature->jacobianDet();
+ const double_array& basis = _quadrature->basis();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
// Get density at quadrature points for this cell
- const double* density = _material->calcDensity(*cellIter, patch, numQuadPts);
+ const std::vector<double_array>& density = _material->calcDensity();
// Compute action for cell
// Compute action for element body forces
if (!grav.isNull()) {
- const read_section_type::value_type* gravCell =
- grav->restrict(patch, cell);
+ const real_section_type::value_type* gravCell =
+ mesh->restrict(grav, cell);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad];
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
@@ -133,8 +144,9 @@
// Assemble cell contribution into field
- b->updateAdd(patch, *cellIter, _cellVector);
+ mesh->updateAdd(b, *cellIter, _cellVector);
} // for
+ */
} // integrateConstant
// ----------------------------------------------------------------------
@@ -143,104 +155,120 @@
pylith::feassemble::ImplicitElasticity::integrateResidual(
const ALE::Obj<real_section_type>& b,
const ALE::Obj<real_section_type>& dispT,
- const ALE::Obj<real_section_type>& coordinates)
+ const ALE::Obj<Mesh>& mesh)
{ // integrateResidual
assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(!b.isNull());
+ assert(!dispT.isNull());
+ assert(!mesh.isNull());
// Get information about section
- const topology_type::patch_type patch = 0;
- const ALE::Obj<topology_type>& topology = dispT->getTopology();
- const ALE::Obj<topology_type::label_sequence>& cells =
- topology->heightStratum(patch, 0);
- const topology_type::label_sequence::iterator cellsEnd = cells->end();
+ const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
- // Allocate vector for cell values (if necessary)
- _initCellVector();
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
// Get cell geometry information that doesn't depend on cell
const int numQuadPts = _quadrature->numQuadPts();
- const double* quadWts = _quadrature->quadWts();
- const int numBasis = _quadrature->numCorners();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
const int cellDim = _quadrature->cellDim();
- for (topology_type::label_sequence::iterator cellIter=cells->begin();
+ // Allocate vector for cell values (if necessary)
+ _initCellVector();
+
+ // Allocate vector for total strain
+ int tensorSize = 0;
+ if (1 == cellDim)
+ tensorSize = 1;
+ else if (2 == cellDim)
+ tensorSize = 3;
+ else if (3 == cellDim)
+ tensorSize = 6;
+ else
+ throw std::logic_error("Tensor size not implemented for given cellDim.");
+ std::vector<double_array> totalStrain(numQuadPts);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ totalStrain[iQuad].resize(tensorSize);
+ totalStrain[iQuad] = 0.0;
+ } // for
+
+ // Loop over cells
+ for (Mesh::label_sequence::iterator cellIter=cells->begin();
cellIter != cellsEnd;
++cellIter) {
// Compute geometry information for current cell
- _quadrature->computeGeometry(coordinates, *cellIter);
+ _quadrature->computeGeometry(mesh, coordinates, *cellIter);
+ // Set cell data in material
+ _material->initCellData(*cellIter, numQuadPts);
+
// Reset element vector to zero
_resetCellVector();
// Restrict input fields to cell
const real_section_type::value_type* dispTCell =
- dispT->restrict(patch, *cellIter);
+ mesh->restrict(dispT, *cellIter);
// Get cell geometry information that depends on cell
- const double* basis = _quadrature->basis();
- const double* basis = _quadrature->basisDeriv();
- const double* jacobianDet = _quadrature->jacobianDet();
+ const double_array& basis = _quadrature->basis();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+ if (cellDim != spaceDim)
+ throw std::logic_error("Not implemented yet.");
+
// Compute B(transpose) * sigma, first computing strains
if (1 == cellDim) {
// Compute total strains
- const int stressSize = _material->stressSize();
- assert(1 == stressSize);
- const int strainSize = stressSize * numQuadPts;
- double* totalStrain = (strainSize > 0) ? new double[strainSize] : 0;
- memset(totalStrain, 0, strainSize*sizeof(double));
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis)
totalStrain[iQuad] += basisDeriv[iQ+iBasis] * dispTCell[iBasis];
- } // for
- const double* stress =
- _material->calcStress(*cellIter, totalStrain, numQuadPts,
- spaceDim);
+ const std::vector<double_array>& stress =
+ _material->calcStress(totalStrain);
// Compute elastic action
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const int iStress = iQuad*stressSize;
- const double s11 = stress[iStress ];
+ const double s11 = stress[iQuad][0];
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const int iBlock = iBasis * spaceDim;
const double N1 = wt*basisDeriv[iQ+iBasis*cellDim ];
_cellVector[iBlock ] -= N1*s11;
} // for
} // for
- delete[] totalStrain; totalStrain = 0;
PetscErrorCode err = PetscLogFlops(numQuadPts*(1+numBasis*5));
if (err)
throw std::runtime_error("Logging PETSc flops failed.");
+
} else if (2 == cellDim) {
// Compute total strains
- const int stressSize = _material->stressSize();
- assert(3 == stressSize);
- const int strainSize = stressSize * numQuadPts;
- double* totalStrain = (strainSize > 0) ? new double[strainSize] : 0;
- memset(totalStrain, 0, strainSize*sizeof(double));
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- totalStrain[iQuad ] +=
+ totalStrain[iQuad][0] +=
basisDeriv[iQ+iBasis ] * dispTCell[iBasis ];
- totalStrain[iQuad+1] +=
+ totalStrain[iQuad][1] +=
basisDeriv[iQ+iBasis+1] * dispTCell[iBasis+1];
- totalStrain[iQuad+2] +=
+ totalStrain[iQuad][2] +=
0.5 * (basisDeriv[iQ+iBasis+1] * dispTCell[iBasis ] +
basisDeriv[iQ+iBasis ] * dispTCell[iBasis+1]);
} // for
} // for
- const double* stress =
- _material->calcStress(*cellIter, totalStrain, numQuadPts,
- spaceDim);
+ const std::vector<double_array>& stress =
+ _material->calcStress(totalStrain);
+
// Compute elastic action
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const int iStress = iQuad*stressSize;
- const double s11 = stress[iStress ];
- const double s22 = stress[iStress+1];
- const double s12 = stress[iStress+2];
+ const double s11 = stress[iQuad][0];
+ const double s22 = stress[iQuad][1];
+ const double s12 = stress[iQuad][2];
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const int iBlock = iBasis * spaceDim;
const double N1 = wt*basisDeriv[iQ+iBasis*cellDim ];
@@ -252,44 +280,40 @@
err = PetscLogFlops(numQuadPts*(1+numBasis*(8+2+9)));
if (err)
throw std::runtime_error("Logging PETSc flops failed.");
+
} else if (3 == cellDim) {
// Compute total strains
- const int stressSize = _material->stressSize();
- assert(6 == stressSize);
- const int strainSize = stressSize*numQuadPts;
- double* totalStrain = (strainSize > 0) ? new double[strainSize] : 0;
- memset(totalStrain, 0, strainSize*sizeof(double));
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- totalStrain[iQuad ] +=
+ totalStrain[iQuad][0] +=
basisDeriv[iQ+iBasis ] * dispTCell[iBasis ];
- totalStrain[iQuad+1] +=
+ totalStrain[iQuad][1] +=
basisDeriv[iQ+iBasis+1] * dispTCell[iBasis+1];
- totalStrain[iQuad+2] +=
+ totalStrain[iQuad][2] +=
basisDeriv[iQ+iBasis+2] * dispTCell[iBasis+2];
- totalStrain[iQuad+3] +=
+ totalStrain[iQuad][3] +=
0.5 * (basisDeriv[iQ+iBasis+1] * dispTCell[iBasis ] +
basisDeriv[iQ+iBasis ] * dispTCell[iBasis+1]);
- totalStrain[iQuad+4] +=
+ totalStrain[iQuad][4] +=
0.5 * (basisDeriv[iQ+iBasis+2] * dispTCell[iBasis+1] +
basisDeriv[iQ+iBasis+1] * dispTCell[iBasis+2]);
- totalStrain[iQuad+5] +=
+ totalStrain[iQuad][5] +=
0.5 * (basisDeriv[iQ+iBasis+2] * dispTCell[iBasis ] +
basisDeriv[iQ+iBasis ] * dispTCell[iBasis+2]);
} // for
} // for
- const double* stress =
- _material->calcStress(*cellIter, totalStrain, numQuadPts, spaceDim);
+ const std::vector<double_array>& stress =
+ _material->calcStress(totalStrain);
+
// Compute elastic action
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const int iStress = iQuad*stressSize;
- const double s11 = stress[iStress ];
- const double s22 = stress[iStress+1];
- const double s33 = stress[iStress+2];
- const double s12 = stress[iStress+3];
- const double s23 = stress[iStress+4];
- const double s13 = stress[iStress+5];
+ const double s11 = stress[iQuad][0];
+ const double s22 = stress[iQuad][1];
+ const double s33 = stress[iQuad][2];
+ const double s12 = stress[iQuad][3];
+ const double s23 = stress[iQuad][4];
+ const double s13 = stress[iQuad][5];
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const int iBlock = iBasis * spaceDim;
@@ -318,46 +342,59 @@
pylith::feassemble::ImplicitElasticity::integrateJacobian(
PetscMat* mat,
const ALE::Obj<real_section_type>& dispT,
- const ALE::Obj<real_section_type>& coordinates)
+ const ALE::Obj<Mesh>& mesh)
{ // integrateJacobian
assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(0 != _mat);
+ assert(!dispT.isNull());
+ assert(!mesh.isNull());
// Get information about section
- const topology_type::patch_type patch = 0;
- const ALE::Obj<topology_type>& topology = dispT->getTopology();
- const ALE::Obj<topology_type::label_sequence>& cells =
- topology->heightStratum(patch, 0);
- const topology_type::label_sequence::iterator cellsEnd = cells->end();
+ const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+
// Get parameters used in integration.
const double dt = _dt;
+ assert(dt > 0);
- // Allocate vector for cell values (if necessary)
- _initCellMatrix();
-
// Get cell geometry information that doesn't depend on cell
const int numQuadPts = _quadrature->numQuadPts();
- const double* quadWts = _quadrature->quadWts();
+ const double_array& quadWts = _quadrature->quadWts();
const int numBasis = _quadrature->numCorners();
const int spaceDim = _quadrature->spaceDim();
- for (topology_type::label_sequence::iterator cellIter=cells->begin();
+
+ // Allocate vector for cell values (if necessary)
+ _initCellMatrix();
+
+ // Loop over cells
+ for (Mesh::label_sequence::iterator cellIter=cells->begin();
cellIter != cellsEnd;
++cellIter) {
// Compute geometry information for current cell
- _quadrature->computeGeometry(coordinates, *cellIter);
+ _quadrature->computeGeometry(mesh, coordinates, *cellIter);
+ // Set cell data in material
+ _material->initCellData(*cellIter, numQuadPts);
+
// Reset element vector to zero
_resetCellMatrix();
// Get cell geometry information that depends on cell
- const double* basis = _quadrature->basis();
- const double* basisDeriv = _quadrature->basisDeriv();
- const double* jacobianDet = _quadrature->jacobianDet();
+ const double_array& basis = _quadrature->basis();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
// Compute Jacobian for cell, specific for each geometry type
if (cellDim != spaceDim)
throw std::logic_error("Not implemented yet.")
+ // Need to finish fixing from here***************************
// 1D Case
if (1 == cellDim) {
assert(1 == numElasticConsts);
More information about the cig-commits
mailing list