From c9d1f46d998f5a3df5ac88b762eba0be3b5296af Mon Sep 17 00:00:00 2001 From: Christophe Geuzaine <cgeuzaine@ulg.ac.be> Date: Thu, 15 Mar 2018 18:35:36 +0100 Subject: [PATCH] Update microstrip.pro --- Electrostatics/microstrip.pro | 648 +++++++++++++++++----------------- 1 file changed, 324 insertions(+), 324 deletions(-) diff --git a/Electrostatics/microstrip.pro b/Electrostatics/microstrip.pro index f194945..e3524fe 100644 --- a/Electrostatics/microstrip.pro +++ b/Electrostatics/microstrip.pro @@ -1,324 +1,324 @@ -/* ------------------------------------------------------------------- - Tutorial 1 : electrostatic field of a microstrip - - Features: - - Physical and abstract regions - - Scalar FunctionSpace with Dirichlet constraint - - Galerkin term for stiffness matrix - - To compute the solution in a terminal: - getdp microstrip -solve EleSta_v - getdp microstrip -pos Map - getdp microstrip -pos Cut - - To compute the solution interactively from the Gmsh GUI: - File > Open > microstrip.pro - Run (button at the bottom of the left panel) - ------------------------------------------------------------------- */ - -/* In this first tutorial we consider the calculation of the electric field - given a static distribution of electric potential. This corresponds to an - "electrostatic" physical model, obtained by combining the time-invariant - Maxwell-Ampere equation (Curl e = 0, with e the electric field) with Gauss' - law (Div d = rho, with d the displacement field and rho the charge density) - and the dielectric constitutive law (d = epsilon e, with epsilon the - dielectric permittivity). - - Since Curl e = 0, e can be derived from a scalar electric potential v, such - that e = -Grad v. Plugging this potential in Gauss' law and using the - constitutive law leads to a scalar Poisson equation in terms of the electric - potential: -Div(epsilon Grad v) = rho. - - We consider here the special case where rho = 0, to model a conducting - microstrip line on top of a dielectric substrate. A Dirichlet boundary - condition sets the potential to 1 mV on the boundary of the microstrip line - (called "Electrode" below) and 0 V on the ground. A homogeneous Neumann - boundary condition (zero flux of the displacement field, i.e. n.d = 0) is - imposed on the left boundary of the domain to account for the symmetry of the - problem, as well as on the top and right boundaries that truncate the - modelling domain. */ - -Group { - /* One starts by giving explicit meaningful names to the Physical regions - defined in the "microstrip.msh" mesh file. We only use 2 volume regions and - 2 surface regions in this model. */ - - Air = Region[101]; - Diel1 = Region[111]; - - Ground = Region[120]; - Electrode = Region[121]; - - /* We now define abstract regions to be used below in the definition of the - scalar electric potential formulation: - - Vol_Ele : volume where -Div(epsilon Grad v) = 0 is solved - Sur_Neumann_Ele: surface where non homogeneous Neumann boundary conditions - (on n.d == epsilon n.Grad v) are imposed - - Vol_xxx groups contain only volume elements of the mesh (triangles here). - Sur_xxx groups contain only surface elements of the mesh (lines here). - - Since there are no non-homogeneous Neumann conditions in the model, - Sur_Neumann_Ele is defined as empty. - */ - - Vol_Ele = Region[ {Air, Diel1} ]; - Sur_Neumann_Ele = Region[ {} ]; -} - -Function { - /* Material laws (here the dielectric permittivity) are defined as piecewise - functions (note the square brackets), in terms of the above defined - physical regions. */ - - eps0 = 8.854187818e-12; - epsilon[Air] = 1. * eps0; - epsilon[Diel1] = 9.8 * eps0; -} - -Constraint { - /* As for material laws, the Dirichlet boundary condition is defined - piecewise. The constraint "Dirichlet_Ele" is invoked in the FunctionSpace - below. */ - - { Name Dirichlet_Ele; Type Assign; - Case { - { Region Ground; Value 0.; } - { Region Electrode; Value 1.e-3; } - } - } -} - -Group{ - /* The domain of definition of a FunctionSpace lists all regions on which a - field is defined. - - Contrary to Vol_xxx and Sur_xxx regions, which contain only volume or - surface regions, resp., domains of definition Dom_xxx regions may contain - both volume and surface regions. Hence the use of the prefixes Vol_, Sur_ - and Dom_ to avoid confusions.*/ - - Dom_Hgrad_v_Ele = Region[ {Vol_Ele, Sur_Neumann_Ele} ]; -} - -FunctionSpace { - /* The function space in which we shall pick the electric scalar potential "v" - solution is definied by - - a domain of definition (the "Support": "Dom_Hgrad_v_Ele") - - a type ("Form0" means scalar field) - - a set of scalar basis functions ("BF_Node" means nodal basis functions) - - a set of entities to which the basis functions are associated ("Entity": - here all the nodes of the domain of definition "NodesOf[All]") - - a constraint (here the Dirichlet boundary conditions) - - The FE expansion of the unknown field "v" reads - - v(x,y) = Sum_k vn_k sn_k(x,y) - - where the "vn_k" are the nodal values (connectors) and "sn_k(x,y)" the - nodal basis functions. Not all connectors are unknowns of the FE problem, - due to the "Constraint", which assigns particular values to the nodes of - the Ground and Electrode regions. GetDP deals with that automatically on - basis of the definition of the FunctionSpace. */ - - { Name Hgrad_v_Ele; Type Form0; - BasisFunction { - { Name sn; NameOfCoef vn; Function BF_Node; - Support Dom_Hgrad_v_Ele; Entity NodesOf[ All ]; } - // using "NodesOf[All]" instead of "NodesOf[Dom_Hgrad_v_Ele]" is an - // optimization, which allows GetDP to not explicitly build the list of - // all nodes - } - Constraint { - { NameOfCoef vn; EntityType NodesOf; NameOfConstraint Dirichlet_Ele; } - } - } -} - -Jacobian { - /* Jacobians are used to specify the mapping between elements in the mesh and - the reference elements (defined in standardized unit cells) over which - integration is performed. "Vol" represents the classical 1-to-1 mapping - between elements of identical geometrical dimension, i.e. in this case a - reference triangle/quadrangle onto triangles/quadrangles in the z=0 plane - (2D <-> 2D). "Sur" would be used to map the reference triangle/quadrangle - onto triangles/quadrangles in a 3D space (2D <-> 3D), or to map the - reference line segment onto segments in 2D space (1D <-> 2D). "Lin" would - be used to map the reference line segment onto segments in 3D space (1D <-> - 3D). */ - - { Name Vol ; - Case { - { Region All ; Jacobian Vol ; } - } - } -} - -Integration { - /* A Gauss quadrature rule with 4 points is used for all integrations below. */ - - { Name Int ; - Case { { Type Gauss ; - Case { { GeoElement Triangle ; NumberOfPoints 4 ; } - { GeoElement Quadrangle ; NumberOfPoints 4 ; } } - } - } - } -} - -Formulation { - /* The syntax of the Formulation section is a harder nut to crack. So let's - deal with it carefully. - - A GetDP Formulation encodes a so-called weak formulation of the original - partial differential equation, i.e. of -Div(epsilon Grad v) = 0. This weak - formulation involves finding v such that - - (-Div(epsilon Grad v) , v')_Vol_Ele = 0 - - holds for all so-called "test-functions" v', where (.,.)_D denotes an inner - product over the domain D. If the test-functions v' are differentiable, - integration by parts using Green's identity leads to finding v such that - - (epsilon Grad v, Grad v')_Vol_Ele + (epsilon n.Grad v, v')_Bnd_Vol_Ele = 0 - - holds for all v', where Bnd_Vol_Ele is the boundary of Vol_Ele. In our - microstrip example this surface term vanishes, as either there is no test - function v' (on the Dirichlet boundary) or the "epsilon n.Grad v" is zero - (on the homogeneous Neumann boundary). We are thus eventually looking for - functions v in the function space Hgrad_v_Ele such that - - (epsilon Grad v, Grad v')_Vol_Ele = 0 - - holds for all v'. Finally, our choice here is to use a Gakerkin method, - where the test functions v' are the same basis functions ("sn_k") as the - ones used to interpolate the unknown potential v. - - The "Integral" statement in the Formulation is a symbolic representation of - this weak formulation. It has got 4 semicolon separated arguments: - * the density [.,.] to be integrated (note the square brackets instead of - the parentheses), with the test-functions (always) after the comma - * the domain of integration, - * the Jacobian of the transformation reference element <-> real element, - * the integration method to be used. - - In the density, braces around a quantity (such as "{v}") indicate that this - quantity belongs to a FunctionSpace. Differential operators can be applied - within braces (such as "{Grad v}"); in particular the symbol "d" represents - the exterior derivative, and it is a synonym of "Grad" when applied to a - scalar function (such as "{d v}"). - - What can be a bit confusing is that the two comma-separated terms of the - density are not interpreted exactly in the same way. Let us unravel this in - detail. - - As the Galerkin method uses as test functions the same basis functions - "sn_k" as for the unknown field "v", the second term in the density should - be something like this [ ... , basis_functions_of {d v} ]. However, since - the second term is always devoted to test functions, the operator - "basis_functions_of" would always be there. It can therefore be made - implicit, and, according to the GetDP syntax, it is omitted. So, one writes - simply [ ... , {d v} ]. - - The first term, on the other hand, can contain a much wider variety of - expressions than the second one. In our case it should be expressed in - terms of the FE expansion of "v" at the present system solution, i.e. when - the coefficients vn_k in the expansion of "v = Sum_k vn_k sn_k" are - unknown. This is indicated by prefixing the braces with "Dof", which leads - to the following density: - - [ epsilon[] * Dof{d v} , {d v} ], - - a so-called bilinear term that will contribute to the stiffness matrix of - the electrostatic problem at hand. - - Another option, which would not work here, is to evaluate the first - argument with the last available already computed solution, i.e. simply - perform the interpolation with known coefficients vn_k. For this the Dof - prefix operator would be omitted and we would have: - - [ epsilon[] * {d v} , {d v} ], - - a so-called linear term that will contribute to the right-hand side of the - linear system. - - Both choices are commonly used in different contexts, and we shall come - back on this often in subsequent tutorials. */ - { Name Electrostatics_v; Type FemEquation; - Quantity { - { Name v; Type Local; NameOfSpace Hgrad_v_Ele; } - } - Equation { - Integral { [ epsilon[] * Dof{d v} , {d v} ]; - In Vol_Ele; Jacobian Vol; Integration Int; } - } - } -} - -Resolution { - { Name EleSta_v; - System { - { Name Sys_Ele; NameOfFormulation Electrostatics_v; } - } - Operation { - Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele]; - } - } -} - -/* Post-processing is done in two parts. - - The first part defines, in terms of the Formulation, which itself refers to - the FunctionSpace, a number of quantities that can be evaluated at the - postprocessing stage. The three quantities defined here are: - - the scalar vector potential, - - the electric field, - - the electric displacement. - - The second part consists in defining groups of post-processing operations, - which can be invoked separately. The first group is invoked by default when - Gmsh is run interactively. Each Operation specifies - - a quantity to be eveluated, - - the region on which the evaluation is done, - - the name of the output file. - The generated post-processing files are automatically displayed by Gmsh if - the "Merge result automatically" option is enabled (which is the default). */ - -PostProcessing { - { Name EleSta_v; NameOfFormulation Electrostatics_v; - Quantity { - { Name v; Value { - Local { [ {v} ]; In Dom_Hgrad_v_Ele; Jacobian Vol; } - } - } - { Name e; Value { - Local { [ -{d v} ]; In Dom_Hgrad_v_Ele; Jacobian Vol; } - } - } - { Name d; Value { - Local { [ -epsilon[] * {d v} ]; In Dom_Hgrad_v_Ele; Jacobian Vol; } - } - } - } - } -} - -e = 1.e-7; // tolerance to ensure that the cut is inside the simulation domain -h = 2.e-3; // vertical position of the cut - -PostOperation { - { Name Map; NameOfPostProcessing EleSta_v; - Operation { - Print [ v, OnElementsOf Dom_Hgrad_v_Ele, File "mStrip_v.pos" ]; - Print [ e, OnElementsOf Dom_Hgrad_v_Ele, File "mStrip_e.pos" ]; - Print [ e, OnLine {{e,h,0}{14.e-3,h,0}}{60}, File "Cut_e.pos" ]; - } - } - { Name Cut; NameOfPostProcessing EleSta_v; - // same cut as above, with more points and exported in raw text format - Operation { - Print [ e, OnLine {{e,e,0}{14.e-3,e,0}} {500}, Format TimeTable, File "Cut_e.txt" ]; - } - } -} +/* ------------------------------------------------------------------- + Tutorial 1 : electrostatic field of a microstrip + + Features: + - Physical and abstract regions + - Scalar FunctionSpace with Dirichlet constraint + - Galerkin term for stiffness matrix + + To compute the solution in a terminal: + getdp microstrip -solve EleSta_v + getdp microstrip -pos Map + getdp microstrip -pos Cut + + To compute the solution interactively from the Gmsh GUI: + File > Open > microstrip.pro + Run (button at the bottom of the left panel) + ------------------------------------------------------------------- */ + +/* In this first tutorial we consider the calculation of the electric field + given a static distribution of electric potential. This corresponds to an + "electrostatic" physical model, obtained by combining the time-invariant + Faraday equation (Curl e = 0, with e the electric field) with Gauss' + law (Div d = rho, with d the displacement field and rho the charge density) + and the dielectric constitutive law (d = epsilon e, with epsilon the + dielectric permittivity). + + Since Curl e = 0, e can be derived from a scalar electric potential v, such + that e = -Grad v. Plugging this potential in Gauss' law and using the + constitutive law leads to a scalar Poisson equation in terms of the electric + potential: -Div(epsilon Grad v) = rho. + + We consider here the special case where rho = 0, to model a conducting + microstrip line on top of a dielectric substrate. A Dirichlet boundary + condition sets the potential to 1 mV on the boundary of the microstrip line + (called "Electrode" below) and 0 V on the ground. A homogeneous Neumann + boundary condition (zero flux of the displacement field, i.e. n.d = 0) is + imposed on the left boundary of the domain to account for the symmetry of the + problem, as well as on the top and right boundaries that truncate the + modelling domain. */ + +Group { + /* One starts by giving explicit meaningful names to the Physical regions + defined in the "microstrip.msh" mesh file. We only use 2 volume regions and + 2 surface regions in this model. */ + + Air = Region[101]; + Diel1 = Region[111]; + + Ground = Region[120]; + Electrode = Region[121]; + + /* We now define abstract regions to be used below in the definition of the + scalar electric potential formulation: + + Vol_Ele : volume where -Div(epsilon Grad v) = 0 is solved + Sur_Neumann_Ele: surface where non homogeneous Neumann boundary conditions + (on n.d == epsilon n.Grad v) are imposed + + Vol_xxx groups contain only volume elements of the mesh (triangles here). + Sur_xxx groups contain only surface elements of the mesh (lines here). + + Since there are no non-homogeneous Neumann conditions in the model, + Sur_Neumann_Ele is defined as empty. + */ + + Vol_Ele = Region[ {Air, Diel1} ]; + Sur_Neumann_Ele = Region[ {} ]; +} + +Function { + /* Material laws (here the dielectric permittivity) are defined as piecewise + functions (note the square brackets), in terms of the above defined + physical regions. */ + + eps0 = 8.854187818e-12; + epsilon[Air] = 1. * eps0; + epsilon[Diel1] = 9.8 * eps0; +} + +Constraint { + /* As for material laws, the Dirichlet boundary condition is defined + piecewise. The constraint "Dirichlet_Ele" is invoked in the FunctionSpace + below. */ + + { Name Dirichlet_Ele; Type Assign; + Case { + { Region Ground; Value 0.; } + { Region Electrode; Value 1.e-3; } + } + } +} + +Group{ + /* The domain of definition of a FunctionSpace lists all regions on which a + field is defined. + + Contrary to Vol_xxx and Sur_xxx regions, which contain only volume or + surface regions, resp., domains of definition Dom_xxx regions may contain + both volume and surface regions. Hence the use of the prefixes Vol_, Sur_ + and Dom_ to avoid confusions.*/ + + Dom_Hgrad_v_Ele = Region[ {Vol_Ele, Sur_Neumann_Ele} ]; +} + +FunctionSpace { + /* The function space in which we shall pick the electric scalar potential "v" + solution is definied by + - a domain of definition (the "Support": "Dom_Hgrad_v_Ele") + - a type ("Form0" means scalar field) + - a set of scalar basis functions ("BF_Node" means nodal basis functions) + - a set of entities to which the basis functions are associated ("Entity": + here all the nodes of the domain of definition "NodesOf[All]") + - a constraint (here the Dirichlet boundary conditions) + + The FE expansion of the unknown field "v" reads + + v(x,y) = Sum_k vn_k sn_k(x,y) + + where the "vn_k" are the nodal values (connectors) and "sn_k(x,y)" the + nodal basis functions. Not all connectors are unknowns of the FE problem, + due to the "Constraint", which assigns particular values to the nodes of + the Ground and Electrode regions. GetDP deals with that automatically on + basis of the definition of the FunctionSpace. */ + + { Name Hgrad_v_Ele; Type Form0; + BasisFunction { + { Name sn; NameOfCoef vn; Function BF_Node; + Support Dom_Hgrad_v_Ele; Entity NodesOf[ All ]; } + // using "NodesOf[All]" instead of "NodesOf[Dom_Hgrad_v_Ele]" is an + // optimization, which allows GetDP to not explicitly build the list of + // all nodes + } + Constraint { + { NameOfCoef vn; EntityType NodesOf; NameOfConstraint Dirichlet_Ele; } + } + } +} + +Jacobian { + /* Jacobians are used to specify the mapping between elements in the mesh and + the reference elements (defined in standardized unit cells) over which + integration is performed. "Vol" represents the classical 1-to-1 mapping + between elements of identical geometrical dimension, i.e. in this case a + reference triangle/quadrangle onto triangles/quadrangles in the z=0 plane + (2D <-> 2D). "Sur" would be used to map the reference triangle/quadrangle + onto triangles/quadrangles in a 3D space (2D <-> 3D), or to map the + reference line segment onto segments in 2D space (1D <-> 2D). "Lin" would + be used to map the reference line segment onto segments in 3D space (1D <-> + 3D). */ + + { Name Vol ; + Case { + { Region All ; Jacobian Vol ; } + } + } +} + +Integration { + /* A Gauss quadrature rule with 4 points is used for all integrations below. */ + + { Name Int ; + Case { { Type Gauss ; + Case { { GeoElement Triangle ; NumberOfPoints 4 ; } + { GeoElement Quadrangle ; NumberOfPoints 4 ; } } + } + } + } +} + +Formulation { + /* The syntax of the Formulation section is a harder nut to crack. So let's + deal with it carefully. + + A GetDP Formulation encodes a so-called weak formulation of the original + partial differential equation, i.e. of -Div(epsilon Grad v) = 0. This weak + formulation involves finding v such that + + (-Div(epsilon Grad v) , v')_Vol_Ele = 0 + + holds for all so-called "test-functions" v', where (.,.)_D denotes an inner + product over the domain D. If the test-functions v' are differentiable, + integration by parts using Green's identity leads to finding v such that + + (epsilon Grad v, Grad v')_Vol_Ele + (epsilon n.Grad v, v')_Bnd_Vol_Ele = 0 + + holds for all v', where Bnd_Vol_Ele is the boundary of Vol_Ele. In our + microstrip example this surface term vanishes, as either there is no test + function v' (on the Dirichlet boundary) or the "epsilon n.Grad v" is zero + (on the homogeneous Neumann boundary). We are thus eventually looking for + functions v in the function space Hgrad_v_Ele such that + + (epsilon Grad v, Grad v')_Vol_Ele = 0 + + holds for all v'. Finally, our choice here is to use a Gakerkin method, + where the test functions v' are the same basis functions ("sn_k") as the + ones used to interpolate the unknown potential v. + + The "Integral" statement in the Formulation is a symbolic representation of + this weak formulation. It has got 4 semicolon separated arguments: + * the density [.,.] to be integrated (note the square brackets instead of + the parentheses), with the test-functions (always) after the comma + * the domain of integration, + * the Jacobian of the transformation reference element <-> real element, + * the integration method to be used. + + In the density, braces around a quantity (such as "{v}") indicate that this + quantity belongs to a FunctionSpace. Differential operators can be applied + within braces (such as "{Grad v}"); in particular the symbol "d" represents + the exterior derivative, and it is a synonym of "Grad" when applied to a + scalar function (such as "{d v}"). + + What can be a bit confusing is that the two comma-separated terms of the + density are not interpreted exactly in the same way. Let us unravel this in + detail. + + As the Galerkin method uses as test functions the same basis functions + "sn_k" as for the unknown field "v", the second term in the density should + be something like this [ ... , basis_functions_of {d v} ]. However, since + the second term is always devoted to test functions, the operator + "basis_functions_of" would always be there. It can therefore be made + implicit, and, according to the GetDP syntax, it is omitted. So, one writes + simply [ ... , {d v} ]. + + The first term, on the other hand, can contain a much wider variety of + expressions than the second one. In our case it should be expressed in + terms of the FE expansion of "v" at the present system solution, i.e. when + the coefficients vn_k in the expansion of "v = Sum_k vn_k sn_k" are + unknown. This is indicated by prefixing the braces with "Dof", which leads + to the following density: + + [ epsilon[] * Dof{d v} , {d v} ], + + a so-called bilinear term that will contribute to the stiffness matrix of + the electrostatic problem at hand. + + Another option, which would not work here, is to evaluate the first + argument with the last available already computed solution, i.e. simply + perform the interpolation with known coefficients vn_k. For this the Dof + prefix operator would be omitted and we would have: + + [ epsilon[] * {d v} , {d v} ], + + a so-called linear term that will contribute to the right-hand side of the + linear system. + + Both choices are commonly used in different contexts, and we shall come + back on this often in subsequent tutorials. */ + { Name Electrostatics_v; Type FemEquation; + Quantity { + { Name v; Type Local; NameOfSpace Hgrad_v_Ele; } + } + Equation { + Integral { [ epsilon[] * Dof{d v} , {d v} ]; + In Vol_Ele; Jacobian Vol; Integration Int; } + } + } +} + +Resolution { + { Name EleSta_v; + System { + { Name Sys_Ele; NameOfFormulation Electrostatics_v; } + } + Operation { + Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele]; + } + } +} + +/* Post-processing is done in two parts. + + The first part defines, in terms of the Formulation, which itself refers to + the FunctionSpace, a number of quantities that can be evaluated at the + postprocessing stage. The three quantities defined here are: + - the scalar vector potential, + - the electric field, + - the electric displacement. + + The second part consists in defining groups of post-processing operations, + which can be invoked separately. The first group is invoked by default when + Gmsh is run interactively. Each Operation specifies + - a quantity to be eveluated, + - the region on which the evaluation is done, + - the name of the output file. + The generated post-processing files are automatically displayed by Gmsh if + the "Merge result automatically" option is enabled (which is the default). */ + +PostProcessing { + { Name EleSta_v; NameOfFormulation Electrostatics_v; + Quantity { + { Name v; Value { + Local { [ {v} ]; In Dom_Hgrad_v_Ele; Jacobian Vol; } + } + } + { Name e; Value { + Local { [ -{d v} ]; In Dom_Hgrad_v_Ele; Jacobian Vol; } + } + } + { Name d; Value { + Local { [ -epsilon[] * {d v} ]; In Dom_Hgrad_v_Ele; Jacobian Vol; } + } + } + } + } +} + +e = 1.e-7; // tolerance to ensure that the cut is inside the simulation domain +h = 2.e-3; // vertical position of the cut + +PostOperation { + { Name Map; NameOfPostProcessing EleSta_v; + Operation { + Print [ v, OnElementsOf Dom_Hgrad_v_Ele, File "mStrip_v.pos" ]; + Print [ e, OnElementsOf Dom_Hgrad_v_Ele, File "mStrip_e.pos" ]; + Print [ e, OnLine {{e,h,0}{14.e-3,h,0}}{60}, File "Cut_e.pos" ]; + } + } + { Name Cut; NameOfPostProcessing EleSta_v; + // same cut as above, with more points and exported in raw text format + Operation { + Print [ e, OnLine {{e,e,0}{14.e-3,e,0}} {500}, Format TimeTable, File "Cut_e.txt" ]; + } + } +} -- GitLab