From fae5ff6b96a3b26c736852eaef7fc9698639af9d Mon Sep 17 00:00:00 2001 From: kmo Date: Thu, 20 Oct 2011 16:56:07 +0000 Subject: [PATCH] Extensions for axisymmetric formulations git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1270 e10b68d5-8a6e-419e-a041-bce267b0401d --- Apps/FiniteDefElasticity/Material/cons2d.f | 26 +++-- Apps/FiniteDefElasticity/Material/plas2d.f | 26 +++-- Apps/FiniteDefElasticity/NeoHookeElasticity.C | 2 +- Apps/FiniteDefElasticity/NeoHookeMaterial.C | 14 +-- .../FiniteDefElasticity/NonlinearElasticity.C | 6 +- .../NonlinearElasticityFbar.C | 60 ++++++---- .../NonlinearElasticityTL.C | 43 +++++--- .../NonlinearElasticityTL.h | 5 +- .../NonlinearElasticityUL.C | 104 ++++++++++++------ .../NonlinearElasticityUL.h | 14 ++- .../NonlinearElasticityULMX.C | 24 ++-- .../NonlinearElasticityULMixed.C | 60 ++++++---- Apps/FiniteDefElasticity/PlasticMaterial.C | 15 ++- Apps/FiniteDefElasticity/SIMFiniteDefEl.C | 2 +- Apps/FiniteDefElasticity/main_NonLinEl.C | 7 +- 15 files changed, 261 insertions(+), 147 deletions(-) diff --git a/Apps/FiniteDefElasticity/Material/cons2d.f b/Apps/FiniteDefElasticity/Material/cons2d.f index 1dbfb468..4528c239 100644 --- a/Apps/FiniteDefElasticity/Material/cons2d.f +++ b/Apps/FiniteDefElasticity/Material/cons2d.f @@ -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,15 +27,19 @@ C return end if 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) - if (nsig .eq. 4) then - Sig(1:4) = S3D(1:4) + if (n .eq. 4) then + Cst = C3D(1:4,1:4) else - Sig(1:2) = S3D(1:2) - Sig(3) = S3D(4) + 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 + Sig(1:2) = S3D(1:2) + Sig(3) = S3D(4) end if C return diff --git a/Apps/FiniteDefElasticity/Material/plas2d.f b/Apps/FiniteDefElasticity/Material/plas2d.f index b2c28169..ce92f83a 100644 --- a/Apps/FiniteDefElasticity/Material/plas2d.f +++ b/Apps/FiniteDefElasticity/Material/plas2d.f @@ -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,15 +34,19 @@ C return end if 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) - if (nsig .eq. 4) then - Sig(1:4) = S3D + if (n .eq. 4) then + Cst = C3D(1:4,1:4) else - Sig(1:2) = S3D(1:2) - Sig(3) = S3D(4) + 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 + Sig(1:2) = S3D(1:2) + Sig(3) = S3D(4) end if C return diff --git a/Apps/FiniteDefElasticity/NeoHookeElasticity.C b/Apps/FiniteDefElasticity/NeoHookeElasticity.C index 9c0d6475..07fffeaf 100644 --- a/Apps/FiniteDefElasticity/NeoHookeElasticity.C +++ b/Apps/FiniteDefElasticity/NeoHookeElasticity.C @@ -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(); diff --git a/Apps/FiniteDefElasticity/NeoHookeMaterial.C b/Apps/FiniteDefElasticity/NeoHookeMaterial.C index 7fdd8dc3..bc29953f 100644 --- a/Apps/FiniteDefElasticity/NeoHookeMaterial.C +++ b/Apps/FiniteDefElasticity/NeoHookeMaterial.C @@ -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 diff --git a/Apps/FiniteDefElasticity/NonlinearElasticity.C b/Apps/FiniteDefElasticity/NonlinearElasticity.C index e612f3df..b136b8e4 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticity.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticity.C @@ -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 diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityFbar.C b/Apps/FiniteDefElasticity/NonlinearElasticityFbar.C index 8fce0c04..a2e9ef39 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityFbar.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityFbar.C @@ -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); - Fi.fill(F.ptr()); + 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); - Fi.fill(F.ptr()); + 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(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); } diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityTL.C b/Apps/FiniteDefElasticity/NonlinearElasticityTL.C index 42150a90..82999b31 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityTL.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityTL.C @@ -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 - eKm->multiply(Bmat,CB,true,false,true); // EK += B^T * CB + 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 diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityTL.h b/Apps/FiniteDefElasticity/NonlinearElasticityTL.h index 070b26a3..da8dd782 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityTL.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticityTL.h @@ -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 diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityUL.C b/Apps/FiniteDefElasticity/NonlinearElasticityUL.C index da2b4e3a..971cba40 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityUL.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityUL.C @@ -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); - Fi.fill(F.ptr()); + 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(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 - eKm->multiply(Bmat,CB,true,false,true); // EK += B^T * CB + 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(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(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; } diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityUL.h b/Apps/FiniteDefElasticity/NonlinearElasticityUL.h index 0c2c15af..743bd705 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityUL.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticityUL.h @@ -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 diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C b/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C index 42e852c2..80db6829 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C @@ -155,22 +155,26 @@ bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt, // Evaluate the deformation gradient, Fp, at previous configuration const_cast(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(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); - 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); + 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) diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C index b04a7c29..1044474f 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C @@ -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); - Fi.fill(F.ptr()); + 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 @@ -448,14 +457,16 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt, myMats->A[Kup](nsd*(a-1)+i,b) += dNdx(a,i)*fe.N2(b) * dVup; } - 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[Ktt].outer_product(fe.N2,fe.N2*Dmat(7,7),true); // += N2*N2^T*D77 + 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(&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); } diff --git a/Apps/FiniteDefElasticity/PlasticMaterial.C b/Apps/FiniteDefElasticity/PlasticMaterial.C index ced2e15c..84356c7f 100644 --- a/Apps/FiniteDefElasticity/PlasticMaterial.C +++ b/Apps/FiniteDefElasticity/PlasticMaterial.C @@ -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(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, diff --git a/Apps/FiniteDefElasticity/SIMFiniteDefEl.C b/Apps/FiniteDefElasticity/SIMFiniteDefEl.C index f3519a4a..e0dccd51 100644 --- a/Apps/FiniteDefElasticity/SIMFiniteDefEl.C +++ b/Apps/FiniteDefElasticity/SIMFiniteDefEl.C @@ -45,7 +45,7 @@ SIMFiniteDefEl2D::SIMFiniteDefEl2D (const std::vector& 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); diff --git a/Apps/FiniteDefElasticity/main_NonLinEl.C b/Apps/FiniteDefElasticity/main_NonLinEl.C index 2517ec58..0a4d91e5 100644 --- a/Apps/FiniteDefElasticity/main_NonLinEl.C +++ b/Apps/FiniteDefElasticity/main_NonLinEl.C @@ -61,6 +61,7 @@ extern std::vector 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] <<" [-dense|-spr|-superlu[]|-samg|-petsc]\n " - <<" [-lag] [-spec] [-nGauss ] [-2D[pstrain]] [-UL|-MX [

]" + <<" [-lag|-spec] [-nGauss ] [-2D[pstrain|axis]] [-UL|-MX [

]" <<"|-[M|m]ixed|-Fbar ]\n [-vtf [-nviz ]" <<" [-nu ] [-nv ] [-nw ]] [-hdf5]\n " <<" [-saveInc ] [-skip2nd] [-dumpInc [raw]]"