Integrate a scalar product

dolfinx
1

Hello dolfinx users,

I want to compute the surface average of a normalised scalar product \frac{1}{S}\int\frac{\left|\left<u,v\right>\right|}{\left||u\right||\left||v\right||}dS over a region of the mesh in dolfinx. I’m doing something wrong though, because I’m always getting something larger than one.

import ufl
import numpy as np
import dolfinx as dfx
from mpi4py.MPI import COMM_WORLD as comm

# Shortforms
n=10
mesh = dfx.mesh.create_unit_square(comm, n, n)
def save(name,f):
	with dfx.io.XDMFFile(comm, name+".xdmf", "w") as xdmf:
		xdmf.write_mesh(mesh)
		xdmf.write_function(f)

# Create vectors (placeholder)
FS_v = dfx.fem.FunctionSpace(mesh,ufl.VectorElement("CG",mesh.ufl_cell(),2))
v1, v2 = dfx.fem.Function(FS_v), dfx.fem.Function(FS_v)
v1.interpolate(lambda x: x.reshape(2,-1))
def v2_f(x):
	x[1,:]=x[0,:]-x[1,:]
	return x.reshape(2,-1)
v2.interpolate(v2_f)
save("v1",v1)
save("v2",v2)

# Integration area (triangle)
x0,x1,x2=(0,.5),(1,1),(1,0)
def line(x,x1,x2): return x[1]<(x2[1]-x1[1])*(x[0]-x1[0])/(x2[0]-x1[0])+x1[1]
def kernel(x):
	res = np.zeros_like(x[0])
	res[line(x,x0,x1)]=1
	res[line(x,x0,x2)]=0
	return res

FS_e = dfx.fem.FunctionSpace(mesh,ufl.FiniteElement("CG",mesh.ufl_cell(),2))
k = dfx.fem.Function(FS_e)
k.interpolate(kernel)
s=dfx.fem.assemble_scalar(dfx.fem.form(k*ufl.dx)).real # Scaling factor for weighting surface
if s==0: s=1 # If the kernel has no support over this proc, don't do anything

def dot(u,v=None):
	if v is None: v=u
	return ufl.real(ufl.dot(u,v))
# Compute dot product and orientation
ae,re = dfx.fem.Function(FS_e),dfx.fem.Function(FS_e)
ae.interpolate(dfx.fem.Expression(ufl.sqrt(dot(v1,v2)**2)*k,FS_e.element.interpolation_points()))
ae.x.array[:]=np.positive(np.nan_to_num(ae.x.array[:])/s)
ae.x.scatter_forward()
save("pdt",ae)
re.interpolate(dfx.fem.Expression(ae/ufl.sqrt(dot(v1))/ufl.sqrt(dot(v2)),FS_e.element.interpolation_points()))
re.x.array[:]=np.nan_to_num(re.x.array[:])
re.x.scatter_forward()
save("alignement",re)

# Compute and communicate integrals
ae,re=dfx.fem.assemble_scalar(dfx.fem.form(ae*ufl.dx)).real,dfx.fem.assemble_scalar(dfx.fem.form(re*ufl.dx)).real
ae,re=comm.gather(ae),comm.gather(re)
if comm.rank==0:
	ae,re = sum(ae), sum(re)
	print("a=",ae,"re=",re,flush=True)

Even though everything here is real, my use case would be simpler if I could run it in complex mode. Right now I get the wrong result re\approx7 indifferently in real of complex. My use case is also parallel (I ran this on 8 procs).

2

Found my own mistake should’ve summed up the s instead of dividing by proc. Below is a working code :

import ufl
import numpy as np
import dolfinx as dfx
from mpi4py.MPI import COMM_WORLD as comm

# Shortforms
n=10
mesh = dfx.mesh.create_unit_square(comm, n, n)
def save(name,f):
	with dfx.io.XDMFFile(comm, name+".xdmf", "w") as xdmf:
		xdmf.write_mesh(mesh)
		xdmf.write_function(f)

# Create vectors (placeholder)
FS_v = dfx.fem.FunctionSpace(mesh,ufl.VectorElement("CG",mesh.ufl_cell(),2))
v1, v2 = dfx.fem.Function(FS_v), dfx.fem.Function(FS_v)
v1.interpolate(lambda x: x.reshape(2,-1))
def v2_f(x):
	x[1,:]=x[0,:]-x[1,:]
	return x.reshape(2,-1)
v2.interpolate(v2_f)
save("v1",v1)
save("v2",v2)

# Integration area (triangle)
x0,x1,x2=(0,.5),(1,1),(1,0)
def line(x,x1,x2): return x[1]<(x2[1]-x1[1])*(x[0]-x1[0])/(x2[0]-x1[0])+x1[1]
def kernel(x):
	res = np.zeros_like(x[0])
	res[line(x,x0,x1)]=1
	res[line(x,x0,x2)]=0
	return res

FS_e = dfx.fem.FunctionSpace(mesh,ufl.FiniteElement("CG",mesh.ufl_cell(),2))
k = dfx.fem.Function(FS_e)
k.interpolate(kernel)
s=dfx.fem.assemble_scalar(dfx.fem.form(k*ufl.dx)).real # Scaling factor for weighting surface
if s==0: s=1 # If the kernel has no support over this proc, don't do anything

def dot(u,v=None):
	if v is None: v=u
	return ufl.real(ufl.dot(u,v))
# Compute dot product and orientation
ae,re = dfx.fem.Function(FS_e),dfx.fem.Function(FS_e)
ae.interpolate(dfx.fem.Expression(ufl.sqrt(dot(v1,v2)**2)*k,FS_e.element.interpolation_points()))
ae.x.array[:]=np.positive(np.nan_to_num(ae.x.array[:]))
ae.x.scatter_forward()
save("pdt",ae)
re.interpolate(dfx.fem.Expression(ae/ufl.sqrt(dot(v1))/ufl.sqrt(dot(v2)),FS_e.element.interpolation_points()))
re.x.array[:]=np.nan_to_num(re.x.array[:])
re.x.scatter_forward()
save("alignement",re)

# Compute and communicate integrals
ae,re=dfx.fem.assemble_scalar(dfx.fem.form(ae*ufl.dx)).real,dfx.fem.assemble_scalar(dfx.fem.form(re*ufl.dx)).real
ae,re,s=comm.gather(ae),comm.gather(re),comm.gather(s)
if comm.rank==0:
        s=sum(s)
	ae,re = sum(ae)/s, sum(re)/s
	print("a=",ae,"re=",re,flush=True)