This page was generated from unit-9.1-C++FE/CppExtension.ipynb.
9.1 Implementation of Finite ElementsΒΆ
This tutorial shows how to implement our own finite elements in C++, and how to use them within the NGSolve language. We implement first order and second order triangular finite elements.
Finite elements implement the basis functions: myElement.hpp myElement.cpp
Differential operators define the mapping from coefficients to function values: myDiffOp.hpp
Finite element spaces implement the enumeration of degrees of freedom, and creation of elements: myFESpace.hpp myFESpace.cpp
Everything is included from the main file mymodule.cpp
We read the main cpp-File to a string, and let NGSolve call the compiler, and load the new library as a Python module:
[1]:
from ngsolve.fem import CompilePythonModule
from pathlib import Path
m = CompilePythonModule(Path('mymodule.cpp'), init_function_name='mymodule')
called mymodule
called ExportMyFESpace
[2]:
from netgen.occ import unit_square
from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2, quad_dominated=False))
We can now create an instance of our own finite element space:
[3]:
fes = m.MyFESpace(mesh, secondorder=True, dirichlet="left|bottom|top")
Constructor of MyFESpace
Flags = secondorder = 1
secondorder
dirichlet = 0: 1
1: 3
2: 4
You have chosen second order elements
Update MyFESpace, #vert = 39, #edge = 94
and use it within NGSolve such as the builtin finite element spaces:
[4]:
print ("ndof = ", fes.ndof)
ndof = 133
[5]:
gfu = GridFunction(fes)
gfu.Set(x*y)
Draw (gfu)
Draw (grad(gfu)[0], mesh);
and solve the standard problem:
[6]:
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
f = LinearForm(1*v*dx).Assemble()
gfu.vec.data = a.mat.Inverse(fes.FreeDofs())*f.vec
Draw (gfu);
Draw basis functions:
[7]:
gfu.vec[:] = 0
gfu.vec[mesh.nv-3] = 1
gfu.vec[fes.ndof-1] = 1
Draw (gfu, order=2);
Documentation provided in the DocInfo structure is available in Python - help. Look for 'secondorder'.
[8]:
help (m.MyFESpace)
Help on class MyFESpace:
class MyFESpace(ngsolve.comp.FESpace)
| My own FESpace.
|
| My own FESpace provides first and second order triangular elements.
|
| Keyword arguments can be:
|
| order: int = 1
| order of finite element space
| complex: bool = False
| Set if FESpace should be complex
| dirichlet: regexpr
| Regular expression string defining the dirichlet boundary.
| More than one boundary can be combined by the | operator,
| i.e.: dirichlet = 'top|right'
| dirichlet_bbnd: regexpr
| Regular expression string defining the dirichlet bboundary,
| i.e. points in 2D and edges in 3D.
| More than one boundary can be combined by the | operator,
| i.e.: dirichlet_bbnd = 'top|right'
| dirichlet_bbbnd: regexpr
| Regular expression string defining the dirichlet bbboundary,
| i.e. points in 3D.
| More than one boundary can be combined by the | operator,
| i.e.: dirichlet_bbbnd = 'top|right'
| definedon: Region or regexpr
| FESpace is only defined on specific Region, created with mesh.Materials('regexpr')
| or mesh.Boundaries('regexpr'). If given a regexpr, the region is assumed to be
| mesh.Materials('regexpr').
| dim: int = 1
| Create multi dimensional FESpace (i.e. [H1]^3)
| dgjumps: bool = False
| Enable discontinuous space for DG methods, this flag is needed for DG methods,
| since the dofs have a different coupling then and this changes the sparsity
| pattern of matrices.
| autoupdate: bool = False
| Automatically update on a change to the mesh.
| low_order_space: bool = True
| Generate a lowest order space together with the high-order space,
| needed for some preconditioners.
| hoprolongation: bool = False
| Create high order prolongation operators,
| only available for H1 and L2 on simplicial meshes
| order_policy: ORDER_POLICY = ORDER_POLICY.OLDSTYLE
| CONSTANT .. use the same fixed order for all elements,
| NODAL ..... use the same order for nodes of same shape,
| VARIABLE ... use an individual order for each edge, face and cell,
| OLDSTYLE .. as it used to be for the last decade
| print: bool = False
| (historic) print some output into file set by 'SetTestoutFile'
| secondorder: bool = False
| Use second order basis functions
|
| Method resolution order:
| MyFESpace
| ngsolve.comp.FESpace
| ngsolve.comp.NGS_Object
| pybind11_builtins.pybind11_object
| builtins.object
|
| Methods defined here:
|
| GetNVert(...)
| GetNVert(self: .MyFESpace) -> int
|
| return number of vertices
|
| __getstate__(...)
| __getstate__(self: ngsolve.comp.FESpace) -> tuple
|
| __init__(...)
| __init__(self: .MyFESpace, mesh: ngsolve.comp.Mesh, **kwargs) -> None
|
| __setstate__(...)
| __setstate__(self: .MyFESpace, arg0: tuple) -> None
|
| ----------------------------------------------------------------------
| Static methods defined here:
|
| __flags_doc__(...) from builtins.PyCapsule
| __flags_doc__() -> dict
|
| ----------------------------------------------------------------------
| Data descriptors defined here:
|
| __dict__
|
| ----------------------------------------------------------------------
| Data and other attributes defined here:
|
| __pybind11_module_local_v4_gcc_libstdcpp_cxxabi1016__ = <capsule objec...
|
| ----------------------------------------------------------------------
| Methods inherited from ngsolve.comp.FESpace:
|
| ApplyM(...)
| ApplyM(self: ngsolve.comp.FESpace, vec: ngsolve.la.BaseVector, rho: ngsolve.fem.CoefficientFunction = None, definedon: ngsolve.comp.Region = None) -> None
|
| Apply mass-matrix. Available only for L2-like spaces
|
| ConvertL2Operator(...)
| ConvertL2Operator(self: ngsolve.comp.FESpace, l2space: ngsolve.comp.FESpace) -> BaseMatrix
|
| CouplingType(...)
| CouplingType(self: ngsolve.comp.FESpace, dofnr: int) -> ngsolve.comp.COUPLING_TYPE
|
|
| Get coupling type of a degree of freedom.
|
| Parameters:
|
| dofnr : int
| input dof number
|
| CreateDirectSolverCluster(...)
| CreateDirectSolverCluster(self: ngsolve.comp.FESpace, **kwargs) -> list
|
| CreateSmoothingBlocks(...)
| CreateSmoothingBlocks(self: ngsolve.comp.FESpace, **kwargs) -> pyngcore.pyngcore.Table_I
|
|
| Create table of smoothing blocks for block-Jacobi/block-Gauss-Seidel preconditioners.
|
| Every table entry describes the set of dofs belonging a Jacobi/Gauss-Seidel block.
|
| Paramters:
|
| blocktype: string | [ string ] | int
| describes blocktype.
| string form ["vertex", "edge", "face", "facet", "vertexedge", ....]
| or list of strings for combining multiple blocktypes
| int is for backward compatibility with old style blocktypes
|
| condense: bool = False
| boolexclude dofs eliminated by static condensation
|
| Elements(...)
| Elements(*args, **kwargs)
| Overloaded function.
|
| 1. Elements(self: ngsolve.comp.FESpace, VOL_or_BND: ngsolve.comp.VorB = <VorB.VOL: 0>) -> ngsolve.comp.FESpaceElementRange
|
|
| Returns an iterable range of elements.
|
| Parameters:
|
| VOL_or_BND : ngsolve.comp.VorB
| input VOL, BND, BBND,...
|
|
|
| 2. Elements(self: ngsolve.comp.FESpace, arg0: ngsolve.comp.Region) -> Iterator[ngsolve.comp.FESpaceElement]
|
| FinalizeUpdate(...)
| FinalizeUpdate(self: ngsolve.comp.FESpace) -> None
|
| finalize update
|
| FreeDofs(...)
| FreeDofs(self: ngsolve.comp.FESpace, coupling: bool = False) -> pyngcore.pyngcore.BitArray
|
|
|
| Return BitArray of free (non-Dirichlet) dofs\n
| coupling=False ... all free dofs including local dofs\n
| coupling=True .... only element-boundary free dofs
|
| Parameters:
|
| coupling : bool
| input coupling
|
| GetDofNrs(...)
| GetDofNrs(*args, **kwargs)
| Overloaded function.
|
| 1. GetDofNrs(self: ngsolve.comp.FESpace, ei: ngsolve.comp.ElementId) -> tuple
|
|
|
| Parameters:
|
| ei : ngsolve.comp.ElementId
| input element id
|
|
|
| 2. GetDofNrs(self: ngsolve.comp.FESpace, ni: ngsolve.comp.NodeId) -> tuple
|
|
|
| Parameters:
|
| ni : ngsolve.comp.NodeId
| input node id
|
| GetDofs(...)
| GetDofs(self: ngsolve.comp.FESpace, region: ngsolve.comp.Region) -> pyngcore.pyngcore.BitArray
|
|
| Returns all degrees of freedom in given region.
|
| Parameters:
|
| region : ngsolve.comp.Region
| input region
|
| GetFE(...)
| GetFE(self: ngsolve.comp.FESpace, ei: ngsolve.comp.ElementId) -> object
|
|
| Get the finite element to corresponding element id.
|
| Parameters:
|
| ei : ngsolve.comp.ElementId
| input element id
|
| GetOrder(...)
| GetOrder(self: ngsolve.comp.FESpace, nodeid: ngsolve.comp.NodeId) -> int
|
| return order of node.
| by now, only isotropic order is supported here
|
| GetTrace(...)
| GetTrace(self: ngsolve.comp.FESpace, arg0: ngsolve.comp.FESpace, arg1: ngsolve.la.BaseVector, arg2: ngsolve.la.BaseVector, arg3: bool) -> None
|
| GetTraceTrans(...)
| GetTraceTrans(self: ngsolve.comp.FESpace, arg0: ngsolve.comp.FESpace, arg1: ngsolve.la.BaseVector, arg2: ngsolve.la.BaseVector, arg3: bool) -> None
|
| HideAllDofs(...)
| HideAllDofs(self: ngsolve.comp.FESpace, component: object = <ngsolve.ngstd.DummyArgument>) -> None
|
| set all visible coupling types to HIDDEN_DOFs (will be overwritten by any Update())
|
| InvM(...)
| InvM(self: ngsolve.comp.FESpace, rho: ngsolve.fem.CoefficientFunction = None) -> BaseMatrix
|
| Mass(...)
| Mass(self: ngsolve.comp.FESpace, rho: ngsolve.fem.CoefficientFunction = None, definedon: Optional[ngsolve.comp.Region] = None) -> BaseMatrix
|
| ParallelDofs(...)
| ParallelDofs(self: ngsolve.comp.FESpace) -> ngsolve.la.ParallelDofs
|
| Return dof-identification for MPI-distributed meshes
|
| Prolongation(...)
| Prolongation(self: ngsolve.comp.FESpace) -> ngmg::Prolongation
|
| Return prolongation operator for use in multi-grid
|
| Range(...)
| Range(self: ngsolve.comp.FESpace, arg0: int) -> ngsolve.la.DofRange
|
| deprecated, will be only available for ProductSpace
|
| SetCouplingType(...)
| SetCouplingType(*args, **kwargs)
| Overloaded function.
|
| 1. SetCouplingType(self: ngsolve.comp.FESpace, dofnr: int, coupling_type: ngsolve.comp.COUPLING_TYPE) -> None
|
|
| Set coupling type of a degree of freedom.
|
| Parameters:
|
| dofnr : int
| input dof number
|
| coupling_type : ngsolve.comp.COUPLING_TYPE
| input coupling type
|
|
|
| 2. SetCouplingType(self: ngsolve.comp.FESpace, dofnrs: ngsolve.ngstd.IntRange, coupling_type: ngsolve.comp.COUPLING_TYPE) -> None
|
|
| Set coupling type for interval of dofs.
|
| Parameters:
|
| dofnrs : Range
| range of dofs
|
| coupling_type : ngsolve.comp.COUPLING_TYPE
| input coupling type
|
| SetDefinedOn(...)
| SetDefinedOn(self: ngsolve.comp.FESpace, region: ngsolve.comp.Region) -> None
|
|
| Set the regions on which the FESpace is defined.
|
| Parameters:
|
| region : ngsolve.comp.Region
| input region
|
| SetOrder(...)
| SetOrder(*args, **kwargs)
| Overloaded function.
|
| 1. SetOrder(self: ngsolve.comp.FESpace, element_type: ngsolve.fem.ET, order: int) -> None
|
|
|
| Parameters:
|
| element_type : ngsolve.fem.ET
| input element type
|
| order : object
| input polynomial order
|
|
| 2. SetOrder(self: ngsolve.comp.FESpace, nodeid: ngsolve.comp.NodeId, order: int) -> None
|
|
|
| Parameters:
|
| nodeid : ngsolve.comp.NodeId
| input node id
|
| order : int
| input polynomial order
|
| SolveM(...)
| SolveM(self: ngsolve.comp.FESpace, vec: ngsolve.la.BaseVector, rho: ngsolve.fem.CoefficientFunction = None, definedon: ngsolve.comp.Region = None) -> None
|
|
| Solve with the mass-matrix. Available only for L2-like spaces.
|
| Parameters:
|
| vec : ngsolve.la.BaseVector
| input right hand side vector
|
| rho : ngsolve.fem.CoefficientFunction
| input CF
|
| TestFunction(...)
| TestFunction(self: ngsolve.comp.FESpace) -> object
|
| Return a proxy to be used as a testfunction for :any:`Symbolic Integrators<symbolic-integrators>`
|
| TnT(...)
| TnT(self: ngsolve.comp.FESpace) -> tuple[object, object]
|
| Return a tuple of trial and testfunction
|
| TraceOperator(...)
| TraceOperator(self: ngsolve.comp.FESpace, tracespace: ngsolve.comp.FESpace, average: bool) -> BaseMatrix
|
| TrialFunction(...)
| TrialFunction(self: ngsolve.comp.FESpace) -> object
|
| Return a proxy to be used as a trialfunction in :any:`Symbolic Integrators<symbolic-integrators>`
|
| Update(...)
| Update(self: ngsolve.comp.FESpace) -> None
|
| update space after mesh-refinement
|
| UpdateDofTables(...)
| UpdateDofTables(self: ngsolve.comp.FESpace) -> None
|
| update dof-tables after changing polynomial order distribution
|
| __eq__(...)
| __eq__(self: ngsolve.comp.FESpace, space: ngsolve.comp.FESpace) -> bool
|
| __mul__(...)
| __mul__(self: ngsolve.comp.FESpace, arg0: ngsolve.comp.FESpace) -> ngcomp::CompoundFESpace
|
| __pow__(...)
| __pow__(self: ngsolve.comp.FESpace, arg0: int) -> ngcomp::CompoundFESpaceAllSame
|
| __str__(...)
| __str__(self: ngsolve.comp.FESpace) -> str
|
| __timing__(...)
| __timing__(self: ngsolve.comp.FESpace) -> object
|
| ----------------------------------------------------------------------
| Static methods inherited from ngsolve.comp.FESpace:
|
| __special_treated_flags__(...) from builtins.PyCapsule
| __special_treated_flags__() -> dict
|
| ----------------------------------------------------------------------
| Readonly properties inherited from ngsolve.comp.FESpace:
|
| autoupdate
|
| components
| deprecated, will be only available for ProductSpace
|
| couplingtype
|
| dim
| multi-dim of FESpace
|
| globalorder
| query global order of space
|
| is_complex
|
| loembedding
|
| lospace
|
| mesh
| mesh on which the FESpace is created
|
| ndof
| number of degrees of freedom
|
| ndofglobal
| global number of dofs on MPI-distributed mesh
|
| type
| type of finite element space
|
| ----------------------------------------------------------------------
| Data and other attributes inherited from ngsolve.comp.FESpace:
|
| __hash__ = None
|
| ----------------------------------------------------------------------
| Readonly properties inherited from ngsolve.comp.NGS_Object:
|
| __memory__
|
| flags
|
| ----------------------------------------------------------------------
| Data descriptors inherited from ngsolve.comp.NGS_Object:
|
| name
|
| ----------------------------------------------------------------------
| Static methods inherited from pybind11_builtins.pybind11_object:
|
| __new__(*args, **kwargs) from pybind11_builtins.pybind11_type
| Create and return a new object. See help(type) for accurate signature.
Exercises:
Extend MyFESpace by the following 1st and 2nd order elements:
1D finite elements (ET_SEGM), as needed for boundary conditions
quadrilateral elements (ET_QUAD), use geom.GenerateMesh(quad_dominated=True)
tetrahedral elements (ET_TET), test it for 3D domains
Next, implement 3rd order elements
[ ]: