How to calculate the error on the skeleton of mesh

Main Forum Old-Forum - do not post here
1

#mesh
from netgen.geom2d import SplineGeometry
geo = SplineGeometry()
geo.AddRectangle( (0, 0), (1, 1))
mesh = Mesh( geo.GenerateMesh(maxh=1/16))

velocity = CoefficientFunction(UN.components[0:2])
err_u = sqrt (Integrate ( (velocity-exact_u)*(velocity-exact_u), mesh))

Now，I want to calculate the error on the skeleton of mesh， namely ，the L^2 error of mass flux jump [v_h] \cdot n
#tex
[tex] \sum_{F \in \mathcal{F}{h}} \frac{1}{h{F}}\left| \llbracket {v}{h} \rrbracket \cdot {n}{F}|{L^{2}(F)}^{2}[/tex]
#\mathcal{F}{h} stand for all facets of mesh

2

Hi Yongbin,

The simplest way to compute this quantity is by defining a bilinear form a so that a(u,u) coincides with your (sum of) integral(s). Then you Apply the bilinear form to your discrete solution u and take the inner product of the result with u again.

Best,
Christoph

3

Hi，Christoph，

Thank you very much for your reply. Sorry，I may not fully understand what you mean. This is my code for Stokes . I don’t know how to calculate in my code. In addition，I also wonder， is there other way to calculate？

Thank you again！

Yongbin

4

Hi Yongbin,

Here is the code fitting your script. I use hybrid dg jumps here as you use a hybrid DG formulation:

a_jumps = BilinearForm(X)
a_jumps += (ubar-uu)*(vbar-vv) * dx(element_boundary = True)
AUNvec = UN.vec.CreateVector()
a_jumps.Apply(UN.vec,AUNvec)
jump_norm = sqrt(InnerProduct(AUNvec,UN.vec))
print(jump_norm)

Alternatively, if you are interested in the standard DG jumps the code looks as follows:

uuOther = CoefficientFunction((ux.Other(), uy.Other()))
vvOther = CoefficientFunction((vx.Other(), vy.Other()))
a_jumps = BilinearForm(X)
a_jumps += (uuOther-uu)*(vvOther-vv) * dx(skeleton = True)
AUNvec = UN.vec.CreateVector()
a_jumps.Apply(UN.vec,AUNvec)
jump_norm = sqrt(InnerProduct(AUNvec,UN.vec))
print(jump_norm)

Alternatives? I think this is the most straight-forward way to do it. As of now there is no “Integrate” equivalent that would allow you to integrate the skeleton-terms directly. However, I think the above approach is also very natural.

Best,
Christoph

5

I have successfully solved this problem according to your reply, thank you very much for your reply!

Yongbin Han

6

Hi Yongbin,

often one needs these terms in residual-type error estimators. Then you need them element by element.

This is now also supported:

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

    elerr = Integrate ( h*(n*(lam*grad(gfu)-(lam*grad(gfu)).Other()))**2*dx(element_boundary=True), mesh, element_wise=True)

full example is attached.

You can integrate elementwise, on every element boundary, and use also values from the neighbouring elements.

This is a quite new feature and needs some more testing. Faces at the domain boundary are skipped, at the moment.

Joachim

Attachment: residual_ee.py

7

Dear Professor !

Thank you very much for your advice and guidance. Let me have a further understanding of this problem.

Yongbin Han