Modify uflacs to support quadhex elements

#77 Merged at 2aadf3f
Repository
Branch
ivan/uflacs-quad-hex-support
Repository
Branch
master
Author
  1. Ivan Yashchuk
Reviewers
Description

Modifications so that dofmap is generated not based on topological dimension but based on the cell shape. This PR together with https://bitbucket.org/fenics-project/ffc/pull-requests/73 https://bitbucket.org/fenics-project/fiat/pull-requests/38 makes following code work correctly

from dolfin import *
mesh_hex = UnitHexMesh(mpi_comm_world(), 10, 10, 10)
assemble(Expression('x[0]', degree = 1) * dx(mesh_hex))  # Returns 0.5
assemble(Expression('x[1]', degree = 1) * dx(mesh_hex))  # Returns 0.5
assemble(Expression('x[2]', degree = 1) * dx(mesh_hex))  # Returns 0.5

interpolate(Expression('...', degree = ...), V) works correctly both for quadrilateral and hexahedron case. Can be checked by

from dolfin import *
parameters["form_compiler"]["representation"] = "uflacs"
mesh = UnitHexMesh(mpi_comm_world(), 1, 1, 1)
V = FunctionSpace(mesh, "Lagrange", 2)
u = interpolate(Expression('x[0]*x[0] + 2*x[1]*x[1] + 3*x[2]*x[2]', degree=2), V)
vertex_values = u.compute_vertex_values()

for i, x in enumerate(mesh.coordinates()):
    print('vertex %d: vertex_values[%d] = %g \t point(%s)' %
        (i, i, vertex_values[i], x))

It returns

vertex 0: vertex_values[0] = 0   point([ 0.  0.  0.])
vertex 1: vertex_values[1] = 1   point([ 1.  0.  0.])
vertex 2: vertex_values[2] = 2   point([ 0.  1.  0.])
vertex 3: vertex_values[3] = 3   point([ 1.  1.  0.])
vertex 4: vertex_values[4] = 3   point([ 0.  0.  1.])
vertex 5: vertex_values[5] = 4   point([ 1.  0.  1.])
vertex 6: vertex_values[6] = 5   point([ 0.  1.  1.])
vertex 7: vertex_values[7] = 6   point([ 1.  1.  1.])

Which is correct for the given function.

  • Commit status

Comments (13)

    1. Ivan Yashchuk author

      It was named differently in other places. For example in representation.py line 215 it is cell_shape and then line 705 it is cellname. I don't know which option is more descriptive.

        1. Martin Sandve Alnæs

          cellname is used in ufl, shape in ufc, it's crashing at some point anyway.

          Trying to clean it up partially will only propagate the inconsistency to some other place in the code.

          Myself I used the ufl cellname but didn't care to touch the internal ffc ir.

          Feel free to change it back and forth if you think that's fun ;)

  1. Jack S. Hale

    I'm happy for this to go in. Clearly the mapping functionality is broken with respect to higher-order cells e.g. affine_weights and needs revisiting.

  2. Martin Sandve Alnæs

    Comments on affine_weights are still relevant as long as that function is used anywhere (whatever you rename it to). That function is a hack to inject an affine basis into the code generation for affine mesh mappings. Anywhere it's used there is an assumption of affine mesh.

    1. Jan Blechta

      Ok, UFLACS already supports non-affine geometry, but maybe only in integral generation and not element generation (is it true). So should quad and hex go through a different codepath?

      1. Martin Sandve Alnæs

        All element code should go through new codepaths. The split into separate reference basis + transform is all about that. Getting rid of any place that uses affine_weights is part of the completion of that work.