Hi! I am having a trouble with a code and here is the problem:
I want to solve the Poisson equation with a zero Dirichlet BC on a square and a disk. From the FEniCS (dolfin-version 2019.2.0.dev0) tutorial site, I imported the code for the square, modified the EBC (essential boundary condition) on the square and executed the script. Everything went smooth. Here is the code:
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)
# Define boundary condition
u_D = Expression('0', degree=2)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, bc)
Then I created a disk in gmsh (version 4.7.1), meshed it, created the .msh file and imported that to FEniCS and kept the remaining code (that I used for the square) unchanged, except for some modifications to accommodate the new geometry and executed it. I encountered an error that says I am integrating tensors which is not permissible (this error didn’t show up in the square case although I kept the variational form unchanged). Here is the code followed by the error:
from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
import numpy
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh
import meshio
msh = meshio.read("circle.msh")
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
#ds_custom1 = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=1)
#ds_custom2 = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)
#print(assemble(1*ds_custom(1)) ,assemble(1*ds_custom(2)))
# Create mesh and define function space
V = VectorFunctionSpace(mesh, "Lagrange", 1)
# Definining the normal to each of the faces
n=FacetNormal(mesh)
#ds = Measure("ds", subdomain_data=boundary_markers)
u_D = Constant((0.0,0.0))
bc = DirichletBC(V, u_D, mf, 1)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, bc)
Error message:
Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2, 2) and free indices with labels ().
Traceback (most recent call last):
File "poisson_test.py", line 62, in <module>
a = dot(grad(u), grad(v))*dx
File "/usr/lib/python3/dist-packages/ufl/measure.py", line 406, in __rmul__
error("Can only integrate scalar expressions. The integrand is a "
File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2, 2) and free indices with ()
I don’t quite get what the problem is. Any help or suggestion would be much appreciated.
