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

Finite deformation plasticity solver added

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@864 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo 2011-03-15 17:21:47 +00:00 committed by Knut Morten Okstad
parent 324c863956
commit 9a34254c7f
17 changed files with 612 additions and 106 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -56,14 +56,8 @@ C Bmod - Material property (Bulk modulus)
C Emod - Material property (Youngs modulus)
C Pratio - Material property (Poissons ratio)
C Smod - Material property (Shear modulus)
C GLe(6) - Green-Lagrangse strains
C PK2(6) - 2nd Piola-Kirchhoff stresses (material stresses)
C Cmt(6,6) - Material constitutive tensor
C incomp - logical variable
C = .true. if incompressible material ): Pratio = 0.5
C istrt - Start state:
C = 0 : Elastic
C = 1 : Last solution
C
C PRINT SWITCH:
C ipsw = 0 Gives no print
@ -93,7 +87,7 @@ C
real*8 detF, Engy, F(3,*), pMAT(*), Sig(*), Cst(*), HV(*)
C
real*8 Bmod, Emod, Pratio, Smod
logical incomp, state
logical incomp
C
include 'include/feninc/const.h'
C
@ -193,7 +187,7 @@ C
C Finite stretch plasticity model
C
call plas3d(ipsw, iwr, iter, lfirst, ntm, pMAT, detF, F(1,1),
& F(1,4), HV(1), HV(7), HV(8), Sig, Cst, state, ierr)
& F(1,4), HV(1), HV(7), HV(8), Sig, Cst, ierr)
if (ierr .lt. 0) go to 7000
C
Engy = zero

37
Apps/FiniteDefElasticity/Material/plas2d.f Normal file
View File

@ -0,0 +1,37 @@
subroutine plas2d (ipsw,iwr,iter,lfirst,ntm,pmat,detF,
& 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(*)
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)
C
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)
Sig(1:2) = S3D(1:2)
Sig(3) = S3D(4)
C
return
end

24
Apps/FiniteDefElasticity/Material/plas3d.f
View File

@ -1,5 +1,5 @@
subroutine plas3d (ipsw, iwr, iter, lfirst, ntm, pMAT, detF,
& Fn1, Fn, be, Epp, Epl, Sig, Cst, state, ierr)
& Fn1, Fn, be, Epp, Epl, Sig, Cst, ierr)
C
C ----------------------------------------------------------------------
C
@ -68,7 +68,6 @@ C Epp - Cumulative plastic strain at t_n+1
C Epl(*) - Plastic Strain for Hardening at t_n+1
C Sig(6) - Cauchy stresses (spatial stresses)
C Cst(6,6) - Spatial constitutive tensor
C state -
C ierr - Error flag
C = 1 : Illegal number of element nodes
C
@ -93,6 +92,7 @@ C yield - yield function value
C yfunct - yield function
C flg = .false. only yield function
C flg = .true. yield function and derivatives
C state -
C
C PRINT SWITCH:
C ipsw = 0 Gives no print
@ -121,7 +121,6 @@ C
integer ipsw, iwr, istrt, iter, lfirst, ntm, ierr
real*8 pMAT(*), detF, Fn1(3,3), Fn(3,3),
& Epp, be(*), Epl(*), Sig(*), Cst(6,6)
logical state
C
integer i, j, a, b, rot, it, itmax, p1(6), p2(6), ierrl, iop
real*8 Bmod, Smod, F(3,3), Finv(3,3), detFr, detr
@ -135,7 +134,7 @@ C
real*8 dtde(3,3), ff(6,6), Cstm(6,6), Eppn
real*8 nn_t(3,3), nn(3), aa(3), bb(3)
real*8 tolb, tolc, d_el, xx, f112, f123, yield, yfunct
logical conv
logical conv, state
C
include 'include/feninc/const.h'
C
@ -193,21 +192,14 @@ C
C
C Check state for iterations
C
state = .true.
dyld = zero
if (iter .eq. 0) then ! First iteration in step
C
if(istrt.eq.0) then ! Elastic state requested
if (istrt .eq. 0) then ! Elastic state requested
state = .false.
dyld = zero
else ! Last state requested
state = .true.
dyld = 1.0d-08*Y0
dyld = 1.0d-08*Y0
endif
C
else ! Not first iteration in step
C
state = .true.
dyld = zero
C
endif
C
C KINEMATIC COMPUTATIONS (Incremental Lagrangian formulation)
@ -538,8 +530,6 @@ C
else
C
C ELASTIC STEP ( only tangent computation )
C
state = .false. ! Indicate elastic on return
C
dtde(1,1) = one2 * Bmod + two3 * Smod
dtde(1,2) = one2 * Bmod - one3 * Smod

11
Apps/FiniteDefElasticity/NeoHookeMaterial.C
View File

@ -71,7 +71,7 @@ void NeoHookeMaterial::print (std::ostream& os) const
bool NeoHookeMaterial::evaluate (Matrix& C, SymmTensor& sigma, double& U,
const Vec3& X, const Tensor& F,
const Vec3&, const Tensor& F,
const SymmTensor&, char iop,
const TimeDomain*) const
{
@ -86,16 +86,17 @@ bool NeoHookeMaterial::evaluate (Matrix& C, SymmTensor& sigma, double& U,
size_t ndim = F.dim();
size_t ncmp = ndim*(ndim+1)/2;
C.resize(ncmp,ncmp);
int ipsw = INT_DEBUG > 1 ? 9 : 0;
int ierr = 0;
int ierr = -99;
#ifdef USE_FTNMAT
// Invoke the FORTRAN routine for Neo-Hookean hyperelastic material models.
if (ndim == 2)
cons2d_(ipsw,0,6,0,mTYP,mVER,0,3,J,F.ptr(),pmat,
cons2d_(INT_DEBUG,0,6,0,mTYP,mVER,0,3,J,F.ptr(),pmat,
&U,U,sigma.ptr(),C.ptr(),ierr);
else
cons3d_(ipsw,0,6,0,mTYP,mVER,0,6,J,F.ptr(),pmat,
cons3d_(INT_DEBUG,0,6,0,mTYP,mVER,0,6,J,F.ptr(),pmat,
&U,U,sigma.ptr(),C.ptr(),ierr);
#else
std::cerr <<" *** NeoHookeMaterial::evaluate: Not included."<< std::endl;
#endif
if (iop == 2)

7
Apps/FiniteDefElasticity/NonlinearElasticityUL.h
View File

@ -16,17 +16,14 @@
#include "Elasticity.h"
struct TimeDomain;
/*!
\brief Class representing the integrand of the nonlinear elasticity problem.
\details This class implements an Updated Lagrangian formulation. It inherits
most of the Elasticity methods, but reimplements the \a kinematics method
for calculating the deformation gradient and the associated Green-Lagrange
strain tensor, and the method \a constitutive for calculating the tangent
constitutive matrix and the associated stress tensor. The \a evalInt and
\a evalBou methods are also reimplemented to account for the updated geometry.
strain tensor. The \a evalInt and \a evalBou methods are also reimplemented
to account for the updated geometry.
*/
class NonlinearElasticityUL : public Elasticity

102
Apps/FiniteDefElasticity/PlasticMaterial.C Normal file
View File

@ -0,0 +1,102 @@
// $Id$
//==============================================================================
//!
//! \file PlasticMaterial.C
//!
//! \date Mar 16 2011
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Plasticity material models.
//!
//==============================================================================
#include "PlasticMaterial.h"
#include "TimeDomain.h"
#ifdef USE_FTNMAT
extern "C" {
//! \brief Interface to 2D 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* be, const double& Epp, const double* Epl,
const double* Sig, double* Cst, int& ierr);
//! \brief Interface to 3D plastic material routines (FORTRAN-77 code).
void plas3d_(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* be, const double& Epp, const double* Epl,
const double* Sig, double* Cst, int& ierr);
}
#ifndef INT_DEBUG
#define INT_DEBUG 0
#endif
#endif
PlasticMaterial::PlasticMaterial (const PlasticPrm* prm, unsigned short int n)
: pMAT(prm->pMAT), Fp(n)
{
memset(HVc,0,sizeof(HVc));
memset(HVp,0,sizeof(HVc));
}
void PlasticPrm::print (std::ostream& os) const
{
std::cout <<"PlasticMaterial: pMAT =";
for (size_t i = 0; i < pMAT.size(); i++)
std::cout <<" "<< pMAT[i];
if (rho > 0.0)
std::cout <<"\nPlasticMaterial: rho ="<< rho;
std::cout << std::endl;
}
bool PlasticMaterial::evaluate (Matrix& C, SymmTensor& sigma, double&,
const Vec3&, const Tensor& F, const SymmTensor&,
char, const TimeDomain* prm) const
{
double J = F.det();
if (J == 0.0)
{
std::cerr <<" *** PlasticMaterial::evaluate: "
<<" 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();
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));
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);
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);
#if INT_DEBUG > 0
std::cout <<"PlasticMaterial::sigma =\n"<< sigma;
std::cout <<"PlasticMaterial::C ="<< C;
#endif
#else
std::cerr <<" *** PlasticMaterial::evaluate: Not included."<< std::endl;
#endif
return ierr == 0;
}
void PlasticMaterial::updateHistoryVars ()
{
memcpy(HVp,HVc,sizeof(HVc));
}

87
Apps/FiniteDefElasticity/PlasticMaterial.h Normal file
View File

@ -0,0 +1,87 @@
// $Id$
//==============================================================================
//!
//! \file PlasticMaterial.h
//!
//! \date Mar 16 2011
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Plasticity material models.
//!
//==============================================================================
#ifndef _PLASTIC_MATERIAL_H
#define _PLASTIC_MATERIAL_H
#include "MaterialBase.h"
#include "Tensor.h"
/*!
\brief Class containing parameters of a plasticity material model.
*/
class PlasticPrm : public Material
{
public:
//! \brief Constructor initializing the material parameters.
PlasticPrm(const RealArray& p, double dens = 0.0) : pMAT(p), rho(dens) {}
//! \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 rho; }
//! \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
double rho; //!< Mass density
friend class PlasticMaterial;
};
/*!
\brief Class representing a plasticity material model.
*/
class PlasticMaterial : public Material
{
public:
//! \brief Constructor initializing the material parameters.
PlasticMaterial(const PlasticPrm* prm, unsigned short int n = 0);
//! \brief Empty destructor.
virtual ~PlasticMaterial() {}
//! \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[in] F Deformation gradient at current point
//! \param[in] prm Nonlinear solution algorithm parameters
virtual bool evaluate(Matrix& C, SymmTensor& sigma, double&,
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
double HVc[10]; //!< History variables, current configuration
double HVp[10]; //!< History variables, previous configuration
Tensor Fp; //!< Deformation gradient, previous configuration
};
#endif

143
Apps/FiniteDefElasticity/PlasticityUL.C Normal file
View File

@ -0,0 +1,143 @@
// $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
{
material->print(os);
std::cout <<"PlasticityUL: Updated Lagrangian formulation"<< std::endl;
}
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 Normal file
View File

@ -0,0 +1,72 @@
// $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

78
Apps/FiniteDefElasticity/SIMFiniteDefEl.C
View File

@ -15,6 +15,8 @@
#include "LinIsotropic.h"
#include "LinearMaterial.h"
#include "NeoHookeMaterial.h"
#include "PlasticMaterial.h"
#include "PlasticityUL.h"
#include "NonlinearElasticityULMixed.h"
#include "NonlinearElasticityULMX.h"
#include "NeoHookeElasticity.h"
@ -30,6 +32,9 @@ 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);
@ -58,7 +63,7 @@ bool SIMFiniteDefEl2D::parse (char* keyWord, std::istream& is)
{
if (!strncasecmp(keyWord,"ISOTROPIC",9))
{
int nmat = atoi(keyWord+10);
int nmat = atoi(keyWord+9);
if (myPid == 0)
std::cout <<"\nNumber of isotropic materials: "<< nmat << std::endl;
@ -68,6 +73,7 @@ bool SIMFiniteDefEl2D::parse (char* keyWord, std::istream& is)
int code = atoi(strtok(cline," "));
if (code > 0)
this->setPropertyType(code,Property::MATERIAL,mVec.size());
double E = atof(strtok(NULL," "));
double nu = atof(strtok(NULL," "));
double rho = atof(strtok(NULL," "));
@ -81,6 +87,38 @@ bool SIMFiniteDefEl2D::parse (char* keyWord, std::istream& is)
<< E <<" "<< nu <<" "<< rho <<" ("<< matVer <<")"<< std::endl;
}
}
else if (!strncasecmp(keyWord,"PLASTIC",7))
{
int nmat = atoi(keyWord+7);
if (myPid == 0)
std::cout <<"\nNumber of plastic materials: "<< nmat << std::endl;
char* cline = 0;
for (int i = 0; i < nmat && (cline = utl::readLine(is)); i++)
{
int code = atoi(strtok(cline," "));
if (code > 0)
this->setPropertyType(code,Property::MATERIAL,mVec.size());
RealArray pMAT;
double rho = atof(strtok(NULL," "));
while ((cline = strtok(NULL," ")))
{
pMAT.push_back(rho);
rho = atof(cline);
}
mVec.push_back(new PlasticPrm(pMAT,rho));
if (myPid == 0)
{
std::cout <<"\tMaterial code "<< code <<":";
for (size_t i = 0; i < pMAT.size(); i++)
std::cout <<" "<< pMAT[i];
std::cout <<" rho = "<< rho << std::endl;
}
}
}
else
return this->SIMLinEl2D::parse(keyWord,is);
@ -103,6 +141,9 @@ SIMFiniteDefEl3D::SIMFiniteDefEl3D (bool checkRHS,
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();
@ -137,7 +178,7 @@ bool SIMFiniteDefEl3D::parse (char* keyWord, std::istream& is)
{
if (!strncasecmp(keyWord,"ISOTROPIC",9))
{
int nmat = atoi(keyWord+10);
int nmat = atoi(keyWord+9);
if (myPid == 0)
std::cout <<"\nNumber of isotropic materials: "<< nmat << std::endl;
@ -147,6 +188,7 @@ bool SIMFiniteDefEl3D::parse (char* keyWord, std::istream& is)
int code = atoi(strtok(cline," "));
if (code > 0)
this->setPropertyType(code,Property::MATERIAL,mVec.size());
double E = atof(strtok(NULL," "));
double nu = atof(strtok(NULL," "));
double rho = atof(strtok(NULL," "));
@ -160,6 +202,38 @@ bool SIMFiniteDefEl3D::parse (char* keyWord, std::istream& is)
<< E <<" "<< nu <<" "<< rho <<" ("<< matVer <<")"<< std::endl;
}
}
else if (!strncasecmp(keyWord,"PLASTIC",7))
{
int nmat = atoi(keyWord+7);
if (myPid == 0)
std::cout <<"\nNumber of plastic materials: "<< nmat << std::endl;
char* cline = 0;
for (int i = 0; i < nmat && (cline = utl::readLine(is)); i++)
{
int code = atoi(strtok(cline," "));
if (code > 0)
this->setPropertyType(code,Property::MATERIAL,mVec.size());
RealArray pMAT;
double rho = atof(strtok(NULL," "));
while ((cline = strtok(NULL," ")))
{
pMAT.push_back(rho);
rho = atof(cline);
}
mVec.push_back(new PlasticPrm(pMAT,rho));
if (myPid == 0)
{
std::cout <<"\tMaterial code "<< code <<":";
for (size_t i = 0; i < pMAT.size(); i++)
std::cout <<" "<< pMAT[i];
std::cout <<" rho = "<< rho << std::endl;
}
}
}
else
return this->SIMLinEl3D::parse(keyWord,is);

3
Apps/FiniteDefElasticity/SIMFiniteDefEl.h
View File

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

5
src/ASM/IntegrandBase.h
View File

@ -55,6 +55,11 @@ public:
//! \brief Defines the solution mode before the element assembly is started.
virtual void setMode(SIM::SolutionMode) {}
//! \brief Initializes the integrand for a new integration loop.
//! \details This method is invoked once before starting the numerical
//! integration over the entire spatial domain.
virtual void initIntegration(const TimeDomain&) {}
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
//! \param[in] X0 Cartesian coordinates of the element center

11
src/ASM/TimeDomain.h
View File

@ -1,4 +1,4 @@
// $Id: TimeDomain.h,v 1.1 2010-06-01 09:40:25 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file TimeDomain.h
@ -21,12 +21,13 @@
struct TimeDomain
{
double t; //!< Current time (or pseudo time)
double dt; //!< Current time increment (or load parameter increment)
int it; //!< Current iteration within current time/load step
double t; //!< Current time (or pseudo time, load parameter)
double dt; //!< Current time (or load parameter) increment
int it; //!< Current iteration within current time/load step
bool first; //!< If \e true, this is the first load/time step
//! \brief Default constructor.
TimeDomain() { t = dt = 0.0; it = 0; }
TimeDomain() { t = dt = 0.0; it = 0; first = true; }
};
#endif

15
src/LinAlg/SystemMatrix.h
View File

@ -76,6 +76,11 @@ public:
//! \brief Initializes the vector assuming it is properly dimensioned.
virtual void init(real value = real(0)) = 0;
//! \brief Begins communication step needed in parallel vector assembly.
virtual bool beginAssembly() { return true; }
//! \brief Ends communication step needed in parallel vector assembly.
virtual bool endAssembly() { return true; }
//! \brief Multiplication with a scalar.
virtual void mult(real alpha) = 0;
@ -146,6 +151,11 @@ public:
//! \brief Initializes the vector to a given scalar value.
virtual void init(real value = real(0)) { this->fill(value); }
//! \brief Begins communication step needed in parallel vector assembly.
virtual bool beginAssembly() { return true; }
//! \brief Ends communication step needed in parallel vector assembly.
virtual bool endAssembly() { return true; }
//! \brief Multiplication with a scalar.
virtual void mult(real alpha) { this->operator*=(alpha); }
@ -212,6 +222,11 @@ public:
//! \brief Initializes the matrix to zero assuming it is properly dimensioned.
virtual void init() = 0;
//! \brief Begins communication step needed in parallel matrix assembly.
virtual bool beginAssembly() { return true; }
//! \brief Ends communication step needed in parallel matrix assembly.
virtual bool endAssembly() { return true; }
//! \brief Adds an element stiffness matrix into the system stiffness matrix.
//! \param[in] eM The element stiffness matrix
//! \param[in] sam Auxilliary data describing the FE model topology,

17
src/SIM/NonLinSIM.C
View File

@ -1,4 +1,4 @@
// $Id: NonLinSIM.C,v 1.28 2011-02-05 18:27:02 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file NonLinSIM.C
@ -18,11 +18,8 @@
#include <sstream>
NonLinSIM::NonLinSIM (SIMbase* sim)
NonLinSIM::NonLinSIM (SIMbase* sim) : model(sim), nBlock(0)
{
model = sim;
nBlock = 0;
#if SP_DEBUG > 2
msgLevel = 100; // prints the linear solution vector if size < 100
#else
@ -37,7 +34,7 @@ NonLinSIM::NonLinSIM (SIMbase* sim)
convTol = 1.0e-6;
divgLim = 1.0;
eta = 0.0;
iteNorm = ENERGY;
iteNorm = ENERGY;
}
@ -159,6 +156,7 @@ bool NonLinSIM::solveStep (SolvePrm& param, SIM::SolutionMode mode)
if (!this->updateConfiguration(param))
return false;
param.time.first = false;
model->setMode(SIM::RECOVERY);
return this->solutionNorms(param.time);
@ -264,12 +262,9 @@ NonLinSIM::ConvStatus NonLinSIM::checkConvergence (SolvePrm& param)
static double prevNorm = 0.0;
static int nIncrease = 0;
double norm=1.f, enorm, resNorm, linsolNorm;
double enorm, resNorm, linsolNorm;
model->iterationNorms(linsol,residual,enorm,resNorm,linsolNorm);
if (iteNorm == ENERGY)
norm = enorm;
else if (iteNorm == L2)
norm = resNorm;
double norm = iteNorm == ENERGY ? enorm : resNorm;
if (param.iter == 0)
{

59
src/SIM/NonLinSIM.h
View File

@ -1,4 +1,4 @@
// $Id: NonLinSIM.h,v 1.15 2011-01-28 15:28:54 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file NonLinSIM.h
@ -16,9 +16,8 @@
#include "SIMinput.h"
#include "SIMenums.h"
#include "TimeDomain.h"
#include "MatVec.h"
#include "SIMparameters.h"
#include "MatVec.h"
class SIMbase;
@ -39,27 +38,23 @@ public:
//! \brief The destructor frees the dynamically allocated FE model object.
virtual ~NonLinSIM();
/*!
\brief A class for nonlinear solution parameters.
*/
//! \brief A class for nonlinear solution parameters.
class SolvePrm : public SIMparameters
{
public:
int maxit; //!< Maximum number of iterations
int nupdat; //!< Number of iterations with updated tangent
double convTol; //!< Relative convergence tolerance
double divgLim; //!< Relative divergence limit
double refNorm; //!< Reference energy norm used in convergence checks
double alpha; //!< Iteration acceleration parameter (line search)
double eta; //!< Line searce tolerance
SolvePrm()
{
refNorm = alpha = 1;
}
public:
//! \brief Default constructor.
SolvePrm() : refNorm(1.0), alpha(1.0) {}
int maxit; //!< Maximum number of iterations
int nupdat; //!< Number of iterations with updated tangent
double convTol; //!< Relative convergence tolerance
double divgLim; //!< Relative divergence limit
double refNorm; //!< Reference energy norm used in convergence checks
double alpha; //!< Iteration acceleration parameter (line search)
double eta; //!< Line search tolerance
};
//! \brief Initializes solution parameters object with values read from file.
//! \brief Initializes the solution parameters with values read from file.
virtual void init(SolvePrm& param);
//! \brief Advances the time/load step one step forward.
virtual bool advanceStep(SolvePrm& param);
@ -75,19 +70,15 @@ public:
//! \param[in] mode Solution mode to use for this step
bool solveStep(SolvePrm& param, SIM::SolutionMode mode = SIM::STATIC);
//! \brief Enum describing the norm used for convergence tests in the
//! Newton-Raphson iterations
enum ERRORNORM {
L2,
ENERGY
};
//! \brief Enum describing the norm used for convergence checks.
enum CNORM { L2, ENERGY };
//! \brief Set the current norm in used in iteration error estimates
//! \brief Sets the current norm in used in iteration error estimates.
//! param[in] norm The wanted error norm
void setErrorNorm(ERRORNORM norm) { iteNorm = norm; }
void setErrorNorm(CNORM norm) { iteNorm = norm; }
//! \brief Computes and prints some solution norm quantities.
//! \param[in] time Parameters for nonlinear/time-dependent simulations.
//! \param[in] time Parameters for nonlinear/time-dependent simulations
//! \param[in] compName Solution name to be used in the norm output
bool solutionNorms(const TimeDomain& time,
const char* compName = "displacement");
@ -102,7 +93,7 @@ public:
//! \brief Dumps the primary solution to given stream for inspection.
//! \param[in] iStep Load/time step identifier
//! \param[in] time Current time/load parameter
//! \param[in] os The output streeam to write to.
//! \param[in] os The output streeam to write to
//! \param[in] withID If \e true, write node ID and coordinates too
void dumpStep(int iStep, double time, std::ostream& os,
bool withID = true) const;
@ -127,9 +118,9 @@ private:
SIMbase* model; //!< The isogeometric FE model
protected:
Vector linsol; //!< Linear solution vector
Vector residual; //!< Residual force vector
std::vector<Vector> solution; //!< Total solution vectors
Vector linsol; //!< Linear solution vector
Vector residual; //!< Residual force vector
Vectors solution; //!< Total solution vectors
private:
// Nonlinear solution algorithm parameters
@ -141,7 +132,7 @@ private:
double eta; //!< Line search tolerance
int maxit; //!< Maximum number of iterations in a time/load step
int nupdat; //!< Number of iterations with updated tangent
ERRORNORM iteNorm;//!< The norm used to measure the residual
CNORM iteNorm; //!< The norm type used to measure the residual
// Post-processing attributes
int nBlock; //!< Running VTF result block counter

37
src/SIM/SIMbase.C
View File

@ -13,7 +13,6 @@
#include "SIMbase.h"
#include "ASMbase.h"
#include "PETScMatrix.h"
#ifdef PARALLEL_PETSC
#include "SAMpatchPara.h"
#include "petscksp.h"
@ -94,10 +93,8 @@ SIMbase::~SIMbase ()
bool SIMbase::parse (char* keyWord, std::istream& is)
{
if (!strncasecmp(keyWord,"LINEARSOLVER",12))
{
if (!mySolParams) mySolParams = new LinSolParams();
mySolParams->read(is,atoi(keyWord+12));
}
this->readLinSolParams(is,atoi(keyWord+12));
else if (!strncasecmp(keyWord,"PARTITIONING",12))
{
int nproc = atoi(keyWord+12);
@ -129,6 +126,7 @@ bool SIMbase::parse (char* keyWord, std::istream& is)
return false;
}
}
else if (isalpha(keyWord[0]))
std::cerr <<" *** SIMbase::parse: Unknown keyword: "<< keyWord << std::endl;
@ -410,6 +408,9 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol,
{
PROFILE1("Element assembly");
if (!myProblem) return false;
myProblem->initIntegration(time);
myEqSys->init(newLHSmatrix);
bool ok = true;
@ -1257,26 +1258,26 @@ void SIMbase::readLinSolParams (std::istream& is, int npar)
bool SIMbase::finalizeAssembly (bool newLHSmatrix)
{
// Communication of matrix and vector assembly (for PETSC only)
SystemMatrix* A = myEqSys->getMatrix();
SystemVector* b = myEqSys->getVector();
if (A && b && A->getType() == SystemMatrix::PETSC)
if (A && newLHSmatrix)
{
// Communication of matrix and vector assembly (for PETSC only)
#ifdef HAS_PETSC
if (newLHSmatrix)
{
if (!static_cast<PETScMatrix*>(A)->beginAssembly()) return false;
if (!static_cast<PETScMatrix*>(A)->endAssembly()) return false;
}
if (!static_cast<PETScVector*>(b)->beginAssembly()) return false;
if (!static_cast<PETScVector*>(b)->endAssembly()) return false;
if (!A->beginAssembly()) return false;
if (!A->endAssembly()) return false;
#if SP_DEBUG > 3
std::cout <<"\nSystem coefficient matrix:"<< *A;
#endif
}
SystemVector* b = myEqSys->getVector();
if (b)
{
if (!b->beginAssembly()) return false;
if (!b->endAssembly()) return false;
#if SP_DEBUG > 3
if (A && newLHSmatrix) std::cout <<"\nSystem coefficient matrix:"<< *A;
if (b) std::cout <<"\nSystem right-hand-side vector:"<< *b;
std::cout <<"\nSystem right-hand-side vector:"<< *b;
#endif
}
return true;
}