Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

Integral on interface between two domains

More
5 years 5 months ago #737 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
Attachments:
More
5 years 5 months ago #744 by joachim
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
Code:
for el in mesh.Elements(BND): if el.mat == "gammaSD": for edge in el.edges: interface_indicator[edge.nr] = True

Joachim
More
5 years 5 months ago - 5 years 5 months ago #748 by rhebergens
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
Code:
dgjumps=True
is added to the FESpace. Setting
Code:
dgjumps=False
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
Attachments:
Last edit: 5 years 5 months ago by rhebergens.
More
5 years 5 months ago #749 by joachim
Hi Sander,

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"))
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
More
5 years 5 months ago #750 by rhebergens
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
More
5 years 5 months ago #751 by joachim
overseeing the Trace() operator WAS a common pitfall.
I just added a test whether the trial and test-functions match with the form.
Joachim
Time to create page: 0.106 seconds