DOLFIN
DOLFIN C++ interface
Loading...
Searching...
No Matches
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
34namespace 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>
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
A Cell is a MeshEntity of topological codimension 0.
Definition Cell.h:43
Definition MultiMesh.h:51
const std::vector< unsigned int > cut_cells(std::size_t part) const
Definition MultiMesh.cpp:126
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
double compute_volume() const
Corresponding function for volume.
Definition MultiMesh.cpp:350
std::shared_ptr< const Mesh > part(std::size_t i) const
Definition MultiMesh.cpp:112
bool is_built() const
Check whether multimesh has been built.
Definition MultiMesh.h:354
void build(std::size_t quadrature_order=2)
Build multimesh.
Definition MultiMesh.cpp:256
std::vector< Point > Simplex
A simplex is a list of points.
Definition MultiMesh.h:58
double compute_area() const
Definition MultiMesh.cpp:307
std::pair< std::vector< double >, std::vector< double > > quadrature_rule
Structure storing a quadrature rule.
Definition MultiMesh.h:55
~MultiMesh()
Destructor.
Definition MultiMesh.cpp:102
static Parameters default_parameters()
Default parameter values.
Definition MultiMesh.h:360
std::size_t num_parts() const
Definition MultiMesh.cpp:107
void add(std::shared_ptr< const Mesh > mesh)
Definition MultiMesh.cpp:249
const std::vector< unsigned int > & uncut_cells(std::size_t part) const
Definition MultiMesh.cpp:119
std::shared_ptr< const BoundingBoxTree > bounding_box_tree_boundary(std::size_t part) const
Definition MultiMesh.cpp:243
void mark_covered(std::size_t part, const std::vector< unsigned int > &cells)
Definition MultiMesh.cpp:153
const std::vector< unsigned int > & covered_cells(std::size_t part) const
Definition MultiMesh.cpp:146
MultiMesh()
Create empty multimesh.
Definition MultiMesh.cpp:45
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
const std::map< unsigned int, std::vector< quadrature_rule > > & quadrature_rules_interface(std::size_t part) const
Definition MultiMesh.cpp:213
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::map< unsigned int, std::vector< std::vector< double > > > & facet_normals(std::size_t part) const
Definition MultiMesh.cpp:229
const std::map< unsigned int, quadrature_rule > & quadrature_rules_cut_cells(std::size_t part) const
Definition MultiMesh.cpp:180
const std::map< unsigned int, std::vector< quadrature_rule > > & quadrature_rules_overlap(std::size_t part) const
Definition MultiMesh.cpp:197
std::vector< std::size_t > IncExcKey
Key to identify polyhedra.
Definition MultiMesh.h:64
std::shared_ptr< const BoundingBoxTree > bounding_box_tree(std::size_t part) const
Definition MultiMesh.cpp:236
void auto_cover(std::size_t p, const Point &point)
Definition MultiMesh.cpp:457
void clear()
Clear multimesh.
Definition MultiMesh.cpp:294
Definition Parameters.h:95
void add(std::string key)
Definition Parameters.h:128
Definition Point.h:41
This class defines quadrature rules for simplices.
Definition SimplexQuadrature.h:37
Common base class for DOLFIN variables.
Definition Variable.h:36
Definition adapt.h:30