Commit 214fd0e0 authored by Christophe Geuzaine's avatar Christophe Geuzaine

work-in-progress benchmarks have moved to https://gitlab.onelab.info/doc/sandbox

parent ac3356da
Pipeline #2120 canceled with stage
This diff is collapsed.
Include "cavitySchrod2d_data.geo";
Solver.AutoMesh = 0;
If(MPI_Size == 1) // sequential meshing
start = 0;
end = N_DOM-1;
EndIf
If(MPI_Size > 1) // parallel meshing
start = MPI_Rank;
end = MPI_Rank;
EndIf
For idom In {start:end}
Delete Model;
ct = Cos(theta);
st = Sin(theta);
Point(1) = {idom*dDom*ct, idom*dDom*st, 0., LC};
myExtrudedLine[] = Extrude {-DY*st, DY*ct, 0} {
Point{1}; Layers{DY/LC};
};
myExtrudedSurface[] = Extrude {dDom*ct, dDom*st, 0} {
Line{myExtrudedLine[1]}; Layers{dDom/LC}; Recombine;
};
lateralSides[] = {};
For i In {2:3}
lateralSides += myExtrudedSurface[i];
EndFor
pmlLeft[] = Extrude {-dBb*ct, -dBb*st, 0} {
Line{myExtrudedLine[1]}; Layers{dBb/LC}; Recombine;
};
pmlRight[] = Extrude {dBb*ct, dBb*st, 0} {
Line{myExtrudedSurface[0]}; Layers{dBb/LC}; Recombine;
};
pmlLeftSides[] = {};
For i In {2:3}
pmlLeftSides += pmlLeft[i];
EndFor
pmlRightSides[] = {};
For i In {2:3}
pmlRightSides += pmlRight[i];
EndFor
Physical Line(-((idom+1)*1000+10)) = myExtrudedLine[1]; // left face
Physical Line(((idom+1)*1000+20)) = myExtrudedSurface[0]; // right face
Physical Line(-((idom+1)*1000+202)) = lateralSides[]; // lateral shell
Physical Surface(((idom+1)*1000+200)) = myExtrudedSurface[1];
Physical Line(-((idom+1)*1000+1)) = pmlLeft[0]; // left face
Physical Line(((idom+1)*1000+102)) = pmlLeftSides[]; // lateral shell
Physical Surface(((idom+1)*1000+100)) = pmlLeft[1];
Physical Line(((idom+1)*1000+4)) = pmlRight[0]; // right face
Physical Line(-((idom+1)*1000+302)) = pmlRightSides[]; // lateral shell
Physical Surface(((idom+1)*1000+300)) = pmlRight[1];
If(StrCmp(OnelabAction, "check")) // only mesh if not in onelab check mode
Printf("Meshing waveguide subdomain %g...", idom);
Mesh 2;
CreateDir Str(DIR);
Save StrCat(MSH_NAME, Sprintf("%g.msh", idom));
Printf("Done.");
EndIf
EndFor
BoundingBox {0, DX, 0, DY, 0, DZ};
Include "cavitySchrod2d_data.geo";
DefineConstant[ // allows to set these from outside
// type of walls
WALLS = {1, Name "Input/05Walls",
Choices {0="Transparent", 1="Metallic"}},
// excitation mode
MODE = {2, Name "Input/05m"}, // y
// transmission boundary condition
TC_TYPE = {3, Name "Input/01Transmission condition",
Choices {0="Order 0", /*1="Order 2", 2="Pade (OSRC)",*/ 3="PML"}},
NP_OSRC = 4,
// parameters for the DDM iterative solver
SOLVER = "gmres", // bcgs, gmsh_pcleft, ...
TOL = 1e-6,
MAXIT = 1000,
RESTART = MAXIT
];
Function {
I[] = Complex[0, 1];
// local interface coordinates for the delta functions
For idom In {0:N_DOM-1}
For jdom In {0:1}
P~{idom}~{jdom}[] = (-X[]*Sin[theta]+Y[]*Cos[theta]);
Q~{idom}~{jdom}[] = Z[];
EndFor
EndFor
P[] = (-X[]*Sin[theta]+Y[]*Cos[theta]);
deltaT = 1.;
c[] = 1.;
om[] = 1.;
R[] = Sqrt[(X[]-DX/2.)^2+(Y[]-DY/2.)^2];
V[] = 2*Fabs[R[]]^2;
F[] = (4*I[]/deltaT + 2*V[]);
Kn[] = 4*I[]/deltaT - 2*V[];
k[] = 1.; // will be used to build the TC
Psi_n[] = Exp[-Fabs[R[]]^2]; // from outside
// parameter for ABC
kInf[] = k[];//BETA_M[];
alphaBT[] = 0; //1/(2*R_EXT) - I[]/(8*k*R_EXT^2*(1+I[]/(k*R_EXT)));
betaBT[] = 0; // -1/(2*I[]*k); //- 1/(2*I[]*k*(1+I[]/(k*R_EXT)));
// parameter for 0th order TC
kDtN[] = k[];
// parameters for 2nd order TC
// J.-F. Lee
kmax[] = Pi/LC ;
delt[] = Sqrt[kmax[]^2-k[]^2]/Sqrt[k[]^2];
Coef_Lee1[] = 1/(1 + I[]*delt[]);
Coef_Lee2[] = -Coef_Lee1[];
// OO2 Gander 2002, pp. 46-47
xsimin = 0;
xsimax = Pi / LC;
deltak[] = Pi;
alphastar[] = I[] * ((k[]^2 - xsimin^2) * (k[]^2 - (k[]-deltak[])^2))^(1/4);
betastar[] = ((xsimax^2 - k[]^2) * ((k[]+deltak[])^2 - k[]^2))^(1/4);
a[] = - (alphastar[] * betastar[] - k[]^2) / (alphastar[] + betastar[]);
b[] = - 1 / (alphastar[] + betastar[]);
// parameters for Pade-type TC
kepsI = 0.;
keps[] = k[]*(1+kepsI*I[]);
theta_branch = Pi/4;
// parameters for PML TC are defined piecewise on groups (see below)
}
Group{
For idom In {0:N_DOM-1}
Omega~{idom} = Region[((idom+1)*1000+200)];
GammaD0~{idom} = Region[{}];
GammaD~{idom} = Region[{}];
Pml~{idom}~{0} = Region[{}];
Pml~{idom}~{1} = Region[{}];
PmlD0~{idom}~{0} = Region[{}];
PmlD0~{idom}~{1} = Region[{}];
PmlInf~{idom}~{0} = Region[{}];
PmlInf~{idom}~{1} = Region[{}];
If(WALLS == 1)
GammaD0~{idom} += Region[ ((idom+1)*1000+202) ];
EndIf
GammaInf~{idom} = Region[{}];
If(WALLS == 0)
GammaInf~{idom} += Region[ ((idom+1)*1000+202) ];
EndIf
GammaN~{idom} = Region[{}];
If (idom == 0)
Sigma~{idom}~{0} = Region[{}];
Sigma~{idom}~{1} = Region[ ((idom+1)*1000+20) ];
GammaD0~{idom} += Region[ ((idom+1)*1000+10) ];
Pml~{idom}~{1} += Region[ ((idom+1)*1000+300) ];
If(WALLS == 1)
PmlD0~{idom}~{1} += Region[ ((idom+1)*1000+302) ];
EndIf
If(WALLS == 0)
PmlInf~{idom}~{1} += Region[ ((idom+1)*1000+302) ];
EndIf
PmlInf~{idom}~{1} += Region[ ((idom+1)*1000+4) ];
EndIf
If (idom == N_DOM-1)
Sigma~{idom}~{0} = Region[ ((idom+1)*1000+10) ];
Sigma~{idom}~{1} = Region[{}];
GammaD0~{idom} += Region[ ((idom+1)*1000+20) ];
GammaD~{idom} = Region[{}];
Pml~{idom}~{0} += Region[ ((idom+1)*1000+100) ];
If(WALLS == 1)
PmlD0~{idom}~{0} += Region[ ((idom+1)*1000+102) ];
EndIf
If(WALLS == 0)
PmlInf~{idom}~{0} += Region[ ((idom+1)*1000+102) ];
EndIf
PmlInf~{idom}~{0} += Region[ ((idom+1)*1000+1) ];
EndIf
If (idom >= 1 && idom < N_DOM-1)
Sigma~{idom}~{0} = Region[ ((idom+1)*1000+10) ];
Sigma~{idom}~{1} = Region[ ((idom+1)*1000+20) ];
GammaD~{idom} = Region[{}];
Pml~{idom}~{0} += Region[ ((idom+1)*1000+100) ];
Pml~{idom}~{1} += Region[ ((idom+1)*1000+300) ];
If(WALLS == 1)
PmlD0~{idom}~{0} += Region[ ((idom+1)*1000+102) ];
PmlD0~{idom}~{1} += Region[ ((idom+1)*1000+302) ];
EndIf
If(WALLS == 0)
PmlInf~{idom}~{0} += Region[ ((idom+1)*1000+102) ];
PmlInf~{idom}~{1} += Region[ ((idom+1)*1000+302) ];
EndIf
PmlInf~{idom}~{0} += Region[ ((idom+1)*1000+1) ];
PmlInf~{idom}~{1} += Region[ ((idom+1)*1000+4) ];
EndIf
Sigma~{idom} = Region[{Sigma~{idom}~{0}, Sigma~{idom}~{1}}] ;
BndSigma~{idom}~{0} = Region[{}];
BndSigma~{idom}~{1} = Region[{}];
BndSigma~{idom} = Region[{BndSigma~{idom}~{0}, BndSigma~{idom}~{1}}] ;
BndGammaInf~{idom}~{0} = Region[{}];
BndGammaInf~{idom}~{1} = Region[{}];
BndGammaInf~{idom} = Region[{BndGammaInf~{idom}~{0}, BndGammaInf~{idom}~{1}}] ;
EndFor
}
Function{
// definitions for parallel (MPI) runs:
ListOfDom = {} ; // the domains that I'm in charge of
ListOfField = {}; // my fields
ListOfNeighborField = {}; // my neighbors
// this describes a layered (1-d like) decomposition
// +------+------+------+---...---+------+
// field: | 0|1 2|3 4|5 2N-4|2N-3 |
// idom: | 0 | 1 | 2 | | N-1 |
// +------+------+------+---...---+------+
For idom In {0:N_DOM-1}
If (idom % MPI_Size == MPI_Rank)
If(idom == 0)
// my fields
myFieldLeft = {};
myFieldRight = {0};
// fields to exchange with
exchangeFieldLeft = {};
exchangeFieldRight = {1};
// as many "blocks" as I have fields
ListOfNeighborField += 1;
ListOfNeighborField += exchangeFieldRight();
EndIf
If(idom == N_DOM-1)
myFieldLeft = {2*idom-1};
myFieldRight = {};
exchangeFieldLeft = {2*(idom-1)};
exchangeFieldRight = {};
ListOfNeighborField += 1;
ListOfNeighborField += exchangeFieldLeft();
EndIf
If(idom > 0 && idom < N_DOM-1)
myFieldLeft = {2*idom-1};
myFieldRight = {2*idom};
exchangeFieldLeft = {2*(idom-1)};
exchangeFieldRight = {2*idom+1};
// 2 "blocks"
ListOfNeighborField += 1;
ListOfNeighborField += exchangeFieldLeft();
ListOfNeighborField += 1;
ListOfNeighborField += exchangeFieldRight();
EndIf
ListOfDom += idom;
ListOfField += {myFieldLeft(), myFieldRight()};
If(ANALYSIS == 0)
g_in~{idom}~{0}[Sigma~{idom}~{0}] = ComplexScalarField[XYZ[]]{exchangeFieldLeft()};
g_in~{idom}~{1}[Sigma~{idom}~{1}] = ComplexScalarField[XYZ[]]{exchangeFieldRight()};
EndIf
If(ANALYSIS == 1)
g_in~{idom}~{0}[Sigma~{idom}~{0}] = ComplexVectorField[XYZ[]]{exchangeFieldLeft()};
g_in~{idom}~{1}[Sigma~{idom}~{1}] = ComplexVectorField[XYZ[]]{exchangeFieldRight()};
EndIf
EndIf
EndFor
/*
MPI_Printf["ListOfDom = ", ListOfDom()];
MPI_Printf["ListOfField = ", ListOfField()];
MPI_Printf["ListOfNeighborField = ", ListOfNeighborField()];
*/
}
Function{
// parameters for PML TC
xSigmaList = {};
thetaList = {};
For i In {0:N_DOM}
xSigmaList += i*dDom;
thetaList += theta;
EndFor
For ii In {0: N_DOM-1}
idom = ii;
xSigma~{idom}~{0} = xSigmaList(idom);
xSigma~{idom}~{1} = xSigmaList(idom+1);
theta~{idom}[] = thetaList(idom);
EndFor
For i In {0:N_DOM}
distSigma~{i}[] = X[]*Cos[theta] + Y[]*Sin[theta] - xSigmaList(i);
EndFor
For ii In {0: N_DOM-1}
idom = ii;
For jdom In {0:1}
kPml~{idom}~{jdom}[] = k[ Vector[xSigma~{idom}~{jdom},Y[],Z[]] ] ;
EndFor
// Bermudez damping functions
SigmaX[Pml~{idom}~{1}] = distSigma~{idom+1}[] > dTr ? c[]/(dBb-distSigma~{idom+1}[]) : 0;
SigmaX[Pml~{idom}~{0}] = -distSigma~{idom}[] > dTr ? c[]/Fabs[(dBb+distSigma~{idom}[])] : 0 ;
EndFor
SigmaY[] = 0.;
SigmaZ[] = 0.;
Kx[] = Complex[1, SigmaX[]/om[]];
Ky[] = Complex[1, SigmaY[]/om[]];
Kz[] = Complex[1, SigmaZ[]/om[]];
D[] = TensorDiag[Ky[]*Kz[]/Kx[], Kx[]*Kz[]/Ky[], Kx[]*Ky[]/Kz[]];
}
If(ANALYSIS == 0)
Include "Schrodinger.pro" ;
EndIf
If(ANALYSIS == 1)
Include "Maxwell.pro" ;
EndIf
DefineConstant[
// default getdp parameters for onelab
R_ = {"DDM", Name "GetDP/1ResolutionChoices", Visible 0},
C_ = {"-solve -v 3 -bin -ksp_monitor", Name "GetDP/9ComputeCommand", Visible 0},
P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}
];
DefineConstant[ // allows to set these from outside
// Analysis type
ANALYSIS = {0, Name "Input/00Type of analysis",
Choices {0="Helmholtz"}},
// wavenumber
WAVENUMBER = {30, Name "Input/0Wavenumber"},
LAMBDA = {2*Pi/WAVENUMBER, Name "Input/1Wavelength", ReadOnly 1},
// number of points per wavelength
N_LAMBDA = {20, Name "Input/2Points per wavelength"},
// dimensions of the waveguide
DX = {1, Name "Input/X dimension"},
DY = {1, Name "Input/Y dimension"},
DZ = {1, Name "Input/Z dimension"},
// number of subdmains in the DDM
N_DOM = {2, Name "Input/04Number of subdomains"},
// base msh filename
MSH_BASE_NAME = "mesh_",
// directory for output files
DIR = "out/"
];
LC = LAMBDA/N_LAMBDA;
// prefix for (split) mesh files (one for each partition)
MSH_NAME = StrCat(DIR, MSH_BASE_NAME) ;
nLayersTr = 1;
nLayersPml = 5;//*meshFactor;
dTr = nLayersTr*LC;
dPml = nLayersPml*LC;
dBb = (nLayersPml+nLayersTr)*LC;
theta = 0.;
dDom = DX / N_DOM;
Include "cavitySchrod2d_data.geo";
For idom In {0:N_DOM-1}
Merge Sprintf("out/u_%g.pos", idom);
EndFor
Combine ElementsFromVisibleViews;
View[0].Visible = 0;
For idom In {0:N_DOM-1}
Merge Sprintf("out/f_%g.pos", idom);
EndFor
Combine ElementsFromVisibleViews;
View[0].Visible = 1;
// Plugin(MathEval).Expression0 = "v0-w0";
// Plugin(MathEval).View = 0;
// Plugin(MathEval).OtherView = 1;
// Plugin(MathEval).ForceInterpolation = 0;
// Plugin(MathEval).Run;
// File "LaplacianNeumann.geo".
// We include the file containing the numbering of the geometry.
// This is usefull at the end of this file, and used to "synchronise" GMSH and GetDP
Include "LaplacianNeumann.param";
//Caracteristic length of the finite elements (reffinement is also possible after the mesh is built):
lc = 0.05;
// This parameter could be placed for instance in "param.geo", to separate more easyly the geometry
// and the discretization parameters.
// The parameters of the border of the domain :
x_max = 1;
x_min = 0;
y_max= 1;
y_min = 0;
//Creation of the 4 angle points of the domain Omega (=square)
p1 = newp; Point(p1) = {x_min,y_min,0,lc};
p2 = newp; Point(p2) = {x_min,y_max,0,lc};
p3 = newp; Point(p3) = {x_max,y_max,0,lc};
p4 = newp; Point(p4) = {x_max,y_min,0,lc};
// Remarks:
// -"newp" is a GMSH function that give the first available number for describing a point.
// For any other entity, like Line, Surface, etc. We recommand the use of "newreg" (see below).
// - By default, GMSH create a 3D domain. The z-coordinate must always be precised.
//The four edges of the square
L1 = newreg; Line(L1) = {p1,p2};
L2 = newreg; Line(L2) = {p2,p3};
L3 = newreg; Line(L3) = {p3,p4};
L4 = newreg; Line(L4) = {p4,p1};
// Line Loop (= boundary of the square)
Bound = newreg; Line Loop(Bound) = {L1,L2,L3,L4};
//Surface of the square
SurfaceOmega = newreg; Plane Surface(SurfaceOmega) = {Bound};
// To conclude, we define the physical entities, that is "what GetDP could see/use".
// "Omega" is a number imported from the file "param.geo".
Physical Surface(Omega) = {SurfaceOmega};
// Do not forget to let a blank line at the end, this could make GMSH crash...
// File "LaplacianNeumann.param"
//Numbers that caracterise the interior of the square (Omega) and its boundary (Gama):
Omega = 1000;
// Three remarks on these numbers :
// - They are arbitrary choosen.
// - They are placed in a separated file to be readable by both GMSH and GetDP.
// - "Gamma" is a special word used by GMSH/GetDP, that is why the boundary is named "Gama", with one "m"...
// Do not forget to let a blank line at the end, this could make GMSH crash...
// File "LaplacianNeumann.pro"
// As in the .geo file, we include the file containing the numbering of the geometry.
Include "LaplacianNeumann.param";
// Group
//======
//We now build the "Groups", that is, the geometrical entities, the different domains of computation.
//Here we only need the interior of the open domain Omega
Group{
Omega = Region[{Omega}];
}
// Function
//=========
Function{
// Pi is a special value in GetDp (=3.1415...)
Coef = Pi;
// Definition of the source function f(x,y)
f[] = (1+2*Coef*Coef)*Cos[Coef*X[]]*Cos[Coef*Y[]];
}
/*
Remark:
- the argument (as "x" and "y") are not writen between the bracket [].
Indeed, between the bracket is writen the domain of definition of the function.
- The argument "x" and "y" are here designed by the GetDP inner functions X[] and Y[], which give
respectively the x-coordinate and the y-coordinate of a considered point.
- To define a function globaly (i.e: not only on a subdomain), we write:
f[] = ...
- In our example, we could also have written
f[Omega] = (1+2*Coef*Coef)*Cos[Coef*X[]]*Cos[Coef*Y[]];
*/
//Jacobian
//========
Jacobian {
{ Name JVol ;
Case {
{ Region All ; Jacobian Vol ; }
}
}
{ Name JSur ;
Case {
{ Region All ; Jacobian Sur ; }
}
}
{ Name JLin ;
Case {
{ Region All ; Jacobian Lin ; }
}
}
}
/*
Remark: roughly speacking, we make use of...:
- Jacobian "Vol" as far as the integration domain is of same dimension than the problem (e.g: calculation in a 3D domain, in a 3D problem,
a 2D domain (surface) in a 2 problem, a 1D domain (line) in a 1D problem)
- Jacobian "Sur" when the domain of integration has one dimension less than the global problem (e.g:
calculation on a surface (2D) in a 3D problem, calculation on a line (1D) in a 2D problem).
- Jacobian "Lin" when the domain of integration has 2 dimensions less than the problem. That is,
for example, a calculation on a line (1D) in a 3D problem.
- Here, we just have define some Jacobian, you will see later that we just use the Jacobian "JVol"
*/
//Integration (parameters)
//=======================
Integration {
{ Name I1 ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 4 ; }
{ GeoElement Triangle ; NumberOfPoints 6 ; }
{ GeoElement Quadrangle ; NumberOfPoints 7 ; }
{ GeoElement Tetrahedron ; NumberOfPoints 15 ; }
{ GeoElement Hexahedron ; NumberOfPoints 34 ; }
}
}
}
}
}
// There is no Constrain because of the Neumann boundary condition
// We go directly to the FuntionSpace
//FunctionSpace
//=============
FunctionSpace{
{ Name Vh; Type Form0;
BasisFunction{
{Name wn; NameOfCoef vn; Function BF_Node;
Support Omega; Entity NodesOf[All];}
}
}
}
/*
Explanation:
The space of approximation is called "Vh".
Its type is "Form0", that means "scalar".
The basis functions are called "wn", and the associated coefficients "vn".
In other words, a function "v" of "Vh" can be written as
v(x,y) = sum_{n} vn.wn(x,y)
The functions "wn" are nodale functions P1 ("BF_Node"), supported on the domain Omega ("Support Omega").
*/
//(Weak) Formulation
//==================
Formulation{
{Name LaplacianNeumann; Type FemEquation;
// We decided to call the formulation "LaplacianNeumann".
// Its type is "FemEquation" ("Finite Element Method")
Quantity{
{Name u; Type Local; NameOfSpace Vh;}
}
// Here, we introduce a quantity "u", which belongs to the function space Vh, defined above.
Equation{
Galerkin{ [Dof{Grad u}, {Grad u}];
In Omega; Jacobian JVol; Integration I1;}
Galerkin{ [ Dof{u}, {u}];
In Omega; Jacobian JVol; Integration I1;}
Galerkin{ [-f[], {u}];
In Omega; Jacobian JVol; Integration I1;}
}
}
}
/*
The variationnal formulation is written between the acolades {}, after "Equation".
Let us first give some vocabulary:
- "Galerkin" : GetDP syntaxic word. This could be translated mathematicaly by "integration" (see below)
- Dof{u}: "Degree Of Freedom". This is used to specify that the quantity is the unknown.
If "Dof" is not written, then "u" is seen as a test function and not as the unknown.
(Be carreful, the unknown and the test functions has the same name in GetDP !
The "Dof" is here to distinguish the unknown to the test-functions)
Now, let us detail the program written between the acolades {}, :
Galerkin{ [Dof{Grad u}, {Grad u}];
In Omega; Jacobian JVol; Integration Int;}
This could be translated mathematicaly by :
Integration on "Omega" of Grad(u).Grad(v)
where "u" is the unknown and "v" a test function.
Note the use of the Jacobian JVol(2D problem and integration on a 2D domain).
Moreover the number of integrations point is given in "I1" (see above "Integration").
The total variationnal formulation can be read as :
search "Dof{u}" in V_h such that, for every "{u}" in V_h,
\int_{\Omega} Grad( Dof{u}) . Grad( {u}) d\Omega + \int_{\Omega} Dof{u}.{u} d\Omega - ...
... - \int_{\Omega} f.{u} d\Omega = 0
(this is exactly the weak formulation of our problem !)
Remarks:
- Between two "Galerkin", a positive sign "+" is implicitly written
- The sum of all integrals "Galerkin" is equal to 0 (do not forget the "minus" sign of the right hand side)
- Why do we use a volumic jacobian even in a 2D problem? The problem is a 2 dimensional problem and
the integral is defined on a 2D domain (Omega). If the integral was writen on, e.g, the boundary,
then the Jacobian "Jsur" would have been used.
*/
// Resolution (of the problem)
//============================
Resolution{
{Name LaplacianNeumann;
// We chose the name LaplacianNeumann for the resolution
// Remark: in GetDP, every entity has a name. The same same can be used for different entities but of different kind, of course
// Here we chose the same name as the formulation, but this is just a choice, no obligation !
System{
{Name Syst; NameOfFormulation LaplacianNeumann;}
}
// A system is linked to a weak formulation
// Here, we only have one weak formulation, which is "LaplacianNeumann"
Operation{
Generate[Syst]; Solve[Syst]; SaveSolution[Syst];
}
/* When we launch GetDP, the program will ask the user to chose a Resolution.
When calling the Resolution "LaplacianNeumann", GetDP will...
- Generate the system
- Solve the system
- Save the solution
Note that GetDP respects the order of the operation !
*/
}
}
//Post Processing
//===============
PostProcessing{
{Name LaplacianNeumann; NameOfFormulation LaplacianNeumann;
Quantity{
{Name u; Value {Local{[{u}];In Omega;Jacobian JVol;}}}
}
}
}
/*
The name of the PostProcessing is LaplacianNeumann.
It call the weak formulation LaplacianNeumann, then take the solution "u" and compute it on all the domain Omega, by the operation:
Value {Local{[{u}];In Carre; Jacobian JVol;}}}
This means that "u" is the interpolate of the solution "Dof{u}" computed in the weak formulation "LaplacianNeumann"
Remark: Again, we chose the same name, but this is just a choice, this has no influence on GetDP resolution.
*/
//Post Operation
//==============
PostOperation{
{Name Map_u; NameOfPostProcessing LaplacianNeumann;
Operation{
Print[u, OnElementsOf Omega, File "u_Neumann.pos"];
}
}
}
/*
The only PostOperation we write is to display "u" introduced in the PostProcessing.
This PostOperation is called "Map_u".
When launching GetDP, it will ask the user to chose :