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

Extensions for axisymmetric formulations

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1270 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo 2011-10-20 16:56:07 +00:00 committed by Knut Morten Okstad
parent 5530f29bef
commit fae5ff6b96
15 changed files with 261 additions and 147 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -1,12 +1,12 @@
subroutine cons2d (ipsw,iter,iwr,lfirst,mTYP,mVER,
& nsig,ndf,detF,F,pMAT,U,Sig,Cst,ierr)
& n,nsig,ndf,detF,F,pMAT,U,Sig,Cst,ierr)
C
C CONS2D: A 2D wrapper on the 3D material routines of FENRIS.
C
implicit none
C
integer ipsw, iter, iwr, lfirst, mTYP, mVER, nsig, ndf, ierr
real*8 detF, F(*), pMAT(*), U, Sig(4), Cst(3,3)
integer ipsw, iter, iwr, lfirst, mTYP, mVER, n, nsig, ndf, ierr
real*8 detF, F(*), pMAT(*), U, Sig(*), Cst(n,n)
real*8 F3D(3,3), S3D(6), C3D(6,6)
C
if (ndf .eq. 2) then
@ -27,10 +27,14 @@ C
return
end if
C
if (n .eq. 4) then
Cst = C3D(1:4,1:4)
else
Cst(1:2,1:2) = C3D(1:2,1:2)
Cst(1:2,3) = C3D(1:2,4)
Cst(3,1:2) = C3D(4,1:2)
Cst(3,3) = C3D(4,4)
end if
if (nsig .eq. 4) then
Sig(1:4) = S3D(1:4)
else

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

@ -1,13 +1,13 @@
subroutine plas2d (ipsw,iwr,iter,lfirst,pmat,nsig,ndf,detF,
subroutine plas2d (ipsw,iwr,iter,lfirst,pmat,n,nsig,ndf,detF,
& Fn1,Fn,be,Epp,Epl,Sig,Cst,ierr)
C
C PLAS2D: A 2D wrapper on the 3D plastic material routine of FENRIS.
C
implicit none
C
integer ipsw, iwr, iter, lfirst, nsig, ndf, ierr
integer ipsw, iwr, iter, lfirst, n, nsig, ndf, ierr
real*8 pmat(*), detF, Fn1(*), Fn(*)
real*8 be(*), Epp, Epl(*), Sig(4), Cst(3,3)
real*8 be(*), Epp, Epl(*), Sig(*), Cst(n,n)
real*8 F3D(3,3), F3P(3,3), S3D(4), C3D(6,6)
C
if (ndf .eq. 2) then
@ -34,10 +34,14 @@ C
return
end if
C
if (n .eq. 4) then
Cst = C3D(1:4,1:4)
else
Cst(1:2,1:2) = C3D(1:2,1:2)
Cst(1:2,3) = C3D(1:2,4)
Cst(3,1:2) = C3D(4,1:2)
Cst(3,3) = C3D(4,4)
end if
if (nsig .eq. 4) then
Sig(1:4) = S3D
else

2
Apps/FiniteDefElasticity/NeoHookeElasticity.C
View File

@ -70,7 +70,7 @@ bool NeoHookeElasticity::kinematics (const Matrix& dNdX,
Tensor& _F, SymmTensor& _E) const
{
// Form the deformation gradient, F
if (!this->NonlinearElasticity::kinematics(dNdX,_F,_E))
if (!this->NonlinearElasticity::kinematics(Vector(),dNdX,0.0,_F,_E))
return false;
else
J = _F.det();

14
Apps/FiniteDefElasticity/NeoHookeMaterial.C
View File

@ -19,8 +19,8 @@ extern "C" {
//! \brief Interface to 2D nonlinear material routines (FORTRAN-77 code).
void cons2d_(const int& ipsw, const int& iter, const int& iwr,
const int& lfirst, const int& mTYP, const int& mVER,
const int& nSig, const int& nDF, const double& detF,
const double* F, const double* pMAT,
const int& n, const int& nSig, const int& nDF,
const double& detF, const double* F, const double* pMAT,
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,
@ -112,17 +112,17 @@ bool NeoHookeMaterial::evaluate (Matrix& C, SymmTensor& sigma, double& U,
}
size_t ndim = sigma.dim();
size_t ncmp = ndim*(ndim+1)/2;
size_t ncmp = sigma.size();
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,sigma.size(),F.dim(),J,F.ptr(),pmat,
U,sigma.ptr(),C.ptr(),ierr);
cons2d_(INT_DEBUG,0,6,0,mTYP,mVER,C.cols(),sigma.size(),F.dim(),
J,F.ptr(),pmat,U,sigma.ptr(),C.ptr(),ierr);
else
cons3d_(INT_DEBUG,0,6,0,mTYP,mVER,0,sigma.size(),0,J,F.ptr(),pmat,&U,
U,sigma.ptr(),C.ptr(),ierr);
cons3d_(INT_DEBUG,0,6,0,mTYP,mVER,0,sigma.size(),0,
J,F.ptr(),pmat,&U,U,sigma.ptr(),C.ptr(),ierr);
#else
std::cerr <<" *** NeoHookeMaterial::evaluate: Not included."<< std::endl;
#endif

6
Apps/FiniteDefElasticity/NonlinearElasticity.C
View File

@ -84,7 +84,7 @@ bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt,
// Evaluate the kinematic quantities, F and E, at this point
Tensor F(nsd);
if (!this->kinematics(fe.dNdX,F,E))
if (!this->kinematics(fe.N,fe.dNdX,X.x,F,E))
return false;
// Evaluate current tangent at this point, that is
@ -148,7 +148,7 @@ bool NonlinearElasticity::evalInt (LocalIntegral*& elmInt,
if (eKg && haveStrains)
// Integrate the geometric stiffness matrix
this->formKG(*eKg,fe.dNdX,S,fe.detJxW);
this->formKG(*eKg,fe.N,fe.dNdX,X.x,S,fe.detJxW);
if (eM)
// Integrate the mass matrix
@ -232,7 +232,7 @@ bool NonlinearElasticity::formStressTensor (const Matrix& dNdX, const Vec3& X,
// Evaluate the kinematic quantities, F and E, at this point
Tensor F(nsd);
if (!this->kinematics(dNdX,F,E))
if (!this->kinematics(Vector(),dNdX,X.x,F,E))
return false;
// Evaluate the 2nd Piola-Kirchhoff stress tensor, S, at this point

56
Apps/FiniteDefElasticity/NonlinearElasticityFbar.C
View File

@ -76,9 +76,9 @@ bool NonlinearElasticityFbar::reducedInt (const FiniteElement& fe) const
VolPtData& ptData = myVolData[iP++];
// Evaluate the deformation gradient, F, at current configuration
Tensor F(nsd);
SymmTensor E(nsd);
if (!this->kinematics(fe.dNdX,F,E))
Tensor F(nDF);
SymmTensor E(nsd,axiSymmetry);
if (!this->kinematics(fe.N,fe.dNdX,0.0,F,E))
return false;
if (E.isZero(1.0e-16))
@ -91,8 +91,15 @@ bool NonlinearElasticityFbar::reducedInt (const FiniteElement& fe) const
{
// Invert the deformation gradient ==> Fi
Matrix Fi(nsd,nsd);
if (nsd == nDF)
Fi.fill(F.ptr());
else
for (unsigned short int i = 1; i <= nsd; i++)
for (unsigned short int j = 1; j <= nsd; j++)
Fi(i,j) = F(i,j);
ptData.J = Fi.inverse();
if (axiSymmetry) ptData.J *= F(3,3);
if (ptData.J == 0.0)
return false;
@ -274,9 +281,9 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
#endif
// Evaluate the deformation gradient, F, at current configuration
Tensor F(nsd);
SymmTensor E(nsd);
if (!this->kinematics(fe.dNdX,F,E))
Tensor F(nDF);
SymmTensor E(nsd,axiSymmetry);
if (!this->kinematics(fe.N,fe.dNdX,X.x,F,E))
return false;
double J, Jbar = 0.0;
@ -285,8 +292,15 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
{
// Invert the deformation gradient ==> Fi
Matrix Fi(nsd,nsd);
if (nDF == nsd)
Fi.fill(F.ptr());
else
for (unsigned short int i = 1; i <= nsd; i++)
for (unsigned short int j = 1; j <= nsd; j++)
Fi(i,j) = F(i,j);
J = Fi.inverse();
if (axiSymmetry) J *= F(3,3);
if (J == 0.0) return false;
if (eKm || iS)
@ -343,7 +357,7 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
}
// Compute modified deformation gradient, Fbar
F *= pow(fabs(Jbar/J),1.0/(double)nsd);
F *= pow(fabs(Jbar/J),1.0/(double)nDF);
#ifdef INT_DEBUG
std::cout <<"NonlinearElasticityFbar::Jbar = "<< Jbar
<<"\nNonlinearElasticityFbar::dMdx ="<< dMdx
@ -351,13 +365,16 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
#endif
// Evaluate the constitutive relation (Jbar is dummy here)
SymmTensor Sig(nsd);
SymmTensor Sig(nsd,axiSymmetry);
if (!material->evaluate(Cmat,Sig,Jbar,X,F,E,lHaveStrains,&prm))
return false;
// Axi-symmetric integration point volume; 2*pi*r*|J|*w
double detJW = (axiSymmetry ? 2.0*M_PI*X.x : 1.0)*fe.detJxW*J;
// Multiply tangent moduli and stresses by integration point volume
Cmat *= fe.detJxW*J;
Sig *= fe.detJxW*J;
Cmat *= detJW;
Sig *= detJW;
if (iS && lHaveStrains)
{
@ -420,11 +437,11 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
if (eM)
// Integrate the mass matrix
this->formMassMatrix(*eM,fe.N,X,J*fe.detJxW);
this->formMassMatrix(*eM,fe.N,X,detJW);
if (eS)
// Integrate the load vector due to gravitation and other body forces
this->formBodyForce(*eS,fe.N,X,J*fe.detJxW);
this->formBodyForce(*eS,fe.N,X,detJW);
return this->getIntegralResult(elmInt);
}
@ -452,9 +469,9 @@ bool ElasticityNormFbar::evalInt (LocalIntegral*& elmInt,
NonlinearElasticityFbar& p = static_cast<NonlinearElasticityFbar&>(myProblem);
// Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
Tensor F(nsd);
SymmTensor E(nsd);
if (!p.kinematics(fe.dNdX,F,E))
Tensor F(p.nDF);
SymmTensor E(p.nDF);
if (!p.kinematics(fe.N,fe.dNdX,X.x,F,E))
return false;
double Jbar = 0.0;
@ -479,16 +496,19 @@ bool ElasticityNormFbar::evalInt (LocalIntegral*& elmInt,
// Compute modified deformation gradient, Fbar
Tensor Fbar(F);
Fbar *= pow(fabs(Jbar/F.det()),1.0/(double)nsd);
Fbar *= pow(fabs(Jbar/F.det()),1.0/(double)p.nDF);
// Compute the strain energy density, U(E) = Int_E (S:Eps) dEps
// and the Cauchy stress tensor, sigma
double U = 0.0;
SymmTensor sigma(nsd,p.material->isPlaneStrain());
SymmTensor sigma(nsd, p.isAxiSymmetric() || p.material->isPlaneStrain());
if (!p.material->evaluate(p.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),
sigma,U,F.det(),fe.detJxW);
sigma,U,F.det(),detJW);
}

41
Apps/FiniteDefElasticity/NonlinearElasticityTL.C
View File

@ -92,7 +92,7 @@ bool NonlinearElasticityTL::evalInt (LocalIntegral*& elmInt,
// and compute the nonlinear strain-displacement matrix, B, from dNdX and F
Tensor F(nsd);
SymmTensor E(nsd), S(nsd);
if (!this->kinematics(fe.dNdX,F,E))
if (!this->kinematics(fe.N,fe.dNdX,X.x,F,E))
return false;
// Evaluate the constitutive relation
@ -104,32 +104,35 @@ bool NonlinearElasticityTL::evalInt (LocalIntegral*& elmInt,
return false;
}
// Axi-symmetric integration point volume; 2*pi*r*|J|*w
const double detJW = axiSymmetry ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW;
if (eKm)
{
// Integrate the material stiffness matrix
CB.multiply(Cmat,Bmat).multiply(fe.detJxW); // CB = C*B*|J|*w
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,fe.dNdX,S,fe.detJxW);
this->formKG(*eKg,fe.N,fe.dNdX,X.x,S,detJW);
if (eM)
// Integrate the mass matrix
this->formMassMatrix(*eM,fe.N,X,fe.detJxW);
this->formMassMatrix(*eM,fe.N,X,detJW);
if (iS && lHaveStrains)
{
// Integrate the internal forces
S *= -fe.detJxW;
S *= -detJW;
if (!Bmat.multiply(S,*iS,true,true)) // ES -= B^T*S
return false;
}
if (eS)
// Integrate the load vector due to gravitation and other body forces
this->formBodyForce(*eS,fe.N,X,fe.detJxW);
this->formBodyForce(*eS,fe.N,X,detJW);
return this->getIntegralResult(elmInt);
}
@ -165,7 +168,7 @@ bool NonlinearElasticityTL::evalBou (LocalIntegral*& elmInt,
// Compute the deformation gradient, F
Tensor F(nsd);
SymmTensor dummy(0);
if (!this->kinematics(fe.dNdX,F,dummy)) return false;
if (!this->kinematics(fe.N,fe.dNdX,X.x,F,dummy)) return false;
// Compute its inverse and determinant, J
double J = F.inverse();
@ -180,24 +183,29 @@ bool NonlinearElasticityTL::evalBou (LocalIntegral*& elmInt,
T[i-1] += t[j-1]*F(j,i);
}
// Axi-symmetric integration point volume; 2*pi*r*|J|*w
const double detJW = axiSymmetry ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW;
// Integrate the force vector
Vector& ES = *eS;
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)*fe.detJxW;
ES(nsd*(a-1)+i) += T[i-1]*fe.N(a)*detJW;
return this->getIntegralResult(elmInt);
}
bool NonlinearElasticityTL::kinematics (const Matrix& dNdX,
Tensor& F, SymmTensor& E) const
bool NonlinearElasticityTL::kinematics (const Vector& N, const Matrix& dNdX,
double r, Tensor& F,
SymmTensor& E) const
{
if (!eV || eV->empty())
{
// Initial state, unit deformation gradient and linear B-matrix
F = 1.0;
return this->Elasticity::kinematics(dNdX,F,E);
E.zero();
return axiSymmetry ? this->formBmatrix(N,dNdX,r) : this->formBmatrix(dNdX);
}
const size_t nenod = dNdX.rows();
@ -216,12 +224,17 @@ bool NonlinearElasticityTL::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*dNd
F = dUdX;
else
if (!dUdX.multiplyMat(*eV,dNdX)) // dUdX = Grad{u} = eV*dNd
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

5
Apps/FiniteDefElasticity/NonlinearElasticityTL.h
View File

@ -62,7 +62,9 @@ public:
protected:
//! \brief Calculates some kinematic quantities at current point.
//! \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[in] F Deformation gradient at current point
//! \param[out] E Green-Lagrange strain tensor at current point
//!
@ -70,7 +72,8 @@ protected:
//! strain-displacement matrix \b B are established. The latter matrix
//! is stored in the mutable class member \a Bmat of the parent class.
//! The B-matrix is formed only when the variable \a formB is true.
virtual bool kinematics(const Matrix& dNdX, Tensor& F, SymmTensor& E) const;
virtual bool kinematics(const Vector& N, const Matrix& dNdX, double r,
Tensor& F, SymmTensor& E) const;
protected:
bool formB; //!< Flag determining whether we need to form the B-matrix

100
Apps/FiniteDefElasticity/NonlinearElasticityUL.C
View File

@ -20,9 +20,15 @@
#include "Tensor.h"
#include "Vec3Oper.h"
#ifndef epsR
//! \brief Zero tolerance for the radial coordinate.
#define epsR 1.0e-16
#endif
NonlinearElasticityUL::NonlinearElasticityUL (unsigned short int n, char lop)
: Elasticity(n), loadOp(lop), plam(-999.9)
NonlinearElasticityUL::NonlinearElasticityUL (unsigned short int n,
bool axS, char lop)
: Elasticity(n,axS), loadOp(lop), plam(-999.9)
{
// Only the current solution is needed
primsol.resize(1);
@ -108,29 +114,47 @@ bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt,
const Vec3& X) const
{
// Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
Tensor F(nsd);
SymmTensor E(nsd);
if (!this->kinematics(fe.dNdX,F,E))
Tensor F(nDF);
SymmTensor E(nsd,axiSymmetry);
if (!this->kinematics(fe.N,fe.dNdX,X.x,F,E))
return false;
// 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 : 0.0;
bool lHaveStrains = !E.isZero(1.0e-16);
if (lHaveStrains)
{
// Invert the deformation gradient ==> Fi
Matrix Fi(nsd,nsd);
if (nDF == nsd)
Fi.fill(F.ptr());
else
for (unsigned short int i = 1; i <= nsd; i++)
for (unsigned short int j = 1; j <= nsd; j++)
Fi(i,j) = F(i,j);
double J = Fi.inverse();
if (axiSymmetry) J *= F(3,3);
if (J == 0.0) return false;
// Scale with J=|F| since we are integrating over current configuration
const_cast<double&>(fe.detJxW) *= J;
detJW *= J;
if (eKm || iS)
{
// Push-forward the basis function gradients to current configuration
dNdx.multiply(fe.dNdX,Fi); // dNdx = dNdX * F^-1
// Compute the small-deformation strain-displacement matrix B from dNdx
if (!this->formBmatrix(dNdx)) return false;
if (axiSymmetry)
{
r += eV->dot(fe.N,0,nsd);
this->formBmatrix(fe.N,dNdx,r);
}
else
this->formBmatrix(dNdx);
#ifdef INT_DEBUG
std::cout <<"NonlinearElasticityUL::dNdx ="<< dNdx;
std::cout <<"NonlinearElasticityUL::B ="<< Bmat;
@ -138,11 +162,16 @@ bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt,
}
}
else if (eKm || iS)
{
// Initial state, no deformation yet
if (!this->formBmatrix(fe.dNdX)) return false;
if (axiSymmetry)
this->formBmatrix(fe.N,fe.dNdX,r);
else
this->formBmatrix(fe.dNdX);
}
// Evaluate the constitutive relation
SymmTensor sigma(nsd);
SymmTensor sigma(nsd,axiSymmetry);
if (eKm || eKg || iS)
{
double U = 0.0;
@ -153,29 +182,29 @@ bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt,
if (eKm)
{
// Integrate the material stiffness matrix
CB.multiply(Cmat,Bmat).multiply(fe.detJxW); // CB = C*B*|J|*w
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,fe.detJxW);
this->formKG(*eKg,fe.N,dNdx,r,sigma,detJW);
if (eM)
// Integrate the mass matrix
this->formMassMatrix(*eM,fe.N,X,fe.detJxW);
this->formMassMatrix(*eM,fe.N,X,detJW);
if (iS && lHaveStrains)
{
// Integrate the internal forces
sigma *= -fe.detJxW;
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,fe.N,X,fe.detJxW);
this->formBodyForce(*eS,fe.N,X,detJW);
return this->getIntegralResult(elmInt);
}
@ -204,12 +233,15 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt,
// Store the traction value for vizualization
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;
unsigned short int i, j;
if (loadOp == 1)
{
// Compute the deformation gradient, F
Tensor F(nsd);
if (!this->formDefGradient(fe.dNdX,F)) return false;
Tensor F(nDF);
if (!this->formDefGradient(fe.N,fe.dNdX,X.x,F)) return false;
// Check for with-rotated pressure load
if (tracFld->isNormalPressure())
@ -228,29 +260,31 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt,
}
else
// Scale with J=|F| since we are integrating over current configuration
const_cast<double&>(fe.detJxW) *= F.det();
detJW *= F.det();
}
// Integrate the force vector
Vector& ES = *eS;
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)*fe.detJxW;
ES(nsd*(a-1)+i) += T[i-1]*fe.N(a)*detJW;
return this->getIntegralResult(elmInt);
}
bool NonlinearElasticityUL::formDefGradient (const Matrix& dNdX,
bool NonlinearElasticityUL::formDefGradient (const Vector& N,
const Matrix& dNdX, double r,
Tensor& F) const
{
static SymmTensor dummy(0);
return this->kinematics(dNdX,F,dummy);
return this->kinematics(N,dNdX,r,F,dummy);
}
bool NonlinearElasticityUL::kinematics (const Matrix& dNdX,
Tensor& F, SymmTensor& E) const
bool NonlinearElasticityUL::kinematics (const Vector& N, const Matrix& dNdX,
double r, Tensor& F,
SymmTensor& E) const
{
if (!eV || eV->empty())
{
@ -268,7 +302,6 @@ 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).
@ -291,6 +324,7 @@ bool NonlinearElasticityUL::kinematics (const Matrix& dNdX,
// Now form the Green-Lagrange strain tensor.
// Note that for the shear terms (i/=j) we actually compute 2*E_ij
// to be consistent with the engineering strain style constitutive matrix.
// TODO: How is this for axisymmetric problems?
for (i = 1; i <= E.dim(); i++)
for (j = 1; j <= i; j++)
{
@ -302,6 +336,8 @@ bool NonlinearElasticityUL::kinematics (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;
#ifdef INT_DEBUG
std::cout <<"NonlinearElasticityUL::eV ="<< *eV;
@ -333,20 +369,23 @@ bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt,
NonlinearElasticityUL& ulp = static_cast<NonlinearElasticityUL&>(myProblem);
// Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
Tensor F(fe.dNdX.cols());
SymmTensor E(fe.dNdX.cols());
if (!ulp.kinematics(fe.dNdX,F,E))
Tensor F(ulp.nDF);
SymmTensor E(ulp.nDF);
if (!ulp.kinematics(fe.N,fe.dNdX,X.x,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;
SymmTensor sigma(E.dim(),ulp.material->isPlaneStrain());
SymmTensor sigma(E.dim(),ulp.isAxiSymmetric()||ulp.material->isPlaneStrain());
if (!ulp.material->evaluate(ulp.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(),fe.detJxW);
return evalInt(getElmNormBuffer(elmInt,6),sigma,U,F.det(),detJW);
}
@ -403,6 +442,9 @@ bool ElasticityNormUL::evalBou (LocalIntegral*& elmInt,
up.push_back(u);
}
getElmNormBuffer(elmInt)[1] += Ux[iP++]*fe.detJxW;
// 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;
return true;
}

14
Apps/FiniteDefElasticity/NonlinearElasticityUL.h
View File

@ -31,8 +31,10 @@ class NonlinearElasticityUL : public Elasticity
public:
//! \brief The default constructor invokes the parent class constructor.
//! \param[in] n Number of spatial dimensions
//! \param[in] axS \e If \e true, and axisymmetric 3D formulation is assumed
//! \param[in] lop Load option (0=on initial length, 1=on updated length)
NonlinearElasticityUL(unsigned short int n = 3, char lop = 0);
NonlinearElasticityUL(unsigned short int n = 3, bool axS = false,
char lop = 0);
//! \brief Empty destructor.
virtual ~NonlinearElasticityUL() {}
@ -78,16 +80,22 @@ public:
virtual NormBase* getNormIntegrand(AnaSol* = 0) const;
//! \brief Calculates some kinematic quantities at current point.
//! \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] E Green-Lagrange strain tensor at current point
virtual bool kinematics(const Matrix& dNdX, Tensor& F, SymmTensor& E) const;
virtual bool kinematics(const Vector& N, const Matrix& dNdX, double r,
Tensor& F, SymmTensor& E) const;
protected:
//! \brief Calculates the deformation gradient at current point.
//! \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 Matrix& dNdX, Tensor& F) const;
bool formDefGradient(const Vector& N, const Matrix& dNdX, double r,
Tensor& F) const;
protected:
mutable Matrix dNdx; //!< Basis function gradients in current configuration

18
Apps/FiniteDefElasticity/NonlinearElasticityULMX.C
View File

@ -155,22 +155,26 @@ bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt,
// Evaluate the deformation gradient, Fp, at previous configuration
const_cast<NonlinearElasticityULMX*>(this)->eV = &eVs[1];
if (!this->formDefGradient(fe.dNdX,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(fe.dNdX,ptData.F))
if (!this->formDefGradient(fe.N,fe.dNdX,X.x,ptData.F))
return false;
if (nsd == 2) // In 2D we always assume plane strain so set F(3,3)=1
if (nDF == 2) // In 2D we always assume plane strain so set F(3,3)=1
ptData.F(3,3) = ptData.Fp(3,3) = 1.0;
// Invert the deformation gradient ==> Fi
Matrix Fi(nsd,nsd);
if (nsd == 3)
Fi.fill(ptData.F.ptr());
else
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;
@ -191,15 +195,15 @@ bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt,
if (Hh)
// Integrate the Hh-matrix
Hh->outer_product(ptData.Phi,ptData.Phi*fe.detJxW,true);
Hh->outer_product(ptData.Phi,ptData.Phi*ptData.detJW,true);
if (eM)
// Integrate the mass matrix
this->formMassMatrix(*eM,fe.N,X,J*fe.detJxW);
this->formMassMatrix(*eM,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*fe.detJxW);
this->formBodyForce(*eS,fe.N,X,J*ptData.detJW);
return this->getIntegralResult(elmInt);
}
@ -409,7 +413,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
#endif
if (eKg)
// Integrate the geometric stiffness matrix
this->formKG(*eKg,pt.dNdx,Sigma,dFmix*pt.detJW);
this->formKG(*eKg,Vector(),pt.dNdx,X.x,Sigma,dFmix*pt.detJW);
}
if (eKm)

56
Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C
View File

@ -312,9 +312,9 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
#endif
// Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
Tensor F(nsd);
SymmTensor E(nsd);
if (!this->kinematics(fe.dN1dX,F,E))
Tensor F(nDF);
SymmTensor E(nsd,axiSymmetry);
if (!this->kinematics(fe.N1,fe.dN1dX,X.x,F,E))
return false;
bool lHaveStrains = !E.isZero(1.0e-16);
@ -329,28 +329,37 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
<<", 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;
// Compute the mixed model deformation gradient, F_bar
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
if (nDF == 2) // In 2D we always assume plane strain so set F(3,3)=1
Fbar(3,3) = r1;
else if (nsd != 3)
return false;
// Invert the deformation gradient ==> Fi
Matrix Fi(nsd,nsd);
if (nDF == nsd)
Fi.fill(F.ptr());
else
for (unsigned short int i = 1; i <= nsd; i++)
for (unsigned short int j = 1; j <= nsd; j++)
Fi(i,j) = F(i,j);
double J = Fi.inverse();
if (axiSymmetry) J *= F(3,3);
if (J == 0.0) return false;
// Push-forward the basis function gradients to current configuration
dNdx.multiply(fe.dN1dX,Fi); // dNdx = dNdX * F^-1
// Compute the mixed integration point volume
double dVol = Theta*fe.detJxW;
double dVup = J*fe.detJxW;
double dVol = Theta*detJW;
double dVup = J*detJW;
#if INT_DEBUG > 0
std::cout <<"NonlinearElasticityULMixed::dNdX ="<< fe.dN1dX;
@ -363,14 +372,14 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
if (eM)
// Integrate the mass matrix
this->formMassMatrix(*eM,fe.N1,X,J*fe.detJxW);
this->formMassMatrix(*eM,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*fe.detJxW);
this->formBodyForce(*eS,fe.N1,X,J*detJW);
// Evaluate the constitutive relation
SymmTensor Sig(3), Sigma(nsd);
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;
@ -410,7 +419,7 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
Sigma = Sig + (Mpres-Bpres);
// Integrate the geometric stiffness matrix
this->formKG(*eKg,dNdx,Sigma,dVol);
this->formKG(*eKg,fe.N1,dNdx,r,Sigma,dVol);
}
#ifdef USE_FTNMAT
@ -449,13 +458,15 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
}
myMats->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*(-fe.detJxW),true); // -= N2*N2^T*|J|
myMats->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 (!this->Elasticity::formBmatrix(dNdx))
return false;
if (axiSymmetry)
this->formBmatrix(fe.N1,dNdx,r);
else
this->formBmatrix(dNdx);
// Integrate the internal forces
Sigma *= -dVol;
@ -463,8 +474,8 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
return false;
// Integrate the volumetric change and pressure forces
myMats->b[Rt].add(fe.N2,(Press - Bpres)*fe.detJxW); // += N2*(p-pBar)*|J|
myMats->b[Rp].add(fe.N2,(Theta - J)*fe.detJxW); // += N2*(Theta-J)*|J|
myMats->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|
#if INT_DEBUG > 4
std::cout <<"NonlinearElasticityULMixed::Sigma*dVol =\n"<< Sigma;
@ -541,9 +552,9 @@ bool ElasticityNormULMixed::evalIntMx (LocalIntegral*& elmInt,
ulp = static_cast<NonlinearElasticityULMixed*>(&myProblem);
// Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
Tensor F(fe.dN1dX.cols());
SymmTensor E(F.dim());
if (!ulp->kinematics(fe.dN1dX,F,E))
Tensor F(ulp->nDF);
SymmTensor E(ulp->nDF);
if (!ulp->kinematics(fe.N1,fe.dN1dX,X.x,F,E))
return false;
// Evaluate the volumetric change field
@ -564,8 +575,11 @@ bool ElasticityNormULMixed::evalIntMx (LocalIntegral*& elmInt,
if (!ulp->material->evaluate(ulp->Cmat,Sig,U,X,ulp->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,fe.detJxW);
return evalInt(getElmNormBuffer(elmInt,6),Sig,U,Theta,detJW);
}

15
Apps/FiniteDefElasticity/PlasticMaterial.C
View File

@ -19,8 +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 double* pMAT,
const int& nSig, const int& nDF, const double& detF,
const double* Fn1, const double* Fn,
const int& n, const int& nSig, const int& nDF,
const double& detF, 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).
@ -143,6 +143,8 @@ bool PlasticMaterial::evaluate (Matrix& C, SymmTensor& sigma, double& U,
const SymmTensor& eps, char iop,
const TimeDomain* prm, const Tensor* Fpf) const
{
C.resize(sigma.size(),sigma.size());
if (iAmIntegrating)
{
if (!prm)
@ -262,10 +264,6 @@ bool PlasticMaterial::PlasticPoint::evaluate (Matrix& C,
// 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
@ -275,8 +273,9 @@ bool PlasticMaterial::PlasticPoint::evaluate (Matrix& C,
for (int i = 0; i < 10; i++) std::cout <<" "<< HVc[i];
std::cout << std::endl;
#endif
if (ndim == 2)
plas2d_(INT_DEBUG,6,prm.it,prm.first,&pMAT.front(),sigma.size(),F.dim(),J,
if (sigma.dim() == 2)
plas2d_(INT_DEBUG,6,prm.it,prm.first,&pMAT.front(),
C.cols(),sigma.size(),F.dim(),J,
F.ptr(),Fp.ptr(),HVc,HVc[6],HVc+7,sigma.ptr(),C.ptr(),ierr);
else
plas3d_(INT_DEBUG,6,prm.it,prm.first,6,6,&pMAT.front(),J,

2
Apps/FiniteDefElasticity/SIMFiniteDefEl.C
View File

@ -45,7 +45,7 @@ SIMFiniteDefEl2D::SIMFiniteDefEl2D (const std::vector<int>& options)
myProblem = new NonlinearElasticityULMX(2,pOrd);
break;
case SIM::UPDATED_LAGRANGE:
myProblem = new NonlinearElasticityUL(2);
myProblem = new NonlinearElasticityUL(2,axiSymmetry);
break;
case SIM::TOTAL_LAGRANGE:
myProblem = new NonlinearElasticityTL(2);

7
Apps/FiniteDefElasticity/main_NonLinEl.C
View File

@ -61,6 +61,7 @@ extern std::vector<int> mixedDbgEl; //!< List of elements for additional output
\arg -fixDup : Resolve co-located nodes by merging them into a single node
\arg -2D : Use two-parametric simulation driver (plane stress)
\arg -2Dpstrain : Use two-parametric simulation driver (plane strain)
\arg -2Daxisymm : Use two-parametric simulation driver (axi-symmetric solid)
\arg -UL : Use updated Lagrangian formulation with nonlinear material
\arg -MX \a pord : Mixed formulation with internal discontinuous pressure
\arg -mixed : Mixed formulation with continuous pressure and volumetric change
@ -163,8 +164,10 @@ int main (int argc, char** argv)
checkRHS = true;
else if (!strcmp(argv[i],"-fixDup"))
fixDup = true;
else if (!strcmp(argv[i],"-2Dpstrain"))
else if (!strncmp(argv[i],"-2Dpstra",8))
twoD = SIMLinEl2D::planeStrain = true;
else if (!strncmp(argv[i],"-2Daxi",6))
twoD = SIMLinEl2D::axiSymmetry = true;
else if (!strncmp(argv[i],"-2D",3))
twoD = true;
else if (!strcmp(argv[i],"-UL"))
@ -206,7 +209,7 @@ int main (int argc, char** argv)
{
std::cout <<"usage: "<< argv[0]
<<" <inputfile> [-dense|-spr|-superlu[<nt>]|-samg|-petsc]\n "
<<" [-lag] [-spec] [-nGauss <n>] [-2D[pstrain]] [-UL|-MX [<p>]"
<<" [-lag|-spec] [-nGauss <n>] [-2D[pstrain|axis]] [-UL|-MX [<p>]"
<<"|-[M|m]ixed|-Fbar <nvp>]\n [-vtf <format> [-nviz <nviz>]"
<<" [-nu <nu>] [-nv <nv>] [-nw <nw>]] [-hdf5]\n "
<<" [-saveInc <dtSave>] [-skip2nd] [-dumpInc <dtDump> [raw]]"