Hermes2D  3.0
element.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 "element.h"
17 #include "refmap.h"
18 #include "global.h"
19 #include "forms.h"
20 
21 namespace Hermes
22 {
23  namespace Hermes2D
24  {
26  {
27  assert(type == HERMES_TYPE_VERTEX);
28  return ref <= 3 && !bnd;
29  }
30 
31  void Node::ref_element(Element* e)
32  {
33  if (type == HERMES_TYPE_EDGE)
34  {
35  // store the element pointer in a free slot of 'elem'
36  if (elem[0] == nullptr)
37  elem[0] = e;
38  else
39  {
40  if (elem[1] == nullptr)
41  elem[1] = e;
42  else
43  throw Hermes::Exceptions::Exception("No free slot 'elem'");
44  }
45  }
46  ref++;
47  }
48 
49  void Node::unref_element(HashTable* ht, Element* e)
50  {
51  if (type == HERMES_TYPE_VERTEX)
52  {
53  if (!--ref) ht->remove_vertex_node(id);
54  }
55  else
56  {
57  // remove the element from the array 'elem'
58  if (elem[0] == e) elem[0] = nullptr;
59  else if (elem[1] == e) elem[1] = nullptr;
60 
61  if (!--ref) ht->remove_edge_node(id);
62  }
63  }
64 
66  {
67  for (unsigned int i = 0; i < nvert; i++)
68  {
69  vn[i]->ref_element();
70  en[i]->ref_element(this);
71  }
72  this->calc_area();
73  this->calc_diameter();
74  }
75 
76  bool Element::get_edge_orientation(int ie) const
77  {
78  return !(this->vn[ie]->id < this->vn[this->next_vert(ie)]->id);
79  }
80 
82  {
83  for (unsigned int i = 0; i < nvert; i++)
84  {
85  vn[i]->unref_element(ht);
86  en[i]->unref_element(ht, this);
87  }
88  this->calc_area();
89  this->calc_diameter();
90  }
91 
92  Element::Element() : visited(false), area(0.0), diameter(0.0), center_set(false)
93  {
94  };
95 
96  bool Element::hsplit() const
97  {
98  if (active)
99  return false;
100  return sons[0] != nullptr;
101  }
102 
103  bool Element::vsplit() const
104  {
105  if (active)
106  return false;
107  return sons[2] != nullptr;
108  }
109 
110  bool Element::bsplit() const
111  {
112  if (active)
113  return false;
114  return sons[0] != nullptr && sons[2] != nullptr;
115  }
116 
118  {
119  Element** elem = en[ie]->elem;
120  if (elem[0] == this)
121  return elem[1];
122  if (elem[1] == this)
123  return elem[0];
124  assert(0);
125  return nullptr;
126  }
127 
128  void Element::calc_area(bool precise_for_curvature)
129  {
130  // First some basic arithmetics.
131  double ax, ay, bx, by;
132  ax = vn[1]->x - vn[0]->x;
133  ay = vn[1]->y - vn[0]->y;
134  bx = vn[2]->x - vn[0]->x;
135  by = vn[2]->y - vn[0]->y;
136 
137  this->area = 0.5*(ax*by - ay*bx);
138  if (is_quad())
139  {
140  ax = bx; ay = by;
141  bx = vn[3]->x - vn[0]->x;
142  by = vn[3]->y - vn[0]->y;
143 
144  this->area = area + 0.5*(ax*by - ay*bx);
145  }
146 
147  // Either the basic approximation is fine.
148  if (!this->is_curved() || !precise_for_curvature)
149  return;
150  // Or we want to capture the curvature precisely.
151  else
152  {
153  // Utility data.
154  RefMap refmap_curv;
155  RefMap refmap_straight;
156  double3* tan;
157 
158  double x_center, y_center;
159  this->get_center(x_center, y_center);
160 
161  for (unsigned char isurf = 0; isurf < this->nvert; isurf++)
162  {
163  // 0 - prepare data structures.
164  int eo = g_quad_2d_std.get_edge_points(isurf, this->get_mode() == HERMES_MODE_TRIANGLE ? g_max_tri : g_max_quad, this->get_mode());
165  unsigned char np = g_quad_2d_std.get_num_points(eo, this->get_mode());
166  double* x_curv = new double[np];
167  double* y_curv = new double[np];
168  double* x_straight = new double[np];
169  double* y_straight = new double[np];
170 
171  // 1 - get the x,y coordinates for the curved element.
172  refmap_curv.set_active_element(this);
173  GeomSurf<double> geometry;
174  init_geom_surf_allocated(geometry, &refmap_curv, isurf, this->en[isurf]->marker, eo, tan);
175  memcpy(x_curv, geometry.x, np*sizeof(double));
176  memcpy(y_curv, geometry.y, np*sizeof(double));
177 
178  // 2. - act if there was no curvature
179  CurvMap* cm_temp = this->cm;
180  this->cm = nullptr;
181  refmap_straight.set_active_element(this);
182  init_geom_surf_allocated(geometry, &refmap_straight, isurf, this->en[isurf]->marker, eo, tan);
183  memcpy(x_straight, geometry.x, np*sizeof(double));
184  memcpy(y_straight, geometry.y, np*sizeof(double));
185 
186  // 3. - compare the two, get the updated area.
187  double previous_distance;
188  for (int i = 0; i < np; i++)
189  {
190  // Distance between the curved and straight edges.
191  double distance_i = std::sqrt(std::pow(x_straight[i] - x_curv[i], 2.0) + std::pow(y_straight[i] - y_curv[i], 2.0));
192 
193  // Add to- or Subtract from- the area (depends on the curvature and we shall decide based on distance from the element center).
194  double distance_from_center_curved = std::pow(x_center - x_curv[i], 2.0) + std::pow(y_center - y_curv[i], 2.0);
195  double distance_from_center_straight = std::pow(x_center - x_straight[i], 2.0) + std::pow(y_center - y_straight[i], 2.0);
196  bool add = distance_from_center_curved > distance_from_center_straight;
197 
198  // Calculate now the area delta.
199  // It depends on the integration point number etc.
200  double area_delta;
201  if (i == 0)
202  {
203  double distance_along_edge = std::sqrt(std::pow(x_straight[i] - this->vn[isurf]->x, 2.0) + std::pow(y_straight[i] - this->vn[isurf]->y, 2.0));
204  area_delta = distance_i * distance_along_edge * 0.5;
205  }
206  if (i > 0 && i < np - 1)
207  {
208  double distance_along_edge = std::sqrt(std::pow(x_straight[i] - x_straight[i - 1], 2.0) + std::pow(y_straight[i] - y_straight[i - 1], 2.0));
209  area_delta = 0.5*(distance_i + previous_distance) * distance_along_edge;
210  }
211  if (i == np - 1)
212  {
213  double distance_along_edge = std::sqrt(std::pow(x_straight[i] - this->vn[(isurf + 1) % this->nvert]->x, 2.0) + std::pow(y_straight[i] - this->vn[(isurf + 1) % this->nvert]->y, 2.0));
214  area_delta = distance_i * distance_along_edge * 0.5;
215  }
216 
217  if (add)
218  area += area_delta;
219  else
220  area -= area_delta;
221 
222  previous_distance = distance_i;
223  }
224 
225  // 4. - re-add the curvature.
226  this->cm = cm_temp;
227 
228  // clean up
229  delete[] x_curv;
230  delete[] y_curv;
231  delete[] x_straight;
232  delete[] y_straight;
233  }
234  }
235  }
236 
237  void Element::get_center(double& x, double& y)
238  {
239  if (!this->center_set)
240  {
241 #pragma omp critical
242  {
243  if (!this->center_set)
244  {
245  // x - coordinate
246  this->x_center = this->vn[0]->x + this->vn[1]->x + this->vn[2]->x;
247  this->y_center = this->vn[0]->y + this->vn[1]->y + this->vn[2]->y;
248  if (this->is_quad())
249  {
250  this->x_center += this->vn[3]->x;
251  this->x_center = this->x_center / 4.0;
252  this->y_center += this->vn[3]->y;
253  this->y_center = this->y_center / 4.0;
254  }
255  else
256  {
257  this->x_center = this->x_center / 3.0;
258  this->y_center = this->y_center / 3.0;
259  }
260  x = this->x_center;
261  y = this->y_center;
262  this->center_set = true;
263  }
264  }
265  }
266  else
267  {
268  x = this->x_center;
269  y = this->y_center;
270  return;
271  }
272  }
273 
275  {
276  double max, l;
277  if (is_triangle())
278  {
279  max = 0.0;
280  for (int i = 0; i < 3; i++)
281  {
282  int j = next_vert(i);
283  l = sqr(vn[i]->x - vn[j]->x) + sqr(vn[i]->y - vn[j]->y);
284  if (l > max)
285  max = l;
286  }
287  }
288  else
289  {
290  max = sqr(vn[0]->x - vn[2]->x) + sqr(vn[0]->y - vn[2]->y);
291  l = sqr(vn[1]->x - vn[3]->x) + sqr(vn[1]->y - vn[3]->y);
292  if (l > max)
293  max = l;
294  }
295  this->diameter = sqrt(max);
296  }
297 
298  Node* get_vertex_node(Node* v1, Node* v2)
299  {
300  // initialize the new_ Node
301  Node* newnode = new Node();
302  newnode->type = HERMES_TYPE_VERTEX;
303  newnode->ref = 0;
304  newnode->bnd = 0;
305  newnode->p1 = -9999;
306  newnode->p2 = -9999;
307  newnode->x = (v1->x + v2->x) * 0.5;
308  newnode->y = (v1->y + v2->y) * 0.5;
309 
310  return newnode;
311  }
312 
313  Node* get_edge_node()
314  {
315  // initialize the new_ Node
316  Node* newnode = new Node();
317  newnode->type = HERMES_TYPE_EDGE;
318  newnode->ref = 0;
319  newnode->bnd = 0;
320  newnode->p1 = -9999;
321  newnode->p2 = -9999;
322  newnode->marker = 0;
323  newnode->elem[0] = newnode->elem[1] = nullptr;
324 
325  return newnode;
326  }
327  }
328 }
double x
vertex node coordinates
Definition: element.h:63
unsigned type
0 = vertex node; 1 = edge node
Definition: element.h:52
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
int marker
element marker
Definition: element.h:144
HERMES_API void init_geom_surf_allocated(GeomSurf< double > &geom, RefMap *rm, unsigned char isurf, int marker, const int order, double3 *&tan)
Init element geometry for surface integrals.
Definition: forms.cpp:444
Stores one element of a mesh.
Definition: element.h:107
T y[H2D_MAX_INTEGRATION_POINTS_COUNT]
y-coordinates[in physical domain].
Definition: forms.h:46
T x[H2D_MAX_INTEGRATION_POINTS_COUNT]
x-coordinates[in physical domain].
Definition: forms.h:44
Element * elem[2]
elements sharing the edge node
Definition: element.h:70
CurvMap * cm
curved mapping, nullptr if not curvilinear
Definition: element.h:188
Common definitions for Hermes2D.
Stores one node of a mesh.
Definition: element.h:45
unsigned char nvert
number of vertices (3 or 4)
Definition: element.h:222
void ref_all_nodes()
Internal.
Definition: element.cpp:65
void unref_all_nodes(HashTable *ht)
Internal.
Definition: element.cpp:81
double area
Serves for saving the once calculated area of this element.
Definition: element.h:190
Node * en[H2D_MAX_NUMBER_EDGES]
edge node pointers
Definition: element.h:138
int p1
parent id numbers
Definition: element.h:75
unsigned bnd
1 = boundary node; 0 = inner node
Definition: element.h:54
Element * get_neighbor(int ie) const
Definition: element.cpp:117
void get_center(double &x, double &y)
Returns the center of gravity.
Definition: element.cpp:237
bool is_constrained_vertex() const
Returns true if the (vertex) node is constrained.
Definition: element.cpp:25
Represents the reference mapping.
Definition: refmap.h:40
unsigned ref
the number of elements using the node
Definition: element.h:50
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 diameter
Serves for saving the once calculated diameter of this element.
Definition: element.h:196
void calc_diameter()
Calculates the diameter.
Definition: element.cpp:274
Stores and searches node tables.
Definition: hash.h:39
void calc_area(bool precise_for_curvature=false)
This takes much longer.
Definition: element.cpp:128
virtual void set_active_element(Element *e)
Definition: refmap.cpp:70
Node * vn[H2D_MAX_NUMBER_VERTICES]
vertex node pointers
Definition: element.h:134