# 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 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.

Thanks!

# My setup

I’m running fenics from the quay.io/fenicsproject/stable:current 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>

/usr/local/lib/python3.6/dist-packages/dolfin/fem/interpolation.py 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/function.py in interpolate(self, u)
363             self._cpp_object.interpolate(u._cpp_object)
364         else:
--> 365             self._cpp_object.interpolate(u)
366
367     def compute_vertex_values(self, mesh=None):

/usr/local/lib/python3.6/dist-packages/dolfin/function/expression.py in wrapped_eval(self, values, x)
51         # Wrap eval functions
52         def wrapped_eval(self, values, x):
---> 53             self.user_expression.eval(values, x)
54
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
29

<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/function.py in __call__(self, *args, **kwargs)
335
336         # The actual evaluation
--> 337         self._cpp_object.eval(values, x)
338
339         # If scalar return statement, return scalar value.

RuntimeError:

*** -------------------------------------------------------------------------
*** 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],
mesh_size=mesh_size
)
simple_mesh = geom.generate_mesh()

lines, triangles, _  = simple_mesh.cells
mesh = Mesh()
editor = MeshEditor()
editor.open(mesh, 'triangle', 2, 2)
editor.init_vertices(len(simple_mesh.points))
editor.init_cells(len(triangles.data))
for i, point in enumerate(simple_mesh.points):
for i, cell in enumerate(triangles.data):
editor.close()

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

# 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()
editor.open(curve_mesh, "interval", tdim, gdim)
n = mesh_points.shape[0] - 1
editor.init_vertices(n + 1)
editor.init_cells(n)

for vi, v in enumerate(mesh_points):

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)

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

# In[10]:

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

Parameters
----------

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):
try:
result = self._function(x[0], self._surface(x[0]))
except Exception as e:
print(x)
raise e
value[:] = result

def value_shape(self):
return (2,)

# In[11]: