# Support for quadrilateral and hexahedral meshes

#73 Declined
Repository
fenics-project
Branch
master
Reviewers
Description

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

With this PR + the corresponding PR in FIAT (https://bitbucket.org/fenics-project/fiat/pull-requests/39/) 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
A.ident_local(dof_set)
A.apply('insert')

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]
b.set_local(b_values)
b.apply('insert')

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 https://bitbucket.org/fenics-project/ffc/pull-requests/82 and https://bitbucket.org/fenics-project/dolfin/pull-requests/371 topological `DirichletBC` will work correctly for quads and hexes

• Commit status

# Comments (13)

1. 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. Please, next time do not prefix/postfix PR title with an author's name and a project name. That is not usual.

3. 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/codegeneration.py", line 69, in generate_code
for ir in ir_finite_elements]
File "/vivid/home/jan/dev/fenics-master/src/ffc/ffc/codegeneration.py", line 69, in <listcomp>
for ir in ir_finite_elements]
File "/vivid/home/jan/dev/fenics-master/src/ffc/ffc/codegeneration.py", 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/generator.py", line 131, in generate_snippets
value = method(*args)
File "/vivid/home/jan/dev/fenics-master/src/ffc/ffc/uflacs/backends/ufc/finite_element.py", line 405, in interpolate_vertex_values
element_cellname = ir["evaluate_basis"]["cellname"]
TypeError: string indices must be integers
[...]
```
1. 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.

1. Don't implement it - it's deprecated. Just implement the reference version.

4. The following script adds some plotting and `DirichletBC` usage:

```from dolfin import *
import numpy as np
from matplotlib import pyplot

import dijitso
dijitso.set_log_level("debug")

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)
warp(mesh.coordinates())
print(mesh.coordinates())

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)
warp(mesh.coordinates())
print(mesh.coordinates())

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')
pyplot.show()
```

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

5. So I see that some of the implementation details improved, but I would still like to see a FIAT/FFC refactoring over these pull requests.

6. author

Okay, I will need to move flattening of entities from FFC to a wrapper class in FIAT

1. 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

1. You have to modify `bitbucket-pipelines.yml` to use the correct branch.