Zero matrix

  • Dori
  • New Member
  • New Member
More
1 year 2 months 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
1 year 2 months 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.158 seconds