OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Ongoing work on plasticity solver. Including calculation of stresses at separate result points.

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@872 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo
2011-03-24 17:37:35 +00:00
committed by Knut Morten Okstad
parent be161ffc21
commit f4d6909e7d
19 changed files with 361 additions and 408 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Apps/FiniteDefElasticity/LinearMaterial.C
View File

@@ -43,7 +43,7 @@ bool LinearMaterial::evaluate (Matrix& C, SymmTensor& sigma, double& U,
// Push-forward the constitutive matrix to current configuration
size_t i, j, k, l;
size_t n = F.dim();
size_t n = eps.dim();
size_t m = C.rows();
Matrix T(m,m), Ctmp;

12
Apps/FiniteDefElasticity/LinearMaterial.h
View File

@@ -18,13 +18,21 @@
/*!
\brief Class representing a isotropic linear elastic material model.
\brief Class representing a general linear elastic material model.
\details This class is a wrapper for any linear-elastic material model,
when it is to be used in a nonlinear finite element formulation based on
the updated Lagrangian formulation. The \a evaluate method of this class
just invokes the corresponding method of the object pointed to by the
\a material member, and then transforms the resulting constitutive matrix
and stress tensor to the current updated reference frame, based on the
supplied deformation tensor, \a F.
*/
class LinearMaterial : public Material
{
public:
//! \brief The constructor initializes the material properties pointer.
//! \param[in] Pointer to material model for linear elastic analysis.
LinearMaterial(const Material* mat) : material(mat) {}
//! \brief The destructor deletes the material properties object.
virtual ~LinearMaterial() { delete material; }
@@ -53,7 +61,7 @@ public:
char iop = 1, const TimeDomain* prm = 0) const;
private:
const Material* material; //!< Pointer to actual material properties object
const Material* material; //!< Linear-elastic material properties object
};
#endif

38
Apps/FiniteDefElasticity/Material/cons2d.f
View File

@@ -1,29 +1,27 @@
subroutine cons2d (ipsw,iter,iwr,lfirst,mTYP,mVER,
& detF,F,Fp,pmat,HV,U,Sig,Cst,ierr)
subroutine cons2d (ipsw,iter,iwr,lfirst,mTYP,mVER,nHV,nTM,
& detF,nDF,F,pMAT,HV,U,Sig,Cst,ierr)
C
C CONS2D: 2D wrapper on the 3D material of FENRIS.
C
implicit none
C
integer ipsw, iter, iwr, lfirst, mTYP, mVER, ierr
real*8 detF, F(4), Fp(4), pmat(*), HV, U, Sig(3), Cst(3,3)
real*8 F3D(3,3), F3P(3,3), S3D(6), C3D(6,6)
integer ipsw, iter, iwr, lfirst, mTYP, mVER, nHV, nTM, nDF, ierr
real*8 detF, F(*), pMAT(*), HV(*), U, Sig(3), Cst(3,3)
real*8 F3D(3,3), S3D(6), C3D(6,6)
C
F3D = 0.0D0
F3D(1,1) = F(1)
F3D(2,1) = F(2)
F3D(1,2) = F(3)
F3D(2,2) = F(4)
F3D(3,3) = 1.0D0
F3P = 0.0D0
F3P(1,1) = Fp(1)
F3P(2,1) = Fp(2)
F3P(1,2) = Fp(3)
F3P(2,2) = Fp(4)
F3P(3,3) = 1.0D0
C
call cons3d (ipsw,iter,iwr,lfirst,mTYP,mVER,
& detF,F3D,F3P,pmat,HV,U,S3D,C3D,ierr)
if (ndf .eq. 2) then
F3D = 0.0D0
F3D(1,1) = F(1)
F3D(2,1) = F(2)
F3D(1,2) = F(3)
F3D(2,2) = F(4)
F3D(3,3) = 1.0D0
call cons3d (ipsw,iter,iwr,lfirst,mTYP,mVER,nHV,nTM,
& detF,F3D,pMAT,HV,U,S3D,C3D,ierr)
else
call cons3d (ipsw,iter,iwr,lfirst,mTYP,mVER,nHV,nTM,
& detF,F,pMAT,HV,U,S3D,C3D,ierr)
end if
C
Cst(1:2,1:2) = C3D(1:2,1:2)
Cst(1:2,3) = C3D(1:2,4)

40
Apps/FiniteDefElasticity/Material/plas2d.f
View File

@@ -1,30 +1,34 @@
subroutine plas2d (ipsw,iwr,iter,lfirst,ntm,pmat,detF,
& Fn1,Fn,be,Epp,Epl,Sig,Cst,ierr)
& ndf,Fn1,Fn,be,Epp,Epl,Sig,Cst,ierr)
C
C PLAS2D: 2D wrapper on the 3D plastic material of FENRIS.
C
implicit none
C
integer ipsw, iwr, iter, lfirst, ntm, ierr
real*8 detF, Fn1(4), Fn(4), pmat(*), be(*)
integer ipsw, iwr, iter, lfirst, ntm, ndf, ierr
real*8 detF, Fn1(*), Fn(*), pmat(*), be(*)
real*8 Epp, Epl(*), Sig(3), Cst(3,3)
real*8 F3D(3,3), F3P(3,3), S3D(6), C3D(6,6)
C
F3D = 0.0D0
F3D(1,1) = Fn(1)
F3D(2,1) = Fn(2)
F3D(1,2) = Fn(3)
F3D(2,2) = Fn(4)
F3D(3,3) = 1.0D0
F3P = 0.0D0
F3P(1,1) = Fn1(1)
F3P(2,1) = Fn1(2)
F3P(1,2) = Fn1(3)
F3P(2,2) = Fn1(4)
F3P(3,3) = 1.0D0
C
call plas3d (ipsw,iwr,iter,lfirst,ntm,pmat,detF,
& F3P,F3D,be,Epp,Epl,S3D,C3D,ierr)
if (ndf .eq. 2) then
F3D = 0.0D0
F3D(1,1) = Fn(1)
F3D(2,1) = Fn(2)
F3D(1,2) = Fn(3)
F3D(2,2) = Fn(4)
F3D(3,3) = 1.0D0
F3P = 0.0D0
F3P(1,1) = Fn1(1)
F3P(2,1) = Fn1(2)
F3P(1,2) = Fn1(3)
F3P(2,2) = Fn1(4)
F3P(3,3) = 1.0D0
call plas3d (ipsw,iwr,iter,lfirst,ntm,pmat,detF,
& F3P,F3D,be,Epp,Epl,S3D,C3D,ierr)
else
call plas3d (ipsw,iwr,iter,lfirst,ntm,pmat,detF,
& Fn,Fn1,be,Epp,Epl,S3D,C3D,ierr)
end if
C
Cst(1:2,1:2) = C3D(1:2,1:2)
Cst(1:2,3) = C3D(1:2,4)

13
Apps/FiniteDefElasticity/NeoHookeMaterial.C
View File

@@ -20,14 +20,15 @@ extern "C" {
void cons2d_(const int& ipsw, const int& iter, const int& iwr,
const int& lfirst, const int& mTYP, const int& mVER,
const int& nHV, const int& nTM, const double& detF,
const double* F, const double* pMAT, double* HV, double& Engy,
const double* Sig, double* Cst, int& ierr);
const int& nDF,
const double* F, const double* pMAT, double* HV,
double& Engy, const double* Sig, double* Cst, int& ierr);
//! \brief Interface to 3D nonlinear material routines (FORTRAN-77 code).
void cons3d_(const int& ipsw, const int& iter, const int& iwr,
const int& lfirst, const int& mTYP, const int& mVER,
const int& nHV, const int& nTM, const double& detF,
const double* F, const double* pMAT, double* HV, double& Engy,
const double* Sig, double* Cst, int& ierr);
const double* F, const double* pMAT, double* HV,
double& Engy, const double* Sig, double* Cst, int& ierr);
}
#ifndef INT_DEBUG
#define INT_DEBUG 0
@@ -83,14 +84,14 @@ bool NeoHookeMaterial::evaluate (Matrix& C, SymmTensor& sigma, double& U,
return false;
}
size_t ndim = F.dim();
size_t ndim = sigma.dim();
size_t ncmp = ndim*(ndim+1)/2;
C.resize(ncmp,ncmp);
int ierr = -99;
#ifdef USE_FTNMAT
// Invoke the FORTRAN routine for Neo-Hookean hyperelastic material models.
if (ndim == 2)
cons2d_(INT_DEBUG,0,6,0,mTYP,mVER,0,3,J,F.ptr(),pmat,
cons2d_(INT_DEBUG,0,6,0,mTYP,mVER,0,3,J,F.dim(),F.ptr(),pmat,
&U,U,sigma.ptr(),C.ptr(),ierr);
else
cons3d_(INT_DEBUG,0,6,0,mTYP,mVER,0,6,J,F.ptr(),pmat,

5
Apps/FiniteDefElasticity/NeoHookeMaterial.h
View File

@@ -19,6 +19,11 @@
/*!
\brief Class representing a Neo-Hookean hyperelastic material model.
\details This class is a wrapper for the FORTRAN material routines of FENRIS,
implemented in by the subroutine CONS3D and subroutines called from CONS3D.
Which material routine to use is governed by the \a mVER parameter.
Only isotropic material properties are supported by this class.
\note In 2D, only plane strain is supported by this class.
*/
class NeoHookeMaterial : public LinIsotropic

57
Apps/FiniteDefElasticity/NonlinearElasticityUL.C
View File

@@ -15,6 +15,7 @@
#include "MaterialBase.h"
#include "ElmMats.h"
#include "ElmNorm.h"
#include "TimeDomain.h"
#include "Tensor.h"
#include "Vec3Oper.h"
@@ -85,7 +86,22 @@ void NonlinearElasticityUL::setMode (SIM::SolutionMode mode)
}
bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt, double detJW,
void NonlinearElasticityUL::initIntegration (const TimeDomain& prm)
{
if (material)
material->initIntegration(prm);
}
void NonlinearElasticityUL::initResultPoints ()
{
if (material)
material->initResultPoints();
}
bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt,
const TimeDomain& prm, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
@@ -128,7 +144,7 @@ bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt, double detJW,
if (eKm || eKg || iS)
{
double U = 0.0;
if (!material->evaluate(Cmat,sigma,U,X,F,E, (eKg || iS) && lHaveStrains))
if (!material->evaluate(Cmat,sigma,U,X,F,E,(eKg || iS),&prm))
return false;
}
@@ -226,7 +242,7 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt, double detJW,
bool NonlinearElasticityUL::formDefGradient (const Matrix& dNdX,
Tensor& F) const
{
SymmTensor dummy(0);
static SymmTensor dummy(0);
return this->kinematics(dNdX,F,dummy);
}
@@ -250,6 +266,7 @@ bool NonlinearElasticityUL::kinematics (const Matrix& dNdX,
return false;
}
// Compute the deformation gradient, [F] = [I] + [dudX] = [I] + [dNdX]*[u],
// and the Green-Lagrange strains, E_ij = 0.5(F_ij+F_ji+F_ki*F_kj).
@@ -257,12 +274,17 @@ bool NonlinearElasticityUL::kinematics (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
F = dUdX;
else
if (!dUdX.multiplyMat(*eV,dNdX)) // dUdX = Grad{u} = eV*dNdX
return false;
unsigned short int i, j, k;
if (dUdX.rows() < F.dim())
F.zero();
// Cannot use operator= here, in case F is of higher dimension than dUdX
for (i = 1; i <= dUdX.rows(); i++)
for (j = 1; j <= dUdX.cols(); j++)
F(i,j) = dUdX(i,j);
// Now form the Green-Lagrange strain tensor.
// Note that for the shear terms (i/=j) we actually compute 2*E_ij
@@ -294,7 +316,14 @@ NormBase* NonlinearElasticityUL::getNormIntegrand (AnaSol*) const
}
bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt, double detJW,
void ElasticityNormUL::initIntegration (const TimeDomain& prm)
{
problem.initIntegration(prm);
}
bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt,
const TimeDomain& prm, double detJW,
const Vector&, const Matrix& dNdX,
const Vec3& X) const
{
@@ -312,18 +341,12 @@ bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt, double detJW,
Matrix C;
SymmTensor S(E.dim());
double U = 0.0;
if (!ulp->material->evaluate(C,S,U,X,F,E,2))
if (!ulp->material->evaluate(C,S,U,X,F,E,2,&prm))
return false;
if (U == 0.0)
{
// Integrate the energy norm a(u^h,u^h) = Int_Omega0 (S:E) dV0
const RealArray& sigma = S;
const RealArray& epsil = E;
for (size_t i = 0; i < sigma.size(); i++)
U += sigma[i]*epsil[i];
}
// Integrate the energy norm a(u^h,u^h) = Int_Omega0 (S:E) dV0
if (U == 0.0) U = S.innerProd(E);
pnorm[0] += U*detJW;
return true;
}

24
Apps/FiniteDefElasticity/NonlinearElasticityUL.h
View File

@@ -43,13 +43,21 @@ public:
//! \param[in] mode The solution mode to use
virtual void setMode(SIM::SolutionMode mode);
//! \brief Initializes the integrand for a new integration loop.
//! \param[in] prm Nonlinear solution algorithm parameters
virtual void initIntegration(const TimeDomain& prm);
//! \brief Initializes the integrand for a result point loop.
virtual void initResultPoints();
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] prm Nonlinear solution algorithm parameters
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& prm, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
@@ -115,6 +123,10 @@ public:
//! \brief Empty destructor.
virtual ~ElasticityNormUL() {}
//! \brief Initializes the integrand for a new integration loop.
//! \param[in] prm Nonlinear solution algorithm parameters
virtual void initIntegration(const TimeDomain& prm);
//! \brief Initializes current element for numerical integration (mixed).
//! \param[in] MNPC1 Nodal point correspondance for the displacement field
virtual bool initElement(const std::vector<int>& MNPC1,
@@ -132,25 +144,29 @@ public:
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] prm Nonlinear solution algorithm parameters
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& prm, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Evaluates the integrand at an interior point (mixed).
//! \param elmInt The local integral object to receive the contributions
//! \param[in] prm Nonlinear solution algorithm parameters
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N1 Basis function values for the displacement field
//! \param[in] dN1dX Basis function gradients for the displacement field
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& prm, double detJW,
const Vector& N1, const Vector&,
const Matrix& dN1dX, const Matrix&,
const Vec3& X) const
{
return this->evalInt(elmInt,detJW,N1,dN1dX,X);
return this->evalInt(elmInt,prm,detJW,N1,dN1dX,X);
}
//! \brief Evaluates the integrand at a boundary point (mixed).

40
Apps/FiniteDefElasticity/NonlinearElasticityULMX.C
View File

@@ -15,6 +15,7 @@
#include "MaterialBase.h"
#include "ElmMats.h"
#include "ElmNorm.h"
#include "TimeDomain.h"
#include "Utilities.h"
#include "Vec3Oper.h"
@@ -133,19 +134,13 @@ bool NonlinearElasticityULMX::initElement (const std::vector<int>& MNPC,
size_t nPm = utl::Pascal(p,nsd);
if (Hh) Hh->resize(nPm,nPm,true);
// Extract the element level displacements of previous increment
int ierr = 0;
if (!primsol[1].empty())
if ((ierr = utl::gather(MNPC,nsd,primsol[1],myMats->b[2])))
std::cerr <<" *** NonlinearElasticityULMX::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
// The other element matrices are initialized by the parent class method
return this->NonlinearElasticityUL::initElement(MNPC) && ierr == 0;
return this->NonlinearElasticityUL::initElement(MNPC);
}
bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt, double detJW,
bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt,
const TimeDomain&, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
@@ -172,7 +167,9 @@ bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt, double detJW,
// Invert the deformation gradient ==> Fi
Matrix Fi(nsd,nsd);
Fi.fill(ptData.F.ptr());
for (unsigned short int i = 1; i <= nsd; i++)
for (unsigned short int j = 1; j <= nsd; j++)
Fi(i,j) = ptData.F(i,j);
double J = Fi.inverse();
if (J == 0.0) return false;
@@ -207,7 +204,8 @@ bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt, double detJW,
}
bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt)
bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
const TimeDomain& prm)
{
if (!iS && !eKm && !eKg) return true;
if (!Hh) return false;
@@ -365,21 +363,20 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt)
// Evaluate the constitutive relation
double U = 0.0;
Sig.push_back(SymmTensor(3));
if (!material->evaluate(Cmat,Sig.back(),U,pt.X,pt.F,Sig.back()))
if (!material->evaluate(Cmat,Sig.back(),U,pt.X,pt.F,Sig.back(),true,&prm))
return false;
#ifdef USE_FTNMAT
// Project the constitutive matrix for the mixed model
D[iP].resize(7,7);
int ipsw = INT_DEBUG > 1 ? 9 : 0;
pcst3d_(Cmat.ptr(),D[iP].ptr(),ipsw,6);
pcst3d_(Cmat.ptr(),D[iP].ptr(),INT_DEBUG,6);
#else
std::cerr <<" *** NonlinearElasticityULMX::finalizeElement: "
<<" Does not work when compiled without USE_FTNMAT"<< std::endl;
return false;
#endif
if (U > 0.0)
if (!Sig.back().isZero())
{
// Integrate mean pressure
lHaveStress = true;
@@ -417,8 +414,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt)
if (eKm)
{
#ifdef USE_FTNMAT
int ipsw = INT_DEBUG > 1 ? 9 : 0;
mdma3d_(Bpres,Press,Sig[iP].ptr(),D[iP].ptr(),ipsw,6);
mdma3d_(Bpres,Press,Sig[iP].ptr(),D[iP].ptr(),INT_DEBUG,6);
// Integrate the material stiffness matrix
D[iP] *= dFmix*pt.detJW;
@@ -473,16 +469,18 @@ bool ElasticityNormULMX::initElement (const std::vector<int>& MNPC,
}
bool ElasticityNormULMX::evalInt (LocalIntegral*& elmInt, double detJW,
bool ElasticityNormULMX::evalInt (LocalIntegral*& elmInt,
const TimeDomain& prm, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
NonlinearElasticityULMX& elp = static_cast<NonlinearElasticityULMX&>(problem);
return elp.evalInt(elmInt,detJW,N,dNdX,X);
return elp.evalInt(elmInt,prm,detJW,N,dNdX,X);
}
bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt)
bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
const TimeDomain& prm)
{
NonlinearElasticityULMX& elp = static_cast<NonlinearElasticityULMX&>(problem);
if (!elp.Hh) return false;
@@ -574,7 +572,7 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt)
// Compute the strain energy density
double U = 0.0;
if (!elp.material->evaluate(C,Sig,U,pt.X,pt.F,Sig,false))
if (!elp.material->evaluate(C,Sig,U,pt.X,pt.F,Sig,false,&prm))
return false;
// Integrate strain energy

14
Apps/FiniteDefElasticity/NonlinearElasticityULMX.h
View File

@@ -63,6 +63,7 @@ public:
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] prm Nonlinear solution algorithm parameters
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
@@ -72,7 +73,8 @@ 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, double detJW,
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& prm, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
@@ -82,7 +84,8 @@ public:
//! modes. All data needed during the integration has been stored internally
//! in the \a myData member by the \a evalInt method.
//! \param elmInt The local integral object to receive the contributions
virtual bool finalizeElement(LocalIntegral*&);
//! \param[in] prm Nonlinear solution algorithm parameters
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
@@ -128,6 +131,7 @@ public:
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] prm Nonlinear solution algorithm parameters
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
@@ -137,7 +141,8 @@ 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, double detJW,
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& prm, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
@@ -146,7 +151,8 @@ public:
//! 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
virtual bool finalizeElement(LocalIntegral*& elmInt);
//! \param[in] prm Nonlinear solution algorithm parameters
virtual bool finalizeElement(LocalIntegral*& elmInt, const TimeDomain& prm);
//! \brief Returns which integrand to be used.
virtual int getIntegrandType() const { return 4; }

18
Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C
View File

@@ -190,7 +190,7 @@ const Vector& NonlinearElasticityULMixed::MixedElmMats::getRHSVector () const
NonlinearElasticityULMixed::NonlinearElasticityULMixed (unsigned short int n)
: NonlinearElasticityUL(n), Fbar(n), Dmat(7,7)
: NonlinearElasticityUL(n), Fbar(3), Dmat(7,7)
{
if (myMats) delete myMats;
myMats = new MixedElmMats();
@@ -290,7 +290,8 @@ bool NonlinearElasticityULMixed::initElementBou (const std::vector<int>& MNPC1,
}
bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, double detJW,
bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt,
const TimeDomain& prm, double detJW,
const Vector& N1,
const Vector& N2,
const Matrix& dNdX1,
@@ -298,8 +299,9 @@ bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, double detJW,
const Vec3& X) const
{
// Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
Tensor F(nsd);
SymmTensor E(nsd);
if (!this->kinematics(dNdX1,Fbar,E))
if (!this->kinematics(dNdX1,F,E))
return false;
bool lHaveStrains = !E.isZero();
@@ -317,10 +319,8 @@ bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, double detJW,
double dVol = Theta*detJW;
// Compute the mixed model deformation gradient, F_bar
double r1 = pow(fabs(Theta/Fbar.det()),1.0/3.0);
Matrix Fi(nsd,nsd);
Fi.fill(Fbar.ptr());
Fbar = F; // notice that F_bar always has dimension 3
double r1 = pow(fabs(Theta/F.det()),1.0/3.0);
Fbar *= r1;
if (nsd == 2) // In 2D we always assume plane strain so set F(3,3)=1
@@ -329,6 +329,8 @@ bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, double detJW,
return false;
// Invert the deformation gradient ==> Fi
Matrix Fi(nsd,nsd);
Fi.fill(F.ptr());
double J = Fi.inverse();
if (J == 0.0) return false;
@@ -351,7 +353,7 @@ bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, double detJW,
// Evaluate the constitutive relation
SymmTensor Eps(3), Sig(3), Sigma(nsd);
double U, Bpres = 0.0, Mpres = 0.0;
if (!material->evaluate(Cmat,Sig,U,X,Fbar,Eps=E,lHaveStrains))
if (!material->evaluate(Cmat,Sig,U,X,Fbar,Eps=E,lHaveStrains,&prm))
return false;
#ifdef USE_FTNMAT

4
Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h
View File

@@ -70,13 +70,15 @@ public:
//! \brief Evaluates the mixed field problem integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] prm Nonlinear solution algorithm parameters
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N1 Basis function values, field 1 (displacements)
//! \param[in] N2 Basis function values, field 2 (pressure and volume change)
//! \param[in] dN1dX Basis function gradients, field 1
//! \param[in] dN2dX Basis function gradients, field 2
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& prm, double detJW,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X) const;

145
Apps/FiniteDefElasticity/PlasticMaterial.C
View File

@@ -19,7 +19,8 @@ extern "C" {
//! \brief Interface to 2D elasto-plastic material routines (FORTRAN-77 code).
void plas2d_(const int& ipsw, const int& iwr, const int& iter,
const int& lfirst, const int& ntm, const double* pMAT,
const double& detF, const double* Fn1, const double* Fn,
const double& detF, const int& ndF,
const double* Fn1, const double* Fn,
const double* be, const double& Epp, const double* Epl,
const double* Sig, double* Cst, int& ierr);
//! \brief Interface to 3D elasto-plastic material routines (FORTRAN-77 code).
@@ -36,10 +37,10 @@ extern "C" {
/*!
See the Fortran routine Material/plas3d.f for interpretation of \a pMAT.
See the Fortran the routine Material/plas3d.f for interpretation of \a pMAT.
*/
PlasticPrm::PlasticPrm (const RealArray& p) : pMAT(p)
PlasticMaterial::PlasticMaterial (const RealArray& p) : pMAT(p), iP1(0), iP2(0)
{
if (pMAT.size() < 11) pMAT.resize(11,0.0);
@@ -57,13 +58,13 @@ PlasticPrm::PlasticPrm (const RealArray& p) : pMAT(p)
else if (nu < 0.5)
{
// Calculate the Bulk and Shear moduli from E and nu
Smod = Emod / (2.0 + 2.0*nu);
Bmod = Emod / (3.0 - 6.0*nu);
Smod = Emod / (2.0 + 2.0*nu);
}
else
{
Smod = Emod / 3.0;
Bmod = 0.0;
Smod = Emod / 3.0;
}
pMAT[0] = Emod;
@@ -73,7 +74,17 @@ PlasticPrm::PlasticPrm (const RealArray& p) : pMAT(p)
}
void PlasticPrm::print (std::ostream& os) const
PlasticMaterial::~PlasticMaterial ()
{
size_t i;
for (i = 0; i < itgPoints.size(); i++)
delete itgPoints[i];
for (i = 0; i < resPoints.size(); i++)
delete resPoints[i];
}
void PlasticMaterial::print (std::ostream& os) const
{
std::cout <<"PlasticMaterial: pMAT =";
for (size_t i = 0; i < pMAT.size(); i++)
@@ -82,17 +93,103 @@ void PlasticPrm::print (std::ostream& os) const
}
PlasticMaterial::PlasticMaterial (const PlasticPrm* prm, unsigned short int n)
: pMAT(prm->pMAT), Fp(n)
void PlasticMaterial::initIntegration (const TimeDomain& prm)
{
memset(HVc,0,sizeof(HVc));
memset(HVp,0,sizeof(HVc));
#if INT_DEBUG > 0
std::cout <<"PlasticMaterial::initIntegration: "<< iP1 << std::endl;
#endif
iP1 = 0;
iAmIntegrating = true;
if (prm.it > 0 || prm.first) return;
for (size_t i = 0; i < itgPoints.size(); i++)
itgPoints[i]->updateHistoryVars();
#if INT_DEBUG > 0
std::cout <<"PlasticMaterial::initIntegration: History updated "
<< itgPoints.size() << std::endl;
#endif
}
void PlasticMaterial::initResultPoints ()
{
#if INT_DEBUG > 0
std::cout <<"PlasticMaterial::initResultPoints: "<< iP1 << std::endl;
#endif
iP2 = 0;
iAmIntegrating = false;
for (size_t i = 0; i < resPoints.size(); i++)
resPoints[i]->updateHistoryVars();
#if INT_DEBUG > 0
std::cout <<"PlasticMaterial::initResultPoints: History updated "
<< resPoints.size() << std::endl;
#endif
}
bool PlasticMaterial::evaluate (Matrix& C, SymmTensor& sigma, double&,
const Vec3&, const Tensor& F, const SymmTensor&,
char, const TimeDomain* prm) const
{
if (iAmIntegrating)
{
if (!prm)
{
std::cerr <<" *** PlasticMaterial::evaluate: No time domain."<< std::endl;
return false;
}
while (itgPoints.size() <= iP1)
itgPoints.push_back(new PlasticPoint(this,F.dim()));
if (prm->it == 0 && !prm->first)
itgPoints[iP1]->Fp = F;
return itgPoints[iP1++]->evaluate(C,sigma,F,*prm);
}
else // Result evaluation
{
while (resPoints.size() <= iP2)
resPoints.push_back(new PlasticPoint(this,F.dim()));
if (!resPoints[iP2]->evaluate(C,sigma,F,TimeDomain(0,false)))
return false;
// Assume only one evaluation per increment; always update Fp
resPoints[iP2++]->Fp = F;
return true;
}
}
PlasticMaterial::PlasticPoint::PlasticPoint (const PlasticMaterial* prm,
unsigned short int n)
: pMAT(prm->pMAT), Fp(n)
{
// Initialize the history variables
memset(HVc,0,sizeof(HVc));
memset(HVp,0,sizeof(HVp));
Fp = HVp[0] = HVp[1] = HVp[2] = 1.0;
}
void PlasticMaterial::PlasticPoint::updateHistoryVars ()
{
// Update history variables with values of the new converged solution
memcpy(HVp,HVc,sizeof(HVc));
}
bool PlasticMaterial::PlasticPoint::evaluate (Matrix& C,
SymmTensor& sigma,
const Tensor& F,
const TimeDomain& prm) const
{
double J = F.det();
if (J == 0.0)
@@ -101,25 +198,23 @@ bool PlasticMaterial::evaluate (Matrix& C, SymmTensor& sigma, double&,
<<" Singular/zero deformation gradient\n"<< F;
return false;
}
else if (!prm)
{
std::cerr <<" *** PlasticMaterial::evaluate: No time domain."<< std::endl;
return false;
}
size_t ndim = F.dim();
// Restore history variables from the previous, converged configuration
memcpy(const_cast<double*>(HVc),HVp,sizeof(HVc));
size_t ndim = sigma.dim();
size_t ncmp = ndim*(ndim+1)/2;
C.resize(ncmp,ncmp);
int ierr = -99;
#ifdef USE_FTNMAT
// Invoke the FORTRAN routine for plasticity material models.
memcpy(const_cast<double*>(HVc),HVp,sizeof(HVc));
// Invoke the FORTRAN routine for plasticity material models
if (ndim == 2)
plas2d_(INT_DEBUG,6,prm->it,prm->first,4,&pMAT.front(),J,F.ptr(),Fp.ptr(),
HVc,HVc[6],HVc+7,sigma.ptr(),C.ptr(),ierr);
plas2d_(INT_DEBUG,6,prm.it,prm.first,4,&pMAT.front(),J,F.dim(),
F.ptr(),Fp.ptr(),HVc,HVc[6],HVc+7,sigma.ptr(),C.ptr(),ierr);
else
plas3d_(INT_DEBUG,6,prm->it,prm->first,6,&pMAT.front(),J,F.ptr(),Fp.ptr(),
HVc,HVc[6],HVc+7,sigma.ptr(),C.ptr(),ierr);
plas3d_(INT_DEBUG,6,prm.it,prm.first,6,&pMAT.front(),J,
F.ptr(),Fp.ptr(),HVc,HVc[6],HVc+7,sigma.ptr(),C.ptr(),ierr);
#if INT_DEBUG > 0
std::cout <<"PlasticMaterial::sigma =\n"<< sigma;
std::cout <<"PlasticMaterial::C ="<< C;
@@ -130,9 +225,3 @@ bool PlasticMaterial::evaluate (Matrix& C, SymmTensor& sigma, double&,
return ierr == 0;
}
void PlasticMaterial::updateHistoryVars ()
{
memcpy(HVp,HVc,sizeof(HVc));
}

120
Apps/FiniteDefElasticity/PlasticMaterial.h
View File

@@ -17,55 +17,84 @@
#include "MaterialBase.h"
#include "Tensor.h"
/*!
\brief Class containing parameters of an elasto-plastic material model.
*/
class PlasticPrm : public Material
{
public:
//! \brief Constructor initializing the material parameters.
PlasticPrm(const RealArray&);
//! \brief Empty destructor.
virtual ~PlasticPrm() {}
//! \brief Prints out material parameters to the given output stream.
virtual void print(std::ostream&) const;
//! \brief Evaluates the mass density at current point.
virtual double getMassDensity(const Vec3&) const { return pMAT[3]; }
//! \brief Dummy method (should not be invoked).
virtual bool evaluate(Matrix&, SymmTensor&, double&,
const Vec3&, const Tensor&, const SymmTensor&,
char, const TimeDomain*) const { return false; }
private:
RealArray pMAT; //!< Material property parameters
friend class PlasticMaterial;
};
class PlasticMaterial;
/*!
\brief Class representing an elasto-plastic material point.
\brief Class representing an history-dependent elasto-plastic material model.
\details This class is a wrapper for the plasticity routines of FENRIS,
implemented in by the FORTRAN subroutine PLAS3D.
The plasticity models have history variables at the integration points that
need to be "remembered" from one iteration/increment to the next. This is
maintained inside this class as a vector of integration point objects, and
a global counter keeping track of which integration point we are calculating.
This counter is initialized in the beginning of each iteration, and then
incremented by each invokation of the \a evaluate method. Therefore, it is
paramount that this method is invoked only once per integration point per
iteration, and in the same order in every iteration. Otherwise it will
output incorrect results.
A separate "integration point" vector is dedicated for results points.
This is needed because the results points (for visualization, etc.) are not
the same as the integration points used in the tangent evaluation.
*/
class PlasticMaterial : public Material
{
//! \brief Class representing an elasto-plastic material point.
class PlasticPoint
{
public:
//! \brief Constructor initializing the material parameters.
//! \param[in] prm Pointer to actual material model object.
//! \param[in] n Number of space dimensions
PlasticPoint(const PlasticMaterial* prm, unsigned short int n);
//! \brief Updates the internal history variables after convergence.
void updateHistoryVars();
//! \brief Evaluates the constitutive relation at this point.
//! \param[out] C Constitutive matrix at current point
//! \param[out] sigma Stress tensor at current point
//! \param[in] F Deformation gradient at current point
//! \param[in] prm Nonlinear solution algorithm parameters
bool evaluate(Matrix& C, SymmTensor& sigma,
const Tensor& F, const TimeDomain& prm) const;
private:
const RealArray& pMAT; //!< Material property parameters
double HVc[10]; //!< History variables, current configuration
double HVp[10]; //!< History variables, previous configuration
public:
Tensor Fp; //!< Deformation gradient, previous configuration
};
public:
//! \brief Constructor initializing the material parameters.
//! \param[in] prm Pointer to actual material parameters object.
//! \param[in] n Number of space dimensions
PlasticMaterial(const PlasticPrm* prm, unsigned short int n);
//! \brief Empty destructor.
virtual ~PlasticMaterial() {}
PlasticMaterial(const RealArray&);
//! \brief The destructor frees the dynamically allocated data.
virtual ~PlasticMaterial();
//! \brief Prints out material parameters to the given output stream.
virtual void print(std::ostream&) const;
//! \brief Initializes the material model for a new integration loop.
//! \param[in] prm Nonlinear solution algorithm parameters
//!
//! \details This method is invoked once before starting the numerical
//! integration over the entire spatial domain, and is used to reset the
//! global integration point counter, and to update the history variables.
virtual void initIntegration(const TimeDomain& prm);
//! \brief Initializes the material model for a new result point loop.
virtual void initResultPoints();
//! \brief Evaluates the mass density at current point.
virtual double getMassDensity(const Vec3&) const { return pMAT[3]; }
//! \brief Evaluates the constitutive relation at an integration point.
//! \brief Evaluates the constitutive relation at current integration point.
//! \param[out] C Constitutive matrix at current point
//! \param[out] sigma Stress tensor at current point
//! \param[in] F Deformation gradient at current point
@@ -74,18 +103,17 @@ public:
const Vec3&, const Tensor& F, const SymmTensor&,
char, const TimeDomain* prm) const;
//! \brief Returns a reference to the previous deformation gradient.
Tensor& defGrad() { return Fp; }
//! \brief Updates the internal history variables after convergence.
void updateHistoryVars();
private:
const RealArray& pMAT; //!< Material property parameters
friend class PlasticPoint;
double HVc[10]; //!< History variables, current configuration
double HVp[10]; //!< History variables, previous configuration
Tensor Fp; //!< Deformation gradient, previous configuration
RealArray pMAT; //!< Material property parameters
mutable size_t iP1; //!< Global integration point counter
mutable size_t iP2; //!< Global result point counter
bool iAmIntegrating; //!< Flag indicating integration or result evaluation
mutable std::vector<PlasticPoint*> itgPoints; //!< Integration point data
mutable std::vector<PlasticPoint*> resPoints; //!< Result point data
};
#endif

144
Apps/FiniteDefElasticity/PlasticityUL.C
View File

@@ -1,144 +0,0 @@
// $Id$
//==============================================================================
//!
//! \file PlasticityUL.C
//!
//! \date Mar 16 2011
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for finite deformation plasticity problems.
//!
//==============================================================================
#include "PlasticityUL.h"
#include "PlasticMaterial.h"
#include "ElmMats.h"
#include "TimeDomain.h"
PlasticityUL::PlasticityUL (unsigned short int n, char lop)
: NonlinearElasticityUL(n,lop)
{
iP = 0;
// Both the current and previous solutions are needed
primsol.resize(2);
}
PlasticityUL::~PlasticityUL ()
{
for (size_t i = 0; i < pBuf.size(); i++)
delete pBuf[i];
}
void PlasticityUL::print (std::ostream& os) const
{
std::cout <<"PlasticityUL: Finite deformation plasticity"<< std::endl;
this->NonlinearElasticityUL::print(os);
}
void PlasticityUL::initIntegration (const TimeDomain& prm)
{
iP = 0;
if (prm.it == 0 && !prm.first)
for (size_t i = 0; i < pBuf.size(); i++)
pBuf[i]->updateHistoryVars();
}
bool PlasticityUL::evalInt (LocalIntegral*& elmInt,
const TimeDomain& prm, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
while (pBuf.size() <= iP)
pBuf.push_back(new PlasticMaterial(static_cast<PlasticPrm*>(material),nsd));
PlasticMaterial* ptData = pBuf[iP++];
if (prm.it == 0)
{
// Evaluate the deformation gradient, Fp, at previous configuration
const_cast<PlasticityUL*>(this)->eV = &myMats->b[2];
if (!this->formDefGradient(dNdX,ptData->defGrad()))
return false;
}
// Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
Tensor F(nsd);
SymmTensor E(nsd);
const_cast<PlasticityUL*>(this)->eV = &myMats->b[1];
if (!this->kinematics(dNdX,F,E))
return false;
bool lHaveStrains = !E.isZero();
if (lHaveStrains)
{
// Invert the deformation gradient ==> Fi
Matrix Fi(nsd,nsd);
Fi.fill(F.ptr());
double J = Fi.inverse();
if (J == 0.0) return false;
// Scale with J=|F| since we are integrating over current configuration
detJW *= J;
if (eKm || iS)
{
// Push-forward the basis function gradients to current configuration
dNdx.multiply(dNdX,Fi); // dNdx = dNdX * F^-1
// Compute the small-deformation strain-displacement matrix B from dNdx
if (!this->formBmatrix(dNdx)) return false;
#if INT_DEBUG > 0
std::cout <<"PlasticityUL::dNdx ="<< dNdx;
std::cout <<"PlasticityUL::B ="<< Bmat;
#endif
}
}
else if (eKm || iS)
// Initial state, no deformation yet
if (!this->formBmatrix(dNdX)) return false;
// Evaluate the constitutive relation
SymmTensor sigma(nsd);
if (eKm || eKg || iS)
{
double U = 0.0;
if (!ptData->evaluate(Cmat,sigma,U,X,F,E,true,&prm))
return false;
}
if (eKm)
{
// Integrate the material stiffness matrix
CB.multiply(Cmat,Bmat).multiply(detJW); // CB = C*B*|J|*w
eKm->multiply(Bmat,CB,true,false,true); // EK += B^T * CB
}
if (eKg && lHaveStrains)
// Integrate the geometric stiffness matrix
this->formKG(*eKg,dNdx,sigma,detJW);
if (eM)
// Integrate the mass matrix
this->formMassMatrix(*eM,N,X,detJW);
if (iS && lHaveStrains)
{
// Integrate the internal forces
sigma *= -detJW;
if (!Bmat.multiply(sigma,*iS,true,true)) // ES -= B^T*sigma
return false;
}
if (eS)
// Integrate the load vector due to gravitation and other body forces
this->formBodyForce(*eS,N,X,detJW);
return this->getIntegralResult(elmInt);
}

72
Apps/FiniteDefElasticity/PlasticityUL.h
View File

@@ -1,72 +0,0 @@
// $Id$
//==============================================================================
//!
//! \file PlasticityUL.h
//!
//! \date Sep 21 2010
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integrand implementations for finite deformation plasticity problems.
//!
//==============================================================================
#ifndef _PLASTICITY_UL_H
#define _PLASTICITY_UL_H
#include "NonlinearElasticityUL.h"
class PlasticMaterial;
struct TimeDomain;
/*!
\brief Class representing the integrand of the nonlinear plasticity problem.
\details This class implements an Updated Lagrangian formulation. It inherits
most of the NonlinearElasticityUL methods, but reimplements \a evalInt using
an internal buffer to store history variables in every integration point
between the iterations.
*/
class PlasticityUL : public NonlinearElasticityUL
{
public:
//! \brief The default constructor invokes the parent class constructor.
//! \param[in] n Number of spatial dimensions
//! \param[in] lop Load option (0=on initial length, 1=on updated length)
PlasticityUL(unsigned short int n = 3, char lop = 0);
//! \brief The destructor frees the dynamically allocated integration buffers.
virtual ~PlasticityUL();
//! \brief Prints out problem definition to the given output stream.
virtual void print(std::ostream& os) const;
//! \brief Initializes the integrand for a new integration loop.
//! \param[in] prm Nonlinear solution algorithm parameters
//!
//! \details This method is invoked once before starting the numerical
//! integration over the entire spatial domain, and is used to reset the
//! internal integration point counter, and to update the history variables.
virtual void initIntegration(const TimeDomain& prm);
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] prm Nonlinear solution algorithm parameters
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& prm, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Returns a null pointer for solution norm evaluation.
virtual NormBase* getNormIntegrand(AnaSol* = 0) const { return 0; }
private:
mutable size_t iP; //!< Integration point counter
mutable std::vector<PlasticMaterial*> pBuf; //!< Integration point buffers
};
#endif

11
Apps/FiniteDefElasticity/SIMFiniteDefEl.C
View File

@@ -16,7 +16,6 @@
#include "LinearMaterial.h"
#include "NeoHookeMaterial.h"
#include "PlasticMaterial.h"
#include "PlasticityUL.h"
#include "NonlinearElasticityULMixed.h"
#include "NonlinearElasticityULMX.h"
#include "NeoHookeElasticity.h"
@@ -32,9 +31,6 @@ SIMFiniteDefEl2D::SIMFiniteDefEl2D (const std::vector<int>& options)
switch (form)
{
case SIM::PLASTICITY:
myProblem = new PlasticityUL(2);
break;
case SIM::MIXED_QnQn1:
nf[1] = 2; // continuous volumetric change and pressure fields
myProblem = new NonlinearElasticityULMixed(2);
@@ -109,7 +105,7 @@ bool SIMFiniteDefEl2D::parse (char* keyWord, std::istream& is)
RealArray pMAT;
while ((cline = strtok(NULL," ")))
pMAT.push_back(atof(cline));
mVec.push_back(new PlasticPrm(pMAT));
mVec.push_back(new PlasticMaterial(pMAT));
if (myPid == 0)
{
@@ -150,9 +146,6 @@ SIMFiniteDefEl3D::SIMFiniteDefEl3D (bool checkRHS,
switch (form)
{
case SIM::PLASTICITY:
myProblem = new PlasticityUL();
break;
case SIM::MIXED_QnQn1:
nf[1] = 2; // continuous volumetric change and pressure fields
myProblem = new NonlinearElasticityULMixed();
@@ -233,7 +226,7 @@ bool SIMFiniteDefEl3D::parse (char* keyWord, std::istream& is)
RealArray pMAT;
while ((cline = strtok(NULL," ")))
pMAT.push_back(atof(cline));
mVec.push_back(new PlasticPrm(pMAT));
mVec.push_back(new PlasticMaterial(pMAT));
if (myPid == 0)
{

3
Apps/FiniteDefElasticity/SIMFiniteDefEl.h
View File

@@ -28,8 +28,7 @@ namespace SIM
NEOHOOKE_IV = 5,
UPDATED_LAGRANGE = 6,
MIXED_QnPn1 = 7,
MIXED_QnQn1 = 8,
PLASTICITY = 9
MIXED_QnQn1 = 8
};
};

7
Apps/FiniteDefElasticity/main_NonLinEl.C
View File

@@ -51,7 +51,6 @@
\arg -2Dpstrain : Use two-parametric simulation driver (plane strain)
\arg -Tensor : Use tensorial total Lagrangian formulation (slow...)
\arg -UL : Use updated Lagrangian formulation with nonlinear material
\arg -PL : Use placticity solution driver
\arg -MX \a pord : Mixed formulation with internal discontinuous pressure
\arg -mixed : Mixed formulation with continuous pressure and volumetric change
\arg -lag : Use Lagrangian basis functions instead of splines/NURBS
@@ -146,8 +145,6 @@ int main (int argc, char** argv)
twoD = SIMLinEl2D::planeStrain = true;
else if (!strncmp(argv[i],"-2D",3))
twoD = true;
else if (!strcmp(argv[i],"-PL"))
form = SIM::PLASTICITY;
else if (!strcmp(argv[i],"-UL"))
{
if (form < SIM::UPDATED_LAGRANGE)
@@ -183,8 +180,8 @@ int main (int argc, char** argv)
<<" [-lag] [-spec] [-2D[pstrain]] [-UL] [-MX [<p>]|-mixed]"
<<" [-nGauss <n>]\n [-vtf <format>] [-nviz <nviz>]"
<<" [-nu <nu>] [-nv <nv>] [-nw <nw>]\n "
<<" [-saveInc <dtSave>] [-dumpInc <dtDump> [raw]]\n "
<<" [-ignore <p1> <p2> ...] [-fixDup] [-checkRHS] [-check]\n";
<<" [-saveInc <dtSave>] [-skip2nd] [-dumpInc <dtDump> [raw]]\n "
<<" [-ignore <p1> <p2> ...] [-fixDup] [-checkRHS] [-check]\n";
return 0;
}