Hi,
I’m trying to implement a DG-BEM coupling method. Therefore i would like to iterate over all the elements on the boundary to set up my coupling matrices. Everything works fine if I use a H^1 space, but if I use L^2 the resulting list is always empty.
To illustrate what I mean, I’ve attached a short sample script. If I change the space to H1 there, everthing works as expected.
https://ngsolve.org/media/kunena/attachments/922/test.py
Attachment: test.py
Hi Alexander,
The space L2 does not provide Surface Elements. However, we have SurfaceL2 to obtain an L2 space that lives only on the surface.
Best,
Christoph
Hi,
I guess what you need is the L2 function to be evaluated “at” the boundary. You can do that by iterating over mesh boundary elements and evaluating the L2 function (which will then be evaluated at the adjacent volume element) at your desired points.
Is there no way to get information on which degrees of freedom are attached to a surface element? Otherwise this sounds like a really unclean solution with lots of iterating to generate zeros.
To be more specific, what I would like to do is set up a sparse matrix representing the trace operator, but between the two libraries.
For the H1(and Maxwell) case, i iterate over the surface elements and for each element, i construct a local Dof map by evaluating the shapeset in the lagrange points (as Bem++ uses lagrange based elements). Then i use the local2global maps of both libraries to construct the global map.
While some solution based on a SurfaceL2 might be possible, I would like to keep my discretization of the surface values in the BEM++ format (because thats where the more involved computations are happening).
On a different, but related note: Is it possible in NGsolve to assemble a Mass matrix where test and trial functions are different FESpace? (one surface L^2 and one global L^2) BilinearForm seems to only allow test and trial space to be the same.
Thanks for your help,
Alex
BilinearForms allow different test and trial spaces
a = BilinearForm(trialspace=fes1, testspace=fes2)but combining surface and volume L2 spaces does not make much sense, as they do not have any elements in common (surface l2 is only defined on surface elements, L2 only defined on volume elements)
You can project the function to the surface L2 by using:
sgf.Set(BoundaryFromVolumeCF(gf), BND)and then do the dof mapping using the l2 surface elements.
Hi again,
Perhaps the following works for you?
fesf = FacetFESpace(mesh,order=…, dirichlet=“.*”)
fesf = Compress(fes,active_dofs=~fesf.FreeDofs())
m = Bilinearform(trialspace=fesf, testspace=fesl2)
m += fesf.TrialFunction() * fesl2.TestFunction() * ds(skeleton=True,BND)
Best,
Christoph
Hi,
using (a slighlty modified version) of your code produces a segfault
“malloc(): invalid size (unsorted)”.
https://ngsolve.org/media/kunena/attachments/922/test.py
@christopher:
BoundaryFromVolumeCF seems to roughly do what I want, except that I don’t have just a single GridFunction I have to transfer but I would really like to build the whole transfer matrix. (As doing the whole projection thing either for every dof or for every iteration in a solver seems overkill).
Maybe I can steal enough code from BoundaryFromVolumeCF to hack together a small C+±part that does what I want.
thanks for your help,
Alex
Attachment: test_2020-04-30.py
Do you only need the operator L2 → surface L2, or also surface L2 → L2?
The first is easy to do, the second is more difficult. I think for the second some C++ code is needed.
from ngsolve import *
import netgen.geom2d as geom2d
mesh = Mesh(geom2d.unit_square.GenerateMesh(maxh=0.5))
L = L2(mesh, order=0)
F = FacetFESpace(mesh, order=0, dirichlet=".*")
S = SurfaceL2(mesh, order=0)
E1 = comp.ConvertOperator(spacea=L, spaceb=F)
E2 = comp.ConvertOperator(spacea=F, spaceb=S, vb=BND)
E = E2 @ E1
print(E)This will give you the desired operator L2->surface L2 as a sparse matrix.
The operator “E1” will average on the interior facets, but we do not care about that as we are only interested in the boundary.
(I did have to add a single line of C++ code in one place to make this work, I will put up a merge request for this).
I do have an idea for the other direction, but it is a bit ugly and needs some C++ code.
You can check out the “set_dualshapes” branch from NGSolve and try the attached script.
For the Volume->Surface operator I did have to fix something in some code I wrote a while ago. That fix should come to the master soon.
For the Surface->Volume operator I did have to make one change to the “Discontinuous” spaces that Is kind of hacky, so I am not convinced it will come to the master branch. Also, it only works if you choose the volume L2 space to be at least of order 1.
Also, you will only get an exact extension of the surface space for functions that have an exact extension.
For example, if you have a triangle in a corner, and your surface function is 1 on one edge connected to it and 0 on the other, you can not have any volume l2 function with this trace, so you will get some interpolation of this (it will average the value in the corner).
If your surface functions are continuous, you could simply use an H1 space for the surface and then map the H1 function you get to an L2 function. No hacks needed for that.
Best,
Lukas
Attachment: vol2surf.py
Hi Lukas,
thanks that solution looks interesting. I think the direction L^2->surface L^2 is sufficient.
The one other thing I would need for my method is the same map L^2->surface L^2 but evaluating the normal derivative. Is this possible in any way?
Currently not quite, but it should be easy. Something like this is on my list of things to implement anyways, give me like an hour.
Best,
Lukas
Okay, that was easy. Should work with the “set_dualshapes” branch. This needs another small change I am putting up a merge request for
As before, we convert to a FacetFESpace, “F”, but this time instead of the L2 function itself, we use it’s normal derivative:
V2S_1 = comp.ConvertOperator(spacea=L, spaceb=F, trial_cf=InnerProduct(L.TrialFunction().Deriv(), specialcf.normal(mesh.dim)))Then we go to the boundary as before:
V2S_2 = comp.ConvertOperator(spacea=F, spaceb=S, vb=BND)Best, Lukas
Attachment: vol2surf2.py
Thanks Lukas, your code works like a charm!
No problem, glad to hear it!