Hermes2D  3.0
multimesh_dg_neighbor_tree.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 "discrete_problem/dg/multimesh_dg_neighbor_tree.h"
17 
18 namespace Hermes
19 {
20  namespace Hermes2D
21  {
22  template<typename Scalar>
23  void MultimeshDGNeighborTree<Scalar>::process_edge(NeighborSearch<Scalar>** neighbor_searches, unsigned char num_neighbor_searches, unsigned int& num_neighbors, bool*& processed)
24  {
25  MultimeshDGNeighborTreeNode root(nullptr, 0);
26 
27  build_multimesh_tree(&root, neighbor_searches, num_neighbor_searches);
28 
29  // Update all NeighborSearches according to the multimesh tree.
30  // After this, all NeighborSearches in neighbor_searches should have the same count
31  // of neighbors and proper set of transformations
32  // for the central and the neighbor element(s) alike.
33  // Also check that every NeighborSearch has the same number o f neighbor elements.
34  num_neighbors = 0;
35  for (unsigned int i = 0; i < num_neighbor_searches; i++)
36  {
37  NeighborSearch<Scalar>* ns = neighbor_searches[i];
38  update_neighbor_search(ns, &root);
39 
40  if (num_neighbors == 0)
41  num_neighbors = ns->n_neighbors;
42 
43  if (ns->n_neighbors != num_neighbors)
44  throw Hermes::Exceptions::Exception("Num_neighbors of different NeighborSearches not matching in assemble_one_state().");
45  }
46 
47  processed = new bool[num_neighbors];
48 
49  for (unsigned int neighbor_i = 0; neighbor_i < num_neighbors; neighbor_i++)
50  {
51  // If the active segment has already been processed (when the neighbor element was assembled), it is skipped.
52  // We test all neighbor searches, because in the case of intra-element edge, the neighboring (the same as central) element
53  // will be marked as visited, even though the edge was not calculated.
54  processed[neighbor_i] = true;
55  for (unsigned int i = 0; i < num_neighbor_searches; i++)
56  {
57  if (!neighbor_searches[i]->neighbors.at(neighbor_i)->visited)
58  {
59  processed[neighbor_i] = false;
60  break;
61  }
62  }
63  }
64  }
65 
66  template<typename Scalar>
68  {
69  for (unsigned int i = 0; i < number; i++)
70  {
71  NeighborSearch<Scalar>* ns = neighbor_searches[i];
72  if (ns->n_neighbors == 1 && (ns->central_transformations_size == 0 || ns->central_transformations[0]->num_levels == 0))
73  continue;
74  for (unsigned int j = 0; j < ns->n_neighbors; j++)
75  if (j < ns->central_transformations_alloc_size && ns->central_transformations[j])
76  insert_into_multimesh_tree(root, ns->central_transformations[j]->transf, ns->central_transformations[j]->num_levels);
77  }
78  }
79 
80  template<typename Scalar>
81  void MultimeshDGNeighborTree<Scalar>::insert_into_multimesh_tree(MultimeshDGNeighborTreeNode* node, unsigned int* transformations, unsigned int transformation_count)
82  {
83  // If we are already in the leaf.
84  if (transformation_count == 0)
85  return;
86  // Both sons are null. We have to add a new_ Node. Let us do it for the left sone of node.
87  if (node->get_left_son() == nullptr && node->get_right_son() == nullptr)
88  {
89  node->set_left_son(new MultimeshDGNeighborTreeNode(node, transformations[0]));
90  insert_into_multimesh_tree(node->get_left_son(), transformations + 1, transformation_count - 1);
91  }
92  // At least the left son is not null (it is impossible only for the right one to be not null, because
93  // the left one always gets into the tree first, as seen above).
94  else
95  {
96  // The existing left son is the right one to continue through.
97  if (node->get_left_son()->get_transformation() == transformations[0])
98  insert_into_multimesh_tree(node->get_left_son(), transformations + 1, transformation_count - 1);
99  // The right one also exists, check that it is the right one, or return an error.
100  else if (node->get_right_son())
101  {
102  if (node->get_right_son()->get_transformation() == transformations[0])
103  insert_into_multimesh_tree(node->get_right_son(), transformations + 1, transformation_count - 1);
104  else
105  throw Hermes::Exceptions::Exception("More than two possible sons in insert_into_multimesh_tree().");
106  }
107  // If the right one does not exist and the left one was not correct, create a right son and continue this way.
108  else
109  {
110  node->set_right_son(new MultimeshDGNeighborTreeNode(node, transformations[0]));
111  insert_into_multimesh_tree(node->get_right_son(), transformations + 1, transformation_count - 1);
112  }
113  }
114  }
115 
116  template<typename Scalar>
117  std::vector<std::vector<unsigned int>*> MultimeshDGNeighborTree<Scalar>::get_multimesh_neighbors_transformations(MultimeshDGNeighborTreeNode* multimesh_tree)
118  {
119  // Initialize the vector.
120  std::vector<std::vector<unsigned int>*> running_transformations;
121  // Prepare the first neighbor's vector.
122  running_transformations.push_back(new std::vector<unsigned int>);
123  // Fill the vector.
124  traverse_multimesh_tree(multimesh_tree, running_transformations);
125  return running_transformations;
126  }
127 
128  template<typename Scalar>
129  void MultimeshDGNeighborTree<Scalar>::traverse_multimesh_tree(MultimeshDGNeighborTreeNode* node,
130  std::vector<std::vector<unsigned int>*>& running_transformations)
131  {
132  // If we are in the root.
133  if (node->get_transformation() == 0)
134  {
135  if (node->get_left_son())
136  traverse_multimesh_tree(node->get_left_son(), running_transformations);
137  if (node->get_right_son())
138  traverse_multimesh_tree(node->get_right_son(), running_transformations);
139  // Delete the vector prepared by the last accessed leaf.
140  delete running_transformations.back();
141  running_transformations.pop_back();
142  return;
143  }
144  // If we are in a leaf.
145  if (node->get_left_son() == nullptr && node->get_right_son() == nullptr)
146  {
147  // Create a vector for the new_ neighbor.
148  std::vector<unsigned int>* new_neighbor_transformations = new std::vector < unsigned int > ;
149  // Copy there the whole path except for this leaf.
150  for (unsigned int i = 0; i < running_transformations.back()->size(); i++)
151  new_neighbor_transformations->push_back((*running_transformations.back())[i]);
152  // Insert this leaf into the current running transformation, thus complete it.
153  running_transformations.back()->push_back(node->get_transformation());
154  // Make the new_neighbor_transformations the current running transformation.
155  running_transformations.push_back(new_neighbor_transformations);
156  return;
157  }
158  else
159  {
160  running_transformations.back()->push_back(node->get_transformation());
161  if (node->get_left_son())
162  traverse_multimesh_tree(node->get_left_son(), running_transformations);
163  if (node->get_right_son())
164  traverse_multimesh_tree(node->get_right_son(), running_transformations);
165  running_transformations.back()->pop_back();
166  return;
167  }
168  return;
169  }
170 
171  template<typename Scalar>
172  void MultimeshDGNeighborTree<Scalar>::update_neighbor_search(NeighborSearch<Scalar>* ns, MultimeshDGNeighborTreeNode* multimesh_tree)
173  {
174  // This has to be done, because we pass ns by reference and the number of neighbors is changing.
175  unsigned int num_neighbors = ns->get_num_neighbors();
176 
177  for (int i = 0; i < num_neighbors; i++)
178  {
179  // Find the node corresponding to this neighbor in the tree.
180  MultimeshDGNeighborTreeNode* node;
181  if (i < ns->central_transformations_alloc_size && ns->central_transformations[i])
182  node = find_node(ns->central_transformations[i]->transf, ns->central_transformations[i]->num_levels, multimesh_tree);
183  else
184  node = multimesh_tree;
185 
186  // Update the NeighborSearch.
187  int added = update_ns_subtree(ns, node, i);
188  if (added >= 0)
189  {
190  i--;
191  num_neighbors--;
192  }
193  }
194  }
195 
196  template<typename Scalar>
197  MultimeshDGNeighborTreeNode* MultimeshDGNeighborTree<Scalar>::find_node(unsigned int* transformations,
198  unsigned int transformation_count,
199  MultimeshDGNeighborTreeNode* node)
200  {
201  // If there are no transformations left.
202  if (transformation_count == 0)
203  return node;
204  else
205  {
206  if (node->get_left_son())
207  {
208  if (node->get_left_son()->get_transformation() == transformations[0])
209  return find_node(transformations + 1, transformation_count - 1, node->get_left_son());
210  }
211  if (node->get_right_son())
212  {
213  if (node->get_right_son()->get_transformation() == transformations[0])
214  return find_node(transformations + 1, transformation_count - 1, node->get_right_son());
215  }
216  }
217  // We always should be able to empty the transformations array.
218  throw
219  Hermes::Exceptions::Exception("Transformation of a central element not found in the multimesh tree.");
220  return nullptr;
221  }
222 
223  template<typename Scalar>
224  int MultimeshDGNeighborTree<Scalar>::update_ns_subtree(NeighborSearch<Scalar>* ns,
225  MultimeshDGNeighborTreeNode* node, unsigned int ith_neighbor)
226  {
227  int current_count = ns->get_num_neighbors();
228 
229  // No subtree => no work.
230  // Also check the assertion that if one son is null, then the other too.
231  if (node->get_left_son() == nullptr)
232  {
233  if (node->get_right_son())
234  throw Hermes::Exceptions::Exception("Only one son (right) not null in MultimeshDGNeighborTree<Scalar>::update_ns_subtree.");
235  return -1;
236  }
237 
238  // Key part.
239  // Begin with storing the info about the current neighbor.
240  Element* neighbor = ns->neighbors[ith_neighbor];
241  typename NeighborSearch<Scalar>::NeighborEdgeInfo edge_info = ns->neighbor_edges[ith_neighbor];
242 
243  // Initialize the vector for central transformations->
244  std::vector<std::vector<unsigned int>*> running_central_transformations;
245  // Prepare the first new_ neighbor's vector. Push back the current transformations (in case of GO_DOWN neighborhood).
246  running_central_transformations.push_back(new std::vector<unsigned int>);
247  if (ith_neighbor < ns->central_transformations_alloc_size && ns->central_transformations[ith_neighbor])
248  ns->central_transformations[ith_neighbor]->copy_to(running_central_transformations.back());
249 
250  // Initialize the vector for neighbor transformations->
251  std::vector<std::vector<unsigned int>*> running_neighbor_transformations;
252  // Prepare the first new_ neighbor's vector. Push back the current transformations (in case of GO_UP/NO_TRF neighborhood).
253  running_neighbor_transformations.push_back(new std::vector<unsigned int>);
254  if (ith_neighbor < ns->neighbor_transformations_alloc_size && ns->neighbor_transformations[ith_neighbor])
255  ns->neighbor_transformations[ith_neighbor]->copy_to(running_neighbor_transformations.back());
256 
257  // Delete the current neighbor.
258  ns->delete_neighbor(ith_neighbor);
259 
260  // Move down the subtree.
261  if (node->get_left_son())
262  traverse_multimesh_subtree(node->get_left_son(), running_central_transformations,
263  running_neighbor_transformations, edge_info, ns->active_edge,
264  ns->central_el->get_mode());
265  if (node->get_right_son())
266  traverse_multimesh_subtree(node->get_right_son(), running_central_transformations,
267  running_neighbor_transformations, edge_info, ns->active_edge,
268  ns->central_el->get_mode());
269 
270  // Delete the last neighbors' info (this is a dead end, caused by the function traverse_multimesh_subtree.
271  delete running_central_transformations.back();
272  running_central_transformations.pop_back();
273  delete running_neighbor_transformations.back();
274  running_neighbor_transformations.pop_back();
275 
276  // Insert new_ neighbors.
277  for (unsigned int i = 0; i < running_central_transformations.size(); i++)
278  {
279  ns->neighbors.push_back(neighbor);
280  ns->neighbor_edges.push_back(edge_info);
281 
282  if ((ns->n_neighbors >= ns->central_transformations_alloc_size) || !ns->central_transformations[ns->n_neighbors])
283  ns->add_central_transformations(new typename NeighborSearch<Scalar>::Transformations, ns->n_neighbors);
284 
285  if ((ns->n_neighbors >= ns->neighbor_transformations_alloc_size) || !ns->neighbor_transformations[ns->n_neighbors])
286  ns->add_neighbor_transformations(new typename NeighborSearch<Scalar>::Transformations, ns->n_neighbors);
287 
288  ns->central_transformations[ns->n_neighbors]->copy_from(*running_central_transformations[i]);
289  ns->neighbor_transformations[ns->n_neighbors]->copy_from(*running_neighbor_transformations[i]);
290 
291  ns->n_neighbors++;
292  }
293 
294  for (unsigned int i = 0; i < running_central_transformations.size(); i++)
295  delete running_central_transformations[i];
296  for (unsigned int i = 0; i < running_neighbor_transformations.size(); i++)
297  delete running_neighbor_transformations[i];
298 
299  // Return the number of neighbors added/deleted.
300  return ns->get_num_neighbors() - current_count;
301  }
302 
303  template<typename Scalar>
304  void MultimeshDGNeighborTree<Scalar>::traverse_multimesh_subtree(MultimeshDGNeighborTreeNode* node,
305  std::vector<std::vector<unsigned int>*>& running_central_transformations,
306  std::vector<std::vector<unsigned int>*>& running_neighbor_transformations,
307  const typename NeighborSearch<Scalar>::NeighborEdgeInfo& edge_info, const int& active_edge, const int& mode)
308  {
309  // If we are in a leaf.
310  if (node->get_left_son() == nullptr && node->get_right_son() == nullptr)
311  {
312  // Create vectors for the new_ neighbor.
313  std::vector<unsigned int>* new_neighbor_central_transformations = new std::vector < unsigned int > ;
314  std::vector<unsigned int>* new_neighbor_neighbor_transformations = new std::vector < unsigned int > ;
315 
316  // Copy there the whole path except for this leaf.
317  for (unsigned int i = 0; i < running_central_transformations.back()->size(); i++)
318  new_neighbor_central_transformations->push_back((*running_central_transformations.back())[i]);
319  for (unsigned int i = 0; i < running_neighbor_transformations.back()->size(); i++)
320  new_neighbor_neighbor_transformations->push_back((*running_neighbor_transformations.back())[i]);
321 
322  // Insert this leaf into the current running central transformation, thus complete it.
323  running_central_transformations.back()->push_back(node->get_transformation());
324 
325  // Make the new_neighbor_central_transformations the current running central transformation.
326  running_central_transformations.push_back(new_neighbor_central_transformations);
327 
328  // Take care of the neighbor transformation.
329  // Insert appropriate info from this leaf into the current running neighbor transformation, thus complete it.
330  if (mode == HERMES_MODE_TRIANGLE)
331  if ((active_edge == 0 && node->get_transformation() == 0) ||
332  (active_edge == 1 && node->get_transformation() == 1) ||
333  (active_edge == 2 && node->get_transformation() == 2))
334  running_neighbor_transformations.back()->push_back((!edge_info.orientation ? edge_info.local_num_of_edge : (edge_info.local_num_of_edge + 1) % 3));
335  else
336  running_neighbor_transformations.back()->push_back((edge_info.orientation ? edge_info.local_num_of_edge : (edge_info.local_num_of_edge + 1) % 3));
337  // Quads.
338  else
339  if ((active_edge == 0 && (node->get_transformation() == 0 || node->get_transformation() == 6)) ||
340  (active_edge == 1 && (node->get_transformation() == 1 || node->get_transformation() == 4)) ||
341  (active_edge == 2 && (node->get_transformation() == 2 || node->get_transformation() == 7)) ||
342  (active_edge == 3 && (node->get_transformation() == 3 || node->get_transformation() == 5)))
343  running_neighbor_transformations.back()->push_back((!edge_info.orientation ? edge_info.local_num_of_edge : (edge_info.local_num_of_edge + 1) % H2D_MAX_NUMBER_EDGES));
344  else
345  running_neighbor_transformations.back()->push_back((edge_info.orientation ? edge_info.local_num_of_edge : (edge_info.local_num_of_edge + 1) % H2D_MAX_NUMBER_EDGES));
346 
347  // Make the new_neighbor_neighbor_transformations the current running neighbor transformation.
348  running_neighbor_transformations.push_back(new_neighbor_neighbor_transformations);
349 
350  return;
351  }
352  else
353  {
354  // Insert this leaf into the current running central transformation, thus complete it.
355  running_central_transformations.back()->push_back(node->get_transformation());
356 
357  // Insert appropriate info from this leaf into the current running neighbor transformation, thus complete it.
358  // Triangles.
359  if (mode == HERMES_MODE_TRIANGLE)
360  if ((active_edge == 0 && node->get_transformation() == 0) ||
361  (active_edge == 1 && node->get_transformation() == 1) ||
362  (active_edge == 2 && node->get_transformation() == 2))
363  running_neighbor_transformations.back()->push_back((!edge_info.orientation ? edge_info.local_num_of_edge : (edge_info.local_num_of_edge + 1) % 3));
364  else
365  running_neighbor_transformations.back()->push_back((edge_info.orientation ? edge_info.local_num_of_edge : (edge_info.local_num_of_edge + 1) % 3));
366  // Quads.
367  else
368  if ((active_edge == 0 && (node->get_transformation() == 0 || node->get_transformation() == 6)) ||
369  (active_edge == 1 && (node->get_transformation() == 1 || node->get_transformation() == 4)) ||
370  (active_edge == 2 && (node->get_transformation() == 2 || node->get_transformation() == 7)) ||
371  (active_edge == 3 && (node->get_transformation() == 3 || node->get_transformation() == 5)))
372  running_neighbor_transformations.back()->push_back((!edge_info.orientation ? edge_info.local_num_of_edge : (edge_info.local_num_of_edge + 1) % H2D_MAX_NUMBER_EDGES));
373  else
374  running_neighbor_transformations.back()->push_back((edge_info.orientation ? edge_info.local_num_of_edge : (edge_info.local_num_of_edge + 1) % H2D_MAX_NUMBER_EDGES));
375 
376  // Move down.
377  if (node->get_left_son())
378  traverse_multimesh_subtree(node->get_left_son(), running_central_transformations, running_neighbor_transformations,
379  edge_info, active_edge, mode);
380  if (node->get_right_son())
381  traverse_multimesh_subtree(node->get_right_son(), running_central_transformations, running_neighbor_transformations,
382  edge_info, active_edge, mode);
383 
384  // Take this transformation out.
385  running_central_transformations.back()->pop_back();
386  running_neighbor_transformations.back()->pop_back();
387  return;
388  }
389  return;
390  }
391 
392  template class HERMES_API MultimeshDGNeighborTree < double > ;
393  template class HERMES_API MultimeshDGNeighborTree < std::complex<double> > ;
394  }
395 }
Definition: adapt.h:24
static void process_edge(NeighborSearch< Scalar > **neighbor_searches, unsigned char num_neighbor_searches, unsigned int &num_neighbors, bool *&processed)
The main method, for the passed neighbor searches, it will process all multi-mesh neighbor consolidat...
Transformations ** central_transformations
Array of transformations of the central element to each neighbor.
#define H2D_MAX_NUMBER_EDGES
A maximum number of edges of an element.
Definition: global.h:31
This class characterizes a neighborhood of a given edge in terms of adjacent elements and provides me...