[Getdp] Implementation of PML for linear elastic waves

[Julien de Rosny]

Hello,

I am trying to implement PML for elastic waves (but without success up 
to now). I am starting with the classical coordinate transformation 
introduced by Berenger

Because only D1 and D2 operator are implemented in getdp, I tried (but I 
am not a specialist in finite elements) to use also curl operator and 
perform a  transformation to express the PML with the D1 D2 curl operators

To that end I changed the Formulation proposed in the getdep 
Elasticity_3D.pro example. See bellow.
A compressional wave is absorbed at the PML, but not a transverse one 
that only involves c66.

To workout the C_m_xx matrix, I used Maxima to express  the stiffness 
tensor in the D1 D2 and curl basis.

Any idea why it is is not properly working?

Thanks.

Regards.

Julien de Rosny
Insitut Langevin
1 rue Jussieu
75005 PARIS



Formulation {
   { Name Mec3D_u ; Type FemEquation ;
     Quantity {
       { Name u ; Type Local ; NameOfSpace H_u_Mec3D ; }
          }
     Equation {
        Galerkin { [L[]* C_m_11[] * Dof{D1 u} , {D1 u} ] ;
                   In Domain_Disp ; Jacobian Vol ; Integration GradGrad ; }
        Galerkin { [L[]*  C_m_22[] * Dof{D2 u} , {D2 u} ] ;
                   In Domain_Disp ; Jacobian Vol ; Integration GradGrad ; }

        Galerkin { [ L[]*C_m_32[] * Dof{D2 u} , {Curl u} ] ;
                   In Domain_Disp ; Jacobian Vol ; Integration GradGrad ; }
        Galerkin { [ L[]*C_m_23[] * Dof{Curl u} , {D2 u} ] ;
                   In Domain_Disp ; Jacobian Vol ; Integration GradGrad ; }
        Galerkin { [ L[]*C_m_33[] * Dof{Curl u} , {Curl u} ] ;
                   In Domain_Disp ; Jacobian Vol ; Integration GradGrad ; }

        Galerkin { DtDt [ rho[] * Dof{u} , {u} ];
                   In Domain_Disp ; Jacobian Vol ; Integration GradGrad ; }
        Galerkin { [ -F[] , {u} ] ;
                    In Domain_Force; Jacobian Sur; Integration GradGrad; }

                      }
   }
}

Where


     lambda=nu0*E0/((1+nu0)*(1-2*nu0)); //coef lamé
     mu=E0/(2*(1+nu0)); //coef lamé - shear modulus
     K=E0/3/(1-2*nu0);

     c11 = lambda+2*mu; c12 = lambda; c66 = mu;


        DampingProfileX[] = ( (X[]>=Lx) || (X[]<=(-Lx)) )? 
1/(Lx+PmlThickness-Fabs[X[]] + epsilon)-1/(PmlThickness + epsilon) : 0;
     DampingProfileY[] = ( (Y[]>=Ly) || (Y[]<=(-Ly)) )? 
1/(Ly+PmlThickness-Fabs[Y[]] + epsilon)-1/(PmlThickness + epsilon) : 0;
     DampingProfileZ[] = ( (Z[]>=Lz) || (Z[]<=(-Lz)) )? 
1/(Lz+PmlThickness-Fabs[Z[]] + epsilon)-1/(PmlThickness + epsilon) : 0;

     d1[] =Complex[1,  DampingProfileX[]];
     d2[] =Complex[1,  DampingProfileY[]];
     d3[] =Complex[1,  DampingProfileZ[]];

   L[] = d1[]*d2[]*d3[];

   C_m_11[   ] = Tensor[ c11/d1[]/d1[], c12/d1[]/d2[], c12/d1[]/d3[],
                            c12/d1[]/d2[], c11/d2[]/d2[], c12/d2[]/d3[],
                            c12/d1[]/d3[], c12/d2[]/d3[], c11/d3[]/d3[] ];

   C_m_22[  ] = Tensor[ c66*(1/d1[]/d2[]+1/d1[]^2/2+1/d2[]^2/2), 0, 0,
                            0, c66*(1/d2[]/d3[]+1/d3[]^2/2+1/d2[]^2/2), 0,
                            0, 0, c66*(1/d1[]/d3[]+1/d3[]^2/2+1/d1[]^2/2) ];

   C_m_32[  ] = Tensor[ 0, c66*(1/d2[]^2-1/d3[]^2)/2, 0,
                            0,                    0, 
c66*(1/d3[]^2-1/d1[]^2)/2,
                            c66*(1/d1[]^2-1/d2[]^2)/2, 0,       0];

    C_m_23[  ] = Tensor[ 0, 0,                             
c66*(1/d1[]^2-1/d2[]^2)/2,
                            c66*(1/d2[]^2-1/d3[]^2)/2, 0,       0,
                                0              , 
c66*(1/d3[]^2-1/d1[]^2)/2,                  0];

    C_m_33[  ] = Tensor[ c66*(-1/d2[]/d3[]+1/d3[]^2/2+1/d2[]^2/2), 0, 0,
                            0, c66*(-1/d1[]/d3[]+1/d3[]^2/2+1/d1[]^2/2), 0,
                            0, 0, 
c66*(-1/d2[]/d1[]+1/d1[]^2/2+1/d2[]^2/2) ];



If the
       Galerkin { [L[]* C_m_11[] * Dof{D1 u} , {D1 u} ] ;
                   In Domain_Disp ; Jacobian Vol ; Integration GradGrad ; }
        Galerkin { [L[]*  C_m_22[] * Dof{D2 u} , {D2 u} ] ;
                   In Domain_Disp ; Jacobian Vol ; Integration GradGrad ; }



*** MAXIMA CODE

PML:zeromatrix(9,9);PML[1,1]:1/l1;PML[2,2]:1/l2;PML[3,3]:1/l3;PML[4,4]:1/l2;PML[5,5]:1/l3;PML[6,6]:1/l1;PML[7,7]:1/l3;PML[8,8]:1/l1;PML[9,9]:1/l2;

C:matrix([c11,c12,c12,0,0,0,0,0,0],[C12,c11,c12,0,0,0,0,0,0],[C12,c12,c11,0,0,0,0,0,0 
],[0,0,0,c66,0,0,c66,0,0],[0,0,0,0,c66,0,0,c66,0],[0,0,0,0,0,c66,0,0,c66],[0,0,0,c66,0,0,c66,0,0],[0,0,0,0,c66,0,0,c66,0],[0,0,0,0,0,c66,0,0,c66]);
P:zeromatrix(9,9);P[1,1]:1;P[2,2]:1;P[3,3]:1;P[4,6]:1;P[4,9]:1;P[5,4]:1;P[5,7]:1;P[6,5]:1;P[6,8]:1;P[7,4]:1;P[7,7]:-1;P[8,5]:1;P[8,8]:-1;P[9,6]:1;P[9,9]:-1;P;
Pi:invert(P);
Cget:expand(P.PML.C.PML.Pi);
C22:submatrix(1,2,3,7,8,9,Cget,1,2,3,7,8,9);C32:submatrix(1,2,3,4,5,6,Cget,1,2,3,7,8,9);C23:submatrix(1,2,3,7,8,9,Cget,1,2,3,4,5,6);C33:submatrix(1,2,3,4,5,6,Cget,1,2,3,4,5,6);