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:
- Generate a mesh over 0≤ z ≤\zeta(x) from the arrays
x, zetausingpygmshand theMeshEditor. - Then I solve the Laplace problem and project the gradient of the computed solution.
- Then I interpolate this gradient (
grad_phiin the code) at a mesh containing only the points at the curve (x, \zeta(x)), built using theMeshEditor. - Finally I define a 1D Mesh and use a
fenics.UserExpressionto 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>
----> 1 grad_phi_surface = fenics.interpolate(grad_phi_expression, W_vector_space)
/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
***
*** fenics-support@googlegroups.com
***
*** 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:
geom.add_polygon(
[(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):
editor.add_vertex(i, point[:-1])
for i, cell in enumerate(triangles.data):
editor.add_cell(i, cell)
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)
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()
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):
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.set_allow_extrapolation(True)
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))
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,)
grad_phi_curve.set_allow_extrapolation(True)
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[ ]: