Integrate over some but not all dimensions

1

Hi,

I have managed to get fenics to solve a 2d problem for me. Mesh and Function space are defined below. As a diagnostic it would be very helpful to have the 1d quantity that is given by integrating over ONE of the two dimensions. (The solution quickly falls off towards the boundaries in that dimension and the integral is well defined.)

# 2d phasespace
phasespace = RectangleMesh(Point(0., -vmax), Point(Lx, vmax), Nx, Nv, "left/right")
# Define function space for phasespace density
V2d = FunctionSpace(phasespace, 'P', degree)

I have found 2d integration on cross section of 3d mesh which basically says that this is not possible and Integrating a 2D function along a single dimension - FEniCS Q&A from the old Q/A site that gives a solution, but when I try it the line

self.u_2d=u_2d

leads to a infinite recursion and error messages like

  [Previous line repeated 494 more times]
  File "/usr/lib/python3/dist-packages/dolfin/function/expression.py", line 405, in __getattr__
    if name == 'user_parameters':
RecursionError: maximum recursion depth exceeded in comparison

Given that the post is from 2014 I assume that internals in fenics have changed about three times since then…

Any advice how to go about this?

2

For this specific error message, see for instance: RecursionError: maximum recursion depth exceeded in the ft11_magnetostatics.py example of the "Solving PDEs in Python - The FEniCS tutorial Part 1" - #2 by alechammond

3

Thanks for that pointer. The updated example from back then now looks like

#!/usr/bin/env python3

from dolfin import *
import numpy
import matplotlib.pyplot as plt

#2D domain to define u(x,y)
m_2D=RectangleMesh(Point(0.,0.), Point(1.,1.), 10, 10)

V_2D=FunctionSpace(m_2D, "CG",1)

U_2D_e=Expression('x[0]+x[1]', degree=2)

U_2D_f=interpolate(U_2D_e,V_2D)

#1D integration domain
m_1D=IntervalMesh(20,0.,1.)
V_1D=FunctionSpace(m_1D, "CG",1)

#an Expression class to evaluate u(x,y) on the 1D domain (range y, x fixed)

class my1DExpression(UserExpression):
    def __init__(self,u_2d,x):
        super().__init__(x)
        self.u_2d=u_2d
        self.x=x
        self._vx= numpy.array([0.])
        self._pt_x= numpy.array([0.,0.])
    def eval(self, values, x):
        self._pt_x[0]= self.x
        self._pt_x[1]= x[0]
        self.u_2d.eval(self._vx,self._pt_x)
        values[0] = self._vx[0]

#An exemple int_0^1 x+y dy=x+1/2
for x in numpy.linspace(0.,1.,10):
    Uy=interpolate(my1DExpression(U_2D_f,x=x),V_1D)
    intUy=assemble(Uy*dx)
    print(intUy,x+0.5,intUy-(x+0.5))

#it is possible to define the u(x) function
class UxExpression(UserExpression):
    def __init__(self,u_2d):
        super().__init__()
        self.u_2d=u_2d
        self.x=x
    def eval(self, values, x):
        _Uy=interpolate(my1DExpression(U_2D_f,x=x[0]),V_1D)
        values[0] = assemble(_Uy*dx)

plot(U_2D_f)
plt.savefig("U_2D_f.png")
plt.close()

#1D for ux
m_1Dx=IntervalMesh(20,0.,1.)
V_1Dx=FunctionSpace(m_1Dx, "CG",1)

Ux_f=interpolate(UxExpression(U_2D_f),V_1Dx)
plot(Ux_f)
plt.savefig("Ux_f.png")

and seems to work as intended. Incorporating this into my existing code will take a while. One thing you can maybe answer in the mean time: is there a good way to combine the two helper classes?

4

Since I have added that I am getting a lot of

WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.

messages. Any advice on fixing that?

5

It’s just a warning, so you shouldn’t care about that. However, if you want to avoid the warning try this:

class my1DExpression(UserExpression):
    def __init__(self,u_2d,x):
        super().__init__(x)
        self.u_2d=u_2d
        self.x=x
        self._vx= numpy.array([0.])
        self._pt_x= numpy.array([0.,0.])
    def eval(self, values, x):
        self._pt_x[0]= self.x
        self._pt_x[1]= x[0]
        self.u_2d.eval(self._vx,self._pt_x)
        values[0] = self._vx[0]
    def value_shape(self):
        return ()
1 Like
6

yeah, that helped, thanks.

8

Is there a way to use this method to integrate a 3D function in only ONE dimension?
Any minimum example would help.
Thanks!

9

After a while I need that as well, which required me to actually understand what the 2d example does. I ended up using something like this:


#!/usr/bin/env python3

from fenics import *
import numpy as np
import matplotlib.pyplot as plt

# for testing we are going to use the function f(x,y,z) = 3/16*x*y**2 + z**3
# and integrate it over y and z to get
# F(x) = \int_{-2}^{2} \int_{0}^{1} f(x,y,z) dz dy = x + 1

# 3D mesh on which we define f(x,y,z)
m_3D = BoxMesh(Point(0.,-2.,0.), Point(5,2.,1.), 10, 10, 10)
V_3D = FunctionSpace(m_3D, "CG", 2)

# Expression for f(x,y,z)
f_exp = Expression('3.*x[0]*x[1]*x[1]/16. + x[2]*x[2]*x[2]', degree=2)
# Interpolate onto 3d space
f = interpolate(f_exp, V_3D)


#2D subspace over which we integrate
m_2D = RectangleMesh(Point(-2.,0.), Point(2.,1.), 10, 10)
V_2D = FunctionSpace(m_2D, "CG", 2)

# an Expression class to evaluate f(x,y,z) at a fixed position in x,
# with y and z from the subspace
class my1DExpression(UserExpression):
    def __init__(self,f_3d,xpos):
        super().__init__(xpos)
        self.f_3d = f_3d
        self.xpos = xpos
        self.value = np.array([0.])
        self.point = np.array([0.,0.,0.])
    def eval(self, values, x):
        # x[] is a location in the 2d subspace, as a two element numpy array
        self.point[0] = self.xpos
        self.point[1] = x[0]
        self.point[2] = x[1]
        self.f_3d.eval(self.value,self.point)
        values[0] = self.value[0]
    def value_shape(self):
        return ()

#Check that the integral matches x+1
for x in np.linspace(0.,5.,11):
    fx = interpolate(my1DExpression(f,xpos=x),V_2D)
    F = assemble(fx*dx)
    print(F,x+1.,F-(x+1.))

#it is possible to define a UserExpression for the F(x) function
class FxExpression(UserExpression):
    def __init__(self,f_3d):
        super().__init__()
        self.f_3d = f_3d
    def eval(self, values, x):
        _Fx = interpolate(my1DExpression(self.f_3d,xpos=x[0]),V_2D)
        values[0] = assemble(_Fx*dx)
    def value_shape(self):
        return ()

# 1D space where F(x) lives
m_1D = IntervalMesh(11, 0., 5.)
V_1D = FunctionSpace(m_1D, "CG", 2)

Fx = interpolate(FxExpression(f), V_1D)

plot(Fx)
plt.savefig("Fx.png")
1 Like