NGSTrefftz

https://paulst.github.io/NGSTrefftz/
Admin
NGSTrefftz is an add-on to NGSolve that provides a framework and finite element spaces for Trefftz...
Tokens:41,086
Snippets:258
Trust Score:9.7
Update:1 month ago
Show doc for...
Tokens:
Raw
### Run NGSTrefftz examples with Docker

Source: https://paulst.github.io/NGSTrefftz/readme.html

Commands to clone the repository, build the Docker image, and start a Jupyter environment to run the example notebooks locally.

```bash
git clone https://github.com/PaulSt/NGSTrefftz
cd NGSTrefftz && docker build -t ngstrefftz_jupyter .
docker run -p 8888:8888 ngstrefftz_jupyter
```

--------------------------------

### Import NGSTrefftz and set thread count

Source: https://paulst.github.io/NGSTrefftz/notebooks/qtelliptic.html

Imports ngsolve, ngstrefftz, and netgen.occ, then sets the number of threads to 4. This is the initial setup for the quasi-Trefftz DG example.

```python
from ngsolve import *
from ngstrefftz import *
from netgen.occ import *
SetNumThreads(4)

```

--------------------------------

### Define problem parameters and mesh

Source: https://paulst.github.io/NGSTrefftz/notebooks/embTrefftz-adv.html

Setup of the velocity field, exact solution, right-hand side, and the mesh for the unit square.

```python
b = CF((-sin(y),cos(x)))
exact = sin(3*x)*sin(3*y)
rhs = b*CF((exact.Diff(x),exact.Diff(y)))

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
fes = L2(mesh, order=4,  dgjumps=True)
u,v = fes.TnT()
```

--------------------------------

### Import NGSTrefftz and NGSolve modules

Source: https://paulst.github.io/NGSTrefftz/notebooks/helmholtz.html

Imports the necessary modules for the Helmholtz Trefftz-DG example.

```python
from ngsolve import *
from ngstrefftz import *
from netgen.occ import *
```

--------------------------------

### Install NGSTrefftz via pip

Source: https://paulst.github.io/NGSTrefftz/readme.html

The recommended method for installing the NGSTrefftz package.

```bash
pip install ngstrefftz
```

--------------------------------

### Import libraries and compute exact solution

Source: https://paulst.github.io/NGSTrefftz/notebooks/embTrefftz-wave.html

Imports NGSolve, ngstrefftz, and Netgen OCC, and defines the exact solution for the wave equation as sinusoidal functions with a wave speed c=1. This setup is necessary to later compare the numerical error.

```python
from ngsolve import *
from ngstrefftz import *
from netgen.occ import *
from ngsolve.TensorProductTools import MakeTensorProductMesh, SegMesh
from ngsolve import *
from ngsolve.webgui import Draw

```

```python
from numpy import pi
c=1
sq = sqrt(2.0)
sig_ex = CoefficientFunction((-cos(pi*x)*sin(pi*y)*sin(pi*z*c*sq)/sq,
    -sin(pi*x)*cos(pi*y)*sin(pi*z*c*sq)/sq))
v_ex = sin(pi*x)*sin(pi*y)*cos(pi*z*c*sq)*c
gD = v_ex

```

--------------------------------

### Run the test suite

Source: https://paulst.github.io/NGSTrefftz/contrib.html

Execute this command within your build directory to verify the project's test suite.

```bash
make test
```

--------------------------------

### Define parameters, exact solution, and boundary condition

Source: https://paulst.github.io/NGSTrefftz/notebooks/helmholtz.html

Sets the frequency omega, defines the exact plane wave solution and its gradient, computes the Robin boundary condition, and sets a small epsilon value.

```python
omega=1
n = specialcf.normal(2)
exact = exp(1j*sqrt(0.5)*(x+y))
gradexact = CoefficientFunction((sqrt(0.5)*1j*exact, sqrt(0.5)*1j*exact))
bndc = gradexact*n + 1j*omega*exact
eps = 10**-7
```

--------------------------------

### Set up mesh and visualize basis functions

Source: https://paulst.github.io/NGSTrefftz/notebooks/embTrefftz-stokes.html

Sets up the mesh and finite element spaces for the Stokes problem, calls the TrefftzStokesEmb function, and visualizes the basis functions using the Draw function.

```python
order = 3
mesh = Mesh(unit_square.GenerateMesh(maxh=1))

V = VectorL2(mesh, order=order, dgjumps=True)
Q = L2(mesh, order=order - 1, dgjumps=True)
fes = V * Q
PP,uf = TrefftzStokesEmb(fes)

gfshow = EmbeddedBasisFunctionsAsMultiDim(PP,fes)
Draw (gfshow.components[0], mesh, interpolate_multidim=False, animate=False, autoscale=True, vectors = {"grid_size" : 20}, settings={"subdivision":20})
```

--------------------------------

### Assemble SIP-DG bilinear and linear forms

Source: https://paulst.github.io/NGSTrefftz/notebooks/embTrefftz-poi.html

Sets up the SIP-DG formulation for the Poisson equation, including jump and mean terms for the skeleton and boundary.

```python
exact = sin(x)*sin(y)*sin(z)
rhs = 3*sin(x)*sin(y)*sin(z)
bndc=exact
alpha = 4
n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size
u = fes.TrialFunction()
v = fes.TestFunction()

jump_u = (u-u.Other())*n
jump_v = (v-v.Other())*n
mean_dudn = 0.5 * (grad(u)+grad(u.Other()))
mean_dvdn = 0.5 * (grad(v)+grad(v.TestFunction().Other()))

a = BilinearForm(fes,symmetric=True)
aterms = grad(u)*grad(v) * dx \
    +alpha*order**2/h*jump_u*jump_v * dx(skeleton=True) \
    +(-mean_dudn*jump_v-mean_dvdn*jump_u) * dx(skeleton=True) \
    +alpha*order**2/h*u*v * ds(skeleton=True) \
    +(-n*grad(u)*v-n*grad(v)*u)* ds(skeleton=True)
a += aterms

f = LinearForm(fes)
fterms = rhs*v * dx \
     +alpha*order**2/h*bndc*v * ds(skeleton=True) \
     +(-n*grad(v)*bndc)* ds(skeleton=True)
f += fterms
with TaskManager():
    a.Assemble()
    f.Assemble()
```

--------------------------------

### Generate Circular Obstacle Mesh

Source: https://paulst.github.io/NGSTrefftz/notebooks/qtelliptic.html

Creates a geometry with a circular obstacle and generates a mesh for the circular obstacle example.

```python
h = 0.05
order = 3
wp = WorkPlane().RectangleC(2,2) \
    .Circle(0,0,0.25).Reverse()
geo = wp.Face()
mesh = Mesh(OCCGeometry(geo,dim=2).GenerateMesh(maxh=h))
mesh.Curve(order)
mesh.ngmesh.SetBCName(0,"top")
mesh.ngmesh.SetBCName(1,"bottom")
mesh.ngmesh.SetBCName(2,"left")
mesh.ngmesh.SetBCName(3,"right")
mesh.ngmesh.SetBCName(4,"innercircle")
```

--------------------------------

### Quasi-Trefftz method with tents

Source: https://paulst.github.io/NGSTrefftz/notebooks/twavetents.html

Demonstrates the quasi-Trefftz method for non-constant wavespeed using tent-pitching.

```python
order = 4
delta=500
dt=0.02
ee = exp(-(x)**2*delta )
initc=((ee,-2*(x)*delta * ee, 0, 0))
bdd = CoefficientFunction(( 0,0,0,0))
mesh = Mesh(unit_square.GenerateMesh(maxh = 0.05))
for i in range(0,len(mesh.GetBoundaries())):
    mesh.ngmesh.SetBCName(i,"neumann")

ts = TentSlab(mesh, method="edge")
ts.SetMaxWavespeed(y+1)
ts.PitchTents(dt=dt, local_ct=True)
TT=TWave(order,ts,CoefficientFunction(y+1))
TT.SetInitial(initc)
TT.SetBoundaryCF(bdd)

V = L2(mesh,order=order)
u = GridFunction(V,"u")
gfshow = GridFunction(V, multidim=0)
t=0
SetNumThreads(4)
with TaskManager():
    while t < 0.8:
        TT.Propagate()
        t += dt
        TT.GetWave(u)
        gfshow.AddMultiDimComponent(u.vec)
Draw (gfshow,mesh,interpolate_multidim=True,deformation=True,animate=True,autoscale=False,min=0,max=0.3)
```

```text
WebGLScene
```

--------------------------------

### Get particular solution for right-hand side

Source: https://paulst.github.io/NGSTrefftz/notebooks/qtelliptic.html

Calls fes.GetParticularSolution(rhs) to obtain the particular approximate solution u_f in the quasi-Trefftz space.

```python
with TaskManager():
    uf = fes.GetParticularSolution(rhs)
```

--------------------------------

### Assemble and solve for exact solution

Source: https://paulst.github.io/NGSTrefftz/notebooks/tunfitted.html

Sets up the right-hand side for the exact solution u=exp(x)sin(y), assembles the system, solves, and prints the L2 error and number of elements.

```python
exact = exp(x)*sin(y)
coef_g = exact

f = LinearForm(Vh)
f += (- grad(v) * nh * coef_g + lambda_nitsche / h * coef_g * v) * ds

f.Assemble()
a.Assemble()

nd = Vh.ndof
gfu.vec.data = a.mat.Inverse(Vh.FreeDofs(), 'sparsecholesky') * f.vec

l2error = sqrt(Integrate((gfu - exact)**2 * dx, mesh))
print(l2error)
print(mesh.ne)
```

--------------------------------

### Define manufactured exact solution and coefficients

Source: https://paulst.github.io/NGSTrefftz/notebooks/qtelliptic.html

Defines the exact solution, diffusion tensor K, advection beta, reaction sigma, and right-hand side rhs for the manufactured example.

```python
exact = sin(pi*(x+y))
K = CF((1+x+y,0,0,1+x+y),dims=(2,2))
beta = CF((1,0))
sigma = 3/(1+x+y)
rhs = -sum( (K*CF((exact.Diff(x),exact.Diff(y))))[i].Diff(var) for i,var in enumerate([x,y])) + beta*CF((exact.Diff(x),exact.Diff(y))) + sigma*exact
```

--------------------------------

### Assemble and solve the system

Source: https://paulst.github.io/NGSTrefftz/notebooks/qtwave.html

Assembles the bilinear and linear forms, solves the linear system, and draws the solution.

```python
with TaskManager():
    a.Assemble()
    f.Assemble()

gfu = GridFunction(fes, name="uDG")
gfu.vec.data = a.mat.Inverse()*f.vec
Draw(gfu,clipping={"z":-1, "dist":-.33, "function":True})
```

--------------------------------

### Assemble and solve

Source: https://paulst.github.io/NGSTrefftz/notebooks/twave.html

Assembles the system, solves for the grid function, and draws the result.

```python
with TaskManager():
    a.Assemble()
    f.Assemble()

gfu = GridFunction(fes, name="uDG")
gfu.vec.data = a.mat.Inverse()*f.vec
Draw(gfu,clipping={"z":-1, "dist":-.11, "function":True})
```

```text
WebGLScene
```

--------------------------------

### Solve DG system and compute error

Source: https://paulst.github.io/NGSTrefftz/notebooks/qtelliptic.html

Sets up boundary conditions, calls dgell, solves the linear system, and prints the quasi-Trefftz DG error.

```python
Dbndc = exact
Dbnd = '.*'
Nbndc = 0
Nbnd = ''

a,f = dgell(fes,K,beta,sigma,Dbndc=Dbndc,Dbnd=Dbnd,Nbnd=Nbnd,Nbndc=Nbndc,rhs=rhs,uf=uf)
gfu = GridFunction(fes)

with TaskManager():
    gfu.vec.data = a.mat.Inverse()*f.vec
print("quasi-Trefftz DG error:",sqrt(Integrate((gfu-exact)**2, mesh)))
```

```text
quasi-Trefftz DG error: 0.000440725293956973
```

--------------------------------

### Apply clang-format to C++ files

Source: https://paulst.github.io/NGSTrefftz/contrib.html

Use this command to find all C++ header and source files in the src directory and apply the project's clang-format style.

```bash
find src/ -iname *.hpp -o -iname *.cpp | xargs clang-format -style='file:.clang-format' -i
```

--------------------------------

### Build NGSTrefftz from source

Source: https://paulst.github.io/NGSTrefftz/readme.html

Steps to clone the repository and build the project using cmake and make.

```bash
git clone --recursive https://github.com/PaulSt/NGSTrefftz
mkdir ./NGSTrefftz/make && cd ./NGSTrefftz/make
cmake ../ && make install
```

--------------------------------

### Define boundary conditions

Source: https://paulst.github.io/NGSTrefftz/notebooks/embTrefftz-helm.html

Sets up exact solution and boundary condition coefficients for the Helmholtz problem.

```python
omega=1
n = specialcf.normal(2)
exact = exp(1j*sqrt(0.5)*(x+y))
gradexact = CoefficientFunction((sqrt(0.5)*1j*exact, sqrt(0.5)*1j*exact))
bndc = gradexact*n - 1j*omega*exact
eps = 10**-7
```

--------------------------------

### Solving the acoustic wave equation with tents

Source: https://paulst.github.io/NGSTrefftz/notebooks/twavetents.html

Sets up and propagates the acoustic wave equation using the Trefftz-DG method on the tent-pitched mesh.

```python
for i in range(0,len(mesh.GetBoundaries())):
    mesh.ngmesh.SetBCName(i,"neumann")

order=4
delta = 500
ee = exp(-(x+1)**2*delta )
initc=((ee,-2*(x+1)*delta * ee, 0, 0))
bdd = CoefficientFunction(( 0,0,0,0))

TT=TWave(order,ts,CoefficientFunction(wavespeed))
TT.SetInitial(initc)
TT.SetBoundaryCF(bdd)

V = L2(mesh,order=order)
u = GridFunction(V,"u")
gfshow = GridFunction(V, multidim=0)
t=0
SetNumThreads(4)
with TaskManager():
    while t < 1.5:
        TT.Propagate()
        t += dt
        TT.GetWave(u)
        gfshow.AddMultiDimComponent(u.vec)
Draw (gfshow,mesh,interpolate_multidim=True,deformation=True,animate=True,autoscale=False,min=0,max=0.3)
```

```text
WebGLScene
```

--------------------------------

### TrefftzTents.SetInitial()

Source: https://paulst.github.io/NGSTrefftz/docu.html

Sets the initial condition for the tent-pitching method.

```APIDOC
## TrefftzTents.SetInitial()

### Description
Sets the initial condition for the tent-pitching method.

### Method
SetInitial

### Parameters
No parameters documented.

### Returns
No return value documented.
```

--------------------------------

### Create quasi-Trefftz space with trefftzfespace

Source: https://paulst.github.io/NGSTrefftz/notebooks/qtelliptic.html

Initializes a mesh and a quasi-Trefftz finite element space for the diffusion-advection-reaction equation using the 'qtelliptic' equation keyword.

```python
h = 0.03;
order = 3;
mesh = Mesh(unit_square.GenerateMesh(maxh=h))
fes = trefftzfespace(mesh,order=order,eq="qtelliptic")
```

--------------------------------

### Import ngsolve, ngstrefftz, and netgen.occ

Source: https://paulst.github.io/NGSTrefftz/notebooks/laplace.html

Import required modules for building the mesh, Trefftz space, and finite element formulation.

```python
from ngsolve import *
from ngstrefftz import *
from netgen.occ import *

```

--------------------------------

### TrefftzTents.Propagate()

Source: https://paulst.github.io/NGSTrefftz/docu.html

Propagates the solution in the tent-pitching method.

```APIDOC
## TrefftzTents.Propagate()

### Description
Propagates the solution in the tent-pitching method.

### Method
Propagate

### Parameters
No parameters documented.

### Returns
No return value documented.
```

--------------------------------

### Define test case

Source: https://paulst.github.io/NGSTrefftz/notebooks/twave.html

Sets the exact solution and Dirichlet boundary conditions for the wave equation.

```python
from numpy import pi
sq = sqrt(2.0)
sig_ex = CoefficientFunction((-cos(pi*x)*sin(pi*y)*sin(pi*z*c*sq)/sq,
    -sin(pi*x)*cos(pi*y)*sin(pi*z*c*sq)/sq))
v_ex = sin(pi*x)*sin(pi*y)*cos(pi*z*c*sq)*c
gD = v_ex
```

--------------------------------

### Set up Trefftz FE space with restriction

Source: https://paulst.github.io/NGSTrefftz/notebooks/tunfitted.html

Creates a Trefftz space for the Laplace equation with DG jumps, restricts it to elements with negative level set, and defines trial/test functions and a grid function.

```python
ci = CutInfo(mesh, lsetp1)
els_hasneg = ci.GetElementsOfType(HASNEG)
els_if = ci.GetElementsOfType(IF)
# Facets with parts inside the domain
facets_dg = GetFacetsWithNeighborTypes(mesh, a=els_hasneg, b=els_hasneg)

Vhbase = trefftzfespace(mesh, order=order, eq='laplace', dgjumps=True)
Vh = Restrict(Vhbase, els_hasneg)
u, v = Vh.TnT()
gfu = GridFunction(Vh)
```

--------------------------------

### Pitching tents on the mesh

Source: https://paulst.github.io/NGSTrefftz/notebooks/twavetents.html

Configures the tent slab and pitches tents based on the specified wavespeed and time step.

```python
dt = 0.05   # tent slab height (in time)
wavespeed = 1
ts = TentSlab(mesh, method="edge")
ts.SetMaxWavespeed(wavespeed)
ts.PitchTents(dt=dt, local_ct=True)
print("max slope", ts.MaxSlope())
print("n tents", ts.GetNTents())
```

```text
max slope 1.0
n tents 705
```

--------------------------------

### Import libraries and set parameters

Source: https://paulst.github.io/NGSTrefftz/notebooks/embTrefftz-stokes.html

Imports the required libraries and sets the number of threads to 4. Defines the viscosity nu as 1.

```python
from ngsolve import *
from ngstrefftz import *
from netgen.occ import *
from ngsolve.webgui import Draw
SetNumThreads(4)
nu=1
```

--------------------------------

### Create a feature branch

Source: https://paulst.github.io/NGSTrefftz/contrib.html

Use this command to create and switch to a new branch for your feature implementation.

```bash
git checkout -b <your-feature-name>
```

--------------------------------

### Visualize results

Source: https://paulst.github.io/NGSTrefftz/notebooks/embTrefftz-ns.html

Commands to visualize the velocity magnitude and pressure field of the computed solution.

```python
Draw(Norm(uh), mesh, "speed");
Draw(ph - Integrate(ph, mesh) / Integrate(1, mesh), mesh, "pressure");
```

--------------------------------

### Assemble bilinear and linear forms

Source: https://paulst.github.io/NGSTrefftz/notebooks/qtwave.html

Defines the bilinear form for the quasi-Trefftz DG method and the linear form with boundary conditions.

```python
a = BilinearForm(fes)
dttau = CoefficientFunction(tuple([grad(wtau)[i,D] for i in range(D)]))
gtau = sum([grad(wtau)[i*(D+2)] for i in range(D)])
dtw =grad(wtau)[(D+1)*(D+1)-1]
gw=CoefficientFunction(tuple([grad(wtau)[D,i] for i in range(D)]))
a += (- v*(gtau + pow(c,-2)*dtw) ) * dx
a += ( - sig*(dttau + gw ) ) * dx
a += spacelike * ( pow(c,-2)*IfPos(n_t,v,vo)*jump_wt + IfPos(n_t,sig,sigo)*jump_taut) * dx(skeleton=True)
a += timelike * ( mean_v*jump_taux + mean_sig*jump_wx + alpha*jump_vx*jump_wx + beta*jump_sigx*jump_taux ) * dx(skeleton=True )
a += ( pow(c,-2)*v*w + sig*tau) * IfPos(n_t,1,0)*spacelike * ds(skeleton=True)
a += ( sig*n_x*w + alpha*v*w ) * timelike * ds(skeleton=True)

f = LinearForm(fes)
f += ( pow(c,-2)*v_ex*w + sig_ex*tau ) * IfPos(n_t,0,1)*spacelike * ds(skeleton=True)
f += ( gD * (alpha*w - tau*n_x) ) * timelike *  ds(skeleton=True)
```

--------------------------------

### Oseen Trefftz space construction with NGSolve

Source: https://paulst.github.io/NGSTrefftz/notebooks/embTrefftz-ns.html

Sets up the polynomial trial and test spaces, defines the Oseen operator with a fixed convection field, and creates the embedded Trefftz space. Prints the degrees of freedom for the full polynomial space and the embedded Trefftz space.

```python
mesh = Mesh(unit_square.GenerateMesh(maxh=0.45))
k = 3
nu = 1.0
wn = CF((1, 0))

V = VectorL2(mesh, order=k, dgjumps=True)
Q = L2(mesh, order=k - 1, dgjumps=True)
base = V * Q

R = VectorL2(mesh, order=k - 2) * L2(mesh, order=k - 1)
u, p = base.TrialFunction()[0:2]
v, q = R.TestFunction()[0:2]

op = (-nu * VLap(u, mesh.dim) + grad(p) + grad(u) * wn) * v * dx(bonus_intorder=16)
op += div(u) * q * dx(bonus_intorder=16)

emb = TrefftzEmbedding(op)
tfes = EmbeddedTrefftzFES(emb)

print(f"full Oseen polynomial dofs: {base.ndof}")
print(f"embedded Oseen Trefftz dofs: {tfes.ndof}")

```

```text
full Oseen polynomial dofs: 156
embedded Oseen Trefftz dofs: 84

```

--------------------------------

### Create mesh and Trefftz finite element spaces

Source: https://paulst.github.io/NGSTrefftz/notebooks/helmholtz.html

Creates a unit square mesh and two Trefftz finite element spaces for the Helmholtz equation and its conjugate, both of order 3 with complex and DG jumps enabled.

```python
mesh = Mesh(unit_square.GenerateMesh(maxh=.3))
fes = trefftzfespace(mesh,order=3,eq="helmholtz",complex=True,dgjumps=True)
fes2 = trefftzfespace(mesh,order=3,eq="helmholtzconj",complex=True,dgjumps=True)
```

--------------------------------

### TreffzTents.Propagate

Source: https://paulst.github.io/NGSTrefftz/docu.html

Solves the tent slab for the current time step.

```APIDOC
## TreffzTents.Propagate()

### Description
Solve tent slab.

### Method
Propagate
```

--------------------------------

### Visualize solution

Source: https://paulst.github.io/NGSTrefftz/notebooks/embTrefftz-poi.html

Draws the computed solution using the web GUI.

```python
from ngsolve.webgui import Draw
Draw(BoundaryFromVolumeCF(gfu),mesh,'gfu')
```

--------------------------------

### Assemble and solve the Trefftz-DG system

Source: https://paulst.github.io/NGSTrefftz/notebooks/helmholtz.html

Defines the DG bilinear form and linear form for the Helmholtz equation, assembles them, solves for the finite element solution, and computes the L2 error against the exact solution.

```python
h = specialcf.mesh_size
alpha = 1/(omega*h)
beta = omega*h
delta = omega*h

u = fes.TrialFunction()
v = fes.TestFunction()
if fes2 is not None:
    v = fes2.TestFunction()
jump = lambda u: (u-u.Other())*n
mean = lambda u: 0.5 * ((u)+(u.Other()))

a = BilinearForm(fes,fes2)
a += grad(u)*(grad(v))*dx - omega**2*u*(v)*dx

a += -(jump(u)*(mean(grad(v)))+mean(grad(u))*jump(v)) * dx(skeleton=True)
a += -1/(omega*1j)*beta*jump(grad(u))*(jump(grad(v))) * dx(skeleton=True)
a += omega*1j*alpha*jump(u)*jump(v) * dx(skeleton=True)

a += -delta*(u*(grad(v))*n+grad(u)*n*(v)) * ds(skeleton=True)
a += -1/(omega*1j)*delta*(grad(u)*n)*((grad(v))*n) * ds(skeleton=True)
a += omega*1j*(1-delta)*u*(v) * ds(skeleton=True)

f = LinearForm(fes2)
f += -1/(omega*1j)*delta*bndc*(grad(v))*n*ds(skeleton=True)
f += (1-delta)*bndc*(v)*ds(skeleton=True)

with TaskManager():
    a.Assemble()
    f.Assemble()

gfu = GridFunction(fes)
with TaskManager():
    gfu.vec.data = a.mat.Inverse() * f.vec

error = sqrt(Integrate((gfu-exact)*Conj(gfu-exact), mesh).real)
print("error ",error)
```

--------------------------------

### Initialize NGSTrefftz environment

Source: https://paulst.github.io/NGSTrefftz/notebooks/embTrefftz-poi.html

Imports necessary modules and sets the number of threads for parallel computation.

```python
from ngsolve import *
from ngstrefftz import *
from netgen.occ import *
SetNumThreads(4)
```

--------------------------------

### TrefftzTents.SetInitial

Source: https://paulst.github.io/NGSTrefftz/docu.html

Sets the initial condition for the wave equation.

```APIDOC
## TrefftzTents.SetInitial(CoefficientFunction)

### Description
Set initial condition.

### Method
SetInitial

### Parameters
- **arg0** (ngsolve.fem.CoefficientFunction) – Coefficient function describing the initial condition.
```

--------------------------------

### Define the exact solution and boundary condition

Source: https://paulst.github.io/NGSTrefftz/notebooks/laplace.html

Sets the exact solution exp(x)*sin(y) as the boundary condition for the problem.

```python
exact = exp(x)*sin(y)
bndc = exact

```

--------------------------------

### Define finite element space and Trefftz-DG formulation coefficients

Source: https://paulst.github.io/NGSTrefftz/notebooks/embTrefftz-wave.html

Sets up the full polynomial finite element space on the mesh. Defines the trial/test functions, the jump and mean operators for the dimensionless framework, and the mesh size, normal, and parameters for the Trefftz-DG bilinear and linear forms. The assembly is performed in parallel with a task manager.

```python
order=3
X = L2(mesh,order=order,dgjumps=True)
fes=FESpace([X for _ in range(mesh.dim)])

```

```python
D = fes.mesh.dim - 1

vsig = fes.TrialFunction()
wtau = fes.TestFunction()
vsigo = CoefficientFunction(tuple([vsig[i].Other() for i in  range(D+1)]))
wtauo = CoefficientFunction(tuple([wtau[i].Other() for i in  range(D+1)]))

v = vsig[D]
sig = CoefficientFunction(tuple([vsig[i] for i in  range(D)]))
w = wtau[D]
tau = CoefficientFunction(tuple([wtau[i] for i in  range(D)]))

vo = vsigo[D]
sigo = CoefficientFunction(tuple([vsigo[i] for i in  range(D)]))
wo = wtauo[D]
tauo = CoefficientFunction(tuple([wtauo[i] for i in  range(D)]))

h = specialcf.mesh_size
n = specialcf.normal(D+1)
n_t = n[D]/Norm(n)
n_x = CoefficientFunction( tuple([n[i]/Norm(n) for i in  range(D)]) )

mean_v = 0.5*(v+vo)
mean_w = 0.5*(w+wo)
mean_sig = 0.5*(sig+sigo)
mean_tau = 0.5*(tau+tauo)

jump_vx = ( v - vo ) * n_x
jump_wx = ( w - wo ) * n_x
jump_sigx = (( sig - sigo ) * n_x)
jump_taux = (( tau - tauo ) * n_x)

jump_vt = ( v - vo ) * n_t
jump_wt = ( w - wo ) * n_t
jump_sigt = ( sig - sigo ) * n_t
jump_taut = ( tau - tauo ) * n_t

timelike = n_x*n_x # n_t=0
spacelike = n_t**2 # n_x=0

#params
alpha=0.5
beta=0.5

a = BilinearForm(fes)
#space like faces, w/o x jump ASSUME TENSOR MESH
a += spacelike * ( pow(c,-2)*IfPos(n_t,v,vo)*jump_wt + IfPos(n_t,sig,sigo)*jump_taut) * dx(skeleton=True)
#time like faces
a += timelike * ( mean_v*jump_taux + mean_sig*jump_wx + alpha*jump_vx*jump_wx + beta*jump_sigx*jump_taux ) * dx(skeleton=True )        #t=T (or *x)
a += ( pow(c,-2)*v*w + sig*tau) *  IfPos(n_t,1,0)*spacelike * ds(skeleton=True)
#dirichlet boundary 'timelike'
a += ( sig*n_x*w + alpha*v*w ) * timelike *   ds(skeleton=True)

f = LinearForm(fes)
f += ( pow(c,-2)*v_ex*w + sig_ex*tau ) * IfPos(n_t,0,1)*spacelike *  ds(skeleton=True)
f += ( gD * (alpha*w - tau*n_x) ) * timelike *   ds(skeleton=True)

with TaskManager():
    a.Assemble()
    f.Assemble()

```

--------------------------------

### Define exact solution, mesh, and DG space

Source: https://paulst.github.io/NGSTrefftz/notebooks/embTrefftz.html

Sets up the exact solution, polynomial order, mesh, and L2 space with DG jumps for the Laplace problem.

```python
exact = exp(x+y)*sin(sqrt(2)*z)
order = 4
mesh = Mesh(unit_cube.GenerateMesh(maxh=.2))
fes = L2(mesh, order=order,  dgjumps=True)
```

--------------------------------

### Construct quasi-Trefftz space

Source: https://paulst.github.io/NGSTrefftz/notebooks/qtwave.html

Creates a quasi-Trefftz finite element space for the first-order wave equation using the 'foqtwave' equation and sets the coefficient.

```python
fes = trefftzfespace(mesh,order=3,dgjumps=True,eq="foqtwave")
fes.SetCoeff(c)
```

--------------------------------

### Import required modules and set order

Source: https://paulst.github.io/NGSTrefftz/notebooks/tunfitted.html

Imports ngsolve, xfem, ngstrefftz, and netgen.occ, and sets the polynomial order to 4.

```python
from ngsolve import *
from ngsolve import dx as ngsdx
from xfem import *
from xfem.lsetcurv import *
from ngstrefftz import *
from netgen.occ import *
order=4
```

--------------------------------

### Initialize Mesh and Solve with Trefftz Embedding

Source: https://paulst.github.io/NGSTrefftz/notebooks/embTrefftz-stokes.html

Configures the finite element spaces, initializes the mesh, and solves the Stokes problem using the Trefftz embedding method.

```python
order = 3
mesh = Mesh(unit_square.GenerateMesh(maxh=.1))

V = VectorL2(mesh, order=order, dgjumps=True)
Q = L2(mesh, order=order - 1, dgjumps=True)
fes = V * Q
a,f = StokesDG(fes, 1, rhs)

PP, uf = TrefftzStokesEmb(fes, rhs)
PPT = PP.CreateTranspose()
TA = PPT @ a.mat @ PP
Tgfu = TA.Inverse() * (PPT * (f.vec - a.mat * uf))
gfu = GridFunction(fes)
gfu.vec.data = PP * Tgfu + uf
```

--------------------------------

### Assemble Bilinear and Linear forms

Source: https://paulst.github.io/NGSTrefftz/notebooks/embTrefftz-adv.html

Configuration of the Upwind DG formulation for the advection equation.

```python
bndc=exact
n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size

jump_u = (u-u.Other())*n
jump_v = (v-v.Other())*n

a = BilinearForm(fes,symmetric=False)
a += - u * b * grad(v) * dx
a += b*n*IfPos(b*n,u,u.Other())*v*dx(element_boundary=True)

f = LinearForm(fes)
f += rhs*v * dx \
     - b*n*IfPos(b*n,0,bndc)*v * ds(skeleton=True)
with TaskManager():
    a.Assemble()
    f.Assemble()
```
Tokens:
Raw