Projection between meshes

[aeae]

Hello,
Hopefully someone could explain what I am doing wrong, or whether this is a bug.

I am trying to do some mesh refinement and in this connection need to transfer (e.g. project) already computed results onto a new mesh. This seems not to work with RT elements at all. Consider the following snippet:

from dolfin import *

mesh = UnitSquareMesh(20, 20)

dg0_el = FiniteElement("DG", mesh.ufl_cell(), 0)
dg0_sp = FunctionSpace(mesh,dg0_el)
dg0_f  = Function(dg0_sp)

rt0_el = FiniteElement("RT", mesh.ufl_cell(), 1)
rt0_sp = FunctionSpace(mesh,rt0_el)
rt0_f  = Function(rt0_sp)

mesh1 = refine(mesh)

# This works
dg0_sp1 = FunctionSpace(mesh1,dg0_el)
dg0_f1  = project(dg0_f,dg0_sp1)

# This fails
rt0_sp1 = FunctionSpace(mesh1,rt0_el)
rt0_f1  = project(rt0_f,rt0_sp1)

The second projection fails at the compilation stage with the output below. Any help/workarounds is highly appreciated.

------------------- Start compiler output ------------------------
/tmp/tmpptiopkeb/ffc_form_5b6cfae20d3f11d5f95e3ba07fb7d7fe36556df0.cpp: In member function ‘virtual void ffc_form_5b6cfae20d3f11d5f95e3ba07fb7d7fe36556df0_cell_integral_main_otherwise::tabulate_tensor(double*, const double* const*, const double*, int) const’:
/tmp/tmpptiopkeb/ffc_form_5b6cfae20d3f11d5f95e3ba07fb7d7fe36556df0.cpp:105:18: error: redeclaration of ‘const double J_c0’
     const double J_c0 = coordinate_dofs[0] * FE7_C0_D01_Q3[0][0][0] + coordinate_dofs[2] * FE7_C0_D01_Q3[0][0][1];
                  ^~~~
/tmp/tmpptiopkeb/ffc_form_5b6cfae20d3f11d5f95e3ba07fb7d7fe36556df0.cpp:101:18: note: ‘const double J_c0’ previously declared here
     const double J_c0 = coordinate_dofs[0] * FE7_C0_D01_Q3[0][0][0] + coordinate_dofs[2] * FE7_C0_D01_Q3[0][0][1];
                  ^~~~
/tmp/tmpptiopkeb/ffc_form_5b6cfae20d3f11d5f95e3ba07fb7d7fe36556df0.cpp:106:18: error: redeclaration of ‘const double J_c3’
     const double J_c3 = coordinate_dofs[1] * FE7_C0_D01_Q3[0][0][0] + coordinate_dofs[5] * FE7_C0_D01_Q3[0][0][1];
                  ^~~~
/tmp/tmpptiopkeb/ffc_form_5b6cfae20d3f11d5f95e3ba07fb7d7fe36556df0.cpp:102:18: note: ‘const double J_c3’ previously declared here
     const double J_c3 = coordinate_dofs[1] * FE7_C0_D01_Q3[0][0][0] + coordinate_dofs[5] * FE7_C0_D01_Q3[0][0][1];
                  ^~~~
/tmp/tmpptiopkeb/ffc_form_5b6cfae20d3f11d5f95e3ba07fb7d7fe36556df0.cpp:107:18: error: redeclaration of ‘const double J_c1’
     const double J_c1 = coordinate_dofs[0] * FE7_C0_D01_Q3[0][0][0] + coordinate_dofs[4] * FE7_C0_D01_Q3[0][0][1];
                  ^~~~
/tmp/tmpptiopkeb/ffc_form_5b6cfae20d3f11d5f95e3ba07fb7d7fe36556df0.cpp:103:18: note: ‘const double J_c1’ previously declared here
     const double J_c1 = coordinate_dofs[0] * FE7_C0_D01_Q3[0][0][0] + coordinate_dofs[4] * FE7_C0_D01_Q3[0][0][1];
                  ^~~~
/tmp/tmpptiopkeb/ffc_form_5b6cfae20d3f11d5f95e3ba07fb7d7fe36556df0.cpp:108:18: error: redeclaration of ‘const double J_c2’
     const double J_c2 = coordinate_dofs[1] * FE7_C0_D01_Q3[0][0][0] + coordinate_dofs[3] * FE7_C0_D01_Q3[0][0][1];
                  ^~~~
/tmp/tmpptiopkeb/ffc_form_5b6cfae20d3f11d5f95e3ba07fb7d7fe36556df0.cpp:104:18: note: ‘const double J_c2’ previously declared here
     const double J_c2 = coordinate_dofs[1] * FE7_C0_D01_Q3[0][0][0] + coordinate_dofs[3] * FE7_C0_D01_Q3[0][0][1];
                  ^~~~

-------------------  End compiler output  ------------------------
Compilation failed! Sources, command, and errors have been written to: /home/fenics/shared/jitfailure-ffc_form_5b6cfae20d3f11d5f95e3ba07fb7d7fe36556df0
Traceback (most recent call last):
  File "test_refine.py", line 21, in <module>
    rt0_f1  = project(rt0_f,rt0_sp1)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/projection.py", line 133, in project
    form_compiler_parameters=form_compiler_parameters)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/assembling.py", line 361, in assemble_system
    b_dolfin_form = _create_dolfin_form(b_form, form_compiler_parameters)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/assembling.py", line 58, in _create_dolfin_form
    function_spaces=function_spaces)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/form.py", line 45, in __init__
    mpi_comm=mesh.mpi_comm())
  File "/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py", line 97, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
  File "/usr/local/lib/python3.6/dist-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/usr/local/lib/python3.6/dist-packages/ffc/jitcompiler.py", line 133, in jit_build
    generate=jit_generate)
  File "/usr/local/lib/python3.6/dist-packages/dijitso/jit.py", line 217, in jit
    % err_info['fail_dir'], err_info)
dijitso.jit.DijitsoError: Dijitso JIT compilation failed, see '/home/fenics/shared/jitfailure-ffc_form_5b6cfae20d3f11d5f95e3ba07fb7d7fe36556df0' for details

[aeae]

If anyone is interested out there - this must be a be a bug. A not-so-elegant workaround is to go though the intermediate space, e.g. Vector DG0:

from dolfin import *

mesh = UnitSquareMesh(20, 20)

dg0_el = VectorElement("DG", mesh.ufl_cell(), 0)
dg0_sp = FunctionSpace(mesh,dg0_el)
dg0_f  = Function(dg0_sp)

rt0_el = FiniteElement("RT", mesh.ufl_cell(), 1)
rt0_sp = FunctionSpace(mesh,rt0_el)
rt0_f  = Function(rt0_sp)

mesh1 = refine(mesh)

# This works
dg0_sp1 = FunctionSpace(mesh1,dg0_el)
dg0_f1  = project(dg0_f,dg0_sp1)

# Workaround for RT projection
rt0_sp1 = FunctionSpace(mesh1,rt0_el)
dg0_f   = project(rt0_f,dg0_sp)
rt0_f1  = project(dg0_f,rt0_sp1)

Now this works w/o MPI. In parallel, I get the following (again, any light shed on this is appreciated):

fenics@bec66c438bc8:~/shared$ mpirun -np 2 python3 test_refine.py
Process 0: Computed global bounding box tree with 3 boxes.
Process 1: Computed global bounding box tree with 3 boxes.
Process 1: Building point search tree to accelerate distance queries.
Process 1: Computed bounding box tree with 799 nodes for 400 points.
Process 0: Building point search tree to accelerate distance queries.
Process 0: Computed bounding box tree with 799 nodes for 400 points.
Traceback (most recent call last):
  File "test_refine.py", line 17, in <module>
    dg0_f1  = project(dg0_f,dg0_sp1)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/projection.py", line 133, in project
    form_compiler_parameters=form_compiler_parameters)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/assembling.py", line 382, in assemble_system
    assembler.assemble(A_tensor, b_tensor)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to evaluate function at point.
*** Reason:  The point is not inside the domain. Consider calling "Function::set_allow_extrapolation(true)" on this Function to allow extrapolation.
*** Where:   This error was encountered inside Function.cpp.
*** Process: 0
***
*** DOLFIN version: 2018.1.0
*** Git changeset:  948dc42cc4e06ed9227d0201ad50f94ac94cbf9f
*** -------------------------------------------------------------------------

Traceback (most recent call last):
  File "test_refine.py", line 17, in <module>
    dg0_f1  = project(dg0_f,dg0_sp1)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/projection.py", line 133, in project
    form_compiler_parameters=form_compiler_parameters)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/assembling.py", line 382, in assemble_system
    assembler.assemble(A_tensor, b_tensor)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to evaluate function at point.
*** Reason:  The point is not inside the domain. Consider calling "Function::set_allow_extrapolation(true)" on this Function to allow extrapolation.
*** Where:   This error was encountered inside Function.cpp.
*** Process: 1
***
*** DOLFIN version: 2018.1.0
*** Git changeset:  948dc42cc4e06ed9227d0201ad50f94ac94cbf9f
*** -------------------------------------------------------------------------

[xiaodong_zhang]

Hello, Aeae. Did you solve the problem? I met the similat problem “Building point search tree to accelerate distance queries”. If you solved this, could you please let me know the solution?

[naghamchibli]

Hello,
We have the same problem. Any solution? Thanks!

[dokken]

You need to allow for extrapolation if projections between meshes in legacy dolfin.

In dolfinx, ive sketched out how to do this in

L2 projection of fine mesh solution to a coarse mesh

At some point you have to evaulate the function at the quadrature points of the coarse grid when doing the projection. The following code illlustrates how to compute the L2 projection from a fine grid to a coarse grid by: Get values of the fine function at the quadrature points of the coarse grid. Then do the projection # Illustration of projection from one mesh to another # Author: Jørgen S. Dokken # SPDX License: MIT import ufl import dolfinx from mpi4py import MPI import basix.ufl import d…