Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

Zero matrix

  • Dori
  • New Member
  • New Member
More
11 months 4 weeks ago #4785 by Dori
Zero matrix was created by Dori
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())
More
11 months 4 weeks ago #4787 by joachim
Replied by joachim on topic Zero matrix
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


 
Time to create page: 0.110 seconds