DOLFIN
DOLFIN C++ interface
MultiMesh.h
1 // Copyright (C) 2014-2016 Anders Logg
2 //
3 // This file is part of DOLFIN.
4 //
5 // DOLFIN is free software: you can redistribute it and/or modify
6 // it under the terms of the GNU Lesser General Public License as published by
7 // the Free Software Foundation, either version 3 of the License, or
8 // (at your option) any later version.
9 //
10 // DOLFIN is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU Lesser General Public License for more details.
14 //
15 // You should have received a copy of the GNU Lesser General Public License
16 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
17 //
18 // Modified by August Johansson 2018
19 //
20 // First added: 2014-03-03
21 // Last changed: 2018-04-03
22 
23 #ifndef __MULTI_MESH_H
24 #define __MULTI_MESH_H
25 
26 #include <memory>
27 #include <vector>
28 #include <map>
29 #include <deque>
30 
31 #include <dolfin/common/Variable.h>
32 #include <dolfin/geometry/Point.h>
33 
34 namespace dolfin
35 {
36 
37  // Forward declarations
38  class Cell;
39  class Mesh;
40  class BoundaryMesh;
41  class BoundingBoxTree;
42  class SimplexQuadrature;
43 
49 
50  class MultiMesh : public Variable
51  {
52  public:
53 
55  typedef std::pair<std::vector<double>, std::vector<double> > quadrature_rule;
56 
58  typedef std::vector<Point> Simplex;
59 
61  typedef std::pair<std::vector<Simplex>, std::set<std::size_t>> Polyhedron;
62 
64  typedef std::vector<std::size_t> IncExcKey;
65 
67  MultiMesh();
68 
70  MultiMesh(std::vector<std::shared_ptr<const Mesh>> meshes,
71  std::size_t quadrature_order);
72 
73  //--- Convenience constructors ---
74 
76  MultiMesh(std::shared_ptr<const Mesh> mesh_0,
77  std::size_t quadrature_order);
78 
80  MultiMesh(std::shared_ptr<const Mesh> mesh_0,
81  std::shared_ptr<const Mesh> mesh_1,
82  std::size_t quadrature_order);
83 
85  MultiMesh(std::shared_ptr<const Mesh> mesh_0,
86  std::shared_ptr<const Mesh> mesh_1,
87  std::shared_ptr<const Mesh> mesh_2,
88  std::size_t quadrature_order);
89 
91  ~MultiMesh();
92 
98  std::size_t num_parts() const;
99 
109  std::shared_ptr<const Mesh> part(std::size_t i) const;
110 
122  const std::vector<unsigned int>& uncut_cells(std::size_t part) const;
123 
141  const std::vector<unsigned int> cut_cells(std::size_t part) const;
142 
156  const std::vector<unsigned int>& covered_cells(std::size_t part) const;
157 
165  void mark_covered(std::size_t part, const std::vector<unsigned int>& cells);
166 
178  const std::map<unsigned int,
179  std::vector<std::pair<std::size_t, unsigned int> > >&
180  collision_map_cut_cells(std::size_t part) const;
181 
194  const std::map<unsigned int, quadrature_rule >&
195  quadrature_rules_cut_cells(std::size_t part) const;
196 
211  const quadrature_rule
212  quadrature_rules_cut_cells(std::size_t part, unsigned int cell_index) const;
213 
228  const std::map<unsigned int, std::vector<quadrature_rule> >&
229  quadrature_rules_overlap(std::size_t part) const;
230 
237  // cell (unsigned int)
238  // The cell index
249  const std::vector<quadrature_rule>
250  quadrature_rules_overlap(std::size_t part, unsigned int cell) const;
251 
267  const std::map<unsigned int, std::vector<quadrature_rule> >&
268  quadrature_rules_interface(std::size_t part) const;
269 
270 
293  const std::vector<quadrature_rule>
294  quadrature_rules_interface(std::size_t part,
295  unsigned int cell_index) const;
296 
297 
315  const std::map<unsigned int, std::vector<std::vector<double> > >&
316  facet_normals(std::size_t part) const;
317 
327  std::shared_ptr<const BoundingBoxTree>
328  bounding_box_tree(std::size_t part) const;
329 
340  std::shared_ptr<const BoundingBoxTree>
341  bounding_box_tree_boundary(std::size_t part) const;
342 
348  void add(std::shared_ptr<const Mesh> mesh);
349 
351  void build(std::size_t quadrature_order=2);
352 
354  bool is_built() const { return _is_built; }
355 
357  void clear();
358 
361  {
362  Parameters p("multimesh");
363 
364  //p.add("quadrature_order", 1);
365  p.add("compress_volume_quadrature", false);
366  p.add("compress_interface_quadrature", false);
367 
368  return p;
369  }
370 
371  //--- The functions below are mainly useful for testing/debugging ---
372 
377  double compute_area() const;
378 
380  double compute_volume() const;
381 
383  std::string plot_matplotlib(double delta_z=1,
384  const std::string& filename="") const;
385 
389  void auto_cover(std::size_t p, const Point& point);
390 
391  private:
392 
393  // Flag for whether multimesh has been built
394  bool _is_built;
395 
396  // List of meshes
397  std::vector<std::shared_ptr<const Mesh> > _meshes;
398 
399  // List of boundary meshes
400  std::vector<std::shared_ptr<BoundaryMesh> > _boundary_meshes;
401 
402  // List of bounding box trees for meshes
403  std::vector<std::shared_ptr<BoundingBoxTree> > _trees;
404 
405  // List of bounding box trees for boundary meshes
406  std::vector<std::shared_ptr<BoundingBoxTree> > _boundary_trees;
407 
408  // Cell indices for all uncut cells for all parts. Access data by
409  //
410  // c = _uncut_cells[i][j]
411  //
412  // where
413  //
414  // c = cell index for an uncut cell
415  // i = the part (mesh) number
416  // j = the cell number (in the list of uncut cells)
417  std::vector<std::vector<unsigned int> > _uncut_cells;
418 
419  // Cell indices for all covered cells for all parts. Access data by
420  //
421  // c = _covered_cells[i][j]
422  //
423  // where
424  //
425  // c = cell index for a covered cell
426  // i = the part (mesh) number
427  // j = the cell number (in the list of covered cells)
428  std::vector<std::vector<unsigned int> > _covered_cells;
429 
430  // Developer note 1: The data structures _collision_map_cut_cells
431  // and _quadrature_rules_cut_cells may be changed from maps to
432  // vectors and indexed by the number of the cut cell (in the list
433  // of cut cells), instead of indexed by the local cell index as
434  // here, if we find that this is important for performance.
435  //
436  // Developer note 2: Quadrature points are naturally a part of a
437  // form (or a term in a form) and not a part of a mesh. However,
438  // for now we use a global (to the multimesh) quadrature rule for
439  // all cut cells, for simplicity.
440 
441  // Collision map for cut cells. Access data by
442  //
443  // c = _collision_map_cut_cells[i][j][k]
444  //
445  // where
446  //
447  // c.first = part number for the cutting mesh
448  // c.second = cell index for the cutting cell
449  // i = the part (mesh) number
450  // j = the cell number (local cell index
451  // k = the collision number (in the list of cutting cells)
452  std::vector<std::map<unsigned int,
453  std::vector<std::pair<std::size_t, unsigned int> > > >
454  _collision_maps_cut_cells;
455 
456  // Quadrature rules for cut cells. Access data by
457  //
458  // q = _quadrature_rules_cut_cells[i][j]
459  //
460  // where
461  //
462  // q.first = quadrature weights, array of length num_points
463  // q.second = quadrature points, flattened num_points x gdim array
464  // i = the part (mesh) number
465  // j = the cell number (local cell index)
466  std::vector<std::map<unsigned int, quadrature_rule> >
467  _quadrature_rules_cut_cells;
468 
469  // Quadrature rules for overlap. Access data by
470  //
471  // q = _quadrature_rules_overlap[i][j][k]
472  //
473  // where
474  //
475  // q.first = quadrature weights, array of length num_points
476  // q.second = quadrature points, flattened num_points x gdim array
477  // i = the part (mesh) number
478  // j = the cell number (local cell index)
479  // k = the collision number (in the list of cutting cells)
480  std::vector<std::map<unsigned int, std::vector<quadrature_rule> > >
481  _quadrature_rules_overlap;
482 
483  // Quadrature rules for interface. Access data by
484  //
485  // q = _quadrature_rules_interface[i][j][k]
486  //
487  // where
488  //
489  // q.first = quadrature weights, array of length num_points
490  // q.second = quadrature points, flattened num_points x gdim array
491  // i = the part (mesh) number
492  // j = the cell number (local cell index)
493  // k = the collision number (in the list of cutting cells)
494  std::vector<std::map<unsigned int, std::vector<quadrature_rule> > >
495  _quadrature_rules_interface;
496 
497  // Facet normals for interface. Access data by
498  //
499  // n = _facet_normals_interface[i][j][k][
500  //
501  // where
502  //
503  // n = a flattened array vector of facet normals, one point for
504  // each quadrature point
505  // i = the part (mesh) number
506  // j = the cell number (local cell index)
507  // k = the collision number (in the list of cutting cells)
508  std::vector<std::map<unsigned int, std::vector<std::vector<double> > > >
509  _facet_normals;
510 
511  // Build boundary meshes
512  void _build_boundary_meshes();
513 
514  // Build bounding box trees
515  void _build_bounding_box_trees();
516 
517  // Build collision maps
518  void _build_collision_maps();
519  //void _build_collision_maps_same_topology();
520  //void _build_collision_maps_different_topology();
521 
522  // Build quadrature rules for the cut cells
523  void _build_quadrature_rules_cut_cells(std::size_t quadrature_order);
524 
525  // Build quadrature rules for the overlap
526  void _build_quadrature_rules_overlap(std::size_t quadrature_order);
527 
528  // Build quadrature rules and normals for the interface
529  void _build_quadrature_rules_interface(std::size_t quadrature_order);
530 
531  // Help function to determine if interface intersection is
532  // (exactly) overlapped by a cutting cell
533  bool _is_overlapped_interface(std::vector<Point> simplex,
534  const Cell cut_cell,
535  Point simplex_normal) const;
536 
537  // Add quadrature rule for simplices in polyhedron. Returns the
538  // number of points generated for each simplex.
539  std::size_t _add_quadrature_rule(quadrature_rule& qr,
540  const SimplexQuadrature& sq,
541  const Simplex& simplex,
542  std::size_t gdim,
543  std::size_t quadrature_order,
544  double factor) const;
545 
546  // Add quadrature rule to existing quadrature rule (append dqr to
547  // qr). Returns number of points added.
548  std::size_t _add_quadrature_rule(quadrature_rule& qr,
549  const quadrature_rule& dqr,
550  std::size_t gdim,
551  double factor) const;
552 
553  // Append normal to list of normals npts times
554  void _add_normal(std::vector<double>& normals,
555  const Point& normal,
556  std::size_t npts,
557  std::size_t gdim) const;
558 
559  // Plot multimesh
560  void _plot() const;
561 
562  // Inclusion-exclusion for overlap
563  void _inclusion_exclusion_overlap
564  (std::vector<quadrature_rule>& qr,
565  const SimplexQuadrature& sq,
566  const std::vector<std::pair<std::size_t, Polyhedron> >& initial_polyhedra,
567  std::size_t tdim,
568  std::size_t gdim,
569  std::size_t quadrature_order) const;
570 
571  // Inclusion-exclusion for interface
572  void _inclusion_exclusion_interface
573  (quadrature_rule& qr,
574  std::vector<double>& normals,
575  const SimplexQuadrature& sq,
576  const Simplex& Eij,
577  const Point& facet_normal,
578  const std::vector<std::pair<std::size_t, Polyhedron> >& initial_polygons,
579  std::size_t tdim,
580  std::size_t gdim,
581  std::size_t quadrature_order) const;
582 
583  // Construct and return mapping from boundary facets to full mesh
584  std::vector<std::vector<std::pair<std::size_t, std::size_t> > >
585  _boundary_facets_to_full_mesh(std::size_t part) const;
586 
587  // Impose consistency of _cut_cells, so that only the cells with
588  // a nontrivial interface quadrature rule are classified as cut.
589  void _impose_cut_cell_consistency();
590 
591  // Remove quadrature rule if the sum of the weights is less than a
592  // tolerance
593  static void remove_quadrature_rule(quadrature_rule& qr,
594  double tolerance);
595  };
596 
597 
598 
599 }
600 
601 
602 #endif
void mark_covered(std::size_t part, const std::vector< unsigned int > &cells)
Definition: MultiMesh.cpp:153
Common base class for DOLFIN variables.
Definition: Variable.h:35
~MultiMesh()
Destructor.
Definition: MultiMesh.cpp:102
std::shared_ptr< const BoundingBoxTree > bounding_box_tree(std::size_t part) const
Definition: MultiMesh.cpp:236
bool is_built() const
Check whether multimesh has been built.
Definition: MultiMesh.h:354
Definition: adapt.h:29
std::pair< std::vector< double >, std::vector< double > > quadrature_rule
Structure storing a quadrature rule.
Definition: MultiMesh.h:55
void add(std::string key)
Definition: Parameters.h:128
Definition: Point.h:40
Definition: MultiMesh.h:50
const std::map< unsigned int, quadrature_rule > & quadrature_rules_cut_cells(std::size_t part) const
Definition: MultiMesh.cpp:180
const std::vector< unsigned int > cut_cells(std::size_t part) const
Definition: MultiMesh.cpp:126
std::shared_ptr< const BoundingBoxTree > bounding_box_tree_boundary(std::size_t part) const
Definition: MultiMesh.cpp:243
std::string plot_matplotlib(double delta_z=1, const std::string &filename="") const
Create matplotlib string to plot 2D multimesh (small meshes only)
Definition: MultiMesh.cpp:383
const std::vector< unsigned int > & covered_cells(std::size_t part) const
Definition: MultiMesh.cpp:146
std::shared_ptr< const Mesh > part(std::size_t i) const
Definition: MultiMesh.cpp:112
double compute_area() const
Definition: MultiMesh.cpp:307
std::vector< Point > Simplex
A simplex is a list of points.
Definition: MultiMesh.h:58
const std::vector< unsigned int > & uncut_cells(std::size_t part) const
Definition: MultiMesh.cpp:119
const std::map< unsigned int, std::vector< std::pair< std::size_t, unsigned int > > > & collision_map_cut_cells(std::size_t part) const
Definition: MultiMesh.cpp:173
A Cell is a MeshEntity of topological codimension 0.
Definition: Cell.h:42
std::size_t num_parts() const
Definition: MultiMesh.cpp:107
const std::map< unsigned int, std::vector< quadrature_rule > > & quadrature_rules_interface(std::size_t part) const
Definition: MultiMesh.cpp:213
Definition: Parameters.h:94
std::pair< std::vector< Simplex >, std::set< std::size_t > > Polyhedron
A polyhedron is a list of simplices and the part numbers.
Definition: MultiMesh.h:61
MultiMesh()
Create empty multimesh.
Definition: MultiMesh.cpp:45
This class defines quadrature rules for simplices.
Definition: SimplexQuadrature.h:36
void auto_cover(std::size_t p, const Point &point)
Definition: MultiMesh.cpp:457
std::vector< std::size_t > IncExcKey
Key to identify polyhedra.
Definition: MultiMesh.h:64
double compute_volume() const
Corresponding function for volume.
Definition: MultiMesh.cpp:350
const std::map< unsigned int, std::vector< quadrature_rule > > & quadrature_rules_overlap(std::size_t part) const
Definition: MultiMesh.cpp:197
const std::map< unsigned int, std::vector< std::vector< double > > > & facet_normals(std::size_t part) const
Definition: MultiMesh.cpp:229
void clear()
Clear multimesh.
Definition: MultiMesh.cpp:294
void add(std::shared_ptr< const Mesh > mesh)
Definition: MultiMesh.cpp:249
void build(std::size_t quadrature_order=2)
Build multimesh.
Definition: MultiMesh.cpp:256
static Parameters default_parameters()
Default parameter values.
Definition: MultiMesh.h:360