[cig-commits] r6479 - short/3D/PyLith/trunk/libsrc/feassemble
willic3 at geodynamics.org
willic3 at geodynamics.org
Fri Mar 30 18:24:06 PDT 2007
Author: willic3
Date: 2007-03-30 18:24:06 -0700 (Fri, 30 Mar 2007)
New Revision: 6479
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc
Log:
Incomplete version of ImplicitElasticity.cc.
Contains:
integrateConstant (only contribution at present is body forces)
integrateResidual (this is incomplete and needs strain computation and stress integration
that I will get from Brad's code)
integrateJacobian (computation of stiffness matrix)
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc 2007-03-31 00:39:07 UTC (rev 6478)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc 2007-03-31 01:24:06 UTC (rev 6479)
@@ -57,15 +57,19 @@
} // material
// ----------------------------------------------------------------------
-// Integrate residual term (b) for quasi-static elasticity term for 3-D
-// finite elements. This is the contribution from element internal forces.
+// Integrate portion of RHS that is assumed to remain constant over the
+// entire time step. At present, this only includes body forces.
+// Note: The vector produced by this computation remains unchanged for
+// small strain computations, so it might be reasonable to keep the vector
+// separate from the other contributions to the RHS.
+
void
-pylith::feassemble::ImplicitElasticity::integrateResidual(
+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)
-{ // integrateResidual
+{ // integrateConstant
assert(0 != _quadrature);
// Get information about section
@@ -100,14 +104,10 @@
// Get cell geometry information that depends on cell
const double* basis = _quadrature->basis();
- const double* basisDeriv = _quadrature->basisDeriv();
const double* jacobianDet = _quadrature->jacobianDet();
- // Get material physical properties at quadrature points for this cell
- _material->calcProperties(*cellIter, patch, numQuadPts);
- const double* density = _material->density();
- const double* elasticConsts = _material->elasticConsts();
- const int numElasticConsts = _material->numElasticConsts();
+ // Get density at quadrature points for this cell
+ const double* density = _material->calcDensity(*cellIter, patch, numQuadPts);
// Compute action for cell
@@ -131,170 +131,64 @@
throw std::runtime_error("Logging PETSc flops failed.");
} // if
- // Compute action for element internal forces
-
-
- /** :TODO:
- *
- * If cellDim and spaceDim are different, we need to transform
- * displacements into cellDim, compute action, and transform
- * result back into spaceDim. Can we get this from the inverse of
- * the Jacobian?
- */
- if (cellDim != spaceDim)
- throw std::logic_error("Not implemented yet.");
-
- // Compute action for elastic terms
- if (1 == cellDim) {
- assert(1 == numElasticConsts);
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double C1111 = elasticConsts[iQuad];
- 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];
- _cellVector[iBlock] -= valIJ * dispTCell[jBlock];
- } // for
- } // for
- } // for
- PetscErrorCode err =
- PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*3)));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
- } else if (2 == cellDim) {
- assert(6 == numElasticConsts);
- 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];
- 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 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 Njq = basisDeriv[iQ+jBasis*cellDim+1];
- const double ki0j0 =
- C1111 * Nip * Njp + C1112 * Niq * Njp +
- C1112 * Nip * Njq + C1212 * Niq * Njq;
- const double ki0j1 =
- C1122 * Nip * Njq + C2212 * Niq * Njq +
- C1112 * Nip * Njp + C1212 * Niq * Njp;
- const double ki1j1 =
- C2222 * Niq * Njq + C2212 * Nip * Njq +
- C2212 * Niq * Njp + C1212 * Nip * Njp;
- _cellVector[iBlock ] -=
- ki0j0 * dispTCell[jBlock ] + ki0j1 * dispTCell[jBlock+1];
- _cellVector[iBlock+1] -=
- ki0j1 * dispTCell[jBlock ] + ki1j1 * dispTCell[jBlock+1];
- } // for
- } // for
- } // for
- PetscErrorCode err =
- PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*(2+3*11+2*4))));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
- } else if (3 == cellDim) {
- assert(21 == numElasticConsts);
- 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];
- 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 Niq = wt*basisDeriv[iQ+iBasis*cellDim+1];
- const double Nir = wt*basisDeriv[iQ+iBasis*cellDim+2];
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- const int jBlock = jBasis * spaceDim;
- const double Njp = basisDeriv[iQ+jBasis*cellDim+0];
- const double Njq = basisDeriv[iQ+jBasis*cellDim+1];
- const double Njr = basisDeriv[iQ+jBasis*cellDim+2];
- const double ki0j0 =
- C1111 * Nip * Njp + C1112 * Niq * Njp + C1113 * Nir * Njp +
- C1112 * Nip * Njq + C1212 * Niq * Njq + C1213 * Nir * Njq +
- C1113 * Nip * Njr + C1213 * Niq * Njr + C1313 * Nir * Njr;
- const double ki0j1 =
- C1122 * Nip * Njq + C2212 * Niq * Njq + C2213 * Nir * Njq +
- C1112 * Nip * Njp + C1212 * Niq * Njp + C1213 * Nir * Njp +
- C1123 * Nip * Njr + C1223 * Niq * Njr + C2313 * Nir * Njr;
- const double ki0j2 =
- C1133 * Nip * Njr + C3312 * Niq * Njr + C3313 * Nir * Njr +
- C1123 * Nip * Njq + C1223 * Niq * Njq + C2313 * Nir * Njq +
- C1113 * Nip * Njp + C1213 * Niq * Njp + C1313 * Nir * Njp;
- const double ki1j1 =
- C2222 * Niq * Njq + C2212 * Nip * Njq + C2223 * Nir * Njq +
- C2212 * Niq * Njp + C1212 * Nip * Njp + C1223 * Nir * Njp +
- C2223 * Niq * Njr + C1223 * Nip * Njr + C2323 * Nir * Njr;
- const double ki1j2 =
- C2233 * Niq * Njr + C3312 * Nip * Njr + C3323 * Nir * Njr +
- C2223 * Niq * Njq + C1223 * Nip * Njq + C2323 * Nir * Njq +
- C2213 * Niq * Njp + C1213 * Nip * Njp + C2313 * Nir * Njp;
- const double ki2j2 =
- C3333 * Nir * Njr + C3323 * Niq * Njr + C3313 * Nip * Njr +
- C3323 * Nir * Njq + C2323 * Niq * Njq + C2313 * Nip * Njq +
- C3313 * Nir * Njp + C2313 * Niq * Njp + C1313 * Nip * Njp;
-
- _cellVector[iBlock ] -=
- ki0j0 * dispTCell[jBlock ] +
- ki0j1 * dispTCell[jBlock+1] +
- ki0j2 * dispTCell[jBlock+2];
- _cellVector[iBlock+1] -=
- ki0j1 * dispTCell[jBlock ] +
- ki1j1 * dispTCell[jBlock+1] +
- ki1j2 * dispTCell[jBlock+2];
- _cellVector[iBlock+2] -=
- ki0j2 * dispTCell[jBlock ] +
- ki1j2 * dispTCell[jBlock+1] +
- ki2j2 * dispTCell[jBlock+2];
- } // for
- } // for
- } // for
- PetscErrorCode err =
- PetscLogFlops(numQuadPts*(1+numBasis*(3+numBasis*(3+6*26+3*6))));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
- } // if/else
-
// Assemble cell contribution into field
b->updateAdd(patch, *cellIter, _cellVector);
} // for
} // integrateConstant
// ----------------------------------------------------------------------
-// Compute matrix associated with operator.
+// Compute residual for quasi-static finite elements.
void
-pylith::feassemble::ExplicitElasticity::integrateJacobian(
+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)
+{ // integrateResidual
+ assert(0 != _quadrature);
+
+ // 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();
+
+ // Allocate vector for cell values (if necessary)
+ _initCellVector();
+
+ // 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 int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+
+ for (topology_type::label_sequence::iterator cellIter=cells->begin();
+ cellIter != cellsEnd;
+ ++cellIter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(coordinates, *cellIter);
+
+ // 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* basis = _quadrature->basisDeriv();
+ const double* jacobianDet = _quadrature->jacobianDet();
+
+
+
+
+// ----------------------------------------------------------------------
+// Compute stiffness matrix.
+void
+pylith::feassemble::ImplicitElasticity::integrateJacobian(
PetscMat* mat,
const ALE::Obj<real_section_type>& dispT,
const ALE::Obj<real_section_type>& coordinates)
@@ -310,7 +204,6 @@
// Get parameters used in integration.
const double dt = _dt;
- const double dt2 = dt*dt;
// Allocate vector for cell values (if necessary)
_initCellMatrix();
@@ -331,34 +224,194 @@
// Get cell geometry information that depends on cell
const double* basis = _quadrature->basis();
+ const double* basisDeriv = _quadrature->basisDeriv();
const double* jacobianDet = _quadrature->jacobianDet();
- // Get material physical properties at quadrature points for this cell
- _material->calcProperties(*cellIter, patch, numQuadPts);
- const double* density = _material->density();
+ // Compute Jacobian for cell, specific for each geometry type
+ if (cellDim != spaceDim)
+ throw std::logic_error("Not implemented yet.")
- // Compute Jacobian for cell
+ // 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; */
- // Compute Jacobian for inertial terms
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt =
- quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad] / dt2;
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const int iBlock = (iBasis * spaceDim) * (spaceDim * numBasis);
- const double valI = wt*basis[iQ+iBasis];
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- const int jBlock = (jBasis * spaceDim);
- const double valIJ = valI * basis[iQ+jBasis];
- for (int iDim=0; iDim < spaceDim; ++iDim)
- _cellMatrix[iBlock+jBlock+iDim*spaceDim+iDim] += valIJ;
+ // Get "elasticity" matrix at quadrature points for this cell
+ _material->calcDerivElastic(*cellIter, patch, numQuadPts);
+ const double* elasticConsts = _material->elasticConsts();
+ 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];
+ 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;
+ } // for
+ } // for
+ } // for
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*3)));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ // 2D Case
+ } 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; */
+
+ // Get "elasticity" matrix at quadrature points for this cell
+ _material->calcDerivElastic(*cellIter, patch, numQuadPts);
+ const double* elasticConsts = _material->elasticConsts();
+ 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];
+ 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 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 Njq = basisDeriv[iQ+jBasis*cellDim+1];
+ const double ki0j0 =
+ C1111 * Nip * Njp + C1112 * Niq * Njp +
+ C1112 * Nip * Njq + C1212 * Niq * Njq;
+ const double ki0j1 =
+ C1122 * Nip * Njq + C2212 * Niq * Njq +
+ C1112 * Nip * Njp + C1212 * Niq * Njp;
+ 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;
+ } // for
} // for
} // for
- } // for
- PetscErrorCode err =
- PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(1+spaceDim))));
- if (err)
- throw std::runtime_error("Logging PETSc flops failed.");
-
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*(3*11+4))));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ // 3D Case
+ } 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; */
+
+ // Get "elasticity" matrix at quadrature points for this cell
+ _material->calcDerivElastic(*cellIter, patch, numQuadPts);
+ const double* elasticConsts = _material->elasticConsts();
+ 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];
+ 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 Niq = wt*basisDeriv[iQ+iBasis*cellDim+1];
+ const double Nir = wt*basisDeriv[iQ+iBasis*cellDim+2];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const int jBlock = jBasis * spaceDim;
+ const double Njp = basisDeriv[iQ+jBasis*cellDim+0];
+ const double Njq = basisDeriv[iQ+jBasis*cellDim+1];
+ const double Njr = basisDeriv[iQ+jBasis*cellDim+2];
+ const double ki0j0 =
+ C1111 * Nip * Njp + C1112 * Niq * Njp + C1113 * Nir * Njp +
+ C1112 * Nip * Njq + C1212 * Niq * Njq + C1213 * Nir * Njq +
+ C1113 * Nip * Njr + C1213 * Niq * Njr + C1313 * Nir * Njr;
+ const double ki0j1 =
+ C1122 * Nip * Njq + C2212 * Niq * Njq + C2213 * Nir * Njq +
+ C1112 * Nip * Njp + C1212 * Niq * Njp + C1213 * Nir * Njp +
+ C1123 * Nip * Njr + C1223 * Niq * Njr + C2313 * Nir * Njr;
+ const double ki0j2 =
+ C1133 * Nip * Njr + C3312 * Niq * Njr + C3313 * Nir * Njr +
+ C1123 * Nip * Njq + C1223 * Niq * Njq + C2313 * Nir * Njq +
+ C1113 * Nip * Njp + C1213 * Niq * Njp + C1313 * Nir * Njp;
+ const double ki1j1 =
+ C2222 * Niq * Njq + C2212 * Nip * Njq + C2223 * Nir * Njq +
+ C2212 * Niq * Njp + C1212 * Nip * Njp + C1223 * Nir * Njp +
+ C2223 * Niq * Njr + C1223 * Nip * Njr + C2323 * Nir * Njr;
+ const double ki1j2 =
+ C2233 * Niq * Njr + C3312 * Nip * Njr + C3323 * Nir * Njr +
+ C2223 * Niq * Njq + C1223 * Nip * Njq + C2323 * Nir * Njq +
+ C2213 * Niq * Njp + C1213 * Nip * Njp + C2313 * Nir * Njp;
+ const double ki2j2 =
+ C3333 * Nir * Njr + C3323 * Niq * Njr + C3313 * Nip * Njr +
+ 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;
+ } // for
+ } // for
+ } // for
+ PetscErrorCode 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);
if (err)
More information about the cig-commits
mailing list