Hermes2D  2.0
orderizer.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 "orderizer.h"
17 #include "space.h"
18 #include "refmap.h"
19 #include "linear_data.cpp"
20 
21 namespace Hermes
22 {
23  namespace Hermes2D
24  {
25  namespace Views
26  {
27  // vertices
28  static int* ord_np[2] = { num_vert_tri, num_vert_quad };
29  static double3* ord_tables_tri[2] = { vert_tri0, vert_tri1 };
30  static double3* ord_tables_quad[2] = { vert_quad0, vert_quad1 };
31  static double3** ord_tables[2] = { ord_tables_tri, ord_tables_quad };
32 
33  // triangles
34  static int* num_elem[2] = { num_elem_tri, num_elem_quad};
35  static int3* ord_elem_tri[2] = { elem_tri0, elem_tri1 };
36  static int3* ord_elem_quad[2] = { elem_quad0, elem_quad1 };
37  static int3** ord_elem[2] = { ord_elem_tri, ord_elem_quad };
38 
39  // edges
40  static int* num_edge[2] = { num_edge_tri, num_edge_quad};
41  static int3* ord_edge_tri[2] = { edge_tri0, edge_tri1 };
42  static int3* ord_edge_quad[2] = { edge_quad0, edge_quad1 };
43  static int3** ord_edge[2] = { ord_edge_tri, ord_edge_quad };
44 
45  static class Quad2DOrd : public Quad2D
46  {
47  public:
48 
49  Quad2DOrd()
50  {
51  max_order[0] = max_order[1] = 1;
52  num_tables[0] = num_tables[1] = 2;
53  tables = ord_tables;
54  np = ord_np;
55  };
56 
57  virtual void dummy_fn() {}
58  } quad_ord;
59 
60  Orderizer::Orderizer() : LinearizerBase()
61  {
62  verts = NULL;
63  this->label_size = 0;
64  ltext = NULL;
65  lvert = NULL;
66  lbox = NULL;
67  tris_orders = NULL;
68 
69  label_count = cl1 = cl2 = cl3 = 0;
70 
71  for (int i = 0, p = 0; i <= 10; i++)
72  {
73  for (int j = 0; j <= 10; j++)
74  {
75  assert((unsigned) p < sizeof(buffer)-5);
76  if(i == j)
77  sprintf(buffer + p, "%d", i);
78  else
79  sprintf(buffer + p, "%d|%d", i, j);
80  labels[i][j] = buffer + p;
81  p += strlen(buffer + p) + 1;
82  }
83  }
84  }
85 
86  int Orderizer::add_vertex()
87  {
88  if(this->vertex_count >= this->vertex_size)
89  {
90  this->vertex_size *= 2;
91  verts = (double3*) realloc(verts, sizeof(double3) * vertex_size);
92  this->info = (int4*) realloc(info, sizeof(int4) * vertex_size);
93  }
94  return this->vertex_count++;
95  }
96 
97  void Orderizer::make_vert(int & index, double x, double y, double val)
98  {
99  index = add_vertex();
100  verts[index][0] = x;
101  verts[index][1] = y;
102  verts[index][2] = val;
103  }
104 
105  template<typename Scalar>
106  void Orderizer::process_space(const Space<Scalar>* space)
107  {
108  // sanity check
109  if(space == NULL) throw Hermes::Exceptions::Exception("Space is NULL in Orderizer:process_space().");
110 
111  if(!space->is_up_to_date())
112  throw Hermes::Exceptions::Exception("The space is not up to date.");
113 
114  int type = 1;
115  label_count = 0;
116  vertex_count = 0;
117  triangle_count = 0;
118  edges_count = 0;
119 
120  // estimate the required number of vertices and triangles
121  Mesh* mesh = space->get_mesh();
122  if(mesh == NULL)
123  {
124  throw Hermes::Exceptions::Exception("Mesh is NULL in Orderizer:process_space().");
125  }
126  int nn = mesh->get_num_active_elements();
127  this->vertex_size = std::max(this->vertex_size, 77 * nn);
128  this->triangle_size = std::max(this->triangle_size, 64 * nn);
129  this->edges_size = std::max(this->edges_size, 16 * nn);
130  this->label_size = std::max(this->label_size, nn + 10);
131  cl1 = cl2 = cl3 = label_size;
132 
133  // reuse or allocate vertex, triangle and edge arrays
134  verts = (double3*) realloc(verts, sizeof(double3) * vertex_size);
135  this->tris = (int3*) realloc(this->tris, sizeof(int3) * this->triangle_size);
136  this->tri_markers = (int*) realloc(this->tri_markers, sizeof(int) * this->triangle_size);
137  this->edges = (int2*) realloc(this->edges, sizeof(int2) * this->edges_size);
138  this->edge_markers = (int*) realloc(this->edge_markers, sizeof(int) * this->edges_size);
139  tris_orders = (int*) realloc(tris_orders, sizeof(int) * triangle_size);
140  info = NULL;
141  this->empty = false;
142  lvert = (int*) realloc(lvert, sizeof(int) * label_size);
143  ltext = (char**) realloc(ltext, sizeof(char*) * label_size);
144  lbox = (double2*) realloc(lbox, sizeof(double2) * label_size);
145 
146  int oo, o[6];
147 
148  RefMap refmap;
149  refmap.set_quad_2d(&quad_ord);
150 
151  // make a mesh illustrating the distribution of polynomial orders over the space
152  Element* e;
153  for_all_active_elements(e, mesh)
154  {
155  oo = o[4] = o[5] = space->get_element_order(e->id);
156  for (unsigned int k = 0; k < e->get_nvert(); k++)
157  o[k] = space->get_edge_order(e, k);
158 
159  refmap.set_active_element(e);
160  double* x = refmap.get_phys_x(type);
161  double* y = refmap.get_phys_y(type);
162 
163  double3* pt = quad_ord.get_points(type, e->get_mode());
164  int np = quad_ord.get_num_points(type, e->get_mode());
165  int id[80];
166  assert(np <= 80);
167 
168  int mode = e->get_mode();
169  if(e->is_quad())
170  {
171  o[4] = H2D_GET_H_ORDER(oo);
172  o[5] = H2D_GET_V_ORDER(oo);
173  }
174  make_vert(lvert[label_count], x[0], y[0], o[4]);
175 
176  for (int i = 1; i < np; i++)
177  make_vert(id[i-1], x[i], y[i], o[(int) pt[i][2]]);
178 
179  for (int i = 0; i < num_elem[mode][type]; i++)
180  this->add_triangle(id[ord_elem[mode][type][i][0]], id[ord_elem[mode][type][i][1]], id[ord_elem[mode][type][i][2]], o[4], e->marker);
181 
182  for (int i = 0; i < num_edge[mode][type]; i++)
183  {
184  if(e->en[ord_edge[mode][type][i][2]]->bnd || (y[ord_edge[mode][type][i][0] + 1] < y[ord_edge[mode][type][i][1] + 1]) ||
185  ((y[ord_edge[mode][type][i][0] + 1] == y[ord_edge[mode][type][i][1] + 1]) &&
186  (x[ord_edge[mode][type][i][0] + 1] < x[ord_edge[mode][type][i][1] + 1])))
187  {
188  LinearizerBase::add_edge(id[ord_edge[mode][type][i][0]], id[ord_edge[mode][type][i][1]], e->en[ord_edge[mode][type][i][2]]->marker);
189  }
190  }
191 
192  double xmin = 1e100, ymin = 1e100, xmax = -1e100, ymax = -1e100;
193  for (unsigned int k = 0; k < e->get_nvert(); k++)
194  {
195  if(e->vn[k]->x < xmin) xmin = e->vn[k]->x;
196  if(e->vn[k]->x > xmax) xmax = e->vn[k]->x;
197  if(e->vn[k]->y < ymin) ymin = e->vn[k]->y;
198  if(e->vn[k]->y > ymax) ymax = e->vn[k]->y;
199  }
200  lbox[label_count][0] = xmax - xmin;
201  lbox[label_count][1] = ymax - ymin;
202  ltext[label_count++] = labels[o[4]][o[5]];
203  }
204 
205  refmap.set_quad_2d(&g_quad_2d_std);
206  }
207 
208  void Orderizer::add_triangle(int iv0, int iv1, int iv2, int order, int marker)
209  {
210  int index;
211 #pragma omp critical(realloc_triangles)
212  {
213  if(this->del_slot >= 0) // reuse a slot after a deleted triangle
214  {
215  index = this->del_slot;
216  del_slot = -1;
217  }
218  {
219  if(triangle_count >= triangle_size)
220  {
221  tri_markers = (int*) realloc(tri_markers, sizeof(int) * (triangle_size * 2));
222  tris = (int3*) realloc(tris, sizeof(int3) * (triangle_size * 2));
223  tris_orders = (int*) realloc(tris_orders, sizeof(int) * (triangle_size = triangle_size * 2));
224  }
225  index = triangle_count++;
226 
227  tris[index][0] = iv0;
228  tris[index][1] = iv1;
229  tris[index][2] = iv2;
230  tris_orders[index] = order;
231  tri_markers[index] = marker;
232  }
233  }
234  }
235 
236  void Orderizer::free()
237  {
238  if(verts != NULL)
239  {
240  ::free(verts);
241  verts = NULL;
242  }
243  if(lvert != NULL)
244  {
245  ::free(lvert);
246  lvert = NULL;
247  }
248  if(ltext != NULL)
249  {
250  ::free(ltext);
251  ltext = NULL;
252  }
253  if(lbox != NULL)
254  {
255  ::free(lbox);
256  lbox = NULL;
257  }
258  if(tris_orders != NULL)
259  {
260  ::free(tris_orders);
261  tris_orders = NULL;
262  }
263 
264  LinearizerBase::free();
265  }
266 
267  Orderizer::~Orderizer()
268  {
269  free();
270  }
271 
272  template<typename Scalar>
273  void Orderizer::save_orders_vtk(const Space<Scalar>* space, const char* file_name)
274  {
275  process_space(space);
276 
277  FILE* f = fopen(file_name, "wb");
278  if(f == NULL) throw Hermes::Exceptions::Exception("Could not open %s for writing.", file_name);
279  lock_data();
280 
281  // Output header for vertices.
282  fprintf(f, "# vtk DataFile Version 2.0\n");
283  fprintf(f, "\n");
284  fprintf(f, "ASCII\n\n");
285  fprintf(f, "DATASET UNSTRUCTURED_GRID\n");
286 
287  // Output vertices.
288  fprintf(f, "POINTS %d %s\n", this->vertex_count, "float");
289  for (int i = 0; i < this->vertex_count; i++)
290  {
291  fprintf(f, "%g %g %g\n", this->verts[i][0], this->verts[i][1], 0.0);
292  }
293 
294  // Output elements.
295  fprintf(f, "\n");
296  fprintf(f, "CELLS %d %d\n", this->triangle_count, 4 * this->triangle_count);
297  for (int i = 0; i < this->triangle_count; i++)
298  {
299  fprintf(f, "3 %d %d %d\n", this->tris[i][0], this->tris[i][1], this->tris[i][2]);
300  }
301 
302  // Output cell types.
303  fprintf(f, "\n");
304  fprintf(f, "CELL_TYPES %d\n", this->triangle_count);
305  for (int i = 0; i < this->triangle_count; i++)
306  {
307  fprintf(f, "5\n"); // The "5" means triangle in VTK.
308  }
309 
310  // This outputs double solution values. Look into Hermes2D/src/output/vtk.cpp
311  // for how it is done for vectors.
312  fprintf(f, "\n");
313  fprintf(f, "CELL_DATA %d\n", this->triangle_count);
314  fprintf(f, "SCALARS %s %s %d\n", "Mesh", "float", 1);
315  fprintf(f, "LOOKUP_TABLE %s\n", "default");
316  for (int i = 0; i < this->triangle_count; i++)
317  {
318  fprintf(f, "%i\n", this->tris_orders[i]);
319  }
320 
321  unlock_data();
322  fclose(f);
323  }
324 
325  template<typename Scalar>
326  void Orderizer::save_mesh_vtk(const Space<Scalar>* space, const char* file_name)
327  {
328  process_space(space);
329 
330  FILE* f = fopen(file_name, "wb");
331  if(f == NULL) throw Hermes::Exceptions::Exception("Could not open %s for writing.", file_name);
332  lock_data();
333 
334  // Output header for vertices.
335  fprintf(f, "# vtk DataFile Version 2.0\n");
336  fprintf(f, "\n");
337  fprintf(f, "ASCII\n\n");
338  fprintf(f, "DATASET UNSTRUCTURED_GRID\n");
339 
340  // Output vertices.
341  fprintf(f, "POINTS %d %s\n", this->vertex_count, "float");
342  for (int i = 0; i < this->vertex_count; i++)
343  fprintf(f, "%g %g %g\n", this->verts[i][0], this->verts[i][1], 0.0);
344 
345  // Output elements.
346  fprintf(f, "\n");
347  fprintf(f, "CELLS %d %d\n", this->edges_count, + 3 * this->edges_count);
348  for (int i = 0; i < this->edges_count; i++)
349  fprintf(f, "2 %d %d\n", this->edges[i][0], this->edges[i][1]);
350 
351  // Output cell types.
352  fprintf(f, "\n");
353  fprintf(f, "CELL_TYPES %d\n", this->edges_count);
354 
355  for (int i = 0; i < this->edges_count; i++)
356  fprintf(f, "3\n"); // The "5" means triangle in VTK.
357 
358  // This outputs double solution values. Look into Hermes2D/src/output/vtk.cpp
359  // for how it is done for vectors.
360  fprintf(f, "\n");
361  fprintf(f, "CELL_DATA %d\n", this->edges_count);
362  fprintf(f, "SCALARS %s %s %d\n", "Mesh", "float", 1);
363  fprintf(f, "LOOKUP_TABLE %s\n", "default");
364  for (int i = 0; i < this->edges_count; i++)
365  fprintf(f, "0 \n");
366  unlock_data();
367  fclose(f);
368  }
369 
370  int Orderizer::get_labels(int*& lvert, char**& ltext, double2*& lbox) const
371  {
372  lvert = this->lvert;
373  ltext = this->ltext;
374  lbox = this->lbox;
375  return label_count;
376  }
377 
378  void Orderizer::calc_vertices_aabb(double* min_x, double* max_x, double* min_y, double* max_y) const
379  {
380  if(verts == NULL)
381  throw Exceptions::Exception("Cannot calculate AABB from NULL vertices");
382  calc_aabb(&verts[0][0], &verts[0][1], sizeof(double3), vertex_count, min_x, max_x, min_y, max_y);
383  }
384 
385  double3* Orderizer::get_vertices()
386  {
387  return this->verts;
388  }
389  int Orderizer::get_num_vertices()
390  {
391  return this->vertex_count;
392  }
393 
394  template HERMES_API void Orderizer::save_orders_vtk<double>(const Space<double>* space, const char* file_name);
395  template HERMES_API void Orderizer::save_orders_vtk<std::complex<double> >(const Space<std::complex<double> >* space, const char* file_name);
396  template HERMES_API void Orderizer::save_mesh_vtk<double>(const Space<double>* space, const char* file_name);
397  template HERMES_API void Orderizer::save_mesh_vtk<std::complex<double> >(const Space<std::complex<double> >* space, const char* file_name);
398  template HERMES_API void Orderizer::process_space<double>(const Space<double>* space);
399  template HERMES_API void Orderizer::process_space<std::complex<double> >(const Space<std::complex<double> >* space);
400  }
401  }
402 }