Inner product of vectors is miscalculated

General
1

Hi, I am using legacy fenics and trying to initialize the vector field with unit vectors of random directions.

from dolfin import *
import numpy as np

mesh = Mesh()
mesh = BoxMesh(Point(0,0,0), Point(1,1,1), 1, 1, 1)

x = SpatialCoordinate(mesh)

vector1 = VectorElement("Lagrange", mesh.ufl_cell(), 1)
scalar1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
scalarspace = FunctionSpace(mesh, scalar1)
vectorspace = FunctionSpace(mesh, vector1)

v0 = Function(vectorspace)
num_dofs = vectorspace.dim()
dim = mesh.geometry().dim()
num_nodes = num_dofs // dim
rand_vec = np.random.rand(num_nodes, dim)
norms = np.linalg.norm(rand_vec, axis=1, keepdims=True)
rand_vecs = rand_vec / norms
dofx = vectorspace.sub(0).dofmap().dofs()
dofy = vectorspace.sub(1).dofmap().dofs()
dofz = vectorspace.sub(2).dofmap().dofs()
for i in range(num_nodes):
    v0.vector()[dofx[i]] = rand_vecs[i][0]
    v0.vector()[dofy[i]] = rand_vecs[i][1]
    v0.vector()[dofz[i]] = rand_vecs[i][2]

normv0 = project(dot(v0,v0),scalarspace)
print(normv0.vector()[:])

output_file = XDMFFile("test.xdmf")
    
output_file.parameters["flush_output"] = True
output_file.parameters["functions_share_mesh"] = True

v0.rename("v0", "v0")
normv0.rename("testscal", "testscal")
output_file.write(v0, float(0.0))
output_file.write(normv0, float(0.0))

In paraview, when I open the “test.xdmf” file, the magnitudes of v0 are 1 everywhere in the mesh. But normv0, the norm of v0, is not one. “inner” instead of “dot” also does not work. I cannot figure it out…

2

I am not sure what you expect here. Even if the functions have norm 1 at the vertices, doesn’t mean that this is the case in the entire domain.

Consider for instance:

from mpi4py import MPI

import dolfinx.fem.petsc
from dolfinx import mesh, fem
import numpy as np
import basix.ufl
import ufl

# Define mesh

metadata = {"quadrature_rule": "default"}
# metadata={"quadrature_rule": "vertex"}

nx = 1
domain = mesh.create_unit_interval(MPI.COMM_WORLD, nx)

## Function space and fields
Uel = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1, shape=(2,))
V = fem.functionspace(domain, Uel)
uh = fem.Function(V)
z = 0.2
uh.x.array[0] = z
uh.x.array[1] = np.sqrt(1 - z**2)
uh.x.array[2] = -z
uh.x.array[3] = np.sqrt(1 - z**2)

t = uh.x.array.reshape(-2, 2)
print(np.linalg.norm(t, axis=1))

V0, _ = V.sub(0).collapse()

u = ufl.TrialFunction(V0)
v = ufl.TestFunction(V0)
a = ufl.inner(u, v) * ufl.dx
L = ufl.inner(ufl.dot(uh, uh), v) * ufl.dx(metadata=metadata)

x = fem.Function(V0)
problem = fem.petsc.LinearProblem(
    a,
    L,
    u=x,
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
    },
)
problem.solve()

print(x.x.array)
with dolfinx.io.VTXWriter(domain.comm, "u.bp", [uh]) as bp:
    bp.write(0.0)

with dolfinx.io.VTXWriter(domain.comm, "unorm.bp", [x]) as bp:
    bp.write(0.0)

expr_norm = dolfinx.fem.Expression(ufl.dot(uh, uh), np.array([[0.5]]))
expr = dolfinx.fem.Expression(uh, np.array([[0.5]]))

print(
    expr.eval(domain, np.array([0], dtype=np.int32)),
    expr_norm.eval(domain, np.array([0], dtype=np.int32)),
)

if you use a standard quadrature rule, you evaluate at a place where the norm is not 1.
The code above gives the “right” result if you use a vertex quadrature.
Here I just illustrate this on a single interval mesh, with the output

[1. 1.]
[0.97333333 0.97333333]
[[0.        0.9797959]] [[0.96]]

where the first line is the norm at the vertices, the second line is the projected norm, and the third line is the value at the midpoint of the single cell, and its corresponding norm.

3

Hi dokken, thank you for your reply.

I think, at least, normv0 should be 1 at the vertices, but it is not 1 everywhere.

[1. 1.]
[0.97333333 0.97333333]
[[0. 0.9797959]] [[0.96]]

I understand that ‘vectors’ are interpolated, so the norms would not be 1 and that’s why the third line, ‘[[0. 0.9797959]] [[0.96]]’, is printed. But I still don’t understand why x.x.array is [0.97333333 0.97333333], not [1,1]. Since the vector fields have norm 1 at the nodes, I assumed that the projected norms would also be 1 at the nodes. Why is it [0.97333333 0.97333333] ?

Actually I wanted to assign a randomly oriented normal vector to each element during initialization. But I couldn’t find how to initialize the field by element. So I initialize the vector field at each node, which caused this problem.

4

If you in my code change the quadrature rule to be a vertex quadrature. Ie the integrals are evaluated by summing up the contributions at the the vertices

Ie uncomment this line, you will get that the projected norm will be 1.

A default quadrature evaluates the field at quadrature points internally in the cell, where the norm isn’t 1, thus the approximation is a bit off.

5

You can use a DG 0 space if you want to initialize the field per element

6

I understand now, thank you.
I am not familiar with fenicsx, how can I do that in legacy fenics?

Then does that mean it is impossible to initialize the field per quadrature points (gaussian integration points) using CG space? Is it possible in FEniCSx?

7

I think you can call
dx(metadata={'quadrature_degree': 1, 'quadrature_scheme': 'vertex'})

A CG space does in general not have dofs aligning with the quadrature points.

I think @adeebkor has some examples of collocating dofs and integration points.

1 Like