- Thank you received: 0
Integral on interface between two domains
- rhebergens
- Topic Author
- Offline
- Junior Member
Less
More
6 years 3 months ago #737
by rhebergens
Integral on interface between two domains was created by rhebergens
Hello,
I have a domain that is split in two (using the code by Guosheng:
ngsolve.org/forum/ngspy-forum/114-interface-integrals#491 )
In fluid 1 I define a FacetFESpace Q1 and in fluid 2 I define a FacetFESpace Q2.
Let p1 and q1 be test and trial functions in Q1.
Let p2 and q2 be test and trial functions in Q2.
On the interface between the fluids I compute the integral
100*p1*q1 + 100*p2*q2.
I then print the matrix and I see that I get different matrices depending on how I define Q2. If I define
Q2 = FacetFESpace(mesh, order=1)
I will get, for example, that A(54,54) = 16.6667, but if I define
Q2 = FacetFESpace(mesh, order=1, definedon=mesh.Materials("fluid2"))
I will get that A(54,54) = 0
I would prefer to use the second definition of Q2, but this gives the wrong matrix. Any suggestions how to fix this?
I have attached the code with Q2 defined in lines 52 and 53.
Thanks,
Sander
I have a domain that is split in two (using the code by Guosheng:
ngsolve.org/forum/ngspy-forum/114-interface-integrals#491 )
In fluid 1 I define a FacetFESpace Q1 and in fluid 2 I define a FacetFESpace Q2.
Let p1 and q1 be test and trial functions in Q1.
Let p2 and q2 be test and trial functions in Q2.
On the interface between the fluids I compute the integral
100*p1*q1 + 100*p2*q2.
I then print the matrix and I see that I get different matrices depending on how I define Q2. If I define
Q2 = FacetFESpace(mesh, order=1)
I will get, for example, that A(54,54) = 16.6667, but if I define
Q2 = FacetFESpace(mesh, order=1, definedon=mesh.Materials("fluid2"))
I will get that A(54,54) = 0
I would prefer to use the second definition of Q2, but this gives the wrong matrix. Any suggestions how to fix this?
I have attached the code with Q2 defined in lines 52 and 53.
Thanks,
Sander
Attachments:
6 years 3 months ago #744
by joachim
Replied by joachim on topic Integral on interface between two domains
Hi Sander,
using HDG plus additionally dgjumps looks like an expensive method - do you really want that ?
I can help with the HDG, but will not spend time on the dgjumps.
A note on the side: You can directly iterate over the mesh topology to mark edges etc.
There is a quite new section in the i-tutorials:
ngsolve.org/docu/nightly/i-tutorials/uni...gy/meshtopology.html
In your case
Joachim
using HDG plus additionally dgjumps looks like an expensive method - do you really want that ?
I can help with the HDG, but will not spend time on the dgjumps.
A note on the side: You can directly iterate over the mesh topology to mark edges etc.
There is a quite new section in the i-tutorials:
ngsolve.org/docu/nightly/i-tutorials/uni...gy/meshtopology.html
In your case
Code:
for el in mesh.Elements(BND):
if el.mat == "gammaSD":
for edge in el.edges:
interface_indicator[edge.nr] = True
Joachim
- rhebergens
- Topic Author
- Offline
- Junior Member
Less
More
- Thank you received: 0
6 years 3 months ago - 6 years 3 months ago #748
by rhebergens
Replied by rhebergens on topic Integral on interface between two domains
Hi Joachim,
Thanks for the suggestion regarding the new i-tutorials. I will update my codes.
Regarding the dgjumps, the reason I'm using dgjumps is because of the Stokes-Darcy coupled code by Guosheng. I have attached a slightly modified version. Here
is added to the FESpace. Setting
leeds to errors of the form
catch in AssembleBilinearform 2: SparseMatrixTM::AddElementMatrix: illegal dnums
Traceback (most recent call last):
File "<string>", line 204, in <module>
File "<string>", line 191, in SolveBVP
RuntimeError: SparseMatrixTM::AddElementMatrix: illegal dnumsin Assemble BilinearForm 'biform_from_py'
Do you know why this is? Since this is HDG I would have expected dgjumps not to be necessary.
Regards,
Sander
Thanks for the suggestion regarding the new i-tutorials. I will update my codes.
Regarding the dgjumps, the reason I'm using dgjumps is because of the Stokes-Darcy coupled code by Guosheng. I have attached a slightly modified version. Here
Code:
dgjumps=True
Code:
dgjumps=False
catch in AssembleBilinearform 2: SparseMatrixTM::AddElementMatrix: illegal dnums
Traceback (most recent call last):
File "<string>", line 204, in <module>
File "<string>", line 191, in SolveBVP
RuntimeError: SparseMatrixTM::AddElementMatrix: illegal dnumsin Assemble BilinearForm 'biform_from_py'
Do you know why this is? Since this is HDG I would have expected dgjumps not to be necessary.
Regards,
Sander
Attachments:
Last edit: 6 years 3 months ago by rhebergens.
6 years 3 months ago #749
by joachim
Replied by joachim on topic Integral on interface between two domains
Hi Sander,
I think what you want is this:
The trace of the vector-facet variable can be evaluated on the boundary. It lives already in the tangential space, so you don't need the projection. You have to provide the Trace() operator.
You don't need the dgjumps, which would lead to unnecessarily many matrix entries.
Best, Joachim
I think what you want is this:
Code:
interfaceTerm = 0.5*(1+4*pi*pi)*vbar.Trace()*ubar.Trace()
a += SymbolicBFI (interfaceTerm, definedon=mesh.Boundaries("gammaSD"))
You don't need the dgjumps, which would lead to unnecessarily many matrix entries.
Best, Joachim
- rhebergens
- Topic Author
- Offline
- Junior Member
Less
More
- Thank you received: 0
6 years 3 months ago #750
by rhebergens
Replied by rhebergens on topic Integral on interface between two domains
Hi Joachim,
The .Trace() was indeed the missing piece. It also answered my original question. I attach the now working version with the updates you suggested.
Thanks!
Sander
The .Trace() was indeed the missing piece. It also answered my original question. I attach the now working version with the updates you suggested.
Thanks!
Sander
Attachments:
6 years 3 months ago #751
by joachim
Replied by joachim on topic Integral on interface between two domains
overseeing the Trace() operator WAS a common pitfall.
I just added a test whether the trial and test-functions match with the form.
Joachim
I just added a test whether the trial and test-functions match with the form.
Joachim
Time to create page: 0.116 seconds