Perform elementwise math on solution vectors

Linear Algebra
1

I would like to perform elementwise addition and multiplication on solution vectors, but am not getting anywhere. I am using Fenics 2019.2.0 on Ubuntu 18. The MWE Code (for multiplication) is as follows:

import meshio
from dolfin import *
import os, sys, traceback

gamma_i = 0.2

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mvc, "VolumeRegions")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

cell = tetrahedron
ele_type = FiniteElement('N1curl', cell, 2, variant="integral") # H(curl) element for EM
V = FunctionSpace(mesh, ele_type)

cc = Expression(("cos(gamma * x[2])", "cos(gamma * x[2])", "cos(gamma * x[2])"), degree = 2, gamma = gamma_i)
ss = Expression(("sin(gamma * x[2])", "sin(gamma * x[2])", "sin(gamma * x[2])"), degree = 2, gamma = gamma_i)


s1 = interpolate(ss, V)
c1 = interpolate(cc, V)

ut = Function(V)
ut.vector()[:] = s1.vector()[:] * c1.vector()[:]

fp = File("ProductSum.pvd")
fp << ut

fp = File("Cos.pvd")
fp << c1

fp = File("Sin.pvd")
fp << s1

sys.exit(0)

Mesh File here.

When I print out the files for the Sin and Cos functions, I get the expected pattern.

Cos
Cos1094×803 145 KB

When I plot the product, I get rubbish -

Product
Product1094×803 161 KB

The question is, what am I doing incorrectly in the multiplication?

Cheers!

2

I haven’t used that syntax for vector pointwise multiplication before, although, in simple examples, it seems to be equivalent to directly invoking PETSc’s VecPointwiseDivide. Maybe the root of the problem is in the interpretation of products of Nédélec degrees of freedom, which is unclear to me.

3

Thank you very much @kamensky . You are indeed correct that the problem lies in the interpretation of the Nédélec degrees of freedom. Projecting the Nédélec space onto a Lagrange space and carrying out the component-wise multiplication there fixed the problem.

import meshio
from dolfin import *
import os, sys, traceback

gamma_i = 0.2
w = 22.86

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mvc, "VolumeRegions")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

cell = tetrahedron
ele_type = FiniteElement('N1curl', cell, 2, variant="integral") # H(curl) element for EM
V = FunctionSpace(mesh, ele_type)
VL = VectorFunctionSpace(mesh, "CG", degree=2, dim=3)

cc = Expression(("cos(gamma * x[2])", "cos(gamma * x[2])", "cos(gamma * x[2])"), degree = 2, gamma = gamma_i)
ss = Expression(("0.0", "sin(3.1415926 * x[0] / w)", "0.0"), degree = 2, w = w)

s1 = interpolate(ss, V)
c1 = interpolate(cc, VL)

ut = Function(VL)
ssNew = project(s1, VL)
ut.vector()[:] = ssNew.vector()[:] * c1.vector()[:]

fp = File("ProductSum.pvd")
fp << ut

fp = File("Cos.pvd")
fp << c1

fp = File("Sin.pvd")
fp << s1

sys.exit(0)

Thanks for the tip!
Cheers