Hello everyone,
I have a cut boundary that I have created using level set function (I am using xfem extension of NGsolve). My domain is square and it has a circular cut boundary in it defined by x^2+y^2=0.5^2. After doing my calculations, I want to calculate flux over this circular cut boundary which means calculating a boundary integral. After my previous question I know how to calculate integrals on mesh boundaries but I do not know how to calculate on this cut boundary, I keep getting wrong results or zeros. Analogous to my previous question I tried,
n_ls = grad(level_set_p1)* (-1/grad(level_set_p1).Norm())
VH1= H1(mesh, order=2, dirichlet="wall|inlet|outlet")
gradu00 = GridFunction(VH1)
gradu01 = GridFunction(VH1)
gradu10 = GridFunction(VH1)
gradu11 = GridFunction(VH1)
gradu00.Set(grad(u_h.components[0])[0])
gradu01.Set(grad(u_h.components[0])[1])
gradu10.Set(grad(u_h.components[1])[0])
gradu11.Set(grad(u_h.components[1])[1])
et = GridFunction(VH1)
et.Set(CoefficientFunction(1.0), definedon=level_set_cut)
flux_x = Integrate(levelset_domain = { "levelset" : level_set_p1, "domain_type" : IF},
cf=et*(gradu00*n_ls[0] + gradu01*n_ls[1]), mesh=mesh, order = order)
flux_y = Integrate(levelset_domain = { "levelset" : level_set_p1, "domain_type" : IF},
cf=et*(gradu10*n_ls[0] + gradu11*n_ls[1]), mesh=mesh, order = order)
print ("flux_x:", flux_x)