# 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



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)


and solve the standard problem:

[6]:

u,v = fes.TnT()
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.
|  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
|  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)
|
|      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
|
|  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)
|
|      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)
|
|      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)
|
|      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

[ ]: