Flow123d  master-f44eb46
mh_dofhandler.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 mh_dofhandler.cc
15  * @brief
16  */
17 
18 #include "flow/mh_dofhandler.hh"
20 #include "mesh/mesh.h"
21 #include "mesh/partitioning.hh"
22 #include "mesh/accessors.hh"
23 #include "mesh/range_wrapper.hh"
24 #include "mesh/neighbours.h"
25 #include "system/index_types.hh"
26 #include "system/sys_profiler.hh"
27 
29 :
30 el_4_loc(nullptr),
31 row_4_el(nullptr),
32 side_id_4_loc(nullptr),
33 side_row_4_id(nullptr),
34 edge_4_loc(nullptr),
35 row_4_edge(nullptr),
36 edge_ds(nullptr),
37 el_ds(nullptr),
38 side_ds(nullptr)
39 {}
40 
42 {
43  delete edge_ds;
44  // delete el_ds;
45  delete side_ds;
46 
47  // xfree(el_4_loc);
48  delete [] row_4_el;
49  delete [] side_id_4_loc;
50  delete [] side_row_4_id;
51  delete [] edge_4_loc;
52  delete [] row_4_edge;
53 }
54 
55 
57  mesh_ = mesh;
58  elem_side_to_global.resize(mesh->n_elements() );
59  for (auto ele : mesh->elements_range()) elem_side_to_global[ ele.idx() ].resize(ele->n_sides());
60 
61  unsigned int i_side_global=0;
62  for (auto ele : mesh->elements_range()) {
63  for(unsigned int i_lside=0; i_lside < ele->n_sides(); i_lside++)
64  elem_side_to_global[ ele.idx() ][i_lside] = i_side_global++;
65  }
66 
68 
69 }
70 
71 
72 
73 
74 // ====================================================================================
75 // - compute optimal edge partitioning
76 // - compute appropriate partitioning of elements and sides
77 // - make arrays: *_id_4_loc and *_row_4_id to allow parallel assembly of the MH matrix
78 // ====================================================================================
80 
81  START_TIMER("prepare parallel");
82 
83  LongIdx *loc_part; // optimal (edge,el) partitioning (local chunk)
84  LongIdx *id_4_old; // map from old idx to ids (edge,el)
85  int loc_i;
86 
87  int e_idx;
88 
89 
90  //ierr = MPI_Barrier(PETSC_COMM_WORLD);
91  //ASSERT_PERMANENT(ierr == 0).error("Error in MPI_Barrier.");
92 
93  // row_4_el will be modified so we make a copy of the array from mesh
94  row_4_el = new LongIdx[mesh_->n_elements()];
97  el_ds = mesh_->get_el_ds();
98 
99  //optimal element part; loc. els. id-> new el. numbering
100  Distribution init_edge_ds(DistributionLocalized(), mesh_->n_edges(), PETSC_COMM_WORLD);
101  // partitioning of edges, edge belongs to the proc of his first element
102  // this is not optimal but simple
103  loc_part = new LongIdx[init_edge_ds.lsize()];
104  id_4_old = new LongIdx[mesh_->n_edges()];
105  {
106  loc_i = 0;
107  for(auto edg: mesh_->edge_range())
108  {
109  unsigned int i_edg = edg.idx();
110  // partition
111  e_idx = edg.side(0)->element().idx();
112  if (init_edge_ds.is_local(i_edg)) {
113  // find (new) proc of the first element of the edge
114  loc_part[loc_i++] = el_ds->get_proc(row_4_el[e_idx]);
115  }
116  // id array
117  id_4_old[i_edg] = i_edg;
118  }
119  }
120  Partitioning::id_maps(mesh_->n_edges(), id_4_old, init_edge_ds, loc_part, edge_ds, edge_4_loc, row_4_edge);
121  delete[] loc_part;
122  delete[] id_4_old;
123 
124  // create map from mesh global edge id to new local edge id
125  unsigned int loc_edge_idx=0;
126  for (unsigned int i_el_loc = 0; i_el_loc < el_ds->lsize(); i_el_loc++) {
127  auto ele = mesh_->element_accessor( el_4_loc[i_el_loc] );
128  for (unsigned int i = 0; i < ele->n_sides(); i++) {
129  unsigned int mesh_edge_idx= ele.side(i)->edge_idx();
130  if ( edge_new_local_4_mesh_idx_.count(mesh_edge_idx) == 0 )
131  // new local edge
132  edge_new_local_4_mesh_idx_[mesh_edge_idx] = loc_edge_idx++;
133  }
134  }
135 
136  //optimal side part; loc. sides; id-> new side numbering
137  Distribution init_side_ds(DistributionBlock(), mesh_->n_sides(), PETSC_COMM_WORLD);
138  // partitioning of sides follows elements
139  loc_part = new LongIdx[init_side_ds.lsize()];
140  id_4_old = new LongIdx[mesh_->n_sides()];
141  {
142  int is = 0;
143  loc_i = 0;
144  for (auto ele : mesh_->elements_range())
145  for(SideIter side = ele.side(0); side->side_idx() < ele->n_sides(); ++side) {
146  // partition
147  if (init_side_ds.is_local(is)) {
148  // find (new) proc of the element of the side
149  loc_part[loc_i++] = el_ds->get_proc(
150  row_4_el[ side->element().idx() ]);
151  }
152  // id array
153  id_4_old[is++] = side_dof( side );
154  }
155  }
156 
157  Partitioning::id_maps(mesh_->n_sides(), id_4_old, init_side_ds, loc_part, side_ds,
159  delete [] loc_part;
160  delete [] id_4_old;
161 
162  // convert row_4_id arrays from separate numberings to global numbering of rows
164 
165 }
166 
167 // ========================================================================
168 // to finish row_4_id arrays we have to convert individual numberings of
169 // sides/els/edges to whole numbering of rows. To this end we count shifts
170 // for sides/els/edges on each proc and then we apply them on row_4_id
171 // arrays.
172 // we employ macros to avoid code redundancy
173 // =======================================================================
175  int i, shift;
176  int np = edge_ds->np();
177  int edge_shift[np], el_shift[np], side_shift[np];
178 
179  int edge_n_id = mesh_->n_edges(),
180  el_n_id = mesh_->n_elements(),
181  side_n_id = mesh_->n_sides();
182 
183  // compute shifts on every proc
184  shift = 0; // in this var. we count new starts of arrays chunks
185  for (i = 0; i < np; i++) {
186  side_shift[i] = shift - (side_ds->begin(i)); // subtract actual start of the chunk
187  shift += side_ds->lsize(i);
188  el_shift[i] = shift - (el_ds->begin(i));
189  shift += el_ds->lsize(i);
190  edge_shift[i] = shift - (edge_ds->begin(i));
191  shift += edge_ds->lsize(i);
192  }
193  // apply shifts
194  for (i = 0; i < side_n_id; i++) {
195  LongIdx &what = side_row_4_id[i];
196  if (what >= 0)
197  what += side_shift[side_ds->get_proc(what)];
198  }
199  for (i = 0; i < el_n_id; i++) {
200  LongIdx &what = row_4_el[i];
201  if (what >= 0)
202  what += el_shift[el_ds->get_proc(what)];
203 
204  }
205  for (i = 0; i < edge_n_id; i++) {
206  LongIdx &what = row_4_edge[i];
207  if (what >= 0)
208  what += edge_shift[edge_ds->get_proc(what)];
209  }
210 }
211 
212 
213 unsigned int MH_DofHandler::side_dof(const SideIter side) const {
214  return elem_side_to_global[ side->element().idx() ][ side->side_idx() ];
215 }
216 
217 
218 void MH_DofHandler::set_solution( double time, double * solution) {
219  ASSERT_PTR( solution ).error("Empty solution.\n");
220  mh_solution = solution;
221  time_ = time;
222 }
223 
224 /// temporary replacement for DofHandler accessor, flux through given side
225 double MH_DofHandler::side_flux(const Side &side) const {
226  return mh_solution[ elem_side_to_global[ side.element().idx() ][ side.side_idx() ] ];
227 }
228 
229 /// temporary replacement for DofHandler accessor, scalar (pressure) on edge of the side
230 double MH_DofHandler::side_scalar(const Side &side) const {
231  unsigned int i_edg = side.edge_idx();
232  return mh_solution[ side.mesh()->n_sides() + side.mesh()->n_elements() + i_edg ];
233 }
234 
235 
237  return mh_solution[ mesh_->n_sides() + ele.idx() ];
238 }
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
Definition: asserts.hh:341
bool is_local(unsigned int idx) const
identify local index
unsigned int begin(int proc) const
get starting local index
unsigned int lsize(int proc) const
get local size
unsigned int get_proc(unsigned int idx) const
get processor of the given index
unsigned int np() const
get num of processors
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
double element_scalar(ElementAccessor< 3 > &ele) const
temporary replacement for DofHandler accessor, scalar (pressure) on element
Distribution * el_ds
std::unordered_map< unsigned int, unsigned int > edge_new_local_4_mesh_idx_
vector< vector< unsigned int > > elem_side_to_global
LongIdx * edge_4_loc
LongIdx * row_4_el
void set_solution(double time, double *solution)
void make_row_numberings()
Distribution * edge_ds
LongIdx * row_4_edge
double side_flux(const Side &side) const
temporary replacement for DofHandler accessor, flux through given side
double * mh_solution
void prepare_parallel()
double side_scalar(const Side &side) const
temporary replacement for DofHandler accessor, scalar (pressure) on edge of the side
LongIdx * side_id_4_loc
LongIdx * el_4_loc
unsigned int side_dof(const SideIter side) const
void reinit(Mesh *mesh)
Distribution * side_ds
LongIdx * side_row_4_id
LongIdx * get_el_4_loc() const
Definition: mesh.h:122
Range< Edge > edge_range() const
Return range of edges.
Definition: mesh.cc:125
unsigned int n_edges() const
Definition: mesh.h:114
LongIdx * get_row_4_el() const
Definition: mesh.h:125
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
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
unsigned int n_sides() const
Definition: mesh.cc:308
void id_maps(int n_ids, LongIdx *id_4_old, Distribution *&new_ds, LongIdx *&id_4_loc, LongIdx *&new_4_id)
unsigned int side_idx() const
Returns local index of the side on the element.
Definition: accessors.hh:444
const MeshBase * mesh() const
Returns pointer to the mesh.
Definition: accessors.hh:440
ElementAccessor< 3 > element() const
Returns iterator to the element of the side.
unsigned int edge_idx() const
Returns global index of the edge connected to the side.
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
Definition: index_types.hh:24
Implementation of range helper class.
#define START_TIMER(tag)
Starts a timer with specified tag.