Hermes2D  3.0
traverse.cpp
1 // This file is part of Hermes2D.
2 //
3 // Hermes2D is free software: you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License as published by
5 // the Free Software Foundation, either version 2 of the License, or
6 // (at your option) any later version.
7 //
8 // Hermes2D is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 // GNU General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with Hermes2D. If not, see <http://www.gnu.org/licenses/>.
15 
16 #include "mesh.h"
17 #include "traverse.h"
18 #include "mesh_function.h"
19 
20 namespace Hermes
21 {
22  namespace Hermes2D
23  {
24  static const Rect H2D_UNITY = { 0, 0, ONE, ONE };
25  Traverse::Traverse(int spaces_size) : spaces_size(spaces_size)
26  {
27  }
28 
29  static int get_split_and_sons(Element* e, Rect* cr, Rect* er, int4& sons)
30  {
31  uint64_t hmid = (er->l + er->r) >> 1;
32  uint64_t vmid = (er->t + er->b) >> 1;
33 
34  if (e->bsplit())
35  {
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;
52  else
53  return (sons[0] = 0, sons[1] = 1, sons[2] = 2, sons[3] = 3), 3;
54  }
55  else if (e->hsplit())
56  {
57  if (cr->t <= vmid)
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;
61  else
62  return (sons[0] = sons[1] = 4, sons[2] = sons[3] = 5), 1;
63  }
64  else // e->vsplit()
65  {
66  if (cr->r <= hmid)
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;
70  else
71  return (sons[0] = sons[3] = 6, sons[1] = sons[2] = 7), 2;
72  }
73  }
74 
75  static void move_to_son(Rect* rnew, Rect* rold, int son)
76  {
77  uint64_t hmid = (rold->l + rold->r) >> 1;
78  uint64_t vmid = (rold->t + rold->b) >> 1;
79  if (rnew != rold)
80  memcpy(rnew, rold, sizeof(Rect));
81 
82  switch (son)
83  {
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;
92  }
93  }
94 
95  void Traverse::init_transforms(Traverse::State* s, unsigned char i)
96  {
97  Rect r;
98  memcpy(&r, s->er + i, sizeof(Rect));
99 
100  while (s->cr.l > r.l || s->cr.r < r.r || s->cr.b > r.b || s->cr.t < r.t)
101  {
102  uint64_t hmid = (r.l + r.r) >> 1;
103  uint64_t vmid = (r.t + r.b) >> 1;
104  int son;
105 
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;
114  else assert(0);
115 
116  s->push_transform(son, i, s->is_triangle());
117  move_to_son(&r, &r, son);
118  }
119  }
120 
121  Traverse::State::State()
122  {
123  memset(this, 0, sizeof(Traverse::State));
124  for (int i = 0; i < 4; i++)
125  this->bnd[i] = true;
126  cr = H2D_UNITY;
127  for (int i = 0; i < num; i++)
128  er[i] = H2D_UNITY;
129  isurf = -1;
130  }
131 
132  /*
133  void Traverse::State::operator=(const State * other)
134  {
135  // Delete part.
136  if(e != nullptr)
137  delete [] e;
138  if(sub_idx != nullptr)
139  delete [] sub_idx;
140 
141  this->num = other->num;
142 
143  this->e = new Element*[num];
144  this->sub_idx = new uint64_t[num];
145  memcpy(this->e, other->e, num * sizeof(Element*));
146  memcpy(this->sub_idx, other->sub_idx, num * sizeof(uint64_t));
147  memcpy(this->bnd, other->bnd, 4 * sizeof(bool));
148 
149  this->rep = other->rep;
150  this->visited = other->visited;
151  this->isurf = other->isurf;
152  this->isBnd = other->isBnd;
153  }
154  */
155  Traverse::State* Traverse::State::clone(const Traverse::State * other)
156  {
157  State* state = new State();
158 
159  state->num = other->num;
160 
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));
166 
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;
172 
173  return state;
174  }
175 
176  Traverse::State::~State()
177  {
178  if (e != nullptr)
179  delete[] e;
180  if (sub_idx != nullptr)
181  delete[] sub_idx;
182  }
183 
184  void Traverse::State::push_transform(unsigned char son, unsigned char i, bool is_triangle)
185  {
186  this->sub_idx[i] = (sub_idx[i] << 3) + son + 1;
187 
188  if (is_triangle)
189  {
190  if (son < 3)
191  {
192  switch (son)
193  {
194  case 0: bnd[1] = false; break;
195  case 1: bnd[2] = false; break;
196  case 2: bnd[0] = false; break;
197  }
198  }
199  else
200  {
201  memset(bnd, 0, sizeof(bnd));
202  }
203  }
204  else
205  {
206  if (son != 0 && son != 1 && son != 4 && son != 6 && son != 7)
207  bnd[0] = false;
208  if (son != 1 && son != 2 && son != 7 && son != 4 && son != 5)
209  bnd[1] = false;
210  if (son != 2 && son != 3 && son != 5 && son != 6 && son != 7)
211  bnd[2] = false;
212  if (son != 3 && son != 0 && son != 6 && son != 4 && son != 5)
213  bnd[3] = false;
214  }
215  }
216 
217  uint64_t Traverse::State::get_transform(unsigned char i)
218  {
219  return this->sub_idx[i];
220  }
221 
222  Traverse::State* Traverse::push_state(int* top_by_ref)
223  {
224  int* top_f = (top_by_ref == nullptr) ? &this->top : top_by_ref;
225 
226  if (*top_f >= size)
227  throw Hermes::Exceptions::Exception("Stack overflow. Increase stack size.");
228 
229  if (stack[*top_f].e == nullptr)
230  {
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];
234  }
235 
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;
242 
243  return stack + (*top_f)++;
244  }
245 
246  void Traverse::set_boundary_info(State* s)
247  {
248  Element* e = nullptr;
249  for (int i = 0; i < num; i++)
250  if ((e = s->e[i]) != nullptr) break;
251 
252  if (e->is_triangle())
253  {
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;
257  }
258  else
259  {
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;
265  }
266  }
267 
268  bool Traverse::State::is_triangle()
269  {
270  bool is_triangle = false;
271  for (int j = 0; j < this->num; j++)
272  {
273  if (this->e[j] != nullptr)
274  {
275  if (this->e[j]->is_triangle())
276  is_triangle = true;
277  break;
278  }
279  }
280  return is_triangle;
281  }
282 
283  template<typename Scalar>
284  Traverse::State** Traverse::get_states(std::vector<MeshFunctionSharedPtr<Scalar> > mesh_functions, unsigned int& states_count)
285  {
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);
290  }
291 
292  Traverse::State** Traverse::get_states(std::vector<MeshSharedPtr> meshes, unsigned int& states_count)
293  {
294  return Traverse::get_states(&meshes[0], meshes.size(), states_count);
295  }
296 
297  Traverse::State** Traverse::get_states(MeshSharedPtr* meshes, unsigned short meshes_count, unsigned int& states_count)
298  {
299  // This will be returned.
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);
306 
307  this->begin(num);
308 
309  int id = 0;
310 
311  while (1)
312  {
313  int i, son;
314  // When the stack of states is not empty (i.e. not at the beginning) the function starts here.
315  // If the top state was visited already, we are returning through it:
316  // undo all its transformations, pop it and continue with a non-visited one
317  State* s;
318  while (top > 0 && (s = stack + top - 1)->visited)
319  (top)--;
320 
321  // The stack is empty, take next base element
322  // The process starts here (at the beginning the stack is always empty, i.e. top == 0)
323  if (top <= 0)
324  {
325  // Push the state of a new_ base element.
326  // This function only allocates memory for the new_ state,
327  // with as many Elements* as there are meshes in this stage.
328  s = push_state();
329  s->cr = H2D_UNITY;
330  while (1)
331  {
332  // No more base elements? we're finished.
333  // Id is set to zero at the beginning by the function trav.begin(..).
334  if (id >= meshes[0]->get_num_base_elements())
335  {
336  this->finish();
337  states_count = count;
338  return states;
339  }
340 
341  int nused = 0;
342  // The variable num is the number of meshes in the stage
343  for (i = 0; i < num; i++)
344  {
345  // Retrieve the Element with this id on the i-th mesh.
346  s->e[i] = meshes[i]->get_element(id);
347  if (!s->e[i]->used)
348  {
349  s->e[i] = nullptr;
350  continue;
351  }
352  else
353  {
354  s->rep = s->e[i];
355  s->rep_i = 0;
356  }
357  s->er[i] = H2D_UNITY;
358  nused++;
359  }
360  // If there is any used element in this stage we continue with the calculation
361  // (break this cycle looking for such an element id).
362  if (nused)
363  break;
364  (id)++;
365  }
366 
367  (id)++;
368 
369  if (s->is_triangle())
370  for (i = 0; i < 3; i++)
371  s->bnd[i] = true;
372  }
373 
374  // Entering a new_ state, perform transformations.
375  s->visited = true;
376  for (i = 0; i < num; i++)
377  {
378  // ..where the element is used ..
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);
383  }
384 
385  // Is this the leaf state?
386  bool leaf = true;
387  for (i = 0; i < num; i++)
388  {
389  if (s->e[i] != nullptr && s->e[i]->used)
390  if (!s->e[i]->active)
391  {
392  leaf = false;
393  break;
394  }
395  }
396 
397  // if yes, set boundary flags and return the state
398  if (leaf)
399  {
400  if (count > predictedCount - 1)
401  {
402  predictedCount *= 1.5;
403  states = realloc_with_check<State*>(states, predictedCount);
404  }
405 
406  set_boundary_info(s);
407  s->rep = nullptr;
408  // EXTREMELY IMPORTANT.
409  // The for loop here is NOT through all meshes, but only
410  // through the spaces.
411  // The reason is not to include states that only have elements
412  // on meshes that are not a part of the weak form.
413  for (int j = 0; j < this->spaces_size; j++)
414  if (s->e[j] != nullptr && s->e[j]->used)
415  {
416  s->rep = s->e[j];
417  s->rep_i = j;
418  }
419  if (s->rep)
420  states[count++] = State::clone(s);
421  continue;
422  }
423 
424  // Triangle: push son states
425  if (s->is_triangle())
426  {
427  // Triangle always has 4 sons.
428  for (son = 0; son <= 3; son++)
429  {
430  State* ns = push_state();
431  // For every mesh..
432  for (i = 0; i < num; i++)
433  {
434  // ..if the element is not used.
435  if (s->e[i] == nullptr || !s->e[i]->used)
436  {
437  ns->e[i] = nullptr;
438  }
439  else if (s->e[i]->active)
440  {
441  ns->e[i] = s->e[i];
442  ns->sub_idx[i] = s->sub_idx[i];
443  ns->push_transform(son, i, true);
444  }
445  // ..we move to the son.
446  else
447  {
448  ns->e[i] = s->e[i]->sons[son];
449  // If the son's element is active.
450  if (ns->e[i]->active)
451  ns->sub_idx[i] = 0;
452  }
453  }
454 
455  // Determine boundary flags and positions for the new_ state.
456  if (son < 3)
457  {
458  memcpy(ns->bnd, s->bnd, sizeof(ns->bnd));
459 
460  switch (son)
461  {
462  case 0: ns->bnd[1] = false; break;
463  case 1: ns->bnd[2] = false; break;
464  case 2: ns->bnd[0] = false; break;
465  }
466  }
467  else
468  {
469  memset(ns->bnd, 0, sizeof(ns->bnd));
470  }
471  }
472  }
473  // Quad: this is a little more complicated, same principle, though.
474  else
475  {
476  // Obtain split types and son numbers for the current rectangle on all elements.
477  int4* current_sons = new int4[num];
478  int split = 0;
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]);
482 
483  // Both splits: recur to four sons, similar to triangles.
484  if (split == 3)
485  {
486  for (son = 0; son <= 3; son++)
487  {
488  State* ns = push_state();
489  // Sets the son's "base" rectangle to the correct one.
490  move_to_son(&ns->cr, &s->cr, son);
491 
492  for (i = 0; i < num; i++)
493  {
494  if (s->e[i] == nullptr || !s->e[i]->used)
495  {
496  ns->e[i] = nullptr;
497  }
498  else
499  {
500  if (s->e[i]->active)
501  {
502  ns->e[i] = s->e[i];
503  ns->sub_idx[i] = s->sub_idx[i];
504  ns->push_transform(son, i);
505  }
506  else
507  {
508  ns->e[i] = s->e[i]->sons[current_sons[i][son] & 3];
509  // Sets the son's "current mesh" rectangle correctly.
510  move_to_son(ns->er + i, s->er + i, current_sons[i][son]);
511  if (ns->e[i]->active)
512  ns->sub_idx[i] = 0;
513  }
514  }
515  }
516  }
517  }
518  // V or h split, recur to two sons.
519  else if (split > 0)
520  {
521  int son0 = 4, son1 = 5;
522  if (split == 2) { son0 = 6; son1 = 7; }
523 
524  for (son = son0; son <= son1; son++)
525  {
526  State* ns = push_state();
527  move_to_son(&ns->cr, &s->cr, son);
528 
529  int j = (son == 4 || son == 6) ? 0 : 2;
530  for (i = 0; i < num; i++)
531  {
532  if (s->e[i] == nullptr || !s->e[i]->used)
533  {
534  ns->e[i] = nullptr;
535  }
536  else
537  {
538  if (s->e[i]->active)
539  {
540  ns->e[i] = s->e[i];
541  ns->sub_idx[i] = s->sub_idx[i];
542  ns->push_transform(son, i);
543  }
544  else
545  {
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)
549  ns->sub_idx[i] = 0;
550  }
551  }
552  }
553  }
554  }
555 
556  // No splits, recur to one son.
557  else
558  {
559  State* ns = push_state();
560  memcpy(&ns->cr, &s->cr, sizeof(Rect));
561 
562  for (i = 0; i < num; i++)
563  {
564  if (s->e[i] == nullptr || !s->e[i]->used)
565  {
566  ns->e[i] = nullptr;
567  }
568  else if (s->e[i]->active)
569  {
570  ns->e[i] = s->e[i];
571  memcpy(&ns->er[i], &ns->cr, sizeof(Rect));
572  }
573  else
574  {
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)
578  ns->sub_idx[i] = 0;
579  }
580  }
581  }
582  delete[] current_sons;
583  }
584  }
585  this->finish();
586  }
587 
588  static void testMeshesCompliance(int n, MeshSharedPtr* meshes)
589  {
590  // Test whether all master meshes have the same number of elements.
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().");
595  }
596 
597  static void testMeshesQuality(int n, MeshSharedPtr* meshes)
598  {
599  Element* e;
600  int base_elem_num = meshes[0]->get_num_base_elements();
601  // Test whether areas of corresponding elements are the same.
602  double *areas = new double[base_elem_num];
603  memset(areas, 0, base_elem_num*sizeof(double));
604 
605  // Read base element areas from the first mesh,
606  // Also get minimum element area.
607  int counter = 0;
608  double min_elem_area = 1e30;
609  for_all_base_elements_incl_inactive(e, meshes[0])
610  {
611  if (!e->used)
612  areas[counter] = 0.0;
613  else
614  {
615  areas[counter] = e->area;
616  if (areas[counter] < min_elem_area)
617  min_elem_area = areas[counter];
618  }
619 
620  counter++;
621  }
622  // take one mesh at a time and compare element areas to the areas[] array
623  double tolerance = min_elem_area / 100.;
624 
625  if (min_elem_area < 0)
626  throw Exceptions::ValueException("min_elem_area", 0.0, Hermes::HermesSqrtEpsilon);
627 
628  for (int i = 1; i < n; i++)
629  {
630  counter = 0;
631  for_all_base_elements_incl_inactive(e, meshes[i])
632  {
633  if (e->used)
634  if (fabs(areas[counter] - e->area) > tolerance && areas[counter] > Hermes::HermesSqrtEpsilon)
635  {
636  throw Hermes::Exceptions::Exception("An element is probably too distorted, try different meshing.");
637  }
638  counter++;
639  }
640  }
641  delete[] areas;
642  }
643 
644  void Traverse::begin(int n)
645  {
646  assert(n > 0);
647  num = n;
648 
649  size = 256;
650 
651  stack = new State[size];
652  memset(stack, 0, size * sizeof(State));
653 
654  top = 0;
655  }
656 
657  void Traverse::free_state(Traverse::State* state)
658  {
659  delete[] state->e;
660  delete[] state->er;
661  delete[] state->sub_idx;
662  memset(state, 0, sizeof(Traverse::State));
663  }
664 
665  void Traverse::finish()
666  {
667  if (stack == nullptr)
668  return;
669 
670  for (int i = 0; i < size; i++)
671  if (stack[i].e != nullptr)
672  free_state(stack + i);
673 
674  delete[] stack;
675  stack = nullptr;
676  }
677 
678  uint64_t Traverse::init_idx(Rect* cr, Rect* er)
679  {
680  Rect r;
681  memcpy(&r, er, sizeof(Rect));
682 
683  uint64_t idx = 0;
684  while (cr->l > r.l || cr->r < r.r || cr->b > r.b || cr->t < r.t)
685  {
686  uint64_t hmid = (r.l + r.r) >> 1;
687  uint64_t vmid = (r.t + r.b) >> 1;
688  int son = -1;
689 
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; }
698  else assert(0);
699 
700  idx = (idx << 3) + son + 1;
701  }
702  return idx;
703  }
704 
705  void Traverse::union_recurrent(Rect* cr, Element** e, Rect* er, uint64_t* idx, Element* uni)
706  {
707  int i, j, son;
708 
709  // are we at the bottom?
710  bool leaf = true;
711  for (i = 0; i < num; i++)
712  {
713  if (!e[i]->active)
714  {
715  leaf = false;
716  break;
717  }
718  }
719  // if yes, store the element transformation indices
720  if (leaf)
721  {
722  if (udsize <= uni->id)
723  {
724  if (!udsize) udsize = 1024;
725  while (udsize <= uni->id)
726  udsize *= 2;
727  for (i = 0; i < num; i++)
728  unidata[i] = (UniData*)realloc(unidata[i], udsize * sizeof(UniData));
729  }
730  for (i = 0; i < num; i++)
731  {
732  unidata[i][uni->id].e = e[i];
733  unidata[i][uni->id].idx = idx[i];
734  }
735  return;
736  }
737 
738  // state arrays
739  Element** e_new = new Element*[num];
740  Rect* er_new = new Rect[num];
741  Rect cr_new;
742 
743  int4* sons = new int4[num];
744  uint64_t* idx_new = new uint64_t[num];
745  memcpy(idx_new, idx, num*sizeof(uint64_t));
746 
747  if (uni->is_triangle())
748  {
749  // visit all sons of the triangle
750  unimesh->refine_element_id(uni->id);
751  for (son = 0; son <= 3; son++)
752  {
753  for (i = 0; i < num; i++)
754  {
755  if (e[i]->active)
756  {
757  e_new[i] = e[i];
758  idx_new[i] = (idx[i] << 3) + son + 1;
759  }
760  else
761  e_new[i] = e[i]->sons[son];
762  }
763  union_recurrent(nullptr, e_new, nullptr, idx_new, uni->sons[son]);
764  }
765  }
766  else
767  {
768  // obtain split types and son numbers for the current rectangle on all elements
769  int split = 0;
770  for (i = 0; i < num; i++)
771  if (!e[i]->active)
772  split |= get_split_and_sons(e[i], cr, er + i, sons[i]);
773 
774  // both splits: recur to four sons
775  if (split == 3)
776  {
777  unimesh->refine_element_id(uni->id, 0);
778 
779  for (son = 0; son <= 3; son++)
780  {
781  move_to_son(&cr_new, cr, son);
782  for (i = 0; i < num; i++)
783  {
784  if (e[i]->active)
785  {
786  e_new[i] = e[i];
787  idx_new[i] = (idx[i] << 3) + son + 1;
788  }
789  else
790  {
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]));
795  }
796  }
797  union_recurrent(&cr_new, e_new, er_new, idx_new, uni->sons[son]);
798  }
799  }
800  // v or h split, recur to two sons
801  else if (split > 0)
802  {
803  unimesh->refine_element_id(uni->id, split);
804 
805  int son0 = 4, son1 = 5;
806  if (split == 2) { son0 = 6; son1 = 7; }
807 
808  for (son = son0; son <= son1; son++)
809  {
810  move_to_son(&cr_new, cr, son);
811  j = (son == 4 || son == 6) ? 0 : 2;
812  for (i = 0; i < num; i++)
813  {
814  if (e[i]->active)
815  {
816  e_new[i] = e[i];
817  idx_new[i] = (idx[i] << 3) + son + 1;
818  }
819  else
820  {
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]));
825  }
826  }
827  union_recurrent(&cr_new, e_new, er_new, idx_new, uni->sons[son & 3]);
828  }
829  }
830  // no splits, recur to one son
831  else
832  {
833  memcpy(&cr_new, cr, sizeof(Rect));
834  for (i = 0; i < num; i++)
835  {
836  if (e[i]->active)
837  e_new[i] = e[i];
838  else
839  {
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]));
844  }
845  }
846  union_recurrent(&cr_new, e_new, er_new, idx_new, uni);
847  }
848  }
849 
850  delete[] e_new;
851  delete[] er_new;
852  delete[] sons;
853  delete[] idx_new;
854  }
855 
856  UniData** Traverse::construct_union_mesh(unsigned char n, MeshSharedPtr* meshes, MeshSharedPtr unimesh)
857  {
858  // Initial check.
859  testMeshesCompliance(n, meshes);
860 
861  Traverse traverse(n);
862 
863  // Initialization.
864  traverse.begin(n);
865 
866  int i;
867  Element** e = new Element*[n];
868  Rect* er = new Rect[n];
869  Rect cr;
870 
871  traverse.unimesh = unimesh;
872  unimesh->copy_base(meshes[0]);
873 
874  // Unimesh initialization.
875  traverse.udsize = 0;
876  traverse.unidata = new UniData*[n];
877  memset(traverse.unidata, 0, sizeof(UniData*)* n);
878 
879  uint64_t* idx = new uint64_t[n];
880  memset(idx, 0, n*sizeof(uint64_t));
881 
882  // Calculation.
883  for (int id = 0; id < meshes[0]->get_num_base_elements(); id++)
884  {
885  if (!meshes[0]->get_element(id)->used)
886  continue;
887  for (i = 0; i < n; i++)
888  {
889  e[i] = meshes[i]->get_element(id);
890  static const Rect H2D_UNITY = { 0, 0, ONE, ONE };
891  cr = er[i] = H2D_UNITY;
892  }
893  traverse.union_recurrent(&cr, e, er, idx, unimesh->get_element(id));
894  }
895 
896  delete[] e;
897  delete[] er;
898  delete[] idx;
899 
900  traverse.finish();
901  return traverse.unidata;
902  }
903 
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);
906  }
907 }
::xsd::cxx::tree::id< char, ncname > id
C++ type corresponding to the ID XML Schema built-in type.
Definition: adapt.h:24
State ** get_states(std::vector< MeshSharedPtr > meshes, unsigned int &states_count)
Definition: traverse.cpp:292