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