Integration on predefined Grid

Question

I have given a vector field $$ \vec{u}(x, y) $$ of which I know the values at a non-uniformly distributed set of points. To solve a problem, I need to evaluate a cost function

$$ J = \int_\Omega \sum\limits_{j = 1}^{N} \left|\nabla \vec{u}_j(s)\right|^2 \text{d}s $$

Where $N = 2$ and $\nabla \vec{u}_j(s)$ is the j-th component of $\nabla \vec{u}$. The paper I'm reading suggests using a Delaunay triangulation of the point set and first order polynomial element functions.

As I'm very new to firstly the FEM topic itself, and also to Fenics, I was wondering, how I could formulate my problem in python. I found this thread about numerical integration, however I can't find out how I can enter my known dataset as a expression. I guess I could generate my own Delaunay mesh of my point set using the mesh editor as proposed here. But that's not explaining how to define the function values either.

Basically I'm having something like this in my mind:

mesh = generate_delaunay_mesh_with_editor(my_points)
f = ? # how to enter my cost function?
print assemble(f * ds)

Can you hint me on some part of the documentation to read up? Or propose a way of integrating an arbitrary vector field's gradient?

Comments

Hi, what is $\Omega$ here? Is it the Delaunay triangulation?

---

Yes, it is indeed the triangulation.

Answer 1

Hi, consider

from matplotlib.tri import Triangulation, LinearTriInterpolator
from dolfin import *
import numpy as np


npoints = 20
# 'Measurement' points
x, y = np.tile(np.linspace(0, 5, 10), 2), np.repeat(np.linspace(0, 5, 10), 2)
# Delaunay triangulation for the points (considery meshpy, mshr?)
triang = Triangulation(x, y)
# Turn it into FEniCS mesh
mesh = Mesh()
mesh_vertices = np.c_[x, y]
mesh_cells = triang.triangles

editor = MeshEditor()
editor.open(mesh, 2, 2)
editor.init_vertices(len(mesh_vertices))
editor.init_cells(len(mesh_cells))

for vertex_index, v in enumerate(mesh_vertices): editor.add_vertex(vertex_index, v)

for cell_index, c in enumerate(mesh_cells): editor.add_cell(cell_index, *c)

editor.close()

# Now lets have data at the points
ux, uy = x**2+y**2, x*y

# First option to represent the data: Expression with callback to matplotlib
class TriInterp(Expression):
    def __init__(self, triang, ux, uy, **kwargs):
        self.Ix = LinearTriInterpolator(triang, ux)
        self.Iy = LinearTriInterpolator(triang, uy)

    def eval(self, value, x): 
        value[0] = self.Ix(*x)
        value[1] = self.Iy(*x)

    def value_shape(self): return (2, )

f = TriInterp(triang, ux, uy, degree=1, cell=mesh.ufl_cell())
L = inner(grad(f), grad(f))*dx(domain=mesh)
print assemble(L) 

# For later use it will be faster to represent f as function
V = VectorFunctionSpace(mesh, 'CG', 1)
v = interpolate(f, V)

Comments

This helps me a great deal, thank you very much! It's funny, that the FEniCS notation is actually quite easy to grasp, without really reading too much of the documentation.