Hello!
I want to transfer a solution gridfunction to a new gridfunction which is based on a scaled up mesh.
The orientation and number of elements in the mesh do not change. I am using the ngsolve.Mesh.Scale() attribute to do so.
I tried doing this with the following two approaches:
from ngsolve import *
import netgen.meshing as ngmsh
ngmesh = ngmsh.Mesh()
ngmesh.Load("channel_3bcs.vol")
mesh1 = Mesh(ngmesh)
ngmesh2 = ngmsh.Mesh()
ngmesh2.Load("channel_3bcs.vol")
mesh2 = Mesh(ngmesh2)
mesh2.ngmesh.Scale(2)
fes1 = HDiv(mesh1, order = 3, autoupdate = True)
gfu1 = GridFunction(fes1, autoupdate = True)
fes2 = HDiv(mesh2, order = 3, autoupdate = True)
gfu2 = GridFunction(fes2, autoupdate = True)
gfu1.Set((x, 0))
Method 1: Set()
gfu2.Set(gfu1)
Result: The solution is set only in the region where the original mesh is (a quarter of the scaled up mesh in this case) and the rest of it is set to 0
Method 2: gfu.vec.data
gfu2.vec.data = gfu1.vec.data
Result: This covers the entire scaled up mesh but for HDiv FESpace, the resultant solution is scaled by a factor of the length scale with which the mesh is scale (so halved in this case)
This method works well for L2 and H1 spaces however.
Can you please explain what I’m doing wrong here?
Can you suggest a way i can transfer the solution to a scaled up mesh from a gfu obtained on the original mesh?
Thanks!
Hi,
both approaches seem to work as expected.
Setting with
gfu2.vec.data = gfu1.vec
works for H1 and L2 as the mappings between elements (like reference and physical) are simply by composition. For HDiv-elements the Piola transformation is used, which is not independent of the mesh-size. Have you tried out to set the same function
gfu2.Set((x,0))
and compare the entries of gfu2.vec with gfu1.vec if the difference is exactly the (inverted) factor?
Best,
Michael
You can do what you want to do using the CoordinateTrafo function:
from ngsolve.fem import CoordinateTrafo
trafo = CoordinateTrafo(CF((0.5*x, 0.5*y)), mesh1.Materials(".*"))
gfu2.Set(gfu1(trafo))
Best
Christopher
Hello Christopher!!
Thank you for your reply!!
CoordinateTrafo worked perfectly for HDiv FES! Appreciate your help!
However I am running into this error when working on CompoundFESpace
netgen.libngpy._meshing.NgException: CompoundFESpace does not have an evaluator for VOL!
on the line:
gfu2.Set(gfu1(trafo))
My FESpace has two components: One is HDiv based and the other is L2
With this you can transfer each function. But since for the L2 also setting the vector is enough this would be more efficient:
trafo = CoordinateTrafo(CF((0.5*x, 0.5*y)), mesh1.Materials(".*"))
gfu2.components[0].Set(gfu1.components[0](trafo))
gfu2.components[1].vec.data = gfu1.components[1].vec
assuming HDiv is component 0 and L2 is component 1
Hello!!
On implementing component wise Set() as follows,
gfu2.components[0].Set(gfu1.components[0](trafo))
I am getting
Segmentation fault (core dumped)
error. Any ideas what the issue might be? Appreciate your help a lot, thanks!!
Figured out the Segmentation Fault issue, I was introducing an extra factor in CoordinateTrafo by mistake.
It works now, but the Set() method seems to introduce some inaccuracies. I tried using Set with increased numerical integration order but there seems to be some cap to the accuracy of Set()
I am attaching some screenshots below to show what exactly is happening
This is the Y velocity field which is 0 in the solution, but it gets some non-zero value when Set on a scaled up HDiv based mesh.
Hello @christopher!
PFA the sample code and mesh I am working with for reference. Thanks!!
sample_code.zip (13.0 KB)
With no bonus_intorder option in Set (Default):
With bonus_intorder = 3:
Expected result is (x,0) on the domain, so magnitude ranging from 0 to 1.0
Ah sorry, my bad. The CoordinateTrafo currently only works for scalar functions, I will fix that for future NGSolve releases. For now you can apply it on each vector entry:
from ngsolve.fem import CoordinateTrafo
trafo = CoordinateTrafo(CF((0.5*x, 0.5*y)), mesh1.Materials(".*"))
gfu2.Set((gfu1[0](trafo), gfu1[1](trafo)))
Hello @christopher!! This works perfectly!
Thank you for your response!! 
The interpolation is now working as expected for H(div), the fix for vector-valued functions and CoordinateTrafo is fixed in the latest nightly,
Joachim
