Good evening,
I was finishing a code project in FEniCS, all was working correctly, but at a moment, I think I suppressed a line code, and it doesn’t work anymor …
I have a UFL Exception and I do not understand why.
Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2, 2) and free indices with labels ().
---------------------------------------------------------------------------
UFLException Traceback (most recent call last)
<ipython-input-63-802ab7e871a6> in <module>
3 u, p = TrialFunctions(W)
4
----> 5 a1 = nu*inner(grad(u), grad(v))*dx - inner(p, div(v))*dx - inner(q, div(u))*dx
6 l1 = inner(Constant((0, 0)), v)*dx
/usr/lib/python3/dist-packages/ufl/measure.py in __rmul__(self, integrand)
418 error("Can only integrate scalar expressions. The integrand is a "
419 "tensor expression with value shape %s and free indices with labels %s." %
--> 420 (integrand.ufl_shape, integrand.ufl_free_indices))
421
422 # If we have a tuple of domain ids, delegate composition to
/usr/lib/python3/dist-packages/ufl/log.py in error(self, *message)
170 "Write error message and raise an exception."
171 self._log.error(*message)
--> 172 raise self._exception_type(self._format_raw(*message))
173
174 def begin(self, *message):
UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2, 2) and free indices with labels ().The problem seems to come from my functions space, but I didn’t modify them ! It is surely a very stupid mistake, but I can’t see it anymore.
My complete code is there :
#C0 Modules
from dolfin import *
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt # pour les graphiques
from mshr import * # pour le maillage
from numpy import * # pour la partie aléatoire
#C1 Constantes (adimensionnées)
nu=1E-5 #viscosité
Pg=10+1E-2 #pressions à gauche
Pd=10 #pressions à droite
tol=1E-14
#C2 Géométrie
X=10 #longueur du tuyau
Y=2 #diamètre
# Trous aléatoires
param_poisson = 2
param_exp = 0.5
# Initialisation
nb_trous = random.poisson(param_poisson)
r = []
xcentre = []
ycentre = []
# Paramètre des trous
for i in range(nb_trous):
r.append(random.exponential(param_exp))
xcentre.append(random.uniform(2, 8))
ycentre.append(random.uniform(-1, 1))
print("xcentre",xcentre)
print("ycentre",ycentre)
print("rayon",r)
#C3 Maillage
nx=50
base = Rectangle(Point(0.0, -1.0), Point(X, Y-1.0))
for i in range(nb_trous):
trou = Circle(Point(xcentre[i], ycentre[i]), r[i])
base = base - trou
mesh2 = generate_mesh(base, nx)
plt.figure(1)
plot(mesh2, title="Maillage du rectangle")
#C4 Espace de fonctions
Ve = VectorElement("P",triangle, 2) #Espace de la fonction vectorielle v (vitesse)
Qe = FiniteElement("P",triangle, 1) #Espace de la fonction scalaire p (pression)
W = FunctionSpace(mesh2, MixedElement([Ve, Qe])) #Espace produit pour la résolution
#C5 Conditions aux bords
class Bord(SubDomain): #Bords haut et bas du tuyau, où sont les parois, et bords des obstacles
def inside(self, x, on_boundary):
# Bords haut et bas
haut_bas = ((abs(x[1] - 1) < tol) or (abs(x[1] + 1) < tol))
# Obstacles
tot_obst=False
for i in range(nb_trous):
obst = ((x[0]-xcentre[i])**2 + (x[1]-ycentre[i])**2 < (r[i]+tol)**2)
tot_obst = tot_obst or obst
tot_obst = tot_obst or haut_bas
return on_boundary and tot_obst
class Entree(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (abs(x[0] - 0) < tol)
class Sortie(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (abs(x[0] - 10) < tol)
# Conditions de Dirichlet sur la vitesse
v_0 = (0.125*(Pg-Pd)*nu*3.14/X, 0) #vitesse débitante moyenne de l'écoulement en entrée/sortie
b_fin = Entree()
b_fex = Sortie()
finD = DirichletBC(W.sub(0), v_0, b_fin)
fexD = DirichletBC(W.sub(0), v_0, b_fex)
b_bord = Bord()
bordD = DirichletBC(W.sub(0), (0, 0), b_bord)
# Conditions de Dirichlet sur la pression
b_in = Entree()
b_ex = Sortie()
inD = DirichletBC(W.sub(1), Constant(Pg), b_in)
exD = DirichletBC(W.sub(1), Constant(Pd), b_ex)
bcu = [bordD, inD, exD, finD, fexD] #total des conditions
class BordFlux(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and ((abs(x[0] - 0) < tol) or (abs(x[0] - 10) < tol))
bfN = BordFlux()
boundaries = MeshFunction("size_t", mesh2, mesh2.topology().dim()-1, 0)
bfN.mark(boundaries, 0)
ds = Measure("ds", domain=mesh2, subdomain_data=boundaries)
#C6 Formulation variationnelle
v, q = TestFunctions(W)
u, p = TrialFunctions(W)
a1 = nu*inner(grad(u), grad(v))*dx - inner(p, div(v))*dx - inner(q, div(u))*dx
l1 = inner(Constant((0, 0)), v)*dx
#C7 Résolution de l'approximation variationnelle
up = Function(W)
u, p = split(up)
solve(a1 == l1, up, bcu)
#C8 Tracer la solution numérique (vitesse)
plt.figure(2)
plu=plot(u)
plt.title('Vitesse v')
plt.xlabel('x')
plt.ylabel('y')
plt.colorbar(plu)
plt.figure(3)
plu2=plot(u[0])
plt.title('Vitesse axiale v_x')
plt.xlabel('x')
plt.ylabel('y')
plt.colorbar(plu2)
#C8 Tracer la solution numérique (pression)
plt.figure(2)
plu=plot(p)
plt.title('Pression p')
plt.xlabel('x')
plt.ylabel('y')
plt.colorbar(plu)Problems are at C6 and I think they come from C3 and C4 but I am not sure
I thank you for your help, I am completely getting mad
Hyanoa