Interface Indicator and computing integral on interface of 2 subdomain

More
4 years 5 months ago - 4 years 5 months ago #2733 by dong
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_indicator*ffes.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

Could you please explain what the dom_ind_average is and why we have 1.75 and 1.25 in the code?
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
Attachments:
Last edit: 4 years 5 months ago by dong. Reason: Source of the attached file
More
4 years 5 months ago #2734 by mneunteufel
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
Code:
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

File Attachment:

File Name: 2domain_2020-05-26.py
File Size:4 KB
More
4 years 5 months ago - 4 years 5 months ago #2736 by dong
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.
sourceforge.net/p/ngsolve/discussion/ngs...thread/5bb706f6/a370
Last edit: 4 years 5 months ago by dong.
Time to create page: 0.134 seconds