Flow123d  master-f44eb46
bc_mesh.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file bc_mesh.cc
15  * @ingroup mesh
16  * @brief Mesh construction
17  */
18 
19 
20 #include "system/index_types.hh"
21 #include "mesh/bc_mesh.hh"
22 #include "mesh/accessors.hh"
23 #include "mesh/partitioning.hh"
24 #include "mesh/neighbours.h"
25 #include "mesh/range_wrapper.hh"
26 #include "la/distribution.hh"
27 
28 
29 
30 BCMesh::BCMesh(Mesh* parent_mesh)
31 : parent_mesh_(parent_mesh),
32  local_part_(nullptr)
33 {
35  this->init_element_vector(1); // necessary to allocate 1 element due to elements_range() which otherwise fails for empty BCMesh
36  this->init_node_vector(0);
37  this->nodes_ = parent_mesh_->nodes_;
39 }
40 
41 
43 {
44  if (local_part_!=nullptr) delete local_part_;
45 }
46 
48 {
49  vector<LongIdx> loc_el_ids;
50 
51  // loop local parent elements and their sides to find boundary elements
52  for (unsigned int iel = 0; iel<parent_mesh_->get_el_ds()->lsize(); iel++)
53  {
54  auto parent_el = parent_mesh_->element_accessor(parent_mesh_->get_el_4_loc()[iel]);
55  for (unsigned int sid=0; sid<parent_el->n_sides(); sid++)
56  {
57  if (parent_el.side(sid)->is_boundary())
58  {
59  loc_el_ids.push_back(parent_el.side(sid)->cond().element_accessor().idx());
60  }
61  }
62  }
63  this->el_ds = new Distribution(loc_el_ids.size(), PETSC_COMM_WORLD);
64 
65  this->el_4_loc = new LongIdx[loc_el_ids.size()];
66  for (unsigned int i=0; i<loc_el_ids.size(); i++) this->el_4_loc[i] = loc_el_ids[i];
67 
68  this->row_4_el = new LongIdx[n_elements()];
69  vector<LongIdx> row_4_loc_el(n_elements(), 0);
70  for (unsigned int i=0; i<loc_el_ids.size(); i++)
71  row_4_loc_el[loc_el_ids[i]] = i + this->el_ds->begin();
72  MPI_Allreduce(row_4_loc_el.data(), this->row_4_el, n_elements(), MPI_LONG_IDX, MPI_MAX, PETSC_COMM_WORLD);
73 
75 }
76 
77 
79  return parent_mesh_->get_part();
80 }
81 
83  if (local_part_ == nullptr) {
84  local_part_ = new LongIdx[this->n_elements()];
85  unsigned int bc_ele_idx;
86  for (auto ele : parent_mesh_->elements_range())
87  if (ele->boundary_idx_ != NULL)
88  for (unsigned int i=0; i<ele->n_sides(); ++i)
89  if ((int)ele->boundary_idx_[i] != -1) {
90  bc_ele_idx = parent_mesh_->boundary_[ ele->boundary_idx_[i] ].bc_ele_idx_;
91  local_part_[bc_ele_idx] = parent_mesh_->get_local_part()[ele.idx()];
92  }
93  }
94  return local_part_;
95 }
96 
97 
98 std::shared_ptr<EquivalentMeshMap> BCMesh::check_compatible_mesh( Mesh & input_mesh) {
99  return parent_mesh_->check_compatible_mesh(input_mesh);
100 }
101 
102 
103 Boundary BCMesh::boundary(unsigned int) const
104 {
105  ASSERT_PERMANENT( false );
106  return Boundary();
107 }
108 
109 
111 {
112  Neighbour neighbour;
113  EdgeData *edg = nullptr;
114  unsigned int ngh_element_idx;
115  unsigned int last_edge_idx = Mesh::undef_idx;
116 
117  neighbour.mesh_ = this;
118 
120 
121  // pointers to created edges
122  edges.resize(0); // be sure that edges are empty
123 
125  vector<unsigned int> intersection_list; // list of elements in intersection of node element lists
126 
127  // Now we go through all element sides and create edges and neighbours
128  for (auto e : this->elements_range()) {
129  for (unsigned int s=0; s<e->n_sides(); s++)
130  {
131  // skip sides that were already found
132  if (e->edge_idx(s) != Mesh::undef_idx) continue;
133 
134  // Find all elements that share this side.
135  side_nodes.resize(e.side(s)->n_nodes());
136  for (unsigned n=0; n<e.side(s)->n_nodes(); n++) side_nodes[n] = e.side(s)->node(n).idx();
137  intersect_element_lists(side_nodes, intersection_list);
138 
139  bool is_neighbour = find_lower_dim_element(intersection_list, e->dim(), ngh_element_idx);
140 
141  if (is_neighbour) { // edge connects elements of different dimensions
142  // Initialize for the neighbour case.
143  neighbour.elem_idx_ = ngh_element_idx;
144  } else { // edge connects only elements of the same dimension
145  // Initialize for the edge case.
146  last_edge_idx=edges.size();
147  edges.resize(last_edge_idx+1);
148  edg = &( edges.back() );
149  edg->n_sides = 0;
150  edg->side_ = new struct SideIter[ intersection_list.size() ];
151  if (e->dim() > 0)
152  if (intersection_list.size() > max_edge_sides_[e->dim()-1])
153  max_edge_sides_[e->dim()-1] = intersection_list.size();
154 
155  if (intersection_list.size() <= 1) {
156  // outer edge, create boundary object as well
157  edg->n_sides=1;
158  edg->side_[0] = e.side(s);
159  element_vec_[e.idx()].edge_idx_[s] = last_edge_idx;
160 
161  continue; // next side of element e
162  }
163  }
164 
165  // go through the elements connected to the edge or neighbour
166  // setup neigbour or edge
167  for( vector<unsigned int>::iterator isect = intersection_list.begin(); isect!=intersection_list.end(); ++isect) {
168  ElementAccessor<3> elem = this->element_accessor(*isect);
169  for (unsigned int ecs=0; ecs<elem->n_sides(); ecs++) {
170  if (elem->edge_idx(ecs) != Mesh::undef_idx) continue; // ??? This should not happen.
171  SideIter si = elem.side(ecs);
172  if ( same_sides( si, side_nodes) ) {
173  if (is_neighbour) {
174  // create a new edge and neighbour for this side, and element to the edge
175  last_edge_idx=edges.size();
176  edges.resize(last_edge_idx+1);
177  edg = &( edges.back() );
178  edg->n_sides = 1;
179  edg->side_ = new struct SideIter[1];
180  edg->side_[0] = si;
181  element_vec_[elem.idx()].edge_idx_[ecs] = last_edge_idx;
182 
183  neighbour.edge_idx_ = last_edge_idx;
184 
185  vb_neighbours_.push_back(neighbour); // copy neighbour with this edge setting
186  } else {
187  // connect the side to the edge, and side to the edge
188  ASSERT_PTR(edg);
189  edg->side_[ edg->n_sides++ ] = si;
190  ASSERT(last_edge_idx != Mesh::undef_idx);
191  element_vec_[elem.idx()].edge_idx_[ecs] = last_edge_idx;
192  }
193  break; // next element from intersection list
194  }
195  } // search for side of other connected element
196  } // connected elements
197 
198  if (! is_neighbour)
199  ASSERT_EQ( (unsigned int) edg->n_sides, intersection_list.size())(e.input_id())(s).error("Missing edge sides.");
200  } // for element sides
201  } // for elements
202 
203  MessageOut().fmt( "Created {} edges and {} neighbours on boundary mesh.\n", edges.size(), vb_neighbours_.size() );
204 }
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
Definition: asserts.hh:341
std::shared_ptr< EquivalentMeshMap > check_compatible_mesh(Mesh &input_mesh) override
Overwrite Mesh::check_compatible_mesh()
Definition: bc_mesh.cc:98
void make_neighbours_and_edges()
Definition: bc_mesh.cc:110
void init_distribution()
setup distribution of elements and related vectors
Definition: bc_mesh.cc:47
BCMesh(Mesh *parent_mesh)
Definition: bc_mesh.cc:30
~BCMesh() override
Destructor.
Definition: bc_mesh.cc:42
Partitioning * get_part() override
Overwrite Mesh::get_part()
Definition: bc_mesh.cc:78
const LongIdx * get_local_part() override
Overwrite Mesh::get_local_part()
Definition: bc_mesh.cc:82
Mesh * parent_mesh_
Pointer to parent (bulk) mesh.
Definition: bc_mesh.hh:74
Boundary boundary(unsigned int) const override
Definition: bc_mesh.cc:103
LongIdx * local_part_
Distribution of boundary elements to processors.
Definition: bc_mesh.hh:77
unsigned int begin(int proc) const
get starting local index
SideIter * side_
Definition: mesh_data.hh:32
unsigned int n_sides
Definition: mesh_data.hh:29
SideIter side(const unsigned int loc_index)
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
Definition: accessors.hh:215
unsigned int edge_idx(unsigned int edg_idx) const
Return edge_idx of given index.
Definition: elements.h:135
unsigned int n_sides() const
Definition: elements.h:131
vector< Element > element_vec_
Definition: mesh.h:295
shared_ptr< Armor::Array< double > > nodes_
Definition: mesh.h:312
vector< vector< vector< unsigned int > > > side_nodes
Definition: mesh.h:254
void intersect_element_lists(vector< unsigned int > const &nodes_list, vector< unsigned int > &intersection_element_list)
Definition: mesh.cc:552
LongIdx * get_el_4_loc() const
Definition: mesh.h:122
bool same_sides(const SideIter &si, vector< unsigned int > &side_nodes)
Definition: mesh.cc:606
Distribution * el_ds
Parallel distribution of elements.
Definition: mesh.h:334
std::shared_ptr< RegionDB > region_db_
Definition: mesh.h:343
shared_ptr< BidirectionalMap< int > > node_ids_
Maps node ids to indexes into vector node_vec_.
Definition: mesh.h:315
void init_element_vector(unsigned int size)
Initialize element_vec_, set size and reset counters of boundary and bulk elements.
Definition: mesh.cc:1142
unsigned int n_nodes() const
Definition: mesh.h:171
LongIdx * el_4_loc
Index set assigning to local element index its global index.
Definition: mesh.h:332
bool find_lower_dim_element(vector< unsigned int > &element_list, unsigned int dim, unsigned int &element_idx)
Definition: mesh.cc:584
void create_node_element_lists()
Definition: mesh.cc:539
vector< Neighbour > vb_neighbours_
Vector of compatible neighbourings.
Definition: mesh.h:304
LongIdx * row_4_el
Definition: mesh.h:330
std::vector< EdgeData > edges
Vector of MH edges, this should not be part of the geometrical mesh.
Definition: mesh.h:301
void init_node_vector(unsigned int size)
Initialize node_vec_, set size.
Definition: mesh.cc:1152
Distribution * get_el_ds() const
Definition: mesh.h:119
unsigned int n_elements() const
Definition: mesh.h:111
Range< ElementAccessor< 3 > > elements_range() const
Returns range of mesh elements.
Definition: mesh.cc:1174
static const unsigned int undef_idx
Definition: mesh.h:105
unsigned int max_edge_sides_[3]
Maximal number of sides per one edge in the actual mesh (set in make_neighbours_and_edges()).
Definition: mesh.h:307
ElementAccessor< 3 > element_accessor(unsigned int idx) const
Create and return ElementAccessor to element of given idx.
Definition: mesh.cc:866
Definition: mesh.h:362
const LongIdx * get_local_part() override
Definition: mesh.cc:338
Partitioning * get_part() override
Definition: mesh.cc:334
vector< BoundaryData > boundary_
Definition: mesh.h:506
std::shared_ptr< EquivalentMeshMap > check_compatible_mesh(Mesh &input_mesh) override
Definition: mesh.cc:919
MeshBase * mesh_
Pointer to Mesh to which belonged.
Definition: neighbours.h:136
unsigned int elem_idx_
Index of element in Mesh::element_vec_.
Definition: neighbours.h:137
unsigned int edge_idx_
Index of Edge in Mesh.
Definition: neighbours.h:138
Class for the mesh partitioning. This should provide:
Definition: partitioning.hh:52
Support classes for parallel programing.
#define MPI_LONG_IDX
Definition: index_types.hh:30
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
Definition: index_types.hh:24
#define MessageOut()
Macro defining 'message' record of log.
Definition: logger.hh:275
#define MPI_Allreduce(sendbuf, recvbuf, count, datatype, op, comm)
Definition: mpi.h:612
#define MPI_MAX
Definition: mpi.h:197
Implementation of range helper class.