Interface Indicator and computing integral on interface of 2 subdomain

Main Forum Old-Forum - do not post here
1

I’m trying to compute the integral on the set of facets that lie on the interface of two subdomain. I found a code on the Ngsolve forum which set up the interface indicator as follows. The file was attached.

[code]def GetInterfaceIndicator(mesh):
dom_indicator = GridFunction(L2(mesh,order=0))
dom_indicator.Set(CoefficientFunction([1,2]))

ffes = FacetFESpace(mesh,order=0)
dom_ind_average = GridFunction(ffes)
mf = BilinearForm(ffes,symmetric=True)
mf += SymbolicBFI(ffes.TrialFunction()ffes.TestFunction(),element_boundary=True)
mf.Assemble()
ff = LinearForm(ffes)
ff += SymbolicLFI(dom_indicatorffes.TestFunction(),element_boundary=True)
ff.Assemble()
dom_ind_average.vec.data = mf.mat.Inverse() * ff.vec

interface_indicator = BitArray(ffes.ndof)
for i in range(len(dom_ind_average.vec)):
if (dom_ind_average.vec[i] > 1.75) or (dom_ind_average.vec[i] < 1.25):
interface_indicator[i] = False
else:
interface_indicator[i] = True
print(“facets marked as interface facets:”, interface_indicator)
return interface_indicator
[/code]

Could you please explain what the dom_ind_average is and why we have 1.75 and 1.25 in the code?

if (dom_ind_average.vec[i] > 1.75) or (dom_ind_average.vec[i] < 1.25):
Does Ngsolve have any functions about setting the similar interfaces, or any literature where I can find some related information?
Thank you so much.

Updated: Source of the attached file

Attachment: 2domain.py

2

Hi dong,

please post also the link if you refer to a different post.

If you reorder terms of the bilinear form
[tex]\sum_T\int_{\partial T}uh,vh,ds [/tex]
it is
[tex]\sum_E\int_{\partial E}2uh,vh,ds [/tex]
for the inner edges E. The right hand side is
[tex]\sum_E\int_{\partial E}(do_L+do_R),vh,ds [/tex]
where do_L is the value on the left element and do_R the value on the right one.

Thus if you devide by two you get (in strong form)
[tex]uh = 0.5(do_L+du_R) \text{ on every inner edge}[/tex]

and 0.5*(1+1) = 1, 0.5*(2+2) = 2, and 0.5*(1+2) = 1.5 and thus the values 1.25 and 1.75.

However, if you label your boundaries with names you can get the dofs with

fesff = FacetFESpace(mesh,order=0)
interface_indicator2 = BitArray(fesff.ndof)
interface_indicator2[:] = False
interface_indicator2 |= fesff.GetDofs(mesh.Boundaries("2"))

I would highly recommend labelling all boundaries and materials by name. This makes the code more readable and is also less error-prone (if you give them meaningful names), see attached file.

Best
Michael

https://ngsolve.org/media/kunena/attachments/889/2domain.py

Attachment: 2domain_2020-05-26.py

3

Thank you so much for clarification.
I updated the link of the attached file in my post.
It seems to me that it came from an old version of ngsolve forum.
By. the way, here is that link.
https://sourceforge.net/p/ngsolve/discussion/ngs-py/thread/5bb706f6/a370