Support for quadrilateral and hexahedral meshes

#73 Declined
  1. Ivan Yashchuk

FFC now accepts "Q" and "TensorProductElement" element families. Changes made in to handle TensorProductElements Added quad and hex test cases to test/unit/misc/

With this PR + the corresponding PR in FIAT ( following commands work correctly:

assemble(1.0 * dx(mesh_quad))
assemble(1.0 * dx(mesh_hex))
assemble(Expression('x[...]', degree = 1) * dx(mesh_quad))

Simple Poisson equation can be solved on a quad mesh, but boundary conditions need to be applied manually (bc.apply does not work). Example:

from dolfin import *
import numpy as np

mesh = UnitQuadMesh(mpi_comm_world(), 2, 2)

# Build the set of vertices where bcs should be prescribed
node_set = VertexFunction('size_t', mesh, 0)
bc_boundary = CompiledSubDomain('near(x[0], 0)')
bc_boundary.mark(node_set, 1)

# Only keep vertices marked as 1
node_set = [v.index() for v in SubsetIterator(node_set, 1)]

V = FunctionSpace(mesh, "Lagrange", 1)
dof_set = np.array(vertex_to_dof_map(V)[node_set], dtype='intc')
bc_u = Constant(0.0)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx
L = 1.0*v*dx

A = assemble(a)
# Modify A: zero bc row & set diagonal to 1

b = assemble(L)
# Modify b: entry in the bc row is taken from bc_u
bc_values = interpolate(bc_u, V).vector().array()
b_values = b.array()
b_values[dof_set] = bc_values[dof_set]

u = Function(V)
solve(A, u.vector(), b)

This PR takes care of getting the data from FIAT and computing intermediate representation for the code generation for quadhex case

UPD: together with and topological DirichletBC will work correctly for quads and hexes

  • Commit status

Comments (13)

  1. Jan Blechta

    Maybe bc.apply does not work because there is no/wrong nodal interpolation (see the FIAT PR). You can always project to target space before.

    bc_value = ...
    bc_value = project(bc_value, V)
    bc = DirichletBC(V, bc_value, where)
  2. Jan Blechta

    This seems broken to me. With

    from dolfin import *
    mesh = UnitQuadMesh(mpi_comm_world(), 3, 3)
    FunctionSpace(mesh, "Lagrange", 1)

    I get

      File "/vivid/home/jan/dev/fenics-master/src/ffc/ffc/", line 69, in generate_code
        for ir in ir_finite_elements]
      File "/vivid/home/jan/dev/fenics-master/src/ffc/ffc/", line 69, in <listcomp>
        for ir in ir_finite_elements]
      File "/vivid/home/jan/dev/fenics-master/src/ffc/ffc/", line 286, in _generate_finite_element_code
        return ufc_finite_element().generate_snippets(L, ir, parameters)
      File "/vivid/home/jan/dev/fenics-master/src/ffc/ffc/uflacs/backends/ufc/", line 131, in generate_snippets
        value = method(*args)
      File "/vivid/home/jan/dev/fenics-master/src/ffc/ffc/uflacs/backends/ufc/", line 405, in interpolate_vertex_values
        element_cellname = ir["evaluate_basis"]["cellname"]
    TypeError: string indices must be integers
    1. Ivan Yashchuk author

      evaluate_basis is not implemented for tensor product elements. This line should be something like `element_cellname = ir["cellname"]. I've done something wrong while splitting the PR into two parts. I merged master branch here, now it is working.

  3. Jan Blechta

    The following script adds some plotting and DirichletBC usage:

    from dolfin import *
    import numpy as np
    from matplotlib import pyplot
    import dijitso
    N = 64
    warp_factor = 0.5
    def warp(coords):
        for c in coords:
            c += [warp_factor*np.sin(c[1]), 0]
    mesh = UnitQuadMesh(N, N)
    V = FunctionSpace(mesh, "Lagrange", 1)
    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(grad(u), grad(v))*dx
    L = v*dx
    bc_u = interpolate(Expression("0.0", element=V.ufl_element()), V)
    bc = DirichletBC(V, bc_u, lambda x, b: b)
    u = Function(V)
    solve(a == L, u, bc)
    mesh = UnitSquareMesh(N, N)
    W = FunctionSpace(mesh, "P", 1)
    v = Function(W)
    x = u.vector()
    y = v.vector()
    V_v2d = vertex_to_dof_map(V)
    W_v2d = vertex_to_dof_map(W)
    for vert in vertices(mesh):
        y[W_v2d[vert.index()]] = x[V_v2d[vert.index()]][0]
    plot(v, mode='warp')

    The result looks quite like Poisson solution by the eye-bulb error estimate:


    1. Ivan Yashchuk author

      This pull request is not updated yet. I want to use pipelines for running tests and it would just fail on importing stage since there is no FlattenedTensorProduct class in FIAT master