I was trying to test an electrostatic solver with the mixed method described in section 2.2.2 of the FEniCS Book, hoping that it would just be a cut-and-paste job, but I have run into an error that I can't seem to track down:
ufl.log.UFLException: All terms in form must have same rank.
I am hoping that it is something simple that I have just messed up in the "a" Form in the code below ...
I am loading meshes from files I create externally, and have pared them down to the point where they aren't too large. I have also posted a smaller version of the code that also produces the error (data follows the code):
#-------------------------------#
# Construct Mesh, FunctionSpace #
#-------------------------------#
mesh = Mesh("test.xml.gz")
BDM = FunctionSpace(mesh, "BDM", 1)
DG = FunctionSpace(mesh, "DG", 0)
W = BDM * DG
#---------------------#
# Boundary Conditions #
#---------------------#
mvc = MeshValueCollection("size_t", mesh, "test_facets.xml.gz")
dirichlet_f = MeshFunction("size_t", mesh, mvc)
dirichlet_bc = DirichletBC(W.sub(0), Constant((0., 0.)), dirichlet_f, 3)
#----------------------#
# Dielectric Materials #
#----------------------#
epsilon = Function(DG)
materials = MeshFunction("size_t", mesh, "test_cells.xml.gz")
helper = numpy.asarray(materials.array(), dtype=numpy.int32) - 1
dm = DG.dofmap()
for cell in cells(mesh):
helper[dm.cell_dofs(cell.index())] = materials[cell] - 1
dielectrics = (2.2, 5., 10.)
epsilon.vector()[:] = numpy.choose(helper, dielectrics)
#------------#
# Define PDE #
#------------#
(electric_field, voltage) = TrialFunctions(W)
(test1, test2) = TestFunctions(W)
r = Expression("x[1]") # the radial coordinate is the y coordinate
n = FacetNormal(mesh)
# ds = Measure("ds")[dirichlet_f]
a = (dot(test1, electric_field + nabla_grad(voltage)) - \
epsilon*dot(nabla_grad(test2), electric_field))*r*dx + \
epsilon*dot(test1, n)*ds(1)
L = Constant(0.0)*test2*dx
w = Function(W)
solve(a == L, w, dirichlet_bc)
(electric_field, voltage) = w.split()
# Save solution in VTK format
f = File("problem.pvd")
electric_field.rename("Electric_Field", electric_field.name())
voltage.rename("Voltage", voltage.name())
f << (voltage, 0.)
f << (electric_field, 0.)
plot(voltage, interactive=True)
File 1/3 "test.xml.gz" 46 lines of text (without the gzip, obviously):
<?xml version="1.0" encoding="UTF-8"?>
<dolfin xmlns:dolfin="http://fenicsproject.org">
<mesh celltype="triangle" dim="2">
<cells size="21">
<triangle index="0" v0="0" v1="1" v2="2" />
<triangle index="1" v0="1" v1="3" v2="4" />
<triangle index="2" v0="1" v1="2" v2="5" />
<triangle index="3" v0="0" v1="4" v2="6" />
<triangle index="4" v0="3" v1="5" v2="7" />
<triangle index="5" v0="1" v1="3" v2="5" />
<triangle index="6" v0="0" v1="1" v2="4" />
<triangle index="7" v0="8" v1="9" v2="10" />
<triangle index="8" v0="5" v1="7" v2="11" />
<triangle index="9" v0="9" v1="10" v2="12" />
<triangle index="10" v0="2" v1="5" v2="12" />
<triangle index="11" v0="5" v1="11" v2="12" />
<triangle index="12" v0="2" v1="9" v2="12" />
<triangle index="13" v0="10" v1="11" v2="12" />
<triangle index="14" v0="0" v1="13" v2="14" />
<triangle index="15" v0="2" v1="9" v2="13" />
<triangle index="16" v0="8" v1="9" v2="15" />
<triangle index="17" v0="9" v1="13" v2="15" />
<triangle index="18" v0="0" v1="6" v2="14" />
<triangle index="19" v0="0" v1="2" v2="13" />
<triangle index="20" v0="13" v1="14" v2="15" />
</cells>
<vertices size="16">
<vertex index="0" x="2.50000000e+00" y="0.00000000e+00" />
<vertex index="1" x="1.18969262e+00" y="2.06060806e+00" />
<vertex index="2" x="0.00000000e+00" y="0.00000000e+00" />
<vertex index="3" x="8.68240888e-01" y="4.92403877e+00" />
<vertex index="4" x="3.83022222e+00" y="3.21393805e+00" />
<vertex index="5" x="-1.25000000e+00" y="2.16506351e+00" />
<vertex index="6" x="5.00000000e+00" y="0.00000000e+00" />
<vertex index="7" x="-2.50000000e+00" y="4.33012702e+00" />
<vertex index="8" x="-2.50000000e+00" y="-4.33012702e+00" />
<vertex index="9" x="-1.25000000e+00" y="-2.16506351e+00" />
<vertex index="10" x="-4.69846310e+00" y="-1.71010072e+00" />
<vertex index="11" x="-4.69846310e+00" y="1.71010072e+00" />
<vertex index="12" x="-2.37938524e+00" y="-3.10862447e-16" />
<vertex index="13" x="1.18969262e+00" y="-2.06060806e+00" />
<vertex index="14" x="3.83022222e+00" y="-3.21393805e+00" />
<vertex index="15" x="8.68240888e-01" y="-4.92403877e+00" />
</vertices>
</mesh>
</dolfin>
File 2/3 "test_facets.xml.gz" 14 lines of text
<?xml version="1.0" encoding="UTF-8"?>
<dolfin xmlns:dolfin="http://fenicsproject.org">
<mesh_value_collection name="boundary_markers" type="uint" dim="1" size="9">
<value cell_index="16" local_entity="1" value="3" />
<value cell_index="18" local_entity="0" value="3" />
<value cell_index="13" local_entity="2" value="2" />
<value cell_index="7" local_entity="1" value="2" />
<value cell_index="3" local_entity="0" value="1" />
<value cell_index="1" local_entity="0" value="1" />
<value cell_index="8" local_entity="0" value="2" />
<value cell_index="4" local_entity="1" value="1" />
<value cell_index="20" local_entity="0" value="3" />
</mesh_value_collection>
</dolfin>
File 3/3 "test_cells.xml.gz" 28 lines of text
<?xml version="1.0"?>
<dolfin xmlns:dolfin="http://fenicsproject.org">
<mesh_function>
<mesh_value_collection name="cells" type="uint" dim="2" size="21">
<value cell_index="0" local_entity="0" value="1" />
<value cell_index="1" local_entity="0" value="1" />
<value cell_index="2" local_entity="0" value="1" />
<value cell_index="3" local_entity="0" value="1" />
<value cell_index="4" local_entity="0" value="1" />
<value cell_index="5" local_entity="0" value="1" />
<value cell_index="6" local_entity="0" value="1" />
<value cell_index="7" local_entity="0" value="2" />
<value cell_index="8" local_entity="0" value="2" />
<value cell_index="9" local_entity="0" value="2" />
<value cell_index="10" local_entity="0" value="2" />
<value cell_index="11" local_entity="0" value="2" />
<value cell_index="12" local_entity="0" value="2" />
<value cell_index="13" local_entity="0" value="2" />
<value cell_index="14" local_entity="0" value="3" />
<value cell_index="15" local_entity="0" value="3" />
<value cell_index="16" local_entity="0" value="3" />
<value cell_index="17" local_entity="0" value="3" />
<value cell_index="18" local_entity="0" value="3" />
<value cell_index="19" local_entity="0" value="3" />
<value cell_index="20" local_entity="0" value="3" />
</mesh_value_collection>
</mesh_function>
</dolfin>
