specify intrules for 1d mesh

Jiashun
Topic Author
Offline
New Member

3 years 7 months ago #3709 by Jiashun

specify intrules for 1d mesh was created by Jiashun

Hello,

I am trying to generate the lumped mass matrix by setting integration rules.

For a 2d surface mesh in 3d, the following code works well.

Code:

ir = IntegrationRule(points = [(0,0), (1,0), (0,1)], weights = [1/6, 1/6, 1/6] )
ds_lumping = ds(intrules = { TRIG : ir })

However, for a 1d curved mesh in 2d, I guess it should be

Code:

ir = IntegrationRule(points = [(0,0), (1,0)], weights = [1/2, 1/2] )
ds_lumping = ds(intrules = { TRIG : ir })

However, it does not work.

Here is my code. I hope that the final result is a diagonal matrix.

Code:

import netgen.gui
import netgen.meshing as ngm
from ngsolve import *

p = [[1,0,0],[1,1,0],[0,1,0],[0,0,0]]
nr = len(p)
mymesh= ngm.Mesh(dim=2)
pids = []
for i in range(nr):
    pids.append (mymesh.Add (ngm.MeshPoint(ngm.Pnt(p[i]))))  

idx = mymesh.AddRegion("material", dim=1)
for i in range(nr):
    if i<nr-1:
        mymesh.Add(ngm.Element1D([pids[i],pids[i+1]],index=idx))
    else:
        mymesh.Add(ngm.Element1D([pids[i],pids[0]],index=idx))
    mesh = Mesh(mymesh)

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

ir = IntegrationRule(points = [(0,0), (1,0)], weights = [1/2,1/2])
ds_lumping = ds(intrules = { TRIG : ir }, definedon=mesh.Boundaries("material"))

a = BilinearForm(fes)
a += u*v*ds_lumping

a.Assemble()
print(a.mat)

Thank you very much for your time in advance.

mneunteufel
Offline
Premium Member

3 years 7 months ago #3713 by mneunteufel

Replied by mneunteufel on topic specify intrules for 1d mesh

Hi Jiashun,

for 1D Elements (in 2D) you need to specify the integration rule for segments, not for triangles. If you change the code to

Code:

ir = IntegrationRule(points = [(0,0), (1,0)], weights = [1/2,1/2])
ds_lumping = ds(intrules = { SEGM : ir }, definedon=mesh.Boundaries("material"))

then you should get a diagonal mass matrix.

Best
Michael

Jiashun
Topic Author
Offline
New Member

3 years 7 months ago #3714 by Jiashun

Replied by Jiashun on topic specify intrules for 1d mesh

It works! Thanks for your explanation.

Time to create page: 0.095 seconds