[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" ];
}
}
}