[cig-commits] r6502 - short/3D/PyLith/trunk/libsrc/feassemble

[willic3 at geodynamics.org]

Author: willic3
Date: 2007-04-02 11:46:55 -0700 (Mon, 02 Apr 2007)
New Revision: 6502

Modified:
   short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc
Log:
Put in missing parts of integrateResidual taken from ExplicitElasticity.cc.
This includes portions that compute total strain and then use material
routines to compute stress, then compute the (negative) element internal
force, B(transpose) * stress.


Modified: short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc	2007-04-02 16:30:39 UTC (rev 6501)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc	2007-04-02 18:46:55 UTC (rev 6502)
@@ -182,9 +182,136 @@
     const double* basis = _quadrature->basisDeriv();
     const double* jacobianDet = _quadrature->jacobianDet();
 
+    // 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 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);
 
+      // 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  ];
+	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  ] += 
+	    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
+      const double* stress = 
+	_material->calcStress(*cellIter, totalStrain, numQuadPts, 
+			      spaceDim);
+      // 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];
+	for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+	  const int iBlock = iBasis * spaceDim;
+	  const double N1 = wt*basisDeriv[iQ+iBasis*cellDim  ];
+	  const double N2 = wt*basisDeriv[iQ+iBasis*cellDim+1];
+	  _cellVector[iBlock  ] -= N1*s11 + N2*s12;
+	  _cellVector[iBlock+1] -= N1*s12 + N2*s22;
+	} // for
+      } // for
+      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  ] += 
+	    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
+      const double* stress = 
+	_material->calcStress(*cellIter, totalStrain, numQuadPts, spaceDim);
+      // 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];
 
+	for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+	  const int iBlock = iBasis * spaceDim;
+	  const double N1 = wt*basisDeriv[iQ+iBasis*cellDim+0];
+	  const double N2 = wt*basisDeriv[iQ+iBasis*cellDim+1];
+	  const double N3 = wt*basisDeriv[iQ+iBasis*cellDim+2];
+	  _cellVector[iBlock  ] -= N1*s11 + N2*s12 + N3*s13;
+	  _cellVector[iBlock+1] -= N1*s12 + N2*s22 + N3*s23;
+	  _cellVector[iBlock+2] -= N1*s13 + N2*s23 + N3*s33;
+	} // for
+      } // for
+      err = PetscLogFlops(numQuadPts*(1+numBasis*(3+12)));
+      if (err)
+	throw std::runtime_error("Logging PETSc flops failed.");
+    } // if/else
 
+   // Assemble cell contribution into field
+    mesh->updateAdd(b, *cellIter, _cellVector);
+  } // for
+} // integrateResidual
+
+
 // ----------------------------------------------------------------------
 // Compute stiffness matrix.
 void