Hermes2D  3.0
thread_linearizer.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 "thread_linearizer.h"
17 #include "refmap.h"
18 #include "traverse.h"
19 #include "exact_solution.h"
20 #include "api2d.h"
21 
22 //#define DEBUG_LINEARIZER
23 
24 static const int default_allocation_multiplier_vertices = 10;
25 static const int default_allocation_multiplier_triangles = 15;
26 static const int default_allocation_multiplier_edges = 10;
27 
28 static const int default_allocation_minsize_vertices = 10000;
29 static const int default_allocation_minsize_triangles = 15000;
30 static const int default_allocation_minsize_edges = 10000;
31 
32 static const double vertex_relative_tolerance = 0.01;
33 
34 namespace Hermes
35 {
36  namespace Hermes2D
37  {
38  namespace Views
39  {
40  template<typename LinearizerDataDimensions>
41  ThreadLinearizerMultidimensional<LinearizerDataDimensions>::ThreadLinearizerMultidimensional(LinearizerMultidimensional<LinearizerDataDimensions>* linearizer) : criterion(linearizer->criterion)
42  {
43  vertex_size = 0;
44  triangle_size = 0;
45  edges_size = 0;
46 
47  // OpenGL part.
48  triangles = nullptr;
49  edges = nullptr;
50  edge_markers = nullptr;
51 
52  // FileExport part.
53  triangle_indices = nullptr;
54 
55  // Common part.
56  vertices = nullptr;
57  triangle_markers = nullptr;
58  hash_table = nullptr;
59  info = nullptr;
60  }
61 
62  template<typename LinearizerDataDimensions>
63  void ThreadLinearizerMultidimensional<LinearizerDataDimensions>::init_linearizer_data(LinearizerMultidimensional<LinearizerDataDimensions>* linearizer)
64  {
65  this->criterion = linearizer->criterion;
66  this->curvature_epsilon = linearizer->curvature_epsilon;
67  this->user_xdisp = linearizer->user_xdisp;
68  this->user_ydisp = linearizer->user_ydisp;
69  this->dmult = linearizer->dmult;
70  this->linearizerOutputType = linearizer->linearizerOutputType;
71 
72  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
73  {
74  this->component[k] = linearizer->component[k];
75  this->item[k] = linearizer->item[k];
76  this->value_type[k] = linearizer->value_type[k];
77  }
78  }
79 
80  template<typename LinearizerDataDimensions>
81  void ThreadLinearizerMultidimensional<LinearizerDataDimensions>::free()
82  {
83  // OpenGL part.
84  free_with_check(this->triangles, true);
85  free_with_check(this->edges, true);
86  free_with_check(this->edge_markers, true);
87  // FileExport part.
88  free_with_check(this->triangle_indices, true);
89  // Common part.
90  free_with_check(this->vertices, true);
91  free_with_check(this->triangle_markers, true);
92  free_with_check(this->hash_table, true);
93  free_with_check(this->info, true);
94  }
95 
96  template<typename LinearizerDataDimensions>
97  ThreadLinearizerMultidimensional<LinearizerDataDimensions>::~ThreadLinearizerMultidimensional()
98  {
99  free();
100  }
101 
102  template<typename LinearizerDataDimensions>
103  void ThreadLinearizerMultidimensional<LinearizerDataDimensions>::init_processing(MeshFunctionSharedPtr<double>* sln, LinearizerMultidimensional<LinearizerDataDimensions>* linearizer)
104  {
105  this->init_linearizer_data(linearizer);
106 
107  // Functions.
108  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
109  {
110  Solution<double>* solution = dynamic_cast<Solution<double>*>(sln[k].get());
111  if (solution && solution->get_type() == HERMES_SLN)
112  {
113  fns[k] = new Solution<double>();
114  fns[k]->copy(solution);
115  }
116  else
117  fns[k] = sln[k]->clone();
118  fns[k]->set_quad_2d(&g_quad_lin);
119  }
120  if (user_xdisp)
121  {
122  Solution<double>* xdisp_solution = dynamic_cast<Solution<double>*>(linearizer->xdisp.get());
123  if (xdisp_solution && xdisp_solution->get_type() == HERMES_SLN)
124  {
125  fns[LinearizerDataDimensions::dimension] = new Solution<double>();
126  fns[LinearizerDataDimensions::dimension]->copy(linearizer->xdisp);
127  }
128  else
129  fns[LinearizerDataDimensions::dimension] = linearizer->xdisp->clone();
130 
131  fns[LinearizerDataDimensions::dimension]->set_quad_2d(&g_quad_lin);
132  }
133  if (user_ydisp)
134  {
135  Solution<double>* ydisp_solution = dynamic_cast<Solution<double>*>(linearizer->ydisp.get());
136  if (ydisp_solution && ydisp_solution->get_type() == HERMES_SLN)
137  {
138  fns[LinearizerDataDimensions::dimension + (user_xdisp ? 1 : 0)] = new Solution<double>();
139  fns[LinearizerDataDimensions::dimension + (user_xdisp ? 1 : 0)]->copy(linearizer->ydisp);
140  }
141  else
142  fns[LinearizerDataDimensions::dimension + (user_xdisp ? 1 : 0)] = linearizer->ydisp->clone();
143 
144  fns[LinearizerDataDimensions::dimension + (user_xdisp ? 1 : 0)]->set_quad_2d(&g_quad_lin);
145  }
146 
147  // Init storage data & counts.
148  this->reallocate(sln[0]->get_mesh());
149  }
150 
151  template<typename LinearizerDataDimensions>
152  void ThreadLinearizerMultidimensional<LinearizerDataDimensions>::reallocate(MeshSharedPtr mesh)
153  {
154  int number_of_elements = mesh->get_num_elements();
155 
156  this->vertex_size = std::max(default_allocation_multiplier_vertices * number_of_elements, std::max(this->vertex_size, default_allocation_minsize_vertices));
157  this->triangle_size = std::max(default_allocation_multiplier_triangles * number_of_elements, std::max(this->triangle_size, default_allocation_minsize_triangles));
158  this->edges_size = std::max(default_allocation_multiplier_edges * number_of_elements, std::max(this->edges_size, default_allocation_minsize_edges));
159 
160  // Set counts.
161  this->vertex_count = 0;
162  this->triangle_count = 0;
163  this->edges_count = 0;
164 
165  if (this->linearizerOutputType == OpenGL)
166  {
167  this->triangles = realloc_with_check<ThreadLinearizerMultidimensional, typename LinearizerDataDimensions::triangle_t>(this->triangles, this->triangle_size, this);
168  this->edges = realloc_with_check<ThreadLinearizerMultidimensional, typename LinearizerDataDimensions::edge_t>(this->edges, this->edges_size, this);
169  this->edge_markers = realloc_with_check<ThreadLinearizerMultidimensional, int>(this->edge_markers, this->edges_size, this);
170  }
171  else
172  this->triangle_indices = realloc_with_check<ThreadLinearizerMultidimensional, triangle_indices_t>(this->triangle_indices, this->triangle_size, this);
173 
174  this->vertices = realloc_with_check<ThreadLinearizerMultidimensional, typename LinearizerDataDimensions::vertex_t>(this->vertices, this->vertex_size, this);
175  this->triangle_markers = realloc_with_check<ThreadLinearizerMultidimensional, int>(this->triangle_markers, this->triangle_size, this);
176 
177  this->hash_table = malloc_with_check<ThreadLinearizerMultidimensional<LinearizerDataDimensions>, int>(this->vertex_size, this, true);
178  memset(this->hash_table, 0xff, sizeof(int)* this->vertex_size);
179 
180  this->info = malloc_with_check<ThreadLinearizerMultidimensional<LinearizerDataDimensions>, internal_vertex_info_t>(this->vertex_size, this, true);
181  }
182 
183  template<typename LinearizerDataDimensions>
184  void ThreadLinearizerMultidimensional<LinearizerDataDimensions>::deinit_processing()
185  {
186  for (unsigned int j = 0; j < (LinearizerDataDimensions::dimension + (this->user_xdisp ? 1 : 0) + (this->user_ydisp ? 1 : 0)); j++)
187  delete fns[j];
188 
189  free_with_check(this->hash_table, true);
190  free_with_check(this->info, true);
191  }
192 
193  template<typename LinearizerDataDimensions>
194  void ThreadLinearizerMultidimensional<LinearizerDataDimensions>::process_state(Traverse::State* current_state)
195  {
196  this->rep_element = current_state->e[0];
197  this->curved = this->rep_element->is_curved();
198  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
199  {
200  fns[k]->set_active_element(current_state->e[k]);
201  fns[k]->set_transform(current_state->sub_idx[k]);
202  fns[k]->set_quad_order(0, this->item[k]);
203  val[k] = fns[k]->get_values(component[k], value_type[k]);
204  }
205 
206  const double *dx = nullptr;
207  const double *dy = nullptr;
208 
209  if (user_xdisp)
210  {
211  fns[LinearizerDataDimensions::dimension]->set_active_element(current_state->e[LinearizerDataDimensions::dimension]);
212  fns[LinearizerDataDimensions::dimension]->set_transform(current_state->sub_idx[LinearizerDataDimensions::dimension]);
213  fns[LinearizerDataDimensions::dimension]->set_quad_order(0, H2D_FN_VAL);
214  dx = fns[LinearizerDataDimensions::dimension]->get_fn_values();
215  }
216 
217  if (user_ydisp)
218  {
219  fns[LinearizerDataDimensions::dimension + (user_xdisp ? 1 : 0)]->set_active_element(current_state->e[LinearizerDataDimensions::dimension + (user_xdisp ? 1 : 0)]);
220  fns[LinearizerDataDimensions::dimension + (user_xdisp ? 1 : 0)]->set_transform(current_state->sub_idx[LinearizerDataDimensions::dimension + (user_xdisp ? 1 : 0)]);
221  fns[LinearizerDataDimensions::dimension + (user_xdisp ? 1 : 0)]->set_quad_order(0, H2D_FN_VAL);
222  dy = fns[LinearizerDataDimensions::dimension + (user_xdisp ? 1 : 0)]->get_fn_values();
223  }
224 
225  int iv[H2D_MAX_NUMBER_VERTICES];
226  for (unsigned int i = 0; i < this->rep_element->get_nvert(); i++)
227  {
228  double x_disp = fns[0]->get_refmap()->get_phys_x(0)[i];
229  double y_disp = fns[0]->get_refmap()->get_phys_y(0)[i];
230  if (user_xdisp)
231  x_disp += dmult * dx[i];
232  if (user_ydisp)
233  y_disp += dmult * dy[i];
234 
235  double value[LinearizerDataDimensions::dimension];
236  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
237  value[k] = val[k][i];
238  iv[i] = this->get_vertex(-this->rep_element->vn[i]->id, -this->rep_element->vn[i]->id, x_disp, y_disp, value);
239  }
240 
241  // recur to sub-elements
242  if (current_state->e[0]->is_triangle())
243  process_triangle(iv[0], iv[1], iv[2], 0);
244  else
245  process_quad(iv[0], iv[1], iv[2], iv[3], 0);
246 
247 #ifndef DEBUG_LINEARIZER
248  if (this->linearizerOutputType == OpenGL)
249  {
250  for (unsigned int i = 0; i < this->rep_element->get_nvert(); i++)
251  process_edge(iv[i], iv[this->rep_element->next_vert(i)], this->rep_element->en[i]->marker);
252  }
253 #endif
254  }
255 
256  template<typename LinearizerDataDimensions>
257  void ThreadLinearizerMultidimensional<LinearizerDataDimensions>::process_triangle(int iv0, int iv1, int iv2, int level)
258  {
259  const double* values[LinearizerDataDimensions::dimension];
260  double* physical_x;
261  double* physical_y;
262  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
263  {
264  fns[k]->set_quad_order(1, item[k]);
265  values[k] = fns[k]->get_values(component[k], value_type[k]);
266  }
267  unsigned short* vertex_indices = tri_indices[0];
268 
269  if (curved)
270  {
271  // obtain physical element coordinates
272  RefMap* refmap = fns[0]->get_refmap();
273  physical_x = refmap->get_phys_x(1);
274  physical_y = refmap->get_phys_y(1);
275 
276  if (user_xdisp)
277  {
278  fns[LinearizerDataDimensions::dimension]->set_quad_order(1, H2D_FN_VAL);
279  const double* dx = fns[LinearizerDataDimensions::dimension]->get_fn_values();
280  for (int i = 0; i < lin_np_tri[1]; i++)
281  physical_x[i] += dmult*dx[i];
282  }
283  if (user_ydisp)
284  {
285  fns[LinearizerDataDimensions::dimension + (this->user_xdisp ? 1 : 0)]->set_quad_order(1, H2D_FN_VAL);
286  const double* dy = fns[LinearizerDataDimensions::dimension + (this->user_xdisp ? 1 : 0)]->get_fn_values();
287  for (int i = 0; i < lin_np_tri[1]; i++)
288  physical_y[i] += dmult*dy[i];
289  }
290  }
291 
292  // obtain linearized values and coordinates at the midpoints
293  for (int i = 0; i < 2 + LinearizerDataDimensions::dimension; i++)
294  {
295  midval[i][0] = (this->vertices[iv0][i] + this->vertices[iv1][i])*0.5;
296  midval[i][1] = (this->vertices[iv1][i] + this->vertices[iv2][i])*0.5;
297  midval[i][2] = (this->vertices[iv2][i] + this->vertices[iv0][i])*0.5;
298  };
299 
300  // determine whether or not to split the element
301  int split;
302  if (level == MAX_LINEARIZER_DIVISION_LEVEL)
303  split = 0;
304  else
305  {
306  if (this->criterion.adaptive)
307  this->split_decision(split, iv0, iv1, iv2, 0, rep_element->get_mode(), values, physical_x, physical_y, vertex_indices);
308  else
309  split = (level < this->criterion.refinement_level);
310  }
311 
312  // split the triangle if the error is too large, otherwise produce a linear triangle
313  if (split)
314  {
315  if (curved)
316  {
317  for (int i = 0; i < 3; i++)
318  {
319  midval[0][i] = physical_x[vertex_indices[i]];
320  midval[1][i] = physical_y[vertex_indices[i]];
321  }
322  }
323 
324  double values_vertices[5][LinearizerDataDimensions::dimension];
325  for (int v = 0; v < 3; v++)
326  {
327  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
328  values_vertices[v][k] = values[k][vertex_indices[v]];
329  }
330 
331  // obtain mid-edge vertices
332  int mid0 = get_vertex(iv0, iv1, midval[0][0], midval[1][0], values_vertices[0]);
333  int mid1 = get_vertex(iv1, iv2, midval[0][1], midval[1][1], values_vertices[1]);
334  int mid2 = get_vertex(iv2, iv0, midval[0][2], midval[1][2], values_vertices[2]);
335 
336  // recur to sub-elements
337  this->push_transforms(0);
338  process_triangle(iv0, mid0, mid2, level + 1);
339  this->pop_transforms();
340 
341  this->push_transforms(1);
342  process_triangle(mid0, iv1, mid1, level + 1);
343  this->pop_transforms();
344 
345  this->push_transforms(2);
346  process_triangle(mid2, mid1, iv2, level + 1);
347  this->pop_transforms();
348 
349  this->push_transforms(3);
350  process_triangle(mid1, mid2, mid0, level + 1);
351  this->pop_transforms();
352  }
353  else
354  add_triangle(iv0, iv1, iv2, this->rep_element->marker);
355 
356 #ifdef DEBUG_LINEARIZER
357  if (this->linearizerOutputType == OpenGL)
358  {
359  process_edge(iv0, iv1, this->rep_element->en[0]->marker);
360  process_edge(iv1, iv2, this->rep_element->en[1]->marker);
361  process_edge(iv2, iv0, this->rep_element->en[2]->marker);
362  }
363 #endif
364  }
365 
366  template<typename LinearizerDataDimensions>
367  void ThreadLinearizerMultidimensional<LinearizerDataDimensions>::process_quad(int iv0, int iv1, int iv2, int iv3, int level)
368  {
369  const double* values[LinearizerDataDimensions::dimension];
370  double* physical_x;
371  double* physical_y;
372  unsigned short* vertex_indices = quad_indices[0];
373  bool flip = this->quad_flip(iv0, iv1, iv2, iv3);
374  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
375  {
376  fns[k]->set_quad_order(1, item[k]);
377  values[k] = fns[k]->get_values(component[k], value_type[k]);
378  }
379 
380  if (curved)
381  {
382  // obtain physical element coordinates
383  RefMap* refmap = fns[0]->get_refmap();
384  physical_x = refmap->get_phys_x(1);
385  physical_y = refmap->get_phys_y(1);
386 
387  if (user_xdisp)
388  {
389  fns[LinearizerDataDimensions::dimension]->set_quad_order(1, H2D_FN_VAL);
390  const double* dx = fns[LinearizerDataDimensions::dimension]->get_fn_values();
391  for (int i = 0; i < lin_np_tri[1]; i++)
392  physical_x[i] += dmult*dx[i];
393  }
394  if (user_ydisp)
395  {
396  fns[LinearizerDataDimensions::dimension + (this->user_xdisp ? 1 : 0)]->set_quad_order(1, H2D_FN_VAL);
397  const double* dy = fns[LinearizerDataDimensions::dimension + (this->user_xdisp ? 1 : 0)]->get_fn_values();
398  for (int i = 0; i < lin_np_tri[1]; i++)
399  physical_y[i] += dmult*dy[i];
400  }
401  }
402 
403  // obtain linearized values and coordinates at the midpoints
404  for (int i = 0; i < 2 + LinearizerDataDimensions::dimension; i++)
405  {
406  midval[i][0] = (this->vertices[iv0][i] + this->vertices[iv1][i]) * 0.5;
407  midval[i][1] = (this->vertices[iv1][i] + this->vertices[iv2][i]) * 0.5;
408  midval[i][2] = (this->vertices[iv2][i] + this->vertices[iv3][i]) * 0.5;
409  midval[i][3] = (this->vertices[iv3][i] + this->vertices[iv0][i]) * 0.5;
410  midval[i][4] = (midval[i][0] + midval[i][2]) * 0.5;
411  };
412 
413  // the value of the middle point is not the average of the four vertex values, since quad == 2 triangles
414  midval[LinearizerDataDimensions::dimension + 1][4] = flip ?
415  (this->vertices[iv0][LinearizerDataDimensions::dimension + 1] + this->vertices[iv2][LinearizerDataDimensions::dimension + 1]) * 0.5
416  :
417  (this->vertices[iv1][LinearizerDataDimensions::dimension + 1] + this->vertices[iv3][LinearizerDataDimensions::dimension + 1]) * 0.5;
418 
419  // determine whether or not to split the element
420  int split;
421  if (level == MAX_LINEARIZER_DIVISION_LEVEL)
422  split = 0;
423  else
424  {
425  if (this->criterion.adaptive)
426  this->split_decision(split, iv0, iv1, iv2, 0, rep_element->get_mode(), values, physical_x, physical_y, vertex_indices);
427  else
428  split = (level < this->criterion.refinement_level ? 3 : 0);
429  }
430 
431  // split the quad if the error is too large, otherwise produce two linear triangles
432  if (split)
433  {
434  if (curved)
435  {
436  for (int i = 0; i < 5; i++)
437  {
438  midval[0][i] = physical_x[vertex_indices[i]];
439  midval[1][i] = physical_y[vertex_indices[i]];
440  }
441  }
442 
443  // obtain mid-edge and mid-element vertices
444  int mid0, mid1, mid2, mid3, mid4;
445  double values_vertices[5][LinearizerDataDimensions::dimension];
446  for (int v = 0; v < 5; v++)
447  {
448  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
449  values_vertices[v][k] = values[k][vertex_indices[v]];
450  }
451  if (split != 1) mid0 = get_vertex(iv0, iv1, midval[0][0], midval[1][0], values_vertices[0]);
452  if (split != 2) mid1 = get_vertex(iv1, iv2, midval[0][1], midval[1][1], values_vertices[1]);
453  if (split != 1) mid2 = get_vertex(iv2, iv3, midval[0][2], midval[1][2], values_vertices[2]);
454  if (split != 2) mid3 = get_vertex(iv3, iv0, midval[0][3], midval[1][3], values_vertices[3]);
455  if (split == 3) mid4 = get_vertex(mid0, mid2, midval[0][4], midval[1][4], values_vertices[4]);
456 
457  // recur to sub-elements
458  if (split == 3)
459  {
460  this->push_transforms(0);
461  process_quad(iv0, mid0, mid4, mid3, level + 1);
462  this->pop_transforms();
463 
464  this->push_transforms(1);
465  process_quad(mid0, iv1, mid1, mid4, level + 1);
466  this->pop_transforms();
467 
468  this->push_transforms(2);
469  process_quad(mid4, mid1, iv2, mid2, level + 1);
470  this->pop_transforms();
471 
472  this->push_transforms(3);
473  process_quad(mid3, mid4, mid2, iv3, level + 1);
474  this->pop_transforms();
475  }
476  else
477  if (split == 1) // h-split
478  {
479  this->push_transforms(4);
480  process_quad(iv0, iv1, mid1, mid3, level + 1);
481  this->pop_transforms();
482 
483  this->push_transforms(5);
484  process_quad(mid3, mid1, iv2, iv3, level + 1);
485  this->pop_transforms();
486  }
487  else // v-split
488  {
489  this->push_transforms(6);
490  process_quad(iv0, mid0, mid2, iv3, level + 1);
491  this->pop_transforms();
492 
493  this->push_transforms(7);
494  process_quad(mid0, iv1, iv2, mid2, level + 1);
495  this->pop_transforms();
496  }
497  }
498  else
499  {
500  // output two linear triangles,
501  if (!flip)
502  {
503  add_triangle(iv3, iv0, iv1, rep_element->marker);
504  add_triangle(iv1, iv2, iv3, rep_element->marker);
505  }
506  else
507  {
508  add_triangle(iv0, iv1, iv2, rep_element->marker);
509  add_triangle(iv2, iv3, iv0, rep_element->marker);
510  }
511 #ifdef DEBUG_LINEARIZER
512  if (this->linearizerOutputType == OpenGL)
513  {
514  process_edge(iv0, iv1, this->rep_element->en[0]->marker);
515  process_edge(iv1, iv2, this->rep_element->en[1]->marker);
516  process_edge(iv2, iv3, this->rep_element->en[2]->marker);
517  process_edge(iv3, iv0, this->rep_element->en[3]->marker);
518  }
519 #endif
520  }
521  }
522 
523  template<typename LinearizerDataDimensions>
524  void ThreadLinearizerMultidimensional<LinearizerDataDimensions>::split_decision(int& split, int iv0, int iv1, int iv2, int iv3, ElementMode2D mode, const double** values, double* physical_x, double* physical_y, unsigned short* vertex_indices) const
525  {
526  // Initialization.
527  split = 0;
528  bool done = false;
529 
530  // Core of the decision - calculate the approximate error of linearizing the normalized solution
531  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
532  {
533  // Errors (split element - non-split element values) in edge midpoints summed up.
534  double error = fabs(values[k][vertex_indices[0]] - midval[2 + k][0])
535  + fabs(values[k][vertex_indices[1]] - midval[2 + k][1])
536  + fabs(values[k][vertex_indices[2]] - midval[2 + k][2]);
537 
538  // For the quad we have one more midpoint.
539  if (mode == HERMES_MODE_QUAD)
540  error += fabs(values[k][vertex_indices[3]] - midval[2 + k][3]);
541 
542  // Divide by the edge count.
543  error /= (3 + mode);
544 
545  // Relative error.
546  // max_value_approx here is only an approximation - only taking into account the elements being processed by this thread.
547  double relative_error = error / this->max_value_approx;
548 
549  // Split ?
550  // See the header of this method (split_decision) for explanation.
551  // We put 3 here so that it is easier to test 'full split' both for quads && triangles.
552  split = (relative_error > this->criterion.error_tolerance) ? 3 : 0;
553 
554  // Quads - division type
555  if (mode == HERMES_MODE_QUAD && split)
556  {
557  double horizontal_error = fabs(values[k][vertex_indices[1]] - midval[2 + k][1]) + fabs(values[k][vertex_indices[3]] - midval[2 + k][3]);
558  double vertical_error = fabs(values[k][vertex_indices[0]] - midval[2 + k][0]) + fabs(values[k][vertex_indices[2]] - midval[2 + k][2]);
559 
560  // Decide whether to split horizontally or vertically only
561  // If one error is LINEARIZER_DIRECTIONAL_QUAD_REFINEMENT_REQUIREMENT larger than the other.
562  if (horizontal_error > LINEARIZER_DIRECTIONAL_QUAD_REFINEMENT_REQUIREMENT * vertical_error)
563  // h-split
564  split = 1;
565  else if (vertical_error > LINEARIZER_DIRECTIONAL_QUAD_REFINEMENT_REQUIREMENT * horizontal_error)
566  // v-split
567  split = 2;
568  else
569  split = 3;
570  }
571  }
572 
573  // If we are not splitting into four elements alreasdy and we have a curved element, check if we have to split because of the curvature.
574  if (curved && split != 3)
575  {
576  for (int i = 0; i < 3 + mode; i++)
577  {
578  double error = sqr(physical_x[vertex_indices[i]] - midval[0][i])
579  + sqr(physical_y[vertex_indices[i]] - midval[1][i]);
580 
581  double diameter = sqr(fns[0]->get_active_element()->diameter);
582 
583  split = (error / diameter) > this->curvature_epsilon ? 3 : split;
584  }
585  }
586  }
587 
588  template<typename LinearizerDataDimensions>
589  double ThreadLinearizerMultidimensional<LinearizerDataDimensions>::get_max_value(Traverse::State* current_state)
590  {
591  double local_max = 0.;
592  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
593  {
594  fns[k]->set_active_element(current_state->e[k]);
595  fns[k]->set_transform(current_state->sub_idx[k]);
596  fns[k]->set_quad_order(0, this->item[k]);
597  const double* val = fns[k]->get_values(component[k], value_type[k]);
598 
599  for (unsigned int i = 0; i < current_state->e[k]->get_nvert(); i++)
600  {
601  double f = fabs(val[i]);
602  if (f > local_max)
603  local_max = f;
604  }
605  }
606  return local_max;
607  }
608 
609  template<typename LinearizerDataDimensions>
610  bool ThreadLinearizerMultidimensional<LinearizerDataDimensions>::quad_flip(int iv0, int iv1, int iv2, int iv3) const
611  {
612  int a = (this->vertices[iv0][2] > this->vertices[iv1][2]) ? iv0 : iv1;
613  int b = (this->vertices[iv2][2] > this->vertices[iv3][2]) ? iv2 : iv3;
614  a = (this->vertices[a][2] > this->vertices[b][2]) ? a : b;
615  return (a == iv1 || a == iv3);
616  }
617 
618  template<typename LinearizerDataDimensions>
619  void ThreadLinearizerMultidimensional<LinearizerDataDimensions>::push_transforms(int transform)
620  {
621  fns[0]->push_transform(transform);
622 
623  if (user_xdisp)
624  {
625  if (fns[1] != fns[0])
626  fns[1]->push_transform(transform);
627  }
628  if (user_ydisp)
629  {
630  if (fns[this->user_xdisp ? 1 : 2] != fns[0])
631  {
632  if (user_xdisp && fns[2] == fns[1])
633  return;
634  fns[this->user_xdisp ? 1 : 2]->push_transform(transform);
635  }
636  }
637  }
638 
639  template<typename LinearizerDataDimensions>
640  void ThreadLinearizerMultidimensional<LinearizerDataDimensions>::pop_transforms()
641  {
642  fns[0]->pop_transform();
643 
644  if (user_xdisp)
645  {
646  if (fns[1] != fns[0])
647  fns[1]->pop_transform();
648  }
649  if (user_ydisp)
650  {
651  if (fns[this->user_xdisp ? 1 : 2] != fns[0])
652  {
653  if (user_xdisp && fns[2] == fns[1])
654  return;
655  fns[this->user_xdisp ? 1 : 2]->pop_transform();
656  }
657  }
658  }
659 
660  template<typename LinearizerDataDimensions>
661  int ThreadLinearizerMultidimensional<LinearizerDataDimensions>::peek_vertex(int p1, int p2)
662  {
663  // search for a vertex with parents p1, p2
664  int index = hash(p1 > p2 ? p2 : p1, p1 > p2 ? p1 : p2);
665  int i = hash_table[index];
666  while (i >= 0)
667  {
668  if (info[i][0] == p1 && info[i][1] == p2)
669  return i;
670  i = info[i][2];
671  }
672  return -1;
673  }
674 
675  template<typename LinearizerDataDimensions>
676  int ThreadLinearizerMultidimensional<LinearizerDataDimensions>::hash(int p1, int p2)
677  {
679  return (984120265 * p1 + 125965121 * p2) & (vertex_size - 1);
680  }
681 
682  template<typename LinearizerDataDimensions>
683  void ThreadLinearizerMultidimensional<LinearizerDataDimensions>::process_edge(int iv1, int iv2, int marker)
684  {
685  int mid = peek_vertex(iv1, iv2);
686  if (mid != -1)
687  {
688  process_edge(iv1, mid, marker);
689  process_edge(mid, iv2, marker);
690  }
691  else
692  add_edge(iv1, iv2, marker);
693  }
694 
695  template<typename LinearizerDataDimensions>
696  int ThreadLinearizerMultidimensional<LinearizerDataDimensions>::get_vertex(int p1, int p2, double x, double y, double* value)
697  {
698  // search for an existing vertex
699  if (p1 > p2)
700  std::swap(p1, p2);
701  int index = this->hash(p1, p2);
702  int i = 0;
703  if (index < this->vertex_count)
704  {
705  i = this->hash_table[index];
706  while (i >= 0 && i < this->vertex_count)
707  {
708  if ((this->info[i][0] == p1 && this->info[i][1] == p2)
709  && (fabs(x - this->vertices[i][0]) < Hermes::HermesEpsilon)
710  && (fabs(y - this->vertices[i][1]) < Hermes::HermesEpsilon))
711  {
712  bool check_value = true;
713  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
714  {
715  if (fabs(value[k] - this->vertices[i][2 + k] / value[k]) > vertex_relative_tolerance)
716  check_value = false;
717  }
718  if (check_value)
719  return i;
720  }
721  // note that we won't return a vertex with a different value than the required one;
722  // this takes care for discontinuities in the solution, where more vertices
723  // with different values will be created
724  i = info[i][2];
725  }
726  }
727 
728  // if not found, create a new_ one
729  i = add_vertex();
730 
731  this->vertices[i][0] = x;
732  this->vertices[i][1] = y;
733  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
734  this->vertices[i][2 + k] = value[k];
735  this->info[i][0] = p1;
736  this->info[i][1] = p2;
737  this->info[i][2] = hash_table[index];
738  this->hash_table[index] = i;
739  return i;
740  }
741 
742  template<typename LinearizerDataDimensions>
743  int ThreadLinearizerMultidimensional<LinearizerDataDimensions>::add_vertex()
744  {
745  if (this->vertex_count >= this->vertex_size)
746  {
747  int new_vertex_size = std::ceil(this->vertex_size * 1.5);
748 
749  this->vertices = realloc_with_check<ThreadLinearizerMultidimensional, typename LinearizerDataDimensions::vertex_t>(this->vertices, new_vertex_size, this);
750  this->info = realloc_with_check<ThreadLinearizerMultidimensional, internal_vertex_info_t>(this->info, new_vertex_size, this);
751  this->hash_table = realloc_with_check<ThreadLinearizerMultidimensional, int>(this->hash_table, new_vertex_size, this);
752  memset(this->hash_table + this->vertex_size, 0xff, sizeof(int)* (new_vertex_size - this->vertex_size));
753 
754  this->vertex_size = new_vertex_size;
755  }
756  return this->vertex_count++;
757  }
758 
759  template<typename LinearizerDataDimensions>
760  void ThreadLinearizerMultidimensional<LinearizerDataDimensions>::add_edge(int iv1, int iv2, int marker)
761  {
762  if (edges_count >= edges_size)
763  {
764  this->edges_size = std::ceil(this->edges_size * 1.5);
765  this->edges = realloc_with_check<ThreadLinearizerMultidimensional, typename LinearizerDataDimensions::edge_t>(this->edges, this->edges_size, this);
766  this->edge_markers = realloc_with_check<ThreadLinearizerMultidimensional, int>(this->edge_markers, this->edges_size, this);
767  }
768 
769  this->edges[edges_count][0][0] = this->vertices[iv1][0];
770  this->edges[edges_count][0][1] = this->vertices[iv1][1];
771  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
772  this->edges[edges_count][0][2 + k] = this->vertices[iv1][2 + k];
773  this->edges[edges_count][1][0] = this->vertices[iv2][0];
774  this->edges[edges_count][1][1] = this->vertices[iv2][1];
775  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
776  this->edges[edges_count][1][2 + k] = this->vertices[iv2][2 + k];
777  this->edge_markers[edges_count++] = marker;
778  }
779 
780  template<typename LinearizerDataDimensions>
781  void ThreadLinearizerMultidimensional<LinearizerDataDimensions>::add_triangle(int iv0, int iv1, int iv2, int marker)
782  {
783  if (triangle_count >= triangle_size)
784  {
785  this->triangle_size = std::ceil(this->triangle_size * 1.5);
786  if (this->linearizerOutputType == OpenGL)
787  this->triangles = realloc_with_check<ThreadLinearizerMultidimensional, typename LinearizerDataDimensions::triangle_t>(this->triangles, this->triangle_size, this);
788  else
789  this->triangle_indices = realloc_with_check<ThreadLinearizerMultidimensional, triangle_indices_t>(this->triangle_indices, this->triangle_size, this);
790  this->triangle_markers = realloc_with_check<ThreadLinearizerMultidimensional, int>(this->triangle_markers, this->triangle_size, this);
791  }
792 
793  if (this->linearizerOutputType == OpenGL)
794  {
795  this->triangles[triangle_count][0][0] = this->vertices[iv0][0];
796  this->triangles[triangle_count][0][1] = this->vertices[iv0][1];
797  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
798  this->triangles[triangle_count][0][2 + k] = this->vertices[iv0][2 + k];
799  this->triangles[triangle_count][1][0] = this->vertices[iv1][0];
800  this->triangles[triangle_count][1][1] = this->vertices[iv1][1];
801  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
802  this->triangles[triangle_count][1][2 + k] = this->vertices[iv1][2 + k];
803  this->triangles[triangle_count][2][0] = this->vertices[iv2][0];
804  this->triangles[triangle_count][2][1] = this->vertices[iv2][1];
805  for (int k = 0; k < LinearizerDataDimensions::dimension; k++)
806  this->triangles[triangle_count][2][2 + k] = this->vertices[iv2][2 + k];
807  }
808  else
809  {
810  this->triangle_indices[triangle_count][0] = iv0;
811  this->triangle_indices[triangle_count][1] = iv1;
812  this->triangle_indices[triangle_count][2] = iv2;
813  }
814 
815  this->triangle_markers[triangle_count++] = marker;
816  }
817 
818  template class HERMES_API ThreadLinearizerMultidimensional < ScalarLinearizerDataDimensions<double> > ;
819  template class HERMES_API ThreadLinearizerMultidimensional < VectorLinearizerDataDimensions<double> > ;
820  template class HERMES_API ThreadLinearizerMultidimensional < ScalarLinearizerDataDimensions<float> > ;
821  template class HERMES_API ThreadLinearizerMultidimensional < VectorLinearizerDataDimensions<float> > ;
822  }
823  }
824 }
Hermes
Definition: adapt.h:24
Hermes::Hermes2D::Solution::copy
virtual void copy(const MeshFunction< Scalar > *sln)
Copy from sln to this instance.
Definition: solution.cpp:179
Hermes::Hermes2D::Views::internal_vertex_info_t
int3 internal_vertex_info_t
Typedefs used throughout the Linearizer functionality.
Definition: linearizer_utils.h:109
LINEARIZER_DIRECTIONAL_QUAD_REFINEMENT_REQUIREMENT
#define LINEARIZER_DIRECTIONAL_QUAD_REFINEMENT_REQUIREMENT
We refine a quad directionally (horizontally, vertically) only if the error in one direction is this ...
Definition: linearizer_utils.h:46
thread_linearizer.h
File containing ThreadLinearizerMultidimensional class.
H2D_MAX_NUMBER_VERTICES
#define H2D_MAX_NUMBER_VERTICES
A maximum number of vertices of an element.
Definition: global.h:32
Hermes::Hermes2D::Solution::get_type
SolutionType get_type() const
Returns solution type.
Definition: solution.h:143
xml_schema::error
::xsd::cxx::tree::error< char > error
Error condition.
Definition: solution_h2d_xml.h:512
MAX_LINEARIZER_DIVISION_LEVEL
#define MAX_LINEARIZER_DIVISION_LEVEL
Very important constant putting an upper bound on the maximum number of successive element division (...
Definition: linearizer_utils.h:50