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