Plotting DG function in 3D

[Sander Rhebergen]

Hello,

When trying to plot a zero function in 3D from the L2 space I get a segmentation error. However, when I change my space to H1, I can plot just fine. Also, I can plot just fine in 2D (both L2 and H1). I have attached the code below. Any help is very much appreciated.

thanks, Sander

from ngsolve import *

#2D
#from netgen.geom2d import unit_square
#mesh = Mesh(unit_square.GenerateMesh(maxh=1.6))

# 3D
from netgen.csg import *
cube = OrthoBrick( Pnt(0,0,0), Pnt(1,1,1))
geo = CSGeometry()
geo.Add (cube)
mesh = Mesh( geo.GenerateMesh(maxh=1.6))

# Function spaces
V = H1(mesh, order=1)
#V = L2(mesh, order=1)

# 2D
#X    = FESpace([V, V])
#gfu  = GridFunction(X)
#velocity = CoefficientFunction (gfu.components[0:2])

# 3D
X = FESpace([V, V, V])
gfu  = GridFunction(X)
velocity = CoefficientFunction (gfu.components[0:3])

# for visualization
Draw (Norm(velocity), mesh, "velocity")

[Anonymous]

[Anonymous]

[Anonymous]

[schruste]

Dear Sander,

The segmentation fault happens when Netgen tries to draw surface values. It tries to evaluate traces of the L2 space, but as the trace of an L2 space is not well defined in NGSolve the evaluation crashes (I assume this is due to didactical reasons mainly).
You can however use the option 'draw_surface=False' when calling the Draw-method this will prevent Netgen from trying to evaluate L2 traces.

Best,
Christoph

[Sander Rhebergen]

Dear Christoph,

Great - that works. However, now I have a related question. For DG you need numerical fluxes which depend of the values of your function evaluated on the element boundary. How would you in NGSolve then compute this "trace"?

thanks,
Sander

[schruste]

Dear Sander,

The issue that you had is related to the way a function is evaluated, in the previous case the "trace evaluator" is not well defined.

DG+L2FESpaces still work well together because the way you evaluate L2-functions in an NGSolve-DG-formulation is not using the trace evaluation as the surface visualization does:
When defining bilinear forms for assembly or application you will only use the trace evaluator when using "VOL_or_BND=BND" and that would be a bad idea again for L2 spaces. However, there are other "modes" for defining integrals for bilinear forms. If you want to do a DG-type integral over the skeleton you can use skeleton=True (in combination with "VOL_or_BND=BND" it will give you only the boundary facets, in combination with "VOL_or_BND=VOL" it will give you only interior facets (default)), or if you want loops over domain boundary (e.g. for HDG) you can use element_boundary=True. When the corresponding terms are collect during assembly NGSolve will evaluate the shape functions at integration points w.r.t. the integration point within a volume element. This type of evaluation is well defined for most shape functions in NGSolve, especially for L2 functions.

Perhaps an illustrating example:
You can do

SymbolicBFI( u*v, VOL_or_BND=BND)

and

SymbolicBFI( u*v, VOL_or_BND=BND, skeleton=True)

The first version will use an integration point on the boundary facet and ask for evaluation on the surface element (trace evaluation) while the second one will evaluate via volume integration points (will be only on one facet of the volume element however). The (mapped!) coordinates of the integration points will be the same, however, the shape function evaluation will be different as the reference elements are different ones.

In NGSolve for the first version the FESpace is asked for GetSFE (get surface finite element), in the second version the FESpace is only asked for GetFE (get finite element (of the corresponding neighboring volume element)). The first version will not work with L2Space, they do not provide the GetSFE (the trace reasoning), the second version is however fine (essentially only a volume evaluation).

Hope this helped,
Christoph

[Sander Rhebergen]

Dear Christoph,

Thanks again for the explanation. I'll use your suggestions to continue my 3D HDG implementation.

thanks,
Sander

[schruste]

Dear Sander,

Perhaps I confused you more than I intended? I am not sure what you interpreted from my last post. Too make this clear: In NGSolve you can do "standard" DG equally as well as HDG.

Here a simple snippet for an interiori penalty DG formulation:

V = L2(mesh, order=p, flags={ 'dgjumps': True })

u = V.TrialFunction()
v = V.TestFunction()

n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size

a = BilinearForm (V, symmetric=True)
a += SymbolicBFI(grad(u) * grad(v))
a += SymbolicBFI(0.5 * (grad(u) + grad(u.Other())) * n * (v - v.Other()), VOL, skeleton=True)
a += SymbolicBFI(0.5 * (grad(v) + grad(v.Other())) * n * (u - u.Other()), VOL, skeleton=True)
a += SymbolicBFI(10 * p * (p+1) / h * (u - u.Other()) * (v - v.Other()), VOL, skeleton=True)
a += SymbolicBFI(10 * p * (p+1) / h * (u) * (v), BND, skeleton=True)
a.Assemble()

Best,
Christoph

[Sander Rhebergen]

Dear Christoph,

Thanks for the snippet. I'm actually trying to get an HDG code running. Everything is working fine in 2D, but just getting things to work in 3D seems a bit more tricky.

For my velocity spaces I am using standard DG spaces:

V = L2(mesh, order=k)

and for the HDG velocity facet spaces I have for example

VF = FacetFESpace(mesh, order=k, dirichlet=["wall|inlet"])

I can create then a mixed space, for example

X = FESpace([V, V, V, VF, VF, VF])

(I'm leaving out the pressure spaces for simplicity here).

I'm then using

a = BilinearForm(X, flags={"eliminate_internal" : True})
a += SymbolicBFI(element_terms)

to compute element integral terms and

a += SymbolicBFI(element_boundary_terms, element_boundary=True)

to compute element boundary terms. As far as I understood from your previous emails, this should work for HDG, correct?

To set a Dirichlet boundary condition on the velocity trace space, I do for example:

UN = GridFunction(X)
uin = 4.0*x
UN.components[3].Set(uin, definedon=mesh.Boundaries("inlet"))

Most of this is just straightforward extension from what I have working in 2D. Since this is an HDG code, do I still need to include the flags={ 'dgjumps': True } in the definition of the V function space? This was not needed in 2D. It doesn't seem to make much difference if I include it in my 3D code or not (my 3D code isn't working right now so this may or may not be going wrong).

thanks,
Sander

[Anonymous]

Hi Sanders,

wow, you already got well into ngs-py.

you should NOT include flags={ 'dgjumps': True } for the HDG method. This flags reservers space for potential non-zero sparse-matrix entries according to a DG method, what is expensive. If you don't use it, you reserve the correct sparse-matrix pattern for the HDG method.

There should not be any larger difference between the 2D and 3D ngs-py scripts.
The difference comes in when we worry about solvers.

Best,
Joachim

[Sander Rhebergen]

Hello,
Attached I have two codes (2D and 3D) solving the Poisson equation using an interior penalty HDG method. The 2D code works perfectly, but the 3D code unfortunately doesn't. I'm wondering if I'm missing something that I need to set in NGSolve that is different for 3D than it is in 2D. Attached also the geo file with which I make a cube in 3D.
thanks,
Sander

[Joachim Schoeberl]

Hi Sander,
in the dirichlet flag you can use either a list of intergers, or one string treated as regex (e.g. "left|right"). You had a list of a string (which was silently ignored)

FacetFESpace(mesh, order=k, dirichlet="wall")

You can check the Dirichlet dofs via

print (X.FreeDofs())

which prints an array of 0/1, where the 1 stand for unconstrained dofs.

a few notes to the file:
1. you get the inner product of vector-valued functions with grad(u)*grad(v), or alternatively as InnerProduct(grad(u), grad(v))
2. you can model the geometry also within Python. it is convenient, since you have to send only one file.

Joachim

[Sander Rhebergen]

Hi Joachim,

Thank you very much - I was indeed having trouble with the bcs, but that clears it up!

Sander