Hermes2D  2.0
linearizerTest.h
1 void Linearizer::process_triangle(MeshFunction<double>** fns, int iv0, int iv1, int iv2, int level,
2  double* val, double* phx, double* phy, int* idx, bool curved)
3  {
4  double midval[3][3];
5 
6  if(level < LIN_MAX_LEVEL)
7  {
8  int i;
9  if(!(level & 1))
10  {
11  // obtain solution values
12  fns[0]->set_quad_order(1, item);
13  val = fns[0]->get_values(component, value_type);
14  if(auto_max)
15  for (i = 0; i < lin_np_tri[1]; i++)
16  {
17  double v = val[i];
18 #pragma omp critical(max)
19  if(finite(v) && fabs(v) > max)
20  max = fabs(v);
21  }
22  idx = tri_indices[0];
23 
24  if(curved)
25  {
26  // obtain physical element coordinates
27  RefMap* refmap = fns[0]->get_refmap();
28  phx = refmap->get_phys_x(1);
29  phy = refmap->get_phys_y(1);
30 
31  double* dx = NULL;
32  double* dy = NULL;
33  if(this->xdisp != NULL)
34  {
35  fns[1]->set_quad_order(1, H2D_FN_VAL);
36  dx = fns[1]->get_fn_values();
37  }
38  if(this->ydisp != NULL)
39  {
40  fns[2]->set_quad_order(1, H2D_FN_VAL);
41  dy = fns[2]->get_fn_values();
42  }
43  for (i = 0; i < lin_np_tri[1]; i++)
44  {
45  if(this->xdisp != NULL)
46  phx[i] += dmult*dx[i];
47  if(this->ydisp != NULL)
48  phy[i] += dmult*dy[i];
49  }
50 
51  }
52  }
53 
54  // obtain linearized values and coordinates at the midpoints
55  for (i = 0; i < 3; i++)
56  {
57  midval[i][0] = (verts[iv0][i] + verts[iv1][i])*0.5;
58  midval[i][1] = (verts[iv1][i] + verts[iv2][i])*0.5;
59  midval[i][2] = (verts[iv2][i] + verts[iv0][i])*0.5;
60  };
61 
62  // determine whether or not to split the element
63  bool split;
64  if(eps >= 1.0)
65  {
66  // if eps > 1, the user wants a fixed number of refinements (no adaptivity)
67  split = (level + 5 < eps);
68  }
69  else
70  {
71  if(!auto_max && fabs(verts[iv0][2]) > max && fabs(verts[iv1][2]) > max && fabs(verts[iv2][2]) > max)
72  {
73  // do not split if the whole triangle is above the specified maximum value
74  split = false;
75  }
76  else
77  {
78  // calculate the approximate error of linearizing the normalized solution
79  double err = fabs(val[idx[0]] - midval[2][0]) +
80  fabs(val[idx[1]] - midval[2][1]) +
81  fabs(val[idx[2]] - midval[2][2]);
82  split = !finite(err) || err > max*3*eps;
83  }
84 
85  // do the same for the curvature
86  if(curved)
87  {
88  for (i = 0; i < 3; i++)
89  if(sqr(phx[idx[i]] - midval[0][i]) + sqr(phy[idx[i]] - midval[1][i]) > sqr(fns[0]->get_active_element()->get_diameter()*this->get_curvature_epsilon()))
90  {
91  split = true;
92  break;
93  }
94  }
95 
96  // do extra tests at level 0, so as not to miss some functions with zero error at edge midpoints
97  if(level == 0 && !split)
98  {
99  split = (fabs(val[8] - 0.5*(midval[2][0] + midval[2][1])) +
100  fabs(val[9] - 0.5*(midval[2][1] + midval[2][2])) +
101  fabs(val[4] - 0.5*(midval[2][2] + midval[2][0]))) > max*3*eps;
102  }
103  }
104 
105  // split the triangle if the error is too large, otherwise produce a linear triangle
106  if(split)
107  {
108  if(curved)
109  for (i = 0; i < 3; i++)
110  {
111  midval[0][i] = phx[idx[i]];
112  midval[1][i] = phy[idx[i]];
113  }
114 
115  // obtain mid-edge vertices
116  int mid0 = get_vertex(iv0, iv1, midval[0][0], midval[1][0], val[idx[0]]);
117  int mid1 = get_vertex(iv1, iv2, midval[0][1], midval[1][1], val[idx[1]]);
118  int mid2 = get_vertex(iv2, iv0, midval[0][2], midval[1][2], val[idx[2]]);
119 
120  // recur to sub-elements
121  fns[0]->push_transform(0);
122  if(this->xdisp != NULL)
123  if(fns[1] != fns[0])
124  fns[1]->push_transform(0);
125  if(this->ydisp != NULL)
126  if(fns[2] != fns[1])
127  fns[2]->push_transform(0);
128  process_triangle(fns, iv0, mid0, mid2, level + 1, val, phx, phy, tri_indices[1], curved);
129  fns[0]->pop_transform();
130  if(this->xdisp != NULL)
131  if(fns[1] != fns[0])
132  fns[1]->pop_transform();
133  if(this->ydisp != NULL)
134  if(fns[2] != fns[1])
135  fns[2]->pop_transform();
136 
137  fns[0]->push_transform(1);
138  if(this->xdisp != NULL)
139  if(fns[1] != fns[0])
140  fns[1]->push_transform(1);
141  if(this->ydisp != NULL)
142  if(fns[2] != fns[1])
143  fns[2]->push_transform(1);
144  process_triangle(fns, mid0, iv1, mid1, level + 1, val, phx, phy, tri_indices[2], curved);
145  fns[0]->pop_transform();
146  if(this->xdisp != NULL)
147  if(fns[1] != fns[0])
148  fns[1]->pop_transform();
149  if(this->ydisp != NULL)
150  if(fns[2] != fns[1])
151  fns[2]->pop_transform();
152 
153  fns[0]->push_transform(2);
154  if(this->xdisp != NULL)
155  if(fns[1] != fns[0])
156  fns[1]->push_transform(2);
157  if(this->ydisp != NULL)
158  if(fns[2] != fns[1])
159  fns[2]->push_transform(2);
160  process_triangle(fns, mid2, mid1, iv2, level + 1, val, phx, phy, tri_indices[3], curved);
161  fns[0]->pop_transform();
162  if(this->xdisp != NULL)
163  if(fns[1] != fns[0])
164  fns[1]->pop_transform();
165  if(this->ydisp != NULL)
166  if(fns[2] != fns[1])
167  fns[2]->pop_transform();
168 
169  fns[0]->push_transform(3);
170  if(this->xdisp != NULL)
171  if(fns[1] != fns[0])
172  fns[1]->push_transform(3);
173  if(this->ydisp != NULL)
174  if(fns[2] != fns[1])
175  fns[2]->push_transform(3);
176  process_triangle(fns, mid1, mid2, mid0, level + 1, val, phx, phy, tri_indices[4], curved);
177  fns[0]->pop_transform();
178  if(this->xdisp != NULL)
179  if(fns[1] != fns[0])
180  fns[1]->pop_transform();
181  if(this->ydisp != NULL)
182  if(fns[2] != fns[1])
183  fns[2]->pop_transform();
184  return;
185  }
186  }
187 
188  // no splitting: output a linear triangle
189  add_triangle(iv0, iv1, iv2, fns[0]->get_active_element()->marker);
190  }
191 
192  void Linearizer::set_curvature_epsilon(double curvature_epsilon)
193  {
194  this->curvature_epsilon = curvature_epsilon;
195  }
196 
197  double Linearizer::get_curvature_epsilon()
198  {
199  return this->curvature_epsilon;
200  }
201 
202  void Linearizer::process_quad(MeshFunction<double>** fns, int iv0, int iv1, int iv2, int iv3, int level,
203  double* val, double* phx, double* phy, int* idx, bool curved)
204  {
205  double midval[3][5];
206 
207  // try not to split through the vertex with the largest value
208  int a = (verts[iv0][2] > verts[iv1][2]) ? iv0 : iv1;
209  int b = (verts[iv2][2] > verts[iv3][2]) ? iv2 : iv3;
210  a = (verts[a][2] > verts[b][2]) ? a : b;
211  int flip = (a == iv1 || a == iv3) ? 1 : 0;
212 
213  if(level < LIN_MAX_LEVEL)
214  {
215  int i;
216  if(!(level & 1)) // this is an optimization: do the following only every other time
217  {
218  // obtain solution values
219  fns[0]->set_quad_order(1, item);
220  val = fns[0]->get_values(component, value_type);
221  if(auto_max)
222  for (i = 0; i < lin_np_quad[1]; i++)
223  {
224  double v = val[i];
225  if(finite(v) && fabs(v) > max)
226 #pragma omp critical(max)
227  if(finite(v) && fabs(v) > max)
228  max = fabs(v);
229  }
230 
231  // This is just to make some sense.
232  if(fabs(max) < 1E-10)
233  max = 1E-10;
234 
235  idx = quad_indices[0];
236 
237  if(curved)
238  {
239  RefMap* refmap = fns[0]->get_refmap();
240  phx = refmap->get_phys_x(1);
241  phy = refmap->get_phys_y(1);
242 
243  double* dx = NULL;
244  double* dy = NULL;
245 
246  if(this->xdisp != NULL)
247  fns[1]->set_quad_order(1, H2D_FN_VAL);
248  if(this->ydisp != NULL)
249  fns[2]->set_quad_order(1, H2D_FN_VAL);
250  if(this->xdisp != NULL)
251  dx = fns[1]->get_fn_values();
252  if(this->ydisp != NULL)
253  dy = fns[2]->get_fn_values();
254  for (i = 0; i < lin_np_quad[1]; i++)
255  {
256  if(this->xdisp != NULL)
257  phx[i] += dmult*dx[i];
258  if(this->ydisp != NULL)
259  phy[i] += dmult*dy[i];
260  }
261  }
262  }
263 
264  // obtain linearized values and coordinates at the midpoints
265  for (i = 0; i < 3; i++)
266  {
267  midval[i][0] = (verts[iv0][i] + verts[iv1][i]) * 0.5;
268  midval[i][1] = (verts[iv1][i] + verts[iv2][i]) * 0.5;
269  midval[i][2] = (verts[iv2][i] + verts[iv3][i]) * 0.5;
270  midval[i][3] = (verts[iv3][i] + verts[iv0][i]) * 0.5;
271  midval[i][4] = (midval[i][0] + midval[i][2]) * 0.5;
272  };
273 
274  // the value of the middle point is not the average of the four vertex values, since quad == 2 triangles
275  midval[2][4] = flip ? (verts[iv0][2] + verts[iv2][2]) * 0.5 : (verts[iv1][2] + verts[iv3][2]) * 0.5;
276 
277  // determine whether or not to split the element
278  int split;
279  if(eps >= 1.0)
280  {
281  // if eps > 1, the user wants a fixed number of refinements (no adaptivity)
282  split = (level < eps) ? 3 : 0;
283  }
284  else
285  {
286  if(!auto_max && fabs(verts[iv0][2]) > max && fabs(verts[iv1][2]) > max
287  && fabs(verts[iv2][2]) > max && fabs(verts[iv3][2]) > max)
288  {
289  // do not split if the whole quad is above the specified maximum value
290  split = 0;
291  }
292  else
293  {
294  // calculate the approximate error of linearizing the normalized solution
295  double herr = fabs(val[idx[1]] - midval[2][1]) + fabs(val[idx[3]] - midval[2][3]);
296  double verr = fabs(val[idx[0]] - midval[2][0]) + fabs(val[idx[2]] - midval[2][2]);
297  double err = fabs(val[idx[4]] - midval[2][4]) + herr + verr;
298  split = (!finite(err) || err > max*4*eps) ? 3 : 0;
299 
300  // decide whether to split horizontally or vertically only
301  if(level > 0 && split)
302  {
303  if(herr > 5*verr)
304  split = 1; // h-split
305  else if(verr > 5*herr)
306  split = 2; // v-split
307  }
308  }
309 
310  // also decide whether to split because of the curvature
311  if(split != 3 && curved)
312  {
313  double cm2 = sqr(fns[0]->get_active_element()->get_diameter()*this->get_curvature_epsilon());
314  if(sqr(phx[idx[1]] - midval[0][1]) + sqr(phy[idx[1]] - midval[1][1]) > cm2 ||
315  sqr(phx[idx[3]] - midval[0][3]) + sqr(phy[idx[3]] - midval[1][3]) > cm2) split |= 1;
316  if(sqr(phx[idx[0]] - midval[0][0]) + sqr(phy[idx[0]] - midval[1][0]) > cm2 ||
317  sqr(phx[idx[2]] - midval[0][2]) + sqr(phy[idx[2]] - midval[1][2]) > cm2) split |= 2;
318  }
319 
320  // do extra tests at level 0, so as not to miss some functions with zero error at edge midpoints
321  if(level == 0 && !split)
322  {
323  split = ((fabs(val[13] - 0.5*(midval[2][0] + midval[2][1])) +
324  fabs(val[17] - 0.5*(midval[2][1] + midval[2][2])) +
325  fabs(val[20] - 0.5*(midval[2][2] + midval[2][3])) +
326  fabs(val[9] - 0.5*(midval[2][3] + midval[2][0]))) > max*4*eps) ? 3 : 0;
327  }
328  }
329 
330  // split the quad if the error is too large, otherwise produce two linear triangles
331  if(split)
332  {
333  if(curved)
334  for (i = 0; i < 5; i++)
335  {
336  midval[0][i] = phx[idx[i]];
337  midval[1][i] = phy[idx[i]];
338  }
339 
340  // obtain mid-edge and mid-element vertices
341  int mid0, mid1, mid2, mid3, mid4;
342  if(split != 1) mid0 = get_vertex(iv0, iv1, midval[0][0], midval[1][0], val[idx[0]]);
343  if(split != 2) mid1 = get_vertex(iv1, iv2, midval[0][1], midval[1][1], val[idx[1]]);
344  if(split != 1) mid2 = get_vertex(iv2, iv3, midval[0][2], midval[1][2], val[idx[2]]);
345  if(split != 2) mid3 = get_vertex(iv3, iv0, midval[0][3], midval[1][3], val[idx[3]]);
346  if(split == 3) mid4 = get_vertex(mid0, mid2, midval[0][4], midval[1][4], val[idx[4]]);
347 
348  // recur to sub-elements
349  if(split == 3)
350  {
351  fns[0]->push_transform(0);
352  if(this->xdisp != NULL)
353  if(fns[1] != fns[0])
354  fns[1]->push_transform(0);
355  if(this->ydisp != NULL)
356  if(fns[2] != fns[1])
357  fns[2]->push_transform(0);
358  process_quad(fns, iv0, mid0, mid4, mid3, level + 1, val, phx, phy, quad_indices[1], curved);
359  fns[0]->pop_transform();
360  if(this->xdisp != NULL)
361  if(fns[1] != fns[0])
362  fns[1]->pop_transform();
363  if(this->ydisp != NULL)
364  if(fns[2] != fns[1])
365  fns[2]->pop_transform();
366 
367  fns[0]->push_transform(1);
368  if(this->xdisp != NULL)
369  if(fns[1] != fns[0])
370  fns[1]->push_transform(1);
371  if(this->ydisp != NULL)
372  if(fns[2] != fns[1])
373  fns[2]->push_transform(1);
374  process_quad(fns, mid0, iv1, mid1, mid4, level + 1, val, phx, phy, quad_indices[2], curved);
375  fns[0]->pop_transform();
376  if(this->xdisp != NULL)
377  if(fns[1] != fns[0])
378  fns[1]->pop_transform();
379  if(this->ydisp != NULL)
380  if(fns[2] != fns[1])
381  fns[2]->pop_transform();
382 
383  fns[0]->push_transform(2);
384  if(this->xdisp != NULL)
385  if(fns[1] != fns[0])
386  fns[1]->push_transform(2);
387  if(this->ydisp != NULL)
388  if(fns[2] != fns[1])
389  fns[2]->push_transform(2);
390  process_quad(fns, mid4, mid1, iv2, mid2, level + 1, val, phx, phy, quad_indices[3], curved);
391  fns[0]->pop_transform();
392  if(this->xdisp != NULL)
393  if(fns[1] != fns[0])
394  fns[1]->pop_transform();
395  if(this->ydisp != NULL)
396  if(fns[2] != fns[1])
397  fns[2]->pop_transform();
398 
399  fns[0]->push_transform(3);
400  if(this->xdisp != NULL)
401  if(fns[1] != fns[0])
402  fns[1]->push_transform(3);
403  if(this->ydisp != NULL)
404  if(fns[2] != fns[1])
405  fns[2]->push_transform(3);
406  process_quad(fns, mid3, mid4, mid2, iv3, level + 1, val, phx, phy, quad_indices[4], curved);
407  fns[0]->pop_transform();
408  if(this->xdisp != NULL)
409  if(fns[1] != fns[0])
410  fns[1]->pop_transform();
411  if(this->ydisp != NULL)
412  if(fns[2] != fns[1])
413  fns[2]->pop_transform();
414  }
415  else
416  if(split == 1) // h-split
417  {
418  fns[0]->push_transform(4);
419  if(this->ydisp != NULL)
420  if(fns[1] != fns[0])
421  fns[1]->push_transform(4);
422  if(this->ydisp != NULL)
423  if(fns[2] != fns[1])
424  fns[2]->push_transform(4);
425  process_quad(fns, iv0, iv1, mid1, mid3, level + 1, val, phx, phy, quad_indices[5], curved);
426  fns[0]->pop_transform();
427  if(this->xdisp != NULL)
428  if(fns[1] != fns[0])
429  fns[1]->pop_transform();
430  if(this->ydisp != NULL)
431  if(fns[2] != fns[1])
432  fns[2]->pop_transform();
433 
434  fns[0]->push_transform(5);
435  if(this->xdisp != NULL)
436  if(fns[1] != fns[0])
437  fns[1]->push_transform(5);
438  if(this->ydisp != NULL)
439  if(fns[2] != fns[1])
440  fns[2]->push_transform(5);
441  process_quad(fns, mid3, mid1, iv2, iv3, level + 1, val, phx, phy, quad_indices[6], curved);
442  fns[0]->pop_transform();
443  if(this->xdisp != NULL)
444  if(fns[1] != fns[0])
445  fns[1]->pop_transform();
446  if(this->ydisp != NULL)
447  if(fns[2] != fns[1])
448  fns[2]->pop_transform();
449  }
450  else // v-split
451  {
452  fns[0]->push_transform(6);
453  if(this->xdisp != NULL)
454  if(fns[1] != fns[0])
455  fns[1]->push_transform(6);
456  if(this->ydisp != NULL)
457  if(fns[2] != fns[1])
458  fns[2]->push_transform(6);
459  process_quad(fns, iv0, mid0, mid2, iv3, level + 1, val, phx, phy, quad_indices[7], curved);
460  fns[0]->pop_transform();
461  if(this->xdisp != NULL)
462  if(fns[1] != fns[0])
463  fns[1]->pop_transform();
464  if(this->ydisp != NULL)
465  if(fns[2] != fns[1])
466  fns[2]->pop_transform();
467 
468  fns[0]->push_transform(7);
469  if(this->xdisp != NULL)
470  if(fns[1] != fns[0])
471  fns[1]->push_transform(7);
472  if(this->ydisp != NULL)
473  if(fns[2] != fns[1])
474  fns[2]->push_transform(7);
475  process_quad(fns, mid0, iv1, iv2, mid2, level + 1, val, phx, phy, quad_indices[8], curved);
476  fns[0]->pop_transform();
477  if(this->xdisp != NULL)
478  if(fns[1] != fns[0])
479  fns[1]->pop_transform();
480  if(this->ydisp != NULL)
481  if(fns[2] != fns[1])
482  fns[2]->pop_transform();
483  }
484  return;
485  }
486  }
487 
488  // output two linear triangles,
489  if(!flip)
490  {
491  add_triangle(iv3, iv0, iv1, fns[0]->get_active_element()->marker);
492  add_triangle(iv1, iv2, iv3, fns[0]->get_active_element()->marker);
493  }
494  else
495  {
496  add_triangle(iv0, iv1, iv2, fns[0]->get_active_element()->marker);
497  add_triangle(iv2, iv3, iv0, fns[0]->get_active_element()->marker);
498  }
499  }
500 
501  void Linearizer::process_solution(MeshFunction<double>* sln, int item_, double eps)
502  {
503  // Important, sets the current caughtException to NULL.
504  this->caughtException = NULL;
505 
506  lock_data();
507  this->tick();
508 
509  // Initialization of 'global' stuff.
510  this->item = item_;
511  this->eps = eps;
512  // get the component and desired value from item.
513  if(item >= 0x40)
514  {
515  component = 1;
516  this->item >>= 6;
517  }
518  while (!(item & 1))
519  {
520  this->item >>= 1;
521  value_type++;
522  }
523  // reset the item to the value before the circus with component, value_type.
524  this->item = item_;
525 
526  // Initialization of computation stuff.
527  // sizes.
528  this->vertex_size = std::max(100 * sln->get_mesh()->get_num_elements(), std::max(this->vertex_size, 50000));
529  this->triangle_size = std::max(150 * sln->get_mesh()->get_num_elements(), std::max(this->triangle_size, 75000));
530  this->edges_size = std::max(100 * sln->get_mesh()->get_num_elements(), std::max(this->edges_size, 50000));
531  // counts.
532  this->vertex_count = 0;
533  this->triangle_count = 0;
534  this->edges_count = 0;
535  // reuse or allocate vertex, triangle and edge arrays.
536  this->verts = (double3*) realloc(this->verts, sizeof(double3) * this->vertex_size);
537  this->tris = (int3*) realloc(this->tris, sizeof(int3) * this->triangle_size);
538  this->tri_markers = (int*) realloc(this->tri_markers, sizeof(int) * this->triangle_size);
539  this->edges = (int2*) realloc(this->edges, sizeof(int2) * this->edges_size);
540  this->edge_markers = (int*) realloc(this->edge_markers, sizeof(int) * this->edges_size);
541  this->info = (int4*) malloc(sizeof(int4) * this->vertex_size);
542  this->empty = false;
543  // initialize the hash table
544  this->hash_table = (int*) malloc(sizeof(int) * this->vertex_size);
545  memset(this->hash_table, 0xff, sizeof(int) * this->vertex_size);
546 
547  // select the linearization quadratures
548  Quad2D *old_quad, *old_quad_x = NULL, *old_quad_y = NULL;
549  old_quad = sln->get_quad_2d();
550  if(xdisp != NULL)
551  old_quad_x = xdisp->get_quad_2d();
552  if(ydisp != NULL)
553  old_quad_y = ydisp->get_quad_2d();
554 
555  // obtain the solution in vertices, estimate the maximum solution value
556  // meshes.
557  Hermes::vector<const Mesh*> meshes;
558  meshes.push_back(sln->get_mesh());
559  if(xdisp != NULL)
560  meshes.push_back(xdisp->get_mesh());
561  if(ydisp != NULL)
562  meshes.push_back(ydisp->get_mesh());
563 
564  // Parallelization
565  MeshFunction<double>*** fns = new MeshFunction<double>**[Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
566  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
567  {
568  fns[i] = new MeshFunction<double>*[3];
569  fns[i][0] = sln->clone();
570  fns[i][0]->set_refmap(new RefMap);
571  fns[i][0]->set_quad_2d(&g_quad_lin);
572  if(xdisp != NULL)
573  {
574  fns[i][1] = xdisp->clone();
575  //fns[i][1]->set_refmap(new RefMap);
576  fns[i][1]->set_quad_2d(&g_quad_lin);
577  }
578  if(ydisp != NULL)
579  {
580  fns[i][2] = ydisp->clone();
581  //fns[i][2]->set_refmap(new RefMap);
582  fns[i][2]->set_quad_2d(&g_quad_lin);
583  }
584  }
585 
586  Transformable*** trfs = new Transformable**[Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
587  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
588  {
589  trfs[i] = new Transformable*[3];
590  trfs[i][0] = fns[i][0];
591  if(xdisp != NULL)
592  trfs[i][1] = fns[i][1];
593  if(ydisp != NULL)
594  trfs[i][2] = fns[i][2];
595  }
596 
597  Traverse trav_masterMax(true);
598  unsigned int num_states = trav_masterMax.get_num_states(meshes);
599 
600  trav_masterMax.begin(meshes.size(), &(meshes.front()));
601 
602  Traverse* trav = new Traverse[Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
603 
604  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
605  {
606  trav[i].begin(meshes.size(), &(meshes.front()), trfs[i]);
607  trav[i].stack = trav_masterMax.stack;
608  }
609 
610  int state_i;
611 
612 #define CHUNKSIZE 1
613  int num_threads_used = Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads);
614 #pragma omp parallel shared(trav_masterMax) private(state_i) num_threads(num_threads_used)
615  {
616 #pragma omp for schedule(dynamic, CHUNKSIZE)
617  for(state_i = 0; state_i < num_states; state_i++)
618  {
619  try
620  {
621  Traverse::State current_state;
622 #pragma omp critical(get_next_state)
623  current_state = trav[omp_get_thread_num()].get_next_state(&trav_masterMax.top, &trav_masterMax.id);
624 
625  fns[omp_get_thread_num()][0]->set_quad_order(0, this->item);
626  double* val = fns[omp_get_thread_num()][0]->get_values(component, value_type);
627 
628  for (unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
629  {
630  double f = val[i];
631 #pragma omp critical (max)
632  if(this->auto_max && finite(f) && fabs(f) > this->max)
633  this->max = fabs(f);
634  }
635  }
636  catch(Hermes::Exceptions::Exception& e)
637  {
638  if(this->caughtException == NULL)
639  this->caughtException = e.clone();
640  }
641  catch(std::exception& e)
642  {
643  if(this->caughtException == NULL)
644  this->caughtException = new Hermes::Exceptions::Exception(e.what());
645  }
646  }
647  }
648 
649  trav_masterMax.finish();
650  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
651  trav[i].finish();
652  delete [] trav;
653 
654  Traverse trav_master(true);
655  num_states = trav_master.get_num_states(meshes);
656 
657  trav_master.begin(meshes.size(), &(meshes.front()));
658 
659  trav = new Traverse[Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
660 
661  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
662  {
663  trav[i].begin(meshes.size(), &(meshes.front()), trfs[i]);
664  trav[i].stack = trav_master.stack;
665  }
666 
667 #pragma omp parallel shared(trav_master) private(state_i) num_threads(num_threads_used)
668  {
669 #pragma omp for schedule(dynamic, CHUNKSIZE)
670  for(state_i = 0; state_i < num_states; state_i++)
671  {
672  try
673  {
674  Traverse::State current_state;
675 
676 #pragma omp critical (get_next_state)
677  current_state = trav[omp_get_thread_num()].get_next_state(&trav_master.top, &trav_master.id);
678 
679  fns[omp_get_thread_num()][0]->set_quad_order(0, this->item);
680  double* val = fns[omp_get_thread_num()][0]->get_values(component, value_type);
681  if(val == NULL)
682  {
683  delete [] trav;
684  throw Hermes::Exceptions::Exception("Item not defined in the solution.");
685  }
686 
687  if(xdisp != NULL)
688  fns[omp_get_thread_num()][1]->set_quad_order(0, H2D_FN_VAL);
689  if(ydisp != NULL)
690  fns[omp_get_thread_num()][2]->set_quad_order(0, H2D_FN_VAL);
691 
692  double *dx = NULL;
693  double *dy = NULL;
694  if(xdisp != NULL)
695  dx = fns[omp_get_thread_num()][1]->get_fn_values();
696  if(ydisp != NULL)
697  dy = fns[omp_get_thread_num()][2]->get_fn_values();
698 
699  int iv[H2D_MAX_NUMBER_VERTICES];
700  for (unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
701  {
702  double f = val[i];
703  double x_disp = fns[omp_get_thread_num()][0]->get_refmap()->get_phys_x(0)[i];
704  double y_disp = fns[omp_get_thread_num()][0]->get_refmap()->get_phys_y(0)[i];
705  if(this->xdisp != NULL)
706  x_disp += dmult * dx[i];
707  if(this->ydisp != NULL)
708  y_disp += dmult * dy[i];
709 
710  iv[i] = this->get_vertex(-fns[omp_get_thread_num()][0]->get_active_element()->vn[i]->id, -fns[omp_get_thread_num()][0]->get_active_element()->vn[i]->id, x_disp, y_disp, f);
711  }
712 
713  // recur to sub-elements
714  if(current_state.e[0]->is_triangle())
715  process_triangle(fns[omp_get_thread_num()], iv[0], iv[1], iv[2], 0, NULL, NULL, NULL, NULL, current_state.e[0]->is_curved());
716  else
717  process_quad(fns[omp_get_thread_num()], iv[0], iv[1], iv[2], iv[3], 0, NULL, NULL, NULL, NULL, current_state.e[0]->is_curved());
718 
719  for (unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
720  process_edge(iv[i], iv[current_state.e[0]->next_vert(i)], current_state.e[0]->en[i]->marker);
721  }
722  catch(Hermes::Exceptions::Exception& e)
723  {
724  if(this->caughtException == NULL)
725  this->caughtException = e.clone();
726  }
727  catch(std::exception& e)
728  {
729  if(this->caughtException == NULL)
730  this->caughtException = new Hermes::Exceptions::Exception(e.what());
731  }
732  }
733  }
734 
735  trav_master.finish();
736  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
737  {
738  trav[i].finish();
739  for(unsigned int j = 0; j < (1 + (xdisp != NULL? 1 : 0) + (ydisp != NULL ? 1 : 0)); j++)
740  delete fns[i][j];
741  delete [] fns[i];
742  delete [] trfs[i];
743  }
744  delete [] fns;
745  delete [] trfs;
746  delete [] trav;
747 
748  // for contours, without regularization.
749  this->tris_contours = (int3*) realloc(this->tris_contours, sizeof(int3) * this->triangle_count);
750  memcpy(this->tris_contours, this->tris, this->triangle_count * sizeof(int3));
751  triangle_contours_count = this->triangle_count;
752 
753  if(this->caughtException != NULL)
754  {
755  this->unlock_data();
756  ::free(hash_table);
757  ::free(info);
758  throw *(this->caughtException);
759  }
760 
761  // regularize the linear mesh
762  for (int i = 0; i < this->triangle_count; i++)
763  {
764  int iv0 = tris[i][0], iv1 = tris[i][1], iv2 = tris[i][2];
765 
766  int mid0 = peek_vertex(iv0, iv1);
767  int mid1 = peek_vertex(iv1, iv2);
768  int mid2 = peek_vertex(iv2, iv0);
769  if(mid0 >= 0 || mid1 >= 0 || mid2 >= 0)
770  {
771  this->del_slot = i;
772  regularize_triangle(iv0, iv1, iv2, mid0, mid1, mid2, tri_markers[i]);
773  }
774  }
775 
776  find_min_max();
777 
778  this->unlock_data();
779 
780  // select old quadratrues
781  sln->set_quad_2d(old_quad);
782  if(user_xdisp)
783  xdisp->set_quad_2d(old_quad_x);
784  else
785  delete xdisp;
786  if(user_ydisp)
787  ydisp->set_quad_2d(old_quad_y);
788  else
789  delete ydisp;
790 
791  // clean up
792  ::free(hash_table);
793  ::free(info);
794  }