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)
6 if(level < LIN_MAX_LEVEL)
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++)
18 double m = (sqrt(sqr(xval[i]) + sqr(yval[i])));
19 #pragma omp critical(max)
20 if(finite(m) && fabs(m) > max)
28 RefMap* refmap = fns[0]->get_refmap();
29 phx = refmap->get_phys_x(1);
30 phy = refmap->get_phys_y(1);
35 for (i = 0; i < 4; i++)
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;
47 split = (level < eps);
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;
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()))
69 if(level == 0 && !split)
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;
81 for (i = 0; i < 3; i++)
83 midval[0][i] = phx[idx[i]];
84 midval[1][i] = phy[idx[i]];
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]]);
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();
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();
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();
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();
113 add_triangle(iv0, iv1, iv2, fns[0]->get_active_element()->marker);
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)
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;
127 if(level < LIN_MAX_LEVEL)
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++)
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)
147 if(fabs(max) < 1E-10)
150 idx = quad_indices[0];
154 RefMap* refmap = fns[0]->get_refmap();
155 phx = refmap->get_phys_x(1);
156 phy = refmap->get_phys_y(1);
161 for (i = 0; i < 4; i++)
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;
175 split = (level < eps);
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;
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;
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;
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;
202 if(level == 0 && !split)
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;
217 for (i = 0; i < 5; i++)
219 midval[0][i] = phx[idx[i]];
220 midval[1][i] = phy[idx[i]];
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]]);
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();
243 add_triangle(iv3, iv0, iv1, fns[0]->get_active_element()->marker);
244 add_triangle(iv1, iv2, iv3, fns[0]->get_active_element()->marker);
248 add_triangle(iv0, iv1, iv2, fns[0]->get_active_element()->marker);
249 add_triangle(iv2, iv3, iv0, fns[0]->get_active_element()->marker);
253 void Vectorizer::set_curvature_epsilon(
double curvature_epsilon)
255 this->curvature_epsilon = curvature_epsilon;
258 double Vectorizer::get_curvature_epsilon()
260 return this->curvature_epsilon;
263 void Vectorizer::process_solution(MeshFunction<double>* xsln, MeshFunction<double>* ysln,
int xitem_orig,
int yitem_orig,
double eps)
266 this->caughtException = NULL;
269 if(xsln == NULL || ysln == NULL)
270 throw Hermes::Exceptions::Exception(
"One of the solutions is NULL in Vectorizer:process_solution().");
276 this->xitem = xitem_orig;
277 this->yitem = yitem_orig;
283 int nn = xsln->get_mesh()->get_num_elements() + ysln->get_mesh()->get_num_elements();
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));
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);
306 hash_table = (
int*) malloc(
sizeof(
int) * vertex_size);
307 memset(hash_table, 0xff,
sizeof(
int) * vertex_size);
310 Quad2D *old_quad_x, *old_quad_y;
311 old_quad_x = xsln->get_quad_2d();
312 old_quad_y = ysln->get_quad_2d();
314 Hermes::vector<const Mesh*> meshes;
315 meshes.push_back(xsln->get_mesh());
316 meshes.push_back(ysln->get_mesh());
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++)
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);
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++)
334 trfs[i] =
new Transformable*[2];
335 trfs[i][0] = fns[i][0];
336 trfs[i][1] = fns[i][1];
364 Traverse trav_masterMax(
true);
365 unsigned int num_states = trav_masterMax.get_num_states(meshes);
367 trav_masterMax.begin(meshes.size(), &(meshes.front()));
369 Traverse* trav =
new Traverse[
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
371 for(
unsigned int i = 0; i <
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
373 trav[i].begin(meshes.size(), &(meshes.front()), trfs[i]);
374 trav[i].stack = trav_masterMax.stack;
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)
383 #pragma omp for schedule(dynamic, CHUNKSIZE)
384 for(state_i = 0; state_i < num_states; state_i++)
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);
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);
397 for (
unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
401 #pragma omp critical(max)
402 if(fabs(sqrt(fx*fx + fy*fy)) > max)
403 max = fabs(sqrt(fx*fx + fy*fy));
406 catch(Hermes::Exceptions::Exception& e)
408 if(this->caughtException == NULL)
409 this->caughtException = e.clone();
413 if(this->caughtException == NULL)
414 this->caughtException =
new Hermes::Exceptions::Exception(e.what());
419 trav_masterMax.finish();
420 for(
unsigned int i = 0; i <
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
424 Traverse trav_master(
true);
425 num_states = trav_master.get_num_states(meshes);
427 trav_master.begin(meshes.size(), &(meshes.front()));
429 trav =
new Traverse[
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
431 for(
unsigned int i = 0; i <
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
433 trav[i].begin(meshes.size(), &(meshes.front()), trfs[i]);
434 trav[i].stack = trav_master.stack;
437 #pragma omp parallel shared(trav_master) private(state_i) num_threads(num_threads_used)
439 #pragma omp for schedule(dynamic, CHUNKSIZE)
440 for(state_i = 0; state_i < num_states; state_i++)
444 Traverse::State current_state;
446 #pragma omp critical (get_next_state)
447 current_state = trav[omp_get_thread_num()].get_next_state(&trav_master.top, &trav_master.id);
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);
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);
458 for (
unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
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);
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());
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());
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);
475 catch(Hermes::Exceptions::Exception& e)
477 if(this->caughtException == NULL)
478 this->caughtException = e.clone();
482 if(this->caughtException == NULL)
483 this->caughtException =
new Hermes::Exceptions::Exception(e.what());
488 trav_master.finish();
489 for(
unsigned int i = 0; i <
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
492 for(
unsigned int j = 0; j < 2; j++)
502 for (
int i = 0; i < this->triangle_count; i++)
504 int iv0 = tris[i][0], iv1 = tris[i][1], iv2 = tris[i][2];
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)
512 regularize_triangle(iv0, iv1, iv2, mid0, mid1, mid2, tri_markers[i]);
522 xsln->set_quad_2d(old_quad_x);
523 ysln->set_quad_2d(old_quad_y);
526 ::free(this->hash_table);