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.

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

[ ]: