Hermes2D  2.0
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 "linearizer.h"
17 #include "refmap.h"
18 #include "traverse.h"
19 #include "exact_solution.h"
20 #include "api2d.h"
21 
22 namespace Hermes
23 {
24  namespace Hermes2D
25  {
26  namespace Views
27  {
28  double3 lin_pts_0_tri[] =
29  {
30  { -1.0, -1.0, 0.0 },
31  { 1.0, -1.0, 0.0 },
32  { -1.0, 1.0, 0.0 }
33  };
34 
35  double3 lin_pts_0_quad[] =
36  {
37  { -1.0, -1.0, 0.0 },
38  { 1.0, -1.0, 0.0 },
39  { 1.0, 1.0, 0.0 },
40  { -1.0, 1.0, 0.0 }
41  };
42 
43  double3 lin_pts_1_tri[12] =
44  {
45  { 0.0, -1.0, 0.0 }, // 0
46  { 0.0, 0.0, 0.0 }, // 1
47  { -1.0, 0.0, 0.0 }, // 2
48  { -0.5, -1.0, 0.0 }, // 3
49  { -0.5, -0.5, 0.0 }, // 4
50  { -1.0, -0.5, 0.0 }, // 5
51  { 0.5, -1.0, 0.0 }, // 6
52  { 0.5, -0.5, 0.0 }, // 7
53  { 0.0, -0.5, 0.0 }, // 8
54  { -0.5, 0.0, 0.0 }, // 9
55  { -0.5, 0.5, 0.0 }, // 10
56  { -1.0, 0.5, 0.0 } // 11
57  };
58 
59  double3 lin_pts_1_quad[21] =
60  {
61  { 0.0, -1.0, 0.0 }, // 0
62  { 1.0, 0.0, 0.0 }, // 1
63  { 0.0, 1.0, 0.0 }, // 2
64  { -1.0, 0.0, 0.0 }, // 3
65  { 0.0, 0.0, 0.0 }, // 4
66  { -0.5, -1.0, 0.0 }, // 5
67  { 0.0, -0.5, 0.0 }, // 6
68  { -0.5, 0.0, 0.0 }, // 7
69  { -1.0, -0.5, 0.0 }, // 8
70  { -0.5, -0.5, 0.0 }, // 9
71  { 0.5, -1.0, 0.0 }, // 10
72  { 1.0, -0.5, 0.0 }, // 11
73  { 0.5, 0.0, 0.0 }, // 12
74  { 0.5, -0.5, 0.0 }, // 13
75  { 1.0, 0.5, 0.0 }, // 14
76  { 0.5, 1.0, 0.0 }, // 15
77  { 0.0, 0.5, 0.0 }, // 16
78  { 0.5, 0.5, 0.0 }, // 17
79  { -0.5, 1.0, 0.0 }, // 18
80  { -1.0, 0.5, 0.0 }, // 19
81  { -0.5, 0.5, 0.0 } // 20
82  };
83 
84  int quad_indices[9][5] =
85  {
86  { 0, 1, 2, 3, 4 },
87  { 5, 6, 7, 8, 9 }, { 10, 11, 12, 6, 13 },
88  { 12, 14, 15, 16, 17 }, { 7, 16, 18, 19, 20 },
89  { 0, 11, 4, 8, 6 }, { 4, 14, 2, 19, 16 },
90  { 5, 4, 18, 3, 7 }, { 10, 1, 15, 4, 12 }
91  };
92 
93  int tri_indices[5][3] =
94  {
95  { 0, 1, 2 }, { 3, 4, 5 }, { 6, 7, 8 }, { 9, 10, 11 }, { 9, 4, 8 }
96  };
97 
98  int lin_np_tri[2] = { 3, 12 };
99  int lin_np_quad[2] = { 4, 21 };
100  int* lin_np[2] = { lin_np_tri, lin_np_quad };
101 
102  double3* lin_tables_tri[2] = { lin_pts_0_tri, lin_pts_1_tri };
103  double3* lin_tables_quad[2] = { lin_pts_0_quad, lin_pts_1_quad };
104  double3** lin_tables[2] = { lin_tables_tri, lin_tables_quad };
105 
106  Linearizer::Linearizer(bool auto_max) : LinearizerBase(auto_max), dmult(1.0), component(0), value_type(0), curvature_epsilon(1e-3)
107  {
108  verts = NULL;
109  xdisp = NULL;
110  user_xdisp = false;
111  ydisp = NULL;
112  user_ydisp = false;
113  tris_contours = NULL;
114  }
115 
116  void Linearizer::process_triangle(MeshFunction<double>** fns, int iv0, int iv1, int iv2, int level,
117  double* val, double* phx, double* phy, int* idx, bool curved)
118  {
119  double midval[3][3];
120 
121  if(level < LIN_MAX_LEVEL)
122  {
123  int i;
124  if(!(level & 1))
125  {
126  // obtain solution values
127  fns[0]->set_quad_order(1, item);
128  val = fns[0]->get_values(component, value_type);
129  if(auto_max)
130  for (i = 0; i < lin_np_tri[1]; i++)
131  {
132  double v = val[i];
133 #pragma omp critical(max)
134  if(finite(v) && fabs(v) > max)
135  max = fabs(v);
136  }
137  idx = tri_indices[0];
138 
139  if(curved)
140  {
141  // obtain physical element coordinates
142  RefMap* refmap = fns[0]->get_refmap();
143  phx = refmap->get_phys_x(1);
144  phy = refmap->get_phys_y(1);
145 
146  double* dx = NULL;
147  double* dy = NULL;
148  if(this->xdisp != NULL)
149  {
150  fns[1]->set_quad_order(1, H2D_FN_VAL);
151  dx = fns[1]->get_fn_values();
152  }
153  if(this->ydisp != NULL)
154  {
155  fns[2]->set_quad_order(1, H2D_FN_VAL);
156  dy = fns[2]->get_fn_values();
157  }
158  for (i = 0; i < lin_np_tri[1]; i++)
159  {
160  if(this->xdisp != NULL)
161  phx[i] += dmult*dx[i];
162  if(this->ydisp != NULL)
163  phy[i] += dmult*dy[i];
164  }
165 
166  }
167  }
168 
169  // obtain linearized values and coordinates at the midpoints
170  for (i = 0; i < 3; i++)
171  {
172  midval[i][0] = (verts[iv0][i] + verts[iv1][i])*0.5;
173  midval[i][1] = (verts[iv1][i] + verts[iv2][i])*0.5;
174  midval[i][2] = (verts[iv2][i] + verts[iv0][i])*0.5;
175  };
176 
177  // determine whether or not to split the element
178  bool split;
179  if(eps >= 1.0)
180  {
181  // if eps > 1, the user wants a fixed number of refinements (no adaptivity)
182  split = (level + 5 < eps);
183  }
184  else
185  {
186  if(!auto_max && fabs(verts[iv0][2]) > max && fabs(verts[iv1][2]) > max && fabs(verts[iv2][2]) > max)
187  {
188  // do not split if the whole triangle is above the specified maximum value
189  split = false;
190  }
191  else
192  {
193  // calculate the approximate error of linearizing the normalized solution
194  double err = fabs(val[idx[0]] - midval[2][0]) +
195  fabs(val[idx[1]] - midval[2][1]) +
196  fabs(val[idx[2]] - midval[2][2]);
197  split = !finite(err) || err > max*3*eps;
198  }
199 
200  // do the same for the curvature
201  if(!split && curved)
202  {
203  for (i = 0; i < 3; i++)
204  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()))
205  {
206  split = true;
207  break;
208  }
209  }
210 
211  // do extra tests at level 0, so as not to miss some functions with zero error at edge midpoints
212  if(level == 0 && !split)
213  {
214  split = (fabs(val[8] - 0.5*(midval[2][0] + midval[2][1])) +
215  fabs(val[9] - 0.5*(midval[2][1] + midval[2][2])) +
216  fabs(val[4] - 0.5*(midval[2][2] + midval[2][0]))) > max*3*eps;
217  }
218  }
219 
220  // split the triangle if the error is too large, otherwise produce a linear triangle
221  if(split)
222  {
223  if(curved)
224  for (i = 0; i < 3; i++)
225  {
226  midval[0][i] = phx[idx[i]];
227  midval[1][i] = phy[idx[i]];
228  }
229 
230  // obtain mid-edge vertices
231  int mid0 = get_vertex(iv0, iv1, midval[0][0], midval[1][0], val[idx[0]]);
232  int mid1 = get_vertex(iv1, iv2, midval[0][1], midval[1][1], val[idx[1]]);
233  int mid2 = get_vertex(iv2, iv0, midval[0][2], midval[1][2], val[idx[2]]);
234 
235  // recur to sub-elements
236  fns[0]->push_transform(0);
237  if(this->xdisp != NULL)
238  if(fns[1] != fns[0])
239  fns[1]->push_transform(0);
240  if(this->ydisp != NULL)
241  if(fns[2] != fns[1])
242  fns[2]->push_transform(0);
243  process_triangle(fns, iv0, mid0, mid2, level + 1, val, phx, phy, tri_indices[1], curved);
244  fns[0]->pop_transform();
245  if(this->xdisp != NULL)
246  if(fns[1] != fns[0])
247  fns[1]->pop_transform();
248  if(this->ydisp != NULL)
249  if(fns[2] != fns[1])
250  fns[2]->pop_transform();
251 
252  fns[0]->push_transform(1);
253  if(this->xdisp != NULL)
254  if(fns[1] != fns[0])
255  fns[1]->push_transform(1);
256  if(this->ydisp != NULL)
257  if(fns[2] != fns[1])
258  fns[2]->push_transform(1);
259  process_triangle(fns, mid0, iv1, mid1, level + 1, val, phx, phy, tri_indices[2], curved);
260  fns[0]->pop_transform();
261  if(this->xdisp != NULL)
262  if(fns[1] != fns[0])
263  fns[1]->pop_transform();
264  if(this->ydisp != NULL)
265  if(fns[2] != fns[1])
266  fns[2]->pop_transform();
267 
268  fns[0]->push_transform(2);
269  if(this->xdisp != NULL)
270  if(fns[1] != fns[0])
271  fns[1]->push_transform(2);
272  if(this->ydisp != NULL)
273  if(fns[2] != fns[1])
274  fns[2]->push_transform(2);
275  process_triangle(fns, mid2, mid1, iv2, level + 1, val, phx, phy, tri_indices[3], curved);
276  fns[0]->pop_transform();
277  if(this->xdisp != NULL)
278  if(fns[1] != fns[0])
279  fns[1]->pop_transform();
280  if(this->ydisp != NULL)
281  if(fns[2] != fns[1])
282  fns[2]->pop_transform();
283 
284  fns[0]->push_transform(3);
285  if(this->xdisp != NULL)
286  if(fns[1] != fns[0])
287  fns[1]->push_transform(3);
288  if(this->ydisp != NULL)
289  if(fns[2] != fns[1])
290  fns[2]->push_transform(3);
291  process_triangle(fns, mid1, mid2, mid0, level + 1, val, phx, phy, tri_indices[4], curved);
292  fns[0]->pop_transform();
293  if(this->xdisp != NULL)
294  if(fns[1] != fns[0])
295  fns[1]->pop_transform();
296  if(this->ydisp != NULL)
297  if(fns[2] != fns[1])
298  fns[2]->pop_transform();
299  return;
300  }
301  }
302 
303  // no splitting: output a linear triangle
304  add_triangle(iv0, iv1, iv2, fns[0]->get_active_element()->marker);
305  }
306 
307  void Linearizer::set_curvature_epsilon(double curvature_epsilon)
308  {
309  this->curvature_epsilon = curvature_epsilon;
310  }
311 
312  double Linearizer::get_curvature_epsilon()
313  {
314  return this->curvature_epsilon;
315  }
316 
317  void Linearizer::process_quad(MeshFunction<double>** fns, int iv0, int iv1, int iv2, int iv3, int level,
318  double* val, double* phx, double* phy, int* idx, bool curved)
319  {
320  double midval[3][5];
321 
322  // try not to split through the vertex with the largest value
323  int a = (verts[iv0][2] > verts[iv1][2]) ? iv0 : iv1;
324  int b = (verts[iv2][2] > verts[iv3][2]) ? iv2 : iv3;
325  a = (verts[a][2] > verts[b][2]) ? a : b;
326  int flip = (a == iv1 || a == iv3) ? 1 : 0;
327 
328  if(level < LIN_MAX_LEVEL)
329  {
330  int i;
331  if(!(level & 1)) // this is an optimization: do the following only every other time
332  {
333  // obtain solution values
334  fns[0]->set_quad_order(1, item);
335  val = fns[0]->get_values(component, value_type);
336  if(auto_max)
337  for (i = 0; i < lin_np_quad[1]; i++)
338  {
339  double v = val[i];
340  if(finite(v) && fabs(v) > max)
341 #pragma omp critical(max)
342  if(finite(v) && fabs(v) > max)
343  max = fabs(v);
344  }
345 
346  // This is just to make some sense.
347  if(fabs(max) < 1E-10)
348  max = 1E-10;
349 
350  idx = quad_indices[0];
351 
352  if(curved)
353  {
354  RefMap* refmap = fns[0]->get_refmap();
355  phx = refmap->get_phys_x(1);
356  phy = refmap->get_phys_y(1);
357 
358  double* dx = NULL;
359  double* dy = NULL;
360 
361  if(this->xdisp != NULL)
362  fns[1]->set_quad_order(1, H2D_FN_VAL);
363  if(this->ydisp != NULL)
364  fns[2]->set_quad_order(1, H2D_FN_VAL);
365  if(this->xdisp != NULL)
366  dx = fns[1]->get_fn_values();
367  if(this->ydisp != NULL)
368  dy = fns[2]->get_fn_values();
369  for (i = 0; i < lin_np_quad[1]; i++)
370  {
371  if(this->xdisp != NULL)
372  phx[i] += dmult*dx[i];
373  if(this->ydisp != NULL)
374  phy[i] += dmult*dy[i];
375  }
376  }
377  }
378 
379  // obtain linearized values and coordinates at the midpoints
380  for (i = 0; i < 3; i++)
381  {
382  midval[i][0] = (verts[iv0][i] + verts[iv1][i]) * 0.5;
383  midval[i][1] = (verts[iv1][i] + verts[iv2][i]) * 0.5;
384  midval[i][2] = (verts[iv2][i] + verts[iv3][i]) * 0.5;
385  midval[i][3] = (verts[iv3][i] + verts[iv0][i]) * 0.5;
386  midval[i][4] = (midval[i][0] + midval[i][2]) * 0.5;
387  };
388 
389  // the value of the middle point is not the average of the four vertex values, since quad == 2 triangles
390  midval[2][4] = flip ? (verts[iv0][2] + verts[iv2][2]) * 0.5 : (verts[iv1][2] + verts[iv3][2]) * 0.5;
391 
392  // determine whether or not to split the element
393  int split;
394  if(eps >= 1.0)
395  {
396  // if eps > 1, the user wants a fixed number of refinements (no adaptivity)
397  split = (level < eps) ? 3 : 0;
398  }
399  else
400  {
401  if(!auto_max && fabs(verts[iv0][2]) > max && fabs(verts[iv1][2]) > max
402  && fabs(verts[iv2][2]) > max && fabs(verts[iv3][2]) > max)
403  {
404  // do not split if the whole quad is above the specified maximum value
405  split = 0;
406  }
407  else
408  {
409  // calculate the approximate error of linearizing the normalized solution
410  double herr = fabs(val[idx[1]] - midval[2][1]) + fabs(val[idx[3]] - midval[2][3]);
411  double verr = fabs(val[idx[0]] - midval[2][0]) + fabs(val[idx[2]] - midval[2][2]);
412  double err = fabs(val[idx[4]] - midval[2][4]) + herr + verr;
413  split = (!finite(err) || err > max*4*eps) ? 3 : 0;
414 
415  // decide whether to split horizontally or vertically only
416  if(level > 0 && split)
417  {
418  if(herr > 5*verr)
419  split = 1; // h-split
420  else if(verr > 5*herr)
421  split = 2; // v-split
422  }
423  }
424 
425  // also decide whether to split because of the curvature
426  if(split != 3 && curved)
427  {
428  double cm2 = sqr(fns[0]->get_active_element()->get_diameter()*this->get_curvature_epsilon());
429  if(sqr(phx[idx[1]] - midval[0][1]) + sqr(phy[idx[1]] - midval[1][1]) > cm2 ||
430  sqr(phx[idx[3]] - midval[0][3]) + sqr(phy[idx[3]] - midval[1][3]) > cm2) split |= 1;
431  if(sqr(phx[idx[0]] - midval[0][0]) + sqr(phy[idx[0]] - midval[1][0]) > cm2 ||
432  sqr(phx[idx[2]] - midval[0][2]) + sqr(phy[idx[2]] - midval[1][2]) > cm2) split |= 2;
433  }
434 
435  // do extra tests at level 0, so as not to miss some functions with zero error at edge midpoints
436  if(level == 0 && !split)
437  {
438  split = ((fabs(val[13] - 0.5*(midval[2][0] + midval[2][1])) +
439  fabs(val[17] - 0.5*(midval[2][1] + midval[2][2])) +
440  fabs(val[20] - 0.5*(midval[2][2] + midval[2][3])) +
441  fabs(val[9] - 0.5*(midval[2][3] + midval[2][0]))) > max*4*eps) ? 3 : 0;
442  }
443  }
444 
445  // split the quad if the error is too large, otherwise produce two linear triangles
446  if(split)
447  {
448  if(curved)
449  for (i = 0; i < 5; i++)
450  {
451  midval[0][i] = phx[idx[i]];
452  midval[1][i] = phy[idx[i]];
453  }
454 
455  // obtain mid-edge and mid-element vertices
456  int mid0, mid1, mid2, mid3, mid4;
457  if(split != 1) mid0 = get_vertex(iv0, iv1, midval[0][0], midval[1][0], val[idx[0]]);
458  if(split != 2) mid1 = get_vertex(iv1, iv2, midval[0][1], midval[1][1], val[idx[1]]);
459  if(split != 1) mid2 = get_vertex(iv2, iv3, midval[0][2], midval[1][2], val[idx[2]]);
460  if(split != 2) mid3 = get_vertex(iv3, iv0, midval[0][3], midval[1][3], val[idx[3]]);
461  if(split == 3) mid4 = get_vertex(mid0, mid2, midval[0][4], midval[1][4], val[idx[4]]);
462 
463  // recur to sub-elements
464  if(split == 3)
465  {
466  fns[0]->push_transform(0);
467  if(this->xdisp != NULL)
468  if(fns[1] != fns[0])
469  fns[1]->push_transform(0);
470  if(this->ydisp != NULL)
471  if(fns[2] != fns[1])
472  fns[2]->push_transform(0);
473  process_quad(fns, iv0, mid0, mid4, mid3, level + 1, val, phx, phy, quad_indices[1], curved);
474  fns[0]->pop_transform();
475  if(this->xdisp != NULL)
476  if(fns[1] != fns[0])
477  fns[1]->pop_transform();
478  if(this->ydisp != NULL)
479  if(fns[2] != fns[1])
480  fns[2]->pop_transform();
481 
482  fns[0]->push_transform(1);
483  if(this->xdisp != NULL)
484  if(fns[1] != fns[0])
485  fns[1]->push_transform(1);
486  if(this->ydisp != NULL)
487  if(fns[2] != fns[1])
488  fns[2]->push_transform(1);
489  process_quad(fns, mid0, iv1, mid1, mid4, level + 1, val, phx, phy, quad_indices[2], curved);
490  fns[0]->pop_transform();
491  if(this->xdisp != NULL)
492  if(fns[1] != fns[0])
493  fns[1]->pop_transform();
494  if(this->ydisp != NULL)
495  if(fns[2] != fns[1])
496  fns[2]->pop_transform();
497 
498  fns[0]->push_transform(2);
499  if(this->xdisp != NULL)
500  if(fns[1] != fns[0])
501  fns[1]->push_transform(2);
502  if(this->ydisp != NULL)
503  if(fns[2] != fns[1])
504  fns[2]->push_transform(2);
505  process_quad(fns, mid4, mid1, iv2, mid2, level + 1, val, phx, phy, quad_indices[3], curved);
506  fns[0]->pop_transform();
507  if(this->xdisp != NULL)
508  if(fns[1] != fns[0])
509  fns[1]->pop_transform();
510  if(this->ydisp != NULL)
511  if(fns[2] != fns[1])
512  fns[2]->pop_transform();
513 
514  fns[0]->push_transform(3);
515  if(this->xdisp != NULL)
516  if(fns[1] != fns[0])
517  fns[1]->push_transform(3);
518  if(this->ydisp != NULL)
519  if(fns[2] != fns[1])
520  fns[2]->push_transform(3);
521  process_quad(fns, mid3, mid4, mid2, iv3, level + 1, val, phx, phy, quad_indices[4], curved);
522  fns[0]->pop_transform();
523  if(this->xdisp != NULL)
524  if(fns[1] != fns[0])
525  fns[1]->pop_transform();
526  if(this->ydisp != NULL)
527  if(fns[2] != fns[1])
528  fns[2]->pop_transform();
529  }
530  else
531  if(split == 1) // h-split
532  {
533  fns[0]->push_transform(4);
534  if(this->ydisp != NULL)
535  if(fns[1] != fns[0])
536  fns[1]->push_transform(4);
537  if(this->ydisp != NULL)
538  if(fns[2] != fns[1])
539  fns[2]->push_transform(4);
540  process_quad(fns, iv0, iv1, mid1, mid3, level + 1, val, phx, phy, quad_indices[5], curved);
541  fns[0]->pop_transform();
542  if(this->xdisp != NULL)
543  if(fns[1] != fns[0])
544  fns[1]->pop_transform();
545  if(this->ydisp != NULL)
546  if(fns[2] != fns[1])
547  fns[2]->pop_transform();
548 
549  fns[0]->push_transform(5);
550  if(this->xdisp != NULL)
551  if(fns[1] != fns[0])
552  fns[1]->push_transform(5);
553  if(this->ydisp != NULL)
554  if(fns[2] != fns[1])
555  fns[2]->push_transform(5);
556  process_quad(fns, mid3, mid1, iv2, iv3, level + 1, val, phx, phy, quad_indices[6], curved);
557  fns[0]->pop_transform();
558  if(this->xdisp != NULL)
559  if(fns[1] != fns[0])
560  fns[1]->pop_transform();
561  if(this->ydisp != NULL)
562  if(fns[2] != fns[1])
563  fns[2]->pop_transform();
564  }
565  else // v-split
566  {
567  fns[0]->push_transform(6);
568  if(this->xdisp != NULL)
569  if(fns[1] != fns[0])
570  fns[1]->push_transform(6);
571  if(this->ydisp != NULL)
572  if(fns[2] != fns[1])
573  fns[2]->push_transform(6);
574  process_quad(fns, iv0, mid0, mid2, iv3, level + 1, val, phx, phy, quad_indices[7], curved);
575  fns[0]->pop_transform();
576  if(this->xdisp != NULL)
577  if(fns[1] != fns[0])
578  fns[1]->pop_transform();
579  if(this->ydisp != NULL)
580  if(fns[2] != fns[1])
581  fns[2]->pop_transform();
582 
583  fns[0]->push_transform(7);
584  if(this->xdisp != NULL)
585  if(fns[1] != fns[0])
586  fns[1]->push_transform(7);
587  if(this->ydisp != NULL)
588  if(fns[2] != fns[1])
589  fns[2]->push_transform(7);
590  process_quad(fns, mid0, iv1, iv2, mid2, level + 1, val, phx, phy, quad_indices[8], curved);
591  fns[0]->pop_transform();
592  if(this->xdisp != NULL)
593  if(fns[1] != fns[0])
594  fns[1]->pop_transform();
595  if(this->ydisp != NULL)
596  if(fns[2] != fns[1])
597  fns[2]->pop_transform();
598  }
599  return;
600  }
601  }
602 
603  // output two linear triangles,
604  if(!flip)
605  {
606  add_triangle(iv3, iv0, iv1, fns[0]->get_active_element()->marker);
607  add_triangle(iv1, iv2, iv3, fns[0]->get_active_element()->marker);
608  }
609  else
610  {
611  add_triangle(iv0, iv1, iv2, fns[0]->get_active_element()->marker);
612  add_triangle(iv2, iv3, iv0, fns[0]->get_active_element()->marker);
613  }
614  }
615 
616  void Linearizer::set_displacement(MeshFunction<double>* xdisp, MeshFunction<double>* ydisp, double dmult)
617  {
618  if(xdisp != NULL)
619  {
620  user_xdisp = true;
621  this->xdisp = xdisp;
622  }
623  if(ydisp != NULL)
624  {
625  user_ydisp = true;
626  this->ydisp = ydisp;
627  }
628  this->dmult = dmult;
629  }
630 
631  void Linearizer::process_solution(MeshFunction<double>* sln, int item_, double eps)
632  {
633  // Important, sets the current caughtException to NULL.
634  this->caughtException = NULL;
635 
636  lock_data();
637  this->tick();
638 
639  // Initialization of 'global' stuff.
640  this->item = item_;
641  this->eps = eps;
642  // get the component and desired value from item.
643  if(item >= 0x40)
644  {
645  component = 1;
646  this->item >>= 6;
647  }
648  while (!(item & 1))
649  {
650  this->item >>= 1;
651  value_type++;
652  }
653  // reset the item to the value before the circus with component, value_type.
654  this->item = item_;
655 
656  // Initialization of computation stuff.
657  // sizes.
658  this->vertex_size = std::max(100 * sln->get_mesh()->get_num_elements(), std::max(this->vertex_size, 50000));
659  this->triangle_size = std::max(150 * sln->get_mesh()->get_num_elements(), std::max(this->triangle_size, 75000));
660  this->edges_size = std::max(100 * sln->get_mesh()->get_num_elements(), std::max(this->edges_size, 50000));
661  // counts.
662  this->vertex_count = 0;
663  this->triangle_count = 0;
664  this->edges_count = 0;
665  // reuse or allocate vertex, triangle and edge arrays.
666  this->verts = (double3*) realloc(this->verts, sizeof(double3) * this->vertex_size);
667  this->tris = (int3*) realloc(this->tris, sizeof(int3) * this->triangle_size);
668  this->tri_markers = (int*) realloc(this->tri_markers, sizeof(int) * this->triangle_size);
669  this->edges = (int2*) realloc(this->edges, sizeof(int2) * this->edges_size);
670  this->edge_markers = (int*) realloc(this->edge_markers, sizeof(int) * this->edges_size);
671  this->info = (int4*) malloc(sizeof(int4) * this->vertex_size);
672  this->empty = false;
673  // initialize the hash table
674  this->hash_table = (int*) malloc(sizeof(int) * this->vertex_size);
675  memset(this->hash_table, 0xff, sizeof(int) * this->vertex_size);
676 
677  // select the linearization quadratures
678  Quad2D *old_quad, *old_quad_x = NULL, *old_quad_y = NULL;
679  old_quad = sln->get_quad_2d();
680  if(xdisp != NULL)
681  old_quad_x = xdisp->get_quad_2d();
682  if(ydisp != NULL)
683  old_quad_y = ydisp->get_quad_2d();
684 
685  // obtain the solution in vertices, estimate the maximum solution value
686  // meshes.
687  Hermes::vector<const Mesh*> meshes;
688  meshes.push_back(sln->get_mesh());
689  if(xdisp != NULL)
690  meshes.push_back(xdisp->get_mesh());
691  if(ydisp != NULL)
692  meshes.push_back(ydisp->get_mesh());
693 
694  // Parallelization
695  MeshFunction<double>*** fns = new MeshFunction<double>**[Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
696  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
697  {
698  fns[i] = new MeshFunction<double>*[3];
699  fns[i][0] = sln->clone();
700  fns[i][0]->set_refmap(new RefMap);
701  fns[i][0]->set_quad_2d(&g_quad_lin);
702  if(xdisp != NULL)
703  {
704  fns[i][1] = xdisp->clone();
705  //fns[i][1]->set_refmap(new RefMap);
706  fns[i][1]->set_quad_2d(&g_quad_lin);
707  }
708  if(ydisp != NULL)
709  {
710  fns[i][2] = ydisp->clone();
711  //fns[i][2]->set_refmap(new RefMap);
712  fns[i][2]->set_quad_2d(&g_quad_lin);
713  }
714  }
715 
716  Transformable*** trfs = new Transformable**[Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
717  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
718  {
719  trfs[i] = new Transformable*[3];
720  trfs[i][0] = fns[i][0];
721  if(xdisp != NULL)
722  trfs[i][1] = fns[i][1];
723  if(ydisp != NULL)
724  trfs[i][2] = fns[i][2];
725  }
726 
727  Traverse trav_masterMax(true);
728  unsigned int num_states = trav_masterMax.get_num_states(meshes);
729 
730  trav_masterMax.begin(meshes.size(), &(meshes.front()));
731 
732  Traverse* trav = new Traverse[Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
733 
734  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
735  {
736  trav[i].begin(meshes.size(), &(meshes.front()), trfs[i]);
737  trav[i].stack = trav_masterMax.stack;
738  }
739 
740  int state_i;
741 
742 #define CHUNKSIZE 1
743  int num_threads_used = Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads);
744 #pragma omp parallel shared(trav_masterMax) private(state_i) num_threads(num_threads_used)
745  {
746 #pragma omp for schedule(dynamic, CHUNKSIZE)
747  for(state_i = 0; state_i < num_states; state_i++)
748  {
749  try
750  {
751  Traverse::State current_state;
752 #pragma omp critical(get_next_state)
753  current_state = trav[omp_get_thread_num()].get_next_state(&trav_masterMax.top, &trav_masterMax.id);
754 
755  fns[omp_get_thread_num()][0]->set_quad_order(0, this->item);
756  double* val = fns[omp_get_thread_num()][0]->get_values(component, value_type);
757 
758  for (unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
759  {
760  double f = val[i];
761 #pragma omp critical (max)
762  if(this->auto_max && finite(f) && fabs(f) > this->max)
763  this->max = fabs(f);
764  }
765  }
766  catch(Hermes::Exceptions::Exception& e)
767  {
768  if(this->caughtException == NULL)
769  this->caughtException = e.clone();
770  }
771  catch(std::exception& e)
772  {
773  if(this->caughtException == NULL)
774  this->caughtException = new Hermes::Exceptions::Exception(e.what());
775  }
776  }
777  }
778 
779  trav_masterMax.finish();
780  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
781  trav[i].finish();
782  delete [] trav;
783 
784  Traverse trav_master(true);
785  num_states = trav_master.get_num_states(meshes);
786 
787  trav_master.begin(meshes.size(), &(meshes.front()));
788 
789  trav = new Traverse[Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
790 
791  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
792  {
793  trav[i].begin(meshes.size(), &(meshes.front()), trfs[i]);
794  trav[i].stack = trav_master.stack;
795  }
796 
797 #pragma omp parallel shared(trav_master) private(state_i) num_threads(num_threads_used)
798  {
799 #pragma omp for schedule(dynamic, CHUNKSIZE)
800  for(state_i = 0; state_i < num_states; state_i++)
801  {
802  try
803  {
804  Traverse::State current_state;
805 
806 #pragma omp critical (get_next_state)
807  current_state = trav[omp_get_thread_num()].get_next_state(&trav_master.top, &trav_master.id);
808 
809  fns[omp_get_thread_num()][0]->set_quad_order(0, this->item);
810  double* val = fns[omp_get_thread_num()][0]->get_values(component, value_type);
811  if(val == NULL)
812  {
813  delete [] trav;
814  throw Hermes::Exceptions::Exception("Item not defined in the solution.");
815  }
816 
817  if(xdisp != NULL)
818  fns[omp_get_thread_num()][1]->set_quad_order(0, H2D_FN_VAL);
819  if(ydisp != NULL)
820  fns[omp_get_thread_num()][2]->set_quad_order(0, H2D_FN_VAL);
821 
822  double *dx = NULL;
823  double *dy = NULL;
824  if(xdisp != NULL)
825  dx = fns[omp_get_thread_num()][1]->get_fn_values();
826  if(ydisp != NULL)
827  dy = fns[omp_get_thread_num()][2]->get_fn_values();
828 
829  int iv[H2D_MAX_NUMBER_VERTICES];
830  for (unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
831  {
832  double f = val[i];
833  double x_disp = fns[omp_get_thread_num()][0]->get_refmap()->get_phys_x(0)[i];
834  double y_disp = fns[omp_get_thread_num()][0]->get_refmap()->get_phys_y(0)[i];
835  if(this->xdisp != NULL)
836  x_disp += dmult * dx[i];
837  if(this->ydisp != NULL)
838  y_disp += dmult * dy[i];
839 
840  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);
841  }
842 
843  // recur to sub-elements
844  if(current_state.e[0]->is_triangle())
845  process_triangle(fns[omp_get_thread_num()], iv[0], iv[1], iv[2], 0, NULL, NULL, NULL, NULL, current_state.e[0]->is_curved());
846  else
847  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());
848 
849  for (unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
850  process_edge(iv[i], iv[current_state.e[0]->next_vert(i)], current_state.e[0]->en[i]->marker);
851  }
852  catch(Hermes::Exceptions::Exception& e)
853  {
854  if(this->caughtException == NULL)
855  this->caughtException = e.clone();
856  }
857  catch(std::exception& e)
858  {
859  if(this->caughtException == NULL)
860  this->caughtException = new Hermes::Exceptions::Exception(e.what());
861  }
862  }
863  }
864 
865  trav_master.finish();
866  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
867  {
868  trav[i].finish();
869  for(unsigned int j = 0; j < (1 + (xdisp != NULL? 1 : 0) + (ydisp != NULL ? 1 : 0)); j++)
870  delete fns[i][j];
871  delete [] fns[i];
872  delete [] trfs[i];
873  }
874  delete [] fns;
875  delete [] trfs;
876  delete [] trav;
877 
878  // for contours, without regularization.
879  this->tris_contours = (int3*) realloc(this->tris_contours, sizeof(int3) * this->triangle_count);
880  memcpy(this->tris_contours, this->tris, this->triangle_count * sizeof(int3));
881  triangle_contours_count = this->triangle_count;
882 
883  if(this->caughtException != NULL)
884  {
885  this->unlock_data();
886  ::free(hash_table);
887  ::free(info);
888  throw *(this->caughtException);
889  }
890 
891  // regularize the linear mesh
892  for (int i = 0; i < this->triangle_count; i++)
893  {
894  int iv0 = tris[i][0], iv1 = tris[i][1], iv2 = tris[i][2];
895 
896  int mid0 = peek_vertex(iv0, iv1);
897  int mid1 = peek_vertex(iv1, iv2);
898  int mid2 = peek_vertex(iv2, iv0);
899  if(mid0 >= 0 || mid1 >= 0 || mid2 >= 0)
900  {
901  this->del_slot = i;
902  regularize_triangle(iv0, iv1, iv2, mid0, mid1, mid2, tri_markers[i]);
903  }
904  }
905 
906  find_min_max();
907 
908  this->unlock_data();
909 
910  // select old quadratrues
911  sln->set_quad_2d(old_quad);
912  if(user_xdisp)
913  xdisp->set_quad_2d(old_quad_x);
914  else
915  delete xdisp;
916  if(user_ydisp)
917  ydisp->set_quad_2d(old_quad_y);
918  else
919  delete ydisp;
920 
921  // clean up
922  ::free(hash_table);
923  ::free(info);
924  }
925 
926  void Linearizer::find_min_max()
927  {
928  // find min & max vertex values
929  this->min_val = 1e100;
930  this->max_val = -1e100;
931  for (int i = 0; i < this->vertex_count; i++)
932  {
933  if(finite(verts[i][2]) && verts[i][2] < min_val) min_val = verts[i][2];
934  if(finite(verts[i][2]) && verts[i][2] > max_val) max_val = verts[i][2];
935  }
936  }
937 
938  int Linearizer::get_vertex(int p1, int p2, double x, double y, double value)
939  {
940  // search for an existing vertex
941  if(p1 > p2) std::swap(p1, p2);
942  int index = this->hash(p1, p2);
943  int i = 0;
944  if(index < this->vertex_count)
945  {
946  i = this->hash_table[index];
947  while (i >= 0 && i < this->vertex_count)
948  {
949  if(
950  this->info[i][0] == p1 && this->info[i][1] == p2 &&
951  (value == verts[i][2] || fabs(value - verts[i][2]) < this->max*1e-8) &&
952  (fabs(x - verts[i][0]) < 1e-8) &&
953  (fabs(y - verts[i][1]) < 1e-8)
954  )
955  return i;
956  // note that we won't return a vertex with a different value than the required one;
957  // this takes care for discontinuities in the solution, where more vertices
958  // with different values will be created
959  i = info[i][2];
960  }
961  }
962 
963  // if not found, create a new one
964 #pragma omp critical(realloc_vertices)
965  i = add_vertex();
966  verts[i][0] = x;
967  verts[i][1] = y;
968  verts[i][2] = value;
969  this->info[i][0] = p1;
970  this->info[i][1] = p2;
971  this->info[i][2] = hash_table[index];
972  this->hash_table[index] = i;
973  return i;
974  }
975 
976  int Linearizer::add_vertex()
977  {
978  if(this->vertex_count >= this->vertex_size)
979  {
980  this->vertex_size *= 2;
981  verts = (double3*) realloc(verts, sizeof(double3) * vertex_size);
982  this->info = (int4*) realloc(info, sizeof(int4) * vertex_size);
983  this->hash_table = (int*) realloc(hash_table, sizeof(int) * vertex_size);
984  memset(this->hash_table + this->vertex_size / 2, 0xff, sizeof(int) * this->vertex_size / 2);
985  }
986  return this->vertex_count++;
987  }
988 
989  void Linearizer::free()
990  {
991  if(verts != NULL)
992  {
993  ::free(verts);
994  verts = NULL;
995  }
996  if(tris_contours != NULL)
997  {
998  ::free(tris_contours);
999  tris_contours = NULL;
1000  }
1001 
1002  LinearizerBase::free();
1003  }
1004 
1005  Linearizer::~Linearizer()
1006  {
1007  free();
1008  }
1009 
1010  void Linearizer::save_solution_vtk(MeshFunction<double>* sln, const char* filename, const char *quantity_name,
1011  bool mode_3D, int item, double eps)
1012  {
1013  process_solution(sln, item, eps);
1014 
1015  FILE* f = fopen(filename, "wb");
1016  if(f == NULL) throw Hermes::Exceptions::Exception("Could not open %s for writing.", filename);
1017  lock_data();
1018 
1019  // Output header for vertices.
1020  fprintf(f, "# vtk DataFile Version 2.0\n");
1021  fprintf(f, "\n");
1022  fprintf(f, "ASCII\n\n");
1023  fprintf(f, "DATASET UNSTRUCTURED_GRID\n");
1024 
1025  // Output vertices.
1026  fprintf(f, "POINTS %d %s\n", this->vertex_count, "float");
1027  for (int i = 0; i < this->vertex_count; i++)
1028  {
1029  if(mode_3D == true) fprintf(f, "%g %g %g\n", this->verts[i][0], this->verts[i][1], this->verts[i][2]);
1030  else fprintf(f, "%g %g %g\n", this->verts[i][0], this->verts[i][1], 0.0);
1031  }
1032 
1033  // Output elements.
1034  fprintf(f, "\n");
1035  fprintf(f, "CELLS %d %d\n", this->triangle_count, 4 * this->triangle_count);
1036  for (int i = 0; i < this->triangle_count; i++)
1037  {
1038  fprintf(f, "3 %d %d %d\n", this->tris[i][0], this->tris[i][1], this->tris[i][2]);
1039  }
1040 
1041  // Output cell types.
1042  fprintf(f, "\n");
1043  fprintf(f, "CELL_TYPES %d\n", this->triangle_count);
1044  for (int i = 0; i < this->triangle_count; i++)
1045  {
1046  fprintf(f, "5\n"); // The "5" means triangle in VTK.
1047  }
1048 
1049  // This outputs double solution values.
1050  fprintf(f, "\n");
1051  fprintf(f, "POINT_DATA %d\n", this->vertex_count);
1052  fprintf(f, "SCALARS %s %s %d\n", quantity_name, "float", 1);
1053  fprintf(f, "LOOKUP_TABLE %s\n", "default");
1054  for (int i = 0; i < this->vertex_count; i++)
1055  {
1056  fprintf(f, "%g\n", this->verts[i][2]);
1057  }
1058 
1059  unlock_data();
1060  fclose(f);
1061  }
1062 
1063  void Linearizer::calc_vertices_aabb(double* min_x, double* max_x, double* min_y, double* max_y) const
1064  {
1065  if(verts == NULL)
1066  throw Exceptions::Exception("Cannot calculate AABB from NULL vertices");
1067  calc_aabb(&verts[0][0], &verts[0][1], sizeof(double3), vertex_count, min_x, max_x, min_y, max_y);
1068  }
1069 
1070  double3* Linearizer::get_vertices()
1071  {
1072  return this->verts;
1073  }
1074  int Linearizer::get_num_vertices()
1075  {
1076  return this->vertex_count;
1077  }
1078 
1079  int Linearizer::get_num_contour_triangles()
1080  {
1081  return this->triangle_contours_count;
1082  }
1083 
1084  int3* Linearizer::get_contour_triangles()
1085  {
1086  return this->tris_contours;
1087  }
1088  }
1089  }
1090 }