Hermes2D  2.0
optimum_selector.cpp
1 #include <iostream>
2 #include "global.h"
3 #include "solution.h"
4 #include "discrete_problem.h"
5 #include "quad_all.h"
6 #include "element_to_refine.h"
7 #include "optimum_selector.h"
8 
9 #define H2DST_ANY H2DST_VERTEX | H2DST_HORIZ_EDGE | H2DST_VERT_EDGE | H2DST_TRI_EDGE | H2DST_BUBBLE
10 
11 namespace Hermes
12 {
13  namespace Hermes2D
14  {
15  namespace RefinementSelectors
16  {
17  HERMES_API const char* get_cand_list_str(const CandList cand_list)
18  {
19  switch(cand_list)
20  {
21  case H2D_P_ISO: return "P_ISO";
22  case H2D_P_ANISO: return "P_ANISO";
23  case H2D_H_ISO: return "H_ISO";
24  case H2D_H_ANISO: return "H_ANISO";
25  case H2D_HP_ISO: return "HP_ISO";
26  case H2D_HP_ANISO_H: return "HP_ANISO_H";
27  case H2D_HP_ANISO_P: return "HP_ANISO_P";
28  case H2D_HP_ANISO: return "HP_ANISO";
29  default:
30  throw Hermes::Exceptions::Exception("Invalid adapt type %d.", cand_list);
31  return NULL;
32  }
33  }
34 
35  HERMES_API bool is_hp(const CandList cand_list)
36  {
37  switch(cand_list)
38  {
39  case H2D_P_ISO:
40  case H2D_P_ANISO:
41  case H2D_H_ISO:
42  case H2D_H_ANISO: return false; break;
43  case H2D_HP_ISO:
44  case H2D_HP_ANISO_H:
45  case H2D_HP_ANISO_P:
46  case H2D_HP_ANISO: return true; break;
47  default: throw Hermes::Exceptions::Exception("Invalid adapt type %d.", cand_list); return false;
48  }
49  }
50 
51  HERMES_API bool is_p_aniso(const CandList cand_list)
52  {
53  switch(cand_list)
54  {
55  case H2D_P_ISO: return false;
56  case H2D_P_ANISO: return true;
57  case H2D_H_ISO: return false;
58  case H2D_H_ANISO: return false;
59  case H2D_HP_ISO: return false;
60  case H2D_HP_ANISO_H: return false;
61  case H2D_HP_ANISO_P: return true;
62  case H2D_HP_ANISO: return true;
63  default: throw Hermes::Exceptions::Exception("Invalid adapt type %d.", cand_list); return false;
64  }
65  }
66 
67  template<typename Scalar>
68  OptimumSelector<Scalar>::OptimumSelector(CandList cand_list, double conv_exp, int
69  max_order, Shapeset* shapeset, const Range& vertex_order, const
70  Range& edge_bubble_order) :
71  Selector<Scalar>(max_order),
72  opt_symmetric_mesh(true),
73  opt_apply_exp_dof(false),
74  cand_list(cand_list),
75  conv_exp(conv_exp),
76  shapeset(shapeset)
77  {
78  if(shapeset == NULL)
79  throw Exceptions::NullException(3);
80 
81  num_shapes = new int***[2];
82  for(int i = 0; i < 2; i++)
83  {
84  num_shapes[i] = new int**[H2DRS_MAX_ORDER + 2];
85  for(int j = 0; j < H2DRS_MAX_ORDER + 2; j++)
86  {
87  num_shapes[i][j] = new int*[H2DRS_MAX_ORDER + 2];
88  for(int k = 0; k < H2DRS_MAX_ORDER + 2; k++)
89  {
90  num_shapes[i][j][k] = new int[6];
91  memset(num_shapes[i][j][k], 0, 6*sizeof(int));
92  }
93  }
94  }
95 
96  //build shape indices
97  build_shape_indices(HERMES_MODE_TRIANGLE, vertex_order, edge_bubble_order);
98  build_shape_indices(HERMES_MODE_QUAD, vertex_order, edge_bubble_order);
99  }
100 
101  template<typename Scalar>
102  OptimumSelector<Scalar>::~OptimumSelector()
103  {
104  if(!this->isAClone)
105  {
106  for(int i = 0; i < 2; i++)
107  {
108  for(int j = 0; j < H2DRS_MAX_ORDER + 2; j++)
109  {
110  for(int k = 0; k < H2DRS_MAX_ORDER + 2; k++)
111  {
112  delete [] num_shapes[i][j][k];
113  }
114  delete [] num_shapes[i][j];
115  }
116  delete [] num_shapes[i];
117  }
118  delete [] num_shapes;
119  }
120  }
121 
122  template<typename Scalar>
123  void OptimumSelector<Scalar>::add_bubble_shape_index(int order_h, int order_v, std::map<int, bool>& used_shape_index, Hermes::vector<ShapeInx>& indices, ElementMode2D mode)
124  {
125  int quad_order = H2D_MAKE_QUAD_ORDER(order_h, order_v);
126  const int num_bubbles = shapeset->get_num_bubbles(quad_order, mode);
127  int* bubble_inxs = shapeset->get_bubble_indices(quad_order, mode);
128  for(int j = 0; j < num_bubbles; j++)
129  {
130  int inx_bubble = bubble_inxs[j];
131  if(used_shape_index.find(inx_bubble) == used_shape_index.end())
132  {
133  used_shape_index[inx_bubble] = true;
134  indices.push_back(ShapeInx(order_h, order_v, inx_bubble, H2DST_BUBBLE));
135  for(int order_h_i = order_h+1; order_h_i < H2DRS_MAX_ORDER + 2; order_h_i++)
136  for(int order_v_i = order_v+1; order_v_i < H2DRS_MAX_ORDER + 2; order_v_i++)
137  {
138  num_shapes[mode][order_h_i][order_v_i][H2DSI_BUBBLE]++;
139  num_shapes[mode][order_h_i][order_v_i][H2DSI_ANY]++;
140  }
141  num_shapes[mode][0][0][H2DSI_BUBBLE]++;
142  num_shapes[mode][0][0][H2DSI_ANY]++;
143  }
144  }
145  }
146 
147  template<typename Scalar>
148  void OptimumSelector<Scalar>::build_shape_indices(const ElementMode2D mode, const Range& vertex_order, const Range& edge_bubble_order)
149  {
150  Hermes::vector<ShapeInx> &indices = shape_indices[mode];
151  int* next_order = this->next_order_shape[mode];
152  int& max_shape_inx = this->max_shape_inx[mode];
153  int num_edges = (mode == HERMES_MODE_QUAD) ? 4 : 3;
154  bool &has_vertex = has_vertex_shape[mode];
155  bool &has_edge = has_edge_shape[mode];
156  bool &has_bubble = has_bubble_shape[mode];
157 
158  //cleanup
159  indices.clear();
160  indices.reserve((H2DRS_MAX_ORDER + 1) * (H2DRS_MAX_ORDER + 1));
161  has_vertex = has_edge = has_bubble = false;
162 
163  //get total range of orders
164  Range order_range = Range::make_envelope(vertex_order, edge_bubble_order);
165 
166  //prepare array of used shape indices
167  std::map<int, bool> used_shape_index;
168 
169  //for all orders
170  max_shape_inx = 0;
171  int examined_shape = 0;
172  for(int i = order_range.lower(); i <= order_range.upper(); i++)
173  {
174  //vertex functions
175  if(vertex_order.is_in_closed(i))
176  {
177  for (int j = 0; j < num_edges; j++)
178  {
179  int inx = shapeset->get_vertex_index(j, mode);
180  if(inx >= 0)
181  {
182  used_shape_index[inx] = true;
183  indices.push_back(ShapeInx(1, 1, inx, H2DST_VERTEX));
184  if(mode == HERMES_MODE_QUAD)
185  {
186  for(int order_h_i = 2; order_h_i < H2DRS_MAX_ORDER + 2; order_h_i++)
187  for(int order_v_i = 2; order_v_i < H2DRS_MAX_ORDER + 2; order_v_i++)
188  {
189  num_shapes[mode][order_h_i][order_v_i][H2DSI_VERTEX]++;
190  num_shapes[mode][order_h_i][order_v_i][H2DSI_ANY]++;
191  }
192  num_shapes[mode][0][0][H2DSI_VERTEX]++;
193  num_shapes[mode][0][0][H2DSI_ANY]++;
194  }
195  else
196  {
197  for(int order_h_i = 1; order_h_i < H2DRS_MAX_ORDER + 1; order_h_i++)
198  {
199  num_shapes[mode][order_h_i][0][H2DSI_VERTEX]++;
200  num_shapes[mode][order_h_i][0][H2DSI_ANY]++;
201  }
202  num_shapes[mode][0][0][H2DSI_VERTEX]++;
203  num_shapes[mode][0][0][H2DSI_ANY]++;
204  }
205  has_vertex = true;
206  }
207  }
208  }
209 
210  //edge functions
211  if(edge_bubble_order.is_in_closed(i))
212  {
213  //edge functions
214  if(mode == HERMES_MODE_QUAD)
215  {
216  for (int j = 0; j < num_edges; j++)
217  {
218  int inx = shapeset->get_edge_index(j, 0, i, HERMES_MODE_QUAD);
219  if(inx >= 0)
220  {
221  used_shape_index[inx] = true;
222  if((j&1) == 0)//horizontal edge
223  {
224  indices.push_back(ShapeInx(i, 0, inx, H2DST_HORIZ_EDGE));
225  for(int order_h_i = i+1; order_h_i < H2DRS_MAX_ORDER + 2; order_h_i++)
226  for(int order_v_i = 1; order_v_i < H2DRS_MAX_ORDER + 2; order_v_i++)
227  {
228  num_shapes[mode][order_h_i][order_v_i][H2DST_HORIZ_EDGE]++;
229  num_shapes[mode][order_h_i][order_v_i][H2DSI_ANY]++;
230  }
231  num_shapes[mode][0][0][H2DST_HORIZ_EDGE]++;
232  num_shapes[mode][0][0][H2DSI_ANY]++;
233  }
234  else //vertical edge
235  {
236  indices.push_back(ShapeInx(0, i, inx, H2DST_VERT_EDGE));
237  for(int order_h_i = 1; order_h_i < H2DRS_MAX_ORDER + 2; order_h_i++)
238  for(int order_v_i = i+1; order_v_i < H2DRS_MAX_ORDER + 2; order_v_i++)
239  {
240  num_shapes[mode][order_h_i][order_v_i][H2DSI_VERT_EDGE]++;
241  num_shapes[mode][order_h_i][order_v_i][H2DSI_ANY]++;
242  }
243  num_shapes[mode][0][0][H2DSI_VERT_EDGE]++;
244  num_shapes[mode][0][0][H2DSI_ANY]++;
245  }
246  has_edge = true;
247  }
248  }
249  }
250  else
251  {
252  for (int j = 0; j < num_edges; j++)
253  {
254  int inx = shapeset->get_edge_index(j, 0, i, HERMES_MODE_TRIANGLE);
255  if(inx >= 0)
256  {
257  used_shape_index[inx] = true;
258  indices.push_back(ShapeInx(i, i, inx, H2DST_TRI_EDGE));
259  for(int order_h_i = i+1; order_h_i < H2DRS_MAX_ORDER + 2; order_h_i++)
260  {
261  num_shapes[mode][order_h_i][0][H2DSI_TRI_EDGE]++;
262  num_shapes[mode][order_h_i][0][H2DSI_ANY]++;
263  }
264  num_shapes[mode][0][0][H2DSI_TRI_EDGE]++;
265  num_shapes[mode][0][0][H2DSI_ANY]++;
266  has_edge = true;
267  }
268  }
269  }
270 
271  //bubble functions
272  if(mode == HERMES_MODE_QUAD)
273  {
274  //NOTE: shapeset returns a set of all possible bubble functions and it is not possible to identify which is the smallest
275  // 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
276  // a list that contains a function of poly-order 2/0. Also, order of indices is not given.
277  int order = i;
278  unsigned num_indices_prev = indices.size();
279  for(int order_other = edge_bubble_order.lower(); order_other <= order; order_other++)
280  {
281  add_bubble_shape_index(order, order_other, used_shape_index, indices, mode);
282  add_bubble_shape_index(order_other, order, used_shape_index, indices, mode);
283  }
284 
285  //check if indices were added
286  if(num_indices_prev < indices.size())
287  has_bubble = true;
288  }
289  else { //triangles
290  int order = i;
291  int num_bubbles = shapeset->get_num_bubbles(order, mode);
292  int* bubble_inxs = shapeset->get_bubble_indices(order, mode);
293  for(int j = 0; j < num_bubbles; j++)
294  {
295  int inx_bubble = bubble_inxs[j];
296  if(used_shape_index.find(inx_bubble) == used_shape_index.end())
297  {
298  used_shape_index[inx_bubble] = true;
299  indices.push_back(ShapeInx(order, order, inx_bubble, H2DST_BUBBLE));
300  for(int order_h_i = order+1; order_h_i < H2DRS_MAX_ORDER + 2; order_h_i++)
301  {
302  num_shapes[mode][order_h_i][0][H2DSI_BUBBLE]++;
303  num_shapes[mode][order_h_i][0][H2DSI_ANY]++;
304  }
305  num_shapes[mode][0][0][H2DSI_BUBBLE]++;
306  num_shapes[mode][0][0][H2DSI_ANY]++;
307  has_bubble = true;
308  }
309  }
310  }
311 
312  //store index of the next order
313  next_order[i] = (int)indices.size();
314 
315  //update maximum
316  while(examined_shape < next_order[i])
317  {
318  max_shape_inx = std::max(max_shape_inx, indices[examined_shape].inx);
319  examined_shape++;
320  }
321  }
322  else
323  next_order[i] = (int)indices.size();
324  }
325  }
326 
327  template<typename Scalar>
328  int OptimumSelector<Scalar>::calc_num_shapes(int mode, int order_h, int order_v, int allowed_type_mask)
329  {
330  //test whether the evaluation is necessary
331  bool full_eval = false;
332  if((allowed_type_mask & H2DST_VERTEX) != 0)
333  full_eval |= has_vertex_shape[mode];
334  if((allowed_type_mask & (H2DST_HORIZ_EDGE | H2DST_VERT_EDGE | H2DST_TRI_EDGE)) != 0)
335  full_eval |= has_edge_shape[mode];
336  if((allowed_type_mask & H2DST_BUBBLE) != 0)
337  full_eval |= has_bubble_shape[mode];
338 
339  //evaluate
340  if(full_eval)
341  {
342  Hermes::vector<ShapeInx>& shapes = shape_indices[mode];
343  int num = 0;
344  typename Hermes::vector<ShapeInx>::const_iterator shape = shapes.begin();
345  while (shape != shapes.end())
346  {
347  if(((int)shape->type & allowed_type_mask) != 0)
348  {
349  if((order_h == H2DRS_ORDER_ANY || shape->order_h <= order_h) && (order_v == H2DRS_ORDER_ANY || shape->order_v <= order_v))
350  num++;
351  }
352  shape++;
353  }
354  return num;
355  }
356  else
357  return 0;
358  }
359 
360  template<typename Scalar>
361  void OptimumSelector<Scalar>::append_candidates_split(const int start_quad_order, const int last_quad_order, const int split, bool iso_p)
362  {
363  //check whether end orders are not lower than start orders
364  if(last_quad_order < 0 || start_quad_order < 0)
365  return;
366  if(H2D_GET_H_ORDER(start_quad_order) > H2D_GET_H_ORDER(last_quad_order) || H2D_GET_V_ORDER(start_quad_order) > H2D_GET_V_ORDER(last_quad_order))
367  return;
368 
369  //get number of sons
370  const int num_sons = get_refin_sons(split);
371 
372  //initialize orders
373  int quad_orders[H2D_MAX_ELEMENT_SONS];
374  OrderPermutator quad_perms[H2D_MAX_ELEMENT_SONS];
375  for(int i = 0; i < num_sons; i++)
376  {
377  quad_orders[i] = start_quad_order;
378  quad_perms[i] = OrderPermutator(start_quad_order, last_quad_order, iso_p, &quad_orders[i]);
379  }
380  for(int i = num_sons; i < H2D_MAX_ELEMENT_SONS; i++)
381  quad_orders[i] = 0;
382 
383  //generate permutations of orders
384  bool quit = false;
385  while(!quit)
386  {
387  do { //create permutation of son 0
388  candidates.push_back(Cand(split, quad_orders));
389  } while (quad_perms[0].next());
390 
391  //reset son 0
392  quad_perms[0].reset();
393 
394  //increment orders of other sons
395  int inx_son = 1;
396  while (inx_son < num_sons && !quad_perms[inx_son].next())
397  {
398  quad_perms[inx_son].reset(); //reset order of the son
399  inx_son++;
400  }
401  if(inx_son >= num_sons)
402  quit = true;
403  }
404  }
405 
406  template<typename Scalar>
407  void OptimumSelector<Scalar>::create_candidates(Element* e, int quad_order, int max_ha_quad_order, int max_p_quad_order)
408  {
409  int order_h = H2D_GET_H_ORDER(quad_order), order_v = H2D_GET_V_ORDER(quad_order);
410  int max_p_order_h = H2D_GET_H_ORDER(max_p_quad_order), max_p_order_v = H2D_GET_V_ORDER(max_p_quad_order);
411  int max_ha_order_h = H2D_GET_H_ORDER(max_ha_quad_order), max_ha_order_v = H2D_GET_V_ORDER(max_ha_quad_order);
412  bool tri = e->is_triangle();
413 
414  //clear list of candidates
415  candidates.clear();
416  if(candidates.capacity() < H2DRS_ASSUMED_MAX_CANDS)
417  candidates.reserve(H2DRS_ASSUMED_MAX_CANDS);
418 
419  //generate all P-candidates (start from intention of generating all possible candidates
420  //and restrict it according to the given adapt-type)
421  bool iso_p = false;
422  int start_order_h = std::max(current_min_order, order_h - 1), start_order_v = std::max(current_min_order, order_v - 1);
423  int start_quad_order = H2D_MAKE_QUAD_ORDER(start_order_h, start_order_v);
424  int last_quad_order = H2D_MAKE_QUAD_ORDER(std::min(max_p_order_h, order_h + H2DRS_MAX_ORDER_INC), std::min(max_p_order_v, order_v + H2DRS_MAX_ORDER_INC));
425  switch(cand_list)
426  {
427  case H2D_H_ISO:
428  case H2D_H_ANISO: last_quad_order = start_quad_order; break; //no P-candidates except the original candidate
429  case H2D_P_ISO:
430  case H2D_HP_ISO:
431  case H2D_HP_ANISO_H: iso_p = true; break; //iso change of orders
432  }
433  append_candidates_split(quad_order, last_quad_order, H2D_REFINEMENT_P, tri || iso_p);
434 
435  //generate all H-candidates
436  iso_p = false;
437  start_order_h = std::max(current_min_order, order_h - 1), start_order_v = std::max(current_min_order, order_v - 1);
438  start_quad_order = H2D_MAKE_QUAD_ORDER(start_order_h, start_order_v);
439  last_quad_order = H2D_MAKE_QUAD_ORDER(std::min(max_ha_order_h, std::min(start_order_h + H2DRS_MAX_ORDER_INC, order_h)), std::min(max_ha_order_v, std::min(start_order_v + H2DRS_MAX_ORDER_INC, order_v)));
440  switch(cand_list)
441  {
442  case H2D_H_ISO:
443  case H2D_H_ANISO:
444  last_quad_order = start_quad_order = quad_order; break; //no only one candidate will be created
445  case H2D_P_ISO:
446  case H2D_P_ANISO: last_quad_order = -1; break; //no H-candidate will be generated
447  case H2D_HP_ISO:
448  case H2D_HP_ANISO_H: iso_p = true; break; //iso change of orders
449  }
450  append_candidates_split(start_quad_order, last_quad_order, H2D_REFINEMENT_H, tri || iso_p);
451 
452  //generate all ANISO-candidates
453  if(!tri && e->iro_cache < 8
454  && (cand_list == H2D_H_ANISO || cand_list == H2D_HP_ANISO_H || cand_list == H2D_HP_ANISO))
455  {
456  iso_p = false;
457  int start_quad_order_hz = H2D_MAKE_QUAD_ORDER(order_h, std::max(current_min_order, (order_v - 1)));
458  int last_quad_order_hz = H2D_MAKE_QUAD_ORDER(std::min(max_ha_order_h, order_h + H2DRS_MAX_ORDER_INC), std::min(order_v, H2D_GET_V_ORDER(start_quad_order) + H2DRS_MAX_ORDER_INC));
459  int start_quad_order_vt = H2D_MAKE_QUAD_ORDER(std::max(current_min_order, (order_h - 1)), order_v);
460  int last_quad_order_vt = H2D_MAKE_QUAD_ORDER(std::min(order_h, H2D_GET_H_ORDER(start_quad_order) + H2DRS_MAX_ORDER_INC), std::min(max_ha_order_v, order_v + H2DRS_MAX_ORDER_INC));
461  switch(cand_list)
462  {
463  case H2D_H_ANISO:
464  last_quad_order_hz = start_quad_order_hz = quad_order;
465  last_quad_order_vt = start_quad_order_vt = quad_order;
466  break; //only one candidate will be created
467  case H2D_HP_ANISO_H: iso_p = true; break; //iso change of orders
468  }
469  if(iso_p) { //make orders uniform: take mininmum order since nonuniformity is caused by different handling of orders along directions
470  int order = std::min(H2D_GET_H_ORDER(start_quad_order_hz), H2D_GET_V_ORDER(start_quad_order_hz));
471  start_quad_order_hz = H2D_MAKE_QUAD_ORDER(order, order);
472  order = std::min(H2D_GET_H_ORDER(start_quad_order_vt), H2D_GET_V_ORDER(start_quad_order_vt));
473  start_quad_order_vt = H2D_MAKE_QUAD_ORDER(order, order);
474 
475  order = std::min(H2D_GET_H_ORDER(last_quad_order_hz), H2D_GET_V_ORDER(last_quad_order_hz));
476  last_quad_order_hz = H2D_MAKE_QUAD_ORDER(order, order);
477  order = std::min(H2D_GET_H_ORDER(last_quad_order_vt), H2D_GET_V_ORDER(last_quad_order_vt));
478  last_quad_order_vt = H2D_MAKE_QUAD_ORDER(order, order);
479  }
480  append_candidates_split(start_quad_order_hz, last_quad_order_hz, H2D_REFINEMENT_ANISO_H, iso_p);
481  append_candidates_split(start_quad_order_vt, last_quad_order_vt, H2D_REFINEMENT_ANISO_V, iso_p);
482  }
483  }
484 
485  template<typename Scalar>
486  void OptimumSelector<Scalar>::update_cands_info(CandsInfo& info_h, CandsInfo& info_p, CandsInfo& info_aniso) const
487  {
488  typename Hermes::vector<Cand>::const_iterator cand = candidates.begin();
489  while (cand != candidates.end())
490  {
491  CandsInfo* info = NULL;
492  if(cand->split == H2D_REFINEMENT_H) info = &info_h;
493  else if(cand->split == H2D_REFINEMENT_P) info = &info_p;
494  else if(cand->split == H2D_REFINEMENT_ANISO_H || cand->split == H2D_REFINEMENT_ANISO_V) info = &info_aniso;
495  else { throw Hermes::Exceptions::Exception("Invalid candidate type: %d.", cand->split); };
496 
497  //evaluate elements of candidates
498  const int num_elems = cand->get_num_elems();
499  for(int i = 0; i < num_elems; i++)
500  {
501  int elem_order_h = H2D_GET_H_ORDER(cand->p[i]), elem_order_v = H2D_GET_V_ORDER(cand->p[i]);
502  if(elem_order_h != elem_order_v)
503  info->uniform_orders = false;
504  if(info->min_quad_order < 0 || info->max_quad_order < 0)
505  info->min_quad_order = info->max_quad_order = H2D_MAKE_QUAD_ORDER(elem_order_h, elem_order_v);
506  else
507  {
508  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));
509  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));
510  }
511  }
512 
513  //next candidate
514  cand++;
515  }
516  }
517 
518  template<typename Scalar>
519  void OptimumSelector<Scalar>::evaluate_cands_dof(Element* e, Solution<Scalar>* rsln)
520  {
521  bool tri = e->is_triangle();
522 
523  for (unsigned i = 0; i < candidates.size(); i++)
524  {
525  Cand& c = candidates[i];
526  if(tri) { //triangle
527  switch(c.split)
528  {
529  case H2D_REFINEMENT_H:
530  {
531  const int central = 3; //central triangle
532  c.dofs = 0;
533  for(int j = 0; j < H2D_MAX_ELEMENT_SONS; j++)
534  {
535  c.dofs += num_shapes[HERMES_MODE_TRIANGLE][H2D_GET_H_ORDER(c.p[j]) + 1][H2DRS_ORDER_ANY + 1][H2DSI_ANY];
536  if(j != central)
537  {
538  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; //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
539  }
540  }
541  if(has_vertex_shape[HERMES_MODE_TRIANGLE])
542  c.dofs -= 2*3; // Every vertex function belonging to vertices of the middle triangle is added 3-times, so it has to be deducted 2 times.
543  }
544  break;
545 
546  case H2D_REFINEMENT_P:
547  c.dofs = num_shapes[HERMES_MODE_TRIANGLE][H2D_GET_H_ORDER(c.p[0]) + 1][H2DRS_ORDER_ANY + 1][H2DSI_ANY];
548  break;
549 
550  default:
551  throw Hermes::Exceptions::Exception("Unknown split type \"%d\" at candidate %d (element #%d)", c.split, i, e->id);
552  }
553  }
554  else { //quad
555  switch(c.split)
556  {
557  case H2D_REFINEMENT_H:
558  c.dofs = 0;
559  for(int j = 0; j < H2D_MAX_ELEMENT_SONS; j++)
560  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];
561  for(int j = 0; j < 2; j++) { //shared edge functions
562  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; //shared vertical edge functions: every edge is twice there
563  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; //shared horizontal edge functions: every edge is twice there
564  }
565  if(has_vertex_shape[HERMES_MODE_QUAD])
566  c.dofs -= 4 + 3; //edge vertex + central vertex
567  break;
568 
569  case H2D_REFINEMENT_ANISO_H:
570  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];
571  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];
572  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; //shared edge functions
573  if(has_vertex_shape[HERMES_MODE_QUAD])
574  c.dofs -= 2; //shared vertex functions
575  break;
576 
577  case H2D_REFINEMENT_ANISO_V:
578  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];
579  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];
580  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; //shared edge functions
581  if(has_vertex_shape[HERMES_MODE_QUAD])
582  c.dofs -= 2; //shared vertex functions
583  break;
584 
585  case H2D_REFINEMENT_P:
586  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];
587  break;
588 
589  default:
590  throw Hermes::Exceptions::Exception("Unknown split type \"%d\" at candidate %d", c.split, i);
591  }
592  }
593  }
594  }
595 
596  template<typename Scalar>
597  void OptimumSelector<Scalar>::evaluate_candidates(Element* e, Solution<Scalar>* rsln, double* avg_error, double* dev_error)
598  {
599  evaluate_cands_error(e, rsln, avg_error, dev_error);
600 
601  evaluate_cands_dof(e, rsln);
602 
603  evaluate_cands_score(e);
604  }
605 
606  template<typename Scalar>
607  void OptimumSelector<Scalar>::evaluate_cands_score(Element* e)
608  {
609  //calculate score of candidates
610  Cand& unrefined = candidates[0];
611  const int num_cands = (int)candidates.size();
612  unrefined.score = 0;
613  const double unrefined_dofs_exp = std::pow(unrefined.dofs, conv_exp);
614  for (int i = 1; i < num_cands; i++)
615  {
616  Cand& cand = candidates[i];
617  if(cand.error < unrefined.error && cand.dofs > unrefined.dofs)
618  {
619  double delta_dof_exp = std::pow(cand.dofs - unrefined.dofs, conv_exp);
620  if(opt_apply_exp_dof)
621  delta_dof_exp = std::pow(cand.dofs, conv_exp) - unrefined_dofs_exp;
622  candidates[i].score = (log10(unrefined.error) - log10(cand.error)) / delta_dof_exp;
623  }
624  else
625  candidates[i].score = 0;
626  }
627  }
628 
629  template<typename Scalar>
630  bool OptimumSelector<Scalar>::compare_cand_score(const Cand& a, const Cand& b)
631  {
632  return a.score > b.score;
633  }
634 
635  template<typename Scalar>
636  void OptimumSelector<Scalar>::select_best_candidate(Element* e, const double avg_error, const double dev_error, int* selected_cand, int* selected_h_cand)
637  {
638  // select an above-average candidate with the steepest error decrease
639 
640  //sort according to the score
641  const int num_cands = (int)candidates.size();
642  if(num_cands > 2)
643  std::sort(candidates.begin() + 1, candidates.end(), compare_cand_score);
644 
645  //select best candidate
646  int imax = 1, h_imax = 1;
647  if(opt_symmetric_mesh) { //prefer symmetric mesh
648  //find first valid score that diffres from the next scores
649  while ((imax + 1) < num_cands && std::abs(candidates[imax].score - candidates[imax + 1].score) < H2DRS_SCORE_DIFF_ZERO)
650  {
651  //find the first candidate with a different score
652  Cand& cand_current = candidates[imax];
653  int imax_end = imax + 2;
654  while (imax_end < num_cands && std::abs(cand_current.score - candidates[imax_end].score) < H2DRS_SCORE_DIFF_ZERO)
655  imax_end++;
656 
657  imax = imax_end;
658  }
659 
660  //find valid H-refinement candidate
661  h_imax = imax;
662  while (h_imax < num_cands && candidates[h_imax].split != H2D_REFINEMENT_H)
663  h_imax++;
664  }
665 
666  //make sure the result is valid: index is valid, selected candidate has a valid score
667  if(imax >= num_cands || candidates[imax].score == 0)
668  imax = 0;
669  if(h_imax >= num_cands || candidates[h_imax].score == 0)
670  h_imax = 0;
671 
672  //report result
673  *selected_cand = imax;
674  *selected_h_cand = h_imax;
675  }
676 
677  template<typename Scalar>
678  bool OptimumSelector<Scalar>::select_refinement(Element* element, int quad_order, Solution<Scalar>* rsln, ElementToRefine& refinement)
679  {
680  //make an uniform order in a case of a triangle
681  int order_h = H2D_GET_H_ORDER(quad_order), order_v = H2D_GET_V_ORDER(quad_order);
682  if(element->is_triangle())
683  {
684  order_v = order_h;
685  quad_order = H2D_MAKE_QUAD_ORDER(order_h, order_v); //in a case of a triangle, order_v is zero. Set it to order_h in order to simplify the routines.
686  }
687 
688  //set orders
689  set_current_order_range(element);
690 
691  // To generate always at least the unchanged candidate.
692  if(current_max_order < std::max(order_h, order_v))
693  current_max_order = std::max(order_h, order_v);
694 
695  //build candidates
696  int inx_cand, inx_h_cand;
697  create_candidates(element, quad_order
698  , H2D_MAKE_QUAD_ORDER(current_max_order, current_max_order)
699  , H2D_MAKE_QUAD_ORDER(current_max_order, current_max_order));
700 
701  if(candidates.size() > 1) { //there are candidates to choose from
702  // evaluate candidates (sum partial projection errors, calculate dofs)
703  double avg_error, dev_error;
704  evaluate_candidates(element, rsln, &avg_error, &dev_error);
705 
706  //select candidate
707  select_best_candidate(element, avg_error, dev_error, &inx_cand, &inx_h_cand);
708  }
709  else { //there is not candidate to choose from, select the original candidate
710  inx_cand = 0;
711  inx_h_cand = 0;
712  }
713 
714  //copy result to output
715  Cand& cand = candidates[inx_cand];
716  Cand& cand_h = candidates[inx_h_cand];
717  refinement.split = cand.split;
718  ElementToRefine::copy_orders(refinement.p, cand.p);
719  if(candidates[inx_h_cand].split == H2D_REFINEMENT_H) { //inx_h_cand points to a candidate which is a H-candidate: copy orders
720  ElementToRefine::copy_orders(refinement.q, cand_h.p);
721  }
722  else { //the index is not H-candidate because not candidate was generate: so, fake orders
723  int h_cand_orders[H2D_MAX_ELEMENT_SONS] = { cand_h.p[0], cand_h.p[0], cand_h.p[0], cand_h.p[0] };
724  ElementToRefine::copy_orders(refinement.q, h_cand_orders);
725  }
726 
727  //modify orders in a case of a triangle such that order_v is zero
728  if(element->is_triangle())
729  {
730  for(int i = 0; i < H2D_MAX_ELEMENT_SONS; i++)
731  {
732  if(!(H2D_GET_V_ORDER(refinement.p[i]) == 0 || H2D_GET_H_ORDER(refinement.p[i]) == H2D_GET_V_ORDER(refinement.p[i])))
733  throw Exceptions::Exception("Triangle processed but the resulting order (%d, %d) of son %d is not uniform", H2D_GET_H_ORDER(refinement.p[i]), H2D_GET_V_ORDER(refinement.p[i]), i);
734  refinement.p[i] = H2D_MAKE_QUAD_ORDER(H2D_GET_H_ORDER(refinement.p[i]), 0);
735  if(!(H2D_GET_V_ORDER(refinement.q[i]) == 0 || H2D_GET_H_ORDER(refinement.q[i]) == H2D_GET_V_ORDER(refinement.q[i])))
736  throw Exceptions::Exception("Triangle processed but the resulting q-order (%d, %d) of son %d is not uniform", H2D_GET_H_ORDER(refinement.q[i]), H2D_GET_V_ORDER(refinement.q[i]), i);
737  refinement.q[i] = H2D_MAKE_QUAD_ORDER(H2D_GET_H_ORDER(refinement.q[i]), 0);
738  }
739  }
740 
741  if(inx_cand == 0)
742  return false;
743  else
744  return true;
745  }
746 
747  template<typename Scalar>
748  void OptimumSelector<Scalar>::generate_shared_mesh_orders(const Element* element, const int orig_quad_order, const int refinement, int tgt_quad_orders[H2D_MAX_ELEMENT_SONS], const int* suggested_quad_orders)
749  {
750  if(refinement == H2D_REFINEMENT_P)
751  throw Exceptions::Exception("P-candidate not supported for updating shared orders");
752  const int num_sons = get_refin_sons(refinement);
753  if(suggested_quad_orders != NULL)
754  {
755  for(int i = 0; i < num_sons; i++)
756  tgt_quad_orders[i] = suggested_quad_orders[i];
757  }
758  else
759  {
760  //calculate new quad_orders
761  int quad_order = orig_quad_order;
762  if(cand_list != H2D_H_ISO && cand_list != H2D_H_ANISO) { //H_ISO and H_ANISO has to keep given order
763  int order_h = H2D_GET_H_ORDER(quad_order), order_v = H2D_GET_V_ORDER(quad_order);
764  switch(refinement)
765  {
766  case H2D_REFINEMENT_H:
767  order_h = std::max(1, (order_h + 1)/2);
768  order_v = std::max(1, (order_v + 1)/2);
769  break;
770  case H2D_REFINEMENT_ANISO_H:
771  order_v = std::max(1, 2*(order_v + 1)/3);
772  break;
773  case H2D_REFINEMENT_ANISO_V:
774  order_h = std::max(1, 2*(order_h + 1)/3);
775  break;
776  }
777  if(element->is_triangle())
778  quad_order = H2D_MAKE_QUAD_ORDER(order_h, 0);
779  else
780  quad_order = H2D_MAKE_QUAD_ORDER(order_h, order_v);
781  }
782  for(int i = 0; i < num_sons; i++)
783  tgt_quad_orders[i] = quad_order;
784  }
785  }
786 
787  template<typename Scalar>
788  void OptimumSelector<Scalar>::set_option(const SelOption option, bool enable)
789  {
790  switch(option)
791  {
792  case H2D_PREFER_SYMMETRIC_MESH: opt_symmetric_mesh = enable; break;
793  case H2D_APPLY_CONV_EXP_DOF: opt_apply_exp_dof = enable; break;
794  default: throw Hermes::Exceptions::Exception("Unknown option %d.", (int)option);
795  }
796  }
797 
798  template<typename Scalar>
799  OptimumSelector<Scalar>::Range::Range() : empty_range(true) {}
800 
801  template<typename Scalar>
802  OptimumSelector<Scalar>::Range::Range(const int& lower_bound, const int& upper_bound) : lower_bound(lower_bound), upper_bound(upper_bound), empty_range(lower_bound > upper_bound)
803  {
804  }
805 
806  template<typename Scalar>
807  bool OptimumSelector<Scalar>::Range::empty() const
808  {
809  return empty_range;
810  }
811 
812  template<typename Scalar>
813  const int& OptimumSelector<Scalar>::Range::lower() const
814  {
815  return lower_bound;
816  }
817 
818  template<typename Scalar>
819  const int& OptimumSelector<Scalar>::Range::upper() const
820  {
821  return upper_bound;
822  }
823 
824  template<typename Scalar>
825  bool OptimumSelector<Scalar>::Range::is_in_closed(const typename OptimumSelector<Scalar>::Range& range) const
826  {
827  return (range.lower_bound >= lower_bound && range.upper_bound <= upper_bound);
828  }
829 
830  template<typename Scalar>
831  bool OptimumSelector<Scalar>::Range::is_in_closed(const int& value) const
832  {
833  return (value >= lower_bound && value <= upper_bound);
834  }
835 
836  template<typename Scalar>
837  bool OptimumSelector<Scalar>::Range::is_in_open(const int& value) const
838  {
839  return (value > lower_bound && value < upper_bound);
840  }
841 
842  template<typename Scalar>
843  void OptimumSelector<Scalar>::Range::enlarge_to_include(const int& value)
844  {
845  if(empty_range) {
846  lower_bound = upper_bound = value;
847  empty_range = false;
848  }
849  else {
850  if(lower_bound > value)
851  lower_bound = value;
852  if(upper_bound < value)
853  upper_bound = value;
854  }
855  }
856 
857  template<typename Scalar>
858  typename OptimumSelector<Scalar>::Range OptimumSelector<Scalar>::Range::make_envelope(const typename OptimumSelector<Scalar>::Range& a, const typename OptimumSelector<Scalar>::Range& b)
859  {
860  if(a.empty())
861  return b;
862  else if(b.empty())
863  return a;
864  else
865  return Range(std::min(a.lower(), b.lower()), std::max(a.upper(), b.upper()));
866  }
867 
868  template class HERMES_API OptimumSelector<double>;
869  template class HERMES_API OptimumSelector<std::complex<double> >;
870  }
871  }
872 }