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 !
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 ][i][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)
Thanks lot your reply.
I’ve tried your advice, there are still some problem.
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.
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:
for i in range(mesh.ne):
if ci.GetElementsOfType(IF)[i]:
print(i, ‘is a cut element’)
el = mesh.ngmesh[ElementId(VOL, i)] # here is a problem
p = [list(mesh[el.vertices[i]].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[i][0] +=0.05
print ("points associated: ", p)
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 Modifying element's vertices - Kunena
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 ][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