Hermes2D  3.0
mesh_util.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 "refmap.h"
18 
19 namespace Hermes
20 {
21  namespace Hermes2D
22  {
23  void MeshUtil::assign_curve(Node* en, Curve* curve, int p1, int p2)
24  {
25  // assign the arc to the elements sharing the edge node
26  for (unsigned int node_i = 0; node_i < 2; node_i++)
27  {
28  Element* e = en->elem[node_i];
29  if (e == nullptr)
30  continue;
31 
32  if (e->cm == nullptr)
33  {
34  e->cm = new CurvMap;
35  e->cm->toplevel = true;
36  e->cm->order = 4;
37  }
38 
39  int idx = -1;
40  for (unsigned j = 0; j < e->get_nvert(); j++)
41  {
42  if (e->en[j] == en)
43  {
44  idx = j;
45  break;
46  }
47  }
48  assert(idx >= 0);
49 
50  if (e->vn[idx]->id == p1)
51  e->cm->curves[idx] = curve;
52  else
53  {
54  Curve* curve_rev = MeshUtil::reverse_curve(curve);
55  e->cm->curves[idx] = curve_rev;
56  }
57  }
58  }
59 
61  {
62  if (curve->type == ArcType)
63  {
64  Arc* toReturn = new Arc((Arc*)curve);
65  ((Arc*)toReturn)->angle = -((Arc*)toReturn)->angle;
66  for (int i = 0; i < toReturn->np; i++)
67  {
68  toReturn->pt[((Arc*)curve)->np - 1 - i][0] = ((Arc*)curve)->pt[i][0];
69  toReturn->pt[((Arc*)curve)->np - 1 - i][1] = ((Arc*)curve)->pt[i][1];
70  toReturn->pt[((Arc*)curve)->np - 1 - i][2] = ((Arc*)curve)->pt[i][2];
71  }
72 
73  return toReturn;
74  }
75  else
76  {
77  Nurbs* toReturn = new Nurbs((Nurbs*)curve);
78 
79  for (int i = 0; i < ((Nurbs*)curve)->np; i++)
80  {
81  toReturn->pt[((Nurbs*)curve)->np - 1 - i][0] = ((Nurbs*)curve)->pt[i][0];
82  toReturn->pt[((Nurbs*)curve)->np - 1 - i][1] = ((Nurbs*)curve)->pt[i][1];
83  toReturn->pt[((Nurbs*)curve)->np - 1 - i][2] = ((Nurbs*)curve)->pt[i][2];
84  }
85 
86  for (int i = 0; i < ((Nurbs*)curve)->nk; i++)
87  toReturn->kv[i] = ((Nurbs*)curve)->kv[i];
88  for (int i = ((Nurbs*)curve)->degree + 1; i < ((Nurbs*)curve)->nk - ((Nurbs*)curve)->degree - 1; i++)
89  toReturn->kv[((Nurbs*)curve)->nk - 1 - i] = 1.0 - ((Nurbs*)curve)->kv[i];
90 
91  return toReturn;
92  }
93  }
94 
96  {
97  while (!base->active) // we need to go down to an active element
98  {
99  int son1, son2;
100  get_edge_sons(base, edge, son1, son2);
101  base = base->sons[son1];
102  }
103  return base->en[edge];
104  }
105 
106  int MeshUtil::get_edge_sons(Element* e, int edge, int& son1, int& son2)
107  {
108  assert(!e->active);
109 
110  if (!e->is_triangle())
111  {
112  if (e->sons[2] == nullptr) // horz quad
113  {
114  if (edge == 0 || edge == 2) { son1 = edge >> 1; return 1; }
115  else if (edge == 1) { son1 = 0; son2 = 1; return 2; }
116  else { son1 = 1; son2 = 0; return 2; }
117  }
118  else if (e->sons[0] == nullptr) // vert quad
119  {
120  if (edge == 1 || edge == 3) { son1 = (edge == 1) ? 3 : 2; return 1; }
121  else if (edge == 0) { son1 = 2; son2 = 3; return 2; }
122  else { son1 = 3; son2 = 2; return 2; }
123  }
124  }
125 
126  // triangle or 4-son quad
127  son1 = edge;
128  son2 = e->next_vert(edge);
129  return 2;
130  }
131 
132  Arc* MeshUtil::load_arc(MeshSharedPtr mesh, int id, Node** en, int p1, int p2, double angle, bool skip_check)
133  {
134  Arc* curve = new Arc(angle);
135 
136  *en = mesh->peek_edge_node(p1, p2);
137 
138  if (*en == nullptr)
139  {
140  if (!skip_check)
141  throw Hermes::Exceptions::MeshLoadFailureException("Curve #%d: edge %d-%d does not exist.", id, p1, p2);
142  else
143  return nullptr;
144  }
145 
146  // edge endpoints control points.
147  curve->pt[0][0] = mesh->nodes[p1].x;
148  curve->pt[0][1] = mesh->nodes[p1].y;
149  curve->pt[0][2] = 1.0;
150  curve->pt[2][0] = mesh->nodes[p2].x;
151  curve->pt[2][1] = mesh->nodes[p2].y;
152  curve->pt[2][2] = 1.0;
153 
154  // read the arc angle
155  double a = (180.0 - angle) / 180.0 * M_PI;
156 
157  // generate one inner control point
158  double x = 1.0 / std::tan(a * 0.5);
159  curve->pt[1][0] = 0.5*((curve->pt[2][0] + curve->pt[0][0]) + (curve->pt[2][1] - curve->pt[0][1]) * x);
160  curve->pt[1][1] = 0.5*((curve->pt[2][1] + curve->pt[0][1]) - (curve->pt[2][0] - curve->pt[0][0]) * x);
161  curve->pt[1][2] = Hermes::cos((M_PI - a) * 0.5);
162 
163  return curve;
164  }
165 
166  MeshHashGrid::MeshHashGrid(Mesh* mesh) : mesh_seq(mesh->get_seq())
167  {
168  mesh->calc_bounding_box();
169 
170  // create grid
171  double interval_len_x = (mesh->top_right_x - mesh->bottom_left_x) / GRID_SIZE;
172  double interval_len_y = (mesh->top_right_y - mesh->bottom_left_y) / GRID_SIZE;
173 
174  intervals_x[0] = mesh->bottom_left_x;
175  intervals_y[0] = mesh->bottom_left_y;
176 
177  for (int i = 1; i < GRID_SIZE; i++)
178  {
179  intervals_x[i] = intervals_x[i - 1] + interval_len_x;
180  intervals_y[i] = intervals_y[i - 1] + interval_len_y;
181  }
182 
183  intervals_x[GRID_SIZE] = mesh->top_right_x;
184  intervals_y[GRID_SIZE] = mesh->top_right_y;
185 
186  for (int i = 0; i < GRID_SIZE; i++)
187  {
188  for (int j = 0; j < GRID_SIZE; j++)
189  {
190  m_grid[i][j] = new MeshHashGridElement(intervals_x[i], intervals_y[j], intervals_x[i + 1], intervals_y[j + 1]);
191  }
192  }
193 
194  // assign elements
195  Element *element;
196  double2 p1, p2;
197  int x_min, x_max, y_min, y_max;
198  for_all_active_elements(element, mesh)
199  {
200  elementBoundingBox(element, p1, p2);
201 
202  x_min = 0;
203  while (intervals_x[x_min + 1] < p1[0])
204  x_min++;
205 
206  x_max = GRID_SIZE - 1;
207  while (intervals_x[x_max] > p2[0])
208  x_max--;
209 
210  y_min = 0;
211  while (intervals_y[y_min + 1] < p1[1])
212  y_min++;
213 
214  y_max = GRID_SIZE - 1;
215  while (intervals_y[y_max] > p2[1])
216  y_max--;
217 
218  assert((x_max >= x_min) && (y_max >= y_min));
219 
220  for (int i = x_min; i <= x_max; i++)
221  {
222  for (int j = y_min; j <= y_max; j++)
223  {
224  assert(m_grid[i][j]->belongs(element));
225  m_grid[i][j]->insert(element);
226  }
227  }
228  }
229  }
230 
231  MeshHashGrid::~MeshHashGrid()
232  {
233  for (int i = 0; i < GRID_SIZE; i++)
234  {
235  for (int j = 0; j < GRID_SIZE; j++)
236  {
237  delete m_grid[i][j];
238  }
239  }
240  }
241 
242  MeshHashGridElement::MeshHashGridElement(double lower_left_x, double lower_left_y, double upper_right_x, double upper_right_y, int depth) : lower_left_x(lower_left_x), lower_left_y(lower_left_y), upper_right_x(upper_right_x), upper_right_y(upper_right_y), m_depth(depth), m_active(true), element_count(0)
243  {
244  this->elements = new Element*[MAX_ELEMENTS];
245  for (int i = 0; i < 2; i++)
246  for (int j = 0; j < 2; j++)
247  m_sons[i][j] = nullptr;
248  }
249 
250  MeshHashGridElement::~MeshHashGridElement()
251  {
252  if (elements)
253  delete[] elements;
254  for (int i = 0; i < 2; i++)
255  for (int j = 0; j < 2; j++)
256  if (m_sons[i][j])
257  delete m_sons[i][j];
258  }
259 
260  bool MeshHashGridElement::belongs(Element *element)
261  {
262  double2 p1, p2;
263  MeshHashGrid::elementBoundingBox(element, p1, p2);
264  return ((p1[0] <= upper_right_x) && (p2[0] >= lower_left_x) && (p1[1] <= upper_right_y) && (p2[1] >= lower_left_y));
265  }
266 
267  void MeshHashGridElement::insert(Element *element)
268  {
269  if (m_active)
270  {
271  elements[element_count++] = element;
272 
273  if ((element_count >= MAX_ELEMENTS) && (m_depth < MAX_DEPTH))
274  {
275  m_active = false;
276  double xx[3] = { lower_left_x, (lower_left_x + upper_right_x) / 2., upper_right_x };
277  double yy[3] = { lower_left_y, (lower_left_y + upper_right_y) / 2., upper_right_y };
278  for (int i = 0; i < 2; i++)
279  {
280  double x0 = xx[i];
281  double x1 = xx[i + 1];
282  for (int j = 0; j < 2; j++)
283  {
284  double y0 = yy[j];
285  double y1 = yy[j + 1];
286 
287  assert(m_sons[i][j] == nullptr);
288 
289  m_sons[i][j] = new MeshHashGridElement(x0, y0, x1, y1, m_depth + 1);
290 
291  for (int elem_i = 0; elem_i < this->element_count; elem_i++)
292  {
293  if (m_sons[i][j]->belongs(elements[elem_i]))
294  m_sons[i][j]->insert(elements[elem_i]);
295  }
296  }
297  }
298 
299  delete[] elements;
300  elements = nullptr;
301  }
302  }
303  else
304  {
305  for (int i = 0; i < 2; i++)
306  {
307  for (int j = 0; j < 2; j++)
308  {
309  assert(m_sons[i][j]);
310  if (m_sons[i][j]->belongs(element))
311  m_sons[i][j]->insert(element);
312  }
313  }
314  }
315  }
316 
317  bool MeshHashGridElement::belongs(double x, double y)
318  {
319  return (x >= lower_left_x) && (x <= upper_right_x) && (y <= lower_left_y) && (y >= upper_right_y);
320  }
321 
323  {
324  if (m_active)
325  {
326  for (int elem_i = 0; elem_i < this->element_count; elem_i++)
327  if (RefMap::is_element_on_physical_coordinates(elements[elem_i], x, y))
328  return elements[elem_i];
329 
330  return nullptr;
331  }
332  else
333  {
334  Element* element;
335  for (int i = 0; i < 2; i++)
336  {
337  for (int j = 0; j < 2; j++)
338  {
339  element = m_sons[i][j]->getElement(x, y);
340  if (element)
341  return element;
342  }
343  }
344  return nullptr;
345  }
346  }
347 
348  void MeshHashGrid::elementBoundingBox(Element *element, double2 &p1, double2 &p2)
349  {
350  p1[0] = p2[0] = element->vn[0]->x;
351  p1[1] = p2[1] = element->vn[0]->y;
352 
353  for (int i = 1; i < element->get_nvert(); i++)
354  {
355  double xx = element->vn[i]->x;
356  double yy = element->vn[i]->y;
357  if (xx > p2[0])
358  p2[0] = xx;
359  if (xx < p1[0])
360  p1[0] = xx;
361  if (yy > p2[1])
362  p2[1] = yy;
363  if (yy < p1[1])
364  p1[1] = yy;
365  }
366 
367  if (element->is_curved())
368  {
369  // todo: should be improved
370  double diameter_x = p2[0] - p1[0];
371  p2[0] += diameter_x;
372  p1[0] -= diameter_x;
373 
374  double diameter_y = p2[1] - p1[1];
375  p2[1] += diameter_y;
376  p1[1] -= diameter_y;
377  }
378  }
379 
380  Element* MeshHashGrid::getElement(double x, double y)
381  {
382  int i = 0;
383  while ((i < GRID_SIZE) && (intervals_x[i + 1] < x))
384  i++;
385 
386  int j = 0;
387  while ((j < GRID_SIZE) && (intervals_y[j + 1] < y))
388  j++;
389 
390  // this means that x or y is outside mesh, but it can hapen
391  if ((i >= GRID_SIZE) || (j >= GRID_SIZE))
392  return nullptr;
393  else
394  return m_grid[i][j]->getElement(x, y);
395  }
396 
397  int MeshHashGrid::get_mesh_seq() const
398  {
399  return this->mesh_seq;
400  }
401 
402  MarkerArea::MarkerArea(Mesh *mesh, int marker) : mesh_seq(mesh->get_seq())
403  {
404  area = 0;
405  Element* elem;
406  for_all_active_elements(elem, mesh)
407  {
408  if (elem->marker == marker)
409  {
410  elem->calc_area(true);
411  area += elem->area;
412  }
413  }
414  }
415 
416  int MarkerArea::get_mesh_seq() const
417  {
418  return this->mesh_seq;
419  }
420 
421  double MarkerArea::get_area() const
422  {
423  return this->area;
424  }
425  }
426 }
double x
vertex node coordinates
Definition: element.h:63
int id
node id number
Definition: element.h:48
Definition: adapt.h:24
unsigned char next_vert(unsigned char i) const
Helper functions to obtain the index of the next or previous vertex/edge.
Definition: element.h:202
static Arc * load_arc(MeshSharedPtr mesh, int id, Node **en, int p1, int p2, double angle, bool skip_check=false)
Definition: mesh_util.cpp:132
Stores one element of a mesh.
Definition: element.h:107
Represents a finite element mesh. Typical usage: MeshSharedPtr mesh; Hermes::Hermes2D::MeshReaderH2DX...
Definition: mesh.h:61
Hermes::Hermes2D::Element * getElement(double x, double y)
Return the Element.
Definition: mesh_util.cpp:322
Element * elem[2]
elements sharing the edge node
Definition: element.h:70
CurvMap * cm
curved mapping, nullptr if not curvilinear
Definition: element.h:188
static const unsigned char np
Arc has 3 control points.
Definition: curved.h:61
unsigned short order
current polynomial degree of the refmap approximation
Definition: curved.h:122
Represents one NURBS curve.
Definition: curved.h:77
Stores one node of a mesh.
Definition: element.h:45
static bool is_element_on_physical_coordinates(Element *e, double x, double y, double *x_reference=nullptr, double *y_reference=nullptr)
Definition: refmap.cpp:670
Node * en[H2D_MAX_NUMBER_EDGES]
edge node pointers
Definition: element.h:138
static Curve * reverse_curve(Curve *curve)
Definition: mesh_util.cpp:60
Element * sons[H2D_MAX_ELEMENT_SONS]
son elements (up to four)
Definition: element.h:140
bool active
0 = active, no sons; 1 = inactive (refined), has sons
Definition: element.h:114
double * kv
knot vector
Definition: curved.h:93
static Node * get_base_edge_node(Element *base, int edge)
For internal use.
Definition: mesh_util.cpp:95
double3 * pt
control points and their weights
Definition: curved.h:89
Node * vn[H2D_MAX_NUMBER_VERTICES]
vertex node pointers
Definition: element.h:134
static int get_edge_sons(Element *e, int edge, int &son1, int &son2)
For internal use.
Definition: mesh_util.cpp:106
static void assign_curve(Node *en, Curve *curve, int p1, int p2)
For internal use.
Definition: mesh_util.cpp:23