Zero matrix

NGSolve
1

Hello,

Upon running the following code, using dsf1 or dsf2, I get a matrix of zeroes. However, I am anticipating a non-zero matrix. What can be causing the error?

from ngsolve import *
from netgen.geom2d import SplineGeometry
ngsglobals.msg_level=1

def topBottom():
    geometry = SplineGeometry()
    pnts     = [ (0, 0),  (0.5, 0), (1, 0), (1, 1), (0.5, 1), (0, 1)  ]
    pnums    = [geometry.AppendPoint(*p) for p in pnts]
    lines = [ (pnums[0], pnums[1], "gamma1",  1, 0),
              (pnums[1], pnums[2], "gamma2",  2, 0),
              (pnums[2], pnums[3], "gamma2",  2, 0),
              (pnums[3], pnums[4], "gamma2",  2, 0),
              (pnums[4], pnums[5], "gamma1",  1, 0),
              (pnums[5], pnums[0], "gamma1",  1, 0),
              (pnums[1], pnums[4], "gammaf", 1, 2)]
    for p1,p2,bc,left, right in lines:
        geometry.Append( ["line", p1, p2], bc=bc, leftdomain=left, rightdomain=right)
    geometry.SetMaterial(1, "Domain1")
    geometry.SetMaterial(2, "Domain2")
    return geometry
mesh = Mesh(topBottom().GenerateMesh (maxh=0.1))
n = specialcf.normal(mesh.dim)

V = VectorL2(mesh, order=1, dgjumps=True)
Q = H1(mesh, order=1, definedon=mesh.Boundaries("gammaf"), dirichlet="gamma1|gamma2")
X=FESpace([V,Q])
(u,p),(v,q)=X.TnT()

dsf1=dx(element_boundary=True,definedon=mesh.Boundaries("gammaf"))
dsf2=dx(skeleton=True,definedon=mesh.Boundaries("gammaf"))

a=BilinearForm(X,check_unused=False)
a+=(p*(v-v.Other())*n+q*(u-u.Other())*n)*dsf2

_,_,vals = a.Assemble().mat.COO()
npvals = vals.NumPy()
print('assemble values = ', (npvals*npvals).sum())
2

Hi Dori,

you are integrating over an internal interface, and you want to access functions from the left, the right, and the interface, right ?
This is not fully supported in NGSolve. Internal interfaces are treated as faces inside the domain, and ‘u’ and ‘u.Other()’ refer to the left and right volume elements.

Way out:
Q1 = H1(mesh, order=1)
Q = Compress(Q1, active_dofs=Q1.GetDofs(mesh.Boundaries(“gammaf”)))

Then you can access functions from Q in the volume, but you pay only for dofs at the interface.

  • Joachim
3

Thank you for your reply.
I see that you set the order of Q1 to be 1. Can I choose a higher degree for my problem, because we need a higher-order discretization?

Also, another question: How do we implement the Dirichlet condition on Q at the endpoints of the interface gamma f?

Will this work the way we want:
Q1 = H1(mesh, order=order,dirichlet=“gamma1|gamma2”)

Q = Compress(Q1, active_dofs=Q1.GetDofs(mesh.Boundaries(“gammaf”)))
…
gfu.components[0].Set(Exact_p, definedon=mesh.Boundaries(“gamma1|gamma2”))

Thank you

4

Hi Dori,

yes, setting order to higher order is standard in NGSolve.

Setting inhomogeneous Dirichlet values at the end of the interface works as you suggested. An easier alternativ is to set values at the end of the interface, i.e. at an co-dimension 2 region. You can access it as the boundary of a boundary-region. Here a minimal example:

from ngsolve import *
from netgen.geom2d import SplineGeometry
from ngsolve.webgui import Draw
ngsglobals.msg_level=1

def topBottom():
    geometry = SplineGeometry()
    pnts     = [ (0, 0),  (0.5, 0), (1, 0), (1, 1), (0.5, 1), (0, 1)  ]
    pnums    = [geometry.AppendPoint(*p) for p in pnts]
    lines = [ (pnums[0], pnums[1], "gamma1",  1, 0),
              (pnums[1], pnums[2], "gamma2",  2, 0),
              (pnums[2], pnums[3], "gamma2",  2, 0),
              (pnums[3], pnums[4], "gamma2",  2, 0),
              (pnums[4], pnums[5], "gamma1",  1, 0),
              (pnums[5], pnums[0], "gamma1",  1, 0),
              (pnums[1], pnums[4], "gammaf", 1, 2)]
    for p1,p2,bc,left, right in lines:
        geometry.Append( ["line", p1, p2], bc=bc, leftdomain=left, rightdomain=right)
    geometry.SetMaterial(1, "Domain1")
    geometry.SetMaterial(2, "Domain2")
    return geometry
    
mesh = Mesh(topBottom().GenerateMesh (maxh=0.1))
n = specialcf.normal(mesh.dim)

print (mesh.Boundaries("gammaf").Mask())
print (mesh.Boundaries("gammaf").Boundaries().Mask())


order=3
Q1 = H1(mesh, order=order,dirichlet="gamma1|gamma2")
Q = Compress(Q1, active_dofs=Q1.GetDofs(mesh.Boundaries("gammaf")))

gfu = GridFunction(Q)
# both versions do work:
# gfu.Set(2*y-1, definedon=mesh.Boundaries("gamma1|gamma2"))
gfu.Set(2*y-1, definedon=mesh.Boundaries("gammaf").Boundaries())
Draw (gfu)

Best,
Joachim

5

Thank you for being so helpful.

Another issue I am having with the compressed space is the refine mesh.

I have a function called “solve” that calculates the solution for h=1/4. I then use a for loop to refine the mesh and call the “solve” function again.

I add
Q1.Update()
Q.Update()
but it still not working