Dot product of tensor and vector

Linear Algebra
1

I’m trying to evaluate the dot product of a tensor with the surface normal on a curved boundary.

import dolfin as df
import mshr as ms
import matplotlib.pyplot as plt

circle = ms.Circle(df.Point(0,0), 1.0, 18)
mesh = ms.generate_mesh(circle, 18)
bmesh = df.BoundaryMesh(mesh, 'exterior')

TFS = df.TensorFunctionSpace(bmesh, 'P', 1)
sigma = df.Function(TFS)
sigma.interpolate(df.Expression(
    (('x[0]','x[0]*x[1]'),('x[0]*x[1]','x[1]')), 
    degree=1))

VFS = df.VectorFunctionSpace(bmesh, 'P', 1)


normal = df.interpolate(
        df.Expression(
            ('x[0]/sqrt(x[0]*x[0]+x[1]*x[1])', 
                'x[1]/sqrt(x[0]*x[0]+x[1]*x[1])'), 
        degree=1), VFS)
ex = df.interpolate(df.Constant((1,0)), VFS)
ey = df.interpolate(df.Constant((0,1)), VFS)

force = df.project(df.dot(sigma, normal), VFS)
fx = df.project(df.dot(sigma, ex), VFS)
fy = df.project(df.dot(sigma, ey), VFS)

force_exact = df.interpolate(
        df.Expression(
            ('x[0]*(x[0]+x[1]*x[1])/sqrt(x[0]*x[0]+x[1]*x[1])', 
            'x[1]*(x[0]*x[0]+x[1])/sqrt(x[0]*x[0]+x[1]*x[1])'),
        degree=1), VFS)

fx_exact = df.interpolate(
        df.Expression(('x[0]', 'x[0]*x[1]'), degree=1), VFS)
fy_exact = df.interpolate(
        df.Expression(('x[0]*x[1]', 'x[1]'), degree=1), VFS)

err = force.compute_vertex_values() - force_exact.compute_vertex_values()
err_x = fx.compute_vertex_values() - fx_exact.compute_vertex_values()
err_y = fy.compute_vertex_values() - fy_exact.compute_vertex_values()

fig, (ax1, ax2, ax3) = plt.subplots(1,3, figsize=(12,5))
ax1.hist(err)
ax2.hist(err_x)
ax3.hist(err_y)

plt.show()

This produces the figure

Figure
Figure.png1200×500 7.57 KB

The errors when evaluating the dot product of the tensor with a constant vector are at machine precision while the errors when the tensor is dotted with a non-constant vector (in this case the surface normal of a curved boundary) are way large.

Why is this? Am I doing a mistake or is this the best that one can obtain?

2

Did you already try to increase the resolution of you solution, i.e., try to use a finer mesh or a higher degree for the finite elements (especially for some boundary integration with the normal only a finer mesh will help).
I think that this should already help in getting lower errors.

3

Thanks for the comment.

Yes, the errors go down with mesh refinement…but very slowly.

On the other hand, if I compute the components of the (varying) normal vector at the same mesh resolution, I get near machine precision errors.

import dolfin as df
import mshr as ms
import matplotlib.pyplot as plt

circle = ms.Circle(df.Point(0,0), 1.0, 18)
mesh = ms.generate_mesh(circle, 18)
bmesh = df.BoundaryMesh(mesh, 'exterior')

TFS = df.TensorFunctionSpace(bmesh, 'P', 1)
sigma = df.Function(TFS)
sigma.interpolate(df.Expression(
    (('x[0]','x[0]*x[1]'),('x[0]*x[1]','x[1]')), 
    degree=1))

VFS = df.VectorFunctionSpace(bmesh, 'P', 1)
SFS = df.FunctionSpace(bmesh, 'P', 1)


normal = df.interpolate(
        df.Expression(
            ('x[0]/sqrt(x[0]*x[0]+x[1]*x[1])', 
                'x[1]/sqrt(x[0]*x[0]+x[1]*x[1])'), 
        degree=1), VFS)
ex = df.interpolate(df.Constant((1,0)), VFS)
ey = df.interpolate(df.Constant((0,1)), VFS)

nx = df.project(df.dot(normal, ex), SFS)
ny = df.project(df.dot(normal, ey), SFS)

nx_exact = df.interpolate(
            df.Expression('x[0]/sqrt(x[0]*x[0]+x[1]*x[1])', 
        degree=1), SFS)
ny_exact = df.interpolate(
            df.Expression('x[1]/sqrt(x[0]*x[0]+x[1]*x[1])', 
        degree=1), SFS)

err_x = nx.compute_vertex_values() - nx_exact.compute_vertex_values()
err_y = ny.compute_vertex_values() - ny_exact.compute_vertex_values()

fig, (ax2, ax3) = plt.subplots(1,2, figsize=(10,5))
ax2.hist(err_x)
ax3.hist(err_y)

plt.show()

gives

Figure_1
Figure_1.png1000×500 5.63 KB

Why is the dot product between a tensor and a vector not evaluated to the same precision as the dot product between a vector and another vector?

Perhaps I’m missing something very simple and silly.

PS: I have also checked having ex and ey as Expressions rather than Constants with the same results.