Hermes2D  3.0
optimum_selector.cpp
1 #include <iostream>
2 #include "element_to_refine.h"
3 #include "optimum_selector.h"
4 #include "order_permutator.h"
5 
6 #define H2DST_ANY H2DST_VERTEX | H2DST_HORIZ_EDGE | H2DST_VERT_EDGE | H2DST_TRI_EDGE | H2DST_BUBBLE
7 
8 namespace Hermes
9 {
10  namespace Hermes2D
11  {
12  namespace RefinementSelectors
13  {
14  template<typename Scalar>
16  max_order, Shapeset* shapeset, const Range& vertex_order, const
17  Range& edge_bubble_order) : Selector<Scalar>(shapeset->get_min_order(), max_order), cand_list(cand_list), shapeset(shapeset), dof_score_exponent(1.0)
18  {
19  if (shapeset == nullptr)
20  throw Exceptions::NullException(3);
21 
22  num_shapes = new int***[2];
23  for (int i = 0; i < 2; i++)
24  {
28  num_shapes[i] = new int**[H2D_NUM_SHAPES_SIZE];
29  for (int j = 0; j < H2D_NUM_SHAPES_SIZE; j++)
30  {
31  num_shapes[i][j] = new int*[H2D_NUM_SHAPES_SIZE];
32  for (int k = 0; k < H2D_NUM_SHAPES_SIZE; k++)
33  {
34  num_shapes[i][j][k] = new int[6];
35  memset(num_shapes[i][j][k], 0, 6 * sizeof(int));
36  }
37  }
38  }
39 
40  //build shape indices
41  build_shape_indices(HERMES_MODE_TRIANGLE, vertex_order, edge_bubble_order);
42  build_shape_indices(HERMES_MODE_QUAD, vertex_order, edge_bubble_order);
43  }
44 
45  template<typename Scalar>
47  {
48  for (int i = 0; i < 2; i++)
49  {
50  for (int j = 0; j < H2D_NUM_SHAPES_SIZE; j++)
51  {
52  for (int k = 0; k < H2D_NUM_SHAPES_SIZE; k++)
53  {
54  delete[] num_shapes[i][j][k];
55  }
56  delete[] num_shapes[i][j];
57  }
58  delete[] num_shapes[i];
59  }
60  delete[] num_shapes;
61  }
62 
63  template<typename Scalar>
64  void OptimumSelector<Scalar>::add_bubble_shape_index(int order_h, int order_v, std::map<int, bool>& used_shape_index, std::vector<ShapeInx>& indices, ElementMode2D mode)
65  {
66  int quad_order = H2D_MAKE_QUAD_ORDER(order_h, order_v);
67  const int num_bubbles = shapeset->get_num_bubbles(quad_order, mode);
68  short* bubble_inxs = shapeset->get_bubble_indices(quad_order, mode);
69  for (int j = 0; j < num_bubbles; j++)
70  {
71  int inx_bubble = bubble_inxs[j];
72  if (used_shape_index.find(inx_bubble) == used_shape_index.end())
73  {
74  used_shape_index[inx_bubble] = true;
75  indices.push_back(ShapeInx(order_h, order_v, inx_bubble, H2DST_BUBBLE));
76  for (int order_h_i = order_h + 1; order_h_i < H2D_NUM_SHAPES_SIZE; order_h_i++)
77  for (int order_v_i = order_v + 1; order_v_i < H2D_NUM_SHAPES_SIZE; order_v_i++)
78  {
79  num_shapes[mode][order_h_i][order_v_i][H2DSI_BUBBLE]++;
80  num_shapes[mode][order_h_i][order_v_i][H2DSI_ANY]++;
81  }
82  num_shapes[mode][0][0][H2DSI_BUBBLE]++;
83  num_shapes[mode][0][0][H2DSI_ANY]++;
84  }
85  }
86  }
87 
88  template<typename Scalar>
89  void OptimumSelector<Scalar>::build_shape_indices(const ElementMode2D mode, const Range& vertex_order, const Range& edge_bubble_order)
90  {
91  std::vector<ShapeInx> &indices = shape_indices[mode];
92  int* next_order = this->next_order_shape[mode];
93  int& max_shape_inx = this->max_shape_inx[mode];
94  int num_edges = (mode == HERMES_MODE_QUAD) ? 4 : 3;
95  bool &has_vertex = has_vertex_shape[mode];
96  bool &has_edge = has_edge_shape[mode];
97  bool &has_bubble = has_bubble_shape[mode];
98 
99  //cleanup
100  indices.clear();
101  indices.reserve((H2DRS_MAX_ORDER + 2) * (H2DRS_MAX_ORDER + 2));
102  has_vertex = has_edge = has_bubble = false;
103 
104  //get total range of orders
105  Range order_range = Range::make_envelope(vertex_order, edge_bubble_order);
106 
107  //prepare array of used shape indices
108  std::map<int, bool> used_shape_index;
109 
110  //for all orders
111  max_shape_inx = 0;
112  int examined_shape = 0;
113  for (unsigned char i = order_range.lower(); i <= 10; i++)
114  {
115  if (shapeset->get_space_type() != HERMES_L2_SPACE && shapeset->get_space_type() != HERMES_L2_MARKERWISE_CONST_SPACE)
116  {
117  //vertex functions
118  if (vertex_order.is_in_closed(i))
119  {
120  for (int j = 0; j < num_edges; j++)
121  {
122  int inx = shapeset->get_vertex_index(j, mode);
123  if (inx >= 0)
124  {
125  used_shape_index[inx] = true;
126  indices.push_back(ShapeInx(1, 1, inx, H2DST_VERTEX));
127  if (mode == HERMES_MODE_QUAD)
128  {
129  for (int order_h_i = 2; order_h_i < H2D_NUM_SHAPES_SIZE; order_h_i++)
130  for (int order_v_i = 2; order_v_i < H2D_NUM_SHAPES_SIZE; order_v_i++)
131  {
132  num_shapes[mode][order_h_i][order_v_i][H2DSI_VERTEX]++;
133  num_shapes[mode][order_h_i][order_v_i][H2DSI_ANY]++;
134  }
135  num_shapes[mode][0][0][H2DSI_VERTEX]++;
136  num_shapes[mode][0][0][H2DSI_ANY]++;
137  }
138  else
139  {
140  for (int order_h_i = 2; order_h_i < H2D_NUM_SHAPES_SIZE; order_h_i++)
141  {
142  num_shapes[mode][order_h_i][0][H2DSI_VERTEX]++;
143  num_shapes[mode][order_h_i][0][H2DSI_ANY]++;
144  }
145  num_shapes[mode][0][0][H2DSI_VERTEX]++;
146  num_shapes[mode][0][0][H2DSI_ANY]++;
147  }
148  has_vertex = true;
149  }
150  }
151  }
152 
153  //edge functions
154  if (edge_bubble_order.is_in_closed(i))
155  {
156  //edge functions
157  if (mode == HERMES_MODE_QUAD)
158  {
159  for (int j = 0; j < num_edges; j++)
160  {
161  int inx = shapeset->get_edge_index(j, 0, i, HERMES_MODE_QUAD);
162  if (inx >= 0)
163  {
164  used_shape_index[inx] = true;
165  if ((j & 1) == 0)//horizontal edge
166  {
167  indices.push_back(ShapeInx(i, 0, inx, H2DST_HORIZ_EDGE));
168  for (int order_h_i = i + 1; order_h_i < H2D_NUM_SHAPES_SIZE; order_h_i++)
169  for (int order_v_i = 1; order_v_i < H2D_NUM_SHAPES_SIZE; order_v_i++)
170  {
171  num_shapes[mode][order_h_i][order_v_i][H2DST_HORIZ_EDGE]++;
172  num_shapes[mode][order_h_i][order_v_i][H2DSI_ANY]++;
173  }
174  num_shapes[mode][0][0][H2DST_HORIZ_EDGE]++;
175  num_shapes[mode][0][0][H2DSI_ANY]++;
176  }
177  else //vertical edge
178  {
179  indices.push_back(ShapeInx(0, i, inx, H2DST_VERT_EDGE));
180  for (int order_h_i = 1; order_h_i < H2D_NUM_SHAPES_SIZE; order_h_i++)
181  for (int order_v_i = i + 1; order_v_i < H2D_NUM_SHAPES_SIZE; order_v_i++)
182  {
183  num_shapes[mode][order_h_i][order_v_i][H2DSI_VERT_EDGE]++;
184  num_shapes[mode][order_h_i][order_v_i][H2DSI_ANY]++;
185  }
186  num_shapes[mode][0][0][H2DSI_VERT_EDGE]++;
187  num_shapes[mode][0][0][H2DSI_ANY]++;
188  }
189  has_edge = true;
190  }
191  }
192  }
193  else
194  {
195  for (int j = 0; j < num_edges; j++)
196  {
197  int inx = shapeset->get_edge_index(j, 0, i, HERMES_MODE_TRIANGLE);
198  if (inx >= 0)
199  {
200  used_shape_index[inx] = true;
201  indices.push_back(ShapeInx(i, i, inx, H2DST_TRI_EDGE));
202  for (int order_h_i = i + 1; order_h_i < H2D_NUM_SHAPES_SIZE; order_h_i++)
203  {
204  num_shapes[mode][order_h_i][0][H2DSI_TRI_EDGE]++;
205  num_shapes[mode][order_h_i][0][H2DSI_ANY]++;
206  }
207  num_shapes[mode][0][0][H2DSI_TRI_EDGE]++;
208  num_shapes[mode][0][0][H2DSI_ANY]++;
209  has_edge = true;
210  }
211  }
212  }
213  }
214 
215  //bubble functions
216  if (mode == HERMES_MODE_QUAD)
217  {
218  //NOTE: shapeset returns a set of all possible bubble functions and it is not possible to identify which is the smallest
219  // order of an element which contains this function, e.g., in a case of a Hcurl and an element of an order 1/1, it returns
220  // a list that contains a function of poly-order 2/0. Also, order of indices is not given.
221  int order = i;
222  unsigned num_indices_prev = indices.size();
223  for (int order_other = edge_bubble_order.lower(); order_other <= order; order_other++)
224  {
225  add_bubble_shape_index(order, order_other, used_shape_index, indices, mode);
226  add_bubble_shape_index(order_other, order, used_shape_index, indices, mode);
227  }
228 
229  //check if indices were added
230  if (num_indices_prev < indices.size())
231  has_bubble = true;
232  }
233  else { //triangles
234  int order = i;
235  unsigned short num_bubbles = shapeset->get_num_bubbles(order, mode);
236  short* bubble_inxs = shapeset->get_bubble_indices(order, mode);
237  for (unsigned short j = 0; j < num_bubbles; j++)
238  {
239  int inx_bubble = bubble_inxs[j];
240  if (used_shape_index.find(inx_bubble) == used_shape_index.end())
241  {
242  used_shape_index[inx_bubble] = true;
243  indices.push_back(ShapeInx(order, order, inx_bubble, H2DST_BUBBLE));
244  for (int order_h_i = order + 1; order_h_i < H2D_NUM_SHAPES_SIZE; order_h_i++)
245  {
246  num_shapes[mode][order_h_i][0][H2DSI_BUBBLE]++;
247  num_shapes[mode][order_h_i][0][H2DSI_ANY]++;
248  }
249  num_shapes[mode][0][0][H2DSI_BUBBLE]++;
250  num_shapes[mode][0][0][H2DSI_ANY]++;
251  has_bubble = true;
252  }
253  }
254  }
255 
256  //store index of the next order
257  next_order[i] = (int)indices.size();
258 
259  //update maximum
260  while (examined_shape < next_order[i])
261  {
262  max_shape_inx = std::max(max_shape_inx, indices[examined_shape].inx);
263  examined_shape++;
264  }
265  }
266  else
267  next_order[i] = (int)indices.size();
268  }
269  }
270 
271  template<typename Scalar>
272  int OptimumSelector<Scalar>::calc_num_shapes(int mode, int order_h, int order_v, int allowed_type_mask)
273  {
274  //test whether the evaluation is necessary
275  bool full_eval = false;
276  if ((allowed_type_mask & H2DST_VERTEX) != 0)
277  full_eval |= has_vertex_shape[mode];
278  if ((allowed_type_mask & (H2DST_HORIZ_EDGE | H2DST_VERT_EDGE | H2DST_TRI_EDGE)) != 0)
279  full_eval |= has_edge_shape[mode];
280  if ((allowed_type_mask & H2DST_BUBBLE) != 0)
281  full_eval |= has_bubble_shape[mode];
282 
283  //evaluate
284  if (full_eval)
285  {
286  std::vector<ShapeInx>& shapes = shape_indices[mode];
287  int num = 0;
288  typename std::vector<ShapeInx>::const_iterator shape = shapes.begin();
289  while (shape != shapes.end())
290  {
291  if (((int)shape->type & allowed_type_mask) != 0)
292  {
293  if ((order_h == H2DRS_ORDER_ANY || shape->order_h <= order_h) && (order_v == H2DRS_ORDER_ANY || shape->order_v <= order_v))
294  num++;
295  }
296  shape++;
297  }
298  return num;
299  }
300  else
301  return 0;
302  }
303 
304  template<typename Scalar>
305  void OptimumSelector<Scalar>::append_candidates_split(std::vector<Cand>& candidates, const int start_order, const int last_order, const RefinementType split, bool iso_p)
306  {
307  //check whether end orders are not lower than start orders
308  if (last_order < 0 || start_order < 0)
309  return;
310  if (H2D_GET_H_ORDER(start_order) > H2D_GET_H_ORDER(last_order) || H2D_GET_V_ORDER(start_order) > H2D_GET_V_ORDER(last_order))
311  return;
312 
313  //get number of sons
314  const int num_sons = get_refin_sons(split);
315 
316  //initialize orders
317  unsigned short quad_orders[H2D_MAX_ELEMENT_SONS];
319  for (int i = 0; i < num_sons; i++)
320  {
321  quad_orders[i] = start_order;
322  quad_perms[i] = OrderPermutator(start_order, last_order, iso_p, &quad_orders[i]);
323  }
324  for (int i = num_sons; i < H2D_MAX_ELEMENT_SONS; i++)
325  quad_orders[i] = 0;
326 
327  //generate permutations of orders
328  bool quit = false;
329  while (!quit)
330  {
331  do { //create permutation of son 0
332  candidates.push_back(Cand(split, quad_orders));
333  } while (quad_perms[0].next());
334 
335  //reset son 0
336  quad_perms[0].reset();
337 
338  //increment orders of other sons
339  int inx_son = 1;
340  while (inx_son < num_sons && !quad_perms[inx_son].next())
341  {
342  //reset order of the son
343  quad_perms[inx_son].reset();
344  inx_son++;
345  }
346  if (inx_son >= num_sons)
347  quit = true;
348  }
349  }
350 
351  template<typename Scalar>
352  std::vector<Cand> OptimumSelector<Scalar>::create_candidates(Element* e, int quad_order)
353  {
354  std::vector<Cand> candidates;
355 
356  // Get the current order range.
357  int current_min_order, current_max_order;
358  this->get_current_order_range(e, current_min_order, current_max_order);
359 
360  int order_h = H2D_GET_H_ORDER(quad_order), order_v = H2D_GET_V_ORDER(quad_order);
361 
362  if (current_max_order < std::max(order_h, order_v))
363  current_max_order = std::max(order_h, order_v);
364 
365  bool tri = e->is_triangle();
366 
367  //clear list of candidates
368  candidates.reserve(H2DRS_ASSUMED_MAX_CANDS);
369 
370  //generate all P-candidates (start from intention of generating all possible candidates
371  //and restrict it according to the given adapt-type)
372  bool iso_p = false;
373  int last_order_h = std::min(current_max_order, order_h + H2DRS_MAX_ORDER_INC), last_order_v = std::min(current_max_order, order_v + H2DRS_MAX_ORDER_INC);
374  int last_order = H2D_MAKE_QUAD_ORDER(last_order_h, last_order_v);
375  switch (cand_list)
376  {
377  case H2D_H_ISO:
378  case H2D_H_ANISO: last_order = quad_order; break; //no P-candidates except the original candidate
379  case H2D_P_ISO:
380  case H2D_HP_ISO:
381  case H2D_HP_ANISO_H: iso_p = true; break; //iso change of orders
382  }
383  append_candidates_split(candidates, quad_order, last_order, H2D_REFINEMENT_P, tri || iso_p);
384 
385  //generate all H-candidates
386  iso_p = false;
387  int start_order_h = std::max(current_min_order, order_h - 1), start_order_v = std::max(current_min_order, order_v - 1);
388  int start_order = H2D_MAKE_QUAD_ORDER(start_order_h, start_order_v);
389  last_order_h = std::min(current_max_order, order_h + H2DRS_MAX_ORDER_INC), last_order_v = std::min(current_max_order, order_v + H2DRS_MAX_ORDER_INC);
390  last_order = H2D_MAKE_QUAD_ORDER(last_order_h, last_order_v);
391  switch (cand_list)
392  {
393  case H2D_H_ISO:
394  case H2D_H_ANISO:
395  last_order = start_order = quad_order; break; //no only one candidate will be created
396  case H2D_P_ISO:
397  case H2D_P_ANISO: last_order = -1; break; //no H-candidate will be generated
398  case H2D_HP_ISO:
399  case H2D_HP_ANISO_H: iso_p = true; break; //iso change of orders
400  }
401  append_candidates_split(candidates, start_order, last_order, H2D_REFINEMENT_H, tri || iso_p);
402 
403  //generate all ANISO-candidates
405  if (!tri && e->iro_cache < 8 && (cand_list == H2D_H_ANISO || cand_list == H2D_HP_ANISO_H || cand_list == H2D_HP_ANISO))
406  {
407  iso_p = false;
408  int start_order_hz = H2D_MAKE_QUAD_ORDER(order_h, std::max(current_min_order, (order_v - 1)));
409  int last_order_hz = H2D_MAKE_QUAD_ORDER(std::min(current_max_order, order_h + H2DRS_MAX_ORDER_INC), std::min(order_v, H2D_GET_V_ORDER(start_order) + H2DRS_MAX_ORDER_INC));
410  int start_order_vt = H2D_MAKE_QUAD_ORDER(std::max(current_min_order, (order_h - 1)), order_v);
411  int last_order_vt = H2D_MAKE_QUAD_ORDER(std::min(order_h, H2D_GET_H_ORDER(start_order) + H2DRS_MAX_ORDER_INC), std::min(current_max_order, order_v + H2DRS_MAX_ORDER_INC));
412  switch (cand_list)
413  {
414  case H2D_H_ANISO:
415  last_order_hz = start_order_hz = quad_order;
416  last_order_vt = start_order_vt = quad_order;
417  //only one candidate will be created
418  break;
419  case H2D_HP_ANISO_H: iso_p = true; break; //iso change of orders
420  }
421  if (iso_p) { //make orders uniform: take mininmum order since nonuniformity is caused by different handling of orders along directions
422  int order = std::min(H2D_GET_H_ORDER(start_order_hz), H2D_GET_V_ORDER(start_order_hz));
423  start_order_hz = H2D_MAKE_QUAD_ORDER(order, order);
424  order = std::min(H2D_GET_H_ORDER(start_order_vt), H2D_GET_V_ORDER(start_order_vt));
425  start_order_vt = H2D_MAKE_QUAD_ORDER(order, order);
426 
427  order = std::min(H2D_GET_H_ORDER(last_order_hz), H2D_GET_V_ORDER(last_order_hz));
428  last_order_hz = H2D_MAKE_QUAD_ORDER(order, order);
429  order = std::min(H2D_GET_H_ORDER(last_order_vt), H2D_GET_V_ORDER(last_order_vt));
430  last_order_vt = H2D_MAKE_QUAD_ORDER(order, order);
431  }
432  append_candidates_split(candidates, start_order_hz, last_order_hz, H2D_REFINEMENT_H_ANISO_H, iso_p);
433  append_candidates_split(candidates, start_order_vt, last_order_vt, H2D_REFINEMENT_H_ANISO_V, iso_p);
434  }
435 
436  return candidates;
437  }
438 
439  template<typename Scalar>
440  void OptimumSelector<Scalar>::update_cands_info(std::vector<Cand>& candidates, CandsInfo& info_h, CandsInfo& info_p, CandsInfo& info_aniso) const
441  {
442  for (unsigned short i = 0; i < candidates.size(); i++)
443  {
444  CandsInfo* info = nullptr;
445  Cand& cand = candidates[i];
446  if (cand.split == H2D_REFINEMENT_H) info = &info_h;
447  else if (cand.split == H2D_REFINEMENT_P) info = &info_p;
448  else if (cand.split == H2D_REFINEMENT_H_ANISO_H || cand.split == H2D_REFINEMENT_H_ANISO_V) info = &info_aniso;
449  else { throw Hermes::Exceptions::Exception("Invalid candidate type: %d.", cand.split); };
450 
451  //evaluate elements of candidates
452  const int num_elems = cand.get_num_elems();
453  for (int i = 0; i < num_elems; i++)
454  {
455  int elem_order_h = H2D_GET_H_ORDER(cand.p[i]), elem_order_v = H2D_GET_V_ORDER(cand.p[i]);
456  if (elem_order_h != elem_order_v)
457  info->uniform_orders = false;
458  if (info->min_quad_order < 0 || info->max_quad_order < 0)
459  info->min_quad_order = info->max_quad_order = H2D_MAKE_QUAD_ORDER(elem_order_h, elem_order_v);
460  else
461  {
462  info->min_quad_order = H2D_MAKE_QUAD_ORDER(std::min(H2D_GET_H_ORDER(info->min_quad_order), elem_order_h), std::min(H2D_GET_V_ORDER(info->min_quad_order), elem_order_v));
463  info->max_quad_order = H2D_MAKE_QUAD_ORDER(std::max(H2D_GET_H_ORDER(info->max_quad_order), elem_order_h), std::max(H2D_GET_V_ORDER(info->max_quad_order), elem_order_v));
464  }
465  }
466  }
467  }
468 
469  template<typename Scalar>
470  void OptimumSelector<Scalar>::evaluate_cands_dof(std::vector<Cand>& candidates, Element* e, MeshFunction<Scalar>* rsln)
471  {
472  bool tri = e->is_triangle();
473 
474  for (unsigned i = 0; i < candidates.size(); i++)
475  {
476  Cand& c = candidates[i];
477  if (tri) { //triangle
478  switch (c.split)
479  {
480  case H2D_REFINEMENT_H:
481  {
482  //central triangle
483  const int central = 3;
484  c.dofs = 0;
485  for (int j = 0; j < H2D_MAX_ELEMENT_SONS; j++)
486  {
487  c.dofs += num_shapes[HERMES_MODE_TRIANGLE][H2D_GET_H_ORDER(c.p[j]) + 1][H2DRS_ORDER_ANY + 1][H2DSI_ANY];
488  if (j != central)
489  {
490  //shared edge: since triangle has three edges which are identified by a single order this will find 3 x different edge of a given order
491  c.dofs -= num_shapes[HERMES_MODE_TRIANGLE][std::min(H2D_GET_H_ORDER(c.p[j]), H2D_GET_H_ORDER(c.p[central])) + 1][H2DRS_ORDER_ANY + 1][H2DSI_TRI_EDGE] / 3;
492  }
493  }
494  if (has_vertex_shape[HERMES_MODE_TRIANGLE])
495  // Every vertex function belonging to vertices of the middle triangle is added 3-times, so it has to be deducted 2 times.
496  c.dofs -= 2 * 3;
497  }
498  break;
499 
500  case H2D_REFINEMENT_P:
501  c.dofs = num_shapes[HERMES_MODE_TRIANGLE][H2D_GET_H_ORDER(c.p[0]) + 1][H2DRS_ORDER_ANY + 1][H2DSI_ANY];
502  break;
503 
504  default:
505  throw Hermes::Exceptions::Exception("Unknown split type \"%d\" at candidate %d (element #%d)", c.split, i, e->id);
506  }
507  }
508  else { //quad
509  switch (c.split)
510  {
511  case H2D_REFINEMENT_H:
512  c.dofs = 0;
513  for (int j = 0; j < H2D_MAX_ELEMENT_SONS; j++)
514  c.dofs += num_shapes[HERMES_MODE_QUAD][H2D_GET_H_ORDER(c.p[j]) + 1][H2D_GET_V_ORDER(c.p[j]) + 1][H2DSI_ANY];
515  for (int j = 0; j < 2; j++) { //shared edge functions
516  //shared vertical edge functions: every edge is twice there
517  c.dofs -= num_shapes[HERMES_MODE_QUAD][H2DRS_ORDER_ANY + 1][std::min(H2D_GET_V_ORDER(c.p[2 * j]), H2D_GET_V_ORDER(c.p[2 * j + 1])) + 1][H2DSI_VERT_EDGE] / 2;
518  //shared horizontal edge functions: every edge is twice there
519  c.dofs -= num_shapes[HERMES_MODE_QUAD][std::min(H2D_GET_H_ORDER(c.p[j]), H2D_GET_H_ORDER(c.p[j ^ 3])) + 1][H2DRS_ORDER_ANY + 1][H2DSI_HORIZ_EDGE] / 2;
520  }
521  if (has_vertex_shape[HERMES_MODE_QUAD])
522  //edge vertex + central vertex
523  c.dofs -= 4 + 3;
524  break;
525 
527  c.dofs = num_shapes[HERMES_MODE_QUAD][H2D_GET_H_ORDER(c.p[0]) + 1][H2D_GET_V_ORDER(c.p[0]) + 1][H2DSI_ANY];
528  c.dofs += num_shapes[HERMES_MODE_QUAD][H2D_GET_H_ORDER(c.p[1]) + 1][H2D_GET_V_ORDER(c.p[1]) + 1][H2DSI_ANY];
529  //shared edge functions
530  c.dofs -= num_shapes[HERMES_MODE_QUAD][std::min(H2D_GET_H_ORDER(c.p[0]), H2D_GET_H_ORDER(c.p[1])) + 1][H2DRS_ORDER_ANY + 1][H2DSI_HORIZ_EDGE] / 2;
531  if (has_vertex_shape[HERMES_MODE_QUAD])
532  //shared vertex functions
533  c.dofs -= 2;
534  break;
535 
537  c.dofs = num_shapes[HERMES_MODE_QUAD][H2D_GET_H_ORDER(c.p[0]) + 1][H2D_GET_V_ORDER(c.p[0]) + 1][H2DSI_ANY];
538  c.dofs += num_shapes[HERMES_MODE_QUAD][H2D_GET_H_ORDER(c.p[1]) + 1][H2D_GET_V_ORDER(c.p[1]) + 1][H2DSI_ANY];
539  //shared edge functions
540  c.dofs -= num_shapes[HERMES_MODE_QUAD][H2DRS_ORDER_ANY + 1][std::min(H2D_GET_V_ORDER(c.p[0]), H2D_GET_V_ORDER(c.p[1])) + 1][H2DSI_VERT_EDGE] / 2;
541  if (has_vertex_shape[HERMES_MODE_QUAD])
542  //shared vertex functions
543  c.dofs -= 2;
544  break;
545 
546  case H2D_REFINEMENT_P:
547  c.dofs = num_shapes[HERMES_MODE_QUAD][H2D_GET_H_ORDER(c.p[0]) + 1][H2D_GET_V_ORDER(c.p[0]) + 1][H2DSI_ANY];
548  break;
549 
550  default:
551  throw Hermes::Exceptions::Exception("Unknown split type \"%d\" at candidate %d", c.split, i);
552  }
553  }
554  }
555  }
556 
557  template<typename Scalar>
558  void OptimumSelector<Scalar>::evaluate_candidates(std::vector<Cand>& candidates, Element* e, MeshFunction<Scalar>* rsln)
559  {
560  evaluate_cands_error(candidates, e, rsln);
561 
562  evaluate_cands_dof(candidates, e, rsln);
563 
564  evaluate_cands_score(candidates, e);
565  }
566 
567  template<typename Scalar>
569  {
570  this->dof_score_exponent = exponent;
571  }
572 
573  template<typename Scalar>
574  void OptimumSelector<Scalar>::evaluate_cands_score(std::vector<Cand>& candidates, Element* e)
575  {
576  // Original candidate.
577  Cand& unrefined = candidates[0];
578  // Original candidate score is zero.
579  unrefined.score = 0;
580 
581  for (unsigned short i = 1; i < candidates.size(); i++)
582  {
583  Cand& candidate = candidates[i];
584 
585  // We are only interested in candidates decreasing the error.
586  if (candidate.error < unrefined.error)
587  candidate.score = (log(unrefined.error / candidate.error)) / std::pow(candidate.dofs - unrefined.dofs, this->dof_score_exponent);
588  else
589  candidate.score = 0;
590  }
591  }
592 
593  template<typename Scalar>
595  {
596  return a.score > b.score;
597  }
598 
599  template<typename Scalar>
600  void OptimumSelector<Scalar>::select_best_candidate(std::vector<Cand>& candidates, Element* e, Cand*& best_candidate, Cand* best_candidates_specific_type[4])
601  {
602  //sort according to the score
603  const int num_cands = (int)candidates.size();
604 
605  if (num_cands > 2)
606  std::sort(candidates.begin() + 1, candidates.end(), compare_cand_score);
607 
608  // Overall best candidate.
609  if (candidates[1].score == 0)
610  // This means that the best candidate is the 'do-nothing' one.
611  best_candidate = &candidates[0];
612  else
613  // The one with the highest score is the best.
614  best_candidate = &candidates[1];
615 
616  for (int i = 0; i < num_cands; i++)
617  {
618  if (candidates[i].split == H2D_REFINEMENT_P)
619  {
620  best_candidates_specific_type[H2D_REFINEMENT_P] = &candidates[i];
621  break;
622  }
623  }
624 
625  for (int i = 0; i < num_cands; i++)
626  {
627  if (candidates[i].split == H2D_REFINEMENT_H)
628  {
629  best_candidates_specific_type[H2D_REFINEMENT_H] = &candidates[i];
630  break;
631  }
632  }
633 
634  for (int i = 0; i < num_cands; i++)
635  {
636  if (candidates[i].split == H2D_REFINEMENT_H_ANISO_H)
637  {
638  best_candidates_specific_type[H2D_REFINEMENT_H_ANISO_H] = &candidates[i];
639  break;
640  }
641  }
642 
643  for (int i = 0; i < num_cands; i++)
644  {
645  if (candidates[i].split == H2D_REFINEMENT_H_ANISO_V)
646  {
647  best_candidates_specific_type[H2D_REFINEMENT_H_ANISO_V] = &candidates[i];
648  break;
649  }
650  }
651  }
652 
653  template<typename Scalar>
655  {
656  //make an uniform order in a case of a triangle
657  int order_h = H2D_GET_H_ORDER(quad_order), order_v = H2D_GET_V_ORDER(quad_order);
658  if (element->is_triangle())
659  {
660  order_v = order_h;
661  //in a case of a triangle, order_v is zero. Set it to order_h in order to simplify the routines.
662  quad_order = H2D_MAKE_QUAD_ORDER(order_h, order_v);
663  }
664 
665  //build candidates.
666  std::vector<Cand> candidates = create_candidates(element, quad_order);
667  //there are candidates to choose from
668  Cand* best_candidate;
669  Cand* best_candidates_specific_type[4];
670  memset(best_candidates_specific_type, 0, 4 * sizeof(Cand*));
671 
672  if (candidates.size() > 1)
673  {
674  // evaluate candidates (sum partial projection errors, calculate dofs)
675  evaluate_candidates(candidates, element, rsln);
676 
677  //select candidate
678  select_best_candidate(candidates, element, best_candidate, best_candidates_specific_type);
679  }
680 
681  //there is not candidate to choose from, select the original candidate
682  else
683  {
684  best_candidate = &candidates[0];
685  }
686 
687  if (best_candidate == &candidates[0])
688  return false;
689 
690  //copy result to output
691  refinement.split = best_candidate->split;
692  ElementToRefine::copy_orders(refinement.refinement_polynomial_order, best_candidate->p);
693  for (int i = 1; i < 4; i++)
694  if (best_candidates_specific_type[i] != nullptr)
695  ElementToRefine::copy_orders(refinement.best_refinement_polynomial_order_type[i], best_candidates_specific_type[i]->p);
696 
697  ElementToRefine::copy_errors(refinement.errors, best_candidate->errors);
698 
699  //modify orders in a case of a triangle such that order_v is zero
700  if (element->is_triangle())
701  {
702  for (int i = 0; i < H2D_MAX_ELEMENT_SONS; i++)
703  {
704  assert(H2D_GET_V_ORDER(refinement.refinement_polynomial_order[i]) == 0 || H2D_GET_H_ORDER(refinement.refinement_polynomial_order[i]) == H2D_GET_V_ORDER(refinement.refinement_polynomial_order[i]));
705 
706  refinement.refinement_polynomial_order[i] = H2D_MAKE_QUAD_ORDER(H2D_GET_H_ORDER(refinement.refinement_polynomial_order[i]), 0);
707 
708  for (int poly_order_type_i = 1; poly_order_type_i < 4; poly_order_type_i++)
709  refinement.best_refinement_polynomial_order_type[poly_order_type_i][i] = H2D_MAKE_QUAD_ORDER(H2D_GET_H_ORDER(refinement.best_refinement_polynomial_order_type[poly_order_type_i][i]), 0);
710  }
711  }
712 
713  return true;
714  }
715 
716  template class HERMES_API OptimumSelector < double > ;
717  template class HERMES_API OptimumSelector < std::complex<double> > ;
718  }
719  }
720 }
void set_dof_score_exponent(double exponent)
Set the score DOF exponent.
RefinementType split
A refinement, see the enum RefinementType.
Definition: candidates.h:90
Definition: adapt.h:24
int id
element id number
Definition: element.h:112
P-candidates only. Hermes::Orders are modified non-uniformly.
Definition: candidates.h:49
virtual void select_best_candidate(std::vector< Cand > &candidates, Element *e, Cand *&best_candidate, Cand *best_candidates_specific_type[4])
Sorts and selects the best candidate and the best H-candidate according to the score.
Stores one element of a mesh.
Definition: element.h:107
virtual void evaluate_cands_dof(std::vector< Cand > &candidates, Element *e, MeshFunction< Scalar > *rsln)
Calculates DOF of candidates.
virtual std::vector< Cand > create_candidates(Element *e, int quad_order)
Fill a list of candidates.
void build_shape_indices(const ElementMode2D mode, const Range &vertex_order, const Range &edge_bubble_order)
Builds shape index table OptimumSelector::shape_indices.
bool uniform_orders
True if all elements of all examined candidates have uniform orders.
void update_cands_info(std::vector< Cand > &candidates, CandsInfo &info_h, CandsInfo &info_p, CandsInfo &info_aniso) const
Updates information about candidates. Initial information is provided.
#define H2DRS_MAX_ORDER
A maximum order suported by refinement selectors.
Definition: global.h:104
H-candidates only. Hermes::Orders are not modified.
Definition: candidates.h:50
Represents a function defined on a mesh.
Definition: mesh_function.h:56
double errors[H2D_MAX_ELEMENT_SONS]
Error of this candidate's sons.
Definition: candidates.h:86
double error
Error of this candidate's sons.
Definition: candidates.h:84
RefinementType
Possible refinements of an element.
static bool compare_cand_score(const Cand &a, const Cand &b)
Compares scores. Used to sort scores ascending.
#define H2DRS_MAX_ORDER_INC
Maximum increase of an order in candidates.
Definition: global.h:69
H-, ANISO- and P-candidates. Hermes::Orders are modified non-uniformly.
Definition: candidates.h:55
unsigned short refinement_polynomial_order[H2D_MAX_ELEMENT_SONS]
Encoded orders of sons.
H- and ANISO-candidates only. Hermes::Orders are not modified.
Definition: candidates.h:51
CandList
Predefined list of candidates.
Definition: candidates.h:46
int min_quad_order
Minimum quad order of all elements of all examined candidates.
ANISO-refienement. The element is split along the vertical axis. Quadrilaterals only.
static void copy_orders(unsigned short *dest, const unsigned short *src)
Copies array of orders.
#define H2DRS_ASSUMED_MAX_CANDS
An estimated maximum number of candidates. Used for purpose of reserving space.
Definition: global.h:66
unsigned short iro_cache
Increase in integration order, see RefMap::calc_inv_ref_order()
Definition: element.h:199
H- and P-candidates only. Hermes::Orders are modified uniformly.
Definition: candidates.h:52
RefinementType split
Proposed refinement. Possible values are defined in the enum RefinementType.
#define H2D_GET_H_ORDER(encoded_order)
Macros for combining quad horizontal and vertical encoded_orders.
Definition: global.h:98
OptimumSelector(CandList cand_list, int max_order, Shapeset *shapeset, const Range &vertex_order, const Range &edge_bubble_order)
Constructor.
unsigned short dofs
An estimated number of DOFs.
Definition: candidates.h:88
Should be exactly the same as is the count of enum ShapesetType.
Definition: shapeset.h:95
virtual void evaluate_cands_score(std::vector< Cand > &candidates, Element *e)
Evalutes score of candidates.
Hermes::Order permutator. Generates all permutations of orders from a set defined by a range of order...
virtual bool select_refinement(Element *element, int quad_order, MeshFunction< Scalar > *rsln, ElementToRefine &refinement)
Selects a refinement.
#define H2DRS_ORDER_ANY
Any order. Used as a wildcard to indicate that a given order can by any valid order.
Definition: global.h:71
A selector that chooses an optimal candidates based on a score.
Definition: function.h:32
HERMES_API int get_refin_sons(const RefinementType refin_type)
Returns a maximum number of sons that will be generated if a given refinement is applied.
unsigned short p[H2D_MAX_ELEMENT_SONS]
Encoded orders of sons, see ::H2D_MAKE_QUAD_ORDER. In a case of a triangle, the vertical order is equ...
Definition: candidates.h:92
int max_quad_order
Maximum quad order of all elements of all examined candidates. If less than zero, no candidate is gen...
double errors[H2D_MAX_ELEMENT_SONS]
Error of the selected candidate.
ANISO-refienement. The element is split along the horizontal axis. Quadrilaterals only...
A parent of all refinement selectors. Abstract class.
Definition: function.h:29
#define H2D_MAX_ELEMENT_SONS
Macros.
Definition: global.h:30
int calc_num_shapes(int mode, int order_h, int order_v, int allowed_type_mask)
Returns a number of shapes that may be contained in an element of a given order.
void reset()
Resets permutator to the starting order.
void add_bubble_shape_index(int order_h, int order_v, std::map< int, bool > &used_shape_index, std::vector< ShapeInx > &indices, ElementMode2D mode)
Adds an index (or indices) of a bubble function of a given order if the shape index was not used yet...
H-, ANISO- and P-candidates. Hermes::Orders are modified uniformly.
Definition: candidates.h:53
double score
A score of a candidate: the higher the better. If zero, the score is not valid and a candidate should...
Definition: candidates.h:94
unsigned char get_num_elems() const
Returns a number of elements of a candidate.
Definition: candidates.cpp:95
P-candidates only. Hermes::Orders are modified uniformly.
Definition: candidates.h:48
void append_candidates_split(std::vector< Cand > &candidates, const int start_quad_order, const int last_order, const RefinementType split, bool iso_p)
Appends cancidates of a given refinement and a given range of orders.
unsigned short best_refinement_polynomial_order_type[4][H2D_MAX_ELEMENT_SONS]
void evaluate_candidates(std::vector< Cand > &candidates, Element *e, MeshFunction< Scalar > *rsln)
Calculates error, dofs, and score of candidates.
#define H2D_NUM_SHAPES_SIZE
A maximum order suported by refinement selectors.
Definition: global.h:105