added: multi-threaded support for finite deformation applications

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1420 e10b68d5-8a6e-419e-a041-bce267b0401d

Apps/FiniteDefElasticity/CMakeLists.txt

@@ -78,6 +78,14 @@ IF(NOT "${DISABLE_HDF5}")
   FIND_PACKAGE(HDF5)
 ENDIF(NOT "${DISABLE_HDF5}")
 
+IF(NOT "${DISABLE_OPENMP}")
+  FIND_PACKAGE(OpenMP)
+ENDIF(NOT "${DISABLE_OPENMP}")
+
+IF(OPENMP_FOUND)
+  SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS} -DUSE_OPENMP")
+ENDIF(OPENMP_FOUND)
+
 # Required libraries
 SET(DEPLIBS ${IFEM_LIBRARIES}
   ${GoTrivariate_LIBRARIES} ${GoTools_LIBRARIES}

Apps/FiniteDefElasticity/NonlinearElasticity.C

@@ -14,6 +14,7 @@
 #include "NonlinearElasticity.h"
 #include "MaterialBase.h"
 #include "FiniteElement.h"
+#include "ElmMats.h"
 #include "Utilities.h"
 #include "Profiler.h"
 
@@ -54,7 +55,7 @@ extern "C" {
 
 
 NonlinearElasticity::NonlinearElasticity (unsigned short int n)
-  : NonlinearElasticityTL(n), E(n)
+  : NonlinearElasticityTL(n)
 {
   fullCmat = false;
 }
@@ -76,22 +77,27 @@ void NonlinearElasticity::setMode (SIM::SolutionMode mode)
 }
 
 
-bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt,
+bool NonlinearElasticity::evalInt (LocalIntegral& elmInt,
 				   const FiniteElement& fe,
 				   const Vec3& X) const
 {
   PROFILE3("NonlinearEl::evalInt");
 
+  ElmMats& elMat = static_cast<ElmMats&>(elmInt);
+
   // Evaluate the kinematic quantities, F and E, at this point
+  Matrix Bmat;
   Tensor F(nsd);
-  if (!this->kinematics(fe.N,fe.dNdX,X.x,F,E))
+  SymmTensor E(nsd);
+  if (!this->kinematics(elMat.vec.front(),fe.N,fe.dNdX,X.x,Bmat,F,E))
     return false;
 
   // Evaluate current tangent at this point, that is
   // the incremental constitutive matrix, Cmat, and
   // the 2nd Piola-Kirchhoff stress tensor, S
+  Matrix Cmat;
   SymmTensor S(nsd);
-  if (!this->formTangent(Cmat,S,X,F))
+  if (!this->formTangent(Cmat,S,X,F,E))
     return false;
 
   bool haveStrains = !E.isZero(1.0e-16);
@@ -102,33 +108,33 @@ bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt,
   if (eKm)
   {
     // Integrate the material stiffness matrix
+    Matrix& EM = elMat.A[eKm-1];
 #ifndef NO_FORTRAN
     // Using Fortran routines optimized for symmetric constitutive tensors
     PROFILE4("stiff_TL_");
     if (nsd == 3)
       if (fullCmat)
 	stiff_tl3d_(fe.N.size(),fe.detJxW,fe.dNdX.ptr(),F.ptr(),
-		    Cmat.ptr(),eKm->ptr());
+		    Cmat.ptr(),EM.ptr());
       else
 	stiff_tl3d_isoel_(fe.N.size(),fe.detJxW,fe.dNdX.ptr(),F.ptr(),
-			  Cmat(1,1),Cmat(1,2),Cmat(4,4),eKm->ptr());
+			  Cmat(1,1),Cmat(1,2),Cmat(4,4),EM.ptr());
     else if (nsd == 2)
       if (fullCmat)
 	stiff_tl2d_(fe.N.size(),fe.detJxW,fe.dNdX.ptr(),F.ptr(),
-		    Cmat.ptr(),eKm->ptr());
+		    Cmat.ptr(),EM.ptr());
       else
 	stiff_tl2d_isoel_(fe.N.size(),fe.detJxW,fe.dNdX.ptr(),F.ptr(),
-			  Cmat(1,1),Cmat(1,2),Cmat(3,3),eKm->ptr());
+			  Cmat(1,1),Cmat(1,2),Cmat(3,3),EM.ptr());
     else if (nsd == 1)
       for (a = 1; a <= fe.N.size(); a++)
 	for (b = 1; b <= fe.N.size(); b++)
-	  (*eKm)(a,b) += fe.dNdX(a,1)*F(1,1)*Cmat(1,1)*F(1,1)*fe.dNdX(b,1);
+	  EM(a,b) += fe.dNdX(a,1)*F(1,1)*Cmat(1,1)*F(1,1)*fe.dNdX(b,1);
 #else
     // This is too costly, but is basically what is done in the fortran routines
     PROFILE4("dNdX^t*F^t*C*F*dNdX");
     unsigned short int l, m, n;
     SymmTensor4 D(Cmat,nsd); // fourth-order material tensor
-    Matrix& EM = *eKm;
     for (a = 1; a <= fe.N.size(); a++)
       for (b = 1; b <= fe.N.size(); b++)
 	for (m = 1; m <= nsd; m++)
@@ -148,16 +154,16 @@ bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt,
 
   if (eKg && haveStrains)
     // Integrate the geometric stiffness matrix
-    this->formKG(*eKg,fe.N,fe.dNdX,X.x,S,fe.detJxW);
+    this->formKG(elMat.A[eKg-1],fe.N,fe.dNdX,X.x,S,fe.detJxW);
 
   if (eM)
     // Integrate the mass matrix
-    this->formMassMatrix(*eM,fe.N,X,fe.detJxW);
+    this->formMassMatrix(elMat.A[eM-1],fe.N,X,fe.detJxW);
 
   if (iS && haveStrains)
   {
     // Integrate the internal forces
-    Vector& ES = *iS;
+    Vector& ES = elMat.b[iS-1];
     for (a = 1; a <= fe.N.size(); a++)
       for (k = 1; k <= nsd; k++)
       {
@@ -171,9 +177,9 @@ bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt,
 
   if (eS)
     // Integrate the load vector due to gravitation and other body forces
-    this->formBodyForce(*eS,fe.N,X,fe.detJxW);
+    this->formBodyForce(elMat.b[eS-1],fe.N,X,fe.detJxW);
 
-  return this->getIntegralResult(elmInt);
+  return true;
 }
 
 
@@ -183,18 +189,20 @@ bool NonlinearElasticity::evalSol (Vector& s, const Vector&,
 {
   PROFILE3("NonlinearEl::evalSol");
 
+  // Extract element displacements
+  Vector eV;
   int ierr = 0;
-  if (eV && !primsol.front().empty())
+  if (!primsol.front().empty())
-    if ((ierr = utl::gather(MNPC,nsd,primsol.front(),*eV)))
+    if ((ierr = utl::gather(MNPC,nsd,primsol.front(),eV)))
     {
-      std::cerr <<" *** NonlinearElasticity::evalSol: Detected "
+      std::cerr <<" *** NonlinearElasticity::evalSol: Detected "<< ierr
-		<< ierr <<" node numbers out of range."<< std::endl;
+		<<" node numbers out of range."<< std::endl;
       return false;
     }
 
   // Evaluate the stress state at this point
   SymmTensor Sigma(nsd);
-  if (!this->formStressTensor(dNdX,X,Sigma))
+  if (!this->formStressTensor(eV,dNdX,X,Sigma))
     return false;
 
   // Congruence transformation to local coordinate system at current point
@@ -206,13 +214,13 @@ bool NonlinearElasticity::evalSol (Vector& s, const Vector&,
 }
 
 
-bool NonlinearElasticity::evalSol (Vector& s,
+bool NonlinearElasticity::evalSol (Vector& s, const Vector& eV,
 				   const Matrix& dNdX, const Vec3& X) const
 {
   PROFILE3("NonlinearEl::evalSol");
 
   SymmTensor Sigma(nsd);
-  if (!this->formStressTensor(dNdX,X,Sigma))
+  if (!this->formStressTensor(eV,dNdX,X,Sigma))
     return false;
 
   s = Sigma;
@@ -220,10 +228,11 @@ bool NonlinearElasticity::evalSol (Vector& s,
 }
 
 
-bool NonlinearElasticity::formStressTensor (const Matrix& dNdX, const Vec3& X,
+bool NonlinearElasticity::formStressTensor (const Vector& eV,
+					    const Matrix& dNdX, const Vec3& X,
 					    SymmTensor& S) const
 {
-  if (!eV || eV->empty())
+  if (eV.empty())
   {
     // Initial state (zero stresses)
     S.zero();
@@ -231,18 +240,22 @@ bool NonlinearElasticity::formStressTensor (const Matrix& dNdX, const Vec3& X,
   }
 
   // Evaluate the kinematic quantities, F and E, at this point
+  Matrix B;
   Tensor F(nsd);
-  if (!this->kinematics(Vector(),dNdX,X.x,F,E))
+  SymmTensor E(nsd);
+  if (!this->kinematics(eV,Vector(),dNdX,X.x,B,F,E))
     return false;
 
   // Evaluate the 2nd Piola-Kirchhoff stress tensor, S, at this point
+  Matrix Cmat;
   double U;
   return material->evaluate(Cmat,S,U,X,F,E,2);
 }
 
 
 bool NonlinearElasticity::formTangent (Matrix& Ctan, SymmTensor& S,
-				       const Vec3& X, const Tensor& F) const
+				       const Vec3& X, const Tensor& F,
+				       const SymmTensor& E) const
 {
   double U;
   return material->evaluate(Ctan,S,U,X,F,E, E.isZero(1.0e-16) ? 0 : 2);

Apps/FiniteDefElasticity/NonlinearElasticity.h

@@ -51,7 +51,7 @@ public:
   //! \param elmInt The local integral object to receive the contributions
   //! \param[in] fe Finite element data of current integration point
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe,
+  virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
 		       const Vec3& X) const;
 
   //! \brief Evaluates the secondary solution at a result point.
@@ -66,9 +66,11 @@ public:
 
   //! \brief Evaluates the finite element (FE) solution at an integration point.
   //! \param[out] s The FE stress values at current point
+  //! \param[in] eV Element solution vector
   //! \param[in] dNdX Basis function gradients at current point
   //! \param[in] X Cartesian coordinates of current point
-  virtual bool evalSol(Vector& s, const Matrix& dNdX, const Vec3& X) const;
+  virtual bool evalSol(Vector& s, const Vector& eV,
+		       const Matrix& dNdX, const Vec3& X) const;
 
 protected:
   //! \brief Forms tangential tensorial quantities needed by the evalInt method.
@@ -76,19 +78,22 @@ protected:
   //! \param[out] S 2nd Piola-Kirchhoff stress tensor at current point
   //! \param[in] X Cartesian coordinates of current integration point
   //! \param[in] F Deformation gradient at current integration point
+  //! \param[in] E Green-Lagrange strain tensor at current integration point
   virtual bool formTangent(Matrix& Ctan, SymmTensor& S,
-			   const Vec3& X, const Tensor& F) const;
+			   const Vec3& X, const Tensor& F,
+			   const SymmTensor& E) const;
 
   //! \brief Forms the 2nd Piola-Kirchhoff stress tensor.
+  //! \param[in] eV Element solution vector
   //! \param[in] dNdX Basis function gradients at current integration point
   //! \param[in] X Cartesian coordinates of current integration point
   //! \param[out] S 2nd Piola-Kirchhoff stress tensor at current point
-  virtual bool formStressTensor(const Matrix& dNdX, const Vec3& X,
+  virtual bool formStressTensor(const Vector& eV,
+				const Matrix& dNdX, const Vec3& X,
 				SymmTensor& S) const;
 
 protected:
-  bool        fullCmat; //!< If \e true, assume a full (but symmetric) C-matrix
+  bool fullCmat; //!< If \e true, assume a full (but symmetric) C-matrix
-  mutable SymmTensor E; //!< Green-Lagrange strain tensor
 };
 
 #endif

Apps/FiniteDefElasticity/NonlinearElasticityFbar.C

@@ -15,18 +15,97 @@
 #include "MaterialBase.h"
 #include "FiniteElement.h"
 #include "GaussQuadrature.h"
+#include "ElmMats.h"
+#include "ElmNorm.h"
 #include "Lagrange.h"
 #include "TimeDomain.h"
 #include "Tensor.h"
 
 
+/*!
+  \brief A struct with volumetric sampling point data.
+*/
+
+struct VolPtData
+{
+  double J;    //!< Determinant of current deformation gradient
+  Vector Nr;   //!< Basis function values (for axisymmetric problems)
+  Matrix dNdx; //!< Basis function gradients at current configuration
+  //! \brief Default constructor.
+  VolPtData() { J = 1.0; }
+};
+
+
+/*!
+  \brief Class containing internal data for an Fbar element.
+*/
+
+class FbarElmData
+{
+public:
+  //! \brief Default constructor.
+  FbarElmData() : pbar(0), scale(0.0), iP(0) {}
+  //! \brief Empty destructor.
+  virtual ~FbarElmData() {}
+
+  //! \brief Initializes the element data
+  bool init(unsigned short int nsd, size_t nPt)
+  {
+    pbar = ceil(pow((double)nPt,1.0/(double)nsd)-0.5);
+    if (pow((double)pbar,(double)nsd) == (double)nPt)
+      myVolData.resize(nPt);
+    else
+    {
+      pbar = 0;
+      std::cerr <<" *** FbarElmData::init: Invalid element, "<< nPt
+		<<" volumetric sampling points specified."<< std::endl;
+      return false;
+    }
+    scale = pbar > 1 ? 1.0/GaussQuadrature::getCoord(pbar)[pbar-1] : 1.0;
+    iP = 0;
+    return true;
+  }
+
+  int    pbar;  //!< Polynomial order of the internal volumetric data field
+  double scale; //!< Scaling factor for extrapolation from sampling points
+  size_t iP;    //!< Volumetric sampling point counter
+
+  std::vector<VolPtData> myVolData; //!< Volumetric sampling point data
+};
+
+
+/*!
+  \brief Class collecting element matrices associated with a Fbar FEM problem.
+*/
+
+class FbarMats : public FbarElmData, public ElmMats {};
+
+
+/*!
+  \brief Class representing integrated norm quantities for an Fbar element.
+*/
+
+class FbarNorm : public FbarElmData, public LocalIntegral
+{
+public:
+  //! \brief The constructor initializes the ElmNorm pointer.
+  FbarNorm(LocalIntegral& n) { myNorm = static_cast<ElmNorm*>(&n); }
+  //! \brief Alternative constructor allocating an internal ElmNorm.
+  FbarNorm(size_t n) { myNorm = new ElmNorm(n); }
+  //! \brief The destructor destroys the ElmNorm object.
+  virtual ~FbarNorm() { myNorm->destruct(); }
+
+  //! \brief Returns the LocalIntegral object to assemble into the global one.
+  virtual const LocalIntegral* ref() const { return myNorm; }
+
+  ElmNorm* myNorm; //!< Pointer to the actual element norms object
+};
+
+
 NonlinearElasticityFbar::NonlinearElasticityFbar (unsigned short int n,
 						  bool axS, int nvp)
   : NonlinearElasticityUL(n,axS), npt1(nvp)
 {
-  iP = 0;
-  pbar = 0;
-  scale = 1.0;
 }
 
 
@@ -39,48 +118,104 @@ void NonlinearElasticityFbar::print (std::ostream& os) const
 }
 
 
-bool NonlinearElasticityFbar::initElement (const std::vector<int>& MNPC,
+LocalIntegral* NonlinearElasticityFbar::getLocalIntegral (size_t nen, size_t,
-                                           const Vec3&, size_t nPt)
+							  bool neumann) const
 {
-  iP = 0;
+  ElmMats* result = new FbarMats;
-  pbar = ceil(pow((double)nPt,1.0/(double)nsd)-0.5);
+  switch (m_mode)
-  if (pow((double)pbar,(double)nsd) == (double)nPt)
-    myVolData.resize(nPt);
-  else
   {
-    pbar = 0;
+    case SIM::STATIC:
-    std::cerr <<" *** NonlinearElasticityFbar::initElement: Invalid element, "
+      result->rhsOnly = neumann;
-	      << nPt <<" volumetric sampling points specified."<< std::endl;
+      result->withLHS = !neumann;
-    return false;
+      result->resize(neumann?0:1,1);
-  }
+      break;
-  scale = pbar > 1 ? 1.0/GaussQuadrature::getCoord(pbar)[pbar-1] : 1.0;
-#if INT_DEBUG > 0
-  std::cout <<"\n\n *** Entering NonlinearElasticityFbar::initElement";
-  std::cout <<", nPt = "<< nPt <<", pbar = "<< pbar
-	    <<", scale = "<< scale << std::endl;
-#endif
 
-  return this->NonlinearElasticityUL::initElement(MNPC);
+    case SIM::DYNAMIC:
+      result->rhsOnly = neumann;
+      result->withLHS = !neumann;
+      result->resize(neumann?0:2,1);
+      break;
+
+    case SIM::VIBRATION:
+    case SIM::BUCKLING:
+      result->withLHS = true;
+      result->resize(2,0);
+      break;
+
+    case SIM::STIFF_ONLY:
+    case SIM::MASS_ONLY:
+      result->withLHS = true;
+      result->resize(1,0);
+      break;
+
+    case SIM::RHS_ONLY:
+      result->rhsOnly = true;
+      result->resize(neumann?0:1,1);
+      break;
+
+    case SIM::RECOVERY:
+      result->rhsOnly = true;
+      break;
+
+    default:
+      ;
+  }
+
+  for (size_t i = 0; i < result->A.size(); i++)
+    result->A[i].resize(nsd*nen,nsd*nen);
+
+  if (result->b.size())
+    result->b.front().resize(nsd*nen);
+
+  return result;
 }
 
 
-bool NonlinearElasticityFbar::reducedInt (const FiniteElement& fe,
+bool NonlinearElasticityFbar::initElement (const std::vector<int>& MNPC,
+                                           const Vec3&, size_t nPt,
+					   LocalIntegral& elmInt)
+{
+  FbarMats& fbar = static_cast<FbarMats&>(elmInt);
+
+#if INT_DEBUG > 0
+  std::cout <<"\n\n *** Entering NonlinearElasticityFbar::initElement";
+  std::cout <<", nPt = "<< nPt <<", pbar = "<< fbar.pbar
+	    <<", scale = "<< fbar.scale << std::endl;
+#endif
+  return fbar.init(nsd,nPt) && this->IntegrandBase::initElement(MNPC,elmInt);
+}
+
+
+bool NonlinearElasticityFbar::reducedInt (LocalIntegral& elmInt,
+					  const FiniteElement& fe,
 					  const Vec3& X) const
 {
+  FbarNorm* fbNorm = 0;
+  FbarElmData* fbar = dynamic_cast<FbarMats*>(&elmInt);
+  if (!fbar)
+  {
+    if ((fbNorm = dynamic_cast<FbarNorm*>(&elmInt)))
+      fbar = fbNorm;
+    else
+      return false;
+  }
+
 #if INT_DEBUG > 1
   std::cout <<"NonlinearElasticityFbar::u(red) = "<< fe.u;
   if (nsd > 1) std::cout <<" "<< fe.v;
   if (nsd > 2) std::cout <<" "<< fe.w;
-  std::cout <<" (iP="<< iP+1 <<")\n"
+  std::cout <<" (iP="<< fbar->iP+1 <<")\n"
 	    <<"NonlinearElasticityFbar::dNdX ="<< fe.dNdX;
 #endif
 
-  VolPtData& ptData = myVolData[iP++];
+  const Vector& eV = fbNorm ? fbNorm->myNorm->vec.front() : elmInt.vec.front();
+  VolPtData& ptData = fbar->myVolData[fbar->iP++];
 
   // Evaluate the deformation gradient, F, at current configuration
+  Matrix B;
   Tensor F(nDF);
   SymmTensor E(nsd,axiSymmetry);
-  if (!this->kinematics(fe.N,fe.dNdX,X.x,F,E))
+  if (!this->kinematics(eV,fe.N,fe.dNdX,X.x,B,F,E))
     return false;
 
   if (E.isZero(1.0e-16))
@@ -110,7 +245,7 @@ bool NonlinearElasticityFbar::reducedInt (const FiniteElement& fe,
     // Push-forward the basis function gradients to current configuration
     ptData.dNdx.multiply(fe.dNdX,Fi); // dNdx = dNdX * F^-1
     if (axiSymmetry && X.x > 0.0)
-      ptData.Nr = fe.N * (1.0/(X.x + eV->dot(fe.N,0,nsd)));
+      ptData.Nr = fe.N * (1.0/(X.x + eV.dot(fe.N,0,nsd)));
 
 #ifdef INT_DEBUG
     std::cout <<"NonlinearElasticityFbar::J = "<< ptData.J
@@ -120,7 +255,7 @@ bool NonlinearElasticityFbar::reducedInt (const FiniteElement& fe,
 #endif
   }
 #ifdef INT_DEBUG
-  if (iP == myVolData.size())
+  if (fbar->iP == fbar->myVolData.size())
     std::cout <<"NonlinearElasticityFbar: Volumetric sampling points completed."
 	      <<"\n"<< std::endl;
 #endif
@@ -326,7 +461,7 @@ static void getAmat3D (const Matrix& C, const SymmTensor& Sig, Matrix& A)
 }
 
 
-bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
+bool NonlinearElasticityFbar::evalInt (LocalIntegral& elmInt,
 				       const FiniteElement& fe,
 				       const TimeDomain& prm,
 				       const Vec3& X) const
@@ -338,10 +473,13 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
   std::cout <<"\nNonlinearElasticityFbar::dNdX ="<< fe.dNdX;
 #endif
 
+  FbarMats& fbar = static_cast<FbarMats&>(elmInt);
+
   // Evaluate the deformation gradient, F, at current configuration
+  Matrix B, dNdx;
   Tensor F(nDF);
   SymmTensor E(nsd,axiSymmetry);
-  if (!this->kinematics(fe.N,fe.dNdX,X.x,F,E))
+  if (!this->kinematics(fbar.vec.front(),fe.N,fe.dNdX,X.x,B,F,E))
     return false;
 
   double J, Jbar = 0.0;
@@ -380,23 +518,27 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
 
   // Axi-symmetric integration point volume; 2*pi*r*|J|*w
   double detJW = (axiSymmetry ? 2.0*M_PI*X.x : 1.0)*fe.detJxW*J;
-  double r = axiSymmetry ? X.x + eV->dot(fe.N,0,nsd) : 0.0;
+  double r = axiSymmetry ? X.x + fbar.vec.front().dot(fe.N,0,nsd) : 0.0;
 
-  if (myVolData.size() == 1)
+  Vector M;
+  Matrix dMdx;
+  if (fbar.myVolData.size() == 1)
   {
     // Only one volumetric sampling point (linear element)
-    Jbar = myVolData.front().J;
+    Jbar = fbar.myVolData.front().J;
-    dMdx = myVolData.front().dNdx;
+    dMdx = fbar.myVolData.front().dNdx;
-    if (axiSymmetry) M = myVolData.front().Nr;
+    if (axiSymmetry) M = fbar.myVolData.front().Nr;
   }
-  else if (!myVolData.empty())
+  else if (!fbar.myVolData.empty())
   {
     // Evaluate the Lagrange polynomials extrapolating the volumetric
     // sampling points (assuming a regular distribution over the element)
     Vector Nbar;
-    int pbar3 = nsd > 2 ? pbar : 0;
+    int pbar3 = nsd > 2 ? fbar.pbar : 0;
-    if (!Lagrange::computeBasis(Nbar,pbar,fe.xi*scale,pbar,fe.eta*scale,
+    if (!Lagrange::computeBasis(Nbar,
-				pbar3,fe.zeta*scale)) return false;
+				fbar.pbar,fe.xi*fbar.scale,
+				fbar.pbar,fe.eta*fbar.scale,
+				pbar3,fe.zeta*fbar.scale)) return false;
 #ifdef INT_DEBUG
     std::cout <<"NonlinearElasticityFbar::Nbar ="<< Nbar;
 #endif
@@ -408,10 +550,10 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
     if (axiSymmetry) M.resize(dNdx.rows(),true);
     for (size_t a = 0; a < Nbar.size(); a++)
     {
-      Jbar += myVolData[a].J*Nbar[a];
+      Jbar += fbar.myVolData[a].J*Nbar[a];
-      dMdx.add(myVolData[a].dNdx,myVolData[a].J*Nbar[a]);
+      dMdx.add(fbar.myVolData[a].dNdx,fbar.myVolData[a].J*Nbar[a]);
       if (axiSymmetry)
-	M.add(myVolData[a].Nr,myVolData[a].J*Nbar[a]);
+	M.add(fbar.myVolData[a].Nr,fbar.myVolData[a].J*Nbar[a]);
     }
     dMdx.multiply(1.0/Jbar);
     if (axiSymmetry) M /= Jbar;
@@ -434,6 +576,7 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
 #endif
 
   // Evaluate the constitutive relation (Jbar is dummy here)
+  Matrix Cmat;
   SymmTensor Sig(nsd,axiSymmetry);
   if (!material->evaluate(Cmat,Sig,Jbar,X,F,E,lHaveStrains,&prm))
     return false;
@@ -445,7 +588,7 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
   if (iS && lHaveStrains)
   {
     // Compute and accumulate contribution to internal force vector
-    Vector& ES = *iS;
+    Vector& ES = fbar.b[iS-1];
     size_t a, d;
     unsigned short int i, j;
     for (a = d = 1; a <= dNdx.rows(); a++)
@@ -465,6 +608,7 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
   if (eKm)
   {
     // Compute standard and modified discrete spatial gradient operators
+    Matrix G, Gbar;
     getGradOperator(axiSymmetry ? r : -1.0, fe.N, dNdx, G);
     getGradOperator(axiSymmetry ? 1.0 : -1.0, M, dMdx, Gbar);
 
@@ -496,78 +640,113 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
       }
 
 #ifdef INT_DEBUG
-  std::cout <<"NonlinearElasticityFbar::G ="<< G
+    std::cout <<"NonlinearElasticityFbar::G ="<< G
-	    <<"NonlinearElasticityFbar::Gbar ="<< Gbar
+	      <<"NonlinearElasticityFbar::Gbar ="<< Gbar
-	    <<"NonlinearElasticityFbar::A ="<< A
+	      <<"NonlinearElasticityFbar::A ="<< A
-	    <<"NonlinearElasticityFbar::Q ="<< Q;
+	      <<"NonlinearElasticityFbar::Q ="<< Q;
 #endif
 
     // Compute and accumulate contribution to tangent stiffness matrix.
 
     // First, standard (material and geometric) tangent stiffness terms
-    eKm->multiply(G,CB.multiply(A,G),true,false,true); // EK += G^T * A*G
+    Matrix CB; CB.multiply(A,G);
+    fbar.A[eKm-1].multiply(G,CB,true,false,true); // EK += G^T * A*G
 
     // Then, modify the tangent stiffness for the F-bar terms
-    Gbar -= G;
+    Gbar -= G; CB.multiply(Q,Gbar);
-    eKm->multiply(G,CB.multiply(Q,Gbar),true,false,true); // EK += G^T * Q*Gbar
+    fbar.A[eKm-1].multiply(G,CB,true,false,true); // EK += G^T * Q*Gbar
   }
 
   if (eM)
     // Integrate the mass matrix
-    this->formMassMatrix(*eM,fe.N,X,detJW);
+    this->formMassMatrix(fbar.A[eM-1],fe.N,X,detJW);
 
   if (eS)
     // Integrate the load vector due to gravitation and other body forces
-    this->formBodyForce(*eS,fe.N,X,detJW);
+    this->formBodyForce(fbar.b[eS-1],fe.N,X,detJW);
 
-  return this->getIntegralResult(elmInt);
+  return true;
 }
 
 
-
 NormBase* NonlinearElasticityFbar::getNormIntegrand (AnaSol*) const
 {
   return new ElasticityNormFbar(*const_cast<NonlinearElasticityFbar*>(this));
 }
 
 
-bool ElasticityNormFbar::reducedInt (const FiniteElement& fe,
+LocalIntegral* ElasticityNormFbar::getLocalIntegral (size_t nen, size_t iEl,
-				     const Vec3& X) const
+						     bool neumann) const
 {
-  return myProblem.reducedInt(fe,X);
+  LocalIntegral* result = this->NormBase::getLocalIntegral(nen,iEl,neumann);
+  if (result) return new FbarNorm(static_cast<ElmNorm&>(*result));
+
+  // Element norms are not requested, so allocate an internal object instead
+  // that will delete itself when invoking its destruct method.
+  return new FbarNorm(this->getNoFields());
 }
 
 
-bool ElasticityNormFbar::evalInt (LocalIntegral*& elmInt,
+bool ElasticityNormFbar::initElement (const std::vector<int>& MNPC,
+				      const Vec3&, size_t nPt,
+				      LocalIntegral& elmInt)
+{
+  NonlinearElasticityFbar& p = static_cast<NonlinearElasticityFbar&>(myProblem);
+  FbarNorm& fbar = static_cast<FbarNorm&>(elmInt);
+
+  return fbar.init(p.nsd,nPt) && this->NormBase::initElement(MNPC,*fbar.myNorm);
+}
+
+
+bool ElasticityNormFbar::initElementBou (const std::vector<int>& MNPC,
+					 LocalIntegral& elmInt)
+{
+  return myProblem.initElementBou(MNPC,*static_cast<FbarNorm&>(elmInt).myNorm);
+}
+
+
+bool ElasticityNormFbar::reducedInt (LocalIntegral& elmInt,
+				     const FiniteElement& fe,
+				     const Vec3& X) const
+{
+  return myProblem.reducedInt(elmInt,fe,X);
+}
+
+
+bool ElasticityNormFbar::evalInt (LocalIntegral& elmInt,
 				  const FiniteElement& fe,
 				  const TimeDomain& prm,
 				  const Vec3& X) const
 {
   size_t nsd = fe.dNdX.cols();
   NonlinearElasticityFbar& p = static_cast<NonlinearElasticityFbar&>(myProblem);
+  FbarNorm& fbar = static_cast<FbarNorm&>(elmInt);
 
   // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
+  Matrix B;
   Tensor F(p.nDF);
   SymmTensor E(p.nDF);
-  if (!p.kinematics(fe.N,fe.dNdX,X.x,F,E))
+  if (!p.kinematics(fbar.myNorm->vec.front(),fe.N,fe.dNdX,X.x,B,F,E))
     return false;
 
   double Jbar = 0.0;
-  if (p.myVolData.size() == 1)
+  if (fbar.myVolData.size() == 1)
-    Jbar = p.myVolData.front().J;
+    Jbar = fbar.myVolData.front().J;
-  else if (!p.myVolData.empty())
+  else if (!fbar.myVolData.empty())
   {
     // Evaluate the Lagrange polynomials extrapolating the volumetric
     // sampling points (assuming a regular distribution over the element)
     RealArray Nbar;
-    int pbar3 = nsd > 2 ? p.pbar : 0;
+    int pbar3 = nsd > 2 ? fbar.pbar : 0;
-    if (!Lagrange::computeBasis(Nbar,p.pbar,fe.xi*p.scale,p.pbar,fe.eta*p.scale,
+    if (!Lagrange::computeBasis(Nbar,
-				pbar3,fe.zeta*p.scale)) return false;
+				fbar.pbar,fe.xi*fbar.scale,
+				fbar.pbar,fe.eta*fbar.scale,
+				pbar3,fe.zeta*fbar.scale)) return false;
 
     // Compute modified deformation gradient determinant
     // by extrapolating the volume sampling points
     for (size_t a = 0; a < Nbar.size(); a++)
-      Jbar += p.myVolData[a].J*Nbar[a];
+      Jbar += fbar.myVolData[a].J*Nbar[a];
   }
   else
     Jbar = F.det();
@@ -578,15 +757,23 @@ bool ElasticityNormFbar::evalInt (LocalIntegral*& elmInt,
 
   // Compute the strain energy density, U(E) = Int_E (S:Eps) dEps
   // and the Cauchy stress tensor, sigma
-  double U = 0.0;
+  Matrix Cmat; double U = 0.0;
   SymmTensor sigma(nsd, p.isAxiSymmetric() || p.material->isPlaneStrain());
-  if (!p.material->evaluate(p.Cmat,sigma,U,X,Fbar,E,3,&prm,&F))
+  if (!p.material->evaluate(Cmat,sigma,U,X,Fbar,E,3,&prm,&F))
     return false;
 
   // Axi-symmetric integration point volume; 2*pi*r*|J|*w
   double detJW = p.isAxiSymmetric() ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW;
 
   // Integrate the norms
-  return ElasticityNormUL::evalInt(getElmNormBuffer(elmInt,6),
+  return ElasticityNormUL::evalInt(*fbar.myNorm,sigma,U,F.det(),detJW);
-				   sigma,U,F.det(),detJW);
+}
+
+
+bool ElasticityNormFbar::evalBou (LocalIntegral& elmInt,
+				  const FiniteElement& fe,
+				  const Vec3& X, const Vec3& normal) const
+{
+  return this->ElasticityNormUL::evalBou(*static_cast<FbarNorm&>(elmInt).myNorm,
+					 fe,X,normal);
 }

Apps/FiniteDefElasticity/NonlinearElasticityFbar.h

@@ -15,6 +15,8 @@
 #define _NONLINEAR_ELASTICITY_FBAR_H
 
 #include "NonlinearElasticityUL.h"
+#include "ElmMats.h"
+#include "ElmNorm.h"
 
 
 /*!
@@ -24,16 +26,6 @@
 
 class NonlinearElasticityFbar : public NonlinearElasticityUL
 {
-  //! \brief A struct with volumetric sampling point data.
-  struct VolPtData
-  {
-    double J;    //!< Determinant of current deformation gradient
-    Vector Nr;   //!< Basis function values (for axisymmetric problems)
-    Matrix dNdx; //!< Basis function gradients at current configuration
-    //! \brief Default constructor.
-    VolPtData() { J = 1.0; }
-  };
-
 public:
   //! \brief The default constructor invokes the parent class constructor.
   //! \param[in] n Number of spatial dimensions
@@ -47,23 +39,33 @@ public:
   //! \brief Prints out problem definition to the given output stream.
   virtual void print(std::ostream& os) const;
 
+  //! \brief Returns a local integral contribution object for the given element.
+  //! \param[in] nen Number of nodes on element
+  //! \param[in] neumann Whether or not we are assembling Neumann BC's
+  virtual LocalIntegral* getLocalIntegral(size_t nen, size_t,
+                                          bool neumann) const;
+
   //! \brief Initializes current element for numerical integration.
   //! \param[in] MNPC Matrix of nodal point correspondance for current element
   //! \param[in] nPt Number of volumetric sampling points in this element
+  //! \param elmInt The local integral object for current element
   virtual bool initElement(const std::vector<int>& MNPC,
-			   const Vec3&, size_t nPt);
+			   const Vec3&, size_t nPt,
+			   LocalIntegral& elmInt);
 
   //! \brief Evaluates volumetric sampling point data at an interior point.
+  //! \param elmInt The local integral object to receive the contributions
   //! \param[in] fe Finite element data of current point
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool reducedInt(const FiniteElement& fe, const Vec3& X) const;
+  virtual bool reducedInt(LocalIntegral& elmInt, const FiniteElement& fe,
+			  const Vec3& X) const;
 
   //! \brief Evaluates the integrand at an interior point.
   //! \param elmInt The local integral object to receive the contributions
   //! \param[in] fe Finite element data of current integration point
   //! \param[in] prm Nonlinear solution algorithm parameters
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe,
+  virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
 		       const TimeDomain& prm, const Vec3& X) const;
 
   //! \brief Returns a pointer to an Integrand for solution norm evaluation.
@@ -76,17 +78,7 @@ public:
   virtual int getIntegrandType() const { return 10+npt1; }
 
 private:
-  int    npt1;  //!< Number of volumetric sampling points in each direction
+  int npt1; //!< Number of volumetric sampling points in each direction
-  int    pbar;  //!< Polynomial order of the internal volumetric data field
-  double scale; //!< Scaling factor for extrapolation from sampling points
-
-  mutable std::vector<VolPtData> myVolData; //!< Volumetric sampling point data
-
-  mutable size_t iP;   //!< Volumetric sampling point counter
-  mutable Vector M;    //!< Modified basis function values
-  mutable Matrix dMdx; //!< Modified basis function gradients
-  mutable Matrix G;    //!< Discrete gradient operator
-  mutable Matrix Gbar; //!< Modified discrete gradient operator
 
   friend class ElasticityNormFbar;
 };
@@ -105,19 +97,50 @@ public:
   //! \brief Empty destructor.
   virtual ~ElasticityNormFbar() {}
 
+  //! \brief Returns a local integral contribution object for the given element.
+  //! \param[in] nen Number of nodes on element
+  //! \param[in] iEl Global element number
+  //! \param[in] neumann Whether or not we are assembling Neumann BC's
+  virtual LocalIntegral* getLocalIntegral(size_t nen, size_t iEl,
+                                          bool neumann) const;
+
+  //! \brief Initializes current element for numerical integration.
+  //! \param[in] MNPC Matrix of nodal point correspondance for current element
+  //! \param[in] nPt Number of volumetric sampling points in this element
+  //! \param elmInt The local integral object for current element
+  virtual bool initElement(const std::vector<int>& MNPC,
+			   const Vec3&, size_t nPt,
+			   LocalIntegral& elmInt);
+
+  //! \brief Initializes current element for boundary integration.
+  //! \param[in] MNPC Matrix of nodal point correspondance for current element
+  //! \param elmInt The local integral object for current element
+  virtual bool initElementBou(const std::vector<int>& MNPC,
+			      LocalIntegral& elmInt);
+
   //! \brief Evaluates volumetric sampling point data at an interior point.
+  //! \param elmInt The local integral object to receive the contributions
   //! \param[in] fe Finite element data of current point
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool reducedInt(const FiniteElement& fe, const Vec3& X) const;
+  virtual bool reducedInt(LocalIntegral& elmInt, const FiniteElement& fe,
+			  const Vec3& X) const;
 
   //! \brief Evaluates the integrand at an interior point.
   //! \param elmInt The local integral object to receive the contributions
   //! \param[in] fe Finite element data of current integration point
   //! \param[in] prm Nonlinear solution algorithm parameters
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe,
+  virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
 		       const TimeDomain& prm, const Vec3& X) const;
 
+  //! \brief Evaluates the integrand at a boundary point.
+  //! \param elmInt The local integral object to receive the contributions
+  //! \param[in] fe Finite element data of current integration point
+  //! \param[in] X Cartesian coordinates of current integration point
+  //! \param[in] normal Boundary normal vector at current integration point
+  virtual bool evalBou(LocalIntegral& elmInt, const FiniteElement& fe,
+		       const Vec3& X, const Vec3& normal) const;
+
   //! \brief Returns which integrand to be used.
   virtual int getIntegrandType() const { return myProblem.getIntegrandType(); }
 };

Apps/FiniteDefElasticity/NonlinearElasticityTL.C

@@ -23,6 +23,10 @@
 #define epsR 1.0e-16
 #endif
 
+#ifdef USE_OPENMP
+#include <omp.h>
+#endif
+
 
 NonlinearElasticityTL::NonlinearElasticityTL (unsigned short int n, bool axS)
   : Elasticity(n,axS)
@@ -135,7 +139,7 @@ bool NonlinearElasticityTL::evalInt (LocalIntegral& elmInt,
   Matrix Bmat;
   Tensor F(nDF);
   SymmTensor E(nsd,axiSymmetry), S(nsd,axiSymmetry);
-  if (!this->kinematics(elMat.vec,fe.N,fe.dNdX,X.x,Bmat,F,E))
+  if (!this->kinematics(elMat.vec.front(),fe.N,fe.dNdX,X.x,Bmat,F,E))
     return false;
 
   // Evaluate the constitutive relation
@@ -206,7 +210,10 @@ bool NonlinearElasticityTL::evalBou (LocalIntegral& elmInt,
   Vec3 T = (*tracFld)(X,normal);
 
   // Store the traction value for vizualization
-  if (!T.isZero()) tracVal[X] = T;
+#ifdef USE_OPENMP
+  if (omp_get_max_threads() == 1)
+#endif
+    if (!T.isZero()) tracVal[X] = T;
 
   // Check for with-rotated pressure load
   unsigned short int i, j;
@@ -216,7 +223,8 @@ bool NonlinearElasticityTL::evalBou (LocalIntegral& elmInt,
     Matrix B;
     Tensor F(nDF);
     SymmTensor dummy(0);
-    if (!this->kinematics(elMat.vec,fe.N,fe.dNdX,X.x,B,F,dummy)) return false;
+    if (!this->kinematics(elMat.vec.front(),fe.N,fe.dNdX,X.x,B,F,dummy))
+      return false;
 
     // Compute its inverse and determinant, J
     double J = F.inverse();
@@ -244,12 +252,12 @@ bool NonlinearElasticityTL::evalBou (LocalIntegral& elmInt,
 }
 
 
-bool NonlinearElasticityTL::kinematics (const Vectors& eV,
+bool NonlinearElasticityTL::kinematics (const Vector& eV,
 					const Vector& N, const Matrix& dNdX,
 					double r, Matrix& Bmat, Tensor& F,
 					SymmTensor& E) const
 {
-  if (eV.empty() || eV.front().empty())
+  if (eV.empty())
   {
     // Initial state, unit deformation gradient and linear B-matrix
     F = 1.0;
@@ -262,7 +270,7 @@ bool NonlinearElasticityTL::kinematics (const Vectors& eV,
 
   const size_t nenod = dNdX.rows();
   const size_t nstrc = axiSymmetry ? 4 : nsd*(nsd+1)/2;
-  if (eV.front().size() != nenod*nsd || dNdX.cols() < nsd)
+  if (eV.size() != nenod*nsd || dNdX.cols() < nsd)
   {
     std::cerr <<" *** NonlinearElasticityTL::kinematics: Invalid dimension,"
 	      <<" dNdX("<< nenod <<","<< dNdX.cols() <<")"<< std::endl;
@@ -273,10 +281,10 @@ bool NonlinearElasticityTL::kinematics (const Vectors& eV,
   // and the Green-Lagrange strains, E_ij = 0.5(F_ij+F_ji+F_ki*F_kj).
 
   // Notice that the matrix multiplication method used here treats the
-  // element displacement vector, *eV, as a matrix whose number of columns
+  // element displacement vector, eV, as a matrix whose number of columns
   // equals the number of rows in the matrix dNdX.
   Matrix dUdX;
-  if (!dUdX.multiplyMat(eV.front(),dNdX)) // dUdX = Grad{u} = eV*dNd
+  if (!dUdX.multiplyMat(eV,dNdX)) // dUdX = Grad{u} = eV*dNd
     return false;
 
   unsigned short int i, j, k;
@@ -304,7 +312,7 @@ bool NonlinearElasticityTL::kinematics (const Vectors& eV,
   // Add the unit tensor to F to form the deformation gradient
   F += 1.0;
   // Add the dU/r term to the F(3,3)-term for axisymmetric problems
-  if (axiSymmetry && r > epsR) F(3,3) += eV.front().dot(N,0,nsd)/r;
+  if (axiSymmetry && r > epsR) F(3,3) += eV.dot(N,0,nsd)/r;
 
 #ifdef INT_DEBUG
   std::cout <<"NonlinearElasticityTL::F =\n"<< F;

Apps/FiniteDefElasticity/NonlinearElasticityTL.h

@@ -69,7 +69,7 @@ public:
 
 protected:
   //! \brief Calculates some kinematic quantities at current point.
-  //! \param[in] eV Element solution vectors
+  //! \param[in] eV Element solution vector
   //! \param[in] N Basis function values at current point
   //! \param[in] dNdX Basis function gradients at current point
   //! \param[in] r Radial coordinate of current point
@@ -80,7 +80,7 @@ protected:
   //! \details The deformation gradient \b F and the nonlinear
   //! strain-displacement matrix \b B are established.
   //! The B-matrix is formed only when the variable \a formB is true.
-  virtual bool kinematics(const Vectors& eV,
+  virtual bool kinematics(const Vector& eV,
 			  const Vector& N, const Matrix& dNdX, double r,
 			  Matrix& Bmat, Tensor& F, SymmTensor& E) const;
 

Apps/FiniteDefElasticity/NonlinearElasticityUL.C

@@ -25,6 +25,10 @@
 #define epsR 1.0e-16
 #endif
 
+#ifdef USE_OPENMP
+#include <omp.h>
+#endif
+
 
 NonlinearElasticityUL::NonlinearElasticityUL (unsigned short int n,
 					      bool axS, char lop)
@@ -46,35 +50,20 @@ void NonlinearElasticityUL::print (std::ostream& os) const
 
 void NonlinearElasticityUL::setMode (SIM::SolutionMode mode)
 {
-  if (!myMats) return;
+  m_mode = mode;
-
-  myMats->rhsOnly = false;
   eM = eKm = eKg = 0;
-  iS = eS  = eV  = 0;
+  eS = iS  = 0;
 
   switch (mode)
     {
     case SIM::STATIC:
-      myMats->resize(1,1);
+    case SIM::RHS_ONLY:
-      eKm = &myMats->A[0];
+      eKm = eKg = 1;
-      eKg = &myMats->A[0];
       break;
 
     case SIM::DYNAMIC:
-      myMats->resize(2,1);
+      eKm = eKg = 1;
-      eKm = &myMats->A[0];
+      eM  = 2;
-      eKg = &myMats->A[0];
-      eM  = &myMats->A[1];
-      break;
-
-    case SIM::RHS_ONLY:
-      if (myMats->A.empty())
-	myMats->resize(1,1);
-      else
-        myMats->b.resize(1);
-      eKm = &myMats->A[0];
-      eKg = &myMats->A[0];
-      myMats->rhsOnly = true;
       break;
 
     default:
@@ -82,15 +71,61 @@ void NonlinearElasticityUL::setMode (SIM::SolutionMode mode)
       return;
     }
 
-  // We always need the force+displacement vectors in nonlinear simulations
+  // We always need the force vectors in nonlinear simulations
-  iS = &myMats->b[0];
+  iS = eS = 1;
-  eS = &myMats->b[0];
-  mySols.resize(1);
-  eV = &mySols[0];
   tracVal.clear();
 }
 
 
+LocalIntegral* NonlinearElasticityUL::getLocalIntegral (size_t nen, size_t,
+                                                        bool neumann) const
+{
+  ElmMats* result = new ElmMats;
+  switch (m_mode)
+  {
+    case SIM::STATIC:
+      result->withLHS = !neumann;
+    case SIM::RHS_ONLY:
+      result->rhsOnly = neumann;
+      result->resize(neumann?0:1,1);
+      break;
+
+    case SIM::DYNAMIC:
+      result->withLHS = !neumann;
+      result->rhsOnly = neumann;
+      result->resize(neumann?0:2,1);
+      break;
+
+    case SIM::VIBRATION:
+    case SIM::BUCKLING:
+      result->withLHS = true;
+      result->resize(2,0);
+      break;
+
+    case SIM::STIFF_ONLY:
+    case SIM::MASS_ONLY:
+      result->withLHS = true;
+      result->resize(1,0);
+      break;
+
+    case SIM::RECOVERY:
+      result->rhsOnly = true;
+      break;
+
+    default:
+      ;
+  }
+
+  for (size_t i = 0; i < result->A.size(); i++)
+    result->A[i].resize(nsd*nen,nsd*nen);
+
+  if (result->b.size())
+    result->b.front().resize(nsd*nen);
+
+  return result;
+}
+
+
 void NonlinearElasticityUL::initIntegration (const TimeDomain& prm)
 {
   if (material)
@@ -108,15 +143,18 @@ void NonlinearElasticityUL::initResultPoints (double lambda)
 }
 
 
-bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt,
+bool NonlinearElasticityUL::evalInt (LocalIntegral& elmInt,
 				     const FiniteElement& fe,
 				     const TimeDomain& prm,
 				     const Vec3& X) const
 {
+  ElmMats& elMat = static_cast<ElmMats&>(elmInt);
+
   // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
+  Matrix Bmat, dNdx;
   Tensor F(nDF);
   SymmTensor E(nsd,axiSymmetry);
-  if (!this->kinematics(fe.N,fe.dNdX,X.x,F,E))
+  if (!this->kinematics(elMat.vec.front(),fe.N,fe.dNdX,X.x,Bmat,F,E))
     return false;
 
   // Axi-symmetric integration point volume; 2*pi*r*|J|*w
@@ -149,11 +187,11 @@ bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt,
       // Compute the small-deformation strain-displacement matrix B from dNdx
       if (axiSymmetry)
       {
-	r += eV->dot(fe.N,0,nsd);
+	r += elMat.vec.front().dot(fe.N,0,nsd);
-	this->formBmatrix(fe.N,dNdx,r);
+	this->formBmatrix(Bmat,fe.N,dNdx,r);
       }
       else
-	this->formBmatrix(dNdx);
+	this->formBmatrix(Bmat,dNdx);
 
 #ifdef INT_DEBUG
       std::cout <<"NonlinearElasticityUL::dNdx ="<< dNdx;
@@ -165,12 +203,13 @@ bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt,
   {
     // Initial state, no deformation yet
     if (axiSymmetry)
-      this->formBmatrix(fe.N,fe.dNdX,r);
+      this->formBmatrix(Bmat,fe.N,fe.dNdX,r);
     else
-      this->formBmatrix(fe.dNdX);
+      this->formBmatrix(Bmat,fe.dNdX);
   }
 
   // Evaluate the constitutive relation
+  Matrix Cmat;
   SymmTensor sigma(nsd,axiSymmetry);
   if (eKm || eKg || iS)
   {
@@ -182,35 +221,36 @@ bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt,
   if (eKm)
   {
     // Integrate the material stiffness matrix
+    Matrix CB;
     CB.multiply(Cmat,Bmat).multiply(detJW); // CB = C*B*|J|*w
-    eKm->multiply(Bmat,CB,true,false,true); // EK += B^T * CB
+    elMat.A[eKm-1].multiply(Bmat,CB,true,false,true); // EK += B^T * CB
   }
 
   if (eKg && lHaveStrains)
     // Integrate the geometric stiffness matrix
-    this->formKG(*eKg,fe.N,dNdx,r,sigma,detJW);
+    this->formKG(elMat.A[eKg-1],fe.N,dNdx,r,sigma,detJW);
 
   if (eM)
     // Integrate the mass matrix
-    this->formMassMatrix(*eM,fe.N,X,detJW);
+    this->formMassMatrix(elMat.A[eM-1],fe.N,X,detJW);
 
   if (iS && lHaveStrains)
   {
     // Integrate the internal forces
     sigma *= -detJW;
-    if (!Bmat.multiply(sigma,*iS,true,true)) // ES -= B^T*sigma
+    if (!Bmat.multiply(sigma,elMat.b[iS-1],true,true)) // ES -= B^T*sigma
       return false;
   }
 
   if (eS)
     // Integrate the load vector due to gravitation and other body forces
-    this->formBodyForce(*eS,fe.N,X,detJW);
+    this->formBodyForce(elMat.b[eS-1],fe.N,X,detJW);
 
-  return this->getIntegralResult(elmInt);
+  return true;
 }
 
 
-bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt,
+bool NonlinearElasticityUL::evalBou (LocalIntegral& elmInt,
 				     const FiniteElement& fe,
 				     const Vec3& X, const Vec3& normal) const
 {
@@ -231,7 +271,10 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt,
   Vec3 T = (*tracFld)(X,normal);
 
   // Store the traction value for vizualization
-  if (!T.isZero()) tracVal[X] = T;
+#ifdef USE_OPENMP
+  if (omp_get_max_threads() == 1)
+#endif
+    if (!T.isZero()) tracVal[X] = T;
 
   // Axi-symmetric integration point volume; 2*pi*r*|J|*w
   double detJW = axiSymmetry ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW;
@@ -241,7 +284,8 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt,
   {
     // Compute the deformation gradient, F
     Tensor F(nDF);
-    if (!this->formDefGradient(fe.N,fe.dNdX,X.x,F)) return false;
+    if (!this->formDefGradient(elmInt.vec.front(),fe.N,fe.dNdX,X.x,F))
+      return false;
 
     // Check for with-rotated pressure load
     if (tracFld->isNormalPressure())
@@ -264,29 +308,31 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt,
   }
 
   // Integrate the force vector
-  Vector& ES = *eS;
+  Vector& ES = static_cast<ElmMats&>(elmInt).b[eS-1];
   for (size_t a = 1; a <= fe.N.size(); a++)
     for (i = 1; i <= nsd; i++)
       ES(nsd*(a-1)+i) += T[i-1]*fe.N(a)*detJW;
 
-  return this->getIntegralResult(elmInt);
+  return true;
 }
 
 
-bool NonlinearElasticityUL::formDefGradient (const Vector& N,
+bool NonlinearElasticityUL::formDefGradient (const Vector& eV, const Vector& N,
 					     const Matrix& dNdX, double r,
 					     Tensor& F) const
 {
-  static SymmTensor dummy(0);
+  Matrix B;
-  return this->kinematics(N,dNdX,r,F,dummy);
+  SymmTensor dummy(0);
+  return this->kinematics(eV,N,dNdX,r,B,F,dummy);
 }
 
 
-bool NonlinearElasticityUL::kinematics (const Vector& N, const Matrix& dNdX,
+bool NonlinearElasticityUL::kinematics (const Vector& eV,
-					double r, Tensor& F,
+					const Vector& N, const Matrix& dNdX,
+					double r, Matrix&, Tensor& F,
 					SymmTensor& E) const
 {
-  if (!eV || eV->empty())
+  if (eV.empty())
   {
     // Initial state, unit deformation gradient and zero strains
     F = 1.0;
@@ -295,7 +341,7 @@ bool NonlinearElasticityUL::kinematics (const Vector& N, const Matrix& dNdX,
   }
 
   const size_t nenod = dNdX.rows();
-  if (eV->size() != nenod*nsd || dNdX.cols() < nsd)
+  if (eV.size() != nenod*nsd || dNdX.cols() < nsd)
   {
     std::cerr <<" *** NonlinearElasticityUL::kinematics: Invalid dimension,"
 	      <<" dNdX("<< nenod <<","<< dNdX.cols() <<")"<< std::endl;
@@ -309,7 +355,7 @@ bool NonlinearElasticityUL::kinematics (const Vector& N, const Matrix& dNdX,
   // element displacement vector, *eV, as a matrix whose number of columns
   // equals the number of rows in the matrix dNdX.
   Matrix dUdX;
-  if (!dUdX.multiplyMat(*eV,dNdX)) // dUdX = Grad{u} = eV*dNdX
+  if (!dUdX.multiplyMat(eV,dNdX)) // dUdX = Grad{u} = eV*dNdX
     return false;
 
   unsigned short int i, j, k;
@@ -337,10 +383,10 @@ bool NonlinearElasticityUL::kinematics (const Vector& N, const Matrix& dNdX,
   // Add the unit tensor to F to form the deformation gradient
   F += 1.0;
   // Add the dU/r term to the F(3,3)-term for axisymmetric problems
-  if (axiSymmetry && r > epsR) F(3,3) += eV->dot(N,0,nsd)/r;
+  if (axiSymmetry && r > epsR) F(3,3) += eV.dot(N,0,nsd)/r;
 
 #ifdef INT_DEBUG
-  std::cout <<"NonlinearElasticityUL::eV ="<< *eV;
+  std::cout <<"NonlinearElasticityUL::eV ="<< eV.front();
   std::cout <<"NonlinearElasticityUL::F =\n"<< F;
 #endif
 
@@ -361,7 +407,7 @@ void ElasticityNormUL::initIntegration (const TimeDomain& prm)
 }
 
 
-bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt,
+bool ElasticityNormUL::evalInt (LocalIntegral& elmInt,
 				const FiniteElement& fe,
 				const TimeDomain& prm,
 				const Vec3& X) const
@@ -369,23 +415,24 @@ bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt,
   NonlinearElasticityUL& ulp = static_cast<NonlinearElasticityUL&>(myProblem);
 
   // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
+  Matrix B;
   Tensor F(ulp.nDF);
   SymmTensor E(ulp.nDF);
-  if (!ulp.kinematics(fe.N,fe.dNdX,X.x,F,E))
+  if (!ulp.kinematics(elmInt.vec.front(),fe.N,fe.dNdX,X.x,B,F,E))
     return false;
 
   // Compute the strain energy density, U(E) = Int_E (S:Eps) dEps
   // and the Cauchy stress tensor, sigma
-  double U = 0.0;
+  Matrix Cmat; double U = 0.0;
   SymmTensor sigma(E.dim(),ulp.isAxiSymmetric()||ulp.material->isPlaneStrain());
-  if (!ulp.material->evaluate(ulp.Cmat,sigma,U,X,F,E,3,&prm))
+  if (!ulp.material->evaluate(Cmat,sigma,U,X,F,E,3,&prm))
     return false;
 
   // Axi-symmetric integration point volume; 2*pi*r*|J|*w
   double detJW = ulp.isAxiSymmetric() ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW;
 
   // Integrate the norms
-  return evalInt(getElmNormBuffer(elmInt,6),sigma,U,F.det(),detJW);
+  return evalInt(static_cast<ElmNorm&>(elmInt),sigma,U,F.det(),detJW);
 }
 
 
@@ -414,7 +461,7 @@ bool ElasticityNormUL::evalInt (ElmNorm& pnorm, const SymmTensor& S,
 }
 
 
-bool ElasticityNormUL::evalBou (LocalIntegral*& elmInt,
+bool ElasticityNormUL::evalBou (LocalIntegral& elmInt,
 				const FiniteElement& fe,
 				const Vec3& X, const Vec3& normal) const
 {
@@ -424,7 +471,7 @@ bool ElasticityNormUL::evalBou (LocalIntegral*& elmInt,
   // Evaluate the current surface traction
   Vec3 t = ulp.getTraction(X,normal);
   // Evaluate the current displacement field
-  Vec3 u = ulp.evalSol(fe.N);
+  Vec3 u = ulp.evalSol(elmInt.vec.front(),fe.N);
 
   // Integrate the external energy (path integral)
   if (iP < Ux.size())
@@ -445,6 +492,6 @@ bool ElasticityNormUL::evalBou (LocalIntegral*& elmInt,
   // Axi-symmetric integration point volume; 2*pi*r*|J|*w
   double detJW = ulp.isAxiSymmetric() ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW;
 
-  getElmNormBuffer(elmInt)[1] += Ux[iP++]*detJW;
+  static_cast<ElmNorm&>(elmInt)[1] += Ux[iP++]*detJW;
   return true;
 }

Apps/FiniteDefElasticity/NonlinearElasticityUL.h

@@ -16,6 +16,8 @@
 
 #include "Elasticity.h"
 
+class ElmNorm;
+
 
 /*!
   \brief Class representing the integrand of the nonlinear elasticity problem.
@@ -45,6 +47,12 @@ public:
   //! \param[in] mode The solution mode to use
   virtual void setMode(SIM::SolutionMode mode);
 
+  //! \brief Returns a local integral container for the given element.
+  //! \param[in] nen Number of DOFs on element
+  //! \param[in] neumann Whether or not we are assembling Neumann BC's
+  virtual LocalIntegral* getLocalIntegral(size_t nen, size_t,
+					  bool neumann) const;
+
   //! \brief Initializes the integrand for a new integration loop.
   //! \param[in] prm Nonlinear solution algorithm parameters
   virtual void initIntegration(const TimeDomain& prm);
@@ -57,7 +65,7 @@ public:
   //! \param[in] fe Finite element data of current integration point
   //! \param[in] prm Nonlinear solution algorithm parameters
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe,
+  virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
 		       const TimeDomain& prm, const Vec3& X) const;
 
   //! \brief Evaluates the integrand at a boundary point.
@@ -70,7 +78,7 @@ public:
   //! possibly with-rotated traction fields (non-conservative loads).
   //! For uni-directional (conservative) loads, it is similar to the
   //! \a LinearElasticity::evalBou method.
-  virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe,
+  virtual bool evalBou(LocalIntegral& elmInt, const FiniteElement& fe,
                        const Vec3& X, const Vec3& normal) const;
 
   //! \brief Returns a pointer to an Integrand for solution norm evaluation.
@@ -80,26 +88,26 @@ public:
   virtual NormBase* getNormIntegrand(AnaSol* = 0) const;
 
   //! \brief Calculates some kinematic quantities at current point.
+  //! \param[in] eV Element solution vector
   //! \param[in] N Basis function values at current point
   //! \param[in] dNdX Basis function gradients at current point
   //! \param[in] r Radial coordinate of current point
   //! \param[out] F Deformation gradient at current point
+  //! \param[out] B The strain-displacement matrix at current point
   //! \param[out] E Green-Lagrange strain tensor at current point
-  virtual bool kinematics(const Vector& N, const Matrix& dNdX, double r,
+  virtual bool kinematics(const Vector& eV,
-			  Tensor& F, SymmTensor& E) const;
+                          const Vector& N, const Matrix& dNdX, double r,
+			  Matrix& B, Tensor& F, SymmTensor& E) const;
 
 protected:
   //! \brief Calculates the deformation gradient at current point.
+  //! \param[in] eV Element solution vector
   //! \param[in] N Basis function values at current point
   //! \param[in] dNdX Basis function gradients at current point
   //! \param[in] r Radial coordinate of current point
   //! \param[out] F Deformation gradient at current point
-  bool formDefGradient(const Vector& N, const Matrix& dNdX, double r,
+  bool formDefGradient(const Vector& eV, const Vector& N,
-		       Tensor& F) const;
+		       const Matrix& dNdX, double r, Tensor& F) const;
-
-protected:
-  mutable Matrix dNdx; //!< Basis function gradients in current configuration
-  mutable Matrix CB;   //!< Result of the matrix-matrix product C*B
 
 private:
   char loadOp; //!< Load option
@@ -135,7 +143,7 @@ public:
   //! \param[in] fe Finite element data of current integration point
   //! \param[in] prm Nonlinear solution algorithm parameters
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe,
+  virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
 		       const TimeDomain& prm, const Vec3& X) const;
 
   //! \brief Evaluates the integrand at a boundary point.
@@ -143,7 +151,7 @@ public:
   //! \param[in] fe Finite element data of current integration point
   //! \param[in] X Cartesian coordinates of current integration point
   //! \param[in] normal Boundary normal vector at current integration point
-  virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe,
+  virtual bool evalBou(LocalIntegral& elmInt, const FiniteElement& fe,
 		       const Vec3& X, const Vec3& normal) const;
 
   //! \brief Returns the number of norm quantities.

Apps/FiniteDefElasticity/NonlinearElasticityULMX.C

@@ -15,8 +15,10 @@
 #include "MaterialBase.h"
 #include "FiniteElement.h"
 #include "ElmMats.h"
+#include "ElmNorm.h"
 #include "TimeDomain.h"
 #include "Utilities.h"
+#include "Tensor.h"
 #include "Vec3Oper.h"
 
 #ifdef USE_FTNMAT
@@ -42,13 +44,102 @@ extern "C" {
 #endif
 
 
+/*!
+  \brief A struct with integration point data needed by \a finalizeElement.
+*/
+
+struct ItgPtData
+{
+  Tensor F;     //!< Deformation gradient at current step/iteration
+  Tensor Fp;    //!< Deformation gradient at previous (converged) step
+  Vector Nr;    //!< Basis function values (for axisymmetric problems)
+  Matrix dNdx;  //!< Basis function gradients at current configuration
+  Vector Phi;   //!< Internal modes for the pressure/volumetric-change fields
+  Vec4   X;     //!< Cartesian coordinates of current integration point
+  double detJW; //!< Jacobian determinant times integration point weight
+  //! \brief Default constructor.
+  ItgPtData() : F(3), Fp(3) { detJW = 0.0; }
+};
+
+
+/*!
+  \brief Class containing internal data for a mixed element.
+*/
+
+class MxElmData
+{
+public:
+  //! \brief Default constructor.
+  MxElmData() : Hh(0), iP(0) {}
+  //! \brief Empty destructor.
+  virtual ~MxElmData() {}
+
+  //! \brief Initializes the element data.
+  bool init(const Vec3& Xc, size_t nPt, int p, unsigned short int nsd)
+  {
+    iP = 0;
+    X0 = Xc;
+    myData.resize(nPt);
+    if (Hh)
+    {
+      size_t nPm = utl::Pascal(p,nsd);
+      Hh->resize(nPm,nPm,true);
+      return true;
+    }
+    std::cerr <<" *** MxElmData::init: No Hh matrix defined."<< std::endl;
+    return false;
+  }
+
+  Vec3    X0; //!< Cartesian coordinates of the element center
+  Matrix* Hh; //!< Pointer to element Hh-matrix associated with pressure modes
+  size_t  iP; //!< Local integration point counter
+
+  std::vector<ItgPtData> myData; //!< Local integration point data
+};
+
+
+/*!
+  \brief Class collecting element matrices associated with a mixed FEM problem.
+*/
+
+class MxMats : public MxElmData, public ElmMats {};
+
+
+/*!
+  \brief Class representing integrated norm quantities for a mixed element.
+*/
+
+class MxNorm : public MxElmData, public LocalIntegral
+{
+public:
+  //! \brief The constructor initializes the ElmNorm pointer.
+  MxNorm(LocalIntegral& n)
+  {
+    myNorm = static_cast<ElmNorm*>(&n);
+    Hh = new Matrix();
+  }
+
+  //! \brief Alternative constructor allocating an internal ElmNorm.
+  MxNorm(size_t n)
+  {
+    myNorm = new ElmNorm(n);
+    Hh = new Matrix();
+  }
+
+  //! \brief The destructor destroys the ElmNorm object.
+  virtual ~MxNorm() { myNorm->destruct(); }
+
+  //! \brief Returns the LocalIntegral object to assemble into the global one.
+  virtual const LocalIntegral* ref() const { return myNorm; }
+
+  ElmNorm* myNorm; //!< Pointer to the actual element norms object
+};
+
+
 NonlinearElasticityULMX::NonlinearElasticityULMX (unsigned short int n,
 						  bool axS, int pp)
-  : NonlinearElasticityUL(n,axS)
+  : NonlinearElasticityUL(n,axS), p(pp)
 {
-  p = pp;
-  iP = 0;
-
   // Both the current and previous solutions are needed
   primsol.resize(2);
 }
@@ -63,106 +154,90 @@ void NonlinearElasticityULMX::print (std::ostream& os) const
 }
 
 
-void NonlinearElasticityULMX::setMode (SIM::SolutionMode mode)
+LocalIntegral* NonlinearElasticityULMX::getLocalIntegral (size_t nen, size_t,
+							  bool neumann) const
 {
-  if (!myMats) return;
+  MxMats* result = new MxMats;
 
-  myMats->rhsOnly = false;
+  switch (m_mode)
-  Hh = eM = eKm = eKg = 0;
-  iS = eS = eV  = 0;
-
-  switch (mode)
     {
     case SIM::STATIC:
-      myMats->resize(2,1);
+      result->withLHS = !neumann;
-      eKm = &myMats->A[0];
+    case SIM::RHS_ONLY:
-      eKg = &myMats->A[0];
+      result->rhsOnly = neumann;
-      Hh  = &myMats->A[1];
+      result->resize(neumann?1:2,1);
       break;
 
     case SIM::DYNAMIC:
-      myMats->resize(3,1);
+      result->withLHS = !neumann;
-      eKm = &myMats->A[0];
+      result->rhsOnly = neumann;
-      eKg = &myMats->A[0];
+      result->resize(neumann?1:3,1);
-      eM  = &myMats->A[1];
-      Hh  = &myMats->A[2];
-      break;
-
-    case SIM::RHS_ONLY:
-      if (myMats->A.size() < 2)
-	myMats->resize(2,1);
-      else
-	myMats->b.resize(1);
-      eKm = &myMats->A[0];
-      eKg = &myMats->A[0];
-      Hh  = &myMats->A.back();
-      myMats->rhsOnly = true;
       break;
 
     case SIM::RECOVERY:
-      if (myMats->A.empty())
+      result->rhsOnly = true;
-	myMats->A.resize(1);
+      result->resize(1,0);
-      if (mySols.empty())
+      break;
-	mySols.resize(1);
-      Hh = &myMats->A.back();
-      eV = &mySols[0];
-      myMats->rhsOnly = true;
-      return;
 
     default:
-      this->Elasticity::setMode(mode);
+      return result;
-      return;
     }
 
-  // We always need the force+displacement vectors in nonlinear simulations
+  result->Hh = &result->A.back();
-  iS = &myMats->b[0];
+
-  eS = &myMats->b[0];
+  for (size_t i = 1; i < result->A.size(); i++)
-  mySols.resize(2);
+    result->A[i-1].resize(nsd*nen,nsd*nen);
-  eV = &mySols[0];
+
-  tracVal.clear();
+  if (result->b.size())
+    result->b.front().resize(nsd*nen);
+
+  return result;
 }
 
 
 bool NonlinearElasticityULMX::initElement (const std::vector<int>& MNPC,
-					   const Vec3& Xc, size_t nPt)
+					   const Vec3& Xc, size_t nPt,
+                                           LocalIntegral& elmInt)
 {
 #if INT_DEBUG > 0
   std::cout <<"\n\n *** Entering NonlinearElasticityULMX::initElement";
   std::cout <<", Xc = "<< Xc << std::endl;
 #endif
 
-  iP = 0;
+  return static_cast<MxMats&>(elmInt).init(Xc,nPt,p,nsd) &&
-  X0 = Xc;
+         this->IntegrandBase::initElement(MNPC,elmInt);
-  myData.resize(nPt);
-  size_t nPm = utl::Pascal(p,nsd);
-  if (Hh) Hh->resize(nPm,nPm,true);
-
-  // The other element matrices are initialized by the parent class method
-  return this->NonlinearElasticityUL::initElement(MNPC);
 }
 
 
-bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt,
+bool NonlinearElasticityULMX::evalInt (LocalIntegral& elmInt,
 				       const FiniteElement& fe,
 				       const TimeDomain&, const Vec3& X) const
 {
-  if (mySols.size() < 2) return false;
+  MxElmData* mx = 0;
+  MxNorm* mxNrm = 0;
+  MxMats* mxMat = dynamic_cast<MxMats*>(&elmInt);
+  if (mxMat)
+    mx = mxMat;
+  else if ((mxNrm = dynamic_cast<MxNorm*>(&elmInt)))
+    mx = mxNrm;
+  else
+    return false;
+
+  ItgPtData& ptData = mx->myData[mx->iP++];
+  const Vectors& eV = mxNrm ? mxNrm->myNorm->vec : elmInt.vec;
+  if (eV.size() < 2) return false;
 
-  ItgPtData& ptData = myData[iP++];
-  Vectors& eVs = const_cast<Vectors&>(mySols);
 
 #if INT_DEBUG > 1
   std::cout <<"NonlinearElasticityULMX::dNdX ="<< fe.dNdX;
 #endif
 
   // Evaluate the deformation gradient, Fp, at previous configuration
-  const_cast<NonlinearElasticityULMX*>(this)->eV = &eVs[1];
+  if (!this->formDefGradient(eV[1],fe.N,fe.dNdX,X.x,ptData.Fp))
-  if (!this->formDefGradient(fe.N,fe.dNdX,X.x,ptData.Fp))
     return false;
 
   // Evaluate the deformation gradient, F, at current configuration
-  const_cast<NonlinearElasticityULMX*>(this)->eV = &eVs[0];
+  if (!this->formDefGradient(eV[0],fe.N,fe.dNdX,X.x,ptData.F))
-  if (!this->formDefGradient(fe.N,fe.dNdX,X.x,ptData.F))
     return false;
 
   if (nDF == 2) // In 2D we always assume plane strain so set F(3,3)=1
@@ -184,13 +259,13 @@ bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt,
   // Push-forward the basis function gradients to current configuration
   ptData.dNdx.multiply(fe.dNdX,Fi); // dNdx = dNdX * F^-1
   if (axiSymmetry)
-    ptData.Nr = fe.N * (1.0/(X.x + eV->dot(fe.N,0,nsd)));
+    ptData.Nr = fe.N * (1.0/(X.x + eV.front().dot(fe.N,0,nsd)));
 
   ptData.X.assign(X);
   ptData.detJW = axiSymmetry ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW;
 
   // Evaluate the pressure modes (generalized coordinates)
-  Vec3 Xg = X-X0;
+  Vec3 Xg = X - mx->X0;
   if (nsd == 3)
     utl::Pascal(p,Xg.x,Xg.y,Xg.z,ptData.Phi);
   else if (nsd == 2)
@@ -198,33 +273,36 @@ bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt,
   else
     return false;
 
-  if (Hh)
+  if (mx->Hh)
     // Integrate the Hh-matrix
-    Hh->outer_product(ptData.Phi,ptData.Phi*ptData.detJW,true);
+    mx->Hh->outer_product(ptData.Phi,ptData.Phi*ptData.detJW,true);
 
   if (eM)
     // Integrate the mass matrix
-    this->formMassMatrix(*eM,fe.N,X,J*ptData.detJW);
+    this->formMassMatrix(mxMat->A[eM-1],fe.N,X,J*ptData.detJW);
 
   if (eS)
     // Integrate the load vector due to gravitation and other body forces
-    this->formBodyForce(*eS,fe.N,X,J*ptData.detJW);
+    this->formBodyForce(mxMat->b[eS-1],fe.N,X,J*ptData.detJW);
 
-  return this->getIntegralResult(elmInt);
+  return true;
 }
 
 
-bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
+bool NonlinearElasticityULMX::finalizeElement (LocalIntegral& elmInt,
 					       const TimeDomain& prm)
 {
   if (!iS && !eKm && !eKg) return true;
-  if (!Hh) return false;
 
+  MxMats& mx = static_cast<MxMats&>(elmInt);
+  if (!mx.Hh) return false;
+
+  size_t iP;
 #if INT_DEBUG > 0
   std::cout <<"\n\n *** Entering NonlinearElasticityULMX::finalizeElement\n";
-  for (iP = 1; iP <= myData.size(); iP++)
+  for (iP = 1; iP <= mx.myData.size(); iP++)
   {
-    const ItgPtData& pt = myData[iP-1];
+    const ItgPtData& pt = mx.myData[iP-1];
     std::cout <<"\n     X"   << iP <<" = "<< pt.X;
     std::cout <<"\n     detJ"<< iP <<"W = "<< pt.detJW;
     std::cout <<"\n     F"   << iP <<" =\n"<< pt.F;
@@ -237,14 +315,14 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
   // 1. Eliminate the internal pressure DOFs by static condensation.
 
   size_t i, j;
-  size_t nPM = Hh->rows();
+  size_t nPM = mx.Hh->rows();
-  size_t nEN = myData.front().dNdx.rows();
+  size_t nEN = mx.myData.front().dNdx.rows();
-  size_t nGP = myData.size();
+  size_t nGP = mx.myData.size();
 
   // Invert the H-matrix
-  if (!utl::invert(*Hh)) return false;
+  if (!utl::invert(*mx.Hh)) return false;
 
-  const Matrix& Hi = *Hh;
+  const Matrix& Hi = *mx.Hh;
 
   Matrix Theta(2,nGP), dNdxBar[nGP];
 
@@ -258,7 +336,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
 
     for (iP = 0; iP < nGP; iP++)
     {
-      const ItgPtData& pt = myData[iP];
+      const ItgPtData& pt = mx.myData[iP];
       double dVol = pt.F.det()*pt.detJW;
       h1 += dVol;
       h2 += pt.Fp.det()*pt.detJW;
@@ -292,7 +370,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
     // Integrate the Ji- and G-matrices
     for (iP = 0; iP < nGP; iP++)
     {
-      const ItgPtData& pt = myData[iP];
+      const ItgPtData& pt = mx.myData[iP];
       for (i = 0; i < nPM; i++)
       {
 	double h1 = pt.Phi[i] * pt.F.det() * pt.detJW;
@@ -332,7 +410,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
     for (iP = 0; iP < nGP; iP++)
     {
       i = iP+1;
-      ItgPtData& pt = myData[iP];
+      ItgPtData& pt = mx.myData[iP];
 
       // Calculate Theta = Hj*Phi
       Theta(1,i) = Hj1.dot(pt.Phi);
@@ -349,7 +427,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
   Vector detF(nGP);
   for (iP = 1; iP <= nGP; iP++)
   {
-    ItgPtData& pt = myData[iP-1];
+    ItgPtData& pt = mx.myData[iP-1];
     detF(iP) = pt.F.det();
     pt.F  *= pow(fabs(Theta(1,iP)/detF(iP)),1.0/3.0);
     pt.Fp *= pow(fabs(Theta(2,iP)/pt.Fp.det()),1.0/3.0);
@@ -364,7 +442,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
 
   // 2. Evaluate the constitutive relation and calculate mixed pressure state.
 
-  Matrix D[nGP];
+  Matrix Cmat, D[nGP];
   Vector Sigm(nPM), Hsig;
   std::vector<SymmTensor> Sig;
   Sig.reserve(nGP);
@@ -375,7 +453,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
 #if INT_DEBUG > 0
     std::cout <<"\n   iGP =     "<< iP+1 << std::endl;
 #endif
-    const ItgPtData& pt = myData[iP];
+    const ItgPtData& pt = mx.myData[iP];
 
     // Evaluate the constitutive relation
     double U = 0.0;
@@ -405,7 +483,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
 
   // Divide pressure by reference volume; Hsig = Hi*Sigm
   if (lHaveStress)
-    if (!Hh->multiply(Sigm,Hsig)) return false;
+    if (!mx.Hh->multiply(Sigm,Hsig)) return false;
 
   // 3. Evaluate tangent matrices.
 
@@ -413,7 +491,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
   SymmTensor Sigma(nsd,axiSymmetry);
   for (iP = 0; iP < nGP; iP++)
   {
-    const ItgPtData& pt = myData[iP];
+    const ItgPtData& pt = mx.myData[iP];
     double dFmix = Theta(1,iP+1);
     if (lHaveStress)
     {
@@ -425,7 +503,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
 #endif
       if (eKg)
 	// Integrate the geometric stiffness matrix
-	this->formKG(*eKg,pt.Nr,pt.dNdx,1.0,Sigma,dFmix*pt.detJW);
+	this->formKG(mx.A[eKg-1],pt.Nr,pt.dNdx,1.0,Sigma,dFmix*pt.detJW);
     }
 
     if (eKm)
@@ -437,31 +515,33 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
       D[iP] *= dFmix*pt.detJW;
       if (nsd == 2)
 	acckmx2d_(axiSymmetry,nEN,pt.Nr.ptr(),pt.dNdx.ptr(),dNdxBar[iP].ptr(),
-		  D[iP].ptr(),eKm->ptr());
+		  D[iP].ptr(),mx.A[eKm-1].ptr());
       else
-	acckmx3d_(nEN,pt.dNdx.ptr(),dNdxBar[iP].ptr(),D[iP].ptr(),eKm->ptr());
+	acckmx3d_(nEN,pt.dNdx.ptr(),dNdxBar[iP].ptr(),
+		  D[iP].ptr(),mx.A[eKm-1].ptr());
 #endif
     }
 
     if (iS && lHaveStress)
     {
       // Compute the small-deformation strain-displacement matrix B from dNdx
+      Matrix Bmat;
       if (axiSymmetry)
-	this->formBmatrix(pt.Nr,pt.dNdx,1.0);
+	this->formBmatrix(Bmat,pt.Nr,pt.dNdx,1.0);
       else
-	this->formBmatrix(pt.dNdx);
+	this->formBmatrix(Bmat,pt.dNdx);
 
       // Integrate the internal forces
       Sigma *= -dFmix*pt.detJW;
-      if (!Bmat.multiply(Sigma,*iS,true,true)) // ES -= B^T*sigma
+      if (!Bmat.multiply(Sigma,mx.b[iS-1],true,true)) // ES -= B^T*sigma
 	return false;
     }
   }
 
 #if INT_DEBUG > 0
   std::cout <<"\n *** Leaving NonlinearElasticityULMX::finalizeElement";
-  if (eKm) std::cout <<", eKm ="<< *eKm;
+  if (eKm) std::cout <<", eKm ="<< mx.A[eKm-1];
-  if (iS) std::cout <<" iS ="<< *iS;
+  if (iS) std::cout <<" iS ="<< mx.b[iS-1];
   std::cout << std::endl;
 #endif
   return true;
@@ -474,59 +554,86 @@ NormBase* NonlinearElasticityULMX::getNormIntegrand (AnaSol*) const
 }
 
 
-bool ElasticityNormULMX::initElement (const std::vector<int>& MNPC,
+LocalIntegral* ElasticityNormULMX::getLocalIntegral (size_t nen, size_t iEl,
-				      const Vec3& Xc, size_t nPt)
+						     bool neumann) const
 {
-  NonlinearElasticityULMX& p = static_cast<NonlinearElasticityULMX&>(myProblem);
+  LocalIntegral* result = this->NormBase::getLocalIntegral(nen,iEl,neumann);
+  if (result) return new MxNorm(static_cast<ElmNorm&>(*result));
 
-  p.iP = 0;
+  // Element norms are not requested, so allocate an inter object instead that
-  p.X0 = Xc;
+  // will delete itself when invoking the destruct method.
-  p.myData.resize(nPt);
+  return new MxNorm(this->getNoFields());
-  size_t nPm = utl::Pascal(p.p,p.nsd);
-  if (p.Hh) p.Hh->resize(nPm,nPm,true);
-
-  // The other element matrices are initialized by the parent class method
-  return p.NonlinearElasticityUL::initElement(MNPC);
 }
 
 
-bool ElasticityNormULMX::evalInt (LocalIntegral*& elmInt,
+bool ElasticityNormULMX::initElement (const std::vector<int>& MNPC,
+				      const Vec3& Xc, size_t nPt,
+				      LocalIntegral& elmInt)
+{
+  NonlinearElasticityULMX& p = static_cast<NonlinearElasticityULMX&>(myProblem);
+  MxNorm& mx = static_cast<MxNorm&>(elmInt);
+
+  return mx.init(Xc,nPt,p.p,p.nsd) &&
+         this->NormBase::initElement(MNPC,*mx.myNorm);
+}
+
+
+bool ElasticityNormULMX::initElementBou (const std::vector<int>& MNPC,
+					 LocalIntegral& elmInt)
+{
+  return myProblem.initElementBou(MNPC,*static_cast<MxNorm&>(elmInt).myNorm);
+}
+
+
+bool ElasticityNormULMX::evalInt (LocalIntegral& elmInt,
 				  const FiniteElement& fe,
 				  const TimeDomain& td, const Vec3& X) const
 {
-  return static_cast<NonlinearElasticityULMX&>(myProblem).evalInt(elmInt,
+  NonlinearElasticityULMX& p = static_cast<NonlinearElasticityULMX&>(myProblem);
-								  fe,td,X);
+
+  return p.evalInt(elmInt,fe,td,X);
 }
 
 
-bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
+bool ElasticityNormULMX::evalBou (LocalIntegral& elmInt,
+                                  const FiniteElement& fe,
+                                  const Vec3& X, const Vec3& normal) const
+{
+  return this->ElasticityNormUL::evalBou(*static_cast<MxNorm&>(elmInt).myNorm,
+					 fe,X,normal);
+}
+
+
+bool ElasticityNormULMX::finalizeElement (LocalIntegral& elmInt,
 					  const TimeDomain& prm)
 {
   NonlinearElasticityULMX& p = static_cast<NonlinearElasticityULMX&>(myProblem);
-  if (!p.Hh) return false;
+  MxNorm& mx = static_cast<MxNorm&>(elmInt);
+
+  if (!mx.Hh) return false;
 
   size_t iP;
 #if INT_DEBUG > 0
   std::cout <<"\n\n *** Entering ElasticityNormULMX::finalizeElement\n";
-  for (iP = 1; iP <= p.myData.size(); iP++)
+  for (iP = 1; iP <= mx.myData.size(); iP++)
   {
-    const NonlinearElasticityULMX::ItgPtData& pt = p.myData[iP-1];
+    const NonlinearElasticityULMX::ItgPtData& pt = mx.myData[iP-1];
     std::cout <<"\n     X"   << iP <<" = "<< pt.X;
     std::cout <<"\n     detJ"<< iP <<"W = "<< pt.detJW;
     std::cout <<"\n     F"   << iP <<" =\n"<< pt.F;
     std::cout <<"     dNdx"<< iP <<" ="<< pt.dNdx;
     std::cout <<"     Phi" << iP <<" ="<< pt.Phi;
   }
-  std::cout <<"\n     H ="<< *(p.Hh) << std::endl;
+  std::cout <<"\n     H ="<< *(mx.Hh) << std::endl;
 #endif
 
   // 1. Eliminate the internal pressure DOFs by static condensation.
 
-  size_t nPM = p.Hh->rows();
+  size_t nPM = mx.Hh->rows();
-  size_t nGP = p.myData.size();
+  size_t nGP = mx.myData.size();
 
   // Invert the H-matrix
-  Matrix& Hi = *p.Hh;
+  Matrix& Hi = *mx.Hh;
   if (!utl::invert(Hi)) return false;
 
   Vector Theta(nGP);
@@ -537,7 +644,7 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
 
     double h1 = 0.0;
     for (iP = 0; iP < nGP; iP++)
-      h1 += p.myData[iP].F.det() * p.myData[iP].detJW;
+      h1 += mx.myData[iP].F.det() * mx.myData[iP].detJW;
 
     // All gauss point values are identical
     Theta.front() = h1 * Hi(1,1);
@@ -551,7 +658,7 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
     // Integrate the Ji-matrix
     for (iP = 0; iP < nGP; iP++)
     {
-      const NonlinearElasticityULMX::ItgPtData& pt = p.myData[iP];
+      const ItgPtData& pt = mx.myData[iP];
       for (size_t i = 0; i < nPM; i++)
 	Ji[i] += pt.Phi[i] * pt.F.det() * pt.detJW;
     }
@@ -566,7 +673,7 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
 
     // Calculate Theta = Hj*Phi
     for (iP = 0; iP < nGP; iP++)
-      Theta[iP] = Hj.dot(p.myData[iP].Phi);
+      Theta[iP] = Hj.dot(mx.myData[iP].Phi);
   }
 
 #if INT_DEBUG > 1
@@ -575,15 +682,14 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
 
   // 2. Evaluate the constitutive relation and integrate energy.
 
-  ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt,6);
+  Matrix Cmat;
-
   SymmTensor Sig(3), Eps(3);
   for (iP = 0; iP < nGP; iP++)
   {
 #if INT_DEBUG > 0
     std::cout <<"\n   iGP =     "<< iP+1 << std::endl;
 #endif
-    NonlinearElasticityULMX::ItgPtData& pt = p.myData[iP];
+    ItgPtData& pt = mx.myData[iP];
     Eps.rightCauchyGreen(pt.F); // Green Lagrange strain tensor
     Eps -= 1.0;
     Eps *= 0.5;
@@ -595,11 +701,11 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
     // Compute the strain energy density, U(Eps) = Int_Eps (S:E) dE
     // and the Cauchy stress tensor, Sig
     double U = 0.0;
-    if (!p.material->evaluate(p.Cmat,Sig,U,pt.X,Fbar,Eps,3,&prm,&pt.F))
+    if (!p.material->evaluate(Cmat,Sig,U,pt.X,Fbar,Eps,3,&prm,&pt.F))
       return false;
 
     // Integrate the norms
-    if (!ElasticityNormUL::evalInt(pnorm,Sig,U,Theta[iP],pt.detJW))
+    if (!ElasticityNormUL::evalInt(*mx.myNorm,Sig,U,Theta[iP],pt.detJW))
       return false;
   }
 

Apps/FiniteDefElasticity/NonlinearElasticityULMX.h

@@ -15,7 +15,6 @@
 #define _NONLINEAR_ELASTICITY_UL_MX_H
 
 #include "NonlinearElasticityUL.h"
-#include "Tensor.h"
 
 
 /*!
@@ -26,20 +25,6 @@
 
 class NonlinearElasticityULMX : public NonlinearElasticityUL
 {
-  //! \brief A struct with integration point data needed by \a finalizeElement.
-  struct ItgPtData
-  {
-    Tensor F;     //!< Deformation gradient at current step/iteration
-    Tensor Fp;    //!< Deformation gradient at previous (converged) step
-    Vector Nr;    //!< Basis function values (for axisymmetric problems)
-    Matrix dNdx;  //!< Basis function gradients at current configuration
-    Vector Phi;   //!< Internal modes for the pressure/volumetric-change fields
-    Vec4   X;     //!< Cartesian coordinates of current integration point
-    double detJW; //!< Jacobian determinant times integration point weight
-    //! \brief Default constructor.
-    ItgPtData() : F(3), Fp(3) { detJW = 0.0; }
-  };
-
 public:
   //! \brief The default constructor invokes the parent class constructor.
   //! \param[in] n Number of spatial dimensions
@@ -53,16 +38,20 @@ public:
   //! \brief Prints out problem definition to the given output stream.
   virtual void print(std::ostream& os) const;
 
-  //! \brief Defines the solution mode before the element assembly is started.
+  //! \brief Returns a local integral container for the given element.
-  //! \param[in] mode The solution mode to use
+  //! \param[in] nen Number of nodes on element
-  virtual void setMode(SIM::SolutionMode mode);
+  //! \param[in] neumann Whether or not we are assembling Neumann BC's
+  virtual LocalIntegral* getLocalIntegral(size_t nen, size_t,
+					  bool neumann) const;
 
   //! \brief Initializes current element for numerical integration.
   //! \param[in] MNPC Matrix of nodal point correspondance for current element
   //! \param[in] Xc Cartesian coordinates of the element center
   //! \param[in] nPt Number of integration points in this element
+  //! \param elmInt The local integral object for current element
   virtual bool initElement(const std::vector<int>& MNPC,
-			   const Vec3& Xc, size_t nPt);
+			   const Vec3& Xc, size_t nPt,
+			   LocalIntegral& elmInt);
 
   //! \brief Evaluates the integrand at an interior point.
   //! \param elmInt The local integral object to receive the contributions
@@ -74,17 +63,17 @@ public:
   //! depending on the basis function values in the internal member \a myData.
   //! The actual numerical integration of the tangent stiffness and internal
   //! forces is then performed by the \a finalizeElement method.
-  virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe,
+  virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
 		       const TimeDomain& prm, const Vec3& X) const;
 
   //! \brief Finalizes the element matrices after the numerical integration.
   //! \details This method is used to implement the multiple integration point
   //! loops within an element, required by the mixed formulation with internal
   //! modes. All data needed during the integration has been stored internally
-  //! in the \a myData member by the \a evalInt method.
+  //! in the \a elmInt object by the \a evalInt method.
   //! \param elmInt The local integral object to receive the contributions
   //! \param[in] prm Nonlinear solution algorithm parameters
-  virtual bool finalizeElement(LocalIntegral*& elmInt, const TimeDomain& prm);
+  virtual bool finalizeElement(LocalIntegral& elmInt, const TimeDomain& prm);
 
   //! \brief Returns a pointer to an Integrand for solution norm evaluation.
   //! \note The Integrand object is allocated dynamically and has to be deleted
@@ -96,13 +85,7 @@ public:
   virtual int getIntegrandType() const { return 4; }
 
 private:
-  int  p;  //!< Polynomial order of the internal pressure field
+  int p; //!< Polynomial order of the internal pressure field
-  Vec3 X0; //!< Cartesian coordinates of the element center
-
-  Matrix* Hh; //!< Pointer to element Hh-matrix associated with pressure modes
-
-  mutable size_t                 iP;     //!< Local integration point counter
-  mutable std::vector<ItgPtData> myData; //!< Local integration point data
 
   friend class ElasticityNormULMX;
 };
@@ -121,12 +104,27 @@ public:
   //! \brief Empty destructor.
   virtual ~ElasticityNormULMX() {}
 
+  //! \brief Returns a local integral contribution object for the given element.
+  //! \param[in] nen Number of nodes on element
+  //! \param[in] iEl The element number
+  //! \param[in] neumann Whether or not we are assembling Neumann BC's
+  virtual LocalIntegral* getLocalIntegral(size_t nen, size_t iEl,
+					  bool neumann) const;
+
   //! \brief Initializes current element for numerical integration.
   //! \param[in] MNPC Matrix of nodal point correspondance for current element
   //! \param[in] Xc Cartesian coordinates of the element center
   //! \param[in] nPt Number of integration points in this element
+  //! \param elmInt The local integral object for current element
   virtual bool initElement(const std::vector<int>& MNPC,
-			   const Vec3& Xc, size_t nPt);
+			   const Vec3& Xc, size_t nPt,
+			   LocalIntegral& elmInt);
+
+  //! \brief Initializes current element for boundary integration.
+  //! \param[in] MNPC Matrix of nodal point correspondance for current element
+  //! \param elmInt The local integral object for current element
+  virtual bool initElementBou(const std::vector<int>& MNPC,
+			      LocalIntegral& elmInt);
 
   //! \brief Evaluates the integrand at an interior point.
   //! \param elmInt The local integral object to receive the contributions
@@ -138,16 +136,24 @@ public:
   //! corresponding NonlinearElasticityULMX object, for which this class is
   //! declared as friend such that it can access the data members. The actual
   //! norm integration us then performed by the \a finalizeElement method.
-  virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe,
+  virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
 		       const TimeDomain& prm, const Vec3& X) const;
 
+  //! \brief Evaluates the integrand at a boundary point.
+  //! \param elmInt The local integral object to receive the contributions
+  //! \param[in] fe Finite element data of current integration point
+  //! \param[in] X Cartesian coordinates of current integration point
+  //! \param[in] normal Boundary normal vector at current integration point
+  virtual bool evalBou(LocalIntegral& elmInt, const FiniteElement& fe,
+		       const Vec3& X, const Vec3& normal) const;
+
   //! \brief Finalizes the element norms after the numerical integration.
   //! \details This method is used to implement the multiple integration point
   //! loops within an element required by the mixed formulation. It performs
   //! a subset of the tasks done by NonlinearElasticityULMX::finalizeElement.
   //! \param elmInt The local integral object to receive the contributions
   //! \param[in] prm Nonlinear solution algorithm parameters
-  virtual bool finalizeElement(LocalIntegral*& elmInt, const TimeDomain& prm);
+  virtual bool finalizeElement(LocalIntegral& elmInt, const TimeDomain& prm);
 
   //! \brief Returns which integrand to be used.
   virtual int getIntegrandType() const { return 4; }

Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C

@@ -15,7 +15,7 @@
 #include "MaterialBase.h"
 #include "FiniteElement.h"
 #include "TimeDomain.h"
-#include "ElmMats.h"
+#include "ElmNorm.h"
 #include "Utilities.h"
 #include "Vec3Oper.h"
 
@@ -35,6 +35,11 @@ extern "C" {
   void acckm3d_(const int& nEN, const double* dNdx,
 		const double* D, double* eKt);
 }
+#ifdef USE_OPENMP
+#if INT_DEBUG > 0
+#undef INT_DEBUG
+#endif
+#endif
 #ifndef INT_DEBUG
 #define INT_DEBUG 0
 #endif
@@ -195,18 +200,10 @@ const Vector& NonlinearElasticityULMixed::MixedElmMats::getRHSVector () const
 
 NonlinearElasticityULMixed::NonlinearElasticityULMixed (unsigned short int n,
 							bool axS)
-  : NonlinearElasticityUL(n,axS), Fbar(3), Dmat(7,7)
+  : NonlinearElasticityUL(n,axS)
 {
-  if (myMats) delete myMats;
+  eKm = eKg = Kuu+1;
-  myMats = new MixedElmMats();
+  iS  = eS  = Ru+1;
-  mySols.resize(NSOL);
-
-  eKm = &myMats->A[Kuu];
-  eKg = &myMats->A[Kuu];
-
-  iS = &myMats->b[Ru];
-  eS = &myMats->b[Ru];
-  eV = &mySols[U];
 }
 
 
@@ -221,18 +218,14 @@ void NonlinearElasticityULMixed::print (std::ostream& os) const
 
 void NonlinearElasticityULMixed::setMode (SIM::SolutionMode mode)
 {
+  m_mode = mode;
   switch (mode)
     {
     case SIM::INIT:
     case SIM::STATIC:
-      tracVal.clear();
-      myMats->rhsOnly = false;
-      break;
-
     case SIM::RHS_ONLY:
       tracVal.clear();
     case SIM::RECOVERY:
-      myMats->rhsOnly = true;
       break;
 
     default:
@@ -242,100 +235,111 @@ void NonlinearElasticityULMixed::setMode (SIM::SolutionMode mode)
 }
 
 
+LocalIntegral* NonlinearElasticityULMixed::getLocalIntegral (size_t nen1,
+							     size_t nen2,
+							     size_t,
+							     bool neumann) const
+{
+  const size_t nedof1 = nsd*nen1;
+  const size_t nedof  = nedof1 + 2*nen2;
+
+  ElmMats* result = new MixedElmMats;
+
+  result->rhsOnly = neumann || m_mode == SIM::RECOVERY;
+  result->b[Rres].resize(nedof);
+  result->b[Ru].resize(nedof1);
+  if (!neumann)
+  {
+    result->A[Kuu].resize(nedof1,nedof1);
+    result->A[Kut].resize(nedof1,nen2);
+    result->A[Kup].resize(nedof1,nen2);
+    result->A[Ktt].resize(nen2,nen2);
+    result->A[Ktp].resize(nen2,nen2);
+    result->A[Ktan].resize(nedof,nedof);
+    result->b[Rt].resize(nen2);
+    result->b[Rp].resize(nen2);
+  }
+  if (m_mode == SIM::STATIC)
+    result->withLHS = !neumann;
+
+  return result;
+}
+
+
+#ifndef USE_OPENMP
 static int iP = 0; //!< Local integration point counter for debug output
 std::vector<int> mixedDbgEl; //!< List of elements for additional output
+#endif
 
 
 bool NonlinearElasticityULMixed::initElement (const std::vector<int>& MNPC1,
 					      const std::vector<int>& MNPC2,
-					      size_t n1)
+					      size_t n1, LocalIntegral& elmInt)
 {
+#ifndef USE_OPENMP
   iP = 0;
+#endif
+  if (primsol.front().empty()) return true;
 
-  const size_t nen1 = MNPC1.size();
+  // Extract the element level solution vectors
-  const size_t nen2 = MNPC2.size();
+  elmInt.vec.resize(NSOL);
-  const size_t nedof = nsd*nen1 + 2*nen2;
+  int ierr = utl::gather(MNPC1,nsd,primsol.front(),elmInt.vec[U])
+    +        utl::gather(MNPC2,0,2,primsol.front(),elmInt.vec[T],nsd*n1,n1)
+    +        utl::gather(MNPC2,1,2,primsol.front(),elmInt.vec[P],nsd*n1,n1);
+  if (ierr == 0) return true;
 
-  myMats->A[Kut].resize(nsd*nen1,nen2,true);
+  std::cerr <<" *** NonlinearElasticityULMixed::initElement: Detected "
-  myMats->A[Kup].resize(nsd*nen1,nen2,true);
+	    << ierr/3 <<" node numbers out of range."<< std::endl;
-  myMats->A[Ktt].resize(nen2,nen2,true);
+  return false;
-  myMats->A[Ktp].resize(nen2,nen2,true);
-  myMats->A[Ktan].resize(nedof,nedof,true);
-
-  myMats->b[Rt].resize(nen2,true);
-  myMats->b[Rp].resize(nen2,true);
-  myMats->b[Rres].resize(nedof,true);
-
-  // Extract the element level volumetric change and pressure vectors
-  int ierr = 0;
-  if (!primsol.front().empty())
-    if ((ierr = utl::gather(MNPC2,0,2,primsol.front(),mySols[T],nsd*n1,n1) +
-		utl::gather(MNPC2,1,2,primsol.front(),mySols[P],nsd*n1,n1)))
-      std::cerr <<" *** NonlinearElasticityULMixed::initElement: Detected "
-		<< ierr/2 <<" node numbers out of range."<< std::endl;
-
-  // The other element matrices are initialized by the parent class method
-  return this->NonlinearElasticityUL::initElement(MNPC1) && ierr == 0;
 }
 
 
 bool NonlinearElasticityULMixed::initElementBou (const std::vector<int>& MNPC1,
-						 const std::vector<int>& MNPC2,
+						 const std::vector<int>&,
-						 size_t)
+						 size_t, LocalIntegral& elmInt)
 {
-  const size_t nen1 = MNPC1.size();
+  return this->IntegrandBase::initElementBou(MNPC1,elmInt);
-  const size_t nen2 = MNPC2.size();
-  const size_t nedof = nsd*nen1 + 2*nen2;
-
-  myMats->b[Ru].resize(nsd*nen1,true);
-  myMats->b[Rres].resize(nedof,true);
-
-  // Extract the element level displacement vector
-  int ierr = 0;
-  if (!primsol.front().empty())
-    if ((ierr = utl::gather(MNPC1,nsd,primsol.front(),mySols[U])))
-      std::cerr <<" *** NonlinearElasticityULMixed::initElementBou: Detected "
-                << ierr <<" node numbers out of range."<< std::endl;
-
-  myMats->withLHS = false;
-  return ierr == 0;
 }
 
 
-bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
+bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral& elmInt,
 					    const MxFiniteElement& fe,
 					    const TimeDomain& prm,
 					    const Vec3& X) const
 {
+  ElmMats& elMat = static_cast<ElmMats&>(elmInt);
+
 #if INT_DEBUG > 0
   std::cout <<"\n\n *** Entering NonlinearElasticityULMixed::evalIntMx: iP = "
 	    << ++iP <<", (X,t) = "<< X << std::endl;
 #endif
 
   // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
+  Matrix B;
   Tensor F(nDF);
   SymmTensor E(nsd,axiSymmetry);
-  if (!this->kinematics(fe.N1,fe.dN1dX,X.x,F,E))
+  if (!this->kinematics(elmInt.vec[U],fe.N1,fe.dN1dX,X.x,B,F,E))
     return false;
 
   bool lHaveStrains = !E.isZero(1.0e-16);
 
   // Evaluate the volumetric change and pressure fields
-  double Theta = fe.N2.dot(mySols[T]) + 1.0;
+  double Theta = fe.N2.dot(elmInt.vec[T]) + 1.0;
-  double Press = fe.N2.dot(mySols[P]);
+  double Press = fe.N2.dot(elmInt.vec[P]);
 #if INT_DEBUG > 0
-  std::cout <<"NonlinearElasticityULMixed::b_theta ="<< mySols[T];
+  std::cout <<"NonlinearElasticityULMixed::b_theta ="<< elmInt.vec[T];
-  std::cout <<"NonlinearElasticityULMixed::b_p ="<< mySols[P];
+  std::cout <<"NonlinearElasticityULMixed::b_p ="<< elmInt.vec[P];
   std::cout <<"NonlinearElasticityULMixed::Theta = "<< Theta
 	    <<", Press = "<< Press << std::endl;
 #endif
 
   // Axi-symmetric integration point volume; 2*pi*r*|J|*w
   double detJW = axiSymmetry ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW;
-  double r     = axiSymmetry ? X.x + eV->dot(fe.N1,0,nsd) : 0.0;
+  double r     = axiSymmetry ? X.x + elmInt.vec[U].dot(fe.N1,0,nsd) : 0.0;
 
   // Compute the mixed model deformation gradient, F_bar
-  Fbar = F; // notice that F_bar always has dimension 3
+  Tensor Fbar(3); // notice that F_bar always has dimension 3
+  Fbar = F;
   double r1 = pow(fabs(Theta/F.det()),1.0/3.0);
 
   Fbar *= r1;
@@ -356,6 +360,7 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
   if (J == 0.0) return false;
 
   // Push-forward the basis function gradients to current configuration
+  Matrix dNdx;
   dNdx.multiply(fe.dN1dX,Fi); // dNdx = dNdX * F^-1
 
   // Compute the mixed integration point volume
@@ -373,18 +378,20 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
 
   if (eM)
     // Integrate the mass matrix
-    this->formMassMatrix(*eM,fe.N1,X,J*detJW);
+    this->formMassMatrix(elMat.A[eM-1],fe.N1,X,J*detJW);
 
   if (eS)
     // Integrate the load vector due to gravitation and other body forces
-    this->formBodyForce(*eS,fe.N1,X,J*detJW);
+    this->formBodyForce(elMat.b[eS-1],fe.N1,X,J*detJW);
 
   // Evaluate the constitutive relation
+  Matrix Cmat;
   SymmTensor Sig(3), Sigma(nsd,axiSymmetry);
   double U, Bpres = 0.0, Mpres = 0.0;
   if (!material->evaluate(Cmat,Sig,U,X,Fbar,E,lHaveStrains,&prm))
     return false;
 
+  Matrix Dmat(7,7);
 #ifdef USE_FTNMAT
   // Project the constitutive matrix for the mixed model
   pcst3d_(Cmat.ptr(),Dmat.ptr(),INT_DEBUG,6);
@@ -397,6 +404,7 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
 
   if (lHaveStrains)
   {
+#ifndef USE_OPENMP
     if (prm.it == 0 &&
 	find(mixedDbgEl.begin(),mixedDbgEl.end(),fe.iel) != mixedDbgEl.end())
     {
@@ -413,6 +421,7 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
 	std::cout <<" "<< sigma[i];
       std::cout <<" "<< Sig.vonMises() << std::endl;
     }
+#endif
 
     // Compute mixed strees
     Bpres = Sig.trace()/3.0;
@@ -420,7 +429,7 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
     Sigma = Sig + (Mpres-Bpres);
 
     // Integrate the geometric stiffness matrix
-    this->formKG(*eKg,fe.N1,dNdx,r,Sigma,dVol);
+    this->formKG(elMat.A[eKg-1],fe.N1,dNdx,r,Sigma,dVol);
   }
 
   // Radial shape function for axi-symmetric problems
@@ -434,9 +443,11 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
   // Integrate the material stiffness matrix
   Dmat *= dVol;
   if (nsd == 2)
-    acckm2d_(axiSymmetry,Nr.size(),Nr.ptr(),dNdx.ptr(),Dmat.ptr(),eKm->ptr());
+    acckm2d_(axiSymmetry,Nr.size(),Nr.ptr(),dNdx.ptr(),Dmat.ptr(),
+	     elMat.A[eKm-1].ptr());
   else
-    acckm3d_(fe.N1.size(),dNdx.ptr(),Dmat.ptr(),eKm->ptr());
+    acckm3d_(fe.N1.size(),dNdx.ptr(),Dmat.ptr(),
+	     elMat.A[eKm-1].ptr());
 #endif
 
   // Integrate the volumetric change and pressure tangent terms
@@ -450,55 +461,55 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
     {
       for (i = 1; i <= nsd; i++)
       {
-	myMats->A[Kut](nsd*(a-1)+i,b) += dNdx(a,i)*fe.N2(b) * Dmat(i,7);
+	elMat.A[Kut](nsd*(a-1)+i,b) += dNdx(a,i)*fe.N2(b) * Dmat(i,7);
 	if (nsd == 2)
-	  myMats->A[Kut](nsd*(a-1)+i,b) += dNdx(a,3-i)*fe.N2(b) * Dmat(4,7);
+	  elMat.A[Kut](nsd*(a-1)+i,b) += dNdx(a,3-i)*fe.N2(b) * Dmat(4,7);
 	else if (i < 3)
 	{
 	  j = i + 1;
 	  k = j%3 + 1;
-	  myMats->A[Kut](nsd*(a-1)+i,b) += dNdx(a,3-i)*fe.N2(b) * Dmat(4,7);
+	  elMat.A[Kut](nsd*(a-1)+i,b) += dNdx(a,3-i)*fe.N2(b) * Dmat(4,7);
-	  myMats->A[Kut](nsd*(a-1)+j,b) += dNdx(a,5-j)*fe.N2(b) * Dmat(5,7);
+	  elMat.A[Kut](nsd*(a-1)+j,b) += dNdx(a,5-j)*fe.N2(b) * Dmat(5,7);
-	  myMats->A[Kut](nsd*(a-1)+k,b) += dNdx(a,4-k)*fe.N2(b) * Dmat(6,7);
+	  elMat.A[Kut](nsd*(a-1)+k,b) += dNdx(a,4-k)*fe.N2(b) * Dmat(6,7);
 	}
-	myMats->A[Kup](nsd*(a-1)+i,b) += dNdx(a,i)*fe.N2(b) * dVup;
+	elMat.A[Kup](nsd*(a-1)+i,b) += dNdx(a,i)*fe.N2(b) * dVup;
-      } 
+      }
       if (axiSymmetry)
       {
 	double Nr = fe.N1(a)/r;
-	myMats->A[Kut](nsd*(a-1)+1,b) += Nr*fe.N2(b) * Dmat(3,7);
+	elMat.A[Kut](nsd*(a-1)+1,b) += Nr*fe.N2(b) * Dmat(3,7);
-	myMats->A[Kup](nsd*(a-1)+1,b) += Nr*fe.N2(b) * dVup;
+	elMat.A[Kup](nsd*(a-1)+1,b) += Nr*fe.N2(b) * dVup;
       }
     }
 
-  myMats->A[Ktt].outer_product(fe.N2,fe.N2*Dmat(7,7),true); // += N2*N2^T*D77
+  elMat.A[Ktt].outer_product(fe.N2,fe.N2*Dmat(7,7),true); // += N2*N2^T*D77
-  myMats->A[Ktp].outer_product(fe.N2,fe.N2*(-detJW),true);  // -= N2*N2^T*|J|
+  elMat.A[Ktp].outer_product(fe.N2,fe.N2*(-detJW),true);  // -= N2*N2^T*|J|
 
   if (lHaveStrains)
   {
     // Compute the small-deformation strain-displacement matrix B from dNdx
     if (axiSymmetry)
-      this->formBmatrix(fe.N1,dNdx,r);
+      this->formBmatrix(B,fe.N1,dNdx,r);
     else
-      this->formBmatrix(dNdx);
+      this->formBmatrix(B,dNdx);
 
     // Integrate the internal forces
     Sigma *= -dVol;
-    if (!Bmat.multiply(Sigma,myMats->b[Ru],true,true))  // -= B^T*sigma
+    if (!B.multiply(Sigma,elMat.b[Ru],true,true))  // -= B^T*sigma
       return false;
 
     // Integrate the volumetric change and pressure forces
-    myMats->b[Rt].add(fe.N2,(Press - Bpres)*detJW); // += N2*(p-pBar)*|J|
+    elMat.b[Rt].add(fe.N2,(Press - Bpres)*detJW); // += N2*(p-pBar)*|J|
-    myMats->b[Rp].add(fe.N2,(Theta - J)*detJW);     // += N2*(Theta-J)*|J|
+    elMat.b[Rp].add(fe.N2,(Theta - J)*detJW);     // += N2*(Theta-J)*|J|
 
 #if INT_DEBUG > 4
     std::cout <<"NonlinearElasticityULMixed::Sigma*dVol =\n"<< Sigma;
-    std::cout <<"NonlinearElasticityULMixed::Bmat ="<< Bmat;
+    std::cout <<"NonlinearElasticityULMixed::Bmat ="<< B;
-    std::cout <<"NonlinearElasticityULMixed::Ru("<< iP <<") ="<< myMats->b[Ru];
+    std::cout <<"NonlinearElasticityULMixed::Ru("<< iP <<") ="<< elMat.b[Ru];
 #endif
   }
 
-  return this->getIntegralResult(elmInt);
+  return true;
 }
 
 
@@ -508,7 +519,7 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
   forwarded to the corresponding single-field method of the parent class.
 */
 
-bool NonlinearElasticityULMixed::evalBouMx (LocalIntegral*& elmInt,
+bool NonlinearElasticityULMixed::evalBouMx (LocalIntegral& elmInt,
 					    const MxFiniteElement& fe,
 					    const Vec3& X,
 					    const Vec3& normal) const
@@ -557,7 +568,7 @@ NormBase* NonlinearElasticityULMixed::getNormIntegrand (AnaSol*) const
 }
 
 
-bool ElasticityNormULMixed::evalIntMx (LocalIntegral*& elmInt,
+bool ElasticityNormULMixed::evalIntMx (LocalIntegral& elmInt,
 				       const MxFiniteElement& fe,
 				       const TimeDomain& prm,
 				       const Vec3& X) const
@@ -566,40 +577,44 @@ bool ElasticityNormULMixed::evalIntMx (LocalIntegral*& elmInt,
   ulp = static_cast<NonlinearElasticityULMixed*>(&myProblem);
 
   // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
+  Matrix B;
   Tensor F(ulp->nDF);
   SymmTensor E(ulp->nDF);
-  if (!ulp->kinematics(fe.N1,fe.dN1dX,X.x,F,E))
+  if (!ulp->kinematics(elmInt.vec[U],fe.N1,fe.dN1dX,X.x,B,F,E))
     return false;
 
   // Evaluate the volumetric change field
-  double Theta = fe.N2.dot(ulp->mySols[T]) + 1.0;
+  double Theta = fe.N2.dot(elmInt.vec[T]) + 1.0;
 
   // Compute the mixed model deformation gradient, F_bar
-  ulp->Fbar = F; // notice that F_bar always has dimension 3
+  Tensor Fbar(3); // notice that F_bar always has dimension 3
+  Fbar = F; // notice that F_bar always has dimension 3
   double r1 = pow(fabs(Theta/F.det()),1.0/3.0);
 
-  ulp->Fbar *= r1;
+  Fbar *= r1;
   if (F.dim() == 2) // In 2D we always assume plane strain so set F(3,3)=1
-    ulp->Fbar(3,3) = r1;
+    Fbar(3,3) = r1;
 
   // Compute the strain energy density, U(E) = Int_E (Sig:Eps) dEps
   // and the Cauchy stress tensor, Sig
+  Matrix Cmat;
   double U = 0.0;
   SymmTensor Sig(3);
-  if (!ulp->material->evaluate(ulp->Cmat,Sig,U,X,ulp->Fbar,E,3,&prm,&F))
+  if (!ulp->material->evaluate(Cmat,Sig,U,X,Fbar,E,3,&prm,&F))
     return false;
 
   // Axi-symmetric integration point volume; 2*pi*r*|J|*w
   double detJW = ulp->isAxiSymmetric() ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW;
 
   // Integrate the norms
-  return evalInt(getElmNormBuffer(elmInt,6),Sig,U,Theta,detJW);
+  return evalInt(static_cast<ElmNorm&>(elmInt),Sig,U,Theta,detJW);
 }
 
 
-bool ElasticityNormULMixed::evalBouMx (LocalIntegral*& elmInt,
+bool ElasticityNormULMixed::evalBouMx (LocalIntegral& elmInt,
 				       const MxFiniteElement& fe,
 				       const Vec3& X, const Vec3& normal) const
 {
   return this->ElasticityNormUL::evalBou(elmInt,fe,X,normal);
 }
+

Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h

@@ -57,24 +57,35 @@ public:
   //! \param[in] mode The solution mode to use
   virtual void setMode(SIM::SolutionMode mode);
 
+  //! \brief Returns a local integral container for the given element.
+  //! \param[in] nen1 Number of nodes on element for basis 1
+  //! \param[in] nen2 Number of nodes on element for basis 2
+  //! \param[in] neumann Whether or not we are assembling Neumann BC's
+  virtual LocalIntegral* getLocalIntegral(size_t nen1, size_t nen2, size_t,
+					  bool neumann) const;
+
   //! \brief Initializes current element for numerical integration.
   //! \param[in] MNPC1 Nodal point correspondance for the basis 1
   //! \param[in] MNPC2 Nodal point correspondance for the basis 2
   //! \param[in] n1 Number of nodes in basis 1 on this patch
+  //! \param elmInt The local integral object for current element
   virtual bool initElement(const std::vector<int>& MNPC1,
-                           const std::vector<int>& MNPC2, size_t n1);
+                           const std::vector<int>& MNPC2, size_t n1,
+			   LocalIntegral& elmInt);
   //! \brief Initializes current element for numerical boundary integration.
   //! \param[in] MNPC1 Nodal point correspondance for the basis 1
   //! \param[in] MNPC2 Nodal point correspondance for the basis 2
+  //! \param elmInt The local integral object for current element
   virtual bool initElementBou(const std::vector<int>& MNPC1,
-			      const std::vector<int>& MNPC2, size_t);
+			      const std::vector<int>& MNPC2, size_t,
+			      LocalIntegral& elmInt);
 
   //! \brief Evaluates the mixed field problem integrand at an interior point.
   //! \param elmInt The local integral object to receive the contributions
   //! \param[in] fe Mixed finite element data of current integration point
   //! \param[in] prm Nonlinear solution algorithm parameters
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool evalIntMx(LocalIntegral*& elmInt, const MxFiniteElement& fe,
+  virtual bool evalIntMx(LocalIntegral& elmInt, const MxFiniteElement& fe,
 			 const TimeDomain& prm, const Vec3& X) const;
 
   //! \brief Evaluates the integrand at a boundary point.
@@ -82,7 +93,7 @@ public:
   //! \param[in] fe Mixed finite element data of current integration point
   //! \param[in] X Cartesian coordinates of current integration point
   //! \param[in] normal Boundary normal vector at current integration point
-  virtual bool evalBouMx(LocalIntegral*& elmInt, const MxFiniteElement& fe,
+  virtual bool evalBouMx(LocalIntegral& elmInt, const MxFiniteElement& fe,
 			 const Vec3& X, const Vec3& normal) const;
 
   //! \brief Returns a pointer to an Integrand for solution norm evaluation.
@@ -102,10 +113,6 @@ public:
   //! \param[in] prefix Name prefix for all components
   virtual const char* getField1Name(size_t i, const char* prefix = 0) const;
 
-protected:
-  mutable Tensor Fbar; //!< Mixed model deformation gradient
-  mutable Matrix Dmat; //!< Projected mixed constitutive matrix
-
   friend class ElasticityNormULMixed;
 };
 
@@ -128,7 +135,7 @@ public:
   //! \param[in] fe Mixed finite element data of current integration point
   //! \param[in] prm Nonlinear solution algorithm parameters
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool evalIntMx(LocalIntegral*& elmInt, const MxFiniteElement& fe,
+  virtual bool evalIntMx(LocalIntegral& elmInt, const MxFiniteElement& fe,
 			 const TimeDomain& prm, const Vec3& X) const;
 
   //! \brief Evaluates the integrand at a boundary point (mixed).
@@ -136,7 +143,7 @@ public:
   //! \param[in] fe Mixed finite element data of current integration point
   //! \param[in] X Cartesian coordinates of current integration point
   //! \param[in] normal Boundary normal vector at current integration point
-  virtual bool evalBouMx(LocalIntegral*& elmInt, const MxFiniteElement& fe,
+  virtual bool evalBouMx(LocalIntegral& elmInt, const MxFiniteElement& fe,
 			 const Vec3& X, const Vec3& normal) const;
 };
 

Apps/FiniteDefElasticity/main_NonLinEl.C

@@ -25,7 +25,9 @@
 #include <string.h>
 #include <ctype.h>
 
+#ifndef USE_OPENMP
 extern std::vector<int> mixedDbgEl; //!< List of elements for additional output
+#endif
 
 
 /*!
@@ -146,9 +148,11 @@ int main (int argc, char** argv)
       else
 	--i;
     }
+#ifndef USE_OPENMP
     else if (!strcmp(argv[i],"-dbgElm"))
       while (i < argc-1 && isdigit(argv[i+1][0]))
 	utl::parseIntegers(mixedDbgEl,argv[++i]);
+#endif
     else if (!strcmp(argv[i],"-ignore"))
       while (i < argc-1 && isdigit(argv[i+1][0]))
 	utl::parseIntegers(ignoredPatches,argv[++i]);