BilinearForm incorrect?

JanWesterdiep
Topic Author
Offline
Junior Member

5 years 3 weeks ago - 5 years 3 weeks ago #2008 by JanWesterdiep

BilinearForm incorrect? was created by JanWesterdiep

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:

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:

Code:

[[ 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

Code:

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?

Last edit: 5 years 3 weeks ago by JanWesterdiep.

JanWesterdiep
Topic Author
Offline
Junior Member

5 years 3 weeks ago #2009 by JanWesterdiep

Replied by JanWesterdiep on topic BilinearForm incorrect?

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 ngsolve.org/forum/ngspy-forum/877-segfau...hen-creating-1d-mesh ). Will report back whenever I get more answers!

JanWesterdiep
Topic Author
Offline
Junior Member

5 years 3 weeks ago #2016 by JanWesterdiep

Replied by JanWesterdiep on topic BilinearForm incorrect?

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.

Code:

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))

christopher
Offline
Administrator

5 years 3 weeks ago #2018 by christopher

Replied by christopher on topic BilinearForm incorrect?

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?

JanWesterdiep
Topic Author
Offline
Junior Member

5 years 3 weeks ago #2021 by JanWesterdiep

Replied by JanWesterdiep on topic BilinearForm incorrect?

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!

matthiash
Offline
Administrator

5 years 3 weeks ago #2023 by matthiash

Replied by matthiash on topic BilinearForm incorrect?

Hi Jan,

Which problems did you encounter when building NGSolve from source? Did you follow the instructions in the documentation here ?

Best,
Matthias

Time to create page: 0.099 seconds