2.5. A Slighly More Complex Problem#
Let’s create now a capacitor problem in which we prescribe the potential at the electrodes. The equation that we need to solve is the following:
Find \(u\) such that
\[\begin{align*} -\varepsilon \Delta u &= 0 \quad \text{in } \Omega, \\ u &= 1 \quad \text{on } \text{Electrode 1},\\ u &= -1 \quad \text{on } \text{Electrode 2},\\ n\cdot \nabla u &= 0 \quad \text{on } \text{outer}, \end{align*}\]
With permettivity given by:
\[\begin{align*} \varepsilon = \begin{cases} 1 & \text{in air},\\ 2 & \text{in dielectric}. \end{cases} \end{align*}\]
And domain as in figure below:
from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
square = Circle( (0,0), 2).Face() # create a rectangle with boundary conditions
square.edges.name = "outer"
square.faces.name = "air"
el1 = MoveTo(-0.4, 0.2).Rectangle( 0.8, 0.1).Face() # create a rectangle with boundary conditions
el1.edges.name = "el1"
el1.vertices.name = "el1"
el2 = MoveTo(-0.4, -0.2).Rectangle( 0.8, 0.1).Face() # create a rectangle with boundary conditions
el2.edges.name = "el2"
el2.vertices.name = "el2"
dielec = MoveTo(-0.9, -0.1).Rectangle( 1.8, 0.3).Face()
dielec.faces.name = "dielec"
air = square # subtract the rectangles from the air rectangle
shape = Glue([air - dielec, dielec])
shape =shape - el1 - el2
### adding extra specifications on the shape
#predefined mesh size for the shape
shape.edges["el.*"].maxh = 0.1 # set the mesh size of the edges with the name "el.*" to 0.05
shape.vertices["el.*"].maxh = 0.05 # set the mesh size of the vertices with the name "el.*" to 0.05
Draw(shape); # draw the shape
mesh = Mesh(OCCGeometry(shape, dim = 2).GenerateMesh(maxh=0.5)).Curve(4) # create a mesh from the shape
eps = CoefficientFunction( [1, 2] ) # define the coefficient function # try sin(5*x) instead of 2
Draw(eps, mesh, "eps"); # draw the coefficient function
help(mesh)
Help on Mesh in module ngsolve.comp object:
class Mesh(pybind11_builtins.pybind11_object)
| NGSolve interface to the Netgen mesh. Provides access and functionality
| to use the mesh for finite element calculations.
|
| Parameters:
|
| mesh (netgen.Mesh): a mesh generated from Netgen
|
| Method resolution order:
| Mesh
| pybind11_builtins.pybind11_object
| builtins.object
|
| Methods defined here:
|
| BBBoundaries(...)
| BBBoundaries(self: ngsolve.comp.Mesh, pattern: str) -> ngsolve.comp.Region
|
| Return co dim 3 boundary mesh-region matching the given regex pattern
|
| BBoundaries(...)
| BBoundaries(self: ngsolve.comp.Mesh, pattern: str) -> ngsolve.comp.Region
|
| Return co dim 2 boundary mesh-region matching the given regex pattern
|
| Boundaries(...)
| Boundaries(*args, **kwargs)
| Overloaded function.
|
| 1. Boundaries(self: ngsolve.comp.Mesh, pattern: str) -> ngsolve.comp.Region
|
| Return boundary mesh-region matching the given regex pattern
|
| 2. Boundaries(self: ngsolve.comp.Mesh, bnds: list[int]) -> ngsolve.comp.Region
|
| Generate boundary mesh-region by boundary condition numbers
|
| BoundaryCF(...)
| BoundaryCF(self: ngsolve.comp.Mesh, values: dict, default: ngsolve.fem.CoefficientFunction = None) -> ngsolve.fem.CoefficientFunction
|
| Boundary wise CoefficientFunction.
| First argument is a dict from either boundary names or Region objects to
| CoefficientFunction (-values). Later given names/regions override earlier
| values. Optional last argument (default) is the value for not given boundaries.
| >>> penalty = mesh.BoundaryCF({ "top" : 1e6 }, default=0)
| will create a CF being 1e6 on the top boundary and 0. elsewhere.
|
| BuildRefinementTree(...)
| BuildRefinementTree(self: ngsolve.comp.Mesh) -> pyngcore.pyngcore.Array_S_S
|
| Contains(...)
| Contains(self: ngsolve.comp.Mesh, x: float = 0.0, y: float = 0.0, z: float = 0.0) -> bool
|
| Check if the point (x,y,z) is in the meshed domain (is inside a volume element)
|
| Curve(...)
| Curve(self: ngsolve.comp.Mesh, order: int) -> ngsolve.comp.Mesh
|
| Curve the mesh elements for geometry approximation of given order
|
| Elements(...)
| Elements(self: ngsolve.comp.Mesh, VOL_or_BND: ngsolve.comp.VorB = <VorB.VOL: 0>) -> ngcomp::ElementRange
|
| Return an iterator over elements on VOL/BND
|
| GeoParamCF(...)
| GeoParamCF(self: ngsolve.comp.Mesh) -> ngsolve.fem.CoefficientFunction
|
| GetBBBoundaries(...)
| GetBBBoundaries(self: ngsolve.comp.Mesh) -> tuple
|
| Return list of boundary conditions for co dimension 3
|
| GetBBoundaries(...)
| GetBBoundaries(self: ngsolve.comp.Mesh) -> tuple
|
| Return list of boundary conditions for co dimension 2
|
| GetBoundaries(...)
| GetBoundaries(self: ngsolve.comp.Mesh) -> tuple
|
| Return list of boundary condition names
|
| GetCurveOrder(...)
| GetCurveOrder(self: ngsolve.comp.Mesh) -> int
|
| GetHPElementLevel(...)
| GetHPElementLevel(self: ngsolve.comp.Mesh, ei: ngsolve.comp.ElementId) -> int
|
| THIS FUNCTION IS WIP!
| Return HP-refinement level of element
|
| GetMaterials(...)
| GetMaterials(self: ngsolve.comp.Mesh) -> tuple
|
| Return list of material names
|
| GetNE(...)
| GetNE(self: ngsolve.comp.Mesh, arg0: ngsolve.comp.VorB) -> int
|
| Return number of elements of codimension VorB.
|
| GetPMLTrafo(...)
| GetPMLTrafo(self: ngsolve.comp.Mesh, dom: int = 1) -> ngsolve.comp.pml.PML
|
| Return pml transformation on domain dom
|
| GetPMLTrafos(...)
| GetPMLTrafos(self: ngsolve.comp.Mesh) -> list
|
| Return list of pml transformations
|
| GetParentElement(...)
| GetParentElement(self: ngsolve.comp.Mesh, ei: ngsolve.comp.ElementId) -> ngsolve.comp.ElementId
|
| Return parent element id on refined mesh
|
| GetParentFaces(...)
| GetParentFaces(self: ngsolve.comp.Mesh, fnum: int) -> tuple
|
| Return parent faces
|
| GetParentVertices(...)
| GetParentVertices(self: ngsolve.comp.Mesh, vnum: int) -> tuple
|
| Return parent vertex numbers on refined mesh
|
| GetPeriodicNodePairs(...)
| GetPeriodicNodePairs(self: ngsolve.comp.Mesh, arg0: ngsolve.fem.NODE_TYPE) -> list
|
| returns list of periodic nodes with their identification number as [((master_nr, minion_nr),idnr),...]
|
| GetTrafo(...)
| GetTrafo(self: ngsolve.comp.Mesh, eid: ngsolve.comp.ElementId) -> ngsolve.fem.ElementTransformation
|
| returns element transformation of given element id
|
| LocalHCF(...)
| LocalHCF(self: ngsolve.comp.Mesh) -> ngsolve.fem.CoefficientFunction
|
| MapToAllElements(...)
| MapToAllElements(*args, **kwargs)
| Overloaded function.
|
| 1. MapToAllElements(self: ngsolve.comp.Mesh, arg0: ngsolve.fem.IntegrationRule, arg1: Union[ngsolve.comp.VorB, ngsolve.comp.Region]) -> numpy.ndarray[ngsolve.fem.MeshPoint]
|
| 2. MapToAllElements(self: ngsolve.comp.Mesh, arg0: dict[ngsolve.fem.ET, ngsolve.fem.IntegrationRule], arg1: Union[ngsolve.comp.VorB, ngsolve.comp.Region]) -> numpy.ndarray[ngsolve.fem.MeshPoint]
|
| MaterialCF(...)
| MaterialCF(self: ngsolve.comp.Mesh, values: dict, default: ngsolve.fem.CoefficientFunction = None) -> ngsolve.fem.CoefficientFunction
|
| Domain wise CoefficientFunction.
| First argument is a dict from either material names or Region objects to
| CoefficientFunction (-values). Later given names/regions override earlier
| values. Optional last argument (default) is the value for not given materials.
| >>> sigma = mesh.MaterialCF({ "steel_.*" : 2e6 }, default=0)
| will create a CF being 2e6 on all domains starting with 'steel_' and 0 elsewhere.
|
| Materials(...)
| Materials(*args, **kwargs)
| Overloaded function.
|
| 1. Materials(self: ngsolve.comp.Mesh, pattern: str) -> ngsolve.comp.Region
|
| Return mesh-region matching the given regex pattern
|
| 2. Materials(self: ngsolve.comp.Mesh, domains: list[int]) -> ngsolve.comp.Region
|
| Generate mesh-region by domain numbers
|
| Refine(...)
| Refine(self: ngsolve.comp.Mesh, mark_surface_elements: bool = False, onlyonce: bool = False) -> None
|
| Local mesh refinement based on marked elements, uses element-bisection algorithm
|
| RefineFromTree(...)
| RefineFromTree(self: ngsolve.comp.Mesh, arg0: pyngcore.pyngcore.Array_S_S) -> None
|
| RefineHP(...)
| RefineHP(self: ngsolve.comp.Mesh, levels: int, factor: float = 0.125) -> None
|
| Geometric mesh refinement towards marked vertices and edges, uses factor for placement of new points
|
| Region(...)
| Region(self: ngsolve.comp.Mesh, vb: ngsolve.comp.VorB, pattern: Optional[str] = '.*') -> ngsolve.comp.Region
|
| Return boundary mesh-region matching the given regex pattern
|
| RegionCF(...)
| RegionCF(self: ngsolve.comp.Mesh, VorB: ngsolve.comp.VorB, value: dict, default: ngsolve.fem.CoefficientFunction = None) -> ngsolve.fem.CoefficientFunction
|
| Region wise CoefficientFunction.
| First argument is VorB, defining the co-dimension,
| second argument is a dict from either region names or Region objects to
| CoefficientFunction (-values). Later given names/regions override earlier
| values. Optional last argument (default) is the value for not given regions.
| >>> sigma = mesh.RegionCF(VOL, { "steel_.*" : 2e6 }, default=0)
| will create a CF being 2e6 on all domains starting with 'steel_' and 0 elsewhere.
|
| SetDeformation(...)
| SetDeformation(self: ngsolve.comp.Mesh, gf: ngcomp::GridFunction) -> None
|
| Deform the mesh with the given GridFunction
|
| SetElementOrder(...)
| SetElementOrder(self: ngsolve.comp.Mesh, eid: ngsolve.comp.ElementId, order: int) -> None
|
| For backward compatibility, not recommended to use
|
| SetPML(...)
| SetPML(self: ngsolve.comp.Mesh, pmltrafo: ngsolve.comp.pml.PML, definedon: object) -> None
|
| Set PML transformation on domain
|
| SetRefinementFlag(...)
| SetRefinementFlag(self: ngsolve.comp.Mesh, ei: ngsolve.comp.ElementId, refine: bool) -> None
|
| Set refinementflag for mesh-refinement
|
| SetRefinementFlags(...)
| SetRefinementFlags(self: ngsolve.comp.Mesh, refine: list[bool]) -> None
|
| Set refinementflags for mesh-refinement
|
| SplitElements_Alfeld(...)
| SplitElements_Alfeld(self: ngsolve.comp.Mesh) -> None
|
| UnSetPML(...)
| UnSetPML(self: ngsolve.comp.Mesh, definedon: object) -> None
|
| Unset PML transformation on domain
|
| UnsetDeformation(...)
| UnsetDeformation(self: ngsolve.comp.Mesh) -> None
|
| Unset the deformation
|
| __call__(...)
| __call__(self: ngsolve.comp.Mesh, x: numpy.ndarray[numpy.float64] = 0.0, y: numpy.ndarray[numpy.float64] = 0.0, z: numpy.ndarray[numpy.float64] = 0.0, VOL_or_BND: ngsolve.comp.VorB = <VorB.VOL: 0>) -> object
|
| Get a MappedIntegrationPoint in the point (x,y,z) on the matching volume (VorB=VOL, default) or surface (VorB=BND) element. BBND elements aren't supported
|
| __eq__(...)
| __eq__(self: ngsolve.comp.Mesh, mesh: ngsolve.comp.Mesh) -> bool
|
| __getitem__(...)
| __getitem__(*args, **kwargs)
| Overloaded function.
|
| 1. __getitem__(self: ngsolve.comp.Mesh, arg0: ngsolve.comp.ElementId) -> ngsolve.comp.Ngs_Element
|
| Return Ngs_Element from given ElementId
|
| 2. __getitem__(self: ngsolve.comp.Mesh, arg0: ngsolve.comp.NodeId) -> ngsolve.comp.MeshNode
|
| Return MeshNode from given NodeId
|
| __getstate__(...)
| __getstate__(self: ngsolve.comp.Mesh) -> tuple
|
| __init__(...)
| __init__(*args, **kwargs)
| Overloaded function.
|
| 1. __init__(self: ngsolve.comp.Mesh, ngmesh: netgen.libngpy._meshing.Mesh) -> None
|
| Make an NGSolve-mesh from a Netgen-mesh
|
| 2. __init__(self: ngsolve.comp.Mesh, filename: str, comm: pyngcore.pyngcore.MPI_Comm = <pyngcore.pyngcore.MPI_Comm object at 0x7ff90d8b5ab0>) -> None
|
| Load a mesh from file.
| In MPI-parallel mode the mesh is distributed over the MPI-group given by the communicator (WIP!)
|
| __setstate__(...)
| __setstate__(self: ngsolve.comp.Mesh, arg0: tuple) -> None
|
| nnodes(...)
| nnodes(self: ngsolve.comp.Mesh, arg0: ngsolve.fem.NODE_TYPE) -> int
|
| number of nodes given type
|
| nodes(...)
| nodes(self: ngsolve.comp.Mesh, node_type: ngsolve.fem.NODE_TYPE) -> ngsolve.comp.MeshNodeRange
|
| iterable of mesh nodes of type node_type
|
| ----------------------------------------------------------------------
| Readonly properties defined here:
|
| comm
| MPI-communicator the Mesh lives in
|
| dim
| mesh dimension
|
| edges
| iterable of mesh edges
|
| faces
| iterable of mesh faces
|
| facets
| iterable of mesh facets
|
| levels
| multigrid levels
|
| ne
| number of volume elements
|
| nedge
| number of edges
|
| nface
| number of faces
|
| nfacet
| number of facets
|
| ngmesh
| the Netgen mesh
|
| nv
| number of vertices
|
| vertices
| iterable of mesh vertices
|
| ----------------------------------------------------------------------
| Data descriptors defined here:
|
| __dict__
|
| deformation
| mesh deformation
|
| ----------------------------------------------------------------------
| Data and other attributes defined here:
|
| __hash__ = None
|
| ----------------------------------------------------------------------
| Static methods inherited from pybind11_builtins.pybind11_object:
|
| __new__(*args, **kwargs) class method of pybind11_builtins.pybind11_object
| Create and return a new object. See help(type) for accurate signature.
The bilinear and liner forms now are :
\[\begin{align*} a(u,v) &= \int_{\Omega} \varepsilon \nabla u \cdot \nabla v \, dx,\\ f(v) &= 0. \end{align*}\]
order = 3
fes = H1(mesh, order=order, dirichlet="el.*") # define the finite element space
u ,v = fes.TnT() # define the trial and test functions
a = BilinearForm(fes)
a += eps * grad(u) * grad(v) * dx
a.Assemble()
f = LinearForm(fes)
f.Assemble() # not needed it is zero
<ngsolve.comp.LinearForm at 0x7ff90d8914b0>
We dont have homogeneous Dirichlet boundary conditions, so we need a homogenization technique.
Homogenization technique in a nutshell: Suppose we the solution \(u_{\text{sol}} = u_0 + u_{D}\), where \(u_0\) is the solution of the problem with homogeneous boundary conditions and \( u_{D}\) is the solution
\[\begin{align*}
A(u_{\text{sol}}, v ) = f(v)\\
A(u_0 + u_D, v ) = f(v)\\
A(u_0, v ) = f(v)-A(u_D, v )\\
\end{align*}\]
After the discretization is equivalent to solve the following linear system:
\[\begin{align*}
u_0 = A^{-1}(f-Au_D)
\end{align*}\]
The solution can be found as
\[\begin{align*}
u_\text{sol} = A^{-1}(f-Au_{D}) + u_D
\end{align*}\]
gfu = GridFunction(fes) # define the grid function
electrode = mesh.BoundaryCF({"el1":1, "el2":-1}, default = 0 ) # define the boundary conditions
gfu.Set(electrode, definedon=mesh.Boundaries("el.*")) # set the boundary conditions
Draw(gfu, deformation =True); # draw the grid function
gfu.vec.data += a.mat.Inverse(freedofs = fes.FreeDofs()) * (f.vec - a.mat * gfu.vec) # solve the system
Draw(gfu , deformation =True, scale = 1); # draw the solution
# Electic field
Draw(-eps*grad(gfu) , mesh, "E", min = 0, max = 4); # draw the solution
# if you zoom in you will see the difference in the two materials.
2.5.1. Some visualization tools#
from ngsolve.webgui import Draw, FieldLines, AddFieldLines
fes_flux = HDiv(mesh, order=order-1) # Same as
#fes_flux = H1(mesh, order=order-1, dim=mesh.dim)
# fes_flux = HDiv(mesh, order=order-1)
flux = GridFunction(fes_flux, name="flux")
flux.Set( eps*grad(gfu) )
N = 10
#fac = 0 if mesh.dim == 2 else 1
#p = [(-1+4*i/N,-2+4*j/N,fac * 2*k/N) for i in range(1,2*N) for j in range(1,2*N) for k in range(1,N)]
p = [(i / N, -1 + (2*j)/ N, 0) for i in range(N+1) for j in range(N+1) ]
fieldlines = flux._BuildFieldLines(mesh, p, num_fieldlines=N**3//5, randomized=False, length=0.3)
clipping = { "clipping" : { "y":0, "z":-1} }
Draw(-eps*grad(gfu), mesh, "X", draw_vol=True, draw_surf=True, objects=[fieldlines], \
autoscale=True, min = 0, max = 10, settings={"Objects": {"Surface": False}}, **clipping);