Flow123d  master-f44eb46
output_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 output_mesh.cc
15  * @brief Classes for auxiliary output mesh.
16  */
17 
18 #include "system/index_types.hh"
19 #include "output_mesh.hh"
20 #include "output_element.hh"
21 #include "mesh/mesh.h"
22 #include "mesh/ref_element.hh"
23 #include "mesh/accessors.hh"
24 #include "mesh/node_accessor.hh"
25 #include "mesh/range_wrapper.hh"
26 #include "la/distribution.hh"
27 
28 
29 namespace IT=Input::Type;
30 
32  return IT::Record("OutputMesh", "Parameters of the refined output mesh. [Not impemented]")
33  .declare_key("max_level", IT::Integer(1,20),IT::Default("3"),
34  "Maximal level of refinement of the output mesh.")
35  .declare_key("refine_by_error", IT::Bool(), IT::Default("false"),
36  "Set true for using ``error_control_field``. Set false for global uniform refinement to max_level.")
37  .declare_key("error_control_field",IT::String(), IT::Default::optional(),
38  "Name of an output field, according to which the output mesh will be refined. The field must be a SCALAR one.")
39  .declare_key("refinement_error_tolerance",IT::Double(0.0), IT::Default("0.01"),
40  "Tolerance for element refinement by error. If tolerance is reached, refinement is stopped."
41  "Relative difference between error control field and its linear approximation on element is computed"
42  "and compared with tolerance.")
43  .close();
44 }
45 
47 :
48  orig_mesh_(&mesh),
49  max_level_(0),
50  refine_by_error_(false),
51  refinement_error_tolerance_(0.0),
52  el_ds_(nullptr),
53  node_ds_(nullptr)
54 {
55 }
56 
57 
59 :
60  input_record_(in_rec),
61  orig_mesh_(&mesh),
62  max_level_(input_record_.val<int>("max_level")),
63  refine_by_error_(input_record_.val<bool>("refine_by_error")),
64  refinement_error_tolerance_(input_record_.val<double>("refinement_error_tolerance")),
65  el_ds_(nullptr),
66  node_ds_(nullptr)
67 {
68 }
69 
71 {
72  // Refined mesh creates own special distributions and local to global maps and needs destroy these objects.
73  if ( (el_ds_!=nullptr) && (mesh_type_ == MeshType::refined)) {
74  delete[] el_4_loc_;
75  delete[] node_4_loc_;
76  delete el_ds_;
77  delete node_ds_;
78  }
79 }
80 
82 {
84 // ASSERT(offsets_->n_values() > 0);
85  return OutputElementIterator(OutputElement(0, shared_from_this()));
86 }
87 
89 {
91 // ASSERT(offsets_->n_values() > 0);
92  return OutputElementIterator(OutputElement(offsets_->n_values()-1, shared_from_this()));
93 }
94 
95 
97 {
98  error_control_field_func_ = error_control_field_func;
99 }
100 
102 {
104  return offsets_->n_values()-1;
105 }
106 
108 {
110  return nodes_->n_values();
111 }
112 
114 {
115  unsigned int elm_idx[1];
116  unsigned int node_idx[1];
117  unsigned int region_idx[1];
118  int partition[1];
119  elem_ids_ = std::make_shared< ElementDataCache<unsigned int> >("elements_ids", (unsigned int)1, this->n_elements());
120  node_ids_ = std::make_shared< ElementDataCache<unsigned int> >("node_ids", (unsigned int)1, this->n_nodes());
121  region_ids_ = std::make_shared< ElementDataCache<unsigned int> >("region_ids", (unsigned int)1, this->n_elements());
122  partitions_ = std::make_shared< ElementDataCache<int> >("partitions", (unsigned int)1, this->n_elements());
123  OutputElementIterator it = this->begin();
124  for (unsigned int i = 0; i < this->n_elements(); ++i, ++it) {
125  if (mesh_type_ == MeshType::orig) elm_idx[0] = orig_mesh_->find_elem_id(it->idx());
126  else elm_idx[0] = it->idx() + 1; // GMSH starts at 1
127  elem_ids_->store_value( i, elm_idx );
128 
129  region_idx[0] = orig_mesh_->element_accessor(it->idx()).region().id();
130  region_ids_->store_value( i, region_idx );
131 
132  partition[0] = orig_mesh_->element_accessor(it->idx())->pid();
133  partitions_->store_value( i, partition );
134 
135  std::vector< unsigned int > node_list = it->node_list();
136  for (unsigned int j = 0; j < it->n_nodes(); ++j) {
137  if (mesh_type_ == MeshType::orig) node_idx[0] = orig_mesh_->find_node_id(node_list[j]);
138  else node_idx[0] = node_list[j] + 1; // GMSH starts at 1
139  node_ids_->store_value( node_list[j], node_idx );
140  }
141  }
142 }
143 
144 
146 {
147  return (nodes_ && connectivity_ && offsets_);
148 }
149 
150 
152 {
153  ASSERT( !is_created() ).error("Multiple initialization of OutputMesh!\n");
154 
155  DebugOut() << "Create output submesh containing only local elements.";
156 
157  unsigned int ele_id = 0,
158  offset = 0, // offset of node indices of element in node vector
159  coord_id = 0, // coordinate id in node vector
160  conn_id = 0; // index to connectivity vector
161  ElementAccessor<3> elm;
162 
168 
169  const unsigned int n_local_elements = el_ds_->lsize();
170  unsigned int n_nodes = node_ds_->end( node_ds_->np()-1 );
171  std::vector<unsigned int> local_nodes_map(n_nodes, undef_idx); // map global to local ids of nodes
172  for (unsigned int i=0; i<n_local_nodes_; ++i) local_nodes_map[ node_4_loc_[i] ] = i;
173 
174  orig_element_indices_ = std::make_shared<std::vector<unsigned int>>(n_local_elements);
175  offsets_ = std::make_shared<ElementDataCache<unsigned int>>("offsets", (unsigned int)ElementDataCacheBase::N_SCALAR, n_local_elements+1);
176  auto &offset_vec = *( offsets_->get_data().get() );
177 
178  offset_vec[0] = 0;
179  loc_4_el_ = new LongIdx [ el_ds_->size() ];
180  for (unsigned int loc_el = 0; loc_el < el_ds_->size(); loc_el++) loc_4_el_[loc_el] = -1;
181  for (unsigned int loc_el = 0; loc_el < n_local_elements; loc_el++) {
182  loc_4_el_[ el_4_loc_[loc_el] ] = loc_el;
183  elm = orig_mesh_->element_accessor( el_4_loc_[loc_el] );
184  // increase offset by number of nodes of the simplicial element
185  offset += elm->dim() + 1;
186  offset_vec[ele_id+1] = offset;
187  (*orig_element_indices_)[ele_id] = el_4_loc_[loc_el];
188  ele_id++;
189  }
190 
191  connectivity_ = std::make_shared<ElementDataCache<unsigned int>>("connectivity", (unsigned int)ElementDataCacheBase::N_SCALAR,
192  offset_vec[offset_vec.size()-1]);
193  auto &connectivity_vec = *( connectivity_->get_data().get() );
194  for (unsigned int loc_el = 0; loc_el < n_local_elements; loc_el++) {
195  elm = orig_mesh_->element_accessor( el_4_loc_[loc_el] );
196  for (unsigned int li=0; li<elm->n_nodes(); li++) {
197  ASSERT(local_nodes_map[ elm.node(li).idx() ] != undef_idx)(elm.node(li).idx()).error("Undefined global to local node index!");
198  connectivity_vec[conn_id++] = local_nodes_map[ elm.node(li).idx() ];
199  }
200  }
201 
202  // set coords of nodes
203  nodes_ = std::make_shared<ElementDataCache<double>>("", (unsigned int)ElementDataCacheBase::N_VECTOR, n_local_nodes_);
204  auto &node_vec = *( nodes_->get_data() );
205  for(unsigned int i_node=0; i_node<local_nodes_map.size(); ++i_node) {
206  if (local_nodes_map[i_node]==undef_idx) continue; // skip element if it is not local
207  auto node = *orig_mesh_->node(i_node);
208  coord_id = 3*local_nodes_map[i_node]; // id of first coordinates in node_vec
209  node_vec[coord_id++] = node[0];
210  node_vec[coord_id++] = node[1];
211  node_vec[coord_id] = node[2];
212  }
213 }
214 
215 
216 
218 {
219  std::shared_ptr<ElementDataCache<unsigned int>> global_offsets; // needs for creating serial nodes and connectivity caches on zero process
220  auto elems_n_nodes = get_elems_n_nodes(); // collects number of nodes on each elements (for fill master_mesh_->offsets_)
221  int rank = el_ds_->myp();
222 
223  if (rank==0) {
224  // create serial output mesh, fill offsets cache and orig_element_indices vector
225  unsigned int n_elems = el_ds_->end( el_ds_->np()-1 );
226  master_mesh_ = this->construct_mesh();
227  master_mesh_->orig_element_indices_ = std::make_shared<std::vector<unsigned int>>(n_elems);
228  master_mesh_->offsets_ = std::make_shared<ElementDataCache<unsigned int>>(
229  "offsets", (unsigned int) ElementDataCacheBase::N_SCALAR, n_elems+1);
230  auto &offsets_vec = *( master_mesh_->offsets_->get_data().get() );
231  auto &elems_n_nodes_vec = *( elems_n_nodes->get_data().get() );
232  unsigned int offset=0;
233  offsets_vec[0] = 0;
234  for (unsigned int i=0; i<n_elems; ++i) {
235  offset += elems_n_nodes_vec[i];
236  offsets_vec[i+1] = offset;
237  (*master_mesh_->orig_element_indices_)[i] = i;
238  }
239  global_offsets = master_mesh_->offsets_;
240  }
241 
242  // collects serial caches
243  std::shared_ptr<ElementDataCache<double>> serial_nodes_cache = make_serial_nodes_cache(global_offsets);
244  std::shared_ptr<ElementDataCache<unsigned int>> serial_connectivity_cache = make_serial_connectivity_cache(global_offsets);
245 
246  if (rank==0) {
247  // set serial output mesh caches
248  master_mesh_->connectivity_ = serial_connectivity_cache;
249  master_mesh_->nodes_ = serial_nodes_cache;
250 
251  master_mesh_->mesh_type_ = this->mesh_type_;
252  }
253 }
254 
255 
256 std::shared_ptr<ElementDataCache<unsigned int>> OutputMeshBase::get_elems_n_nodes()
257 {
258  // Compute (locally) number of nodes of each elements
259  ElementDataCache<unsigned int> local_elems_n_nodes("elems_n_nodes", (unsigned int) ElementDataCacheBase::N_SCALAR, offsets_->n_values()-1);
260  auto &local_elems_n_nodes_vec = *( local_elems_n_nodes.get_data().get() );
261  auto &offset_vec = *( offsets_->get_data().get() );
262  for (unsigned int i=0; i<local_elems_n_nodes.n_values(); ++i) local_elems_n_nodes_vec[i] = offset_vec[i+1] - offset_vec[i];
263 
264  // Collect data, set on zero process
265  std::shared_ptr<ElementDataCache<unsigned int>> global_elems_n_nodes;
266  auto gather_cache = local_elems_n_nodes.gather(el_ds_, el_4_loc_);
267  if (el_ds_->myp()==0) global_elems_n_nodes = std::dynamic_pointer_cast< ElementDataCache<unsigned int> >(gather_cache);
268  return global_elems_n_nodes;
269 }
270 
271 
272 
273 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
274 
275 
277 : OutputMeshBase(mesh)
278 {
279  this->mesh_type_ = MeshType::orig;
280 }
281 
283 : OutputMeshBase(mesh, in_rec)
284 {
285  this->mesh_type_ = MeshType::orig;
286 }
287 
288 
290 {
291 }
292 
293 
295 {
296  ASSERT_PERMANENT(0).error("Not implemented yet.");
297 }
298 
300 {
301  ASSERT_PERMANENT(0).error("Not implemented yet.");
302  return false;
303 }
304 
305 
306 std::shared_ptr<OutputMeshBase> OutputMesh::construct_mesh()
307 {
308  return std::make_shared<OutputMesh>(*orig_mesh_);
309 }
310 
311 
312 std::shared_ptr<ElementDataCache<double>> OutputMesh::make_serial_nodes_cache(FMT_UNUSED std::shared_ptr<ElementDataCache<unsigned int>> global_offsets)
313 {
314  std::shared_ptr<ElementDataCache<double>> serial_nodes_cache;
315 
316  // collects nodes_ data (coordinates)
317  auto serial_nodes = nodes_->gather(node_ds_, node_4_loc_);
318 
319  if (el_ds_->myp()==0) serial_nodes_cache = std::dynamic_pointer_cast< ElementDataCache<double> >(serial_nodes);
320  return serial_nodes_cache;
321 }
322 
323 
324 std::shared_ptr<ElementDataCache<unsigned int>> OutputMesh::make_serial_connectivity_cache(std::shared_ptr<ElementDataCache<unsigned int>> global_offsets)
325 {
326  std::shared_ptr<ElementDataCache<unsigned int>> serial_connectivity_cache;
327 
328  // re-number connectivity indices from local to global
329  auto &conn_vec = *( connectivity_->get_data().get() );
330  ElementDataCache<unsigned int> global_conn("connectivity", (unsigned int)1, conn_vec.size()); // holds global indices of nodes
331  auto &global_conn_vec = *( global_conn.get_data().get() );
332  for(unsigned int i=0; i<conn_vec.size(); i++) {
333  global_conn_vec[i] = node_4_loc_[ conn_vec[i] ];
334  }
335 
336  // collects global connectivities
337  auto &local_offset_vec = *( offsets_->get_data().get() );
338  auto global_fix_size_conn = global_conn.element_node_cache_fixed_size(local_offset_vec);
339  auto collective_conn = global_fix_size_conn->gather(el_ds_, el_4_loc_);
340 
341  if (el_ds_->myp()==0) {
342  auto &offset_vec = *( global_offsets->get_data().get() );
343  serial_connectivity_cache = std::dynamic_pointer_cast< ElementDataCache<unsigned int> >( collective_conn->element_node_cache_optimize_size(offset_vec) );
344  }
345  return serial_connectivity_cache;
346 }
347 
348 
349 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
350 
352 : OutputMeshBase(mesh)
353 {
354  this->mesh_type_ = MeshType::discont;
355 }
356 
358 : OutputMeshBase(mesh, in_rec)
359 {
360  this->mesh_type_ = MeshType::discont;
361 }
362 
363 
365 {
366 }
367 
368 
369 template<int dim>
372  const ElementAccessor<spacedim> &ele_acc)
373 {
374  static const unsigned int n_subelements = 1 << dim; //2^dim
375 
376 // The refinement of elements for the output mesh is done using edge splitting
377 // technique (so called red refinement). Since we use this only for better output
378 // visualization of non-polynomial solutions, we do not care for existence of hanging
379 // nodes.
380 // In 2D case, it is straightforward process: find the midpoints of all sides,
381 // connect them and generate 4 triangles. These triangles are congruent and have
382 // equal surface areas.
383 // On the other hand, the 3D case is more complicated. After splitting the
384 // edges, we obtain 4 tetrahedra at the vertices of the original one. The octahedron
385 // that remains in the middle can be subdivided according to one of its three
386 // diagonals. Only the choice of the shortest octahedron diagonal leads to a regular
387 // tetrahedra decomposition. This algorithm originally comes from Bey.
388 // Bey's algorithm (red refinement of tetrahedron):
389 // p.29 https://www5.in.tum.de/pub/Joshi2016_Thesis.pdf
390 // p.108 http://www.bcamath.org/documentos_public/archivos/publicaciones/sergey_book.pdf
391 // https://www.math.uci.edu/~chenlong/iFEM/doc/html/uniformrefine3doc.html#1
392 // J. Bey. Simplicial grid refinement: on Freudenthal's algorithm and the optimal number of congruence classes.
393 // Numer. Math. 85(1):1--29, 2000. p11 Algorithm: RedRefinement3D.
394 // p.4 http://www.vis.uni-stuttgart.de/uploads/tx_vispublications/vis97-grosso.pdf
395 
396  // connectivity of refined element
397  // these arrays are hardwired to the current reference element
398  static const std::vector<unsigned int> conn[] = {
399  {}, //0D
400 
401  //1D:
402  // 0,1 original nodes, 2 is in the middle
403  // get 2 elements
404  {0, 2,
405  2, 1},
406 
407  //2D:
408  // 0,1,2 original nodes
409  // 3,4,5 nodes are in the middle of sides 0,1,2 in the respective order
410  // get 4 elements
411  {0, 3, 4,
412  3, 1, 5,
413  4, 5, 2,
414  3, 5, 4},
415 
416  //3D:
417  // 0,1,2,3 original nodes
418  // 4,5,6,7,8,9 are nodes in the middle of edges 0,1,2,3,4,5 in the respective order
419  // first 4 tetrahedra are at the original nodes
420  // next 4 tetrahedra are from the inner octahedron - 4 choices according to the diagonal
421  {1, 7, 4, 8,
422  7, 2, 5, 9,
423  4, 5, 0, 6,
424  8, 9, 6, 3,
425  7, 4, 8, 9, // 4-9 octahedron diagonal
426  7, 4, 5, 9,
427  4, 8, 9, 6,
428  4, 5, 9, 6},
429 
430  {1, 7, 4, 8,
431  7, 2, 5, 9,
432  4, 5, 0, 6,
433  8, 9, 6, 3,
434  8, 5, 4, 6, // 5-8 octahedron diagonal
435  5, 8, 9, 6,
436  5, 8, 4, 7,
437  8, 5, 9, 7},
438 
439  {1, 7, 4, 8,
440  7, 2, 5, 9,
441  4, 5, 0, 6,
442  8, 9, 6, 3,
443  6, 8, 7, 9, // 6-7 octahedron diagonal
444  8, 6, 7, 4,
445  6, 5, 7, 4,
446  5, 6, 7, 9}
447  };
448 // DBGMSG("level = %d, %d\n", aux_element.level, max_refinement_level_);
449 
450  ASSERT_EQ(dim, aux_element.nodes.size()-1);
451 
452  // if not refining any further, push into final vector
453  if( ! refinement_criterion(aux_element, ele_acc) ) {
454  refinement.push_back(aux_element);
455  return;
456  }
457 
458  std::vector<AuxElement> subelements(n_subelements);
459 
460  const unsigned int n_old_nodes = RefElement<dim>::n_nodes,
461  n_new_nodes = RefElement<dim>::n_lines; // new points are in the center of lines
462 
463  // auxiliary vectors
464  std::vector<Space<spacedim>::Point> nodes = aux_element.nodes;
465 
466  nodes.reserve(n_old_nodes+n_new_nodes);
467 
468  // create new points in the element
469  for(unsigned int e=0; e < n_new_nodes; e++)
470  {
473  nodes.push_back( p / 2.0);
474  //nodes.back().print();
475  }
476 
477  unsigned int diagonal = 0;
478  // find shortest diagonal: [0]:4-9, [1]:5-8 or [2]:6-7
479  if(dim == 3){
480  double min_diagonal = arma::norm(nodes[4]-nodes[9],2);
481  double d = arma::norm(nodes[5]-nodes[8],2);
482  if(d < min_diagonal){
483  min_diagonal = d;
484  diagonal = 1;
485  }
486  d = arma::norm(nodes[6]-nodes[7],2);
487  if(d < min_diagonal){
488  min_diagonal = d;
489  diagonal = 2;
490  }
491  }
492 
493  for(unsigned int i=0; i < n_subelements; i++)
494  {
495  AuxElement& sub_ele = subelements[i];
496  sub_ele.nodes.resize(n_old_nodes);
497  sub_ele.level = aux_element.level+1;
498 
499  // over nodes
500  for(unsigned int j=0; j < n_old_nodes; j++)
501  {
502  unsigned int conn_id = (n_old_nodes)*i + j;
503  sub_ele.nodes[j] = nodes[conn[dim+diagonal][conn_id]];
504  }
505  refine_aux_element<dim>(sub_ele, refinement, ele_acc);
506  }
507 }
508 
509 
510 
511 template void OutputMeshDiscontinuous::refine_aux_element<1>(const OutputMeshDiscontinuous::AuxElement&,std::vector< OutputMeshDiscontinuous::AuxElement >&, const ElementAccessor<spacedim> &);
512 template void OutputMeshDiscontinuous::refine_aux_element<2>(const OutputMeshDiscontinuous::AuxElement&,std::vector< OutputMeshDiscontinuous::AuxElement >&, const ElementAccessor<spacedim> &);
513 template void OutputMeshDiscontinuous::refine_aux_element<3>(const OutputMeshDiscontinuous::AuxElement&,std::vector< OutputMeshDiscontinuous::AuxElement >&, const ElementAccessor<spacedim> &);
514 
515 
517  const ElementAccessor<spacedim> &ele_acc)
518 {
519  // check refinement criteria:
520 
521  //first check max. level
522  bool refine = refinement_criterion_uniform(aux_ele);
523 
524  //if max. level not reached and refinement by error is set
525  if(refine && refine_by_error_)
526  {
527  // compute centre of aux element
528  Space<spacedim>::Point centre({0,0,0});
529  for(auto& v : aux_ele.nodes ) centre += v;
530  centre = centre/aux_ele.nodes.size();
531  return refinement_criterion_error(aux_ele, centre, ele_acc);
532  }
533 
534  return refine;
535 }
536 
538 {
539  return (ele.level < max_level_);
540 }
541 
543  const Space<spacedim>::Point &centre,
544  const ElementAccessor<spacedim> &ele_acc
545  )
546 {
547  ASSERT(error_control_field_func_).error("Error control field not set!");
548 
549  // evaluate at nodes and center in a single call
550  std::vector<double> val_list(ele.nodes.size()+1);
551  Armor::array point_list(spacedim,1,1+ele.nodes.size());
552  point_list.set(0) = centre;
553  unsigned int i=0;
554  for (auto node : ele.nodes) point_list.set(++i) = node;
555  error_control_field_func_(point_list, ele_acc, val_list);
556 
557  //TODO: compute L1 or L2 error using standard quadrature
558 
559  //compare average value at nodes with value at center
560 
561  double average_val = 0.0;
562  for(unsigned int i=1; i<ele.nodes.size()+1; ++i)//(double& v: nodes_val)
563  average_val += val_list[i];
564  average_val = average_val / ele.nodes.size();
565 
566  double diff = std::abs((average_val - val_list[0])/val_list[0]);
567 // DebugOut().fmt("diff: {} {} {}\n", diff, average_val, val_list[0]);
568  return ( diff > refinement_error_tolerance_);
569 
570 }
571 
572 
573 std::shared_ptr<OutputMeshBase> OutputMeshDiscontinuous::construct_mesh()
574 {
575  return std::make_shared<OutputMeshDiscontinuous>(*orig_mesh_);
576 }
577 
578 
579 std::shared_ptr<ElementDataCache<double>> OutputMeshDiscontinuous::make_serial_nodes_cache(std::shared_ptr<ElementDataCache<unsigned int>> global_offsets)
580 {
581  std::shared_ptr<ElementDataCache<double>> serial_nodes_cache;
582 
583  // Create helper cache of discontinuous node data ordering by elements
584  std::shared_ptr< ElementDataCache<double> > discont_node_cache = std::make_shared<ElementDataCache<double>>("",
585  (unsigned int) ElementDataCacheBase::N_VECTOR, this->connectivity_->n_values());
586  auto &discont_node_vec = *( discont_node_cache->get_data().get() );
587  auto &local_nodes_vec = *( this->nodes_->get_data().get() );
588  auto &local_conn_vec = *( this->connectivity_->get_data().get() );
589  auto &local_offset_vec = *( this->offsets_->get_data().get() );
590  unsigned int i_old, i_new;
591  for (unsigned int i_conn=0; i_conn<this->connectivity_->n_values(); ++i_conn) {
592  i_old = local_conn_vec[i_conn] * ElementDataCacheBase::N_VECTOR;
593  i_new = i_conn * ElementDataCacheBase::N_VECTOR;
594  for(unsigned int i = 0; i < ElementDataCacheBase::N_VECTOR; i++) {
595  discont_node_vec[i_new+i] = local_nodes_vec[i_old+i];
596  }
597  }
598  // Collects node data
599  auto fix_size_node_cache = discont_node_cache->element_node_cache_fixed_size(local_offset_vec);
600  auto collect_fix_size_node_cache = fix_size_node_cache->gather(el_ds_, el_4_loc_);
601 
602  if (el_ds_->myp()==0) {
603  auto &offset_vec = *( global_offsets->get_data().get() );
604  serial_nodes_cache = std::dynamic_pointer_cast< ElementDataCache<double> >(collect_fix_size_node_cache->element_node_cache_optimize_size(offset_vec));
605  }
606  return serial_nodes_cache;
607 }
608 
609 
610 std::shared_ptr<ElementDataCache<unsigned int>> OutputMeshDiscontinuous::make_serial_connectivity_cache(std::shared_ptr<ElementDataCache<unsigned int>> global_offsets)
611 {
612  std::shared_ptr<ElementDataCache<unsigned int>> serial_connectivity_cache;
613 
614  if (el_ds_->myp()==0) {
615  auto &offset_vec = *( global_offsets->get_data().get() );
616  serial_connectivity_cache = std::make_shared<ElementDataCache<unsigned int>>(
617  "connectivity", (unsigned int) ElementDataCacheBase::N_SCALAR, offset_vec[offset_vec.size()-1]);
618  auto &conn_vec = *( serial_connectivity_cache->get_data().get() );
619  for (unsigned int i=0; i<conn_vec.size(); ++i) conn_vec[i] = i;
620  }
621  return serial_connectivity_cache;
622 }
623 
624 
626 {
627  ASSERT( !is_created() ).error("Multiple initialization of OutputMesh!\n");
628 
629  DebugOut() << "Create refined discontinuous submesh containing only local elements.";
630  // initial guess of size: n_elements
631  nodes_ = std::make_shared<ElementDataCache<double>>("",(unsigned int)ElementDataCacheBase::N_VECTOR,0);
632  connectivity_ = std::make_shared<ElementDataCache<unsigned int>>("connectivity",(unsigned int) ElementDataCacheBase::N_SCALAR,0);
633  offsets_ = std::make_shared<ElementDataCache<unsigned int>>("offsets",(unsigned int) ElementDataCacheBase::N_SCALAR,0);
634  orig_element_indices_ = std::make_shared<std::vector<unsigned int>>();
635 
636  // index of last node added; set at the end of original ones
637  unsigned int last_offset = 0;
638 
639  auto &node_vec = *( nodes_->get_data().get() );
640  auto &conn_vec = *( connectivity_->get_data().get() );
641  auto &offset_vec = *( offsets_->get_data().get() );
642 
643  node_vec.reserve(4*orig_mesh_->n_nodes());
644  conn_vec.reserve(4*4*orig_mesh_->n_elements());
645  offset_vec.reserve(4*orig_mesh_->n_elements()+4);
646  offset_vec.push_back(0);
647 
648  LongIdx *el_4_loc = orig_mesh_->get_el_4_loc();
649  const unsigned int n_local_elements = orig_mesh_->get_el_ds()->lsize();
650 
651  for (unsigned int loc_el = 0; loc_el < n_local_elements; loc_el++) {
652  auto ele = orig_mesh_->element_accessor( el_4_loc[loc_el] );
653  const unsigned int
654  dim = ele->dim(),
655  ele_idx = ele.idx();
656 
657  AuxElement aux_ele;
658  aux_ele.nodes.resize(ele->n_nodes());
659  aux_ele.level = 0;
660 
661  unsigned int li;
662  for (li=0; li<ele->n_nodes(); li++) {
663  aux_ele.nodes[li] = *ele.node(li);
664  }
665 
666  std::vector<AuxElement> refinement;
667 
668  switch(dim){
669  case 1: this->refine_aux_element<1>(aux_ele, refinement, ele); break;
670  case 2: this->refine_aux_element<2>(aux_ele, refinement, ele); break;
671  case 3: this->refine_aux_element<3>(aux_ele, refinement, ele); break;
672  default: ASSERT_PERMANENT(0).error("Should not happen.\n");
673  }
674 
675  //skip unrefined element
676 // if(refinement.size() < 2) continue;
677  unsigned int node_offset = node_vec.size(),
678  con_offset = conn_vec.size();
679  node_vec.resize(node_vec.size() + (refinement.size() * (dim+1))*spacedim);
680  conn_vec.resize(conn_vec.size() + refinement.size()*(dim+1));
681 // orig_element_indices_->resize(orig_element_indices_->size() + refinement.size()*(dim+1));
682 
683 // DebugOut() << "ref size = " << refinement.size() << "\n";
684  //gather coords and connectivity (in a continous way inside element)
685  for(unsigned int i=0; i < refinement.size(); i++)
686  {
687  last_offset += dim+1;
688  offset_vec.push_back(last_offset);
689  (*orig_element_indices_).push_back(ele_idx);
690  for(unsigned int j=0; j < dim+1; j++)
691  {
692  unsigned int con = i*(dim+1) + j;
693  conn_vec[con_offset + con] = con_offset + con;
694 
695  for(unsigned int k=0; k < spacedim; k++) {
696  node_vec[node_offset + con*spacedim + k] = refinement[i].nodes[j][k];
697  }
698  }
699  }
700  }
701 
702  conn_vec.shrink_to_fit();
703  node_vec.shrink_to_fit();
704  offset_vec.shrink_to_fit();
705 
706  connectivity_->set_n_values(conn_vec.size());
707  nodes_->set_n_values(node_vec.size() / spacedim);
708  offsets_->set_n_values(offset_vec.size());
709 
710  // Create special distributions and arrays of local to global indexes of refined mesh
711  el_ds_ = new Distribution(offset_vec.size()-1, PETSC_COMM_WORLD);
712  node_ds_ = new Distribution(offset_vec[offset_vec.size()-1], PETSC_COMM_WORLD);
714  el_4_loc_ = new LongIdx [ el_ds_->lsize() ];
715  LongIdx global_el_idx = el_ds_->begin();
716  for (unsigned int i=0; i<el_ds_->lsize(); ++i, ++global_el_idx) {
717  el_4_loc_[i] = global_el_idx;
718  }
719  node_4_loc_ = new LongIdx [ node_ds_->lsize() ];
720  LongIdx global_node_idx = node_ds_->begin();
721  for (unsigned int i=0; i<node_ds_->lsize(); ++i, ++global_node_idx) {
722  node_4_loc_[i] = global_node_idx;
723  }
724 
725  mesh_type_ = MeshType::refined;
726 }
727 
728 
730 {
731  master_mesh_ = this->construct_mesh();
732  master_mesh_->offsets_ = this->offsets_;
733  master_mesh_->orig_element_indices_ = this->orig_element_indices_;
734 
735  auto &conn_vec = *( this->connectivity_->get_data().get() );
736  master_mesh_->connectivity_ = std::make_shared<ElementDataCache<unsigned int>>(
737  "connectivity", (unsigned int) ElementDataCacheBase::N_SCALAR, conn_vec.size());
738  auto &master_conn_vec = *( master_mesh_->connectivity_->get_data().get() );
739  for (unsigned int i=0; i<master_conn_vec.size(); ++i) master_conn_vec[i] = i;
740 
741  master_mesh_->nodes_ = std::make_shared<ElementDataCache<double>>(
742  "", (unsigned int) ElementDataCacheBase::N_VECTOR, conn_vec.size());
743  auto &node_vec = *( this->nodes_->get_data().get() );
744  auto &master_node_vec = *( master_mesh_->nodes_->get_data().get() );
745  unsigned int i_own, i_master, j;
746  for (unsigned int i=0; i<conn_vec.size(); ++i) {
747  i_own = conn_vec[i]*ElementDataCacheBase::N_VECTOR;
748  i_master = i*ElementDataCacheBase::N_VECTOR;
749  for (j=0; j<ElementDataCacheBase::N_VECTOR; ++j) {
750  master_node_vec[i_master+j] = node_vec[i_own+j];
751  }
752  }
753 
754  master_mesh_->mesh_type_ = this->mesh_type_;
755 }
#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
ArrayMatSet set(uint index)
Definition: armor.hh:838
unsigned int myp() const
get my processor
unsigned int begin(int proc) const
get starting local index
unsigned int end(int proc) const
get last local index +1
unsigned int lsize(int proc) const
get local size
unsigned int size() const
get global size
unsigned int np() const
get num of processors
NodeAccessor< 3 > node(unsigned int ni) const
Definition: accessors.hh:230
unsigned int n_values() const
std::shared_ptr< ElementDataCacheBase > gather(Distribution *distr, LongIdx *local_to_global) override
Implements ElementDataCacheBase::gather.
std::shared_ptr< ElementDataCacheBase > element_node_cache_fixed_size(std::vector< unsigned int > &offset_vec) override
Implements ElementDataCacheBase::element_node_cache_fixed_size.
CacheData get_data()
Return underlying vector of element data.
unsigned int dim() const
Definition: elements.h:120
unsigned int n_nodes() const
Definition: elements.h:125
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
Class for declaration of the input of type Bool.
Definition: type_base.hh:452
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
static Default optional()
The factory function to make an empty default value which is optional.
Definition: type_record.hh:124
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:534
Class for declaration of the integral input data.
Definition: type_base.hh:483
Record type proxy class.
Definition: type_record.hh:182
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:503
Class for declaration of the input data that are in string format.
Definition: type_base.hh:582
General iterator template. Provides iterator over objects of type Object in some container.
LongIdx * get_el_4_loc() const
Definition: mesh.h:122
unsigned int n_nodes() const
Definition: mesh.h:171
int find_node_id(unsigned int pos) const
Return node id (in GMSH file) of node of given position in node vector.
Definition: mesh.h:156
NodeAccessor< 3 > node(unsigned int idx) const
Create and return NodeAccessor to node of given idx.
Definition: mesh.cc:872
int find_elem_id(unsigned int pos) const
Return element id (in GMSH file) of element of given position in element vector.
Definition: mesh.h:138
Distribution * get_el_ds() const
Definition: mesh.h:119
unsigned int n_elements() const
Definition: mesh.h:111
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
LongIdx * get_node_4_loc() const
Definition: mesh.h:444
unsigned int n_local_nodes() const
Definition: mesh.h:447
Distribution * get_node_ds() const
Definition: mesh.h:441
Base class for Output mesh.
Definition: output_mesh.hh:71
double refinement_error_tolerance_
Tolerance for error criterion refinement.
Definition: output_mesh.hh:200
virtual std::shared_ptr< ElementDataCache< double > > make_serial_nodes_cache(std::shared_ptr< ElementDataCache< unsigned int >> global_offsets)=0
friend class OutputElement
Friend provides access to vectors for element accessor class.
Definition: output_mesh.hh:242
OutputMeshBase(Mesh &mesh)
Constructor. Takes computational mesh as a parameter.
Definition: output_mesh.cc:46
LongIdx * node_4_loc_
Index set assigning to local node index its global index.
Definition: output_mesh.hh:237
std::shared_ptr< ElementDataCache< unsigned int > > connectivity_
Vector maps the nodes to their coordinates in vector nodes_.
Definition: output_mesh.hh:208
void make_serial_master_mesh()
Synchronize parallel data and create serial COLECTIVE output mesh on zero process.
Definition: output_mesh.cc:217
std::function< void(const Armor::array &, const ElementAccessor< spacedim > &, std::vector< double > &)> ErrorControlFieldFunc
Definition: output_mesh.hh:77
virtual std::shared_ptr< OutputMeshBase > construct_mesh()=0
static const Input::Type::Record & get_input_type()
The specification of output mesh.
Definition: output_mesh.cc:31
std::shared_ptr< ElementDataCache< unsigned int > > get_elems_n_nodes()
Compute and return number of nodes for each elements (compute from offsets)
Definition: output_mesh.cc:256
unsigned int n_elements()
Returns number of element.
Definition: output_mesh.cc:101
friend class OutputMeshDiscontinuous
Definition: output_mesh.hh:247
std::shared_ptr< OutputMeshBase > master_mesh_
Definition: output_mesh.hh:227
bool is_created()
Check if nodes_, connectivity_ and offsets_ data caches are created.
Definition: output_mesh.cc:145
LongIdx * el_4_loc_
Index set assigning to local element index its global index.
Definition: output_mesh.hh:234
virtual std::shared_ptr< ElementDataCache< unsigned int > > make_serial_connectivity_cache(std::shared_ptr< ElementDataCache< unsigned int >> global_offsets)=0
Distribution * el_ds_
Parallel distribution of elements.
Definition: output_mesh.hh:236
std::shared_ptr< ElementDataCache< unsigned int > > node_ids_
Vector gets ids of nodes. Data is used in GMSH output.
Definition: output_mesh.hh:213
void create_sub_mesh()
Definition: output_mesh.cc:151
bool refine_by_error_
True, if output mesh is to be refined by error criterion.
Definition: output_mesh.hh:199
friend class OutputMesh
Definition: output_mesh.hh:246
ErrorControlFieldFunc error_control_field_func_
Refinement error control field function (hold value_list function of field).
Definition: output_mesh.hh:196
OutputElementIterator begin()
Gives iterator to the FIRST element of the output mesh.
Definition: output_mesh.cc:81
std::shared_ptr< std::vector< unsigned int > > orig_element_indices_
Vector of element indices in the computational mesh. (Important when refining.)
Definition: output_mesh.hh:203
const unsigned int max_level_
Maximal level of refinement.
Definition: output_mesh.hh:193
OutputElementIterator end()
Gives iterator to the LAST element of the output mesh.
Definition: output_mesh.cc:88
std::shared_ptr< ElementDataCache< unsigned int > > elem_ids_
Vector gets ids of elements. Data is used in GMSH output.
Definition: output_mesh.hh:215
std::shared_ptr< ElementDataCache< int > > partitions_
Vector gets partitions of elements. Data is used in GMSH output.
Definition: output_mesh.hh:219
std::shared_ptr< ElementDataCache< double > > nodes_
Vector of node coordinates. [spacedim x n_nodes].
Definition: output_mesh.hh:206
static const unsigned int spacedim
Shortcut instead of spacedim template. We suppose only spacedim=3 at the moment.
Definition: output_mesh.hh:74
MeshType mesh_type_
Type of OutputMesh.
Definition: output_mesh.hh:198
std::shared_ptr< ElementDataCache< unsigned int > > offsets_
Vector of offsets of node indices of elements. Maps elements to their nodes in connectivity_.
Definition: output_mesh.hh:210
unsigned int n_nodes()
Returns number of nodes.
Definition: output_mesh.cc:107
unsigned int n_local_nodes_
Hold number of local nodes (own + ghost), value is equal with size of node_4_loc array.
Definition: output_mesh.hh:239
Mesh * orig_mesh_
Pointer to the computational mesh.
Definition: output_mesh.hh:190
void set_error_control_field(ErrorControlFieldFunc error_control_field_func)
Selects the error control field computing function of output field set according to input record.
Definition: output_mesh.cc:96
LongIdx * loc_4_el_
Index set assigning to global index its local element index.
Definition: output_mesh.hh:235
virtual ~OutputMeshBase()
Definition: output_mesh.cc:70
Distribution * node_ds_
Parallel distribution of nodes. Depends on elements distribution.
Definition: output_mesh.hh:238
std::shared_ptr< ElementDataCache< unsigned int > > region_ids_
Vector gets ids of regions. Data is used in GMSH output.
Definition: output_mesh.hh:217
void create_id_caches()
Create nodes and elements data caches.
Definition: output_mesh.cc:113
bool refinement_criterion_error(const AuxElement &ele, const Space< spacedim >::Point &centre, const ElementAccessor< spacedim > &ele_acc)
Refinement flag - measures discretisation error according to error control field.
Definition: output_mesh.cc:542
void make_parallel_master_mesh() override
Overrides OutputMeshBase::make_parallel_master_mesh.
Definition: output_mesh.cc:729
void create_refined_sub_mesh() override
Implements OutputMeshBase::create_refined_sub_mesh.
Definition: output_mesh.cc:625
bool refinement_criterion_uniform(const AuxElement &ele)
Refinement flag - checks only maximal level of refinement.
Definition: output_mesh.cc:537
bool refinement_criterion(const AuxElement &ele, const ElementAccessor< spacedim > &ele_acc)
Collects different refinement criteria results.
Definition: output_mesh.cc:516
void refine_aux_element(const AuxElement &aux_element, std::vector< AuxElement > &refinement, const ElementAccessor< spacedim > &ele_acc)
Performs the actual refinement of AuxElement. Recurrent.
Definition: output_mesh.cc:370
std::shared_ptr< ElementDataCache< unsigned int > > make_serial_connectivity_cache(std::shared_ptr< ElementDataCache< unsigned int >> global_offsets) override
Implements OutputMeshBase::make_serial_connectivity_cache.
Definition: output_mesh.cc:610
std::shared_ptr< ElementDataCache< double > > make_serial_nodes_cache(std::shared_ptr< ElementDataCache< unsigned int >> global_offsets) override
Implements OutputMeshBase::make_serial_nodes_cache.
Definition: output_mesh.cc:579
std::shared_ptr< OutputMeshBase > construct_mesh() override
Implements OutputMeshBase::construct_mesh.
Definition: output_mesh.cc:573
std::shared_ptr< ElementDataCache< unsigned int > > make_serial_connectivity_cache(std::shared_ptr< ElementDataCache< unsigned int >> global_offsets) override
Implements OutputMeshBase::make_serial_connectivity_cache.
Definition: output_mesh.cc:324
std::shared_ptr< ElementDataCache< double > > make_serial_nodes_cache(std::shared_ptr< ElementDataCache< unsigned int >> global_offsets) override
Implements OutputMeshBase::make_serial_nodes_cache.
Definition: output_mesh.cc:312
std::shared_ptr< OutputMeshBase > construct_mesh() override
Implements OutputMeshBase::construct_mesh.
Definition: output_mesh.cc:306
bool refinement_criterion()
Definition: output_mesh.cc:299
void create_refined_sub_mesh() override
Implements OutputMeshBase::create_refined_sub_mesh.
Definition: output_mesh.cc:294
static const IdxVector<(InDim >OutDim ? InDim+1 :dim-InDim) > interact(TInteraction< OutDim, InDim > interaction, bool inv=false)
Armor::ArmaVec< double, spacedim > Point
Definition: point.hh:42
Support classes for parallel programing.
const unsigned int undef_idx
Definition: index_types.hh:32
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
Definition: index_types.hh:24
#define DebugOut()
Macro defining 'debug' record of log.
Definition: logger.hh:284
Class OutputElement and its iterator OutputElementIterator on the output mesh.
Classes for auxiliary output mesh.
Iter< OutputElement > OutputElementIterator
Definition: output_mesh.hh:37
#define FMT_UNUSED
Definition: posix.h:75
Implementation of range helper class.
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.
Auxiliary structure defining element of refined output mesh.
Definition: output_mesh.hh:295
std::vector< Space< spacedim >::Point > nodes
Definition: output_mesh.hh:296