Cartesian product of square meshes?

mesh
1

Hello. I would like to solve a PDE with solutions like u\colon (0,1)^2\times (0,1)^2 \to \mathbb{R}. Attempting to use TensorProductMesh in ufl gives the error “Product of vector-valued elements not supported”. So is it at all possible to work with domains in more than 3 dimensions? Are there any known hacks for dealing with a Cartesian product of 2-dimensional meshes?

2

Is there a particular reason for wanting to do a cartiesian product of the meshes, and not make a suitable function space based on VectorElement or TensorElement?
At least in your example, the grids are identical in your cartesian product, is this supposed to be the case in all of your applications?

3

Hi dokken, thanks for the reply. Yes, the intended purpose for this application is to have a 2-dimensional mesh modeling a habitat, and u is a function of two individual’s spatial locations within that habitat. So in all instances the Cartesian product will have identical grids. Also the PDE is of the form

L_x u(x,y) + L_y u(x,y) = c

where c is constant, L_x is a differential operator acting only on the x slot of u and L_y is the same operator acting only on the y slot.

I might be mistaken but I thought TensorElements are used for tensor-valued functions, but here u is a scalar-valued function with tensor arguments. I don’t know how I would go about making a function space containing functions that are evaluated using coordinates like [[x1, x2], [y1, y2]]. If I could do that I’d be able to handle the rest of the implementation, I think.

4

Hi Matt,

I have a similarly structured problem where I have a scalar equation of a 4-D space with the 4-D space being a cartesian product between 2 2-D spaces. Were you ever able to find any solutions to this problem?

Thanks,
Mark

5

Hi Mark,

I do have a partial solution for this by using tensor products, but have not figured out how to formulate the boundary conditions on the product space. For my particular problem the boundary is the diagonal in the product space.

Have you made any progress yourself? Perhaps we could collaborate a bit?

Best,
Matt

6

Hi Matt,

Thanks for the reply. I looked into a workaround by using an operator splitting approach but haven’t really done much besides paper and pencil work. Recently I have just been working on some smaller test cases for my problem where I can just reduce to 3d. Is it straightforward in fenics to make, say a 4-d basis function, by doing a tensor product of 2 2d guys? I tried a little bit but didn’t really have any luck. I’d love to talk about this with you sometime via zoom or whatever is convenient for you. My email is markloveland@utexas.edu if you’re interested in looking into this.

Thanks,
Mark

7

I wouldn’t say it’s straightforward, but I do think it is doable! I actually have an example solving Poisson’s equation on the unit square using vector-valued functions on the unit interval, which is kind of a proof of concept for the 4D case. We should definitely discuss some of this over Zoom but for now I’ll share my example in case someone else finds this post useful in the future.

8

That sounds great, thanks

9

Equation

Ok here I’ll show how to solve Poisson’s equation on the unit square using only a unit interval mesh. The equation is

\begin{align} -\Delta u &= f(x,y)&&\text{on }\Omega=(0,1)^2 \\ u &= 0 &&\text{on }\partial\Omega \end{align}

where f(x,y) = 2x(y-1)(xy-2x+y+2)\exp(x-y). This has the analytic solution

u(x,y) = xy(x-1)(1-y)\exp(x-y)

which we can compare to the final numerical solution.

Analytic calculations

The crux of this approach is, given a basis \{\phi_i\} for functions on \Omega we can construct a basis \{\phi_i\phi_j\} for functions on \Omega\times\Omega.

Product basis

Let the unit interval mesh have n nodes such that x_i = \frac{i}{n-1} and let \{\phi_i\} be degree 1 CG finite elements. The idea here is to use vector-valued functions \tilde u\colon [0,1] \to \mathbb{R}^n to represent the n ordinate nodal values of u on a structured unit square mesh \{(x_i, x_j)\}, i.e. \tilde u(x)_j \approx u(x, x_j). Now let U_{ij} = u(x_i, x_j) and expand \tilde u(x)_j in our basis elements to obtain

\tilde u(x)_j \approx \sum_i U_{ij}\phi_i(x)

and hence

u(x, y) \approx \sum_{j} \tilde u(x)_j \phi_j(y) \approx \sum_{ij} U_{ij}\phi_i(x)\phi_j(y).

This way we can use \{\phi_i\phi_j\} as a finite element basis.

Variational formulation

We quickly compute a(u,v) = L(v) as

\begin{align} a(u,v) &= \int_\Omega\! \nabla u \cdot \nabla v \, dxdy \\ L(v) &= \int_\Omega\! fv\, dxdy \end{align}

Written in terms of the product basis elements the stiffness tensor A_{ij, k\ell} := a(\phi_k\phi_\ell, \hat\phi_i\hat\phi_j) where (after applying Fubini’s theorem)

a(\phi_k\phi_\ell, \hat\phi_i\hat\phi_j) = \int_0^1\! \phi'_k \hat\phi'_i \,dx \int_0^1\! \phi_\ell\hat\phi_j \,dy + \int_0^1\! \phi_k\hat\phi_i\,dx \int_0^1\! \phi'_\ell\hat\phi'_j \,dy

This suggests defining matrices

B_{ij} := \int_0^1\! \phi'_j\hat\phi'_i \,dx \quad\text{and}\quad C_{ij} := \int_0^1\! \phi_j\hat\phi_i\,dx

so we can write A_{ij,k\ell} = B_{ik}C_{j\ell} + C_{ik}B_{j\ell}, or equivalently A = B\otimes C + C\otimes B. Now let’s exploit the fact that f(x,y) = \sum_k X_k(x)Y_k(x) and compute

b_{ij} := L(\hat\phi_i\hat\phi_j ) = \int_\Omega\! f \hat\phi_i \hat\phi_j \,dxdy = \sum_k\int_0^1\! X_k \hat\phi_i\,dx \int_0^1\! Y_k \hat\phi_j\,dy

Hence we need only solve AU = b to obtain the desired solution.

Fenics implementation

We compute the integrals defining A and b using a subspace W of the vector function space V.

import numpy as np
from fenics import *

n = 11 # number of mesh nodes
mesh = UnitIntervalMesh(n-1)
V = VectorFunctionSpace(mesh, 'CG', 1, n)
v2d = vertex_to_dof_map(V)
W = V.sub(0).collapse()

Next let’s compute b as a PETSc Vector using integral forms defined on W. Explicitly, f(x,y) = \sum_{k=0}^3 X_k(x)Y_k(y) where

\begin{align} X_0(x) &= 2x \exp(x) & Y_0(x) &= y(y-1) \exp(-y) \\ X_1(x) &= -4x^2 \exp(x) & Y_1(x) &= (y-1) \exp(-y) \\ X_2(x) &= 2x^2 \exp(x) & Y_2(x) &= y(y-1) \exp(-y) \\ X_3(x) &= 4x \exp(x) & Y_3(x) &= (y-1) \exp(-y) \end{align} 
X = ['2 * x[0] * exp(x[0])', 
     '-4 * x[0] * x[0] * exp(x[0])',
     '2 * x[0] * x[0] * exp(x[0])', 
     '4 * x[0] * exp(x[0])']
X = [Expression(x, element=W.ufl_element()) for x in X]

Y = ['x[0] * (x[0] - 1) * exp(-x[0])', 
     '(x[0] - 1) * exp(-x[0])',
     'x[0] * (x[0] - 1) * exp(-x[0])', 
     '(x[0] - 1) * exp(-x[0])']
Y = [Expression(y, element=W.ufl_element()) for y in Y]

u, v = TrialFunction(W), TestFunction(W)
X_forms = [x * v * dx for x in X]
Y_forms = [y * v * dx for y in Y]

X = [assemble(xf) for xf in X_forms]
Y = [assemble(yf) for yf in Y_forms]
    
b = Vector(mesh.mpi_comm(), n**2)
b.set_local(np.kron(X[0], Y[0])[v2d])
for x, y in zip(*[X[1:], Y[1:]]):
    b.add_local(np.kron(xi, yi)[v2d])

Now we’ll assemble A using forms from W and V

B_form = u.dx(0) * v.dx(0) * dx
C_form = u * v * dx
BW = as_matrix(assemble(B_form).array())
CW = as_matrix(assemble(C_form).array())

u, v = TrialFunction(V), TestFunction(V)
BB = outer(u.dx(0), v.dx(0))
CC = outer(u, v)
a = inner(BB, C) * dx + inner(CC, B) * dx
A = assemble(a)

To handle the vanishing boundary conditions, note we need to define the whole space as the boundary for the first and last components of \tilde u.

class VecBoundary(SubDomain):
    def __init__(self, j, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.j = j

    def inside(self, x, on_boundary):
        if self.j==0 or self.j==n-1:
            on_bdy = x[0]==x[0]
        else:
            on_bdy = on_boundary
        return on_bdy

bcs = [DirichletBC(V.sub(j), 0, VecBoundary(j)) for j in range(n)]
[bc.apply(A, b) for bc in bcs]

Finally we can apply the boundary conditions, solve the linear system for U, and compare it to the analytic solution.

# get matrix representation of solution
u = Function(V)
U = u.vector()
solve(A, U, b)
U_mat = U[v2d].reshape(n,n)
U_mat = np.flip(U_mat, axis=0) # so plot is same as Cartesian coords

# get analytic matrix solution
analytic_u = lambda x,y: x * y * (x - 1) * (1 - y) * np.exp(x - y)
xy = UnitSquareMesh(n-1, n-1).coordinates()
u_mat = np.array([analytic_u(x[0], x[1]) for x in xy]).reshape(n,n)
u_mat = np.flip(u_mat, axis=0) 

print(np.linalg.norm(U_mat - u_mat)) # = 0.001789965687964315

This solution is actually a better estimate than using a unit square mesh with degree 1 CG elements, perhaps because f is better approximated in the product basis? You can view the matrix representations with, say, plt.imshow().