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

Corrected the von Mises stress calculation in 2D and plane strain in particular; must include the sigma_zz component also

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@917 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo
2011-04-26 12:30:53 +00:00
committed by Knut Morten Okstad
parent 67a0212aa4
commit 74c7d26258
13 changed files with 190 additions and 101 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
Apps/FiniteDefElasticity/CMakeLists.txt
View File

@@ -13,6 +13,11 @@ IF(NOT CMAKE_BUILD_TYPE)
SET(CMAKE_BUILD_TYPE Release)
ENDIF(NOT CMAKE_BUILD_TYPE)
SET(IFEM_BUILD_TYPE ${CMAKE_BUILD_TYPE})
IF(${CMAKE_BUILD_TYPE} MATCHES "Nopt")
SET(CMAKE_BUILD_TYPE Debug)
ENDIF(${CMAKE_BUILD_TYPE} MATCHES "Nopt")
# Add local modules
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH}
${PROJECT_SOURCE_DIR}/../../cmake/Modules

3
Apps/FiniteDefElasticity/LinearMaterial.h
View File

@@ -37,6 +37,9 @@ public:
//! \brief The destructor deletes the material properties object.
virtual ~LinearMaterial() { delete material; }
//! \brief Returns \e false if plane stress in 2D.
virtual bool isPlaneStrain() const { return material->isPlaneStrain(); }
//! \brief Prints out material parameters to the given output stream.
virtual void print(std::ostream& os) const { material->print(os); }

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

@@ -1,12 +1,12 @@
subroutine cons2d (ipsw,iter,iwr,lfirst,mTYP,mVER,
& ndf,detF,F,pMAT,U,Sig,Cst,ierr)
& 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, ndf, ierr
real*8 detF, F(*), pMAT(*), U, Sig(3), Cst(3,3)
integer ipsw, iter, iwr, lfirst, mTYP, mVER, nsig, ndf, ierr
real*8 detF, F(*), pMAT(*), U, Sig(4), Cst(3,3)
real*8 F3D(3,3), S3D(6), C3D(6,6)
C
if (ndf .eq. 2) then
@@ -31,8 +31,12 @@ C
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)
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
end

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

@@ -1,13 +1,13 @@
subroutine plas2d (ipsw,iwr,iter,lfirst,pmat,ndf,detF,
subroutine plas2d (ipsw,iwr,iter,lfirst,pmat,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, ndf, ierr
integer ipsw, iwr, iter, lfirst, nsig, ndf, ierr
real*8 pmat(*), detF, Fn1(*), Fn(*)
real*8 be(*), Epp, Epl(*), Sig(3), Cst(3,3)
real*8 be(*), Epp, Epl(*), Sig(4), Cst(3,3)
real*8 F3D(3,3), F3P(3,3), S3D(4), C3D(6,6)
C
if (ndf .eq. 2) then
@@ -38,8 +38,12 @@ C
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)
if (nsig .eq. 4) then
Sig(1:4) = S3D
else
Sig(1:2) = S3D(1:2)
Sig(3) = S3D(4)
end if
C
return
end

8
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& nDF, const double& detF, const double* F,
const double* pMAT,
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,
@@ -91,10 +91,10 @@ bool NeoHookeMaterial::evaluate (Matrix& C, SymmTensor& sigma, double& U,
#ifdef USE_FTNMAT
// Invoke the FORTRAN routine for Neo-Hookean hyperelastic material models.
if (ndim == 2)
cons2d_(INT_DEBUG,0,6,0,mTYP,mVER,F.dim(),J,F.ptr(),pmat,
cons2d_(INT_DEBUG,0,6,0,mTYP,mVER,sigma.size(),F.dim(),J,F.ptr(),pmat,
U,sigma.ptr(),C.ptr(),ierr);
else
cons3d_(INT_DEBUG,0,6,0,mTYP,mVER,0,6,0,J,F.ptr(),pmat,&U,
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;

4
Apps/FiniteDefElasticity/PlasticMaterial.C
View File

@@ -19,7 +19,7 @@ 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& nDF, const double& detF,
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);
@@ -251,7 +251,7 @@ bool PlasticMaterial::PlasticPoint::evaluate (Matrix& C,
std::cout << std::endl;
#endif
if (ndim == 2)
plas2d_(INT_DEBUG,6,prm.it,prm.first,&pMAT.front(),F.dim(),J,
plas2d_(INT_DEBUG,6,prm.it,prm.first,&pMAT.front(),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,

8
src/Integrands/AnalyticSolutions.C
View File

@@ -35,11 +35,11 @@ SymmTensor Hole::evaluate (const Vec3& X) const
double R2 = R <= a ? 1.0 : a*a/(R*R);
double R4 = R <= a ? 1.0 : R2*R2;
SymmTensor sigma(is3D ? 3 : 2);
SymmTensor sigma(is3D ? 3 : 2, true);
sigma(1,1) = F0 * (1.0 - R2*(1.5*C2 + C4) + 1.5*R4*C4);
sigma(2,2) = F0 * ( - R2*(0.5*C2 - C4) - 1.5*R4*C4);
sigma(1,2) = F0 * ( - R2*(0.5*S2 + S4) + 1.5*R4*S4);
if (is3D) sigma(3,3) = F0 * nu*(1.0 - 2.0*R2*C2);
sigma(3,3) = F0 * nu*(1.0 - 2.0*R2*C2);
return sigma;
}
@@ -98,12 +98,12 @@ SymmTensor Lshape::evaluate (const Vec3& X) const
double theta = atan2(y,x);
// Set up the stress tensor in local system
SymmTensor sigma(is3D ? 3 : 2);
SymmTensor sigma(is3D ? 3 : 2, true);
double c0 = F0*lambda*pow(r,lm1);
sigma(1,1) = c0*((2.0-q*lp1)*cos(lm1*theta) - lm1*cos(lm3*theta));
sigma(2,2) = c0*((2.0+q*lp1)*cos(lm1*theta) + lm1*cos(lm3*theta));
sigma(1,2) = c0*( q*lp1 *sin(lm1*theta) + lm1*sin(lm3*theta));
if (is3D) sigma(3,3) = nu * (sigma(1,1)+sigma(2,2));
sigma(3,3) = nu * (sigma(1,1)+sigma(2,2));
// Transform to global coordinates
return sigma.transform(T);

34
src/Integrands/Elasticity.C
View File

@@ -457,7 +457,7 @@ bool Elasticity::evalSol (Vector& s, const Vector&,
return false;
// Calculate the stress tensor through the constitutive relation
SymmTensor sigma(nsd); double U;
SymmTensor sigma(nsd,material->isPlaneStrain()); double U;
if (!material->evaluate(Cmat,sigma,U,X,dUdX,eps))
return false;
@@ -493,7 +493,7 @@ bool Elasticity::evalSol (Vector& s, const Matrix& dNdX, const Vec3& X) const
return false;
// Calculate the stress tensor through the constitutive relation
SymmTensor sigma(nsd); double U;
SymmTensor sigma(nsd,material->isPlaneStrain()); double U;
if (!material->evaluate(Cmat,sigma,U,X,dUdX,eps))
return false;
@@ -513,15 +513,16 @@ bool Elasticity::evalSol (Vector& s, const STensorFunc& asol,
size_t Elasticity::getNoFields () const
{
return nsd*(nsd+1)/2 + 1;
if (nsd == 2 && material->isPlaneStrain())
return 5;
else
return nsd*(nsd+1)/2 + 1;
}
const char* Elasticity::getFieldLabel (size_t i, const char* prefix) const
{
static const char* s1[1] = { "Axial stress" };
static const char* s2[3] = { "s_xx","s_yy","s_xy" };
static const char* s3[6] = { "s_xx","s_yy","s_zz","s_xy","s_yz","s_xz" };
static const char* s[6] = { "s_xx","s_yy","s_zz","s_xy","s_yz","s_xz" };
static std::string label;
if (prefix)
@@ -529,13 +530,14 @@ const char* Elasticity::getFieldLabel (size_t i, const char* prefix) const
else
label.clear();
if (i >= (size_t)(nsd*(nsd+1)/2))
if (nsd == 1)
label += "Axial stress";
else if (i+1 >= this->getNoFields())
label += "von Mises stress";
else switch (nsd) {
case 1: label += s1[i]; break;
case 2: label += s2[i]; break;
case 3: label += s3[i]; break;
}
else if (i == 2 && this->getNoFields() == 4)
label += s[3];
else
label += s[i];
return label.c_str();
}
@@ -599,7 +601,10 @@ bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
// Evaluate the finite element stress field
Vector sigmah, sigma;
if (!problem.evalSol(sigmah,fe.dNdX,X)) return false;
if (!problem.evalSol(sigmah,fe.dNdX,X))
return false;
else if (sigmah.size() == 4)
sigmah.erase(sigmah.begin()+2); // Remove the sigma_zz if plane strain
// Integrate the energy norm a(u^h,u^h)
pnorm[0] += sigmah.dot(Cinv*sigmah)*fe.detJxW;
@@ -618,6 +623,9 @@ bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
{
// Evaluate the analytical stress field
sigma = (*anasol)(X);
if (sigma.size() == 4)
sigma.erase(sigma.begin()+2); // Remove the sigma_zz if plane strain
// Integrate the energy norm a(u,u)
pnorm[2] += sigma.dot(Cinv*sigma)*fe.detJxW;
// Integrate the error in energy norm a(u-u^h,u-u^h)

9
src/Integrands/LinIsotropic.C
View File

@@ -134,16 +134,21 @@ bool LinIsotropic::evaluate (Matrix& C, SymmTensor& sigma, double& U,
if (iop > 0)
{
// Calculate the stress tensor, sigma = C*eps
Vector sig; // Use a local variable to avoid a redim of sigma
if (eps.dim() != sigma.dim())
{
// Account for non-matching tensor dimensions
SymmTensor epsil(sigma.dim());
if (!C.multiply(epsil=eps,sigma))
if (!C.multiply(epsil=eps,sig))
return false;
}
else
if (!C.multiply(eps,sigma))
if (!C.multiply(eps,sig))
return false;
sigma = sig; // Add sigma_zz in case of plane strain
if (!planeStress && nsd == 2 && sigma.size() == 4)
sigma(3,3) = nu * (sigma(1,1)+sigma(2,2));
}
if (iop == 3) // Calculate strain energy density, // U = 0.5*sigma:eps

3
src/Integrands/LinIsotropic.h
View File

@@ -33,6 +33,9 @@ public:
//! \brief Empty destructor.
virtual ~LinIsotropic() {}
//! \brief Returns \e false if plane stress in 2D.
virtual bool isPlaneStrain() const { return !planeStress; }
//! \brief Prints out material parameters to the given output stream.
virtual void print(std::ostream&) const;

3
src/Integrands/MaterialBase.h
View File

@@ -36,6 +36,9 @@ public:
//! \brief Empty destructor.
virtual ~Material() {}
//! \brief Returns \e false if plane stress in 2D.
virtual bool isPlaneStrain() const { return true; }
//! \brief Prints out material parameters to the given output stream.
virtual void print(std::ostream&) const {}

161
src/Utility/Tensor.C
View File

@@ -73,6 +73,9 @@ Tensor& Tensor::operator= (const Tensor& T)
for (t_ind i = 1; i <= ndim; i++)
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)
v[2] = T(3,3);
}
return *this;
@@ -326,10 +329,10 @@ std::ostream& SymmTensor::print (std::ostream& os) const
{
switch (n) {
case 1:
return os << v[0] << std::endl;
return os << v.front() << std::endl;
case 2:
return os << v[0] <<'\n'
<< v[2] <<' '<< v[1] << std::endl;
return os << v.front() <<'\n'
<< v.back() <<' '<< v[1] << std::endl;
case 3:
return os << v[0] <<'\n'
<< v[3] <<' '<< v[1] <<'\n'
@@ -340,12 +343,18 @@ std::ostream& SymmTensor::print (std::ostream& os) const
}
bool SymmTensor::redim (const t_ind dim)
SymmTensor::SymmTensor (const t_ind nsd, bool with33) : Tensor(0)
{
if (n == dim) return false;
this->redim(nsd,with33);
}
const_cast<t_ind&>(n) = dim;
v.resize(n*(n+1)/2,real(0));
bool SymmTensor::redim (const t_ind nsd, bool with33)
{
if (n == nsd) return false;
const_cast<t_ind&>(n) = nsd;
v.resize(nsd == 2 && with33 ? 4 : n*(n+1)/2, real(0));
return true;
}
@@ -355,7 +364,7 @@ SymmTensor::SymmTensor (const std::vector<real>& vec) : Tensor(0)
if (vec.size() > 5)
this->redim(3);
else if (vec.size() > 2)
this->redim(2);
this->redim(2, vec.size() == 4);
else if (vec.size() > 0)
this->redim(1);
else
@@ -367,7 +376,7 @@ SymmTensor::SymmTensor (const std::vector<real>& vec) : Tensor(0)
void SymmTensor::copy (const SymmTensor& T)
{
this->redim(T.n);
this->redim(T.n, T.v.size() == 4);
std::copy(T.v.begin(),T.v.end(),v.begin());
}
@@ -384,14 +393,14 @@ SymmTensor& SymmTensor::transform (const Tensor& T)
real S11, S12, S13, S21, S22, S23, S31, S32, S33;
switch (n) {
case 2:
S11 = T(1,1)*v[0] + T(1,2)*v[2];
S12 = T(1,1)*v[2] + T(1,2)*v[1];
S21 = T(2,1)*v[0] + T(2,2)*v[2];
S22 = T(2,1)*v[2] + T(2,2)*v[1];
S11 = T(1,1)*v.front() + T(1,2)*v.back();
S12 = T(1,1)*v.back() + T(1,2)*v[1];
S21 = T(2,1)*v.front() + T(2,2)*v.back();
S22 = T(2,1)*v.back() + T(2,2)*v[1];
v[0] = S11*T(1,1) + S12*T(1,2);
v[1] = S21*T(2,1) + S22*T(2,2);
v[2] = S11*T(2,1) + S12*T(2,2);
v.front() = S11*T(1,1) + S12*T(1,2);
v[1] = S21*T(2,1) + S22*T(2,2);
v.back() = S11*T(2,1) + S12*T(2,2);
break;
case 3:
@@ -419,29 +428,39 @@ SymmTensor& SymmTensor::transform (const Tensor& T)
real SymmTensor::trace () const
{
if (n == 3)
return v[0] + v[1] + v[2];
else if (n == 2)
return v[0] + v[1];
else if (n == 1)
return v[0];
real t = real(0);
return real(0);
if (n == 3)
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;
}
real SymmTensor::det () const
{
if (n == 3)
return v[0]*(v[1]*v[2] - v[4]*v[4])
- v[3]*(v[3]*v[2] - v[5]*v[4])
+ v[5]*(v[3]*v[4] - v[5]*v[1]);
else if (n == 2)
return (v[0]-v[2])*v[1];
else if (n == 1)
return v[0];
real d = real(0);
return real(0);
if (n == 3)
d = v[0]*(v[1]*v[2] - v[4]*v[4])
- v[3]*(v[3]*v[2] - v[5]*v[4])
+ v[5]*(v[3]*v[4] - v[5]*v[1]);
else if (n == 2)
d = (v.front()-v.back())*v[1];
else if (n == 1)
d = v.front();
if (v.size() == 4)
d *= v[2];
return d;
}
@@ -468,14 +487,21 @@ real SymmTensor::inverse (real tol)
}
else if (n == 2)
{
real T11 = v[0];
real T21 = v[2]; real T22 = v[1];
v[0] = T22 / det;
v[1] = T11 / det;
v[2] = -T21 / det;
real T11 = v.front();
real T21 = v.back(); real T22 = v[1];
v.front()= T22 / det;
v[1] = T11 / det;
v.back() = -T21 / det;
if (v.size() == 4)
{
v[0] *= v[2];
v[1] *= v[2];
v[3] *= v[2];
v[2] = (T11*T22 - T21*T21) / det;
}
}
else if (n == 1)
v[0] = real(1) / det;
v.front() = real(1) / det;
return det;
}
@@ -509,27 +535,44 @@ SymmTensor& SymmTensor::rightCauchyGreen (const Tensor& F)
}
/*!
The von Mises value of a symmetric 3D (stress) tensor is defined as follows:
\f[ s_{\rm vm} = \sqrt{
s_{11}(s_{11}-s_{22}) +
s_{22}(s_{22}-s_{33}) +
s_{33}(s_{33}-s_{11}) + 3(s_{12}^2 + s_{23}^2 + s_{31}^2)}
\f]
*/
real SymmTensor::vonMises () const
{
switch (n) {
case 1:
return v[0];
case 2:
return sqrt(v[0]*v[0] + v[1]*v[1] - v[0]*v[1] + 2.0*v[2]*v[2]);
case 3:
return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2] -
v[0]*v[1] - v[1]*v[2] - v[2]*v[0] +
3.0*(v[3]*v[3] + v[4]*v[4] + v[5]*v[5]));
double vms = 0.0;
if (n == 3)
{
double s11 = v[0];
double s22 = v[1];
double s33 = v[2];
double s12 = v[3];
double s23 = v[4];
double s31 = v[5];
vms = sqrt(s11*(s11-s22) + s22*(s22-s33) + s33*(s33-s11) +
3.0*(s12*s12 + s23*s23 + s31*s31));
}
else if (n == 2)
{
double s11 = v[0];
double s22 = v[1];
double s33 = v.size() == 4 ? v[2] : 0.0;
double s12 = v.size() == 4 ? v[3] : v[2];
vms = sqrt(s11*(s11-s22) + s22*(s22-s33) + s33*(s33-s11) + 3.0*s12*s12);
}
else if (n == 1)
vms = v.front();
return real(0);
return vms;
}
/*!
\brief Multiplication between a scalar and a symmetric tensor.
*/
SymmTensor operator* (real a, const SymmTensor& T)
{
SymmTensor S(T.dim());
@@ -541,10 +584,6 @@ SymmTensor operator* (real a, const SymmTensor& T)
}
/*!
\brief Adding a scaled unit tensor to a symmetric tensor.
*/
SymmTensor operator+ (const SymmTensor& T, real a)
{
SymmTensor S(T);
@@ -552,14 +591,13 @@ SymmTensor operator+ (const SymmTensor& T, real a)
for (Tensor::t_ind i = 0; i < S.n; i++)
S.v[i] += a;
if (S.v.size() == 4)
S.v[2] += a;
return S;
}
/*!
\brief Subtracting a scaled unit tensor from a symmetric tensor.
*/
SymmTensor operator- (const SymmTensor& T, real a)
{
SymmTensor S(T);
@@ -567,6 +605,9 @@ SymmTensor operator- (const SymmTensor& T, real a)
for (Tensor::t_ind i = 0; i < S.n; i++)
S.v[i] -= a;
if (S.v.size() == 4)
S.v[2] -= a;
return S;
}

25
src/Utility/Tensor.h
View File

@@ -89,6 +89,9 @@ public:
//! \brief Returns the dimension of this tensor.
t_ind dim() const { return n; }
//! \brief Returns the size of this tensor.
size_t size() const { return v.size(); }
//! \brief Query whether this tensor is symmetric or not.
bool symmetric() const { return v.size() == (size_t)(n*(n+1)/2); }
@@ -130,23 +133,26 @@ public:
class SymmTensor : public Tensor
{
//! \brief Sets the number of spatial dimensions for the tensor to \a newDim.
//! \brief Resets the number of spatial dimensions of the tensor.
//! \param[in] nsd The new tensor dimension
//! \param[in] with33 If \e true and \a nsd = 2, include the 33 term also
//! \return \e true if the dimension was changed, otherwise \e false
//!
//! \details This method is private because the tensor dimension is not
//! supposed to be changed by the application. It is only for internal use.
bool redim(const t_ind newDim);
bool redim(const t_ind nsd, bool with33 = false);
protected:
//! \brief Returns a 0-based array index for the given tensor indices.
//! \details Symmetric tensors are assumed stored with the following order:
//! s11, s22, s33, s12, s23, s13.
//! \details Symmetric 3D tensors are assumed stored with the following order:
//! s11, s22, s33, s12, s23, s13. Symmetric 2D tensors have the order
//! s11, s22, s12 and if the 33 component is included: s11, s22, s33, s12.
virtual t_ind index(t_ind i, t_ind j) const
{
if (i == j)
return i-1; // diagonal term
else if (n == 2)
return 2; // off-diagonal term (2D)
return v.size()-1; // off-diagonal term (2D)
if (i == j+1 || i+2 == j) std::swap(i,j);
return i+2; // upper triangular term (3D)
@@ -157,7 +163,14 @@ protected:
public:
//! \brief Constructor creating a zero tensor.
SymmTensor(const t_ind nsd) : Tensor(nsd) { v.resize(n*(n+1)/2,real(0)); }
//! \param[in] nsd The tensor dimension (1, 2 or 3)
//! \param[in] with33 If \e true and \a nsd = 2, include the 33 term also
//!
//! \details The combination \a nsd = 2 and \a with33 = \e true results in a
//! 2D tensor with the 33-term as the forth component. This is typically used
//! to represent the symmetric stress tensor in 2D plane strain models,
//! where the \f$\sigma_{zz}\f$ component is nonzero.
SymmTensor(const t_ind nsd, bool with33 = false);
//! \brief Constructor creating a symmetric tensor from a vector.
SymmTensor(const std::vector<real>& vec);
//! \brief Copy constructor.