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