This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Function evaluates with scalar or constant field give different result

0 votes

Hello people,

To make it simple, i have defined a function D_el(x). When i evaluate it with x=C a scalar, i have the expected result. When i evaluate it with x being a field (with x a constant field x=C), i have another (then, wrong) result.

Below is a minimum code that illustrates my problem:

# Impose / is the decimal division (must be at the start of the code)
from __future__ import division

# Dolfin library
from dolfin import *
# Numpy library
import numpy as np

# Current version of FEniCS
print "The current version of FEniCS is",dolfin.dolfin_version()

# Mesh
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(10, 10, 10), 10, 10,10)

# Define function
def D_el(Ce,Temp):
    C=Ce*1e-6
    D=0.00584*exp(-2870/Temp)*((1000*C)**2) - 0.0339*exp(-2920/Temp)*1000*C + 0.129*exp(-3200/Temp)
    return(D) 

# Test function with scalar 
# Input
C=10000
Temp=298.15
# Output
D=D_el(C,Temp)
print D # Correct

# Set function space
V = FunctionSpace(mesh, 'Lagrange', 1)
# Project C on the function space
c_on_V = project(C,V)
plot(c_on_V, title='C',interactive=True) # Value display is correct
# Test function with field
D=D_el(c_on_V,Temp)
plot(D, title='D_on_V',interactive=True) # WRONG !

# Let's get the min and max of the function calculated on the field
D_on_V = project(D,V)
array_=D_on_V.vector().array()
max_=np.max(array_)
print max_ # WRONG !
min_=np.min(array_)
print min_ # WRONG !

I obtain this (all the the three number should be identical)

The current version of FEniCS is 2016.1.0
2.24353357359e-05
Object cannot be plotted directly, projecting to piecewise linears.
-1.6102060254e-05
-1.6102060254e-05

Where is my error?

asked Jun 16, 2017 by fussegli FEniCS Novice (700 points)

Ok, so i have found the error: the second input variable (Temp) must also being projected to V

 # Set function space
V = FunctionSpace(mesh, 'Lagrange', 1)
# Project C on the function space
c_on_V = project(C,V)
plot(c_on_V, title='C',interactive=True) # Value display is correct
# Project Temp on the function space
Temp_on_V = project(Temp,V)
plot(Temp_on_V, title='T',interactive=True) # Value display is correct

# Test function with field
D=D_el(c_on_V,Temp_on_V)
plot(D, title='D_on_V',interactive=True) # Correct !

# Let's get the min and max of the function calculated on the field
D_on_V = project(D,V)
array_=D_on_V.vector().array()
max_=np.max(array_)
print max_ # Correct !
min_=np.min(array_)
print min_ # Correct !

This time, i get:

The current version of FEniCS is 2016.1.0
2.24353357359e-05
Calling FFC just-in-time (JIT) compiler, this may take some time.
Object cannot be plotted directly, projecting to piecewise linears.
Calling FFC just-in-time (JIT) compiler, this may take some time.
2.24353357359e-05
2.24353357359e-05

However, i do not understand why using Temp as a scalar lead to an error. If someone can explain it to me...

...