16 #include "linearizer.h"
19 #include "exact_solution.h"
35 double3 lin_pts_0_quad[] =
43 double3 lin_pts_1_tri[12] =
59 double3 lin_pts_1_quad[21] =
84 int quad_indices[9][5] =
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 }
93 int tri_indices[5][3] =
95 { 0, 1, 2 }, { 3, 4, 5 }, { 6, 7, 8 }, { 9, 10, 11 }, { 9, 4, 8 }
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 };
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 };
106 Linearizer::Linearizer(
bool auto_max) : LinearizerBase(auto_max), dmult(1.0), component(0), value_type(0), curvature_epsilon(1e-3)
113 tris_contours = NULL;
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)
121 if(level < LIN_MAX_LEVEL)
127 fns[0]->set_quad_order(1, item);
128 val = fns[0]->get_values(component, value_type);
130 for (i = 0; i < lin_np_tri[1]; i++)
133 #pragma omp critical(max)
134 if(finite(v) && fabs(v) > max)
137 idx = tri_indices[0];
142 RefMap* refmap = fns[0]->get_refmap();
143 phx = refmap->get_phys_x(1);
144 phy = refmap->get_phys_y(1);
148 if(this->xdisp != NULL)
151 dx = fns[1]->get_fn_values();
153 if(this->ydisp != NULL)
156 dy = fns[2]->get_fn_values();
158 for (i = 0; i < lin_np_tri[1]; i++)
160 if(this->xdisp != NULL)
161 phx[i] += dmult*dx[i];
162 if(this->ydisp != NULL)
163 phy[i] += dmult*dy[i];
170 for (i = 0; i < 3; i++)
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;
182 split = (level + 5 < eps);
186 if(!auto_max && fabs(verts[iv0][2]) > max && fabs(verts[iv1][2]) > max && fabs(verts[iv2][2]) > max)
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;
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()))
212 if(level == 0 && !split)
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;
224 for (i = 0; i < 3; i++)
226 midval[0][i] = phx[idx[i]];
227 midval[1][i] = phy[idx[i]];
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]]);
236 fns[0]->push_transform(0);
237 if(this->xdisp != NULL)
239 fns[1]->push_transform(0);
240 if(this->ydisp != NULL)
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)
247 fns[1]->pop_transform();
248 if(this->ydisp != NULL)
250 fns[2]->pop_transform();
252 fns[0]->push_transform(1);
253 if(this->xdisp != NULL)
255 fns[1]->push_transform(1);
256 if(this->ydisp != NULL)
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)
263 fns[1]->pop_transform();
264 if(this->ydisp != NULL)
266 fns[2]->pop_transform();
268 fns[0]->push_transform(2);
269 if(this->xdisp != NULL)
271 fns[1]->push_transform(2);
272 if(this->ydisp != NULL)
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)
279 fns[1]->pop_transform();
280 if(this->ydisp != NULL)
282 fns[2]->pop_transform();
284 fns[0]->push_transform(3);
285 if(this->xdisp != NULL)
287 fns[1]->push_transform(3);
288 if(this->ydisp != NULL)
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)
295 fns[1]->pop_transform();
296 if(this->ydisp != NULL)
298 fns[2]->pop_transform();
304 add_triangle(iv0, iv1, iv2, fns[0]->get_active_element()->marker);
307 void Linearizer::set_curvature_epsilon(
double curvature_epsilon)
309 this->curvature_epsilon = curvature_epsilon;
312 double Linearizer::get_curvature_epsilon()
314 return this->curvature_epsilon;
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)
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;
328 if(level < LIN_MAX_LEVEL)
335 val = fns[0]->get_values(component, value_type);
337 for (i = 0; i < lin_np_quad[1]; i++)
340 if(finite(v) && fabs(v) > max)
341 #pragma omp critical(max)
342 if(finite(v) && fabs(v) > max)
347 if(fabs(max) < 1E-10)
350 idx = quad_indices[0];
354 RefMap* refmap = fns[0]->get_refmap();
361 if(this->xdisp != NULL)
363 if(this->ydisp != NULL)
365 if(this->xdisp != NULL)
367 if(this->ydisp != NULL)
369 for (i = 0; i < lin_np_quad[1]; i++)
371 if(this->xdisp != NULL)
372 phx[i] += dmult*dx[i];
373 if(this->ydisp != NULL)
374 phy[i] += dmult*dy[i];
380 for (i = 0; i < 3; i++)
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;
390 midval[2][4] = flip ? (verts[iv0][2] + verts[iv2][2]) * 0.5 : (verts[iv1][2] + verts[iv3][2]) * 0.5;
397 split = (level < eps) ? 3 : 0;
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)
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;
416 if(level > 0 && split)
420 else if(verr > 5*herr)
426 if(split != 3 && curved)
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;
436 if(level == 0 && !split)
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;
449 for (i = 0; i < 5; i++)
451 midval[0][i] = phx[idx[i]];
452 midval[1][i] = phy[idx[i]];
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]]);
467 if(this->xdisp != NULL)
470 if(this->ydisp != NULL)
473 process_quad(fns, iv0, mid0, mid4, mid3, level + 1, val, phx, phy, quad_indices[1], curved);
475 if(this->xdisp != NULL)
478 if(this->ydisp != NULL)
483 if(this->xdisp != NULL)
486 if(this->ydisp != NULL)
489 process_quad(fns, mid0, iv1, mid1, mid4, level + 1, val, phx, phy, quad_indices[2], curved);
491 if(this->xdisp != NULL)
494 if(this->ydisp != NULL)
499 if(this->xdisp != NULL)
502 if(this->ydisp != NULL)
505 process_quad(fns, mid4, mid1, iv2, mid2, level + 1, val, phx, phy, quad_indices[3], curved);
507 if(this->xdisp != NULL)
510 if(this->ydisp != NULL)
515 if(this->xdisp != NULL)
518 if(this->ydisp != NULL)
521 process_quad(fns, mid3, mid4, mid2, iv3, level + 1, val, phx, phy, quad_indices[4], curved);
523 if(this->xdisp != NULL)
526 if(this->ydisp != NULL)
534 if(this->ydisp != NULL)
537 if(this->ydisp != NULL)
540 process_quad(fns, iv0, iv1, mid1, mid3, level + 1, val, phx, phy, quad_indices[5], curved);
542 if(this->xdisp != NULL)
545 if(this->ydisp != NULL)
550 if(this->xdisp != NULL)
553 if(this->ydisp != NULL)
556 process_quad(fns, mid3, mid1, iv2, iv3, level + 1, val, phx, phy, quad_indices[6], curved);
558 if(this->xdisp != NULL)
561 if(this->ydisp != NULL)
568 if(this->xdisp != NULL)
571 if(this->ydisp != NULL)
574 process_quad(fns, iv0, mid0, mid2, iv3, level + 1, val, phx, phy, quad_indices[7], curved);
576 if(this->xdisp != NULL)
579 if(this->ydisp != NULL)
584 if(this->xdisp != NULL)
587 if(this->ydisp != NULL)
590 process_quad(fns, mid0, iv1, iv2, mid2, level + 1, val, phx, phy, quad_indices[8], curved);
592 if(this->xdisp != NULL)
595 if(this->ydisp != NULL)
606 add_triangle(iv3, iv0, iv1, fns[0]->get_active_element()->marker);
607 add_triangle(iv1, iv2, iv3, fns[0]->get_active_element()->marker);
611 add_triangle(iv0, iv1, iv2, fns[0]->get_active_element()->marker);
612 add_triangle(iv2, iv3, iv0, fns[0]->get_active_element()->marker);
616 void Linearizer::set_displacement(MeshFunction<double>* xdisp, MeshFunction<double>* ydisp,
double dmult)
631 void Linearizer::process_solution(MeshFunction<double>* sln,
int item_,
double eps)
634 this->caughtException = NULL;
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));
662 this->vertex_count = 0;
663 this->triangle_count = 0;
664 this->edges_count = 0;
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);
674 this->hash_table = (
int*) malloc(
sizeof(
int) * this->vertex_size);
675 memset(this->hash_table, 0xff,
sizeof(
int) * this->vertex_size);
678 Quad2D *old_quad, *old_quad_x = NULL, *old_quad_y = NULL;
679 old_quad = sln->get_quad_2d();
681 old_quad_x = xdisp->get_quad_2d();
683 old_quad_y = ydisp->get_quad_2d();
687 Hermes::vector<const Mesh*> meshes;
688 meshes.push_back(sln->get_mesh());
690 meshes.push_back(xdisp->get_mesh());
692 meshes.push_back(ydisp->get_mesh());
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++)
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);
704 fns[i][1] = xdisp->clone();
706 fns[i][1]->set_quad_2d(&g_quad_lin);
710 fns[i][2] = ydisp->clone();
712 fns[i][2]->set_quad_2d(&g_quad_lin);
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++)
719 trfs[i] =
new Transformable*[3];
720 trfs[i][0] = fns[i][0];
722 trfs[i][1] = fns[i][1];
724 trfs[i][2] = fns[i][2];
727 Traverse trav_masterMax(
true);
728 unsigned int num_states = trav_masterMax.get_num_states(meshes);
730 trav_masterMax.begin(meshes.size(), &(meshes.front()));
732 Traverse* trav =
new Traverse[
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
734 for(
unsigned int i = 0; i <
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
736 trav[i].begin(meshes.size(), &(meshes.front()), trfs[i]);
737 trav[i].stack = trav_masterMax.stack;
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)
746 #pragma omp for schedule(dynamic, CHUNKSIZE)
747 for(state_i = 0; state_i < num_states; state_i++)
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);
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);
758 for (
unsigned int i = 0; i < current_state.e[0]->get_nvert(); i++)
761 #pragma omp critical (max)
762 if(this->auto_max && finite(f) && fabs(f) > this->max)
766 catch(Hermes::Exceptions::Exception& e)
768 if(this->caughtException == NULL)
769 this->caughtException = e.clone();
773 if(this->caughtException == NULL)
774 this->caughtException =
new Hermes::Exceptions::Exception(e.what());
779 trav_masterMax.finish();
780 for(
unsigned int i = 0; i <
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
784 Traverse trav_master(
true);
785 num_states = trav_master.get_num_states(meshes);
787 trav_master.begin(meshes.size(), &(meshes.front()));
789 trav =
new Traverse[
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads)];
791 for(
unsigned int i = 0; i <
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
793 trav[i].begin(meshes.size(), &(meshes.front()), trfs[i]);
794 trav[i].stack = trav_master.stack;
797 #pragma omp parallel shared(trav_master) private(state_i) num_threads(num_threads_used)
799 #pragma omp for schedule(dynamic, CHUNKSIZE)
800 for(state_i = 0; state_i < num_states; state_i++)
804 Traverse::State current_state;
806 #pragma omp critical (get_next_state)
807 current_state = trav[omp_get_thread_num()].get_next_state(&trav_master.top, &trav_master.id);
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);
814 throw Hermes::Exceptions::Exception(
"Item not defined in the solution.");
818 fns[omp_get_thread_num()][1]->set_quad_order(0,
H2D_FN_VAL);
820 fns[omp_get_thread_num()][2]->set_quad_order(0,
H2D_FN_VAL);
825 dx = fns[omp_get_thread_num()][1]->get_fn_values();
827 dy = fns[omp_get_thread_num()][2]->get_fn_values();
830 for (
unsigned int i = 0; i < current_state.e[0]->get_nvert(); 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];
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);
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());
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());
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);
852 catch(Hermes::Exceptions::Exception& e)
854 if(this->caughtException == NULL)
855 this->caughtException = e.clone();
859 if(this->caughtException == NULL)
860 this->caughtException =
new Hermes::Exceptions::Exception(e.what());
865 trav_master.finish();
866 for(
unsigned int i = 0; i <
Hermes2DApi.get_integral_param_value(Hermes::Hermes2D::numThreads); i++)
869 for(
unsigned int j = 0; j < (1 + (xdisp != NULL? 1 : 0) + (ydisp != NULL ? 1 : 0)); j++)
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;
883 if(this->caughtException != NULL)
888 throw *(this->caughtException);
892 for (
int i = 0; i < this->triangle_count; i++)
894 int iv0 = tris[i][0], iv1 = tris[i][1], iv2 = tris[i][2];
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)
902 regularize_triangle(iv0, iv1, iv2, mid0, mid1, mid2, tri_markers[i]);
911 sln->set_quad_2d(old_quad);
913 xdisp->set_quad_2d(old_quad_x);
917 ydisp->set_quad_2d(old_quad_y);
926 void Linearizer::find_min_max()
929 this->min_val = 1e100;
930 this->max_val = -1e100;
931 for (
int i = 0; i < this->vertex_count; i++)
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];
938 int Linearizer::get_vertex(
int p1,
int p2,
double x,
double y,
double value)
941 if(p1 > p2) std::swap(p1, p2);
942 int index = this->hash(p1, p2);
944 if(index < this->vertex_count)
946 i = this->hash_table[index];
947 while (i >= 0 && i < this->vertex_count)
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)
964 #pragma omp critical(realloc_vertices)
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;
976 int Linearizer::add_vertex()
978 if(this->vertex_count >= this->vertex_size)
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);
986 return this->vertex_count++;
989 void Linearizer::free()
996 if(tris_contours != NULL)
998 ::free(tris_contours);
999 tris_contours = NULL;
1002 LinearizerBase::free();
1005 Linearizer::~Linearizer()
1011 bool mode_3D,
int item,
double eps)
1013 process_solution(sln, item, eps);
1015 FILE* f = fopen(filename,
"wb");
1016 if(f == NULL)
throw Hermes::Exceptions::Exception(
"Could not open %s for writing.", filename);
1020 fprintf(f,
"# vtk DataFile Version 2.0\n");
1022 fprintf(f,
"ASCII\n\n");
1023 fprintf(f,
"DATASET UNSTRUCTURED_GRID\n");
1026 fprintf(f,
"POINTS %d %s\n", this->vertex_count,
"float");
1027 for (
int i = 0; i < this->vertex_count; i++)
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);
1035 fprintf(f,
"CELLS %d %d\n", this->triangle_count, 4 * this->triangle_count);
1036 for (
int i = 0; i < this->triangle_count; i++)
1038 fprintf(f,
"3 %d %d %d\n", this->tris[i][0], this->tris[i][1], this->tris[i][2]);
1043 fprintf(f,
"CELL_TYPES %d\n", this->triangle_count);
1044 for (
int i = 0; i < this->triangle_count; i++)
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++)
1056 fprintf(f,
"%g\n", this->verts[i][2]);
1063 void Linearizer::calc_vertices_aabb(
double* min_x,
double* max_x,
double* min_y,
double* max_y)
const
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);
1070 double3* Linearizer::get_vertices()
1074 int Linearizer::get_num_vertices()
1076 return this->vertex_count;
1079 int Linearizer::get_num_contour_triangles()
1081 return this->triangle_contours_count;
1084 int3* Linearizer::get_contour_triangles()
1086 return this->tris_contours;