Discrete/Interpolated Matrix Multiplication with Test/Trial Function in Variational Formulation

variational formulation
1

Is there a way to implement a variational formulation of the following form on the CG1 FEM space:

\frac{1}{k} (u^n - u^{n-1}, \phi)_2 + (\mathcal{I}_h (A u^n), M \mathcal{I}_h (A \phi))_2

for square Matrix Functions A,M and where we used the global Interpolation operator \mathcal{I}_h s. t.

\mathcal{I}_h (A \phi) = \sum_z A(z) \phi(z) \psi_z

for all nodes z with the according shape functions \psi_z.
This formulation should still be linear in the DOFs and therefore (hopefully) be automatically integrable.
Applying the Interpolation operator globally can be achieved using the quadrature options, so the formulation

\frac{1}{k} (u^n - u^{n-1}, \phi) + \int_\Omega \mathcal{I}_h [(A u^n) \cdot M B \phi ] dx

would be feasible, but is not what I am looking for exactly.

Any help is appreciated. :slight_smile:

2

Either I misunderstood completely, or the formulation isn’t linear, considering that there is a product of \phi times itself.

3

Sorry that was an unfortunate typo. Edited it accordingly. \psi_z are the basis or shape functions of the FEM Space whereas \phi is the test function.

4

Let’s start with the simple things: do you know how to write a weighted inner product (MatrixValuedFunction trial, test) in ufl?. The short version: you could interpolate your matrix-valued function into a tensor CG 1 space (say, call it M), and the simply write ufl.inner(M*trial, test).

See for instance dolfinx/python/test/unit/fem/test_bcs.py at v0.8.0 · FEniCS/dolfinx · GitHub for how to define a tensor function space, or many of the online tutorials linked in this forum. The line above would also tell you how to define a function space for a vector field since, I guess, your test and trial would be vector valued (otherwise I don’t understand how A and M could be matrix-valued functions).

I don’t fully understand the rest of the question. You call \mathcal{I} a global interpolation operator, but it seems to me that the result is still local to the cell (because either \phi or \psi would be zero unless they live on the same cell).

5

Thanks for the answer so far! Constructing the matrix is not the problem, but really the application of the nodal interpolation Operator \mathcal{I} (for CG it is global and local, it does not make a difference, so sorry for the confusing notation). I’ll try to make it more practical. Let us consider the Lagrangian CG1 FEM Space on the triangulation of a mesh \Omega = \cup_i K_i with the nodes N_h:

S_h = \{ v \in C^0: v|_K \in P^1 \}.

So lets take v\in S_h. Then it can represented in form of shape functions \psi_z

v = \sum_{z \in N_h} v_z \psi_z \quad s.t. v(z) = v_z , \quad \psi_z (\hat{z}) = \delta_{z, \hat{z}}.

These basis functions span the FEM space, so

span\{\psi_z : z \in N_h\} = S_h.

So we represent the FEM function v in terms of the Lagrangian basis functions \psi_z which are 1 at node z and 0 at all other nodes.
Now I want to manipulate the Trial or Test Function on a nodal level not in the L^2 sense.
So, for a Test, Trial Function u, \phi\in S_h and given Matrix FEM Functions A\in S_h I would like to implement the term

(\mathcal{I} (Au), \phi)_2 = \sum_{i,j} (Au) (z_i) \cdot \phi(z_j) \int_\Omega \psi_{z_i} \psi_{z_j} dx

which is different from

(Au, \phi)_2 = \sum_{i,j,k} A (z_k) u (z_i) \cdot \phi(z_j) \int_\Omega \psi_{z_k} \psi_{z_i} \psi_{z_j} dx.

As far as I understand and as far as I have used FEniCS(x), your suggestion would leave me with the latter term. However, for my particular case the nodal realization of that manipulation really is crucial.

I know that UFL has functions like as_matrix, as_vector, but I have only seen them being used with constants.

6

I constructed a MWE below for an example of the harmonic heat flow into the sphere with the weak formulation

(\partial_t d, \phi) - (d \times \Delta d, d\times \phi) =0

Hereby, I compute one iteration of three different variational formulations.

As matrix we simply consider the cross product in this example.

Case 1 corresponds to the standard L^2 formulation.

(\partial_t d, \phi)_h - (d \times \Delta d, d\times \phi)_2 =0

with

(f,g)_h = \int_\Omega \mathcal{I} ( f\cdot g) dx

This corresponds to

inner(d - d0, td)*dxL - dt*inner(cross(d0,q), cross(d0,td)) *dx

Case 3 corresponds to the case with mass-lumping applied.

(\partial_t d, \phi)_h - (d \times \Delta d, d\times \phi)_h =0 
inner(d - d0, td)*dxL - dt*inner(cross(d0,q) , cross(d0,td)) *dxL

Case 2 corresponds to the above described need of interpolating the Test Function with a matrix.

(\partial_t d, \phi)_h - (\mathcal{I} [d \times \Delta d], \mathcal{I} [d\times \phi ])_2 =0

I try to achieve this using the ufl Function as_vector. Since this is part of the ufl standard repertoire I however expect the same behaviour as in case 1 due to the assembly.

MWE:

import numpy as np

from dolfinx.fem import Function, functionspace, form, assemble_scalar, ElementMetaData
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import CellType, create_box, CellType
from dolfinx.io import XDMFFile
from ufl import dx, grad, inner, div, dot, TrialFunctions, TrialFunction, TestFunction, TestFunctions, lhs, rhs, Measure, as_vector, cross, VectorElement, MixedElement, split

from mpi4py import MPI

def L2error(a):
    error = form(inner(a, a) * dx)
    error_a = assemble_scalar(error)
    return error_a

def ic(x: np.ndarray)-> np.ndarray:
    # x hase shape (dimension, points)
    dim = 3
    values = np.zeros((dim, x.shape[1])) # values is going to be the output
    values[0]= np.sin( 2.0*np.pi*(np.cos(x[0])-np.sin(x[1]) ) )
    values[1]= np.cos( 2.0*np.pi*(np.cos(x[0])-np.sin(x[1]) ) )
    if dim == 3: values[2]=0.0

    return values

def main():
    dim         = 3
    dt          = 0.5 
    dh          = 10    
    domain = create_box(MPI.COMM_WORLD, [np.array([-1, -1,-1]), np.array([1.0, 1.0,1.0])],  [dh,dh,dh], cell_type = CellType.tetrahedron)
        
    el        = VectorElement("Lagrange", domain.ufl_cell(), 1, dim = dim)
    me        = MixedElement([el,el])
    FS        = functionspace(domain, me)

    d, q = TrialFunctions(FS) 
    td, tq = TestFunctions(FS)

    u1 = Function(FS)
    u2 = Function(FS)
    u3 = Function(FS)
    
    u0 = Function(FS)
    d0, q0 = split(u0)

    
    def g_times_f(g,f):
        """
        cross product as matrix
        https://en.wikipedia.org/wiki/Cross_product#Alternative_ways_to_compute
        """
        f = as_vector([
            -g[2]*f[1] + g[1]*f[2],
            g[2]*f[0] - g[0]*f[2],
            -g[1]*f[0]+g[0]*f[1]
            ])

        return f
    

    dxL = Measure("dx", metadata = {"quadrature_rule": "vertex"}) # mass lumping 

    # Example:
    # Harmonic Heat Flow into sphere

    energy_eq = inner( grad(d), grad(tq))*dx - inner(q,tq)*dxL

    # Case 1: L2 case (standard)
    eq1 = inner(d - d0, td)*dxL - dt*inner(cross(d0,q), cross(d0,td)) *dx  + energy_eq

    # Case 2: reformulation with ufl.as_vector
    # Expected Behaviour: Same as L2 case 1
    eq2 = inner(d - d0, td)*dxL - dt*inner(g_times_f(d0,q) , g_times_f(d0,td)) *dx  + energy_eq

    # Case 3: fully mass-lumped
    eq3 = inner(d - d0, td)*dxL - dt*inner(cross(d0,q) , cross(d0,td)) *dxL  + energy_eq

    a1=lhs(eq1)
    b1=rhs(eq1)

    a2=lhs(eq2)
    b2=rhs(eq2)

    a3=lhs(eq3)
    b3=rhs(eq3)

    u0.sub(0).interpolate(ic)

    solver_param = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"} 

    LP1 = LinearProblem(a1,b1 , bcs=[], u=u1, petsc_options=solver_param) 
    LP2 = LinearProblem(a2,b2 , bcs=[], u=u2, petsc_options=solver_param) 
    LP3 = LinearProblem(a3,b3 , bcs=[], u=u3, petsc_options=solver_param) 

    LP1.solve()  
    LP2.solve()
    LP3.solve()

    print("diff between 1 and 2", L2error(u1.sub(0).collapse()-u2.sub(0).collapse()))
    print("diff between 1 and 3", L2error(u1.sub(0).collapse()-u3.sub(0).collapse()))
    print("diff between 2 and 3", L2error(u3.sub(0).collapse()-u2.sub(0).collapse()))



if __name__ == '__main__':
    main()

The output is:

diff between 1 and 2 0.0
diff between 1 and 3 3.8015563928781586
diff between 2 and 3 3.8015563928781586

This shows that using the method as_vector is assembled in the same way as the standard ufl formulation. This was also my expected behaviour.

So I am still wondering if one can construct a Trial Function where the nodal coefficients can be scaled dependent on another function. So for a Test/Trial Function u that can be represented as

u = \sum_z u_z \psi_z,

I instead want to use

\mathcal{I}(Au) = \sum_z (A_z u_z) \psi_z,

as Test/Trial Function in the variational formulation.

I hope that can explain my question and the motivation behind it better. :sweat_smile: