Use of numpy functions with Fenics?

variational formulation
1

Hi everyone,

I am currently trying to use trignometric functions in order to calculate the value of \arccos(n[i]*[0,0,1]) across my mesh, where n[i] is defined as:

i = Index()
n = FacetNormal(mesh)

Once this has been calculated, I would then like to use this to calculate another function R across my entire mesh and then feed this value into an Expression().
However, when I try to run the code: I get this error:

Traceback (most recent call last):
  File "reflectivity.py", line 20, in <module>
    a_ = np.arccos(n[i]*zDir[i])
TypeError: loop of ufunc does not support argument 0 of type IndexSum which has no callable arccos method

I seem to remember that there is an issue using numpy with Fenics, how do I resolve this?
My full code is shown below:

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from ufl import as_tensor
from ufl import Index
import math
from mpi4py import MPI

mesh = UnitSquareMesh(16,16)
Space = FunctionSpace(mesh,'P', 1)
T = project(Expression("x[0]+0.1*x[1]", degree=2), Space)
n1 = 1.0 #refractive index of vacuum
n2 = 2.9421 # refractive index of W at 515nm
i = Index()
n = FacetNormal(mesh)
zDir = as_tensor([0,0,1]) #incoming direction of light (unit vector)
a_ = np.arccos(n[i]*zDir[i])
angle = project(a_)

def R(n1,n2,angle_i):
	angle_t = np.arcsin((n1/n2) * np.sin(angle_i))
	rs = ((n1*np.cos(angle_i) - n2*np.cos(angle_t))/(n1*np.cos(angle_i) + n2*np.cos(angle_t)))**2
	rp = ((n1*np.cos(angle_t) - n2*np.cos(angle_i))/(n1*np.cos(angle_t) + n2*np.cos(angle_i)))**2
	return 1/2 * (rs + rp)

R = R(n1,n2,angle)
L = Expression('(1-R)',degree=2, R=R)
2

Simply use ufl directly:

from fenics import *
import numpy as np
from ufl import as_tensor, Index

mesh = UnitCubeMesh(16,16, 16)
Space = FunctionSpace(mesh,'P', 1)
T = project(Expression("x[0]+0.1*x[1]", degree=2), Space)
n1 = 1.0 #refractive index of vacuum
n2 = 2.9421 # refractive index of W at 515nm
i = Index()
n = FacetNormal(mesh)
zDir = as_tensor([0,0,1]) #incoming direction of light (unit vector)
a = acos(n[i]*zDir[i])

def R(n1,n2,angle_i):
	angle_t = asin((n1/n2) * sin(angle_i))
	rs = ((n1*cos(angle_i) - n2*cos(angle_t))/(n1*cos(angle_i) + n2*cos(angle_t)))**2
	rp = ((n1*cos(angle_t) - n2*cos(angle_i))/(n1*cos(angle_t) + n2*cos(angle_i)))**2
	return 1/2 * (rs + rp)

L = 1-R(n1, n2, a)
print(assemble(L*ds))
3

Ok that makes sense, however, I do not seem to be able to use the Expession() function anymore, which is what I wanted to achieve:

L = Expression('1-R', R = R(n1, n2, a))

Gives the error:

Traceback (most recent call last):
  File "reflectivity.py", line 21, in <module>
    L = Expression('1-R', R = R(n1, n2, a))
  File "/usr/local/fkf/Fenics_openmpi2/Python3/lib64/python3.6/site-packages/dolfin/function/expression.py", line 400, in __init__
    self._cpp_object = jit.compile_expression(cpp_code, params)
  File "/usr/local/fkf/Fenics_openmpi2/Python3/lib64/python3.6/site-packages/dolfin/function/jit.py", line 158, in compile_expression
    expression = compile_class(cpp_data, mpi_comm=mpi_comm)
  File "/usr/local/fkf/Fenics_openmpi2/Python3/lib64/python3.6/site-packages/dolfin/jit/jit.py", line 170, in compile_class
    raise RuntimeError("Unable to compile C++ code with dijitso")
RuntimeError: Unable to compile C++ code with dijitso
4

Why do you want it to be a dolfin.Expression? What I have showed you is how to express the problem using pure ufl syntax, which avoids the interpolation of an expression for each assembly.
As you can observe, this works perfectly well to use in assembly.

5

In my more complex 3D example, I wish to use Expression() in order represent a Gaussian heat source illuminating one face of a 3D mesh, then use this expression inside my weak form in order to account for this heat flux on one face. My desired weak form looks like this:

F = (rho*c/dtime*(T-T0)*V - q[i]*V.dx(i)) * dv + epsilon*sigma*(T**4 - T_a4)*V*da - L*V*da(3)

Where T is the function I want to find across my whole mesh and L is the desired expression I have posted previously. Given that, from my understanding, assemble acts like a summation, I wish to have an expression that has a different value across each node on my mesh, depending on the value of a

6

So I still do not understand why you cannot just use L (as the ufl expression I showed above) directly in your variational form?
i.e.

L = 1-R(n1, n2, a)
F = (rho*c/dtime*(T-T0)*V - q[i]*V.dx(i)) * dv + epsilon*sigma*(T**4 - T_a4)*V*da - L*V*da(3)
7

Ill give it a try on my current compex code, but would L and Expression(L) do functionally the same thing, except one uses ufl and the other uses dolfin?

8

So they do mostly the same thing. One of the core differences is that an expression is interpolated into an appropriate function space (specified by degree or element as kwarg),
while the ufl expression is evaluated at quadrature points. In some cases, the quadrature degree will increase when using a ufl expression (to try to capture highly fluctuating functions).
However, you can fix the quadrature degree in the integration measures, i.e. dx(degree=5)

9

So swapping the Expression term for the pure ufl term has lead to some errors popping up since my full Expression term is defined as a Gaussian:

Laser = Expression('(1-R)*exp(-pow((x[0] - 0), 2)/(w*w)-pow((x[1]-0), 2)/(w*w))', R=R, w = 0.5)

And then switching it to the ufl expression

w = 0.5
Laser = (1-R)*exp(-pow((x[0] - 0), 2)/(w*w)-pow((x[1]-0), 2)/(w*w))

But there is an issue here as x is not defined.

10

Set
x=SpatialCoordinate(mesh)