### 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 ``` -------------------------------- ### 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 ``` -------------------------------- ### 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 ``` -------------------------------- ### 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 ``` -------------------------------- ### 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() ```