Obtain Integration Points for ngsxfem model

NGSolve
1

Dear all,

I tried to obtain integration points for a ngsxfem model which I could then use as a reference in an integral computation wrt. another code. I would be perfect to obtain the integration points element-wise and with additional information whether the respective element is cut or not.

However, I struggled to find the integration points anywhere. In the “Integrate” function in “init.py” in xfem, I found this “Integrate_X_special_args” which could get “ip_container” as an input. Unfortunately, I think this function is not working, but I think the purpose of it is what I am looking for.

Is there a way to obtain integration points for a CutFEM problem?

I would appreciate any help :slight_smile:

As I first model, I wanted to integrate simply a constant function. My model would look like:

#import netgen.gui
# ngsolve stuff
from ngsolve import *
from ngsolve.webgui import *
# basic geometry features (for the background mesh)
from netgen.occ import *
# visualization stuff
from ngsolve.internal import *
# basic xfem functionality
from xfem import *

geo = OCCGeometry(unit_square_shape.Scale((0,0,0),3).Move((-1.5,-1.5,0)), dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=0.5,quad_dominated=True))
levelset = (sqrt(sqrt(x**2+y**2)) - 1.0)
lsetp1 = GridFunction(H1(mesh,order=1))
InterpolateToP1(levelset,lsetp1)
f = CoefficientFunction(1)
order = 2
dx = dCut(levelset=lsetp1, domain_type=NEG, order=order)
integral = Integrate(f * dx, mesh=mesh)

Best,
Michael Loibl

2

Hi Michael ,

I am not sure how to get the integration points when using an integration symbol like dCut, but with the levelset_domain dictionary you can do it like this:

IntPoints = []
Integrate(levelset_domain = { "levelset" : lsetp1, "domain_type" : IF}, cf=0, mesh = mesh, ip_container = IntPoints, order=1)

The coefficientfunction is irrelevant, so I always choose cf=0 for a a fast computation.

Best wishes,
Paul

3

Hey Michael,

yes I see. The general understanding is that we want to somewhat allow for an export of the ngsxfem integration points, whilst we do not regard this as an important core feature of ngsxfem. This is why there is no support for this in the context of the newer Integration Symbols dx/ dCut etc. However, as Paul pointed out, with the older interface of the levelset dict in Integrate, the array can be brought in. The items in the ip_container will be MeshPoints, so .nr shall give the relevant element number. The information about cut topology could be only reconstructed with high efforts, but I would suggest to construct a CutInfo object to obtain it.

Hope this helps, if you need further specific code snippets etc let us know.
Best,
Fabian