Nodal interpolation of inner product of vector valued functions


Considering V=VectorH1(mesh, order=1) and Q=H1(mesh, order=1), I am trying to assemble a bilinear form
B = BilinearForm(trialspace=V, testspace=Q) with the form b(u,q)=\int q I_h(u \cdot u_0). Here, I_h denotes the nodal interpolation into the space Q, and u_0 is a given GridFunction in V.

I tried to compute the interpolation I_h(u \cdot u_0) as in 2.10 of the tutorial (

defining nodal interpolation

def nodal_interpolation(func,fes):
gfu = GridFunction(fes)
u,v = fes.TnT()
vdual = v.Operator(“dual”)
a_int = BilinearForm(fes)
a_int += uvdualdx(element_vb=BBND)
f_int = LinearForm(fes)
f_int += funcvdualdx(element_vb=BBND)
f_int.Assemble() = a_int.mat.Inverse() * f_int.vec
return gfu

u = V.TrialFunction()
q = Q.TestFunction()
w = nodal_interpolation(u_0*u,Q)

However, I got a “segmentation fault”. I guess this is because u is TrialFunction function in V, but in nodal_interpolation f_int is a LinearForm(Q) so that it is not allowed to call a TrialFunction of another space here.

I wonder how to define a nodal_interpolation(func, fes), where “func” contains a TrialFunction function of a space different from “fes”? More specially, what would be best way to define a nodal interpolation of inner product of vector valued TrialFunction and vector valued GridFunction into the space of scalar valued function?

Many thanks in advance.


it depends a bit what you really need …

In the simplest, low order case numerical integration with points in the vertices does the job, use a user-defined IntegrationRule.

In the second part of tutorial 2.10 we define a linear operator to interpolate a function from one space to another. Similar to this, you can implement an interpolation operator of u_0*u to a high order H1-space.

Here is also a related thread:


The operator
u\in V \rightarrow I_h(u\cdot u_0) \in Q
itself can be set up conveniently by

Ih = ConvertOperator(spacea=V, spaceb=Q, trial_cf=InnerProduct(u, u_0))

So if you just want to compute the interpolation you can use this.


Hi Joachim,

Thanks for your reply, and thanks for pointing to the related thread Nodal interpolant available? - Kunena! I think the idea of introducing an auxiliary variable discussed there should also work for my problem.

However, I am wondering if there is an easier way of doing this by any chance. So I have the following two questions:

  1. I am trying to assemble bilinear form B = BilinearForm(trialspace=V, testspace=Q) from the integral:

b(u,q)=\int q I_h(u \cdot u_0).

If the integral of the form \int I_h(something), I can see that numerical integration with points in the vertices does the job. However, in my case, by taking trial function u to be nodal basis \Phi_i in V and test function q to be nodal basis \Psi_j in Q, I can simplify:

\int \Psi_j I_h(\Phi_i \cdot u_0) = \int \Psi_j (\Phi_i \cdot c_i),

where c_i is a constant vector determined by nodal value of u_0 at x_i. Is it possible to use the symbolic integration to compute B_ij = \int \Psi_j (\Phi_i \cdot c_i)?

  1. I also applied the idea of tutorial 2.10 to define an operator nodal_interpolation(func1,func2,fes1,fes2) interpolating func1\cdot func2 into fes2, where func1 and func2 belong to fes1. This works fine for the cases when func1 is GridFunction or CoefficientFunction. However, I need to apply this operator to a TrialFunction func1:

u,v = V.TrialFunction(), V.TestFunction()
p,q = Q.TrialFunction(), Q.TestFunction()
u_old = GridFunction(Vh)
w = nodal_interpolation(u,u_old,V,Q)
b = BilinearForm(trialspace=V, testspace=Q)
b += q * w * dx

but it seems that b.mat has all elements zero, and thus applying this interpolation operator to TrialFunction and using it in symbolic integration does not work. Is it possible to apply an interpolation to a TrialFunction and use it directly in symbolic integration?

Many thanks in advance.


we have implemented now the possibility to interpolate expressions into some finite element space, and use it within matrix assembling.

For example, you can interpolate u*v into the nodal P1 to obtain mass lumping (see attachment).

Hope this helps in your case.

BUT, A WARNING: This is just a quick proof of concepts implementation, and orders of magnitudes slower than the usual NGSolve assembling. If you have a working formulation using the Interpolate - function, you can post the example to discuss an efficient realization.

The feature is available within the available nightly release.

Best, Joachim