How to interpolate at a curve (of a boundary) with a higher-degree-polynomial Lagrange function space?

Given I have a solution to a Laplace problem over a mesh that is generated dynamically, I want to extract the gradient of the solution (and maybe other functions too) at a curve that in this case coincides with the boundary of the domain. This curve is defined as a function of one parameter only and I know that it can be written as (x, \zeta(x)), with \zeta a scalar function. I want to be able to also solve the Laplace equation with Lagrange elements of arbitrary degree, but can’t figure out why it does not work :sweat_smile: for degree > 1.

Until now I have this implementation where I basically follow these steps:

  1. Generate a mesh over 0≤ z ≤\zeta(x) from the arrays x, zeta using pygmsh and the MeshEditor.
  2. Then I solve the Laplace problem and project the gradient of the computed solution.
  3. Then I interpolate this gradient (grad_phi in the code) at a mesh containing only the points at the curve (x, \zeta(x)), built using the MeshEditor.
  4. Finally I define a 1D Mesh and use a fenics.UserExpression to convert this function from a two-parameter function (\mathbb{R}^2 \to \mathbb{R}) into a 1-parameter function ([a,b] \to \mathbb{R}), i.e., a function of the x parameter only.

For degree = 1 everything works well, but when I try degree = 2 I get this error:

Mesh entity index -1 out of range [0, 99] for entity of dimension 1..

Any ideas on how to fix/improve it?

Below is a minimal working example (or the “most minimal” I could think of) that I adapted from the poisson example and my “real” code.


My setup

I’m running fenics from the docker image on a docker installation on macOS Big Sur using JupyterLab.

The trace

RuntimeError                              Traceback (most recent call last)
<ipython-input-11-d61373bff4a6> in <module>
----> 1 grad_phi_surface = fenics.interpolate(grad_phi_expression, W_vector_space)

/usr/local/lib/python3.6/dist-packages/dolfin/fem/ in interpolate(v, V)
     69     # Compute interpolation
     70     if hasattr(v, "_cpp_object"):
---> 71         Pv.interpolate(v._cpp_object)
     72     else:
     73         Pv.interpolate(v)

/usr/local/lib/python3.6/dist-packages/dolfin/function/ in interpolate(self, u)
    363             self._cpp_object.interpolate(u._cpp_object)
    364         else:
--> 365             self._cpp_object.interpolate(u)
    367     def compute_vertex_values(self, mesh=None):

/usr/local/lib/python3.6/dist-packages/dolfin/function/ in wrapped_eval(self, values, x)
     51         # Wrap eval functions
     52         def wrapped_eval(self, values, x):
---> 53             self.user_expression.eval(values, x)
     55         def wrapped_eval_cell(self, values, x, cell):

<ipython-input-10-f26abbf87b9c> in eval(self, value, x)
     25         except Exception as e:
     26             print(x)
---> 27             raise e
     28         value[:] = result

<ipython-input-10-f26abbf87b9c> in eval(self, value, x)
     22     def eval(self, value, x):
     23         try:
---> 24             result = self._function(x[0], self._surface(x[0]))
     25         except Exception as e:
     26             print(x)

/usr/local/lib/python3.6/dist-packages/dolfin/function/ in __call__(self, *args, **kwargs)
    336         # The actual evaluation
--> 337         self._cpp_object.eval(values, x)
    339         # If scalar return statement, return scalar value.


*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error:   Unable to create mesh entity.
*** Reason:  Mesh entity index -1 out of range [0, 99] for entity of dimension 1.
*** Where:   This error was encountered inside MeshEntity.cpp.
*** Process: 0
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

The code

This is a jupyter nbconvert of the original notebook.

#!/usr/bin/env python
# coding: utf-8

# In[1]:

import numpy as np
import fenics
import matplotlib.pyplot as plt
import ZCS

import scipy.interpolate
import pygmsh

from fenics import *

# In[2]:

polynomial_degree = 2

x = np.linspace(0, 10, 100)

zeta = np.cos(x*np.pi/3)+10

# Build the mesh

# In[3]:

mesh_size = 1
with pygmsh.geo.Geometry() as geom:
        [(xi,0) for xi in x] 
         + [(xi, zi) for (xi,zi) in zip(x, zeta)][::-1],
    simple_mesh = geom.generate_mesh()

lines, triangles, _  = simple_mesh.cells
mesh = Mesh()
editor = MeshEditor(), 'triangle', 2, 2)
for i, point in enumerate(simple_mesh.points):
    editor.add_vertex(i, point[:-1])
for i, cell in enumerate(
    editor.add_cell(i, cell)

# Solve the Laplace problem

# In[4]:

V = FunctionSpace(mesh, 'P', polynomial_degree)

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

def surface_boundary(x, on_boundary):
    return on_boundary and x[1]>=1

bc = DirichletBC(V, u_D, surface_boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Project `grad_phi`

# In[5]:

vector_element = fenics.MixedElement([V.ufl_element(), V.ufl_element()])
V_vec = fenics.FunctionSpace(V.mesh(), vector_element)
grad_phi = project(grad(u), V_vec)

# Build a mesh over the curve of the top boundary

# In[6]:

mesh_points = np.array([[xi, zetai] for (xi, zetai) in zip(x, zeta)])

tdim = 1
gdim = 2
curve_mesh = Mesh()
editor = MeshEditor(), "interval", tdim, gdim)
n = mesh_points.shape[0] - 1
editor.init_vertices(n + 1)

for vi, v in enumerate(mesh_points):
    editor.add_vertex(vi, fenics.Point(v))

for ci in range(n):
    editor.add_cell(ci, np.array([ci, ci + 1], dtype="uintp"))

# Interpolate grad_phi over that boundary

# In[7]:

V_curve = fenics.FunctionSpace(curve_mesh, "P", polynomial_degree)
V_curve_vector_element = fenics.MixedElement([V_curve.ufl_element(), V_curve.ufl_element()])
V_curve_vector_space = fenics.FunctionSpace(V_curve.mesh(), V_curve_vector_element)


grad_phi_curve = fenics.interpolate(grad_phi, V_curve_vector_space)

# Build the target 1D mesh and function space

# In[8]:

x_mesh = fenics.IntervalMesh(x.shape[0] - 1, x.min(), x.max())

W = fenics.FunctionSpace(x_mesh, "P", polynomial_degree)
W_vector_element = fenics.MixedElement([W.ufl_element(), W.ufl_element()])
W_vector_space = fenics.FunctionSpace(W.mesh(), W_vector_element)

# Make a degree-1 spline for the surface points that should match the straight lines of the mesh

# In[9]:

surface_lines = scipy.interpolate.interp1d(
    grad_phi_curve.function_space().mesh().coordinates()[:, 0],
    grad_phi_curve.function_space().mesh().coordinates()[:, 1],

# In[10]:

class mapping_2D_to_1D(fenics.UserExpression):
    fenics.UserExpression that represents the mapping x -> function(x, zeta(x))


    function: function
        Function "f" in the formula above,
        over R2 that can be evaluated as f(x,z)
    surface: function
        Scalar function "zeta" in the formula above,
        over R1 that can be evaluated as zeta(x)


    def __init__(self, function, surface, *arg, **kwargs):
        super().__init__(*arg, **kwargs)
        self._function = function
        self._surface = surface

    def eval(self, value, x):
            result = self._function(x[0], self._surface(x[0]))
        except Exception as e:
            raise e
        value[:] = result

    def value_shape(self):
        return (2,)

grad_phi_expression = mapping_2D_to_1D(grad_phi_curve, surface_lines)

# In[11]:

grad_phi_surface = fenics.interpolate(grad_phi_expression, W_vector_space)

# In[ ]: