18 #include "mesh_function.h"
24 static const Rect H2D_UNITY = { 0, 0, ONE, ONE };
25 Traverse::Traverse(
int spaces_size) : spaces_size(spaces_size)
29 static int get_split_and_sons(Element* e, Rect* cr, Rect* er, int4& sons)
31 uint64_t hmid = (er->l + er->r) >> 1;
32 uint64_t vmid = (er->t + er->b) >> 1;
36 if (cr->r <= hmid && cr->t <= vmid)
37 return (sons[0] = sons[1] = sons[2] = sons[3] = 0), 0;
38 else if (cr->l >= hmid && cr->t <= vmid)
39 return (sons[0] = sons[1] = sons[2] = sons[3] = 1), 0;
40 else if (cr->l >= hmid && cr->b >= vmid)
41 return (sons[0] = sons[1] = sons[2] = sons[3] = 2), 0;
42 else if (cr->r <= hmid && cr->b >= vmid)
43 return (sons[0] = sons[1] = sons[2] = sons[3] = 3), 0;
44 else if (cr->r <= hmid)
45 return (sons[0] = sons[1] = 0, sons[2] = sons[3] = 3), 1;
46 else if (cr->l >= hmid)
47 return (sons[0] = sons[1] = 1, sons[2] = sons[3] = 2), 1;
48 else if (cr->t <= vmid)
49 return (sons[0] = sons[3] = 0, sons[1] = sons[2] = 1), 2;
50 else if (cr->b >= vmid)
51 return (sons[0] = sons[3] = 3, sons[1] = sons[2] = 2), 2;
53 return (sons[0] = 0, sons[1] = 1, sons[2] = 2, sons[3] = 3), 3;
58 return (sons[0] = sons[1] = sons[2] = sons[3] = 4), 0;
59 else if (cr->b >= vmid)
60 return (sons[0] = sons[1] = sons[2] = sons[3] = 5), 0;
62 return (sons[0] = sons[1] = 4, sons[2] = sons[3] = 5), 1;
67 return (sons[0] = sons[1] = sons[2] = sons[3] = 6), 0;
68 else if (cr->l >= hmid)
69 return (sons[0] = sons[1] = sons[2] = sons[3] = 7), 0;
71 return (sons[0] = sons[3] = 6, sons[1] = sons[2] = 7), 2;
75 static void move_to_son(Rect* rnew, Rect* rold,
int son)
77 uint64_t hmid = (rold->l + rold->r) >> 1;
78 uint64_t vmid = (rold->t + rold->b) >> 1;
80 memcpy(rnew, rold,
sizeof(Rect));
84 case 0: rnew->r = hmid; rnew->t = vmid;
break;
85 case 1: rnew->l = hmid; rnew->t = vmid;
break;
86 case 2: rnew->l = hmid; rnew->b = vmid;
break;
87 case 3: rnew->r = hmid; rnew->b = vmid;
break;
88 case 4: rnew->t = vmid;
break;
89 case 5: rnew->b = vmid;
break;
90 case 6: rnew->r = hmid;
break;
91 case 7: rnew->l = hmid;
break;
95 void Traverse::init_transforms(Traverse::State* s,
unsigned char i)
98 memcpy(&r, s->er + i,
sizeof(Rect));
100 while (s->cr.l > r.l || s->cr.r < r.r || s->cr.b > r.b || s->cr.t < r.t)
102 uint64_t hmid = (r.l + r.r) >> 1;
103 uint64_t vmid = (r.t + r.b) >> 1;
106 if (s->cr.r <= hmid && s->cr.t <= vmid) son = 0;
107 else if (s->cr.l >= hmid && s->cr.t <= vmid) son = 1;
108 else if (s->cr.l >= hmid && s->cr.b >= vmid) son = 2;
109 else if (s->cr.r <= hmid && s->cr.b >= vmid) son = 3;
110 else if (s->cr.r <= hmid) son = 6;
111 else if (s->cr.l >= hmid) son = 7;
112 else if (s->cr.t <= vmid) son = 4;
113 else if (s->cr.b >= vmid) son = 5;
116 s->push_transform(son, i, s->is_triangle());
117 move_to_son(&r, &r, son);
121 Traverse::State::State()
123 memset(
this, 0,
sizeof(Traverse::State));
124 for (
int i = 0; i < 4; i++)
127 for (
int i = 0; i < num; i++)
155 Traverse::State* Traverse::State::clone(
const Traverse::State * other)
157 State* state =
new State();
159 state->num = other->num;
161 state->e =
new Element*[state->num];
162 state->sub_idx =
new uint64_t[state->num];
163 memcpy(state->e, other->e, state->num *
sizeof(Element*));
164 memcpy(state->sub_idx, other->sub_idx, state->num *
sizeof(uint64_t));
165 memcpy(state->bnd, other->bnd, 4 *
sizeof(
bool));
167 state->rep = other->rep;
168 state->rep_i = other->rep_i;
169 state->visited = other->visited;
170 state->isurf = other->isurf;
171 state->isBnd = other->isBnd;
176 Traverse::State::~State()
180 if (sub_idx !=
nullptr)
184 void Traverse::State::push_transform(
unsigned char son,
unsigned char i,
bool is_triangle)
186 this->sub_idx[i] = (sub_idx[i] << 3) + son + 1;
194 case 0: bnd[1] =
false;
break;
195 case 1: bnd[2] =
false;
break;
196 case 2: bnd[0] =
false;
break;
201 memset(bnd, 0,
sizeof(bnd));
206 if (son != 0 && son != 1 && son != 4 && son != 6 && son != 7)
208 if (son != 1 && son != 2 && son != 7 && son != 4 && son != 5)
210 if (son != 2 && son != 3 && son != 5 && son != 6 && son != 7)
212 if (son != 3 && son != 0 && son != 6 && son != 4 && son != 5)
217 uint64_t Traverse::State::get_transform(
unsigned char i)
219 return this->sub_idx[i];
222 Traverse::State* Traverse::push_state(
int* top_by_ref)
224 int* top_f = (top_by_ref ==
nullptr) ? &this->top : top_by_ref;
227 throw Hermes::Exceptions::Exception(
"Stack overflow. Increase stack size.");
229 if (stack[*top_f].e ==
nullptr)
231 stack[*top_f].e =
new Element*[num];
232 stack[*top_f].er =
new Rect[num];
233 stack[*top_f].sub_idx =
new uint64_t[num];
236 stack[*top_f].visited =
false;
237 stack[*top_f].isurf = -1;
238 memset(stack[*top_f].sub_idx, 0, num *
sizeof(uint64_t));
239 for (
int i = 0; i < 4; i++)
240 stack[*top_f].bnd[i] =
true;
241 stack[*top_f].num = this->num;
243 return stack + (*top_f)++;
246 void Traverse::set_boundary_info(State* s)
248 Element* e =
nullptr;
249 for (
int i = 0; i < num; i++)
250 if ((e = s->e[i]) !=
nullptr)
break;
252 if (e->is_triangle())
254 for (
int i = 0; i < 3; i++)
255 (s->bnd[i] = (s->bnd[i] && e->en[i]->bnd));
256 s->isBnd = s->bnd[0] || s->bnd[1] || s->bnd[2] || e->vn[0]->bnd || e->vn[1]->bnd || e->vn[2]->bnd;
260 s->bnd[0] = s->bnd[0] && (s->cr.b == 0) && e->en[0]->bnd;
261 s->bnd[1] = s->bnd[1] && (s->cr.r == ONE) && e->en[1]->bnd;
262 s->bnd[2] = s->bnd[2] && (s->cr.t == ONE) && e->en[2]->bnd;
263 s->bnd[3] = s->bnd[3] && (s->cr.l == 0) && e->en[3]->bnd;
264 s->isBnd = s->bnd[0] || s->bnd[1] || s->bnd[2] || s->bnd[3] || e->vn[0]->bnd || e->vn[1]->bnd || e->vn[2]->bnd || e->vn[3]->bnd;
268 bool Traverse::State::is_triangle()
270 bool is_triangle =
false;
271 for (
int j = 0; j < this->num; j++)
273 if (this->e[j] !=
nullptr)
275 if (this->e[j]->is_triangle())
283 template<
typename Scalar>
286 std::vector<MeshSharedPtr> meshes;
287 for (
unsigned short i = 0; i < mesh_functions.size(); i++)
288 meshes.push_back(mesh_functions[i]->get_mesh());
289 return this->
get_states(meshes, states_count);
300 int count = 0, predictedCount = 0;
301 this->num = meshes_count;
302 for (
int i = 0; i < meshes_count; i++)
303 if (meshes[i]->get_num_active_elements() > predictedCount)
304 predictedCount = meshes[i]->get_num_active_elements();
305 State** states = malloc_with_check<State*>(predictedCount);
318 while (top > 0 && (s = stack + top - 1)->visited)
334 if (
id >= meshes[0]->get_num_base_elements())
337 states_count = count;
343 for (i = 0; i < num; i++)
346 s->e[i] = meshes[i]->get_element(
id);
357 s->er[i] = H2D_UNITY;
369 if (s->is_triangle())
370 for (i = 0; i < 3; i++)
376 for (i = 0; i < num; i++)
379 if (s->e[i] !=
nullptr && s->e[i]->used)
380 if (s->sub_idx[i] == 0 && s->e[i]->active)
381 if (!s->e[i]->is_triangle())
382 init_transforms(s, i);
387 for (i = 0; i < num; i++)
389 if (s->e[i] !=
nullptr && s->e[i]->used)
390 if (!s->e[i]->active)
400 if (count > predictedCount - 1)
402 predictedCount *= 1.5;
403 states = realloc_with_check<State*>(states, predictedCount);
406 set_boundary_info(s);
413 for (
int j = 0; j < this->spaces_size; j++)
414 if (s->e[j] !=
nullptr && s->e[j]->used)
420 states[count++] = State::clone(s);
425 if (s->is_triangle())
428 for (son = 0; son <= 3; son++)
430 State* ns = push_state();
432 for (i = 0; i < num; i++)
435 if (s->e[i] ==
nullptr || !s->e[i]->used)
439 else if (s->e[i]->active)
442 ns->sub_idx[i] = s->sub_idx[i];
443 ns->push_transform(son, i,
true);
448 ns->e[i] = s->e[i]->sons[son];
450 if (ns->e[i]->active)
458 memcpy(ns->bnd, s->bnd,
sizeof(ns->bnd));
462 case 0: ns->bnd[1] =
false;
break;
463 case 1: ns->bnd[2] =
false;
break;
464 case 2: ns->bnd[0] =
false;
break;
469 memset(ns->bnd, 0,
sizeof(ns->bnd));
477 int4* current_sons =
new int4[num];
479 for (i = 0; i < num; i++)
480 if (s->e[i] !=
nullptr && !s->e[i]->active)
481 split |= get_split_and_sons(s->e[i], &s->cr, s->er + i, current_sons[i]);
486 for (son = 0; son <= 3; son++)
488 State* ns = push_state();
490 move_to_son(&ns->cr, &s->cr, son);
492 for (i = 0; i < num; i++)
494 if (s->e[i] ==
nullptr || !s->e[i]->used)
503 ns->sub_idx[i] = s->sub_idx[i];
504 ns->push_transform(son, i);
508 ns->e[i] = s->e[i]->sons[current_sons[i][son] & 3];
510 move_to_son(ns->er + i, s->er + i, current_sons[i][son]);
511 if (ns->e[i]->active)
521 int son0 = 4, son1 = 5;
522 if (split == 2) { son0 = 6; son1 = 7; }
524 for (son = son0; son <= son1; son++)
526 State* ns = push_state();
527 move_to_son(&ns->cr, &s->cr, son);
529 int j = (son == 4 || son == 6) ? 0 : 2;
530 for (i = 0; i < num; i++)
532 if (s->e[i] ==
nullptr || !s->e[i]->used)
541 ns->sub_idx[i] = s->sub_idx[i];
542 ns->push_transform(son, i);
546 ns->e[i] = s->e[i]->sons[current_sons[i][j] & 3];
547 move_to_son(ns->er + i, s->er + i, current_sons[i][j]);
548 if (ns->e[i]->active)
559 State* ns = push_state();
560 memcpy(&ns->cr, &s->cr,
sizeof(Rect));
562 for (i = 0; i < num; i++)
564 if (s->e[i] ==
nullptr || !s->e[i]->used)
568 else if (s->e[i]->active)
571 memcpy(&ns->er[i], &ns->cr,
sizeof(Rect));
575 ns->e[i] = s->e[i]->sons[current_sons[i][0] & 3];
576 move_to_son(ns->er + i, s->er + i, current_sons[i][0]);
577 if (ns->e[i]->active)
582 delete[] current_sons;
588 static void testMeshesCompliance(
int n, MeshSharedPtr* meshes)
591 int base_elem_num = meshes[0]->get_num_base_elements();
592 for (
int i = 1; i < n; i++)
593 if (base_elem_num != meshes[i]->get_num_base_elements())
594 throw Hermes::Exceptions::Exception(
"Meshes not compatible in Traverse::begin().");
597 static void testMeshesQuality(
int n, MeshSharedPtr* meshes)
600 int base_elem_num = meshes[0]->get_num_base_elements();
602 double *areas =
new double[base_elem_num];
603 memset(areas, 0, base_elem_num*
sizeof(
double));
608 double min_elem_area = 1e30;
609 for_all_base_elements_incl_inactive(e, meshes[0])
612 areas[counter] = 0.0;
615 areas[counter] = e->area;
616 if (areas[counter] < min_elem_area)
617 min_elem_area = areas[counter];
623 double tolerance = min_elem_area / 100.;
625 if (min_elem_area < 0)
626 throw Exceptions::ValueException(
"min_elem_area", 0.0, Hermes::HermesSqrtEpsilon);
628 for (
int i = 1; i < n; i++)
631 for_all_base_elements_incl_inactive(e, meshes[i])
634 if (fabs(areas[counter] - e->area) > tolerance && areas[counter] > Hermes::HermesSqrtEpsilon)
636 throw Hermes::Exceptions::Exception(
"An element is probably too distorted, try different meshing.");
644 void Traverse::begin(
int n)
651 stack =
new State[size];
652 memset(stack, 0, size *
sizeof(State));
657 void Traverse::free_state(Traverse::State* state)
661 delete[] state->sub_idx;
662 memset(state, 0,
sizeof(Traverse::State));
665 void Traverse::finish()
667 if (stack ==
nullptr)
670 for (
int i = 0; i < size; i++)
671 if (stack[i].e !=
nullptr)
672 free_state(stack + i);
678 uint64_t Traverse::init_idx(Rect* cr, Rect* er)
681 memcpy(&r, er,
sizeof(Rect));
684 while (cr->l > r.l || cr->r < r.r || cr->b > r.b || cr->t < r.t)
686 uint64_t hmid = (r.l + r.r) >> 1;
687 uint64_t vmid = (r.t + r.b) >> 1;
690 if (cr->r <= hmid && cr->t <= vmid) { son = 0; r.r = hmid; r.t = vmid; }
691 else if (cr->l >= hmid && cr->t <= vmid) { son = 1; r.l = hmid; r.t = vmid; }
692 else if (cr->l >= hmid && cr->b >= vmid) { son = 2; r.l = hmid; r.b = vmid; }
693 else if (cr->r <= hmid && cr->b >= vmid) { son = 3; r.r = hmid; r.b = vmid; }
694 else if (cr->t <= vmid) { son = 4; r.t = vmid; }
695 else if (cr->b >= vmid) { son = 5; r.b = vmid; }
696 else if (cr->r <= hmid) { son = 6; r.r = hmid; }
697 else if (cr->l >= hmid) { son = 7; r.l = hmid; }
700 idx = (idx << 3) + son + 1;
705 void Traverse::union_recurrent(Rect* cr, Element** e, Rect* er, uint64_t* idx, Element* uni)
711 for (i = 0; i < num; i++)
722 if (udsize <= uni->
id)
724 if (!udsize) udsize = 1024;
725 while (udsize <= uni->
id)
727 for (i = 0; i < num; i++)
728 unidata[i] = (UniData*)realloc(unidata[i], udsize *
sizeof(UniData));
730 for (i = 0; i < num; i++)
732 unidata[i][uni->id].e = e[i];
733 unidata[i][uni->id].idx = idx[i];
739 Element** e_new =
new Element*[num];
740 Rect* er_new =
new Rect[num];
743 int4* sons =
new int4[num];
744 uint64_t* idx_new =
new uint64_t[num];
745 memcpy(idx_new, idx, num*
sizeof(uint64_t));
747 if (uni->is_triangle())
750 unimesh->refine_element_id(uni->id);
751 for (son = 0; son <= 3; son++)
753 for (i = 0; i < num; i++)
758 idx_new[i] = (idx[i] << 3) + son + 1;
761 e_new[i] = e[i]->sons[son];
763 union_recurrent(
nullptr, e_new,
nullptr, idx_new, uni->sons[son]);
770 for (i = 0; i < num; i++)
772 split |= get_split_and_sons(e[i], cr, er + i, sons[i]);
777 unimesh->refine_element_id(uni->id, 0);
779 for (son = 0; son <= 3; son++)
781 move_to_son(&cr_new, cr, son);
782 for (i = 0; i < num; i++)
787 idx_new[i] = (idx[i] << 3) + son + 1;
791 e_new[i] = e[i]->sons[sons[i][son] & 3];
792 move_to_son(&(er_new[i]), er + i, sons[i][son]);
793 if (e_new[i]->active)
794 idx_new[i] = init_idx(&cr_new, &(er_new[i]));
797 union_recurrent(&cr_new, e_new, er_new, idx_new, uni->sons[son]);
803 unimesh->refine_element_id(uni->id, split);
805 int son0 = 4, son1 = 5;
806 if (split == 2) { son0 = 6; son1 = 7; }
808 for (son = son0; son <= son1; son++)
810 move_to_son(&cr_new, cr, son);
811 j = (son == 4 || son == 6) ? 0 : 2;
812 for (i = 0; i < num; i++)
817 idx_new[i] = (idx[i] << 3) + son + 1;
821 e_new[i] = e[i]->sons[sons[i][j] & 3];
822 move_to_son(&(er_new[i]), er + i, sons[i][j]);
823 if (e_new[i]->active)
824 idx_new[i] = init_idx(&cr_new, &(er_new[i]));
827 union_recurrent(&cr_new, e_new, er_new, idx_new, uni->sons[son & 3]);
833 memcpy(&cr_new, cr,
sizeof(Rect));
834 for (i = 0; i < num; i++)
840 e_new[i] = e[i]->sons[sons[i][0] & 3];
841 move_to_son(&(er_new[i]), er + i, sons[i][0]);
842 if (e_new[i]->active)
843 idx_new[i] = init_idx(&cr_new, &(er_new[i]));
846 union_recurrent(&cr_new, e_new, er_new, idx_new, uni);
856 UniData** Traverse::construct_union_mesh(
unsigned char n, MeshSharedPtr* meshes, MeshSharedPtr unimesh)
859 testMeshesCompliance(n, meshes);
861 Traverse traverse(n);
867 Element** e =
new Element*[n];
868 Rect* er =
new Rect[n];
871 traverse.unimesh = unimesh;
872 unimesh->copy_base(meshes[0]);
876 traverse.unidata =
new UniData*[n];
877 memset(traverse.unidata, 0,
sizeof(UniData*)* n);
879 uint64_t* idx =
new uint64_t[n];
880 memset(idx, 0, n*
sizeof(uint64_t));
883 for (
int id = 0;
id < meshes[0]->get_num_base_elements();
id++)
885 if (!meshes[0]->get_element(
id)->used)
887 for (i = 0; i < n; i++)
889 e[i] = meshes[i]->get_element(
id);
890 static const Rect H2D_UNITY = { 0, 0, ONE, ONE };
891 cr = er[i] = H2D_UNITY;
893 traverse.union_recurrent(&cr, e, er, idx, unimesh->get_element(
id));
901 return traverse.unidata;
904 template HERMES_API Traverse::State** Traverse::get_states<double>(std::vector<MeshFunctionSharedPtr<double> > mesh_functions,
unsigned int& states_count);
905 template HERMES_API Traverse::State** Traverse::get_states<std::complex<double> >(std::vector<MeshFunctionSharedPtr<std::complex<double> > > mesh_functions,
unsigned int& states_count);
::xsd::cxx::tree::id< char, ncname > id
C++ type corresponding to the ID XML Schema built-in type.
State ** get_states(std::vector< MeshSharedPtr > meshes, unsigned int &states_count)