MultiMesh.h

// Copyright (C) 2014-2016 Anders Logg
//
// This file is part of DOLFIN.
//
// DOLFIN is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// DOLFIN is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
//
// Modified by August Johansson 2018
//
// First added:  2014-03-03
// Last changed: 2018-04-03

#ifndef __MULTI_MESH_H
#define __MULTI_MESH_H

#include <memory>
#include <vector>
#include <map>
#include <deque>

#include <dolfin/common/Variable.h>
#include <dolfin/geometry/Point.h>

namespace dolfin
{

  // Forward declarations
  class Cell;
  class Mesh;
  class BoundaryMesh;
  class BoundingBoxTree;
  class SimplexQuadrature;


  class MultiMesh : public Variable
  {
  public:

    typedef std::pair<std::vector<double>, std::vector<double> > quadrature_rule;

    typedef std::vector<Point> Simplex;

    typedef std::pair<std::vector<Simplex>, std::set<std::size_t>> Polyhedron;

    typedef std::vector<std::size_t> IncExcKey;

    MultiMesh();

    MultiMesh(std::vector<std::shared_ptr<const Mesh>> meshes,
              std::size_t quadrature_order);

    //--- Convenience constructors ---

    MultiMesh(std::shared_ptr<const Mesh> mesh_0,
              std::size_t quadrature_order);

    MultiMesh(std::shared_ptr<const Mesh> mesh_0,
              std::shared_ptr<const Mesh> mesh_1,
              std::size_t quadrature_order);

    MultiMesh(std::shared_ptr<const Mesh> mesh_0,
              std::shared_ptr<const Mesh> mesh_1,
              std::shared_ptr<const Mesh> mesh_2,
              std::size_t quadrature_order);

    ~MultiMesh();

    std::size_t num_parts() const;

    std::shared_ptr<const Mesh> part(std::size_t i) const;

    const std::vector<unsigned int>& uncut_cells(std::size_t part) const;

    const std::vector<unsigned int> cut_cells(std::size_t part) const;

    const std::vector<unsigned int>& covered_cells(std::size_t part) const;

    void mark_covered(std::size_t part, const std::vector<unsigned int>& cells);

    const std::map<unsigned int,
                   std::vector<std::pair<std::size_t, unsigned int> > >&
    collision_map_cut_cells(std::size_t part) const;

    const std::map<unsigned int, quadrature_rule >&
    quadrature_rules_cut_cells(std::size_t part) const;

    const quadrature_rule
    quadrature_rules_cut_cells(std::size_t part, unsigned int cell_index) const;

    const std::map<unsigned int, std::vector<quadrature_rule> >&
    quadrature_rules_overlap(std::size_t part) const;

    //      cell (unsigned int)
    //          The cell index
    const std::vector<quadrature_rule>
    quadrature_rules_overlap(std::size_t part, unsigned int cell) const;

    const std::map<unsigned int, std::vector<quadrature_rule> >&
    quadrature_rules_interface(std::size_t part) const;


    const std::vector<quadrature_rule>
    quadrature_rules_interface(std::size_t part,
                   unsigned int cell_index) const;


    const std::map<unsigned int, std::vector<std::vector<double> > >&
    facet_normals(std::size_t part) const;

    std::shared_ptr<const BoundingBoxTree>
    bounding_box_tree(std::size_t part) const;

    std::shared_ptr<const BoundingBoxTree>
    bounding_box_tree_boundary(std::size_t part) const;

    void add(std::shared_ptr<const Mesh> mesh);

    void build(std::size_t quadrature_order=2);

    bool is_built() const { return _is_built; }

    void clear();

    static Parameters default_parameters()
    {
      Parameters p("multimesh");

      //p.add("quadrature_order", 1);
      p.add("compress_volume_quadrature", false);
      p.add("compress_interface_quadrature", false);

      return p;
    }

    //--- The functions below are mainly useful for testing/debugging ---

    double compute_area() const;

    double compute_volume() const;

    std::string plot_matplotlib(double delta_z=1,
                const std::string& filename="") const;

    void auto_cover(std::size_t p, const Point& point);

  private:

    // Flag for whether multimesh has been built
    bool _is_built;

    // List of meshes
    std::vector<std::shared_ptr<const Mesh> > _meshes;

    // List of boundary meshes
    std::vector<std::shared_ptr<BoundaryMesh> > _boundary_meshes;

    // List of bounding box trees for meshes
    std::vector<std::shared_ptr<BoundingBoxTree> > _trees;

    // List of bounding box trees for boundary meshes
    std::vector<std::shared_ptr<BoundingBoxTree> > _boundary_trees;

    // Cell indices for all uncut cells for all parts. Access data by
    //
    //     c = _uncut_cells[i][j]
    //
    // where
    //
    //     c = cell index for an uncut cell
    //     i = the part (mesh) number
    //     j = the cell number (in the list of uncut cells)
    std::vector<std::vector<unsigned int> > _uncut_cells;

    // Cell indices for all covered cells for all parts. Access data by
    //
    //     c = _covered_cells[i][j]
    //
    // where
    //
    //     c = cell index for a covered cell
    //     i = the part (mesh) number
    //     j = the cell number (in the list of covered cells)
    std::vector<std::vector<unsigned int> > _covered_cells;

    // Developer note 1: The data structures _collision_map_cut_cells
    // and _quadrature_rules_cut_cells may be changed from maps to
    // vectors and indexed by the number of the cut cell (in the list
    // of cut cells), instead of indexed by the local cell index as
    // here, if we find that this is important for performance.
    //
    // Developer note 2: Quadrature points are naturally a part of a
    // form (or a term in a form) and not a part of a mesh. However,
    // for now we use a global (to the multimesh) quadrature rule for
    // all cut cells, for simplicity.

    // Collision map for cut cells. Access data by
    //
    //     c = _collision_map_cut_cells[i][j][k]
    //
    // where
    //
    //     c.first  = part number for the cutting mesh
    //     c.second = cell index for the cutting cell
    //            i = the part (mesh) number
    //            j = the cell number (local cell index
    //            k = the collision number (in the list of cutting cells)
    std::vector<std::map<unsigned int,
                         std::vector<std::pair<std::size_t, unsigned int> > > >
    _collision_maps_cut_cells;

    // Quadrature rules for cut cells. Access data by
    //
    //     q = _quadrature_rules_cut_cells[i][j]
    //
    // where
    //
    //     q.first  = quadrature weights, array of length num_points
    //     q.second = quadrature points, flattened num_points x gdim array
    //            i = the part (mesh) number
    //            j = the cell number (local cell index)
    std::vector<std::map<unsigned int, quadrature_rule> >
    _quadrature_rules_cut_cells;

    // Quadrature rules for overlap. Access data by
    //
    //     q = _quadrature_rules_overlap[i][j][k]
    //
    // where
    //
    //     q.first  = quadrature weights, array of length num_points
    //     q.second = quadrature points, flattened num_points x gdim array
    //            i = the part (mesh) number
    //            j = the cell number (local cell index)
    //            k = the collision number (in the list of cutting cells)
    std::vector<std::map<unsigned int, std::vector<quadrature_rule> > >
    _quadrature_rules_overlap;

    // Quadrature rules for interface. Access data by
    //
    //     q = _quadrature_rules_interface[i][j][k]
    //
    // where
    //
    //     q.first  = quadrature weights, array of length num_points
    //     q.second = quadrature points, flattened num_points x gdim array
    //            i = the part (mesh) number
    //            j = the cell number (local cell index)
    //            k = the collision number (in the list of cutting cells)
    std::vector<std::map<unsigned int, std::vector<quadrature_rule> > >
    _quadrature_rules_interface;

    // Facet normals for interface. Access data by
    //
    //     n = _facet_normals_interface[i][j][k][
    //
    // where
    //
    //     n = a flattened array vector of facet normals, one point for
    //         each quadrature point
    //     i = the part (mesh) number
    //     j = the cell number (local cell index)
    //     k = the collision number (in the list of cutting cells)
    std::vector<std::map<unsigned int, std::vector<std::vector<double> > > >
    _facet_normals;

    // Build boundary meshes
    void _build_boundary_meshes();

    // Build bounding box trees
    void _build_bounding_box_trees();

    // Build collision maps
    void _build_collision_maps();
    //void _build_collision_maps_same_topology();
    //void _build_collision_maps_different_topology();

    // Build quadrature rules for the cut cells
    void _build_quadrature_rules_cut_cells(std::size_t quadrature_order);

    // Build quadrature rules for the overlap
    void _build_quadrature_rules_overlap(std::size_t quadrature_order);

    // Build quadrature rules and normals for the interface
    void _build_quadrature_rules_interface(std::size_t quadrature_order);

    // Help function to determine if interface intersection is
    // (exactly) overlapped by a cutting cell
    bool _is_overlapped_interface(std::vector<Point> simplex,
                  const Cell cut_cell,
                  Point simplex_normal) const;

    // Add quadrature rule for simplices in polyhedron. Returns the
    // number of points generated for each simplex.
    std::size_t _add_quadrature_rule(quadrature_rule& qr,
                     const SimplexQuadrature& sq,
                     const Simplex& simplex,
                     std::size_t gdim,
                     std::size_t quadrature_order,
                     double factor) const;

    // Add quadrature rule to existing quadrature rule (append dqr to
    // qr). Returns number of points added.
    std::size_t _add_quadrature_rule(quadrature_rule& qr,
                                     const quadrature_rule& dqr,
                                     std::size_t gdim,
                                     double factor) const;

    // Append normal to list of normals npts times
    void _add_normal(std::vector<double>& normals,
                     const Point& normal,
                     std::size_t npts,
                     std::size_t gdim) const;

    // Plot multimesh
    void _plot() const;

    // Inclusion-exclusion for overlap
    void _inclusion_exclusion_overlap
      (std::vector<quadrature_rule>& qr,
       const SimplexQuadrature& sq,
       const std::vector<std::pair<std::size_t, Polyhedron> >& initial_polyhedra,
       std::size_t tdim,
       std::size_t gdim,
       std::size_t quadrature_order) const;

    // Inclusion-exclusion for interface
    void _inclusion_exclusion_interface
      (quadrature_rule& qr,
       std::vector<double>& normals,
       const SimplexQuadrature& sq,
       const Simplex& Eij,
       const Point& facet_normal,
       const std::vector<std::pair<std::size_t, Polyhedron> >& initial_polygons,
       std::size_t tdim,
       std::size_t gdim,
       std::size_t quadrature_order) const;

    // Construct and return mapping from boundary facets to full mesh
    std::vector<std::vector<std::pair<std::size_t, std::size_t> > >
      _boundary_facets_to_full_mesh(std::size_t part) const;

    // Impose consistency of _cut_cells, so that only the cells with
    // a nontrivial interface quadrature rule are classified as cut.
    void _impose_cut_cell_consistency();

    // Remove quadrature rule if the sum of the weights is less than a
    // tolerance
    static void remove_quadrature_rule(quadrature_rule& qr,
                       double tolerance);
  };



}


#endif