[cig-commits] r6727 - short/3D/PyLith/trunk/libsrc/feassemble
willic3 at geodynamics.org
willic3 at geodynamics.org
Mon Apr 30 13:27:57 PDT 2007
Author: willic3
Date: 2007-04-30 13:27:57 -0700 (Mon, 30 Apr 2007)
New Revision: 6727
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc
Log:
Possibly complete version of ImplicitElasticity.cc. Still need to test.
Definitely need to double-check that clearing and update of global stiffness
is done correctly.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc 2007-04-28 06:48:27 UTC (rev 6726)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc 2007-04-30 20:27:57 UTC (rev 6727)
@@ -226,10 +226,9 @@
// Compute B(transpose) * sigma, first computing strains
if (1 == cellDim) {
- // Compute total strains
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis)
- totalStrain[iQuad] += basisDeriv[iQ+iBasis] * dispTCell[iBasis];
+ // Compute total strains and then use these to compute stresses
+ IntegratorElasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
const std::vector<double_array>& stress =
_material->calcStress(totalStrain);
@@ -248,18 +247,9 @@
throw std::runtime_error("Logging PETSc flops failed.");
} else if (2 == cellDim) {
- // Compute total strains
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- totalStrain[iQuad][0] +=
- basisDeriv[iQ+iBasis ] * dispTCell[iBasis ];
- totalStrain[iQuad][1] +=
- basisDeriv[iQ+iBasis+1] * dispTCell[iBasis+1];
- totalStrain[iQuad][2] +=
- 0.5 * (basisDeriv[iQ+iBasis+1] * dispTCell[iBasis ] +
- basisDeriv[iQ+iBasis ] * dispTCell[iBasis+1]);
- } // for
- } // for
+ // Compute total strains and then use these to compute stresses
+ IntegratorElasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
const std::vector<double_array>& stress =
_material->calcStress(totalStrain);
@@ -282,26 +272,9 @@
throw std::runtime_error("Logging PETSc flops failed.");
} else if (3 == cellDim) {
- // Compute total strains
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- totalStrain[iQuad][0] +=
- basisDeriv[iQ+iBasis ] * dispTCell[iBasis ];
- totalStrain[iQuad][1] +=
- basisDeriv[iQ+iBasis+1] * dispTCell[iBasis+1];
- totalStrain[iQuad][2] +=
- basisDeriv[iQ+iBasis+2] * dispTCell[iBasis+2];
- totalStrain[iQuad][3] +=
- 0.5 * (basisDeriv[iQ+iBasis+1] * dispTCell[iBasis ] +
- basisDeriv[iQ+iBasis ] * dispTCell[iBasis+1]);
- totalStrain[iQuad][4] +=
- 0.5 * (basisDeriv[iQ+iBasis+2] * dispTCell[iBasis+1] +
- basisDeriv[iQ+iBasis+1] * dispTCell[iBasis+2]);
- totalStrain[iQuad][5] +=
- 0.5 * (basisDeriv[iQ+iBasis+2] * dispTCell[iBasis ] +
- basisDeriv[iQ+iBasis ] * dispTCell[iBasis+2]);
- } // for
- } // for
+ // Compute total strains and then use these to compute stresses
+ IntegratorElasticity::calcTotalStrain3D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
const std::vector<double_array>& stress =
_material->calcStress(totalStrain);
@@ -346,10 +319,13 @@
{ // integrateJacobian
assert(0 != _quadrature);
assert(0 != _material);
- assert(0 != _mat);
+ assert(0 != mat);
assert(!dispT.isNull());
assert(!mesh.isNull());
+ // Clear stiffness matrix. Not sure if this is the correct way.
+ PetscErrorCode err = MatZeroEntries(mat);
+
// Get information about section
const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
assert(!cells.isNull());
@@ -368,14 +344,30 @@
const double_array& quadWts = _quadrature->quadWts();
const int numBasis = _quadrature->numCorners();
const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
- // Compute Jacobian for cell, specific for each geometry type
if (cellDim != spaceDim)
throw std::logic_error("Not implemented yet.")
// Allocate vector for cell values (if necessary)
_initCellMatrix();
+ // 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;
@@ -394,37 +386,32 @@
const double_array& basisDeriv = _quadrature->basisDerivQuad();
const double_array& jacobianDet = _quadrature->jacobianDetQuad();
- // Need to finish fixing from here***************************
// 1D Case
if (1 == cellDim) {
assert(1 == numElasticConsts);
// Compute strains
- /* for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const int iBlock = iBasis * numBasis;
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- const int jBlock = jBasis; */
+ IntegratorElasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
// Get "elasticity" matrix at quadrature points for this cell
- _material->calcDerivElastic(*cellIter, patch, numQuadPts);
- const double* elasticConsts = _material->elasticConsts();
+ const double_array& elasticConsts = _material->calcDerivElastic(&totalStrain);
const int numElasticConsts = _material->numElasticConsts();
// Compute Jacobian for consistent tangent matrix
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double C1111 = elasticConsts[iQuad];
+ const double C1111 = elasticConsts[iQuad][0];
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const int iBlock = iBasis * spaceDim;
const double valI = wt*basisDeriv[iQ+iBasis]*C1111;
for (int jBasis=0; jBasis < numBasis; ++jBasis] {
const int jBlock = jBasis * spaceDim;
const double valIJ = valI * basisDeriv[iQ+jBasis];
- _cellMatrix[iBlock][jBlock] += valIJ;
+ _cellMatrix[iBlock+jBlock] += valIJ;
} // for
} // for
} // for
- PetscErrorCode err =
+ err =
PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*3)));
if (err)
throw std::runtime_error("Logging PETSc flops failed.");
@@ -433,34 +420,29 @@
} else if (2 == cellDim) {
assert(6 == numElasticConsts);
// Compute strains
- /* for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const int iBlock = iBasis * numBasis;
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- const int jBlock = jBasis; */
+ IntegratorElasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
// Get "elasticity" matrix at quadrature points for this cell
- _material->calcDerivElastic(*cellIter, patch, numQuadPts);
- const double* elasticConsts = _material->elasticConsts();
+ const double_array& elasticConsts = _material->calcDerivElastic(&totalStrain);
const int numElasticConsts = _material->numElasticConsts();
// Compute Jacobian for consistent tangent matrix
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const int iConst = iQuad*numElasticConsts;
- const double C1111 = elasticConsts[iConst+0];
- const double C1122 = elasticConsts[iConst+1];
- const double C1112 = elasticConsts[iConst+2];
- const double C2222 = elasticConsts[iConst+3];
- const double C2212 = elasticConsts[iConst+4];
- const double C1212 = elasticConsts[iConst+5];
+ const double C1111 = elasticConsts[iQuad][0];
+ const double C1122 = elasticConsts[iQuad][1];
+ const double C1112 = elasticConsts[iQuad][2];
+ const double C2222 = elasticConsts[iQuad][3];
+ const double C2212 = elasticConsts[iQuad][4];
+ const double C1212 = elasticConsts[iQuad][5];
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const int iBlock = iBasis * spaceDim;
- const double Nip = wt*basisDeriv[iQ+iBasis*cellDim+0];
+ const double Nip = wt*basisDeriv[iQ+iBasis*cellDim ];
const double Niq = wt*basisDeriv[iQ+iBasis*cellDim+1];
for (int jBasis=0; jBasis < numBasis; ++jBasis) {
const int jBlock = jBasis * spaceDim;
- const double Njp = basisDeriv[iQ+jBasis*cellDim+0];
+ const double Njp = basisDeriv[iQ+jBasis*cellDim ];
const double Njq = basisDeriv[iQ+jBasis*cellDim+1];
const double ki0j0 =
C1111 * Nip * Njp + C1112 * Niq * Njp +
@@ -471,14 +453,14 @@
const double ki1j1 =
C2222 * Niq * Njq + C2212 * Nip * Njq +
C2212 * Niq * Njp + C1212 * Nip * Njp;
- _cellMatrix[iBlock ][jBlock ] += ki0j0;
- _cellMatrix[iBlock ][jBlock+1] += ki0j1;
- _cellMatrix[iBlock+1][jBlock ] += ki0j1;
- _cellMatrix[iBlock+1][jBlock+1] += ki1j1;
+ _cellMatrix[iBlock +jBlock ] += ki0j0;
+ _cellMatrix[iBlock +jBlock+1] += ki0j1;
+ _cellMatrix[iBlock+1+jBlock ] += ki0j1;
+ _cellMatrix[iBlock+1+jBlock+1] += ki1j1;
} // for
} // for
} // for
- PetscErrorCode err =
+ err =
PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*(3*11+4))));
if (err)
throw std::runtime_error("Logging PETSc flops failed.");
@@ -487,42 +469,37 @@
} else if (3 == cellDim) {
assert(21 == numElasticConsts);
// Compute strains
- /* for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const int iBlock = iBasis * numBasis;
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- const int jBlock = jBasis; */
+ IntegratorElasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
+ dispTCell, numBasis);
// Get "elasticity" matrix at quadrature points for this cell
- _material->calcDerivElastic(*cellIter, patch, numQuadPts);
- const double* elasticConsts = _material->elasticConsts();
+ const double_array& elasticConsts = _material->calcDerivElastic(&totalStrain);
const int numElasticConsts = _material->numElasticConsts();
// Compute Jacobian for consistent tangent matrix
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const int iConst = iQuad*numElasticConsts;
- const double C1111 = elasticConsts[iConst+ 0];
- const double C1122 = elasticConsts[iConst+ 1];
- const double C1133 = elasticConsts[iConst+ 2];
- const double C1112 = elasticConsts[iConst+ 3];
- const double C1123 = elasticConsts[iConst+ 4];
- const double C1113 = elasticConsts[iConst+ 5];
- const double C2222 = elasticConsts[iConst+ 6];
- const double C2233 = elasticConsts[iConst+ 7];
- const double C2212 = elasticConsts[iConst+ 8];
- const double C2223 = elasticConsts[iConst+ 9];
- const double C2213 = elasticConsts[iConst+10];
- const double C3333 = elasticConsts[iConst+11];
- const double C3312 = elasticConsts[iConst+12];
- const double C3323 = elasticConsts[iConst+13];
- const double C3313 = elasticConsts[iConst+14];
- const double C1212 = elasticConsts[iConst+15];
- const double C1223 = elasticConsts[iConst+16];
- const double C1213 = elasticConsts[iConst+17];
- const double C2323 = elasticConsts[iConst+18];
- const double C2313 = elasticConsts[iConst+19];
- const double C1313 = elasticConsts[iConst+20];
+ const double C1111 = elasticConsts[iQuad][ 0];
+ const double C1122 = elasticConsts[iQuad][ 1];
+ const double C1133 = elasticConsts[iQuad][ 2];
+ const double C1112 = elasticConsts[iQuad][ 3];
+ const double C1123 = elasticConsts[iQuad][ 4];
+ const double C1113 = elasticConsts[iQuad][ 5];
+ const double C2222 = elasticConsts[iQuad][ 6];
+ const double C2233 = elasticConsts[iQuad][ 7];
+ const double C2212 = elasticConsts[iQuad][ 8];
+ const double C2223 = elasticConsts[iQuad][ 9];
+ const double C2213 = elasticConsts[iQuad][10];
+ const double C3333 = elasticConsts[iQuad][11];
+ const double C3312 = elasticConsts[iQuad][12];
+ const double C3323 = elasticConsts[iQuad][13];
+ const double C3313 = elasticConsts[iQuad][14];
+ const double C1212 = elasticConsts[iQuad][15];
+ const double C1223 = elasticConsts[iQuad][16];
+ const double C1213 = elasticConsts[iQuad][17];
+ const double C2323 = elasticConsts[iQuad][18];
+ const double C2313 = elasticConsts[iQuad][19];
+ const double C1313 = elasticConsts[iQuad][20];
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const int iBlock = iBasis * spaceDim;
const double Nip = wt*basisDeriv[iQ+iBasis*cellDim+0];
@@ -558,26 +535,30 @@
C3323 * Nir * Njq + C2323 * Niq * Njq + C2313 * Nip * Njq +
C3313 * Nir * Njp + C2313 * Niq * Njp + C1313 * Nip * Njp;
- _cellMatrix[iblock ][jBlock ] += ki0j0;
- _cellMatrix[iblock+1][jBlock ] += ki0j1;
- _cellMatrix[iblock+2][jBlock ] += ki0j2;
- _cellMatrix[iblock ][jBlock+1] += ki0j1;
- _cellMatrix[iblock+1][jBlock+1] += ki1j1;
- _cellMatrix[iblock+2][jBlock+1] += ki1j2;
- _cellMatrix[iblock ][jBlock+2] += ki0j2;
- _cellMatrix[iblock+1][jBlock+2] += ki1j2;
- _cellMatrix[iblock+2][jBlock+2] += ki2j2;
+ _cellMatrix[iblock +jBlock ] += ki0j0;
+ _cellMatrix[iblock+1+jBlock ] += ki0j1;
+ _cellMatrix[iblock+2+jBlock ] += ki0j2;
+ _cellMatrix[iblock +jBlock+1] += ki0j1;
+ _cellMatrix[iblock+1+jBlock+1] += ki1j1;
+ _cellMatrix[iblock+2+jBlock+1] += ki1j2;
+ _cellMatrix[iblock +jBlock+2] += ki0j2;
+ _cellMatrix[iblock+1+jBlock+2] += ki1j2;
+ _cellMatrix[iblock+2+jBlock+2] += ki2j2;
} // for
} // for
} // for
- PetscErrorCode err =
+ err =
PetscLogFlops(numQuadPts*(1+numBasis*(3+numBasis*(6*26+9))));
if (err)
throw std::runtime_error("Logging PETSc flops failed.");
} // if/else
- // Assemble cell contribution into field
- err = assembleMatrix(*mat, *cellIter, _cellMatrix, ADD_VALUES);
+ // Assemble cell contribution into field. Not sure if this is correct for
+ // global stiffness matrix.
+ const ALE::Obj<Mesh::order_type>& globalOrder =
+ mesh->getFactory()->getGlobalOrder(mesh, "default", dispT);
+ err = updateOperator(*mat, mesh, dispT, globalOrder,
+ *cellIter, _cellMatrix, ADD_VALUES);
if (err)
throw std::runtime_error("Update to PETSc Mat failed.");
} // for
More information about the cig-commits
mailing list