Changed: Extended the Material::evaluate interface to take a

FiniteElement reference as argument (replacing integration point counter).
Added: Option to define the Young's modulus via a scalar field.

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

Elasticity.C

@@ -503,10 +503,11 @@ Vec3 Elasticity::evalSol (const Vector& eV, const Vector& N) const
 }
 
 
-bool Elasticity::formCinverse (Matrix& Cinv, const Vec3& X) const
+bool Elasticity::formCinverse (Matrix& Cinv, const FiniteElement& fe,
+                               const Vec3& X) const
 {
   SymmTensor dummy(nsd,axiSymmetry); double U;
-  return material->evaluate(Cinv,dummy,U,0,X,dummy,dummy,-1);
+  return material->evaluate(Cinv,dummy,U,fe,X,dummy,dummy,-1);
 }
 
 
@@ -573,7 +574,7 @@ bool Elasticity::evalSol (Vector& s, const Vector& eV, const FiniteElement& fe,
   // Calculate the stress tensor through the constitutive relation
   Matrix Cmat;
   SymmTensor sigma(nsd, axiSymmetry || material->isPlaneStrain()); double U;
-  if (!material->evaluate(Cmat,sigma,U,fe.iGP,X,dUdX,eps))
+  if (!material->evaluate(Cmat,sigma,U,fe,X,dUdX,eps))
     return false;
 
   // Congruence transformation to local coordinate system at current point
@@ -706,7 +707,8 @@ bool ElasticityNorm::evalInt (LocalIntegral& elmInt, const FiniteElement& fe,
 
   // Evaluate the inverse constitutive matrix at this point
   Matrix Cinv;
-  if (!problem.formCinverse(Cinv,X)) return false;
+  if (!problem.formCinverse(Cinv,fe,X))
+    return false;
 
   // Evaluate the finite element stress field
   Vector sigmah, sigma, error;

Elasticity.h

@@ -215,8 +215,9 @@ class Elasticity : public IntegrandBase
   //! \brief Sets up the inverse constitutive matrix at current point.
   //! \param[out] Cinv \f$6\times6\f$-matrix (in 3D) or \f$3\times3\f$-matrix
   //! (in 2D), representing the inverse constitutive tensor
+  //! \param[in] fe Finite element data at current point
   //! \param[in] X Cartesian coordinates of current point
-  bool formCinverse(Matrix& Cinv, const Vec3& X) const;
+  bool formCinverse(Matrix& Cinv, const FiniteElement& fe, const Vec3& X) const;
 
   //! \brief Returns \e true if this is an axial-symmetric problem.
   bool isAxiSymmetric() const { return axiSymmetry; }

LinIsotropic.C

@@ -12,22 +12,31 @@
 //==============================================================================
 
 #include "LinIsotropic.h"
+#include "Field.h"
 #include "Tensor.h"
 #include "Vec3.h"
 
 
 LinIsotropic::LinIsotropic (bool ps, bool ax) : planeStress(ps), axiSymmetry(ax)
 {
   // Default material properties - typical values for steel (SI units)
-  Efunc = 0;
+  Efunc = NULL;
+  Efield = NULL;
   Emod = 2.05e11;
   nu = 0.29;
   rho = 7.85e3;
 }
 
 
 LinIsotropic::LinIsotropic (RealFunc* E, double v, double den, bool ps, bool ax)
-  : Efunc(E), nu(v), rho(den), planeStress(ps), axiSymmetry(ax)
+  : Efunc(E), Efield(NULL), nu(v), rho(den), planeStress(ps), axiSymmetry(ax)
+{
+  Emod = -1.0; // Should not be referenced
+}
+
+
+LinIsotropic::LinIsotropic (Field* E, double v, double den, bool ps, bool ax)
+  : Efunc(NULL), Efield(E), nu(v), rho(den), planeStress(ps), axiSymmetry(ax)
 {
   Emod = -1.0; // Should not be referenced
 }
@@ -81,17 +90,21 @@ void LinIsotropic::print (std::ostream& os) const
   \end{array}\right] \f]
 */
 
-bool LinIsotropic::evaluate (Matrix& C, SymmTensor& sigma, double& U, size_t,
-			     const Vec3& X, const Tensor&,
-			     const SymmTensor& eps, char iop,
-			     const TimeDomain*, const Tensor*) const
+bool LinIsotropic::evaluate (Matrix& C, SymmTensor& sigma, double& U,
+                             const FiniteElement& fe, const Vec3& X,
+                             const Tensor&, const SymmTensor& eps, char iop,
+                             const TimeDomain*, const Tensor*) const
 {
   const size_t nsd = sigma.dim();
   const size_t nst = nsd == 2 && axiSymmetry ? 4 : nsd*(nsd+1)/2;
   C.resize(nst,nst,true);
 
-  // Evaluate the scalar stiffness function, if defined
-  const double E = Efunc ? (*Efunc)(X) : Emod;
+  // Evaluate the scalar stiffness function or field, if defined
+  double E = Emod;
+  if (Efield)
+    E = Efield->valueFE(fe);
+  else if (Efunc)
+    E = (*Efunc)(X);
 
   if (nsd == 1)
   {
@@ -101,7 +114,7 @@ bool LinIsotropic::evaluate (Matrix& C, SymmTensor& sigma, double& U, size_t,
     {
       sigma = eps; sigma *= E;
       if (iop == 3)
-	U = 0.5*sigma(1,1)*eps(1,1);
+        U = 0.5*sigma(1,1)*eps(1,1);
     }
     return true;
   }
@@ -172,11 +185,11 @@ bool LinIsotropic::evaluate (Matrix& C, SymmTensor& sigma, double& U, size_t,
       // Account for non-matching tensor dimensions
       SymmTensor epsil(sigma.dim(), nsd == 2 && axiSymmetry);
       if (!C.multiply(epsil=eps,sig))
-	return false;
+        return false;
     }
     else
       if (!C.multiply(eps,sig))
-	return false;
+        return false;
 
     sigma = sig; // Add sigma_zz in case of plane strain
     if (!planeStress && ! axiSymmetry && nsd == 2 && sigma.size() == 4)

LinIsotropic.h

@@ -16,6 +16,7 @@
 
 #include "MaterialBase.h"
 #include "Function.h"
+#include "Field.h"
 
 
 /*!
@@ -36,18 +37,27 @@ class LinIsotropic : public Material
   //! \param[in] ps If \e true, assume plane stress in 2D
   //! \param[in] ax If \e true, assume 3D axi-symmetric material
   LinIsotropic(double E, double v = 0.0, double densty = 0.0,
-	       bool ps = false, bool ax = false)
-    : Efunc(0), Emod(E), nu(v), rho(densty), planeStress(ps), axiSymmetry(ax) {}
+               bool ps = false, bool ax = false)
+    : Efunc(NULL), Efield(NULL), Emod(E), nu(v), rho(densty),
+      planeStress(ps), axiSymmetry(ax) {}
   //! \brief Constructor initializing the material parameters.
   //! \param[in] E Young's modulus (spatial function)
   //! \param[in] v Poisson's ratio
   //! \param[in] density Mass density
   //! \param[in] ps If \e true, assume plane stress in 2D
   //! \param[in] ax If \e true, assume 3D axi-symmetric material
   LinIsotropic(RealFunc* E, double v = 0.0, double density = 0.0,
-	       bool ps = false, bool ax = false);
+               bool ps = false, bool ax = false);
+  //! \brief Constructor initializing the material parameters.
+  //! \param[in] E Young's modulus (spatial field)
+  //! \param[in] v Poisson's ratio
+  //! \param[in] density Mass density
+  //! \param[in] ps If \e true, assume plane stress in 2D
+  //! \param[in] ax If \e true, assume 3D axi-symmetric material
+  LinIsotropic(Field* E, double v = 0.0, double density = 0.0,
+               bool ps = false, bool ax = false);
   //! \brief The destructor deletes the stiffness function, if defined.
-  virtual ~LinIsotropic() { delete Efunc; }
+  virtual ~LinIsotropic() { delete Efunc; delete Efield; }
 
   //! \brief Returns \e false if plane stress in 2D.
   virtual bool isPlaneStrain() const { return !planeStress; }
@@ -62,24 +72,29 @@ class LinIsotropic : public Material
   //! \param[out] C Constitutive matrix at current point
   //! \param[out] sigma Stress tensor at current point
   //! \param[out] U Strain energy density at current point
+  //! \param[in] fe Finite element quantities at current point
   //! \param[in] X Cartesian coordinates of current point
   //! \param[in] eps Strain tensor at current point
   //! \param[in] iop Calculation option;
   //!  -1 : Calculate the inverse constitutive matrix only,
   //!   0 : Calculate the constitutive matrix only,
   //!   1 : Calculate Cauchy stresses and the constitutive matrix,
   //!   3 : Calculate the strain energy density only.
-  virtual bool evaluate(Matrix& C, SymmTensor& sigma, double& U, size_t,
-			const Vec3& X, const Tensor&, const SymmTensor& eps,
-			char iop = 1, const TimeDomain* = 0,
-			const Tensor* = 0) const;
+  virtual bool evaluate(Matrix& C, SymmTensor& sigma, double& U,
+                        const FiniteElement& fe, const Vec3& X,
+                        const Tensor&, const SymmTensor& eps,
+                        char iop = 1, const TimeDomain* = NULL,
+                        const Tensor* = NULL) const;
 
   //! \brief Returns the function, if any, describing the stiffness variation.
   const RealFunc* getEfunc() const { return Efunc; }
+  //! \brief Returns the field, if any, describing the stiffness variation.
+  const Field* getEfield() const { return Efield; }
 
 protected:
   // Material properties
-  RealFunc* Efunc;    //!< Young's modules (spatial function)
+  RealFunc* Efunc;    //!< Young's modulus (spatial function)
+  Field* Efield;      //!< Young's modulus (spatial field)
   double Emod;        //!< Young's modulus (constant)
   double nu;          //!< Poisson's ratio
   double rho;         //!< Mass density

LinearElasticity.C

@@ -45,7 +45,7 @@ bool LinearElasticity::evalInt (LocalIntegral& elmInt, const FiniteElement& fe,
 
     // Evaluate the constitutive matrix and the stress tensor at this point
     double U;
-    if (!material->evaluate(Cmat,sigma,U,fe.iGP,X,eps,eps))
+    if (!material->evaluate(Cmat,sigma,U,fe,X,eps,eps))
       return false;
   }
 

MaterialBase.C

@@ -0,0 +1,25 @@
+// $Id$
+//==============================================================================
+//!
+//! \file MaterialBase.C
+//!
+//! \date Mar 01 2011
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Base class for material models.
+//!
+//==============================================================================
+
+#include "MaterialBase.h"
+#include "FiniteElement.h"
+
+
+bool Material::evaluate (Matrix& C, SymmTensor& sigma, double& U,
+                         size_t ip, const Vec3& X, const Tensor& F,
+                         const SymmTensor& eps, char iop,
+                         const TimeDomain* prm, const Tensor* Fpf) const
+{
+  FiniteElement fe(0,ip);
+  return this->evaluate(C,sigma,U,fe,X,F,eps,iop,prm,Fpf);
+}

MaterialBase.h

@@ -19,6 +19,8 @@
 class Vec3;
 class Tensor;
 class SymmTensor;
+class FiniteElement;
+class Field;
 struct TimeDomain;
 
 
@@ -50,6 +52,8 @@ class Material
   virtual void initResultPoints() {}
   //! \brief Defines a point location with some special material properties.
   virtual void addSpecialPoint(const Vec3&) {}
+  //! \brief Assigns a scalar field defining the material properties.
+  virtual void assignScalarField(Field*, size_t = 0) {}
 
   //! \brief Evaluates the mass density at current point.
   virtual double getMassDensity(const Vec3&) const { return 0.0; }
@@ -58,7 +62,7 @@ class Material
   //! \param[out] C Constitutive matrix at current point
   //! \param[out] sigma Stress tensor at current point
   //! \param[out] U Strain energy density at current point
-  //! \param[in] ip Global index for current point
+  //! \param[in] fe Finite element quantities at current point
   //! \param[in] X Cartesian coordinates of current point
   //! \param[in] F Deformation gradient at current point
   //! \param[in] eps Strain tensor at current point
@@ -70,10 +74,29 @@ class Material
   //!   3 : Calculate the strain energy density only.
   //! \param[in] prm Nonlinear solution algorithm parameters
   //! \param[in] Fpf Deformation gradient for push-forward transformation
-  virtual bool evaluate(Matrix& C, SymmTensor& sigma, double& U, size_t ip,
-			const Vec3& X, const Tensor& F, const SymmTensor& eps,
-			char iop = 1, const TimeDomain* prm = 0,
-			const Tensor* Fpf = 0) const = 0;
+  virtual bool evaluate(Matrix& C, SymmTensor& sigma, double& U,
+                        const FiniteElement& fe, const Vec3& X,
+                        const Tensor& F, const SymmTensor& eps,
+                        char iop = 1, const TimeDomain* prm = NULL,
+                        const Tensor* Fpf = NULL) const = 0;
+  //! \brief Evaluates the constitutive relation at an integration point.
+  //! \param[out] C Constitutive matrix at current point
+  //! \param[out] sigma Stress tensor at current point
+  //! \param[out] U Strain energy density at current point
+  //! \param[in] ip Global index for current point
+  //! \param[in] X Cartesian coordinates of current point
+  //! \param[in] F Deformation gradient at current point
+  //! \param[in] eps Strain tensor at current point
+  //! \param[in] iop Calculation option
+  //! \param[in] prm Nonlinear solution algorithm parameters
+  //! \param[in] Fpf Deformation gradient for push-forward transformation
+  //!
+  //! \details This is a wrapper used when parametric or local coordinates
+  //! of current point are not needed (or available).
+  bool evaluate(Matrix& C, SymmTensor& sigma, double& U, size_t ip,
+                const Vec3& X, const Tensor& F, const SymmTensor& eps,
+                char iop = 1, const TimeDomain* prm = NULL,
+                const Tensor* Fpf = NULL) const;
 
   //! \brief Returns number of internal result variables of the material model.
   virtual int getNoIntVariables() const { return 0; }