Hermes2D  2.0
vectorizerTest.h
1 void Vectorizer::process_triangle(MeshFunction<double>** fns, int iv0, int iv1, int iv2, int level,
2  double* xval, double* yval, double* phx, double* phy, int* idx, bool curved)
3  {
4  double midval[4][3];
5 
6  if(level < LIN_MAX_LEVEL)
7  {
8  int i;
9  if(!(level & 1))
10  {
11  // obtain solution values and physical element coordinates
12  fns[0]->set_quad_order(1, xitem);
13  fns[1]->set_quad_order(1, yitem);
14  xval = fns[0]->get_values(component_x, value_type_x);
15  yval = fns[1]->get_values(component_y, value_type_y);
16  for (i = 0; i < lin_np_tri[1]; i++)
17  {
18  double m = (sqrt(sqr(xval[i]) + sqr(yval[i])));
19 #pragma omp critical(max)
20  if(finite(m) && fabs(m) > max)
21  max = fabs(m);
22  }
23 
24  idx = tri_indices[0];
25 
26  if(curved)
27  {
28  RefMap* refmap = fns[0]->get_refmap();
29  phx = refmap->get_phys_x(1);
30  phy = refmap->get_phys_y(1);
31  }
32  }
33 
34  // obtain linearized values and coordinates at the midpoints
35  for (i = 0; i < 4; i++)
36  {
37  midval[i][0] = (verts[iv0][i] + verts[iv1][i])*0.5;
38  midval[i][1] = (verts[iv1][i] + verts[iv2][i])*0.5;
39  midval[i][2] = (verts[iv2][i] + verts[iv0][i])*0.5;
40  };
41 
42  // determine whether or not to split the element
43  bool split;
44  if(eps >= 1.0)
45  {
46  //if eps > 1, the user wants a fixed number of refinements (no adaptivity)
47  split = (level < eps);
48  }
49  else
50  {
51  // calculate the approximate error of linearizing the normalized solution
52  double err = fabs(sqrt(sqr(xval[idx[0]]) + sqr(yval[idx[0]])) - sqrt(sqr(midval[2][0]) + sqr(midval[3][0]))) +
53  fabs(sqrt(sqr(xval[idx[1]]) + sqr(yval[idx[1]])) - sqrt(sqr(midval[2][1]) + sqr(midval[3][1]))) +
54  fabs(sqrt(sqr(xval[idx[2]]) + sqr(yval[idx[2]])) - sqrt(sqr(midval[2][2]) + sqr(midval[3][2])));
55  split = !finite(err) || err > max*3*eps;
56 
57  // do the same for the curvature
58  if(curved)
59  {
60  for (i = 0; i < 4; i++)
61  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()))
62  {
63  split = true;
64  break;
65  }
66  }
67 
68  // do extra tests at level 0, so as not to miss some functions with zero error at edge midpoints
69  if(level == 0 && !split)
70  {
71  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])))) +
72  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])))) +
73  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;
74  }
75  }
76 
77  // split the triangle if the error is too large, otherwise produce a linear triangle
78  if(split)
79  {
80  if(curved)
81  for (i = 0; i < 3; i++)
82  {
83  midval[0][i] = phx[idx[i]];
84  midval[1][i] = phy[idx[i]];
85  }
86 
87  // obtain mid-edge vertices
88  int mid0 = get_vertex(iv0, iv1, midval[0][0], midval[1][0], xval[idx[0]], yval[idx[0]]);
89  int mid1 = get_vertex(iv1, iv2, midval[0][1], midval[1][1], xval[idx[1]], yval[idx[1]]);
90  int mid2 = get_vertex(iv2, iv0, midval[0][2], midval[1][2], xval[idx[2]], yval[idx[2]]);
91 
92  // recur to sub-elements
93  fns[0]->push_transform(0); if(fns[1] != fns[0]) fns[1]->push_transform(0);
94  process_triangle(fns, iv0, mid0, mid2, level + 1, xval, yval, phx, phy, tri_indices[1], curved);
95  fns[0]->pop_transform(); if(fns[1] != fns[0]) fns[1]->pop_transform();
96 
97  fns[0]->push_transform(1); if(fns[1] != fns[0]) fns[1]->push_transform(1);
98  process_triangle(fns, mid0, iv1, mid1, level + 1, xval, yval, phx, phy, tri_indices[2], curved);
99  fns[0]->pop_transform(); if(fns[1] != fns[0]) fns[1]->pop_transform();
100 
101  fns[0]->push_transform(2); if(fns[1] != fns[0]) fns[1]->push_transform(2);
102  process_triangle(fns, mid2, mid1, iv2, level + 1, xval, yval, phx, phy, tri_indices[3], curved);
103  fns[0]->pop_transform(); if(fns[1] != fns[0]) fns[1]->pop_transform();
104 
105  fns[0]->push_transform(3); if(fns[1] != fns[0]) fns[1]->push_transform(3);
106  process_triangle(fns, mid1, mid2, mid0, level + 1, xval, yval, phx, phy, tri_indices[4], curved);
107  fns[0]->pop_transform(); if(fns[1] != fns[0]) fns[1]->pop_transform();
108  return;
109  }
110  }
111 
112  // no splitting: output a linear triangle
113  add_triangle(iv0, iv1, iv2, fns[0]->get_active_element()->marker);
114  }
115 
116  void Vectorizer::process_quad(MeshFunction<double>** fns, int iv0, int iv1, int iv2, int iv3, int level,
117  double* xval, double* yval, double* phx, double* phy, int* idx, bool curved)
118  {
119  double midval[4][5];
120 
121  // try not to split through the vertex with the largest value
122  int a = (sqr(verts[iv1][2]) + sqr(verts[iv1][3]) > sqr(verts[iv0][2]) + sqr(verts[iv0][3])) ? iv0 : iv1;
123  int b = (sqr(verts[iv2][2]) + sqr(verts[iv2][3]) > sqr(verts[iv3][2]) + sqr(verts[iv3][3])) ? iv2 : iv3;
124  a = (sqr(verts[a][2]) + sqr(verts[a][3]) > sqr(verts[b][2]) + sqr(verts[b][3])) ? a : b;
125  int flip = (a == iv1 || a == iv3) ? 1 : 0;
126 
127  if(level < LIN_MAX_LEVEL)
128  {
129  int i;
130  if(!(level & 1))
131  {
132  // obtain solution values and physical element coordinates
133  fns[0]->set_quad_order(1, xitem);
134  fns[1]->set_quad_order(1, yitem);
135  xval = fns[0]->get_values(component_x, value_type_x);
136  yval = fns[1]->get_values(component_y, value_type_y);
137  for (i = 0; i < lin_np_quad[1]; i++)
138  {
139  double m = sqrt(sqr(xval[i]) + sqr(yval[i]));
140  if(finite(m) && fabs(m) > max)
141 #pragma omp critical(max)
142  if(finite(m) && fabs(m) > max)
143  max = fabs(m);
144  }
145 
146  // This is just to make some sense.
147  if(fabs(max) < 1E-10)
148  max = 1E-10;
149 
150  idx = quad_indices[0];
151 
152  if(curved)
153  {
154  RefMap* refmap = fns[0]->get_refmap();
155  phx = refmap->get_phys_x(1);
156  phy = refmap->get_phys_y(1);
157  }
158  }
159 
160  // obtain linearized values and coordinates at the midpoints
161  for (i = 0; i < 4; i++)
162  {
163  midval[i][0] = (verts[iv0][i] + verts[iv1][i]) * 0.5;
164  midval[i][1] = (verts[iv1][i] + verts[iv2][i]) * 0.5;
165  midval[i][2] = (verts[iv2][i] + verts[iv3][i]) * 0.5;
166  midval[i][3] = (verts[iv3][i] + verts[iv0][i]) * 0.5;
167  midval[i][4] = (midval[i][0] + midval[i][2]) * 0.5;
168  };
169 
170  // determine whether or not to split the element
171  bool split;
172 
173  //if eps > 1, the user wants a fixed number of refinements (no adaptivity)
174  if(eps >= 1.0)
175  split = (level < eps);
176  else
177  {
178  // the value of the middle point is not the average of the four vertex values, since quad == 2 triangles
179  midval[2][4] = flip ? (verts[iv0][2] + verts[iv2][2]) * 0.5 : (verts[iv1][2] + verts[iv3][2]) * 0.5;
180  midval[3][4] = flip ? (verts[iv0][3] + verts[iv2][3]) * 0.5 : (verts[iv1][3] + verts[iv3][3]) * 0.5;
181 
182  // calculate the approximate error of linearizing the normalized solution
183  double err = fabs(sqrt(sqr(xval[idx[0]]) + sqr(yval[idx[0]])) - sqrt(sqr(midval[2][0]) + sqr(midval[3][0]))) +
184  fabs(sqrt(sqr(xval[idx[1]]) + sqr(yval[idx[1]])) - sqrt(sqr(midval[2][1]) + sqr(midval[3][1]))) +
185  fabs(sqrt(sqr(xval[idx[2]]) + sqr(yval[idx[2]])) - sqrt(sqr(midval[2][2]) + sqr(midval[3][2]))) +
186  fabs(sqrt(sqr(xval[idx[3]]) + sqr(yval[idx[3]])) - sqrt(sqr(midval[2][3]) + sqr(midval[3][3]))) +
187  fabs(sqrt(sqr(xval[idx[4]]) + sqr(yval[idx[4]])) - sqrt(sqr(midval[2][4]) + sqr(midval[3][4])));
188  split = !finite(err) || err > max*40*eps;
189 
190  // do the same for the curvature
191  if(curved && !split)
192  {
193  double cm2 = sqr(fns[0]->get_active_element()->get_diameter()*this->get_curvature_epsilon());
194  if(sqr(phx[idx[1]] - midval[0][1]) + sqr(phy[idx[1]] - midval[1][1]) > cm2 ||
195  sqr(phx[idx[3]] - midval[0][3]) + sqr(phy[idx[3]] - midval[1][3]) > cm2) split = true;
196 
197  if(sqr(phx[idx[0]] - midval[0][0]) + sqr(phy[idx[0]] - midval[1][0]) > cm2 ||
198  sqr(phx[idx[2]] - midval[0][2]) + sqr(phy[idx[2]] - midval[1][2]) > cm2) split = true;
199  }
200 
201  // do extra tests at level 0, so as not to miss functions with zero error at edge midpoints
202  if(level == 0 && !split)
203  {
204  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])))) +
205  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])))) +
206  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])))) +
207  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]))));
208  split = !finite(err) || (err) > max*4*eps;
209  }
210  }
211 
212  // split the quad if the error is too large, otherwise produce two linear triangles
213  if(split)
214  {
215  if(curved)
216  {
217  for (i = 0; i < 5; i++)
218  {
219  midval[0][i] = phx[idx[i]];
220  midval[1][i] = phy[idx[i]];
221  }
222  }
223 
224  // obtain mid-edge and mid-element vertices
225  int mid0 = get_vertex(iv0, iv1, midval[0][0], midval[1][0], xval[idx[0]], yval[idx[0]]);
226  int mid1 = get_vertex(iv1, iv2, midval[0][1], midval[1][1], xval[idx[1]], yval[idx[1]]);
227  int mid2 = get_vertex(iv2, iv3, midval[0][2], midval[1][2], xval[idx[2]], yval[idx[2]]);
228  int mid3 = get_vertex(iv3, iv0, midval[0][3], midval[1][3], xval[idx[3]], yval[idx[3]]);
229  int mid4 = get_vertex(mid0, mid2, midval[0][4], midval[1][4], xval[idx[4]], yval[idx[4]]);
230 
231  // recur to sub-elements
232  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();
233  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();
234  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();
235  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();
236  return;
237  }
238  }
239 
240  // output two linear triangles,
241  if(!flip)
242  {
243  add_triangle(iv3, iv0, iv1, fns[0]->get_active_element()->marker);
244  add_triangle(iv1, iv2, iv3, fns[0]->get_active_element()->marker);
245  }
246  else
247  {
248  add_triangle(iv0, iv1, iv2, fns[0]->get_active_element()->marker);
249  add_triangle(iv2, iv3, iv0, fns[0]->get_active_element()->marker);
250  }
251  }
252 
253  void Vectorizer::set_curvature_epsilon(double curvature_epsilon)
254  {
255  this->curvature_epsilon = curvature_epsilon;
256  }
257 
258  double Vectorizer::get_curvature_epsilon()
259  {
260  return this->curvature_epsilon;
261  }
262 
263  void Vectorizer::process_solution(MeshFunction<double>* xsln, MeshFunction<double>* ysln, int xitem_orig, int yitem_orig, double eps)
264  {
265  // Important, sets the current caughtException to NULL.
266  this->caughtException = NULL;
267 
268  // sanity check
269  if(xsln == NULL || ysln == NULL)
270  throw Hermes::Exceptions::Exception("One of the solutions is NULL in Vectorizer:process_solution().");
271 
272  lock_data();
273  this->tick();
274 
275  // initialization
276  this->xitem = xitem_orig;
277  this->yitem = yitem_orig;
278  this->eps = eps;
279 
280  // estimate the required number of vertices and triangles
281  // (based on the assumption that the linear mesh will be
282  // about four-times finer than the original mesh).
283  int nn = xsln->get_mesh()->get_num_elements() + ysln->get_mesh()->get_num_elements();
284 
285  this->vertex_size = std::max(100 * nn, std::max(this->vertex_size, 50000));
286  this->triangle_size = std::max(150 * nn, std::max(this->triangle_size, 75000));
287  this->edges_size = std::max(100 * nn, std::max(this->edges_size, 50000));
288  //dashes_size = edges_size;
289 
290  vertex_count = 0;
291  triangle_count = 0;
292  edges_count = 0;
293  dashes_count = 0;
294 
295  // reuse or allocate vertex, triangle and edge arrays
296  verts = (double4*) malloc(sizeof(double4) * vertex_size);
297  this->tris = (int3*) realloc(this->tris, sizeof(int3) * this->triangle_size);
298  this->tri_markers = (int*) realloc(this->tri_markers, sizeof(int) * this->triangle_size);
299  this->edges = (int2*) realloc(this->edges, sizeof(int2) * this->edges_size);
300  this->edge_markers = (int*) realloc(this->edge_markers, sizeof(int) * this->edges_size);
301  this->dashes = (int2*) realloc(this->dashes, sizeof(int2) * dashes_size);
302  info = (int4*) malloc(sizeof(int4) * vertex_size);
303  this->empty = false;
304 
305  // initialize the hash table
306  hash_table = (int*) malloc(sizeof(int) * vertex_size);
307  memset(hash_table, 0xff, sizeof(int) * vertex_size);
308 
309  // select the linearization quadrature
310  Quad2D *old_quad_x, *old_quad_y;
311  old_quad_x = xsln->get_quad_2d();
312  old_quad_y = ysln->get_quad_2d();
313 
314  Hermes::vector<const Mesh*> meshes;
315  meshes.push_back(xsln->get_mesh());
316  meshes.push_back(ysln->get_mesh());
317 
318  // Parallelization
319  MeshFunction<double>*** fns = new MeshFunction<double>**[Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
320  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
321  {
322  fns[i] = new MeshFunction<double>*[2];
323  fns[i][0] = xsln->clone();
324  fns[i][0]->set_refmap(new RefMap);
325  fns[i][0]->set_quad_2d(&g_quad_lin);
326  fns[i][1] = ysln->clone();
327  fns[i][1]->set_refmap(new RefMap);
328  fns[i][1]->set_quad_2d(&g_quad_lin);
329  }
330 
331  Transformable*** trfs = new Transformable**[Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
332  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
333  {
334  trfs[i] = new Transformable*[2];
335  trfs[i][0] = fns[i][0];
336  trfs[i][1] = fns[i][1];
337  }
338 
339  // get the component and desired value from item.
340  if(xitem >= 0x40)
341  {
342  component_x = 1;
343  xitem >>= 6;
344  }
345  while (!(xitem & 1))
346  {
347  xitem >>= 1;
348  value_type_x++;
349  }
350  // get the component and desired value from item.
351  if(yitem >= 0x40)
352  {
353  component_y = 1;
354  yitem >>= 6;
355  }
356  while (!(yitem & 1))
357  {
358  yitem >>= 1;
359  value_type_y++;
360  }
361  xitem = xitem_orig;
362  yitem = yitem_orig;
363 
364  Traverse trav_masterMax(true);
365  unsigned int num_states = trav_masterMax.get_num_states(meshes);
366 
367  trav_masterMax.begin(meshes.size(), &(meshes.front()));
368 
369  Traverse* trav = new Traverse[Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
370 
371  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
372  {
373  trav[i].begin(meshes.size(), &(meshes.front()), trfs[i]);
374  trav[i].stack = trav_masterMax.stack;
375  }
376 
377  int state_i;
378 
379 #define CHUNKSIZE 1
380 int num_threads_used = Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads);
381 #pragma omp parallel shared(trav_masterMax) private(state_i) num_threads(num_threads_used)
382  {
383 #pragma omp for schedule(dynamic, CHUNKSIZE)
384  for(state_i = 0; state_i < num_states; state_i++)
385  {
386  try
387  {
388  Traverse::State current_state;
389 #pragma omp critical(get_next_state)
390  current_state = trav[omp_get_thread_num()].get_next_state(&trav_masterMax.top, &trav_masterMax.id);
391 
392  fns[omp_get_thread_num()][0]->set_quad_order(0, xitem);
393  fns[omp_get_thread_num()][1]->set_quad_order(0, yitem);
394  double* xval = fns[omp_get_thread_num()][0]->get_values(component_x, value_type_x);
395  double* yval = fns[omp_get_thread_num()][1]->get_values(component_y, value_type_y);
396 
397  for (unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
398  {
399  double fx = xval[i];
400  double fy = yval[i];
401 #pragma omp critical(max)
402  if(fabs(sqrt(fx*fx + fy*fy)) > max)
403  max = fabs(sqrt(fx*fx + fy*fy));
404  }
405  }
406  catch(Hermes::Exceptions::Exception& e)
407  {
408  if(this->caughtException == NULL)
409  this->caughtException = e.clone();
410  }
411  catch(std::exception& e)
412  {
413  if(this->caughtException == NULL)
414  this->caughtException = new Hermes::Exceptions::Exception(e.what());
415  }
416  }
417  }
418 
419  trav_masterMax.finish();
420  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
421  trav[i].finish();
422  delete [] trav;
423 
424  Traverse trav_master(true);
425  num_states = trav_master.get_num_states(meshes);
426 
427  trav_master.begin(meshes.size(), &(meshes.front()));
428 
429  trav = new Traverse[Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
430 
431  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
432  {
433  trav[i].begin(meshes.size(), &(meshes.front()), trfs[i]);
434  trav[i].stack = trav_master.stack;
435  }
436 
437 #pragma omp parallel shared(trav_master) private(state_i) num_threads(num_threads_used)
438  {
439 #pragma omp for schedule(dynamic, CHUNKSIZE)
440  for(state_i = 0; state_i < num_states; state_i++)
441  {
442  try
443  {
444  Traverse::State current_state;
445 
446 #pragma omp critical (get_next_state)
447  current_state = trav[omp_get_thread_num()].get_next_state(&trav_master.top, &trav_master.id);
448 
449  fns[omp_get_thread_num()][0]->set_quad_order(0, xitem);
450  fns[omp_get_thread_num()][1]->set_quad_order(0, yitem);
451  double* xval = fns[omp_get_thread_num()][0]->get_values(component_x, value_type_x);
452  double* yval = fns[omp_get_thread_num()][1]->get_values(component_y, value_type_y);
453 
454  double* x = fns[omp_get_thread_num()][0]->get_refmap()->get_phys_x(0);
455  double* y = fns[omp_get_thread_num()][0]->get_refmap()->get_phys_y(0);
456 
457  int iv[4];
458  for (unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
459  {
460  double fx = xval[i];
461  double fy = yval[i];
462 
463  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);
464  }
465 
466  // recur to sub-elements
467  if(current_state.e[0]->is_triangle())
468  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());
469  else
470  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());
471 
472  for (unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
473  process_edge(iv[i], iv[current_state.e[0]->next_vert(i)], current_state.e[0]->en[i]->marker);
474  }
475  catch(Hermes::Exceptions::Exception& e)
476  {
477  if(this->caughtException == NULL)
478  this->caughtException = e.clone();
479  }
480  catch(std::exception& e)
481  {
482  if(this->caughtException == NULL)
483  this->caughtException = new Hermes::Exceptions::Exception(e.what());
484  }
485  }
486  }
487 
488  trav_master.finish();
489  for(unsigned int i = 0; i < Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
490  {
491  trav[i].finish();
492  for(unsigned int j = 0; j < 2; j++)
493  delete fns[i][j];
494  delete [] fns[i];
495  delete [] trfs[i];
496  }
497  delete [] fns;
498  delete [] trfs;
499  delete [] trav;
500 
501  // regularize the linear mesh
502  for (int i = 0; i < this->triangle_count; i++)
503  {
504  int iv0 = tris[i][0], iv1 = tris[i][1], iv2 = tris[i][2];
505 
506  int mid0 = peek_vertex(iv0, iv1);
507  int mid1 = peek_vertex(iv1, iv2);
508  int mid2 = peek_vertex(iv2, iv0);
509  if(mid0 >= 0 || mid1 >= 0 || mid2 >= 0)
510  {
511  this->del_slot = i;
512  regularize_triangle(iv0, iv1, iv2, mid0, mid1, mid2, tri_markers[i]);
513  }
514  }
515 
516  find_min_max();
517 
518  //if(verbose_mode) print_hash_stats();
519  unlock_data();
520 
521  // select old quadratrues
522  xsln->set_quad_2d(old_quad_x);
523  ysln->set_quad_2d(old_quad_y);
524 
525  // clean up
526  ::free(this->hash_table);
527  ::free(this->info);
528  }