[Getdp] unexpected "segmentation violation" at Pre-processing

gdemesy at physics.utoronto.ca gdemesy at physics.utoronto.ca
Tue Feb 9 15:35:22 CET 2010


Dear Getdp List,

I am trying to calculate the field diffracted by a grating enlighten  
by a plane wave in the visible range (for 1 polarization case = 2D pb,  
in harmonic regime). In short, it amounts to solve a scalar Helmholtz  
equation together with Bloch conditions, PMLs and an appropriate  
source term located inside the diffractive element (just like we did  
in www.opticsinfobase.org/oe/abstract.cfm?uri=oe-15-26-18089).

I have been struggling for a while on the rather simple attached  
problem. This gmsh/getdp code is a "benchmark problem" for which I  
already have the solution via another FEM code. I am stuck with this  
"Segmentation violation (invalid memory reference)" at Pre-processing  
step.

I have unsuccessfully tried to simplify each step of the model as much  
as possible:
- a geometry generated whithout using "Extrude"
- a rough mesh
- Neumann instead of Bloch conditions
- source = 0
...
My conclusion: I'm missing something elementary somewhere.
Any hint would be greatly appreciated...

Thanks a lot for your time.

Best,
Guillaume Demésy,
Postdoc at University of Toronto,
Sajeev John's Group

----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.

-------------- next part --------------
n_rode     = 10;
n_air      = 1; 
paramaille = 3;
lambda0    = 700e-9;
lc_rode    = lambda0/(paramaille*n_rode);
lc_air     = lambda0/(paramaille*n_air);
d          = 4e-6;
Point(1) = {-d/2, -14e-006, 0, lc_air};
Point(2) = {-d/2, -10.8e-006, 0, lc_air};
Point(3) = {-d/2, -3.1e-006, 0, lc_air};
Point(4) = {-d/2, -1e-007, 0, lc_air};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};

out[]= Extrude{d,0,0}{Line{1};Line{2};Line{3};Layers{d/lc_air};};
Delete{
	 Surface{7,11,15};
	}

Point(11) = {    0,   -7e-6 , 0,lc_rode};
Point(12) = { 5e-7,   -7e-6 , 0,lc_rode};
Point(13) = {-5e-7,   -7e-6 , 0,lc_rode};

Circle(20) = {12, 11, 13};
Circle(21) = {13, 11, 12};

Line Loop(24) = { 1,  6,  -4, -5};           Plane Surface(24) = {24};
Line Loop(27) = { 6,  8, -10, -2, -20, -21}; Plane Surface(27) = {27};
Line Loop(28) = { 3, 14, -12, -10};          Plane Surface(28) = {28};
Line Loop(29) = {21, 20};                    Plane Surface(29) = {29};

Physical Line(101) = {1, 2, 3};       // Bloch_Left
Physical Line(102) = {4, 8, 12};	  // Bloch_Right
Physical Line(110) = {5, 14};		  // Dirichlet
Physical Line(120) = {6, 21, 20, 10}; // Continuity

Physical Surface(1000) = {24};        // PML_bot
Physical Surface(2000) = {27};        // Air
Physical Surface(3000) = {28};        // PML_top
Physical Surface(4000) = {29};        // Rode
-------------- next part --------------
Group {
        // Domains
		Air     = Region[2000];
		Rode    = Region[4000];
		PMLtop  = Region[3000];
		PMLbot  = Region[1000];
        
        // Boundaries
		SurfBlochLeft   = Region[101];
		SurfBlochRight  = Region[102];
		SurfDirichlet   = Region[110];
		SurfContinuity  = Region[120];
        
        Omega = Region[{Air,Rode,PMLtop,PMLbot}];
        // Gamma = Region[{SurfBlochLeft,SurfDirichlet,SurfBlochRight}];
}

Function{
        // Geometric and Physical parameters
		period           = 4.e-6;
		mu0              = 4.e-7*Pi;
		epsilon0         = 8.854187817e-12;
        cel              = 299792458.;
        
        // Incident plane wave = So-called Ancillary Problem
        lambda0          = 700.e-3;
        Freq             = cel/lambda0;
		k0               = 2.*Pi/lambda0;
		theta0           = 0. *Pi/180.;
        Ae               = 1.;
		alpha0           = k0*Sin[theta0];
		beta0            = k0*Cos[theta0];
		ezi[]            = Ae * Complex[ Sin[alpha0*X[]+beta0*Y[]] , -Cos[alpha0*X[]+beta0*Y[]] ];
        deph[]           =      Complex[ Sin[alpha0*period]        , -Cos[alpha0*period]        ];
        
        // Optical properties
        epsilon_air      = 1.; 
        epsilon_rode     = 3.; 
        mu_air           = 1.; 
        mu_rode          = 1.; 
		epsilon[Air]     = epsilon_air;
		epsilon[Rode]    = epsilon_rode;
		epsilon1[Air]    = epsilon_air;
		epsilon1[Rode]   = epsilon_air;
		mu[Air]          = mu_air;
		mu[Rode]         = mu_rode;
        
        // PML
        a_pml            =  1.;
        b_pml            = -1.;
        sx               =  1.;
        sy[]             = Complex[a_pml,b_pml]; 
        sz               = 1.;
        epsilon[PMLtop]  = epsilon_air * TensorDiag[sz*sy[]/sx,sx*sz/sy[],sx*sy[]/sz]; 
        epsilon[PMLbot]  = epsilon_air * TensorDiag[sz*sy[]/sx,sx*sz/sy[],sx*sy[]/sz];
        mu[PMLtop]       = mu_air * TensorDiag[sz*sy[]/sx,sx*sz/sy[],sx*sy[]/sz];
        mu[PMLbot]       = mu_air * TensorDiag[sz*sy[]/sx,sx*sz/sy[],sx*sy[]/sz];
        }

Constraint {
		{Name Dirichlet; Type Assign;
			Case {
				    { Region SurfDirichlet; Value 0.; }
	             }
		}
		{Name Bloch;
			Case {
                    { Region SurfBlochRight; Type LinkCplx ; RegionRef SurfBlochLeft;
                      Coefficient deph[]; Function Vector[$X+period,$Y,$Z] ;
                    }
			     }
		}
}
Include "Jacobian_Lib.pro"
Include "Integration_Lib.pro"
Include "prop_TE.pro"
-------------- next part --------------
Integration {
  { Name Int_1 ;
    Case { 
      { Type Gauss ;
        Case { 
	  { GeoElement Point       ; NumberOfPoints  1 ; }
	  { GeoElement Line        ; NumberOfPoints  4 ; }
	  { GeoElement Triangle    ; NumberOfPoints  6 ; }
	  { GeoElement Quadrangle  ; NumberOfPoints  6 ; }
	  { GeoElement Tetrahedron ; NumberOfPoints  6 ; }
	  { GeoElement Hexahedron  ; NumberOfPoints  8 ; }
	  { GeoElement Prism       ; NumberOfPoints  8 ; } 
	}
      }
    }
  }
}
-------------- next part --------------
Jacobian {
  { Name JVol ;
    Case { 
      { Region All ; Jacobian Vol ; }
    }
  }
  { Name JSur ;
    Case {
      { Region All ; Jacobian Sur ; }
    }
  }
}
-------------- next part --------------
FunctionSpace {
  { Name Hgrad; Type Form0;
    BasisFunction {
      { Name sn; NameOfCoef un; Function BF_Node; 
        Support Region[All]; Entity NodesOf[All]; }
      //{ Name sn2; NameOfCoef un2; Function BF_Node_2E; 
      //  Support Region[All]; Entity EdgesOf[All]; }
    }
    Constraint {
      { NameOfCoef un; EntityType NodesOf ; NameOfConstraint Dirichlet; }
      { NameOfCoef un; EntityType NodesOf ; NameOfConstraint Bloch; }
   }
  }
}

// Formulation 1 (BIDON)
/*
 Formulation {
   { Name helmoltz_scalar; Type FemEquation;
     Quantity {
       { Name u; Type Local; NameOfSpace Hgrad; }
     }
     Equation {
       Galerkin { [ k0^2*Dof{u} , {u} ];
                  In Omega; Integration Int_1; Jacobian JVol;  }
       Galerkin { [ -Dof{Grad u} , {Grad u} ];
                  In Omega; Integration Int_1; Jacobian JVol;  }
        }
    }
 }
 */


// Formulation 2 = sources + PML
Formulation {{Name helmoltz_scalar; Type FemEquation;
    		Quantity {{ Name u; Type Local; NameOfSpace Hgrad;}}    
		Equation {
	      		   Galerkin { [-1/mu[]*Dof{Grad u} , {Grad u}];                   // AIR
                 		In Air; Jacobian JVol; Integration Int_1;  }
	      		   
                   Galerkin { [k0^2*epsilon[]*Dof{u} , {u}]; 
                 		In Air; Jacobian JVol; Integration Int_1;  }
  	      		   
                   Galerkin { [-1/mu[]*Dof{Grad u} , {Grad u}];                   // RODE
                 		In Rode; Jacobian JVol; Integration Int_1;  }
  	      		   
                   Galerkin { [epsilon[]*k0^2*ezi[] , {u}];
                 		In Rode; Jacobian JVol; Integration Int_1;  }

                   Galerkin { [-epsilon1[]*k0^2*ezi[], {u}];
                 		In Rode; Jacobian JVol; Integration Int_1;  }
                   
                   Galerkin { [k0^2*epsilon[]*Dof{u} , {u}];
                 		In Rode; Jacobian JVol; Integration Int_1;  }
  	      		   
                   Galerkin { [-1/mu[]*Dof{Grad u},{Grad u}];                    // PMLtop
                 		In PMLtop; Jacobian JVol; Integration Int_1;  }
  	      		   
                   Galerkin {[k0^2*CompZZ[epsilon[]]*Dof{u} , {u}];
                 		In PMLtop; Jacobian JVol; Integration Int_1;  }
  	      		   
                   Galerkin { [-1/mu[]*Dof{Grad u}, {Grad u}];                   // PMLbot
                 		In PMLbot; Jacobian JVol; Integration Int_1;  }
  	      		   
                   Galerkin {[k0^2*epsilon[]*Dof{u} , {u}]; 
                 		In PMLbot; Jacobian JVol; Integration Int_1;  }
                }
            }
        }


Resolution {
  { Name helmoltz_scalar;
    System {
      { Name S; NameOfFormulation helmoltz_scalar; Type ComplexValue; Frequency Freq;}
    }
    Operation {
      Generate[S] ;
      // Print[S];
      Solve[S] ;
      SaveSolution[S] ;
    }
  }
}

PostProcessing {
    { Name get_Ed; NameOfFormulation helmoltz_scalar;
            Quantity {
                { Name e_d; Value { Local { [ {u} ]; In Omega; } } }
}}}

PostOperation {
{ Name Ed; NameOfPostProcessing get_Ed ;
Operation {
Print [ e_d, OnElementsOf Omega, File "map_Ez_diffacted.pos" ];
}
}
}