Hermes2D  3.0
proj_based_selector.cpp
1 // This file is part of Hermes2D.
2 //
3 // Hermes2D is free software: you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License as published by
5 // the Free Software Foundation, either version 2 of the License, or
6 // (at your option) any later version.
7 //
8 // Hermes2D is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 // GNU General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with Hermes2D. If not, see <http://www.gnu.org/licenses/>.
15 
16 #include "proj_based_selector.h"
17 #include "hcurl_proj_based_selector.h"
18 #include <algorithm>
19 #include "order_permutator.h"
20 #include "algebra/dense_matrix_operations.h"
21 
23 
24 namespace Hermes
25 {
26  namespace Hermes2D
27  {
28  namespace RefinementSelectors
29  {
30  template<typename Scalar>
32  max_order, Shapeset* shapeset, const Range& vertex_order, const
33  Range& edge_bubble_order) :
34  OptimumSelector<Scalar>(cand_list, max_order, shapeset, vertex_order, edge_bubble_order),
35  warn_uniform_orders(false),
36  error_weight_h(H2DRS_DEFAULT_ERR_WEIGHT_H),
37  error_weight_p(H2DRS_DEFAULT_ERR_WEIGHT_P),
38  error_weight_aniso(H2DRS_DEFAULT_ERR_WEIGHT_ANISO)
39  {
40  cached_shape_vals_valid = new bool[2];
42  cached_shape_vals = new TrfShape[2];
43 
44  //clean svals initialization state
46 
47  //clear matrix cache
48  for (int m = 0; m < H2D_NUM_MODES; m++)
49  for (int i = 0; i < H2DRS_MAX_ORDER + 2; i++)
50  for (int k = 0; k < H2DRS_MAX_ORDER + 2; k++)
51  proj_matrix_cache[m][i][k] = nullptr;
52  }
53 
54  template<typename Scalar>
56  {
57  //delete matrix cache
58  for (int m = 0; m < H2D_NUM_MODES; m++)
59  {
60  for (int i = 0; i < H2DRS_MAX_ORDER + 2; i++)
61  for (int k = 0; k < H2DRS_MAX_ORDER + 2; k++)
62  {
63  if (proj_matrix_cache[m][i][k] != nullptr)
64  free_with_check<double*>(proj_matrix_cache[m][i][k], true);
65  }
66  }
67 
68  delete[] cached_shape_vals_valid;
69  delete[] cached_shape_ortho_vals;
70  delete[] cached_shape_vals;
71  }
72 
73  template<typename Scalar>
74  void ProjBasedSelector<Scalar>::set_error_weights(double weight_h, double weight_p, double weight_aniso)
75  {
76  error_weight_h = weight_h;
77  error_weight_p = weight_p;
78  error_weight_aniso = weight_aniso;
79  }
80 
81  template<typename Scalar>
83  {
84  return error_weight_h;
85  }
86 
87  template<typename Scalar>
88  double ProjBasedSelector<Scalar>::get_error_weight_p() const
89  {
90  return error_weight_p;
91  }
92 
93  template<typename Scalar>
94  double ProjBasedSelector<Scalar>::get_error_weight_aniso() const
95  {
96  return error_weight_aniso;
97  }
98 
99  template<typename Scalar>
100  ProjBasedSelector<Scalar>::TrfShapeExp::TrfShapeExp() : num_gip(0), num_expansion(0), values(nullptr) {};
101 
102  template<typename Scalar>
104  {
105  free_with_check(values, true);
106  }
107 
108  template<typename Scalar>
109  void ProjBasedSelector<Scalar>::TrfShapeExp::allocate(int num_expansion, int num_gip)
110  {
111  free_with_check(values, true);
112  values = new_matrix<double>(num_expansion, num_gip);
113  this->num_expansion = num_expansion;
114  this->num_gip = num_gip;
115  }
116 
117  template<typename Scalar>
118  double* ProjBasedSelector<Scalar>::TrfShapeExp::operator[](int inx_expansion)
119  {
120  return values[inx_expansion];
121  }
122 
123  template<typename Scalar>
124  bool ProjBasedSelector<Scalar>::TrfShapeExp::empty() const
125  {
126  return values == nullptr;
127  }
128 
129  template<typename Scalar>
130  void ProjBasedSelector<Scalar>::evaluate_cands_error(std::vector<Cand>& candidates, Element* e, MeshFunction<Scalar>* rsln)
131  {
132  bool tri = e->is_triangle();
133 
134  // find range of orders
135  typename OptimumSelector<Scalar>::CandsInfo info_h, info_p, info_aniso;
136  this->update_cands_info(candidates, info_h, info_p, info_aniso);
137 
138  // calculate squared projection errors of elements of candidates
139  CandElemProjError herr[4], anisoerr[4], perr;
140  calc_projection_errors(e, info_h, info_p, info_aniso, rsln, herr, perr, anisoerr);
141 
142  //evaluate errors and dofs
143  for (unsigned i = 0; i < candidates.size(); i++)
144  {
145  Cand& c = candidates[i];
146  double error_squared = 0.0;
147  if (tri)
148  {
149  switch (c.split)
150  {
151  case H2D_REFINEMENT_H:
152  error_squared = 0.0;
153  for (int j = 0; j < H2D_MAX_ELEMENT_SONS; j++)
154  {
155  int order = H2D_GET_H_ORDER(c.p[j]);
156  c.errors[j] = herr[j][order][order];
157  error_squared += c.errors[j];
158  }
159  //element of a candidate occupies 1/4 of the reference domain defined over a candidate
160  error_squared *= 0.25;
161  break;
162 
163  case H2D_REFINEMENT_P:
164  {
165  int order = H2D_GET_H_ORDER(c.p[0]);
166  c.errors[0] = perr[order][order];
167  error_squared = perr[order][order];
168  }
169  break;
170 
171  default:
172  throw Hermes::Exceptions::Exception("Unknown split type \"%d\" at candidate %d", c.split, i);
173  }
174  }
175  else
176  {
177  //quad
178  switch (c.split)
179  {
180  case H2D_REFINEMENT_H:
181  error_squared = 0.0;
182  for (int j = 0; j < H2D_MAX_ELEMENT_SONS; j++)
183  {
184  int order_h = H2D_GET_H_ORDER(c.p[j]), order_v = H2D_GET_V_ORDER(c.p[j]);
185  c.errors[j] = herr[j][order_h][order_v];
186  error_squared += c.errors[j];
187  }
188  //element of a candidate occupies 1/4 of the reference domain defined over a candidate
189  error_squared *= 0.25;
190  break;
191 
194  {
195  error_squared = 0.0;
196  for (int j = 0; j < 2; j++)
197  {
198  c.errors[j] = anisoerr[(c.split == H2D_REFINEMENT_H_ANISO_H) ? j : j + 2][H2D_GET_H_ORDER(c.p[j])][H2D_GET_V_ORDER(c.p[j])];
199  error_squared += c.errors[j];
200  }
201  //element of a candidate occupies 1/2 of the reference domain defined over a candidate
202  error_squared *= 0.5;
203  }
204  break;
205 
206  case H2D_REFINEMENT_P:
207  {
208  int order_h = H2D_GET_H_ORDER(c.p[0]), order_v = H2D_GET_V_ORDER(c.p[0]);
209  c.errors[0] = perr[order_h][order_v];
210  error_squared = c.errors[0];
211  }
212  break;
213 
214  default:
215  throw Hermes::Exceptions::Exception("Unknown split type \"%d\" at candidate %d", c.split, i);
216  }
217  }
218 
219  //calculate error from squared error
220  c.error = sqrt(error_squared);
221 
222  //apply weights
223  switch (c.split)
224  {
225  case H2D_REFINEMENT_H: c.error *= error_weight_h; break;
228  case H2D_REFINEMENT_P: c.error *= error_weight_p; break;
229  default: throw Hermes::Exceptions::Exception("Unknown split type \"%d\" at candidate %d", c.split, i);
230  }
231  }
232  }
233 
234  template<typename Scalar>
236  {
237  ElementMode2D mode = e->get_mode();
238 
239  // select quadrature, obtain integration points and weights
240  Quad2D* quad = &g_quad_2d_std;
241  rsln->set_quad_2d(quad);
242  double3* gip_points = quad->get_points(H2DRS_INTR_GIP_ORDER, mode);
243  int num_gip_points = quad->get_num_points(H2DRS_INTR_GIP_ORDER, mode);
244 
245  // everything is done on the reference domain
246  Scalar* rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS];
247 
248  Solution<Scalar>* rsln_sln = dynamic_cast<Solution<Scalar>*>(rsln);
249  if (rsln_sln != nullptr)
250  rsln_sln->enable_transform(false);
251 
252  // obtain reference solution values on all four refined sons
253  Element* base_element = rsln->get_mesh()->get_element(e->id);
254 
255  // value on base element.
256  if (base_element->active)
257  {
258  for (int son = 0; son < H2D_MAX_ELEMENT_SONS; son++)
259  {
260  rsln->set_active_element(base_element);
261  rsln->push_transform(son);
263 
264  //obtain precalculated values
265  precalc_ref_solution(son, rsln, e, H2DRS_INTR_GIP_ORDER, rval);
266  }
267  }
268  else
269  {
270  for (int son = 0; son < H2D_MAX_ELEMENT_SONS; son++)
271  {
272  rsln->set_active_element(base_element->sons[son]);
274 
275  //obtain precalculated values
276  precalc_ref_solution(son, rsln, e, H2DRS_INTR_GIP_ORDER, rval);
277  }
278  }
279 
280  //retrieve transformations
281  Trf* trfs = nullptr;
282  int num_noni_trfs = 0;
283  if (mode == HERMES_MODE_TRIANGLE)
284  {
285  trfs = tri_trf;
286  num_noni_trfs = H2D_TRF_TRI_NUM;
287  }
288  else
289  {
290  trfs = quad_trf;
291  num_noni_trfs = H2D_TRF_QUAD_NUM;
292  }
293 
294  // precalculate values of shape functions
295  TrfShape empty_shape_vals;
296 
297  if (!cached_shape_vals_valid[mode])
298  {
299 #pragma omp critical (cached_shape_vals_valid)
300  if (!cached_shape_vals_valid[mode])
301  {
302  precalc_ortho_shapes(gip_points, num_gip_points, trfs, num_noni_trfs, this->shape_indices[mode], this->max_shape_inx[mode], cached_shape_ortho_vals[mode], mode);
303  precalc_shapes(gip_points, num_gip_points, trfs, num_noni_trfs, this->shape_indices[mode], this->max_shape_inx[mode], cached_shape_vals[mode], mode);
304  cached_shape_vals_valid[mode] = true;
305  }
306  }
307 
308  //issue a warning if ortho values are defined and the selected cand_list might benefit from that but it cannot because elements do not have uniform orders
309  if (!warn_uniform_orders && mode == HERMES_MODE_QUAD && !cached_shape_ortho_vals[mode][H2D_TRF_IDENTITY].empty())
310  {
311  warn_uniform_orders = true;
312  if (this->cand_list == H2D_H_ISO || this->cand_list == H2D_H_ANISO || this->cand_list == H2D_P_ISO || this->cand_list == H2D_HP_ISO || this->cand_list == H2D_HP_ANISO_H)
313  {
314  this->warn_if(!info_h.uniform_orders || !info_aniso.uniform_orders || !info_p.uniform_orders, "Possible inefficiency: %s might be more efficient if the input mesh contains elements with uniform orders strictly.", get_cand_list_str(this->cand_list));
315  }
316  }
317 
318  TrfShape& svals = cached_shape_vals[mode];
319  TrfShape& ortho_svals = cached_shape_ortho_vals[mode];
320 
321 #pragma region candidatesEvaluation
322  //H-candidates
323  if (!info_h.is_empty())
324  {
325  if (base_element->active)
326  {
327  Trf* sub_trfs[4] = { &trfs[0], &trfs[1], &trfs[2], &trfs[3] };
328  std::vector<TrfShapeExp>* p_trf_svals[4] = { &svals[0], &svals[1], &svals[2], &svals[3] };
329  std::vector<TrfShapeExp>* p_trf_ortho_svals[4] = { &ortho_svals[0], &ortho_svals[1], &ortho_svals[2], &ortho_svals[3] };
330  for (int son = 0; son < H2D_MAX_ELEMENT_SONS; son++)
331  {
332  int sub_rval[1] = { son };
333  calc_error_cand_element(mode, gip_points, num_gip_points
334  , 1, &base_element, &sub_trfs[son], sub_rval
335  , &p_trf_svals[son], &p_trf_ortho_svals[son]
336  , info_h, herr[son], rval);
337  }
338  }
339  else
340  {
341  Trf* p_trf_identity[1] = { &trfs[H2D_TRF_IDENTITY] };
342  std::vector<TrfShapeExp>* p_trf_svals[1] = { &svals[H2D_TRF_IDENTITY] };
343  std::vector<TrfShapeExp>* p_trf_ortho_svals[1] = { &ortho_svals[H2D_TRF_IDENTITY] };
344  for (int son = 0; son < H2D_MAX_ELEMENT_SONS; son++)
345  {
346  int sub_rval[1] = { son };
347  calc_error_cand_element(mode, gip_points, num_gip_points
348  , 1, &base_element->sons[son], p_trf_identity, sub_rval
349  , p_trf_svals, p_trf_ortho_svals
350  , info_h, herr[son], rval);
351  }
352  }
353  }
354 
355  //ANISO-candidates
356  if (!info_aniso.is_empty())
357  {
358  if (base_element->active)
359  {
360  const int tr[4][2] = { { 0, 1 }, { 3, 2 }, { 0, 3 }, { 1, 2 } };
361  for (int version = 0; version < 4; version++)
362  {
363  Trf* sub_trfs[2] = { &trfs[tr[version][0]], &trfs[tr[version][1]] };
364  Element* sub_domains[2] = { base_element, base_element };
365  int sub_rval[2] = { tr[version][0], tr[version][1] };
366  std::vector<TrfShapeExp>* sub_svals[2] = { &svals[tr[version][0]], &svals[tr[version][1]] };
367  std::vector<TrfShapeExp>* sub_ortho_svals[2] = { &ortho_svals[tr[version][0]], &ortho_svals[tr[version][1]] };
368  calc_error_cand_element(mode, gip_points, num_gip_points
369  , 2, sub_domains, sub_trfs, sub_rval
370  , sub_svals, sub_ortho_svals
371  , info_aniso, anisoerr[version], rval);
372  }
373  }
374  else
375  {
376  //indices of sons for sub-areas
377  const int sons[4][2] = { { 0, 1 }, { 3, 2 }, { 0, 3 }, { 1, 2 } };
378  //indices of ref. domain transformations for sub-areas
379  const int tr[4][2] = { { 6, 7 }, { 6, 7 }, { 4, 5 }, { 4, 5 } };
380  for (int version = 0; version < 4; version++)
381  { // 2 elements for vertical split, 2 elements for horizontal split
382  Trf* sub_trfs[2] = { &trfs[tr[version][0]], &trfs[tr[version][1]] };
383  Element* sub_domains[2] = { base_element->sons[sons[version][0]], base_element->sons[sons[version][1]] };
384  int sub_rval[2] = { sons[version][0], sons[version][1] };
385  std::vector<TrfShapeExp>* sub_svals[2] = { &svals[tr[version][0]], &svals[tr[version][1]] };
386  std::vector<TrfShapeExp>* sub_ortho_svals[2] = { &ortho_svals[tr[version][0]], &ortho_svals[tr[version][1]] };
387  calc_error_cand_element(mode, gip_points, num_gip_points
388  , 2, sub_domains, sub_trfs, sub_rval
389  , sub_svals, sub_ortho_svals
390  , info_aniso, anisoerr[version], rval);
391  }
392  }
393  }
394 
395  //P-candidates
396  if (!info_p.is_empty())
397  {
398  if (base_element->active)
399  {
400  Trf* sub_trfs[4] = { &trfs[0], &trfs[1], &trfs[2], &trfs[3] };
401  int sub_rval[4] = { 0, 1, 2, 3 };
402  std::vector<TrfShapeExp>* sub_svals[4] = { &svals[0], &svals[1], &svals[2], &svals[3] };
403  std::vector<TrfShapeExp>* sub_ortho_svals[4] = { &ortho_svals[0], &ortho_svals[1], &ortho_svals[2], &ortho_svals[3] };
404  Element* sub_domains[4] = { base_element, base_element, base_element, base_element };
405 
406  calc_error_cand_element(mode, gip_points, num_gip_points
407  , 4, sub_domains, sub_trfs, sub_rval
408  , sub_svals, sub_ortho_svals
409  , info_p, perr, rval);
410  }
411  else
412  {
413  Trf* sub_trfs[4] = { &trfs[0], &trfs[1], &trfs[2], &trfs[3] };
414  int sub_rval[4] = { 0, 1, 2, 3 };
415  std::vector<TrfShapeExp>* sub_svals[4] = { &svals[0], &svals[1], &svals[2], &svals[3] };
416  std::vector<TrfShapeExp>* sub_ortho_svals[4] = { &ortho_svals[0], &ortho_svals[1], &ortho_svals[2], &ortho_svals[3] };
417 
418  calc_error_cand_element(mode, gip_points, num_gip_points
419  , 4, base_element->sons, sub_trfs, sub_rval
420  , sub_svals, sub_ortho_svals
421  , info_p, perr, rval);
422  }
423  }
424 #pragma endregion
425 
426  for (int son = 0; son < H2D_MAX_ELEMENT_SONS; son++)
427  this->free_ref_solution_data(son, rval);
428  }
429 
430  template<typename Scalar>
432  , double3* gip_points, int num_gip_points
433  , const int num_sub, Element** sub_domains, Trf** sub_trfs, int* sons
434  , std::vector<TrfShapeExp>** sub_nonortho_svals, std::vector<TrfShapeExp>** sub_ortho_svals
435  , const typename OptimumSelector<Scalar>::CandsInfo& info
436  , CandElemProjError errors_squared, Scalar* rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS]
437  )
438  {
439  //allocate space
440  int max_num_shapes = this->next_order_shape[mode][this->max_order == H2DRS_DEFAULT_ORDER ? H2DRS_MAX_ORDER : this->max_order];
441  Scalar* right_side = new Scalar[max_num_shapes];
442  int* shape_inxs = new int[max_num_shapes];
443  //solver data
444  int* indx = new int[max_num_shapes];
445  //solver data
446  double* d = new double[max_num_shapes];
447  double** proj_matrix = new_matrix<double>(max_num_shapes, max_num_shapes);
448  ProjMatrixCache& proj_matrices = proj_matrix_cache[mode];
449  std::vector<typename OptimumSelector<Scalar>::ShapeInx>& full_shape_indices = this->shape_indices[mode];
450 
451  //check whether ortho-svals are available
452  bool ortho_svals_available = true;
453  for (int i = 0; i < num_sub && ortho_svals_available; i++)
454  ortho_svals_available &= !sub_ortho_svals[i]->empty();
455 
456  // An array of cached right-hand side values.
457  std::vector< ValueCacheItem<Scalar> > nonortho_rhs_cache;
458  std::vector< ValueCacheItem<Scalar> > ortho_rhs_cache;
459  for (int i = 0; i <= this->max_shape_inx[mode]; i++)
460  {
461  nonortho_rhs_cache.push_back(ValueCacheItem<Scalar>());
462  ortho_rhs_cache.push_back(ValueCacheItem<Scalar>());
463  }
464 
465  //calculate for all orders
466  double sub_area_corr_coef = 1.0 / num_sub;
467  OrderPermutator order_perm(info.min_quad_order, info.max_quad_order, mode == HERMES_MODE_TRIANGLE || info.uniform_orders);
468  do
469  {
470  int quad_order = order_perm.get_quad_order();
471  int order_h = H2D_GET_H_ORDER(quad_order), order_v = H2D_GET_V_ORDER(quad_order);
472 
473  //build a list of shape indices from the full list
474  int num_shapes = 0;
475  unsigned int inx_shape = 0;
476  while (inx_shape < full_shape_indices.size())
477  {
478  typename OptimumSelector<Scalar>::ShapeInx& shape = full_shape_indices[inx_shape];
479  if (order_h >= shape.order_h && order_v >= shape.order_v)
480  {
481  if (num_shapes >= max_num_shapes)
482  throw Exceptions::Exception("more shapes than predicted, possible incosistency");
483  shape_inxs[num_shapes] = shape.inx;
484  num_shapes++;
485  }
486  inx_shape++;
487  }
488 
489  //continue only if there are shapes to process
490  if (num_shapes == 0)
491  continue;
492 
493  bool use_ortho = ortho_svals_available && order_perm.get_order_h() == order_perm.get_order_v();
494 
495  //select a cache
496  std::vector< ValueCacheItem<Scalar> >& rhs_cache = use_ortho ? ortho_rhs_cache : nonortho_rhs_cache;
497  std::vector<TrfShapeExp>** sub_svals = use_ortho ? sub_ortho_svals : sub_nonortho_svals;
498 
499  //calculate projection matrix iff no ortho is used
500  if (!use_ortho)
501  {
502  if (!proj_matrices[order_h][order_v])
503  {
504 #pragma omp critical
505  {
506  if (!proj_matrices[order_h][order_v])
507  {
508  proj_matrices[order_h][order_v] = build_projection_matrix(gip_points, num_gip_points, shape_inxs, num_shapes, mode);
509  }
510  }
511  }
512 
513  //copy projection matrix because original matrix will be modified
514  copy_matrix(proj_matrix, proj_matrices[order_h][order_v], num_shapes, num_shapes);
515  }
516 
517  //build right side (fill cache values that are missing)
518  for (int inx_sub = 0; inx_sub < num_sub; inx_sub++)
519  {
520  Element* this_sub_domain = sub_domains[inx_sub];
521  ElemSubTrf this_sub_trf = { sub_trfs[inx_sub], 1 / sub_trfs[inx_sub]->m[0], 1 / sub_trfs[inx_sub]->m[1] };
522  ElemGIP this_sub_gip = { gip_points, num_gip_points };
523  std::vector<TrfShapeExp>& this_sub_svals = *(sub_svals[inx_sub]);
524 
525  for (int k = 0; k < num_shapes; k++)
526  {
527  int shape_inx = shape_inxs[k];
528  ValueCacheItem<Scalar>& shape_rhs_cache = rhs_cache[shape_inx];
529  if (!shape_rhs_cache.is_valid())
530  {
531  TrfShapeExp empty_sub_vals;
532  ElemSubShapeFunc this_sub_shape = { shape_inx, this_sub_svals.empty() ? empty_sub_vals : this_sub_svals[shape_inx] };
533  shape_rhs_cache.set(shape_rhs_cache.get() + evaluate_rhs_subdomain(this_sub_domain, this_sub_gip, sons[inx_sub], this_sub_trf, this_sub_shape, rval));
534  }
535  }
536  }
537 
538  //copy values from cache and apply area correction coefficient
539  for (int k = 0; k < num_shapes; k++)
540  {
541  ValueCacheItem<Scalar>& rhs_cache_value = rhs_cache[shape_inxs[k]];
542  right_side[k] = sub_area_corr_coef * rhs_cache_value.get();
543  rhs_cache_value.mark();
544  }
545 
546  //solve iff no ortho is used
547  if (!use_ortho)
548  {
549  ludcmp(proj_matrix, num_shapes, indx, d);
550  lubksb<double, Scalar>(proj_matrix, num_shapes, indx, right_side);
551  }
552 
553  //calculate error
554  double error_squared = 0;
555  for (int inx_sub = 0; inx_sub < num_sub; inx_sub++)
556  {
557  Element* this_sub_domain = sub_domains[inx_sub];
558  ElemSubTrf this_sub_trf = { sub_trfs[inx_sub], 1 / sub_trfs[inx_sub]->m[0], 1 / sub_trfs[inx_sub]->m[1] };
559  ElemGIP this_sub_gip = { gip_points, num_gip_points };
560  ElemProj elem_proj = { shape_inxs, num_shapes, *(sub_svals[inx_sub]), right_side, quad_order };
561 
562  error_squared += evaluate_error_squared_subdomain(this_sub_domain, this_sub_gip, sons[inx_sub], this_sub_trf, elem_proj, rval);
563  }
564  //apply area correction coefficient
565  errors_squared[order_h][order_v] = error_squared * sub_area_corr_coef;
566  } while (order_perm.next());
567 
568  free_with_check(proj_matrix, true);
569  delete[] right_side;
570  delete[] shape_inxs;
571  delete[] indx;
572  delete[] d;
573  }
574 
575  template class HERMES_API ProjBasedSelector < double > ;
576  template class HERMES_API ProjBasedSelector < std::complex<double> > ;
577  }
578  }
579 }
RefinementType split
A refinement, see the enum RefinementType.
Definition: candidates.h:90
unsigned short get_order_v() const
Returns the current vertical order.
Definition: adapt.h:24
int order_v
A minimal vertical order of an element that can use this shape function.
virtual void push_transform(int son)
Definition: function.cpp:115
HERMES_API Trf quad_trf[H2D_TRF_NUM]
A table of quad sub-subdomain transforms. Only first ::H2D_TRF_QUAD_NUM transformations are valid...
virtual void precalc_shapes(const double3 *gip_points, const int num_gip_points, const Trf *trfs, const int num_noni_trfs, const std::vector< typename OptimumSelector< Scalar >::ShapeInx > &shapes, const int max_shape_inx, TrfShape &svals, ElementMode2D mode)
Calculates values of shape function at GIP for all transformations.
std::vector< TrfShapeExp > TrfShape[H2D_TRF_NUM]
Evaluated shapes for all possible transformations for all points. The first index is a transformation...
int id
element id number
Definition: element.h:112
virtual void precalc_ref_solution(int inx_son, MeshFunction< Scalar > *rsln, Element *element, int intr_gip_order, Scalar *rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS])=0
Returns an array of values of the reference solution at integration points.
Stores one element of a mesh.
Definition: element.h:107
#define H2DRS_DEFAULT_ERR_WEIGHT_P
A default multiplicative coefficient of an error of a P-candidate.
Definition: global.h:74
virtual void set_active_element(Element *e)
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
double error_weight_h
A coefficient that multiplies error of H-candidate. The default value is H2DRS_DEFAULT_ERR_WEIGHT_H.
void calc_error_cand_element(const ElementMode2D mode, double3 *gip_points, int num_gip_points, const int num_sub, Element **sub_domains, Trf **sub_trfs, int *sons, std::vector< TrfShapeExp > **sub_nonortho_svals, std::vector< TrfShapeExp > **sub_ortho_svals, const typename OptimumSelector< Scalar >::CandsInfo &info, CandElemProjError errors_squared, Scalar *rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS])
Calculate projection errors of an element of an candidate considering multiple orders.
virtual void set_quad_2d(Quad2D *quad_2d)
HERMES_API const char * get_cand_list_str(const CandList cand_list)
Returns a string representation of a predefined candidate list.
Definition: candidates.cpp:9
bool is_valid() const
Returns true if value is mared as valid.
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
const int max_order
A maximum allowed order.
Definition: selector.h:77
Integration points in the reference domain of an element of a candidate.
virtual double ** build_projection_matrix(double3 *gip_points, int num_gip_points, const int *shape_inx, const int num_shapes, ElementMode2D)=0
Builds projection matrix using a given set of shapes.
void set_quad_order(unsigned short order, unsigned short mask=H2D_FN_DEFAULT)
Definition: function.cpp:59
double error_weight_aniso
A coefficient that multiplies error of ANISO-candidate. The default value is H2DRS_DEFAULT_ERR_WEIGHT...
bool warn_uniform_orders
True if the selector has already warned about possible inefficiency.
double ** ProjMatrixCache[H2DRS_MAX_ORDER+2][H2DRS_MAX_ORDER+2]
A projection matrix cache type.
virtual void free_ref_solution_data(int inx_son, Scalar *rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS])=0
Frees the data allocated in precalc_ref_solution.
std::vector< ShapeInx > shape_indices[H2D_NUM_MODES]
Shape indices. The first index is a mode (ElementMode2D).
unsigned short get_quad_order() const
Returns the current order in an encoded form.
ProjMatrixCache proj_matrix_cache[H2D_NUM_MODES]
An array of projection matrices.
int max_shape_inx[H2D_NUM_MODES]
A maximum index of a shape function. The first index is a mode (ElementMode2D).
T get() const
Returns the value. It does check the state of the value.
A general projection-based selector.
Definition: function.h:33
#define H2D_NUM_MODES
Internal.
Definition: global.h:35
TrfShape * cached_shape_ortho_vals
Precalculated valus of orthogonalized shape functions.
H- and ANISO-candidates only. Hermes::Orders are not modified.
Definition: candidates.h:51
unsigned short get_order_h() const
Returns the current horizontal order.
CandList
Predefined list of candidates.
Definition: candidates.h:46
double CandElemProjError[H2DRS_MAX_ORDER+2][H2DRS_MAX_ORDER+2]
Error of an element of a candidate for various permutations of orders.
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.
bool * cached_shape_vals_valid
True if shape values were already initialized.
H- and P-candidates only. Hermes::Orders are modified uniformly.
Definition: candidates.h:52
#define H2D_GET_H_ORDER(encoded_order)
Macros for combining quad horizontal and vertical encoded_orders.
Definition: global.h:98
bool next()
Moves to the next permutation of orders.
virtual Scalar evaluate_rhs_subdomain(Element *sub_elem, const ElemGIP &sub_gip, int son, const ElemSubTrf &sub_trf, const ElemSubShapeFunc &sub_shape, Scalar *rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS])=0
Evaluates a value of the right-hande side in a subdomain.
2D transformation.
Definition: transformable.h:29
HERMES_API Trf tri_trf[H2D_TRF_NUM]
A table of triangle sub-subdomain transforms. Only first ::H2D_TRF_TRI_NUM transformations are valid...
Should be exactly the same as is the count of enum ShapesetType.
Definition: shapeset.h:95
#define H2DRS_INTR_GIP_ORDER
An integration order used to integrate while evaluating a candidate.
Definition: global.h:68
TrfShape * cached_shape_vals
Precalculate values of shape functions.
#define H2DRS_DEFAULT_ERR_WEIGHT_H
A default multiplicative coefficient of an error of a H-candidate.
Definition: global.h:73
Hermes::Order permutator. Generates all permutations of orders from a set defined by a range of order...
double error_weight_p
A coefficient that multiplies error of P-candidate. The default value is H2DRS_DEFAULT_ERR_WEIGHT_P.
void enable_transform(bool enable=true)
Definition: solution.cpp:620
virtual double evaluate_error_squared_subdomain(Element *sub_elem, const ElemGIP &sub_gip, int son, const ElemSubTrf &sub_trf, const ElemProj &elem_proj, Scalar *rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS])=0
Evaluates an squared error of a projection of an element of a candidate onto subdomains.
A transformation from a reference domain of a subdomain to a reference domain of an element of a cand...
TrfShapeExp()
A default contructor. Creates an empty instance.
int next_order_shape[H2D_NUM_MODES][H2DRS_MAX_ORDER+1]
An index to the array OptimumSelector::shape_indices of a shape function of the next uniform order...
A selector that chooses an optimal candidates based on a score.
Definition: function.h:32
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...
ANISO-refienement. The element is split along the horizontal axis. Quadrilaterals only...
#define H2D_MAX_ELEMENT_SONS
Macros.
Definition: global.h:30
int order_h
A minimal horizonal order of an element that can use this shape function.
virtual void evaluate_cands_error(std::vector< Cand > &candidates, Element *e, MeshFunction< Scalar > *rsln)
Calculates error of candidates.
void set_error_weights(double weight_h=H2DRS_DEFAULT_ERR_WEIGHT_H, double weight_p=H2DRS_DEFAULT_ERR_WEIGHT_P, double weight_aniso=H2DRS_DEFAULT_ERR_WEIGHT_ANISO)
Sets error weights.
Element * sons[H2D_MAX_ELEMENT_SONS]
son elements (up to four)
Definition: element.h:140
#define H2DRS_DEFAULT_ERR_WEIGHT_ANISO
A default multiplicative coefficient of an error of a ANISO-candidate.
Definition: global.h:75
virtual void precalc_ortho_shapes(const double3 *gip_points, const int num_gip_points, const Trf *trfs, const int num_noni_trfs, const std::vector< typename OptimumSelector< Scalar >::ShapeInx > &shapes, const int max_shape_inx, TrfShape &ortho_svals, ElementMode2D mode)
Calculates values of orthogonalized shape function at GIP for all transformations.
MeshSharedPtr get_mesh() const
Return the mesh.
H-, ANISO- and P-candidates. Hermes::Orders are modified uniformly.
Definition: candidates.h:53
virtual void calc_projection_errors(Element *e, const typename OptimumSelector< Scalar >::CandsInfo &info_h, const typename OptimumSelector< Scalar >::CandsInfo &info_p, const typename OptimumSelector< Scalar >::CandsInfo &info_aniso, MeshFunction< Scalar > *rsln, CandElemProjError herr[H2D_MAX_ELEMENT_SONS], CandElemProjError perr, CandElemProjError anisoerr[H2D_MAX_ELEMENT_SONS])
Calculates projection errors of an elements of candidates for all permutations of orders...
bool active
0 = active, no sons; 1 = inactive (refined), has sons
Definition: element.h:114
void mark(int new_state=H2DRS_VALCACHE_VALID)
Marks a value.
P-candidates only. Hermes::Orders are modified uniformly.
Definition: candidates.h:48
Represents the solution of a PDE.
Definition: api2d.h:35
#define H2DRS_DEFAULT_ORDER
A default order. Used to indicate an unkonwn order or a maximum support order.
Definition: global.h:103
bool is_empty() const
Returns true if there are no candidates.