co-normal value on a 1D surface segment

Main Forum Old-Forum - do not post here
1

Hello, I am trying to hack a co-normal value for 1D surface segment that is needed for a HDG formulation…
Basically, for each line segment on a 1D surface mesh (embedded in a 2D mesh), I want to assign value +1 or -1 to the two nodes of the segment based on its tangential direction… any suggestion to hack this?
In the 2D surface case, we can do Cross(n, t), but it is not available in 1D.

2

Well, I found a not so elegant solution based on orientation of the edge and a SurfaceL2 grid function…

gfn = GridFunction(Compress(SurfaceL2(mesh, order=1, definedon=mesh.Boundaries("interface"))))
    count = 0
    for e in mesh.Elements(BND):
        if e.mat=="interface":
            if e.vertices[0].nr > e.vertices[1].nr:
                gfn.vec[2*count+1] = 1
            else:
                gfn.vec[2*count+1] = -1 
            count += 1
3

Well, I found a not so elegant solution based on orientation of the edge and a SurfaceL2 grid function…

    gfn = GridFunction(Compress(SurfaceL2(mesh, order=1, definedon=mesh.Boundaries("interface"))))
    count = 0
    for e in mesh.Elements(BND):
        if e.mat=="interface":
            if e.vertices[0].nr > e.vertices[1].nr:
                gfn.vec[2*count+1] = 1
            else:
                gfn.vec[2*count+1] = -1 
            count += 1
4

Hi Guosheng,

for a 1D-surface in 2D you can use the specialcf.tangential function for the co-normal vector.

conormal = specialcf.tangential(2)

Best,
Michael

5

Hi Michael,

Well, I want a hybrid mixed method on a 1D surface mesh, for this purpose, I need to have a co-normal on the nodes that changes signs when evaluated from neighboring segments.
Since there is no SurfaceHDiv space for 1D, currently, I am simply using SurfaceL2 space for the (scalar) velocity field, and evaluate the surface divergence as

grad(u).Trace()*t

where t is the tangential vector.
The term that is missing is the transmission condition that enforce back continuity of velocity

qhat*u*co_n*ds(element_boundary=True"interface")

where qhat is a H1-P1 (hybrid) space on the surface and u is the SurfaceL2 velocity space. You can see that, basically, I need the co-normal direction co_n that is like the normal constant in 1D mesh…
The previous snipt is a way to hack such a co-normal coefficient using P1-SurfaceL2, but I am not confident it will work for all meshes.
Any ideas?

Best,
Guosheng

6

Hi Guosheng,

it is possible to imitate the HDivSurface by using an H1-space and equip it with the tangential vector. E.g. a lowest order HDivSurface could be achieved by

fesH = H1(mesh, order=1)
u,v = fesH.TnT()

con = specialcf.tangential(2)
u = u*con
v = v*con

Then, u and v are co-normal continuous functions on the mesh.

Another idea, similar to your code with the edge orientation, is to use two different specialcf.tangential vectors

nel   = specialcf.tangential()
nel_c = specialcf.tangential(consistent=True)

The second one does not change its sign depending from which element you are coming. Thus the product nel*nel_c will be 1 from the one and -1 from the other side. The consistent=True flag is rather new, but it should work also for 1D-surfaces.

Best,
Michael

7

Hi Michael,

Please see the attached minimal example where I use the 1D normal. I need it for the evaluation of the element boundary of surface segments (codim 2). For this purpose, the co-normal shall (consistently) take -1 on the left end point and +1 on the right end point, which was what the gfn trying to do in the code.
I can not make it work using the consistent tangential vector

from ngsolve import *
from ngsolve.meshes import MakeStructured2DMesh
N = 4
mesh = MakeStructured2DMesh(nx=N, ny=1, mapping=lambda x, y: (x, y/N))

gfn = GridFunction(Compress(SurfaceL2(mesh, order=1, definedon=mesh.Boundaries("bottom"))))
count = 0
for e in mesh.Elements(BND):
    if e.mat=="bottom":
        if e.vertices[0].nr > e.vertices[1].nr:
            gfn.vec[2*count+1] = -1
        else:
            gfn.vec[2*count+1] = 1 
        count += 1    

n = specialcf.normal(2)
# HDiv interpolation
gfInterp = GridFunction(HDiv(mesh, order=1))
gfInterp.Set(gfn*n, definedon=mesh.Boundaries("bottom"))
Draw(gfInterp[1], mesh)

## 
V = Compress(H1(mesh, definedon=mesh.Boundaries("bottom")))
u,v = V.TnT()
t = specialcf.tangential(2)
t1 = specialcf.tangential(2, consistent=True)

f = LinearForm(V)
# This works
f += t*t1*v*ds(element_boundary=True, definedon=mesh.Boundaries("bottom"))
# f += gfn*v*ds(element_boundary=True, definedon=mesh.Boundaries("bottom"))
f.Assemble()
print(f.vec)

Best,
Guosheng

8

Hi Guosheng,

yes, you are right, there is a detail missing for the consistent tangential vector for 1D manifolds.

Until this is fixed you can rotate the normal vector, which does not change from element to element.

n = specialcf.normal(2)
t1 = CF( (-n[1],n[0]) )

With this I got the correct results for your minimal example.

Best,
Michael

9

Hi Guosheng,

we can now access Jacobians, and the coordinates from the reference element via special CoefficientFunctions (thx to Michael). This allows to build the co-normal vector.

The Logging - CoefficientFunction is useful to inspect the function evaluations.

  • Joachim
from netgen.geom2d import unit_square
from ngsolve import *
mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))


fes = SurfaceL2(mesh, order=1)
u,v = fes.TnT()
f = LinearForm(fes)

t = specialcf.tangential(2)
jac = specialcf.JacobianMatrix(2,1)[:,0]
xref = specialcf.xref(1)
conormal = (2*xref[0]-1)*jac

from ngsolve.fem import LoggingCF
conormal = LoggingCF (conormal)

f += conormal*t  *v * ds("bottom", element_boundary=True)

f.Assemble()
10

Thanks!
I’ve now moved to a 3D mesh with a embedded 2D surface. I use mixed FEM to discrete a Laplace operator on the surface.
The following test case does not work: mesh is a unit cube, and the “flat surface” is simply its top boundary.
I use HDivSurface/SurfaceL2 spaces to build the finite element pair on “top” boundary, but it seems that I can not use these spaces directly here as the bilinear form assembler complained about inconsistent # of DOFs. Any idea what’s going on?!

from ngsolve.meshes import MakeStructured3DMesh
from ngsolve import *
mesh = MakeStructured3DMesh(hexes=False, nx=10, ny=10, nz=2)

V = Compress(HDivSurface(mesh, order=1, definedon=mesh.Boundaries("top")))
W = Compress(SurfaceL2(mesh, order=0, definedon=mesh.Boundaries("top")))
fes = V*W
(u, p), (v, q) = fes.TnT()
a = BilinearForm(fes)
a += (u.Trace()*v.Trace()-div(u).Trace()*q.Trace()-p.Trace()*div(v).Trace())*ds("top")
a.Assemble()
11

Hi Guosheng,

a definedon-check was missing in the HDivSurface, should be fixed in the next nightly-build.

Best,
Michael

12

Hi Michael,

Thanks for the very quick fix. I will check it out.
Also, I just realized that your co-normal is not exactly what I am looking for… I haven’t express the problem clear enough. In the attached file, I depicted what I want: a scalar function a 1D surface mesh that acts exactly as the 1D normal direction for a plain 1D mesh, that is, given an orientation of the line segment, it should return a -1 on the left end of the segment and +1 on the right end of segment. (I don’t think it should be called a co-normal as it is a lower dimensional quantity)… Since I couldn’t function such a built-in function, the above surfaceL2 gridfunction gfn is a hacker that do the job.

Since the tangential/normal directions are quantities defined on edges, not vertices, the associated value on left/right endpoints of the segments will always be the same, which is not what I want.

*** More background: I use this “co-normal” to reinforce nodal continuity of a SurfaceL2 function:

# u is a surface L2 trial function
# v is a H1 trial function 

# this bilinearform shall return nodal continuity of u (hybridized!)
a += u*v*con*ds(element_boundary=True)

This has been a very fruitful discussion! Thanks again for the quick responses.

Best,
Guosheng

Attachment: note.pdf

13

The surface mixed FEM now works as expected. I want to hybridize the surface mixed FEM, but the discontinuous flag seems not working properly for surface HDiv: I am expecting with the discontinuous flag, the space will be co-normal discontinuous across surface edges, but this is not the case in the attached code

from ngsolve.meshes import MakeStructured3DMesh
from ngsolve import *
mesh = MakeStructured3DMesh(nx=1, ny=1, nz=1)
VC = Compress(HDivSurface(mesh, order=0, RT=True, discontinuous=False, definedon=mesh.Boundaries("top")))
print(VC.ndof, " free: ", VC.FreeDofs())
VD = Compress(HDivSurface(mesh, order=0, RT=True, discontinuous=True, definedon=mesh.Boundaries("top")))
print(VD.ndof, " free: ", VD.FreeDofs())