Hello team!
I’m trying to build a (very simple) BilinearForm with different trial- and test spaces, specifically, two FESpaces where one is specified on a mesh that is a refinement of the other.
As a minimum working example, consider the following code:
import scipy.sparse as sp
from ngsolve.meshes import Make1DMesh
mesh = Make1DMesh(1)
refined_mesh = Make1DMesh(2)
X = H1(mesh, order=2)
Y = H1(refined_mesh, order=2)
u = X.TrialFunction()
v = Y.TestFunction()
A = BilinearForm(trialspace=X, testspace=Y)
A += u * v * dx
A.Assemble()
mat = sp.csr_matrix(A.mat.CSR()).toarray()
print(mat)
This yields the following output:
[[ 0.33333333 0.16666667 -0.04166667]
[ 0.16666667 0.33333333 -0.04166667]
[ 0. 0. 0. ]
[-0.04166667 -0.04166667 0.00833333]
[ 0. 0. 0. ]]
Now, I am fairly sure this matrix is incorrect. Consider the following test
import numpy as np
gfX = GridFunction(X)
gfY = GridFunction(Y)
gfX.vec.FV().NumPy()[0] = 1
gfY.vec.FV().NumPy()[0] = 1
assert np.isclose(mat[0, 0], Integrate(gfX * gfY, mesh_S))
which yields an AssertionError.
Maybe it’s not supported to create BilinearForms with different underlying meshes.
How could I achieve the same goal (i.e., working with two FESpaces with different underlying meshes) without manually integrating every basis function?
As an update: I do realize this is quite a “weird” thing to do. A solution may be to use the prolongation operator from the multigrid module if I manage to recast the definition of Y on a true refinement of the mesh – unfortunately, this may be hard in this 1D case because 1D mesh refinement is not supported (cf Segfault when creating 1D mesh - Kunena). Will report back whenever I get more answers!
Final update: LinearForms do work as intended (presumably because it uses quadrature under the hood), so I was able to “work around” this problem through a bit of craftiness. This code will be horribly slow so unsuitable for production, but it is enough for my proof of concept.
import scipy.sparse as sp
import numpy as np
from ngsolve.meshes import Make1DMesh
mesh = Make1DMesh(1)
mesh_S = Make1DMesh(2)
Y = H1(mesh_S, order=2)
X = H1(mesh, order=2)
v = Y.TestFunction()
gfs = [GridFunction(X) for _ in range(X.ndof)]
mat = np.zeros([Y.ndof, X.ndof])
for i, gf in enumerate(gfs):
gf.vec.FV().NumPy()[i] = 1.0
l = LinearForm(Y)
l += gf * v * dx
l.Assemble()
mat[:, i] = l.vec.FV().NumPy()[:]
print(mat)
gfX = GridFunction(X)
gfY = GridFunction(Y)
gfX.vec.FV().NumPy()[0] = 1
gfY.vec.FV().NumPy()[0] = 1
assert np.isclose(mat[0, 0], Integrate(gfX * gfY, mesh_S))
Hm yes, the mixed bilinearform expects to be defined only on different spaces and not on different meshes. If one is only a refinement of the other that may be working to only hack some parts of the c++ code (but still not easy), but for completely different meshes it is difficult…
Can you explain the application for this?
Hello Christopher,
Thanks for your reply. My application is a mixed higher-order space-time formulation of a parabolic PDE which needs a higher order H(div)-conforming function space, which was exactly our reason for choosing NGSolve.
I was unable to get a “clean” solution working (using a TensorProductFESpace) in Python, and I was unable to build NGSolve from source so I could work in C++, so I resorted to building my necessary matrices in space and time separately, and build the space-time system matrix using Kronecker products. In building one of my space matrices, I hit the bug mentioned above.
Indeed, one mesh is a refinement of the other, so my hope was to somehow get the prolongation operator from coarser to finer mesh and be on my merry way. It seems that this method is not exposed in the Python API (only a ProlongateInline() method, which I never got to work). I was able to work around this (in 2D at least) by creating a number of GridFunctions with setting nested=True and unit vectors as data, refining the mesh, Update()'ing the GridFunctions, and reading the resulting vectors, but this only seemed to work for lowest-order spaces.
I would love to get into the nitty gritty of NGSolve using C++, however I was unable to build it from source on both of my machines (one OSX 10.14, one OSX 10.12). If I manage to fix this problem, my next NGSolve project will work with C++ under the hood 
Sorry for the lengthy reply!
Hi Jan,
Which problems did you encounter when building NGSolve from source? Did you follow the instructions in the documentation here?
Best,
Matthias
Hey Matthias,
Yes I am following this guide. The problem is (for both machines) at step Build on Mac OS X. — NGS-Py 6.2.2302 documentation
Some of the output I am getting is collected in this gist: NGSolve error building from source · GitHub
[please read shell_output first, and the rest afterwards.]
I tried running make with both Tcl and Tcl/Tk installed and uninstalled. Do you have any idea?
Hi Jan,
A first guess: Did you remove the dmg bundle before starting the configuration?
Ha good catch! I certainly did not.
This got me a little further, it is now crashing at project_tk rather than project_tcl: WITHOUT_TK shell output · GitHub
when I brew uninstall tcl-tk and rerun, I get a different error – see WITH_TK* and WITHOUT_TK* in the gist linked above.
That’s interesting, NGSolve should build and link it’s own version of Tcl/Tk, so brew is not necessary at all.
Try to wipe /Applications/Netgen.app and your build directory completely and run cmake from scratch.
Hey Matthias,
Thanks a lot – I got it working by updating my clang to 10.0 for c++17 (required by -i think- pybind11). I now have my own build from source :-)!!