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

Various finite deformation fixes/adjustments. Mainly related to Fbar method

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1116 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo 2011-08-28 13:52:12 +00:00 committed by Knut Morten Okstad
parent 1a093c5b25
commit d798dfcd6c
7 changed files with 74 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
Apps/FiniteDefElasticity/NeoHookeMaterial.C
View File

@ -136,11 +136,22 @@ bool NeoHookeMaterial::evaluate (Matrix& C, SymmTensor& sigma, double& U,
sigma.transform(Fi); // sigma = F^-1 * sigma * F^-t
sigma *= J;
//TODO: Also pull-back the C-matrix (Total Lagrange formulation)
std::cerr <<" *** NeoHookeMaterial::evaluate: Not available for"
<<" Total Lagrangian formulation, sorry."<< std::endl;
return false;
}
else
{
SymmTensor S(sigma); // Make a copy of sigma since it should be Cauchy stress when iop=3
U = S.transform(Fi).innerProd(eps)*J; //TODO: Replace this by proper path integral
SymmTensor S(sigma); // sigma should be Cauchy stress when iop=3
U = S.transform(Fi).innerProd(eps)*J;
//TODO: Replace the above by proper path integral
static bool first = true;
if (first)
std::cerr <<" ** NeoHookeMaterial::evaluate: Path-integration of"
<<" strain energy density is not implemented.\n"
<<" Warning: The strain energy will be incorrect."
<< std::endl;
first = false;
}
}

31
Apps/FiniteDefElasticity/NonlinearElasticityFbar.C
View File

@ -316,15 +316,12 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
{
// Evaluate the Lagrange polynomials extrapolating the volumetric
// sampling points (assuming a regular distribution over the element)
RealArray Nbar;
Vector Nbar;
int pbar3 = nsd > 2 ? pbar : 0;
if (Lagrange::computeBasis(Nbar,pbar,fe.xi*scale,pbar,fe.eta*scale,
pbar3,fe.zeta*scale)) return false;
if (!Lagrange::computeBasis(Nbar,pbar,fe.xi*scale,pbar,fe.eta*scale,
pbar3,fe.zeta*scale)) return false;
#ifdef INT_DEBUG
std::cout <<"NonlinearElasticityFbar::Nbar =";
for (size_t i = 0; i < Nbar.size(); i++) std::cout <<" "<< Nbar[i];
std::cout << std::endl;
std::cout <<"NonlinearElasticityFbar::Nbar ="<< Nbar;
#endif
// Compute modified deformation gradient determinant and basis function
@ -380,8 +377,7 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
if (eKm)
{
// Compute standard and mixed discrete spatial gradient operators
Matrix G, Gbar;
// Compute standard and modified discrete spatial gradient operators
getGradOperator(dNdx,G);
getGradOperator(dMdx,Gbar);
@ -412,11 +408,12 @@ bool NonlinearElasticityFbar::evalInt (LocalIntegral*& elmInt,
<<"NonlinearElasticityFbar::Q ="<< Q;
#endif
// Compute and accumulate contribution to tangent stiffness matrix
// Compute and accumulate contribution to tangent stiffness matrix.
// First, standard (material and geometric) tangent stiffness terms
eKm->multiply(G,CB.multiply(A,G),true,false,true); // EK += G^T * A*G
// Modify the tangent stiffness for the F-bar terms
// Then, modify the tangent stiffness for the F-bar terms
Gbar -= G;
eKm->multiply(G,CB.multiply(Q,Gbar),true,false,true); // EK += G^T * Q*Gbar
}
@ -451,10 +448,10 @@ bool ElasticityNormFbar::evalInt (LocalIntegral*& elmInt,
const TimeDomain& prm,
const Vec3& X) const
{
size_t nsd = fe.dNdX.cols();
NonlinearElasticityFbar& p = static_cast<NonlinearElasticityFbar&>(myProblem);
// Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
size_t nsd = fe.dNdX.cols();
Tensor F(nsd);
SymmTensor E(nsd);
if (!p.kinematics(fe.dNdX,F,E))
@ -484,14 +481,14 @@ bool ElasticityNormFbar::evalInt (LocalIntegral*& elmInt,
Tensor Fbar(F);
Fbar *= pow(fabs(Jbar/F.det()),1.0/(double)nsd);
// Compute the strain energy density, U(E) = Int_E (Sig:Eps) dEps
// and the Cauchy stress tensor, S
// Compute the strain energy density, U(E) = Int_E (S:Eps) dEps
// and the Cauchy stress tensor, sigma
double U = 0.0;
SymmTensor S(nsd,p.material->isPlaneStrain());
if (!p.material->evaluate(p.Cmat,S,U,X,Fbar,E,3,&prm,&F))
SymmTensor sigma(nsd,p.material->isPlaneStrain());
if (!p.material->evaluate(p.Cmat,sigma,U,X,Fbar,E,3,&prm,&F))
return false;
// Integrate the norms
return ElasticityNormUL::evalInt(getElmNormBuffer(elmInt,6),
S,U,F.det(),fe.detJxW);
sigma,U,F.det(),fe.detJxW);
}

13
Apps/FiniteDefElasticity/NonlinearElasticityFbar.h
View File

@ -36,6 +36,7 @@ class NonlinearElasticityFbar : public NonlinearElasticityUL
public:
//! \brief The default constructor invokes the parent class constructor.
//! \param[in] n Number of spatial dimensions
//! \param[in] nvp Number of volumetric sampling points in each direction
NonlinearElasticityFbar(unsigned short int n = 3, int nvp = 1);
//! \brief Empty destructor.
virtual ~NonlinearElasticityFbar() {}
@ -58,11 +59,6 @@ public:
//! \param[in] fe Finite element data of current integration point
//! \param[in] prm Nonlinear solution algorithm parameters
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This method mainly stores the integration point quantities
//! depending on the basis function values in the internal member \a myData.
//! The actual numerical integration of the tangent stiffness and internal
//! forces is then performed by the \a finalizeElement method.
virtual bool evalInt(LocalIntegral*& elmInt, const FiniteElement& fe,
const TimeDomain& prm, const Vec3& X) const;
@ -80,10 +76,13 @@ private:
int pbar; //!< Polynomial order of the internal volumetric data field
double scale; //!< Scaling factor for extrapolation from sampling points
mutable size_t iP; //!< Local sampling point counter
mutable Matrix dMdx; //!< Modified basis function gradients
mutable std::vector<VolPtData> myVolData; //!< Volumetric sampling point data
mutable size_t iP; //!< Volumetric sampling point counter
mutable Matrix dMdx; //!< Modified basis function gradients
mutable Matrix G; //!< Discrete gradient operator
mutable Matrix Gbar; //!< Modified discrete gradient operator
friend class ElasticityNormFbar;
};

9
Apps/FiniteDefElasticity/NonlinearElasticityUL.h
View File

@ -142,6 +142,15 @@ public:
virtual size_t getNoFields() const { return 6; }
protected:
//! \brief Evaluates and accumulates the point-wise norm quantities.
//! \param pnorm The element norms
//! \param[in] S Cauchy stress tensor
//! \param[in] U Strain energy density
//! \param[in] detF Determinant of deformation gradient
//! \param[in] detJxW Jacobian determinant times integration point weight
//!
//! \details This method is used by the virtual \a evalInt method and is
//! separated out such that it also can be reused by sub-classes.
static bool evalInt(ElmNorm& pnorm, const SymmTensor& S,
double U, double detF, double detJxW);

7
Apps/FiniteDefElasticity/PlasticMaterial.C
View File

@ -183,11 +183,14 @@ bool PlasticMaterial::evaluate (Matrix& C, SymmTensor& sigma, double& U,
sigma.transform(Fi); // sigma = F^-1 * sigma * F^-t
sigma *= J;
//TODO: Also pull-back the C-matrix (Total Lagrange formulation)
std::cerr <<" *** PlasticMaterial::evaluate: Not available for"
<<" Total Lagrangian formulation, sorry."<< std::endl;
return false;
}
else
{
SymmTensor S(sigma); // Make a copy of sigma since it should be Cauchy stress when iop=3
S.transform(Fi); // S = F^-1 * S * F^-t
SymmTensor S(sigma); // sigma should be Cauchy stress when iop=3
S.transform(Fi); // S = F^-1 * sigma * F^-t
if (iAmIntegrating)
U = itgPoints[iP1-1]->energyIntegral(S,eps)*J;
else

22
src/LinAlg/SAM.C
View File

@ -286,10 +286,6 @@ bool SAM::getNoDofCouplings (IntVec& nnz) const
{
nnz.clear();
// RUNAR
// Not implemented for constrained equations yet
//if (nceq > 0) return false;
// Find the set of DOF couplings for each free DOF
std::vector<IntSet> dofc;
if (!this->getDofCouplings(dofc))
@ -531,7 +527,12 @@ void SAM::assembleReactions (Vector& reac, const RealArray& eS, int iel) const
bool SAM::getElmEqns (IntVec& meen, int iel, int nedof) const
{
if (iel < 1 || iel > nel) return false;
if (iel < 1 || iel > nel)
{
std::cerr <<"SAM::getElmEqns: Element "<< iel <<" is out of range [1,"
<< nel <<"]"<< std::endl;
return false;
}
int ip = mpmnpc[iel-1];
int nenod = mpmnpc[iel] - ip;
@ -564,7 +565,12 @@ bool SAM::getElmEqns (IntVec& meen, int iel, int nedof) const
bool SAM::getNodeEqns (IntVec& mnen, int inod) const
{
if (inod < 1 || inod > nnod) return false;
if (inod < 1 || inod > nnod)
{
std::cerr <<"SAM::getNodeEqns: Node "<< inod <<" is out of range [1,"
<< nnod <<"]"<< std::endl;
return false;
}
mnen.clear();
mnen.reserve(madof[inod]-madof[inod-1]);
@ -749,7 +755,11 @@ real SAM::getReaction (int dir, const Vector& rf) const
bool SAM::getNodalReactions (int inod, const Vector& rf, Vector& nrf) const
{
if (inod < 1 || inod > nnod)
{
std::cerr <<"SAM::getNodalReactions: Node "<< inod <<" is out of range [1,"
<< nnod <<"]"<< std::endl;
return false;
}
bool haveRF = false;
int ip = madof[inod-1]-1;

16
src/Utility/Tensor.C
View File

@ -74,7 +74,7 @@ Tensor& Tensor::operator= (const Tensor& T)
for (t_ind j = (this->symmetric() ? i : 1); j <= ndim; j++)
v[this->index(i,j)] = T(i,j);
if (v.size() == 4 && T.v.size() >= 4)
if (this->symmetric() && v.size() == 4 && T.v.size() >= 4)
v[2] = T(3,3);
}
@ -106,6 +106,9 @@ Tensor& Tensor::operator= (real val)
for (i = j = 0; i < n; i++, j += inc)
v[j] = val;
if (inc == 1 && v.size() == 4)
v[2] = val;
return *this;
}
@ -130,6 +133,9 @@ Tensor& Tensor::operator+= (real val)
for (i = j = 0; i < n; i++, j += inc)
v[j] += val;
if (inc == 1 && v.size() == 4)
v[2] += val;
return *this;
}
@ -154,6 +160,9 @@ Tensor& Tensor::operator-= (real val)
for (i = j = 0; i < n; i++, j += inc)
v[j] -= val;
if (inc == 1 && v.size() == 4)
v[2] -= val;
return *this;
}
@ -455,16 +464,13 @@ real SymmTensor::trace () const
{
real t = real(0);
if (n == 3)
if (n == 3 || v.size() == 4)
t = v[0] + v[1] + v[2];
else if (n == 2)
t = v[0] + v[1];
else if (n == 1)
t = v[0];
if (v.size() == 4)
t += v[2];
return t;
}