Computing function map (coordinates : values) over a bounding box fails in MPI

1

Hi community,

I am trying to reduce computational time for a specific code script by creating a map for a mesh - coordinates : values - over an independent overlapping bounding box. The following code script does that. It works perfectly fine in serial. However the code hangs up when running in parallel.

If you play around with the bounding box, you observe that it works well in MPI with the bounding box coinciding with all the mesh partitions. The hang up is observed if any mesh partition doesn’t see the bounding box at all.

I am using fenics 2016.1.0 to run the script using docker on an Ubuntu 18.04 terminal. I am running the code using the command mpirun.mpich -n 4 python untitled.py. (Note: Put both untitled.py and untitled.cpp in the same directory)

untitled.py

from dolfin import *

mesh1 = UnitSquareMesh(11, 11) 
V1 = VectorFunctionSpace(mesh1, 'CG', 2)
u1 = Function(V1)
u1 = interpolate(Expression((("x[0]", "x[1]")), degree = 2),V1)

fem_code = open('untitled.cpp', 'r').read()
compiled_fem_module = compile_extension_module(code=fem_code)

compiled_fem_module.delta(u1)

untitled.cpp

#include <dolfin/mesh/Vertex.h>

namespace dolfin
{

  bool in_bounding_box(const std::vector<double>& point,
                                            const std::vector<double>& bounding_box,
                                            const double tol)
  {
    // Return false if bounding box is empty
    if (bounding_box.empty())
      return false; 

    const std::size_t gdim = point.size();
    dolfin_assert(bounding_box.size() == 2*gdim);
    for (std::size_t i = 0; i < gdim; ++i)
    {
      if (!(point[i] >= (bounding_box[i] - tol)
        && point[i] <= (bounding_box[gdim + i] + tol)))
      {
        return false;
      }
    }

    return true;   
  }

std::map<std::vector<double>, std::vector<double>>
  tabulate_coordinates_to_values(const Function& u0, std::vector<double>& x_x)
  {
  	std::map<std::vector<double>, std::vector<double>>
    	coords_to_values;        

    // Extract mesh
    const FunctionSpace& V0 =  *u0.function_space();
    const Mesh& mesh0 = *u0.function_space()->mesh();
    const std::size_t gdim0 = mesh0.geometry().dim();

    // Create arrays used to evaluate points on mesh0
    std::vector<double> node_coord(gdim0);
  	std::vector<double> node_val(u0.value_size());
  	Array<double> _node_coord(gdim0, node_coord.data());
  	Array<double> _node_val(u0.value_size(), node_val.data());

    for (VertexIterator vert(mesh0); !vert.end(); ++vert)    
    {      
      for(int i = 0; i < gdim0; i++)
        node_coord[i] = vert->x(i);
      
      if(in_bounding_box(node_coord, x_x, 1e-12))  
      {       
         u0.eval(_node_val, _node_coord);
      
         coords_to_values.insert
          (std::make_pair(node_coord, node_val));
      }
    }

	return coords_to_values;
  }

  void delta(const Function& u0) 
  {
    // Function u0
    const FunctionSpace& V0 =  *u0.function_space();
    const Mesh& mesh0 = *V0.mesh();
    const std::size_t gdim0 = mesh0.geometry().dim();  

    // Create arrays used to evaluate points on mesh0
    std::vector<double> x(gdim0);
    std::vector<double> val(u0.value_size());
    Array<double> _x(gdim0, x.data());
    Array<double> _val(u0.value_size(), val.data()); 

    //  Bounding box
    std::vector<double> bounding_box = {0.65, 0.65, 0.9, 0.9};

    // Map for u0 encompassing the region of the bounding box
    std::map<std::vector<double>, std::vector<double>> 
          coords_to_values = tabulate_coordinates_to_values(u0, bounding_box);

    // .................. Post processing ....................    

    for(const auto &map : coords_to_values)
    {
      std::copy(map.first.begin(), map.first.end(), x.begin());
      std::copy(map.second.begin(), map.second.end(), val.begin());
      std::cout << _x << "                " << _val << "\n";        
    }      
          
  }
}

Pls. suggest a workaround to this problem.