Commit ac3356da authored by Christophe Geuzaine's avatar Christophe Geuzaine

benchmarks have moved to https://gitlab.onelab.info/doc/models

parent 3bf42ab2
Pipeline #2119 passed with stage
in 9 minutes 56 seconds
Model developed by [[User:Axel Modave|A. Modave]]. Onelab wrapping by [[User:Francois Henrotte|F. Henrotte]].
Academic example for eigenvalue problems.
[[Category:GetDP]]
[[Category:Eigenvalue problem]]
<!--
To run the example, open main.pro with Gmsh.
-->
//========================================================
// GetDP code for simulation of EIGENVALUES PROBLEMS
// - Helmholtz equation (scalar/vector and 1D/2D/3D)
// - Nodal/Edge finite-elements (Form0/Form1 for u)
// Contributors: A. Itagi (original version, 2007),
// B. Kubicek (minor modif), A. Modave (major modif)
//========================================================
Group {
DefineGroup[ Dom, Bnd, Tot ] ;
}
Jacobian {
{ Name Jac ;
Case {
{ Region Bnd ; Jacobian Sur ; }
{ Region Dom ; Jacobian Vol ; }
}
}
}
Integration {
{ Name I1 ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 4 ; }
{ GeoElement Triangle ; NumberOfPoints 4 ; }
{ GeoElement Quadrangle ; NumberOfPoints 4 ; }
{ GeoElement Tetrahedron ; NumberOfPoints 4 ; }
{ GeoElement Hexahedron ; NumberOfPoints 6 ; }
{ GeoElement Prism ; NumberOfPoints 9 ; }
}
}
}
}
}
Constraint {
{ Name uConstraint ; Type Assign ;
Case {
{ Region Bnd ; Value 0. ; }
}
}
}
FunctionSpace {
If (FLAG_FORM==0)
{ Name uSpace ; Type Form0 ;
BasisFunction {
{ Name sn ; NameOfCoef en ; Function BF_Node ; Support Tot ; Entity NodesOf[All] ; }
}
Constraint {
{ NameOfCoef en ; EntityType NodesOf ; NameOfConstraint uConstraint ; }
}
}
EndIf
If (FLAG_FORM==1)
{ Name uSpace ; Type Form1 ;
BasisFunction {
{ Name sn ; NameOfCoef en ; Function BF_Edge ; Support Tot ; Entity EdgesOf[All] ; }
}
Constraint {
{ NameOfCoef en ; EntityType EdgesOf ; NameOfConstraint uConstraint ; }
}
}
EndIf
If (FLAG_FORM==10)
{ Name uSpace ; Type Form1P ;
BasisFunction {
{ Name sn ; NameOfCoef en ; Function BF_PerpendicularEdge ; Support Tot ; Entity NodesOf[All] ; }
}
Constraint {
{ NameOfCoef en ; EntityType NodesOf ; NameOfConstraint uConstraint ; }
}
}
EndIf
}
Formulation {
{ Name Form ; Type FemEquation ;
Quantity {
{ Name u ; Type Local ; NameOfSpace uSpace ; }
}
Equation {
Galerkin { DtDtDof[ Dof{u} , {u} ] ;
In Dom ; Integration I1 ; Jacobian Jac ; }
Galerkin { [ Dof{d u} , {d u} ] ;
In Dom ; Integration I1 ; Jacobian Jac ; }
}
}
}
Resolution {
{ Name Reso ;
System {
{ Name A ; NameOfFormulation Form ;
Type ComplexValue ;
}
}
Operation {
CreateDir["output/"] ;
GenerateSeparate[A] ;
EigenSolve[A,NbEigenvalues,EigenvalShiftRe,EigenvalShiftIm] ;
SaveSolutions[A] ;
}
}
}
PostProcessing {
{ Name PostPro ; NameOfFormulation Form ;
Quantity {
{ Name u ; Value{ Local{ [{u}] ; In Dom ; Jacobian Jac ; } } }
}
}
}
PostOperation {
{ Name PostOp ; NameOfPostProcessing PostPro ;
Operation {
For n In {0:(NbEigenvalues-1)}
Print [ u, OnElementsOf Dom, TimeStep{n}, File StrCat["output/eigenvector",Sprintf("%g",n),".pos"]];
EndFor
SendMergeFileRequest[ "main.pro.opt" ] ;
}
}
}
DefineConstant[
R_ = {"Reso", Name "GetDP/1ResolutionChoices", Visible 0},
P_ = {"PostOp", Name "GetDP/2PostOperationChoices", Visible 0, ReadOnly 1},
C_ = {"-solve -pos -v2", Name "GetDP/9ComputeCommand", Visible 0}
] ;
//========================================================
// Benchmark "Academic eigenvalues problems"
// File: GMSH geometry (cuboid)
//========================================================
DefineConstant[
FLAG_MESH = {2, Highlight "Black",
Name StrCat[OnelabParamMesh,"1Shape of cells"],
Choices {1="Tetrahedral cells", 2="Hexahedral cells"} },
res = { 0.1, Min 0.01, Max 1, Step 0.01, Visible (FLAG_MESH==1),
Name StrCat[OnelabParamMesh,"2Characteristic length of cells (length of borders = 1)"]},
resHexa = { 30, Min 1, Max 100, Step 1, Visible (FLAG_MESH==2),
Name StrCat[OnelabParamMesh,"3Number of cells along each border"]}
];
p1[] += newp ; Point(newp) = {-0.5,-0.5,-0.5} ;
p1[] += newp ; Point(newp) = { 0.5,-0.5,-0.5} ;
p1[] += newp ; Point(newp) = { 0.5, 0.5,-0.5} ;
p1[] += newp ; Point(newp) = {-0.5, 0.5,-0.5} ;
p2[] += newp ; Point(newp) = {-0.5,-0.5, 0.5} ;
p2[] += newp ; Point(newp) = { 0.5,-0.5, 0.5} ;
p2[] += newp ; Point(newp) = { 0.5, 0.5, 0.5} ;
p2[] += newp ; Point(newp) = {-0.5, 0.5, 0.5} ;
l1[] += newl ; Line(newl) = {p1[0], p1[1]} ;
l1[] += newl ; Line(newl) = {p1[1], p1[2]} ;
l1[] += newl ; Line(newl) = {p1[2], p1[3]} ;
l1[] += newl ; Line(newl) = {p1[3], p1[0]} ;
l2[] += newl ; Line(newl) = {p2[0], p2[1]} ;
l2[] += newl ; Line(newl) = {p2[1], p2[2]} ;
l2[] += newl ; Line(newl) = {p2[2], p2[3]} ;
l2[] += newl ; Line(newl) = {p2[3], p2[0]} ;
lL[] += newl ; Line(newl) = {p1[0], p2[0]} ;
lL[] += newl ; Line(newl) = {p1[1], p2[1]} ;
lL[] += newl ; Line(newl) = {p1[2], p2[2]} ;
lL[] += newl ; Line(newl) = {p1[3], p2[3]} ;
ll1 = newll ; Line Loop(newll) = {-l1[]} ;
s1 = news ; Plane Surface(news) = {ll1} ;
ll2 = newll ; Line Loop(newll) = {l2[]} ;
s2 = news ; Plane Surface(news) = {ll2} ;
llL[] += newll ; Line Loop(newll) = {l1[0], lL[1], -l2[0], -lL[0]} ;
llL[] += newll ; Line Loop(newll) = {l1[1], lL[2], -l2[1], -lL[1]} ;
llL[] += newll ; Line Loop(newll) = {l1[2], lL[3], -l2[2], -lL[2]} ;
llL[] += newll ; Line Loop(newll) = {l1[3], lL[0], -l2[3], -lL[3]} ;
sL[] += news ; Plane Surface(news) = {llL[0]} ;
sL[] += news ; Plane Surface(news) = {llL[1]} ;
sL[] += news ; Plane Surface(news) = {llL[2]} ;
sL[] += news ; Plane Surface(news) = {llL[3]} ;
sl = newsl ; Surface Loop(newsl) = {s1, s2, sL[]} ;
v = newv ; Volume(newv) = {sl} ;
Physical Surface(BND) = {s1, s2, sL[]} ;
Physical Volume(DOM) = {v} ;
If (FLAG_MESH==1)
Mesh.CharacteristicLengthMax = res ;
EndIf
If (FLAG_MESH==2)
Transfinite Line {l1[0], l1[2], l2[0], l2[2]} = resHexa+1 ;
Transfinite Line {l1[1], l1[3], l2[1], l2[3]} = resHexa+1 ;
Transfinite Line {lL[]} = resHexa+1 ;
Transfinite Surface {s1, s2, sL[]} ;
Recombine Surface {s1, s2, sL[]} ;
Transfinite Volume {v} ;
Recombine Volume {v} ;
EndIf
//========================================================
// Benchmark "Academic eigenvalues problems"
// File: GMSH geometry (disk)
//========================================================
DefineConstant[
res = { 0.01, Min 0.001, Max 1, Step 0.001,
Name StrCat[OnelabParamMesh,"2Characteristic length of cells (radius of domain = 1)"]}
];
Mesh.CharacteristicLengthMax = res ;
p[] += newp ; Point(newp) = { 0 , 0 , 0} ;
p[] += newp ; Point(newp) = { 0.5, 0 , 0} ;
p[] += newp ; Point(newp) = { 0 , 0.5, 0} ;
p[] += newp ; Point(newp) = {-0.5, 0 , 0} ;
p[] += newp ; Point(newp) = { 0 ,-0.5, 0} ;
l[] += newl ; Circle(newl) = {p[1], p[0], p[2]} ;
l[] += newl ; Circle(newl) = {p[2], p[0], p[3]} ;
l[] += newl ; Circle(newl) = {p[3], p[0], p[4]} ;
l[] += newl ; Circle(newl) = {p[4], p[0], p[1]} ;
ll = newll ; Line Loop(newll) = {l[]} ;
s = news ; Plane Surface(news) = {ll} ;
Physical Line(BND) = {l[]} ;
Physical Surface(DOM) = {s} ;
//========================================================
// Benchmark "Academic eigenvalues problems"
// File: GMSH geometry (line)
//========================================================
DefineConstant[
res = { 200, Min 1, Max 10000, Step 1,
Name StrCat[OnelabParamMesh,"1Number of cells"]}
];
p[] += newp ; Point(newp) = {-0.5, 0, 0} ;
p[] += newp ; Point(newp) = { 0.5, 0, 0} ;
l = newl ; Line(newl) = {p[0],p[1]} ;
Physical Point(BND) = {p[]} ;
Physical Line(DOM) = {l} ;
Transfinite Line {l} = res+1 ;
//========================================================
// Benchmark "Academic eigenvalues problems"
// File: GMSH geometry (square)
//========================================================
DefineConstant[
FLAG_MESH = {2, Highlight "Black",
Name StrCat[OnelabParamMesh,"1Shape of cells"],
Choices {1="Triangular cells", 2="Rectangular cells"} },
res = { 0.05, Min 0.001, Max 1, Step 0.001, Visible (FLAG_MESH==1),
Name StrCat[OnelabParamMesh,"2Characteristic length of cells (length of borders = 1)"]},
resRect = { 200, Min 1, Max 1000, Step 1, Visible (FLAG_MESH==2),
Name StrCat[OnelabParamMesh,"3Number of cells along each border"]}
];
p[] += newp ; Point(newp) = {-0.5, -0.5, 0} ;
p[] += newp ; Point(newp) = { 0.5, -0.5, 0} ;
p[] += newp ; Point(newp) = { 0.5, 0.5, 0} ;
p[] += newp ; Point(newp) = {-0.5, 0.5, 0} ;
l[] += newl ; Line(newl) = {p[0], p[1]} ;
l[] += newl ; Line(newl) = {p[1], p[2]} ;
l[] += newl ; Line(newl) = {p[2], p[3]} ;
l[] += newl ; Line(newl) = {p[3], p[0]} ;
ll = newll ; Line Loop(newll) = {l[]} ;
s = news ; Plane Surface(news) = {ll} ;
Physical Line(BND) = {l[]} ;
Physical Surface(DOM) = {s} ;
If (FLAG_MESH==1)
Mesh.CharacteristicLengthMax = res ;
EndIf
If (FLAG_MESH==2)
Transfinite Line {l[0], l[2]} = resRect+1 ;
Transfinite Line {l[1], l[3]} = resRect+1 ;
Transfinite Surface {s} ;
Recombine Surface {s} ;
EndIf
//========================================================
// Benchmark "Academic eigenvalues problems"
// File: Parameters
//========================================================
ExtGnuplot = ".dat";
ExtGmsh = ".pos";
DirOutput = "output/";
OnelabParamMesh = "Input/1Mesh parameter(s)/";
OnelabParamSolver = "Input/2Resolution parameters/";
// GmshOption "Reset",
FLAG_GEOM = DefineNumber[1, Name "Input/01Geometry", Highlight "Black",
GmshOption "Reset",
Choices {1="Linear domain [1D]",
2="Squared domain [2D]",
3="Circular domain [2D]",
4="Cuboidal domain [3D]"}];
NbEigenvalues = DefineNumber[ 4, Min 1, Max 100, Step 1,
Name StrCat[OnelabParamSolver,"3Number of eigenvalues to compute"] ] ;
EigenvalShiftRe = DefineNumber[ 0, Min 0, Max 1000, Step 0.001,
Name StrCat[OnelabParamSolver,"4Spectral shift (real part)"] ] ;
EigenvalShiftIm = DefineNumber[ 0, Min 0, Max 1000, Step 0.001,
Name StrCat[OnelabParamSolver,"5Spectral shift (imaginary part)"] ] ;
If (FLAG_GEOM==1)
LinkGeo = "geomLine.geo" ;
DIM = 1 ;
EndIf
If (FLAG_GEOM==2)
LinkGeo = "geomSquare.geo" ;
DIM = 2 ;
EndIf
If (FLAG_GEOM==3)
LinkGeo = "geomDisk.geo" ;
DIM = 2 ;
EndIf
If (FLAG_GEOM==4)
LinkGeo = "geomCuboid.geo" ;
DIM = 3 ;
EndIf
BND = 1000;
DOM = 2000;
//========================================================
// Benchmark "Academic eigenvalues problems"
// File: GMSH geometry (choice of the cavity)
//========================================================
Include "main.dat" ;
Include Str[LinkGeo] ;
Solver.AutoShowLastStep = 0;
//Solver.AutoCheck = 0;
Solver.AutoShowViews = 2;
Solver.AutoMesh = 2;
General.ScaleX = 1;
General.ScaleY = 1;
General.ScaleZ = 1;
General.TrackballQuaternion0 = 0;
General.TrackballQuaternion1 = 0;
General.TrackballQuaternion2 = 0;
General.TrackballQuaternion3 = 1;
FLAG_GEOM = DefineNumber[0, Name "Input/01Geometry"];
Mesh.SurfaceEdges=0; // do not display the mesh
k = PostProcessing.NbViews - 1;
If (FLAG_GEOM == 1)
General.TrackballQuaternion0 = 0;
General.TrackballQuaternion1 = 0;
General.TrackballQuaternion2 = 0;
General.TrackballQuaternion3 = 1;
General.ScaleX = 0.5;
General.ScaleY = 0.5;
General.ScaleZ = 0.5;
For k In {0 : PostProcessing.NbViews - 1}
View[k].Name = Sprintf("%g th mode", k);
View[k].Axes = 2;
View[k].Type = 2;
View[k].RaiseZ = 0;
View[k].Visible=1;
EndFor
EndIf
If (FLAG_GEOM == 2)
General.ScaleX = 0.5;
General.ScaleY = 0.5;
General.ScaleZ = 0.5;
General.TrackballQuaternion0 = 0.5237405039445426;
General.TrackballQuaternion1 = 0.08505041242690579;
General.TrackballQuaternion2 = 0.2061022875509479;
General.TrackballQuaternion3 = 0.8221825581585859;
For k In {0 : PostProcessing.NbViews - 1}
View[k].Name = Sprintf("%g th mode", k);
View[k].Axes = 0;
View[k].Type = 0;
View[k].RaiseZ = 1;
View[k].Visible=0;
EndFor
View[PostProcessing.NbViews - 1].Visible=1;
EndIf
//========================================================
// Benchmark "Academic eigenvalues problems"
// File: GetDP simulation (choice of the problem)
//========================================================
Include "main.dat" ;
FLAG_EQN = DefineNumber[1, Name "Input/02Problem", Highlight "Black",
Choices {1="Scalar Helmholtz equation", 2="Vector Helmholtz equation"}];
If (FLAG_EQN==1)
FLAG_FORM = 0 ; // Nodal elements - Form0
EndIf
If ((FLAG_EQN==2) && ((DIM==1) || (DIM==2)))
FLAG_FORM = 10 ; // Edge elements - Form1P (with rotation for 1D/2D)
EndIf
If ((FLAG_EQN==2) && (DIM==3))
FLAG_FORM = 1 ; // Edge elements - Form1 (only for 3D)
EndIf
Group {
Bnd = Region[{BND}] ;
Dom = Region[{DOM}] ;
Tot = Region[{Bnd,Dom}] ;
}
Include "formulation.pro" ;
FLAG_GEOM = DefineNumber[0, Name "Input/01Geometry"];
Mesh.SurfaceEdges=0; // do not display the mesh
k = PostProcessing.NbViews - 1;
If (FLAG_GEOM == 1)
For k In {0 : PostProcessing.NbViews - 1}
View[k].Name = Sprintf("%g th mode", k);
View[k].Axes = 2;
View[k].Type = 2;
View[k].RaiseZ = 0;
View[k].Visible=1;
EndFor
EndIf
If (FLAG_GEOM == 2 || FLAG_GEOM == 3)
General.ScaleX = 0.5;
General.ScaleY = 0.5;
General.ScaleZ = 0.5;
General.TrackballQuaternion0 = 0.5237405039445426;
General.TrackballQuaternion1 = 0.08505041242690579;