Using @OrthotropicBehaviour<Pipe> and ComputeStiffnessTensor<>

[lmersel]

Hello,

I am using mfront to define an new energy decomposition for a phase fracture modeling, and for any anisotropic elastic materials (ref Nguyen et al. ) . It is similar of what Jeremy Bleyer did with Amor and Miehe decomposition coupling with fenics. At the end, I would like to compare the three formulation in dynamic modelling.

The first step of the law is to provide as an input an anisotropic ( or at least orthotopic) stiffness matrix to mfront (the variable C). Then I have to compute the eigenvalues and and eigenprojectors of C. At the end, I should transform the stiffness matrix into C^(1/2) and C^(-1/2) like in the joined image.

The matrix C is considered constant so the decomposition should carry out only once.
I started to code this step with 3 blocks :
@OrthotropicBehaviour<Pipe>;
@ComputeStiffnessTensor<UnAltered> {};
@InitLocalVariables{}

At the compilation I got this error :
In file included from PhaseFieldDisplacementNguyenSplit-generic.cxx:27: PhaseFieldDisplacementNguyenSplit.mfront: In instantiation of 'void tfel::material::PhaseFieldDisplacementNguyenSplit<hypothesis, Numeric Type, false>::initialize() [with tfel::material::ModellingHypothesis::Hypothesis hypothesis = tfel::material::ModellingHypothesis::AXISYM METRICALGENERALISEDPLANESTRAIN; NumericType = double]': /stck/lmersel/CODE_PROGRAMME/MFront_MGIS/TFEL-4.0/install/include/MFront/GenericBehaviour/Integrate.hxx:318:19: required from 'int mfro nt::gb::integrate(mfront_gb_BehaviourData&, typename Behaviour::SMFlag, tfel::material::OutOfBoundsPolicy) [with Behaviour = tfel::materi al::PhaseFieldDisplacementNguyenSplit<tfel::material::ModellingHypothesis::AXISYMMETRICALGENERALISEDPLANESTRAIN, double, false>; mfront_g b_BehaviourData = mfront_gb_BehaviourData; typename Behaviour::SMFlag = tfel::material::TangentOperatorTraits<tfel::material::MechanicalB ehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR>::SMFlag]' PhaseFieldDisplacementNguyenSplit-generic.cxx:618:145: required from here PhaseFieldDisplacementNguyenSplit.mfront:21:32: error: 'using StiffnessTensor = tfel::config::Types<1, double, false>::StiffnessTensor' { aka 'struct tfel::math::st2tost2<1, double>'} has no member named 'computeEigenVectors' compilation terminated due to -Wfatal-errors. make: *** [PhaseFieldDisplacementNguyenSplit-generic.o] Error 1 callMake: can't build target 'all' libraries building went wrong
I am new in mfront so I propably did something wrong. I attached the .mfront file (in .txt format).

Thank you very much for the help.

Kind regards,

Lamia
PhaseFieldDisplacementNguyenSplit.txt

[thelfer]

Hi Lamia, Unfortunately, there is no function to compute a eigenvalues of a fourth order tensor so far in MFront. But, in the orthotropic case, I think that we can use a small trick. Let me think of it.

[lmersel]

With some help, I finally updated the code in fixing some errors.
PhaseFieldDisplacementNguyenSplit.txt

The new issue is about the product between the tensor 'Cp1_2' and the strain, thats logical.
PhaseFieldDisplacementNguyenSplit.mfront:143:32: error: no match for 'operator*' (operand types are 'std::tuple_element<0, std::tuple<tfel::math::tmatrix<3, 3, tfel::math::stensor<1, double> >, tfel::math::tmatrix<3, 3, tfel::math::stensor<1, double> > > >::type' {aka 'tfel::math::tmatrix<3, 3, tfel::math::stensor<1, double> >'} and 'const StrainStensor' {aka 'const tfel::math::stensor<1, double>'}) compilation terminated due to -Wfatal-errors. make: *** [PhaseFieldDisplacementNguyenSplit-generic.o] Error 1
I keep fixing the code :)

Best,

Lamia

[thelfer]

Dear Lamia,

The returned value from the diagonalize function was wrong, causing the issue.
See this corrected file: PhaseFieldDisplacementNguyenSplit.zip

However, I am pretty sure that you made some mistake in the computation of the square root (sqrt(x)) of the stiffness stensor C because you don't take the square root of the eigen values.

Best,

Thomas

[lmersel]

Dear Thomas,

Thank you for your answer. Indeed, you are right. I made some mistakes and corrected them. Thanks !
For the negative square root of the stiffness, I used power<-1,2>(x).
I have a small question in the calculation of the positive and the negative square root of C. For the case N=1, can i express the tensor like that
if constexpr (N == 1) { // pos and neg square root of stiffness tensor C auto Cp1_2 = sqrt(vp(0)) * (ns[0] ^ ns[0]) ; auto Cn1_2 = 1/sqrt(vp(0)) * (ns[0] ^ ns[0]); return {vp, ns, Cp1_2, Cn1_2};
Do I need to add sqrt(vp(1)) * (ns[1] ^ ns[1]) + sqrt(vp(2))*(ns[2] ^ ns[2]);, since it is 1D ?
Same question for the 2D case, since the case can be reduced like

(in the kelvon basis E1 =e1^e1, E2 = e2^e2, E3 = (e1^e2 + e2 ^e1)/2
and his spetral decomposition

(ai : the eigen vector)

Thank for your help,

Best,

Lamia

[thelfer]

Do I need to add sqrt(vp(1)) * (ns[1] ^ ns[1]) + sqrt(vp(2))*(ns[2] ^ ns[2]);, since it is 1D ?

Yes. There are three eigenvalues.

To be honest, I am now pretty sure that this code may not work in 1D or 2D because there is no reason that the Kelvin modes respects the "space" symmetry, but the ^ operator forces this (for instance, the result of this operator is forced to be diagonal in 1D).

[lmersel]

Thank for your remark. Indeed I am agreed. I am thinking about setting aside the 1D case and focusing on 2D et 3D, since I will never work in 1D. Otherwise, I could compute separatly the 1D case.

If I use the plane strain assumption (N=2), in small strain, I should provide a vector epsilon of size 4x1( e = [e00, e11, e22, sqrt(2)*e01]^T ). Thus the transformed strain, 'et' that I want to compute is the product of et = C^(1/2) : e with C^(1/2) a 4x4 matrix.
C = [ [C00, C01, C02, sqrt(2)*C03],
[C01, C11, C12, sqrt(2)*C13],
[C02, C12, C22, sqrt(2)*C23],
[ sqrt(2)*C03, sqrt(2)*C13, sqrt(2)C23, 2C33]]

In 2D or 3D, the operator '^' used shouldn't provide a diagonal matrix right?
.
In the article I read, they use a generic formula called the Sylvester's formula to compute the eigenprojector wi^wi.

Thanks,

Lamia