modification of CutFEM nodes

Dear ngsolve/ngsxfem community,

I want to modify the position of the nodes that are associated with the cut elements which intersect with level-set interface in the X/Cut-FEM, then update the modified mesh for my following calculation.

I see there is a way to modify the points (nodes) of the existed mesh on the side ngmesh. However, to modify the points, we should first get the information on which elements need to be modified. What I know is that in ngsxfem, CutInfo can be used to get such information, for example using ci.GetElementsOfType(IF). But it is in type BitArray. I do not know how to use that ...

So, could you give me some hint on doing so, thanks a lot !


Shaoqi

Hi Shaoqi,

the index in the BitArray corresponds to the element number:

from netgen.geom2d import unit_square
from ngsolve import *
from xfem import *

mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))
lset = (x-0.5)**2 + (y-0.5)**2 - 0.2**2

lsetp1 = GridFunction(H1(mesh, order=1))
InterpolateToP1(lset, lsetp1)
ci = CutInfo(mesh, lsetp1)

for el in mesh.Elements(VOL):
    print(el.nr)
    print ("vertices: ", el.vertices)   # get vertices of an element
    print ("edges: ", el.edges)

print('xxxxx')

for i in range(mesh.ne):
    if ci.GetElementsOfType(IF)[ i ]:
        print(i, 'is a cut element')
        el = mesh.ngmesh[ElementId(VOL, i)]    
        print ("vertices: ", el.vertices)   # get vertices of an element
        print ("edges: ", el.edges)

Best wishes,
Henry

Hi Henry,

Thanks lot your reply.
I've tried your advice, there are still some problem.
1. The object in mesh.Elements(VOL) is read only, right? So, we have to pass to mesh.ngmesh, where all the stuff can be modified.
However, this code mesh.ngmesh[ElementId(VOL, i)] seems not work: the argument not support.
2. To modify the position of the nodes, the ng.mesh.Points() should be used? or there are other functions?

Here is my simple test code (I just want move the the nodes associated with cut elements vertically for example:

map = lambda x,y: (x*0.7-0.5, y/10)
mesh = MakeStructured2DMesh(quads=False, nx=5, ny=2, mapping=map)

levelset = x-0.
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order-1, levelset=levelset)
deformation = lsetmeshadap.CalcDeformation(levelset)
lsetp1 = lsetmeshadap.lset_p1

ci = CutInfo(mesh, lsetp1)

for i in range(mesh.ne):
if ci.GetElementsOfType(IF):
print(i, 'is a cut element')
el = mesh.ngmesh[ElementId(VOL, i)] # here is a problem
p = [list(mesh[el.vertices].point) for i in range(3)]
print ("vertices: ", el.vertices) # get vertices of an element
print ("points associated: ", p)
for i in range(3):
p[0] +=0.05
print ("points associated: ", p)


Thanks again !

Hi Shaoqi,

1. Yes, from the ngsolve side, the mesh is read only. From the netgen side you can move points if you don't change the mesh topology. See  ngsolve.org/forum/ngspy-forum/1372-modifying-elements-vertices

2. As you are working with the netgen mesh, rather than the NGSolve mesh, you need to work with the netgen types. The following should do what you're after

vertices = [ ]

for i, el in enumerate(mesh.ngmesh.Elements2D()):
    if ci.GetElementsOfType(IF) [ i ]:
        print(el.points)
        
        for p in el.points:
            if p not in vertices:
                vertices.append(p)
                mp = mesh.ngmesh[p]
                mp[1] += 0.01
                mesh.ngmesh[p] = mp

mesh.ngmesh.Update()
Draw(mesh)

Hi Henry,

Your code works! Thanks a lot!

Shaoqi