Hello,
There is a behavior that I don’t understand when taking the average of two functions, I produced the following minimal working example:
from __future__ import print_function
from fenics import *
from mshr import *
# Create mesh
mesh = generate_mesh(Rectangle(Point(0, 0), Point(1.0, 1.0)), 8)
P_omega_n = VectorElement( 'P', triangle, 3 )
P_z_n = FiniteElement( 'P', triangle, 1 )
element = MixedElement( [P_omega_n, P_z_n] )
Q = FunctionSpace(mesh, element)
Q_omega_n = Q.sub(0).collapse()
Q_z_n= Q.sub(1).collapse()
omega_n_1 = Function( Q_omega_n )
z_dummy = Function( Q_z_n )
psi = Function( Q )
class omega_Expression( UserExpression ):
def eval(self, values, x):
values[0] = 1.0
values[1] = 1.0
def value_shape(self):
return (2,)
omega_n, z_n = split( psi )
omega_n_0 = Function( Q_omega_n )
z_n_0 = Function( Q_z_n )
omega_n_0 = interpolate( omega_Expression( element=Q_omega_n.ufl_element() ), Q_omega_n )
omega_n_12 = (omega_n + omega_n_1) / 2.0
# set initial profiles
omega_n_1 = interpolate( omega_Expression( element=Q_omega_n.ufl_element() ), Q_omega_n )
assigner = FunctionAssigner( Q, [Q_omega_n, Q_z_n] )
# Time-stepping
for step in range( 2 ):
#for step = 1, omega_n should contain omega_n_0
XDMFFile('check_omega_n_step' + str( step ) + '.xdmf' ).write( project(omega_n, Q_omega_n) )
#for step = 1, omega_n_1 should contain omega_n_0
XDMFFile('check_omega_n_1_step' + str( step ) + '.xdmf' ).write( omega_n_1 )
#for step = 1, omega_n_12 should containg the average between omega_n and omega_n_1, i.e., omega_n_0, but it does not
XDMFFile('check_omega_n_12_step' + str( step ) + '.xdmf' ).write( project(omega_n_12, Q_omega_n) )
#write omega_n_0, z_n_0 into psi
assigner.assign(psi, [omega_n_0, z_n_0])
# copy psi into omega_n_1, z_dummy
omega_n_1, z_dummy = psi.split( deepcopy=True )
If run it with pyhon3 mwe.py
and look at the output files at step = 1
, I see that omega_n
and omega_n_1
are both equal to the vector (1,1), while omega_n_12
, which has been defined as their average in omega_n_12 = (omega_n + omega_n_1) / 2.0
, is equal to (1/2,1/2) (so it differs from their average!), see screenshots attached.
May you please explain to me what is going on ? If I put the definition omega_n_12 = (omega_n + omega_n_1) / 2.0
in the time loop, right after for step in range( 2 ):
, this behavior goes away.
Thank you for your help.