1 void Linearizer::process_triangle(MeshFunction<double>** fns,
int iv0,
int iv1,
int iv2,
int level,
2 double* val,
double* phx,
double* phy,
int* idx,
bool curved)
6 if(level < LIN_MAX_LEVEL)
12 fns[0]->set_quad_order(1, item);
13 val = fns[0]->get_values(component, value_type);
15 for (i = 0; i < lin_np_tri[1]; i++)
18 #pragma omp critical(max)
19 if(finite(v) && fabs(v) > max)
27 RefMap* refmap = fns[0]->get_refmap();
28 phx = refmap->get_phys_x(1);
29 phy = refmap->get_phys_y(1);
33 if(this->xdisp != NULL)
35 fns[1]->set_quad_order(1, H2D_FN_VAL);
36 dx = fns[1]->get_fn_values();
38 if(this->ydisp != NULL)
40 fns[2]->set_quad_order(1, H2D_FN_VAL);
41 dy = fns[2]->get_fn_values();
43 for (i = 0; i < lin_np_tri[1]; i++)
45 if(this->xdisp != NULL)
46 phx[i] += dmult*dx[i];
47 if(this->ydisp != NULL)
48 phy[i] += dmult*dy[i];
55 for (i = 0; i < 3; i++)
57 midval[i][0] = (verts[iv0][i] + verts[iv1][i])*0.5;
58 midval[i][1] = (verts[iv1][i] + verts[iv2][i])*0.5;
59 midval[i][2] = (verts[iv2][i] + verts[iv0][i])*0.5;
67 split = (level + 5 < eps);
71 if(!auto_max && fabs(verts[iv0][2]) > max && fabs(verts[iv1][2]) > max && fabs(verts[iv2][2]) > max)
79 double err = fabs(val[idx[0]] - midval[2][0]) +
80 fabs(val[idx[1]] - midval[2][1]) +
81 fabs(val[idx[2]] - midval[2][2]);
82 split = !finite(err) || err > max*3*eps;
88 for (i = 0; i < 3; i++)
89 if(sqr(phx[idx[i]] - midval[0][i]) + sqr(phy[idx[i]] - midval[1][i]) > sqr(fns[0]->get_active_element()->get_diameter()*this->get_curvature_epsilon()))
97 if(level == 0 && !split)
99 split = (fabs(val[8] - 0.5*(midval[2][0] + midval[2][1])) +
100 fabs(val[9] - 0.5*(midval[2][1] + midval[2][2])) +
101 fabs(val[4] - 0.5*(midval[2][2] + midval[2][0]))) > max*3*eps;
109 for (i = 0; i < 3; i++)
111 midval[0][i] = phx[idx[i]];
112 midval[1][i] = phy[idx[i]];
116 int mid0 = get_vertex(iv0, iv1, midval[0][0], midval[1][0], val[idx[0]]);
117 int mid1 = get_vertex(iv1, iv2, midval[0][1], midval[1][1], val[idx[1]]);
118 int mid2 = get_vertex(iv2, iv0, midval[0][2], midval[1][2], val[idx[2]]);
121 fns[0]->push_transform(0);
122 if(this->xdisp != NULL)
124 fns[1]->push_transform(0);
125 if(this->ydisp != NULL)
127 fns[2]->push_transform(0);
128 process_triangle(fns, iv0, mid0, mid2, level + 1, val, phx, phy, tri_indices[1], curved);
129 fns[0]->pop_transform();
130 if(this->xdisp != NULL)
132 fns[1]->pop_transform();
133 if(this->ydisp != NULL)
135 fns[2]->pop_transform();
137 fns[0]->push_transform(1);
138 if(this->xdisp != NULL)
140 fns[1]->push_transform(1);
141 if(this->ydisp != NULL)
143 fns[2]->push_transform(1);
144 process_triangle(fns, mid0, iv1, mid1, level + 1, val, phx, phy, tri_indices[2], curved);
145 fns[0]->pop_transform();
146 if(this->xdisp != NULL)
148 fns[1]->pop_transform();
149 if(this->ydisp != NULL)
151 fns[2]->pop_transform();
153 fns[0]->push_transform(2);
154 if(this->xdisp != NULL)
156 fns[1]->push_transform(2);
157 if(this->ydisp != NULL)
159 fns[2]->push_transform(2);
160 process_triangle(fns, mid2, mid1, iv2, level + 1, val, phx, phy, tri_indices[3], curved);
161 fns[0]->pop_transform();
162 if(this->xdisp != NULL)
164 fns[1]->pop_transform();
165 if(this->ydisp != NULL)
167 fns[2]->pop_transform();
169 fns[0]->push_transform(3);
170 if(this->xdisp != NULL)
172 fns[1]->push_transform(3);
173 if(this->ydisp != NULL)
175 fns[2]->push_transform(3);
176 process_triangle(fns, mid1, mid2, mid0, level + 1, val, phx, phy, tri_indices[4], curved);
177 fns[0]->pop_transform();
178 if(this->xdisp != NULL)
180 fns[1]->pop_transform();
181 if(this->ydisp != NULL)
183 fns[2]->pop_transform();
189 add_triangle(iv0, iv1, iv2, fns[0]->get_active_element()->marker);
192 void Linearizer::set_curvature_epsilon(
double curvature_epsilon)
194 this->curvature_epsilon = curvature_epsilon;
197 double Linearizer::get_curvature_epsilon()
199 return this->curvature_epsilon;
202 void Linearizer::process_quad(MeshFunction<double>** fns,
int iv0,
int iv1,
int iv2,
int iv3,
int level,
203 double* val,
double* phx,
double* phy,
int* idx,
bool curved)
208 int a = (verts[iv0][2] > verts[iv1][2]) ? iv0 : iv1;
209 int b = (verts[iv2][2] > verts[iv3][2]) ? iv2 : iv3;
210 a = (verts[a][2] > verts[b][2]) ? a : b;
211 int flip = (a == iv1 || a == iv3) ? 1 : 0;
213 if(level < LIN_MAX_LEVEL)
219 fns[0]->set_quad_order(1, item);
220 val = fns[0]->get_values(component, value_type);
222 for (i = 0; i < lin_np_quad[1]; i++)
225 if(finite(v) && fabs(v) > max)
226 #pragma omp critical(max)
227 if(finite(v) && fabs(v) > max)
232 if(fabs(max) < 1E-10)
235 idx = quad_indices[0];
239 RefMap* refmap = fns[0]->get_refmap();
240 phx = refmap->get_phys_x(1);
241 phy = refmap->get_phys_y(1);
246 if(this->xdisp != NULL)
247 fns[1]->set_quad_order(1, H2D_FN_VAL);
248 if(this->ydisp != NULL)
249 fns[2]->set_quad_order(1, H2D_FN_VAL);
250 if(this->xdisp != NULL)
251 dx = fns[1]->get_fn_values();
252 if(this->ydisp != NULL)
253 dy = fns[2]->get_fn_values();
254 for (i = 0; i < lin_np_quad[1]; i++)
256 if(this->xdisp != NULL)
257 phx[i] += dmult*dx[i];
258 if(this->ydisp != NULL)
259 phy[i] += dmult*dy[i];
265 for (i = 0; i < 3; i++)
267 midval[i][0] = (verts[iv0][i] + verts[iv1][i]) * 0.5;
268 midval[i][1] = (verts[iv1][i] + verts[iv2][i]) * 0.5;
269 midval[i][2] = (verts[iv2][i] + verts[iv3][i]) * 0.5;
270 midval[i][3] = (verts[iv3][i] + verts[iv0][i]) * 0.5;
271 midval[i][4] = (midval[i][0] + midval[i][2]) * 0.5;
275 midval[2][4] = flip ? (verts[iv0][2] + verts[iv2][2]) * 0.5 : (verts[iv1][2] + verts[iv3][2]) * 0.5;
282 split = (level < eps) ? 3 : 0;
286 if(!auto_max && fabs(verts[iv0][2]) > max && fabs(verts[iv1][2]) > max
287 && fabs(verts[iv2][2]) > max && fabs(verts[iv3][2]) > max)
295 double herr = fabs(val[idx[1]] - midval[2][1]) + fabs(val[idx[3]] - midval[2][3]);
296 double verr = fabs(val[idx[0]] - midval[2][0]) + fabs(val[idx[2]] - midval[2][2]);
297 double err = fabs(val[idx[4]] - midval[2][4]) + herr + verr;
298 split = (!finite(err) || err > max*4*eps) ? 3 : 0;
301 if(level > 0 && split)
305 else if(verr > 5*herr)
311 if(split != 3 && curved)
313 double cm2 = sqr(fns[0]->get_active_element()->get_diameter()*this->get_curvature_epsilon());
314 if(sqr(phx[idx[1]] - midval[0][1]) + sqr(phy[idx[1]] - midval[1][1]) > cm2 ||
315 sqr(phx[idx[3]] - midval[0][3]) + sqr(phy[idx[3]] - midval[1][3]) > cm2) split |= 1;
316 if(sqr(phx[idx[0]] - midval[0][0]) + sqr(phy[idx[0]] - midval[1][0]) > cm2 ||
317 sqr(phx[idx[2]] - midval[0][2]) + sqr(phy[idx[2]] - midval[1][2]) > cm2) split |= 2;
321 if(level == 0 && !split)
323 split = ((fabs(val[13] - 0.5*(midval[2][0] + midval[2][1])) +
324 fabs(val[17] - 0.5*(midval[2][1] + midval[2][2])) +
325 fabs(val[20] - 0.5*(midval[2][2] + midval[2][3])) +
326 fabs(val[9] - 0.5*(midval[2][3] + midval[2][0]))) > max*4*eps) ? 3 : 0;
334 for (i = 0; i < 5; i++)
336 midval[0][i] = phx[idx[i]];
337 midval[1][i] = phy[idx[i]];
341 int mid0, mid1, mid2, mid3, mid4;
342 if(split != 1) mid0 = get_vertex(iv0, iv1, midval[0][0], midval[1][0], val[idx[0]]);
343 if(split != 2) mid1 = get_vertex(iv1, iv2, midval[0][1], midval[1][1], val[idx[1]]);
344 if(split != 1) mid2 = get_vertex(iv2, iv3, midval[0][2], midval[1][2], val[idx[2]]);
345 if(split != 2) mid3 = get_vertex(iv3, iv0, midval[0][3], midval[1][3], val[idx[3]]);
346 if(split == 3) mid4 = get_vertex(mid0, mid2, midval[0][4], midval[1][4], val[idx[4]]);
351 fns[0]->push_transform(0);
352 if(this->xdisp != NULL)
354 fns[1]->push_transform(0);
355 if(this->ydisp != NULL)
357 fns[2]->push_transform(0);
358 process_quad(fns, iv0, mid0, mid4, mid3, level + 1, val, phx, phy, quad_indices[1], curved);
359 fns[0]->pop_transform();
360 if(this->xdisp != NULL)
362 fns[1]->pop_transform();
363 if(this->ydisp != NULL)
365 fns[2]->pop_transform();
367 fns[0]->push_transform(1);
368 if(this->xdisp != NULL)
370 fns[1]->push_transform(1);
371 if(this->ydisp != NULL)
373 fns[2]->push_transform(1);
374 process_quad(fns, mid0, iv1, mid1, mid4, level + 1, val, phx, phy, quad_indices[2], curved);
375 fns[0]->pop_transform();
376 if(this->xdisp != NULL)
378 fns[1]->pop_transform();
379 if(this->ydisp != NULL)
381 fns[2]->pop_transform();
383 fns[0]->push_transform(2);
384 if(this->xdisp != NULL)
386 fns[1]->push_transform(2);
387 if(this->ydisp != NULL)
389 fns[2]->push_transform(2);
390 process_quad(fns, mid4, mid1, iv2, mid2, level + 1, val, phx, phy, quad_indices[3], curved);
391 fns[0]->pop_transform();
392 if(this->xdisp != NULL)
394 fns[1]->pop_transform();
395 if(this->ydisp != NULL)
397 fns[2]->pop_transform();
399 fns[0]->push_transform(3);
400 if(this->xdisp != NULL)
402 fns[1]->push_transform(3);
403 if(this->ydisp != NULL)
405 fns[2]->push_transform(3);
406 process_quad(fns, mid3, mid4, mid2, iv3, level + 1, val, phx, phy, quad_indices[4], curved);
407 fns[0]->pop_transform();
408 if(this->xdisp != NULL)
410 fns[1]->pop_transform();
411 if(this->ydisp != NULL)
413 fns[2]->pop_transform();
418 fns[0]->push_transform(4);
419 if(this->ydisp != NULL)
421 fns[1]->push_transform(4);
422 if(this->ydisp != NULL)
424 fns[2]->push_transform(4);
425 process_quad(fns, iv0, iv1, mid1, mid3, level + 1, val, phx, phy, quad_indices[5], curved);
426 fns[0]->pop_transform();
427 if(this->xdisp != NULL)
429 fns[1]->pop_transform();
430 if(this->ydisp != NULL)
432 fns[2]->pop_transform();
434 fns[0]->push_transform(5);
435 if(this->xdisp != NULL)
437 fns[1]->push_transform(5);
438 if(this->ydisp != NULL)
440 fns[2]->push_transform(5);
441 process_quad(fns, mid3, mid1, iv2, iv3, level + 1, val, phx, phy, quad_indices[6], curved);
442 fns[0]->pop_transform();
443 if(this->xdisp != NULL)
445 fns[1]->pop_transform();
446 if(this->ydisp != NULL)
448 fns[2]->pop_transform();
452 fns[0]->push_transform(6);
453 if(this->xdisp != NULL)
455 fns[1]->push_transform(6);
456 if(this->ydisp != NULL)
458 fns[2]->push_transform(6);
459 process_quad(fns, iv0, mid0, mid2, iv3, level + 1, val, phx, phy, quad_indices[7], curved);
460 fns[0]->pop_transform();
461 if(this->xdisp != NULL)
463 fns[1]->pop_transform();
464 if(this->ydisp != NULL)
466 fns[2]->pop_transform();
468 fns[0]->push_transform(7);
469 if(this->xdisp != NULL)
471 fns[1]->push_transform(7);
472 if(this->ydisp != NULL)
474 fns[2]->push_transform(7);
475 process_quad(fns, mid0, iv1, iv2, mid2, level + 1, val, phx, phy, quad_indices[8], curved);
476 fns[0]->pop_transform();
477 if(this->xdisp != NULL)
479 fns[1]->pop_transform();
480 if(this->ydisp != NULL)
482 fns[2]->pop_transform();
491 add_triangle(iv3, iv0, iv1, fns[0]->get_active_element()->marker);
492 add_triangle(iv1, iv2, iv3, fns[0]->get_active_element()->marker);
496 add_triangle(iv0, iv1, iv2, fns[0]->get_active_element()->marker);
497 add_triangle(iv2, iv3, iv0, fns[0]->get_active_element()->marker);
501 void Linearizer::process_solution(MeshFunction<double>* sln,
int item_,
double eps)
504 this->caughtException = NULL;
528 this->vertex_size = std::max(100 * sln->get_mesh()->get_num_elements(), std::max(this->vertex_size, 50000));
529 this->triangle_size = std::max(150 * sln->get_mesh()->get_num_elements(), std::max(this->triangle_size, 75000));
530 this->edges_size = std::max(100 * sln->get_mesh()->get_num_elements(), std::max(this->edges_size, 50000));
532 this->vertex_count = 0;
533 this->triangle_count = 0;
534 this->edges_count = 0;
536 this->verts = (double3*) realloc(this->verts,
sizeof(double3) * this->vertex_size);
537 this->tris = (int3*) realloc(this->tris,
sizeof(int3) * this->triangle_size);
538 this->tri_markers = (
int*) realloc(this->tri_markers,
sizeof(
int) * this->triangle_size);
539 this->edges = (int2*) realloc(this->edges,
sizeof(int2) * this->edges_size);
540 this->edge_markers = (
int*) realloc(this->edge_markers,
sizeof(
int) * this->edges_size);
541 this->info = (int4*) malloc(
sizeof(int4) * this->vertex_size);
544 this->hash_table = (
int*) malloc(
sizeof(
int) * this->vertex_size);
545 memset(this->hash_table, 0xff,
sizeof(
int) * this->vertex_size);
548 Quad2D *old_quad, *old_quad_x = NULL, *old_quad_y = NULL;
549 old_quad = sln->get_quad_2d();
551 old_quad_x = xdisp->get_quad_2d();
553 old_quad_y = ydisp->get_quad_2d();
557 Hermes::vector<const Mesh*> meshes;
558 meshes.push_back(sln->get_mesh());
560 meshes.push_back(xdisp->get_mesh());
562 meshes.push_back(ydisp->get_mesh());
565 MeshFunction<double>*** fns =
new MeshFunction<double>**[
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
566 for(
unsigned int i = 0; i <
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
568 fns[i] =
new MeshFunction<double>*[3];
569 fns[i][0] = sln->clone();
570 fns[i][0]->set_refmap(
new RefMap);
571 fns[i][0]->set_quad_2d(&g_quad_lin);
574 fns[i][1] = xdisp->clone();
576 fns[i][1]->set_quad_2d(&g_quad_lin);
580 fns[i][2] = ydisp->clone();
582 fns[i][2]->set_quad_2d(&g_quad_lin);
586 Transformable*** trfs =
new Transformable**[
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
587 for(
unsigned int i = 0; i <
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
589 trfs[i] =
new Transformable*[3];
590 trfs[i][0] = fns[i][0];
592 trfs[i][1] = fns[i][1];
594 trfs[i][2] = fns[i][2];
597 Traverse trav_masterMax(
true);
598 unsigned int num_states = trav_masterMax.get_num_states(meshes);
600 trav_masterMax.begin(meshes.size(), &(meshes.front()));
602 Traverse* trav =
new Traverse[
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
604 for(
unsigned int i = 0; i <
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
606 trav[i].begin(meshes.size(), &(meshes.front()), trfs[i]);
607 trav[i].stack = trav_masterMax.stack;
613 int num_threads_used =
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads);
614 #pragma omp parallel shared(trav_masterMax) private(state_i) num_threads(num_threads_used)
616 #pragma omp for schedule(dynamic, CHUNKSIZE)
617 for(state_i = 0; state_i < num_states; state_i++)
621 Traverse::State current_state;
622 #pragma omp critical(get_next_state)
623 current_state = trav[omp_get_thread_num()].get_next_state(&trav_masterMax.top, &trav_masterMax.id);
625 fns[omp_get_thread_num()][0]->set_quad_order(0, this->item);
626 double* val = fns[omp_get_thread_num()][0]->get_values(component, value_type);
628 for (
unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
631 #pragma omp critical (max)
632 if(this->auto_max && finite(f) && fabs(f) > this->max)
636 catch(Hermes::Exceptions::Exception& e)
638 if(this->caughtException == NULL)
639 this->caughtException = e.clone();
643 if(this->caughtException == NULL)
644 this->caughtException =
new Hermes::Exceptions::Exception(e.what());
649 trav_masterMax.finish();
650 for(
unsigned int i = 0; i <
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
654 Traverse trav_master(
true);
655 num_states = trav_master.get_num_states(meshes);
657 trav_master.begin(meshes.size(), &(meshes.front()));
659 trav =
new Traverse[
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
661 for(
unsigned int i = 0; i <
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
663 trav[i].begin(meshes.size(), &(meshes.front()), trfs[i]);
664 trav[i].stack = trav_master.stack;
667 #pragma omp parallel shared(trav_master) private(state_i) num_threads(num_threads_used)
669 #pragma omp for schedule(dynamic, CHUNKSIZE)
670 for(state_i = 0; state_i < num_states; state_i++)
674 Traverse::State current_state;
676 #pragma omp critical (get_next_state)
677 current_state = trav[omp_get_thread_num()].get_next_state(&trav_master.top, &trav_master.id);
679 fns[omp_get_thread_num()][0]->set_quad_order(0, this->item);
680 double* val = fns[omp_get_thread_num()][0]->get_values(component, value_type);
684 throw Hermes::Exceptions::Exception(
"Item not defined in the solution.");
688 fns[omp_get_thread_num()][1]->set_quad_order(0, H2D_FN_VAL);
690 fns[omp_get_thread_num()][2]->set_quad_order(0, H2D_FN_VAL);
695 dx = fns[omp_get_thread_num()][1]->get_fn_values();
697 dy = fns[omp_get_thread_num()][2]->get_fn_values();
700 for (
unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
703 double x_disp = fns[omp_get_thread_num()][0]->get_refmap()->get_phys_x(0)[i];
704 double y_disp = fns[omp_get_thread_num()][0]->get_refmap()->get_phys_y(0)[i];
705 if(this->xdisp != NULL)
706 x_disp += dmult * dx[i];
707 if(this->ydisp != NULL)
708 y_disp += dmult * dy[i];
710 iv[i] = this->get_vertex(-fns[omp_get_thread_num()][0]->get_active_element()->vn[i]->
id, -fns[omp_get_thread_num()][0]->get_active_element()->vn[i]->
id, x_disp, y_disp, f);
714 if(current_state.e[0]->is_triangle())
715 process_triangle(fns[omp_get_thread_num()], iv[0], iv[1], iv[2], 0, NULL, NULL, NULL, NULL, current_state.e[0]->is_curved());
717 process_quad(fns[omp_get_thread_num()], iv[0], iv[1], iv[2], iv[3], 0, NULL, NULL, NULL, NULL, current_state.e[0]->is_curved());
719 for (
unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
720 process_edge(iv[i], iv[current_state.e[0]->next_vert(i)], current_state.e[0]->en[i]->marker);
722 catch(Hermes::Exceptions::Exception& e)
724 if(this->caughtException == NULL)
725 this->caughtException = e.clone();
729 if(this->caughtException == NULL)
730 this->caughtException =
new Hermes::Exceptions::Exception(e.what());
735 trav_master.finish();
736 for(
unsigned int i = 0; i <
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
739 for(
unsigned int j = 0; j < (1 + (xdisp != NULL? 1 : 0) + (ydisp != NULL ? 1 : 0)); j++)
749 this->tris_contours = (int3*) realloc(this->tris_contours,
sizeof(int3) * this->triangle_count);
750 memcpy(this->tris_contours, this->tris, this->triangle_count *
sizeof(int3));
751 triangle_contours_count = this->triangle_count;
753 if(this->caughtException != NULL)
758 throw *(this->caughtException);
762 for (
int i = 0; i < this->triangle_count; i++)
764 int iv0 = tris[i][0], iv1 = tris[i][1], iv2 = tris[i][2];
766 int mid0 = peek_vertex(iv0, iv1);
767 int mid1 = peek_vertex(iv1, iv2);
768 int mid2 = peek_vertex(iv2, iv0);
769 if(mid0 >= 0 || mid1 >= 0 || mid2 >= 0)
772 regularize_triangle(iv0, iv1, iv2, mid0, mid1, mid2, tri_markers[i]);
781 sln->set_quad_2d(old_quad);
783 xdisp->set_quad_2d(old_quad_x);
787 ydisp->set_quad_2d(old_quad_y);