Support for quadrilateral and hexahedral meshes

#73 Declined
Repository
fenics-project
Branch
master
Author
  1. Ivan Yashchuk
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. 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/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. 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
    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:

    figure_1-1.png

    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