- Thank you received: 1
Interface Indicator and computing integral on interface of 2 subdomain
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.
Could you please explain what the dom_ind_average is and why we have 1.75 and 1.25 in the code?
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
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):
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
- mneunteufel
- Offline
- Premium Member
Less
More
- Thank you received: 59
4 years 5 months ago #2734
by mneunteufel
Replied by mneunteufel on topic Interface Indicator and computing integral on interface of 2 subdomain
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
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
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
Attachments:
4 years 5 months ago - 4 years 5 months ago #2736
by dong
Replied by dong on topic Interface Indicator and computing integral on interface of 2 subdomain
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
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