scattererTmatrix.pro

///////////////////////////////
// Author : Guillaume Demesy //
// scattererTmatrix.pro      //
///////////////////////////////

Include "scattererTmatrix_data.geo";
myDir = "resT/";

Group {
  Scat_In        = Region[1];
  Scat_Out       = Region[2];
  Domain         = Region[{Scat_In,Scat_Out}];
  PMLs           = Region[3];
  All_domains    = Region[{Scat_In,Scat_Out,PMLs}];
  SurfInt        = Region[{10}];
  SurfDirichlet  = Region[{20}];
}

Function{
  I[]         = Complex[0,1];
  avoid_sing  = 1.e-14;
  mu0         = 4*Pi*100.0*nm;
  epsilon0    = 8.854187817e-3*nm;
  cel         = 1.0/Sqrt[epsilon0 * mu0];
  Freq        = cel/lambda0;
  omega0      = 2.0*Pi*Freq;
  k_Out       = 2.0*Pi*Sqrt[epsr_Out_re]/lambda0;
  R3D[]   = Sqrt[X[]^2+Y[]^2+Z[]^2];
  R2D[]   = Sqrt[X[]^2+Y[]^2];
  cos_theta[] = Z[]/(R3D[]+avoid_sing);
  theta[]     = Acos[cos_theta[]];
  phi[]       = Atan2[-Y[],-X[]]+Pi;
  cart2sph[]  =   Tensor[X[]/R3D[], Y[]/R3D[], Z[]/R3D[],       
                         X[]*Z[]/(R3D[]*R2D[]), Y[]*Z[]/(R3D[]*R2D[]), -(X[]^2+Y[]^2)/(R3D[]*R2D[]),
                        -Y[]/R2D[], X[]/R2D[],0];
  R[] = Vector[X[]/R3D[], Y[]/R3D[], Z[]/R3D[]];

  a_pml =  1.;
  b_pml = -siwt;

  sr[]         = Complex[a_pml,b_pml];
  sphi[]       = Complex[1,0];
  stheta[]     = Complex[1,0];
  r_tild[]     = r_pml_in + (R3D[] - r_pml_in) * sr[];
  theta_tild[] = theta[];
  pml_tens_temp1[]  = TensorDiag[(r_tild[]/R3D[])^2 * sphi[]*stheta[]/sr[],
                                 (r_tild[]/R3D[])   * sr[]*stheta[]/sphi[],
                                 (r_tild[]/R3D[])   * sphi[]*sr[]/stheta[]];
  pml_tens_temp2[]  = Rotate[pml_tens_temp1[],0,-theta[]-Pi/2,0];
  pml_tens[]        = Rotate[pml_tens_temp2[],0,0,-phi[]];

  epsilonr_In[]       = Complex[epsr_In_re  , epsr_In_im];
  epsilonr_Out[]      = Complex[epsr_Out_re , epsr_Out_im];

  epsilonr[Scat_In]   = epsilonr_In[]  * TensorDiag[1.,1.,1.];
  epsilonr[Scat_Out]  = epsilonr_Out[] * TensorDiag[1.,1.,1.];
  epsilonr[PMLs]      = epsilonr_Out[] * pml_tens[];
                          
  epsilonr1[Scat_In]  = epsilonr_Out[] * TensorDiag[1.,1.,1.];
  epsilonr1[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];                
  epsilonr1[PMLs]     = epsilonr_Out[] * pml_tens[];
                          
  mur[Scat_In]        = TensorDiag[1.,1.,1.];
  mur[Scat_Out]       = TensorDiag[1.,1.,1.];
  mur[PMLs]           = pml_tens[];

  For pe In {1:p_max}
    ne = Floor[Sqrt[pe]];
    me = ne*(ne+1) - Floor[pe];
    Mnm_source~{pe}[] = Mnm[1,ne,me,XYZ[],k_Out];
    Nnm_source~{pe}[] = Nnm[1,ne,me,XYZ[],k_Out];
    Mnm_out~{pe}[]    = Mnm[3,ne,me,XYZ[],k_Out];
    Nnm_out~{pe}[]    = Nnm[3,ne,me,XYZ[],k_Out];
    Xnm~{pe}[]        = Xnm[ne,me,X[],Y[],Z[]];
    Ynm~{pe}[]        = Ynm[ne,me,X[],Y[],Z[]];
    Znm~{pe}[]        = Znm[ne,me,X[],Y[],Z[]];
    Ynm_r~{pe}[]      = Ynm[ne,me,X[],Y[],Z[]] * R[];
    source_M~{pe}[]   = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Mnm_source~{pe}[];
    source_N~{pe}[]   = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm_source~{pe}[];
    SphHankelOutgoing_n~{pe}[]   = (JnSph[ne  ,k_Out*r_pml_in]- siwt * I[]*YnSph[ne  ,k_Out*r_pml_in]) ;
    SphHankelOutgoing_nm1~{pe}[] = (JnSph[ne-1,k_Out*r_pml_in]- siwt * I[]*YnSph[ne-1,k_Out*r_pml_in]) ;
    dRicattiBessel~{pe}[]        = (k_Out * r_pml_in * (SphHankelOutgoing_nm1~{pe}[]-((ne+1)/((k_Out*r_pml_in))) * SphHankelOutgoing_n~{pe}[]) + SphHankelOutgoing_n~{pe}[]);
    normalize_fhnm_X~{pe}[]      = 1/SphHankelOutgoing_n~{pe}[];
    normalize_fenm_Z~{pe}[]      = k_Out*r_pml_in/dRicattiBessel~{pe}[];
    normalize_fenm_Y~{pe}[]      = k_Out*r_pml_in/(SphHankelOutgoing_n~{pe}[]*Sqrt[ne*(ne+1)]);
  EndFor  
}

Constraint {
  {Name Dirichlet; Type Assign;
   Case {
     { Region SurfDirichlet; Value 0.; }
   }
  }
}

Jacobian {
  { Name JVol ;
    Case {
      { Region All ; Jacobian Vol ; }
    }
  }
  { Name JSur ;
    Case {
      { Region All ; Jacobian Sur ; }
    }
  }
}

Integration {
  { Name Int_1 ;
    Case { 
      { Type Gauss ;
        Case { 
          { GeoElement Point        ; NumberOfPoints   1 ; }
          { GeoElement Line2        ; NumberOfPoints   4 ; }
          { GeoElement Triangle     ; NumberOfPoints  12 ; }
          { GeoElement Triangle2    ; NumberOfPoints  12 ; }
          { GeoElement Tetrahedron  ; NumberOfPoints  17 ; }
          { GeoElement Tetrahedron2 ; NumberOfPoints  17 ; }
        }
      }
    }
  }
}

FunctionSpace {
  { Name Hcurl; Type Form1;
    BasisFunction {
      { Name sn; NameOfCoef un; Function BF_Edge;
        Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
      { Name sn2; NameOfCoef un2; Function BF_Edge_2E;
        Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
        If (is_FEM_o2==1)
          { Name sn3; NameOfCoef un3; Function BF_Edge_3F_b;
            Support Region[{All_domains,SurfInt}]; Entity FacetsOf[All]; }
          { Name sn4; NameOfCoef un4; Function BF_Edge_3F_c;
            Support Region[{All_domains,SurfInt}]; Entity FacetsOf[All]; }
          { Name sn5; NameOfCoef un5; Function BF_Edge_4E;
            Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
         EndIf
    }
    Constraint {
      { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint Dirichlet; }
      { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
      If (is_FEM_o2==1)
        { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
        { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint Dirichlet; }  
        { NameOfCoef un5; EntityType EdgesOf  ; NameOfConstraint Dirichlet; }  
      EndIf
    }
  }
}

Formulation {
  {Name VPWMN_helmholtz_vector; Type FemEquation;
   Quantity {
     { Name u; Type Local; NameOfSpace Hcurl;}
    }
    Equation {
      Galerkin { [-1/mur[]*Dof{Curl u}             , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1;  }
      Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u}     ]; In All_domains; Jacobian JVol; Integration Int_1;  }
      Galerkin { [ (omega0/cel)^2*(epsilonr[]-epsilonr1[])*
                      ($isN ? Nnm[1,$NE,$ME,XYZ[],k_Out] : Mnm[1,$NE,$ME,XYZ[],k_Out]) 
                                                    , {u} ]; In Scat_In; Jacobian JVol; Integration Int_1;}      
    }
  }
}

Resolution {
  { Name res_VPWall_helmholtz_vector;
    System {
      { Name T; NameOfFormulation VPWMN_helmholtz_vector; Type ComplexValue; }
  }
  Operation {
    CreateDir[Str[myDir]];        
    For pe In {1:p_max}
      Evaluate[ $isN = 0 ];
      Evaluate[ $PE  = pe ];
      Evaluate[ $NE  = Floor[Sqrt[$PE]] ];
      Evaluate[ $ME  = $NE*($NE+1) - Floor[$PE] ];
      If (pe==1)
        Generate[T];
        Solve[T];
      EndIf
      GenerateRHS[T];
      SolveAgain[T];
      PostOperation[VPWM_postop~{pe}];
      Evaluate[$isN=1];
      GenerateRHS[T];
      SolveAgain[T];
      PostOperation[VPWN_postop~{pe}];
    EndFor
    }
  }
}

PostProcessing {
  For pe In {1:p_max}
    { Name VPWM_postpro~{pe}; NameOfFormulation VPWMN_helmholtz_vector; NameOfSystem T;
      Quantity {
        { Name E_scat           ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
        { Name H_scat           ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
        { Name Mnm_source~{pe}  ; Value { Local { [   Mnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }

        For po In {1:p_max}
          { Name feM~{pe}~{po} ; Value { Integral { [ (1/r_pml_in)^2 *normalize_fenm_Z~{po}[]*  {u}*Conj[Znm~{po}[]]  ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
          { Name fhM~{pe}~{po} ; Value { Integral { [ (1/r_pml_in)^2 *normalize_fhnm_X~{po}[]*  {u}*Conj[Xnm~{po}[]]  ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
        EndFor
      }
    }
    { Name VPWN_postpro~{pe}; NameOfFormulation VPWMN_helmholtz_vector; NameOfSystem T;
      Quantity {
        { Name E_scat            ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
        { Name H_scat            ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
        { Name Nnm_source~{pe}   ; Value {    Local { [   Nnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }
        For po In {1:p_max}
          { Name feN~{pe}~{po}  ; Value { Integral { [ 1/(r_pml_in^2*epsilonr_Out[]) *normalize_fenm_Z~{po}[]*  {u}*Conj[Znm~{po}[]]  ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
          { Name fhN~{pe}~{po}  ; Value { Integral { [ 1/(r_pml_in^2*epsilonr_Out[]) *normalize_fhnm_X~{po}[]*  {u}*Conj[Xnm~{po}[]]  ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
        EndFor
      }
    }
  EndFor
}

PostOperation {
  For pe In {1:p_max}
    {Name VPWM_postop~{pe}; NameOfPostProcessing VPWM_postpro~{pe} ;
      Operation {
        If (flag_plotcuts==1)
          Print [ E_scat , OnGrid
                {(r_pml_in-1*nm)*Sin[$B]*Cos[$C], (r_pml_in-1*nm)*Sin[$B]*Sin[$C], (r_pml_in-1*nm)*Cos[$B]} 
                {(r_pml_in-1*nm), {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)}, 
                {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} }, 
                File StrCat[myDir,StrCat["E_scat_onsphere_cart_M",Sprintf["%g.pos",pe]]],
                Name StrCat["E_scat_onsphere_cart_M",Sprintf["%g",pe]]];
        EndIf
        For po In {1:p_max}
          Print[feM~{pe}~{po}[SurfInt] , OnGlobal, Format Table, File StrCat[myDir,StrCat[StrCat["feM_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
          Print[fhM~{pe}~{po}[SurfInt] , OnGlobal, Format Table, File StrCat[myDir,StrCat[StrCat["fhM_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
        EndFor
      }
    }
    {Name VPWN_postop~{pe}; NameOfPostProcessing VPWN_postpro~{pe} ;
      Operation {
        If (flag_plotcuts==1)
          Print [ E_scat , OnGrid
                {(r_pml_in-1*nm)*Sin[$B]*Cos[$C], (r_pml_in-1*nm)*Sin[$B]*Sin[$C], (r_pml_in-1*nm)*Cos[$B]} 
                {(r_pml_in-1*nm), {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)}, 
                {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} }, 
                File StrCat[myDir,StrCat["E_scat_onsphere_cart_N",Sprintf["%g.pos",pe]]],
                Name StrCat["E_scat_onsphere_cart_N",Sprintf["%g",pe]]];
        EndIf
        For po In {1:p_max}
          Print[feN~{pe}~{po}[SurfInt] , OnGlobal, Format Table, File StrCat[myDir,StrCat[StrCat["feN_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
          Print[fhN~{pe}~{po}[SurfInt] , OnGlobal, Format Table, File StrCat[myDir,StrCat[StrCat["fhN_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
        EndFor      
      }
    }
  EndFor
}

DefineConstant[
  R_ = {"res_VPWall_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
  C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
  P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];