How to integrate a single basis function over a given domain?

1

Hello,

I am searching for a function or a way in FEniCS to compute the norm:

||\psi_{E}||_{L^{2}(\Omega)}^{2} = \int_{\Omega} \psi_{E}^{2} \, dx,

where \psi_{E} is a Nedelec basis function on an edge of choice. I need this information for a sharp error estimate.

To compute a norm of a solution of a Maxwell equation, I would solve the equation and then compute a norm of the solution, for example like this:

from fenics import *
import numpy as np

mesh = UnitCubeMesh(2,2,2)
NE = FiniteElement(‘N1curl’,tetrahedron,1)
V = FunctionSpace(mesh, NE)

f = Constant((1,1,1))

u = TrialFunction(V)
v = TestFunction(V)

a = assemble(inner(curl(u), curl(v))*dx)
b = assemble(inner(f,v)*dx)

a = a.array()
b=b.get_local()
x=np.linalg.lstsq(a,b,rcond=None)[0]
c.vector().set_local(x)
norm_x = assemble(inner(c,c)*dx)

Now, I would like to compute a norm of a basis Nedelec function \psi_{E} in H(curl). My idea would give the norm of x = \sum_{E} c_{E} \psi_{E}, where c_{E} are constants. However, I would like to extract a single basis function \psi_E and compute its norm. To my knowledge, I don’t need to solve the Maxwell equation to do that.

I have found how to evaluate basis function at a given point here: https://fenicsproject.org/qa/8627/how-to-evaluate-basis-function-at-points-of-cell/

However, I have found no solution how to integrate a single basis function.

Does someone have an idea?

2

Since I am not sure how matrices in FEniCS are built (especially when a sparse structure is used), I am puzzled whether my idea could work for my purpose:

mesh_d = UnitCubeMesh(1,1,1)
mesh_d.init()
print(mesh_d.num_edges())
NE = FiniteElement(‘N1curl’,tetrahedron,1)
V = FunctionSpace(mesh_d, NE)
w1 = TrialFunction(V)
tau1 = TestFunction(V)

a = assemble(inner(w1,tau1)*dx).array()
#print(a)
print(a.shape)

for E in edges(mesh_d):
print(a[E.index(),E.index()])

Are the entries on the diagonal of the matrix already the norms
$$||\psi_{E}||{L^{2}(\Omega)}^{2} = \int{\Omega} \psi_{E}^{2} , dx?$$

Since the mesh has 19 edges in my example and since the matrix has the shape (19,19), it seems to me, that the norms for each \psi_E are simply placed on the diagonal of the matrix. Is someone familiar with the output of ‘assemble’ and the underlying structure of the matrix? Do the indices in the matrix ‘assemble(inner(w1,tau1)*dx).array()’ coincide with E.index[] of the mesh or is there a special mapping in between?

3

Since there are so many views on my question, I decided to share my code and my solution.

The idea ist based on the github code:

and the paper

Computational homogenization of time-harmonic Maxwell’s equations
Patrick Henning and Anna Persson

Here is the piece of code for computation in 2 D:

from fenics import *
import numpy as np
from dof2edge2 import *

mesh = UnitSquareMesh(1,1)

print(“Lagrange basis functions”)

LA = FiniteElement(‘Lagrange’,triangle,1)
Q = FunctionSpace(mesh,LA)
w = TrialFunction(Q)
tau = TestFunction(Q)

L = assemble(w * tau * dx).array()
v_d=vertex_to_dof_map(Q)

psi=np.zeros(mesh.num_vertices())
for vertt in vertices(mesh):
psi[vertt.index()] = L[v_d[vertt.index()],v_d[vertt.index()]]

print("||phi_0||_omega^{2}",psi)

print(“manual solution”)

phi_basic = Function(Q)

e_vector = np.zeros(mesh.num_vertices())

for k in range(mesh.num_vertices()):
e_vector[v_d[k]]=1
phi_basic.vector().set_local(e_vector)
e_vector[v_d[k]]=0
sol = assemble(phi_basic * phi_basic * dx)
solely = assemble(phi_basic*dx) # integral over a basis function
print("||phi_0||_omega^{2}",sol)
print("\int_Omega phi_0 dx",solely)

print(“Nedelec basis functions”)

NE = FiniteElement(‘N1curl’,triangle,1)
Q1 = FunctionSpace(mesh,NE)
w1 = TrialFunction(Q1)
tau1 = TestFunction(Q1)

L = assemble(inner(w1,tau1)*dx).array()

d2e, e2d = dof2edge(mesh)
psi1=np.zeros(mesh.num_edges())
for E in edges(mesh):
psi1[E.index()] = L[e2d[E.index()],e2d[E.index()]]

print("||phi_0||_omega^{2}",psi1)

print(“manual solution”)

phi_basic1 = Function(Q1)

e_vector = np.zeros(mesh.num_edges())

for k in range(mesh.num_edges()):
e_vector[e2d[k]]=1
phi_basic1.vector().set_local(e_vector)
e_vector[e2d[k]]=0
sol = assemble(inner(phi_basic1,phi_basic1)*dx)
assemble((phi_basic1[0]*Constant(1))*dx +(phi_basic1[1]*Constant(1))*dx )
print("||phi_0||_omega^{2}",sol)

The code in dof2edge2.py constructs a mapping between the numbering of the edges and degrees of freedom.
dof2edge2.py: (Code of Anna Persson)

from fenics import *
import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse as sps

print(" Computes two mappings: dof to edge and edge to dof")

def dof2edge(mesh):

V = FunctionSpace(mesh, 'N1curl', 1)
N = V.dim()

# Cell to edge map
cell2edges = V.mesh().topology()(2, 1)
# Cell to dof map
cell2dofs = V.dofmap().cell_dofs

dof2edge = np.zeros(mesh.num_edges(), dtype="int")
edge2dof = np.zeros(mesh.num_edges(), dtype="int")
# Iterate over cells, associating the edges to the dofs for that cell
for c in range(mesh.num_cells()):
    # get the global edge numbers for this cell
    c_edges = cell2edges(c)
    # get the global dof numbers for this cell
    c_dofs = cell2dofs(c)
    # associate the edge numbers to the corresponding dof numbers and vice versa
    dof2edge[c_dofs] = c_edges
    edge2dof[c_edges] = c_dofs

return dof2edge, edge2dof

I hope this helps other people to solve similar problems.