Hermes2D  2.0
mesh.h
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 #ifndef __H2D_MESH_H
17 #define __H2D_MESH_H
18 
19 #include "../global.h"
20 #include "curved.h"
21 #include "hash.h"
22 #include "../mixins2d.h"
23 
24 namespace Hermes
25 {
26  namespace Hermes2D
27  {
28  enum
29  {
30  HERMES_TYPE_VERTEX = 0,
31  HERMES_TYPE_EDGE = 1
32  };
33 
34  class Element;
35  class HashTable;
36 
37  template<typename Scalar> class Space;
38  template<typename Scalar> class KellyTypeAdapt;
39  struct MItem;
40  struct Rect;
41  extern unsigned g_mesh_seq;
42 
43  namespace RefinementSelectors
44  {
45  template<typename Scalar> class Selector;
46  template<typename Scalar> class POnlySelector;
47  template<typename Scalar> class HOnlySelector;
48  template<typename Scalar> class OptimumSelector;
49  template<typename Scalar> class ProjBasedSelector;
50  template<typename Scalar> class H1ProjBasedSelector;
51  template<typename Scalar> class L2ProjBasedSelector;
52  template<typename Scalar> class HcurlProjBasedSelector;
53  }
54 
55  namespace Views
56  {
57  class MeshView;
58  }
59 
68  struct HERMES_API Node
69  {
70  int id;
71  unsigned ref:29;
72  unsigned type:1;
73  unsigned bnd:1;
74  unsigned used:1;
75 
76  union
77  {
78  struct
79  {
80  double x, y;
81  };
82  struct
83  {
84  int marker;
85  Element* elem[2];
86  };
87  };
88 
89  int p1, p2;
90  Node* next_hash;
91 
93  bool is_constrained_vertex() const;
94 
95  void ref_element(Element* e = NULL);
96  void unref_element(HashTable* ht, Element* e = NULL);
97  };
98 
120  class HERMES_API Element
121  {
122  public:
123  Element();
124  int id;
125  bool active;
126  bool used;
127  Element* parent;
128  bool visited;
129 
132  double get_area();
133 
136  double get_diameter();
137 
139  void get_center(double& x, double& y);
140 
141  Node* vn[H2D_MAX_NUMBER_VERTICES];
142  union
143  {
144  Node* en[H2D_MAX_NUMBER_EDGES];
145  Element* sons[H2D_MAX_ELEMENT_SONS];
146  };
147 
148  int marker;
149 
150  // returns the edge orientation. This works for the unconstrained edges.
151  int get_edge_orientation(int ie) const;
152  ElementMode2D get_mode() const;
153 
154  bool is_triangle() const;
155  bool is_quad() const;
156  bool is_curved() const;
157  int get_nvert() const;
158 
159  bool hsplit() const;
160  bool vsplit() const;
161  bool bsplit() const;
162 
163  protected:
164  CurvMap* cm;
165 
166  bool areaCalculated;
168  double area;
169 
171  bool diameterCalculated;
173  double diameter;
174 
176  int iro_cache;
177 
179  int next_vert(int i) const;
180  int prev_vert(int i) const;
181 
184  Element* get_neighbor(int ie) const;
185 
187  void ref_all_nodes();
189  void unref_all_nodes(HashTable* ht);
190  private:
191  unsigned nvert:30;
192 
193  friend class Mesh;
194  friend class MeshReader;
195  friend class MeshReaderH2D;
196  friend class MeshReaderH1DXML;
197  friend class MeshReaderH2DXML;
198  friend class PrecalcShapeset;
199  template<typename Scalar> friend class Space;
200  template<typename Scalar> friend class Adapt;
201  template<typename Scalar> friend class H1Space;
202  template<typename Scalar> friend class HcurlSpace;
203  template<typename Scalar> friend class HdivSpace;
204  template<typename Scalar> friend class L2Space;
205  template<typename Scalar> friend class KellyTypeAdapt;
206  template<typename Scalar> friend class DiscreteProblem;
207  template<typename Scalar> friend class Solution;
208  template<typename Scalar> friend class NeighborSearch;
209  template<typename Scalar> friend class Filter;
210  template<typename Scalar> friend class MeshFunction;
211  template<typename Scalar> friend class Global;
212  template<typename Scalar> friend class RefinementSelectors::Selector;
213  template<typename Scalar> friend class RefinementSelectors::POnlySelector;
214  template<typename Scalar> friend class RefinementSelectors::HOnlySelector;
215  template<typename Scalar> friend class RefinementSelectors::OptimumSelector;
216  template<typename Scalar> friend class RefinementSelectors::ProjBasedSelector;
217  template<typename Scalar> friend class RefinementSelectors::H1ProjBasedSelector;
218  template<typename Scalar> friend class RefinementSelectors::L2ProjBasedSelector;
219  template<typename Scalar> friend class RefinementSelectors::HcurlProjBasedSelector;
220  friend class Views::ScalarView;
221  friend class Views::Linearizer;
222  friend class RefMap;
223  friend class Traverse;
224  friend class Transformable;
225  friend class CurvMap;
226  friend class Views::Orderizer;
227  friend class Views::Vectorizer;
228  friend bool is_twin_nurbs(Element* e, int i);
229  friend int rtb_criterion(Element* e);
230  friend CurvMap* create_son_curv_map(Element* e, int son);
231  };
232 
246  class HERMES_API Mesh : public HashTable, public Hermes::Hermes2D::Mixins::StateQueryable
247  {
248  public:
249  Mesh();
250  virtual ~Mesh();
251 
254  void init(int size = H2D_DEFAULT_HASH_SIZE);
255 
257  virtual bool isOkay() const;
258  inline std::string getClassName() const { return "Mesh"; }
259 
261  bool rescale(double x_ref, double y_ref);
262 
264  void copy(const Mesh* mesh);
265 
267  void copy_base(Mesh* mesh);
268 
270  void copy_converted(Mesh* mesh);
271 
273  void create(int nv, double2* verts, int nt, int3* tris, std::string* tri_markers,
274  int nq, int4* quads, std::string* quad_markers, int nm, int2* mark, std::string* boundary_markers);
275 
277  void free();
278 
280  Element* get_element(int id) const;
281 
283  int get_num_elements() const;
284 
286  int get_num_base_elements() const;
287 
289  int get_num_used_base_elements() const;
290 
292  int get_num_active_elements() const;
293 
295  int get_max_element_id() const;
296 
298  int get_num_vertex_nodes() const;
299 
301  int get_num_edge_nodes() const;
302 
309  void refine_element_id(int id, int refinement = 0);
310 
313  void refine_all_elements(int refinement = 0, bool mark_as_initial = false);
314 
322  void refine_by_criterion(int (*criterion)(Element* e), int depth = 1, bool mark_as_initial = false);
323 
326  void refine_towards_vertex(int vertex_id, int depth = 1, bool mark_as_initial = false);
327 
332  void refine_towards_boundary(Hermes::vector<std::string> markers, int depth = 1, bool aniso = true, bool mark_as_initial = false);
333  void refine_towards_boundary(std::string marker, int depth = 1, bool aniso = true, bool mark_as_initial = false);
334 
336  void refine_in_area(std::string marker, int depth = 1, bool mark_as_initial = false);
338  void refine_in_areas(Hermes::vector<std::string> markers, int depth = 1, bool mark_as_initial = false);
339 
348  int* regularize(int n);
349 
352  void unrefine_element_id(int id);
353 
359  void unrefine_all_elements(bool keep_initial_refinements = true);
360 
362  Element* get_element_fast(int id) const;
363 
365  unsigned get_seq() const;
366 
368  void set_seq(unsigned seq);
369 
371  class HERMES_API ReferenceMeshCreator
372  {
373  public:
380  ReferenceMeshCreator(Mesh* coarse_mesh, int refinement = 0);
381 
384  virtual Mesh* create_ref_mesh();
385 
386  private:
388  Mesh* coarse_mesh;
389  int refinement;
390  };
391 
392  private:
394  void initial_single_check();
395  static void initial_multimesh_check(Hermes::vector<Mesh*> meshes);
396 
398  int get_edge_sons(Element* e, int edge, int& son1, int& son2) const;
399 
403  void convert_quads_to_triangles();
404 
406  void convert_to_base();
407 
408  void refine_element_to_quads_id(int id);
409 
410  void refine_triangle_to_quads(Mesh* mesh, Element* e, Element** elems_out = NULL);
411 
412  void refine_element_to_triangles_id(int id);
413 
414  void refine_quad_to_triangles(Element* e);
415 
420  void refine_quad_to_quads(Element* e, int refinement = 0);
421 
422  void convert_element_to_base_id(int id);
423  void convert_triangles_to_base(Element* e);
424  void convert_quads_to_base(Element* e);
425 
426  Array<Element> elements;
427  int nactive;
428  unsigned seq;
429 
430  int nbase, ntopvert;
431  int ninitial;
432 
433  void unrefine_element_internal(Element* e);
434 
439  Nurbs* reverse_nurbs(Nurbs* nurbs);
440  Node* get_base_edge_node(Element* base, int edge);
441 
442  int* parents;
443  int parents_size;
444 
445  int get_edge_degree(Node* v1, Node* v2);
446  void assign_parent(Element* e, int i);
447  void regularize_triangle(Element* e);
448  void regularize_quad(Element* e);
449  void flatten();
450 
451  class HERMES_API MarkersConversion
452  {
453  public:
454  MarkersConversion();
455 
459  int min_marker_unused;
460 
464  void insert_marker(int internal_marker, std::string user_marker);
465 
467  struct StringValid
468  {
469  StringValid();
470  StringValid(std::string marker, bool valid);
471  std::string marker;
472  bool valid;
473  };
474 
476  struct IntValid
477  {
478  IntValid();
479  IntValid(int marker, bool valid);
480  int marker;
481  bool valid;
482  };
483 
486  StringValid get_user_marker(int internal_marker) const;
487 
489  IntValid get_internal_marker(std::string user_marker) const;
490 
491  enum MarkersConversionType {
492  HERMES_ELEMENT_MARKERS_CONVERSION = 0,
493  HERMES_BOUNDARY_MARKERS_CONVERSION = 1
494  };
495 
496  virtual MarkersConversionType get_type() const = 0;
497 
498  protected:
501  std::map<int, std::string> conversion_table;
502 
505  std::map<std::string, int> conversion_table_inverse;
506 
507  friend class Space<double>;
508  friend class Space<std::complex<double> >;
509  friend class Mesh;
510  };
511 
514  class HERMES_API CurvedException : public Hermes::Exceptions::Exception
515  {
516  public:
519  CurvedException(int elementId);
520  CurvedException(const CurvedException & e);
521 
522  int getElementId() const;
523  private:
524  int elementId;
525  };
526 
527  class ElementMarkersConversion : public MarkersConversion
528  {
529  public:
530  ElementMarkersConversion();
531  virtual MarkersConversionType get_type() const;
532  };
533 
534  class BoundaryMarkersConversion : public MarkersConversion
535  {
536  public:
537  BoundaryMarkersConversion();
538  virtual MarkersConversionType get_type() const;
539  };
540 
541  ElementMarkersConversion element_markers_conversion;
542  BoundaryMarkersConversion boundary_markers_conversion;
543 
544  friend class MeshReaderH2D;
545  friend class MeshReaderH2DXML;
546  friend class MeshReaderH1DXML;
547  friend class MeshReaderExodusII;
548  friend class DiscreteProblem<double>;
549  friend class DiscreteProblem<std::complex<double> >;
550  friend class WeakForm<double>;
551  friend class WeakForm<std::complex<double> >;
552  template<typename Scalar> friend class Adapt;
553  friend class KellyTypeAdapt<double>;
554  template<typename Scalar> friend class Global;
555  friend class KellyTypeAdapt<std::complex<double> >;
556  template<typename Scalar> friend class Solution;
557  template<typename Scalar> friend class Filter;
558  template<typename Scalar> friend class MeshFunction;
559  friend class RefMap;
560  friend class Traverse;
561  friend class Transformable;
562  friend class Curved;
563  friend class Views::MeshView;
564  template<typename Scalar> friend class RefinementSelectors::Selector;
565  template<typename Scalar> friend class RefinementSelectors::POnlySelector;
566  template<typename Scalar> friend class RefinementSelectors::HOnlySelector;
567  template<typename Scalar> friend class RefinementSelectors::OptimumSelector;
568  template<typename Scalar> friend class RefinementSelectors::ProjBasedSelector;
569  template<typename Scalar> friend class RefinementSelectors::H1ProjBasedSelector;
570  template<typename Scalar> friend class RefinementSelectors::L2ProjBasedSelector;
571  template<typename Scalar> friend class RefinementSelectors::HcurlProjBasedSelector;
572  friend class PrecalcShapeset;
573  template<typename Scalar> friend class Space;
574  template<typename Scalar> friend class H1Space;
575  template<typename Scalar> friend class HcurlSpace;
576  template<typename Scalar> friend class HdivSpace;
577  template<typename Scalar> friend class L2Space;
578  friend class Views::ScalarView;
579  friend class Views::Orderizer;
580  public:
581  ElementMarkersConversion &get_element_markers_conversion();
582  BoundaryMarkersConversion &get_boundary_markers_conversion();
583 
584  /* node and son numbering on a triangle:
585 
586  -Triangle to triangles refinement
587 
588  vn[2] vn[2]
589 
590  * *
591 
592  / \ / \
593  / \ / \
594  / \ / \
595  / \ / son[2]\
596  / \ /_________\
597  en[2] / \ en[1] vn[0] * * vn[1]
598  * * vn[1] *-----------* vn[0]
599  / \ vn[2] * \ / * vn[2]
600  / \ / \ \ son[3]/ / \
601  / \ / \ \ / / \
602  / \ / \ \ / / \
603  / \ / son[0]\ \ / /son[1] \
604  / \ / \ * / \
605  *-------------*-------------* *-----------* *-----------*
606  vn[0] vn[1] vn[2] vn[0] vn[1]
607  vn[0] en[0] vn[1]
608 
609  -Triangle to quads refinement
610 
611  vn[2] vn[2]
612 
613  * *
614  / \ / \
615  / \ / \
616  / \ / \
617  / \ vn[3] * son[2]* vn[1]
618  / \ vn[3] * \ / * vn[2]
619  en[2] * * en[1] / \ \ / / \
620  / \ / \ vn[0] / \
621  / \ / \ * / \
622  / \ / \ / \
623  / * \ / vn[2] * * vn[3] \
624  / \ / | | \
625  / \ / son[0] | | son[1] \
626  / \ / | | \
627  *-------------*-------------* *-------------* *-------------*
628  vn[0] vn[1] vn[0] vn[1]
629  vn[0] en[0] vn[1]
630 
631  node and son numbering on a quad: refinement '0':
632 
633  vn[3] en[2] vn[2] vn[3] vn[2] vn[3] vn[2]
634 
635  *-------------*-------------* *------------* *------------*
636  | | | | | |
637  | | | | | |
638  | | | son[3] | | son[2] |
639  | | | | | |
640  | | | vn[1]| |vn[0] |
641  | | vn[0] *------------* *------------* vn[1]
642  en[3] * * en[1] vn[3] *------------* *------------* vn[2]
643  | | | vn[2]| |vn[3] |
644  | | | | | |
645  | | | son[0] | | son[1] |
646  | | | | | |
647  | | | | | |
648  | | *------------* *------------*
649  *-------------*-------------*
650  vn[0] vn[1] vn[0] vn[1]
651  vn[0] en[0] vn[1]
652 
653  refinement '1': refinement '2':
654 
655  vn[3] vn[2] vn[3] vn[2] vn[3] vn[2]
656 
657  *---------------------------* *------------* *------------*
658  | | | | | |
659  | | | | | |
660  | son[1] | | | | |
661  | | | | | |
662  | | | | | |
663  vn[0] *---------------------------* vn[1] | | | |
664  vn[3] *---------------------------* vn[2] | son[2] | | son[3] |
665  | | | | | |
666  | | | | | |
667  | son[0] | | | | |
668  | | | | | |
669  | | | | | |
670  *---------------------------* *------------* *------------*
671 
672  vn[0] vn[1] vn[0] vn[1] vn[0] vn[1]
673  */
674 
675  Element* create_quad(int marker, Node* v0, Node* v1, Node* v2, Node* v3, CurvMap* cm, int id = -1);
676  Element* create_triangle(int marker, Node* v0, Node* v1, Node* v2, CurvMap* cm, int id = -1);
677  void refine_element(Element* e, int refinement);
678 
681  Hermes::vector<std::pair<unsigned int, int> > refinements;
682 
688  void refine_quad(Element* e, int refinement, Element** sons_out = NULL);
689  void refine_triangle_to_triangles(Element* e, Element** sons = NULL);
690 
692  static double vector_length(double a_1, double a_2);
693 
695  static bool same_line(double p_1, double p_2, double q_1, double q_2, double r_1, double r_2);
696 
698  static bool is_convex(double a_1, double a_2, double b_1, double b_2);
699  static void check_triangle(int i, Node *&v0, Node *&v1, Node *&v2);
700  static void check_quad(int i, Node *&v0, Node *&v1, Node *&v2, Node *&v3);
701  };
702 
703  static Node* get_edge_node();
704  static Node* get_vertex_node(Node* v1, Node* v2);
705 
707  #define for_all_elements(e, mesh) \
708  for (int _id = 0, _max = (mesh)->get_max_element_id(); _id < _max; _id++) \
709  if(((e) = (mesh)->get_element_fast(_id))->used)
710 
711  #define for_all_base_elements(e, mesh) \
712  for (int _id = 0; _id < (mesh)->get_num_base_elements(); _id++) \
713  if(((e) = (mesh)->get_element_fast(_id))->used)
714 
715  #define for_all_base_elements_incl_inactive(e, mesh) \
716  for (int _id = 0; _id < (mesh)->get_num_base_elements(); _id++) \
717  if(((e) = (mesh)->get_element_fast(_id))->used || !((e) = (mesh)->get_element_fast(_id))->used)
718 
719  #define for_all_active_elements(e, mesh) \
720  for (int _id = 0, _max = (mesh)->get_max_element_id(); _id < _max; _id++) \
721  if(((e) = (mesh)->get_element_fast(_id))->used) \
722  if((e)->active)
723 
724  #define for_all_inactive_elements(e, mesh) \
725  for (int _id = 0, _max = (mesh)->get_max_element_id(); _id < _max; _id++) \
726  if(((e) = (mesh)->get_element_fast(_id))->used) \
727  if(!(e)->active)
728 
729  #define for_all_nodes(n, mesh) \
730  for (int _id = 0, _max = (mesh)->get_max_node_id(); _id < _max; _id++) \
731  if(((n) = (mesh)->get_node(_id))->used)
732 
733  #define for_all_vertex_nodes(n, mesh) \
734  for (int _id = 0, _max = (mesh)->get_max_node_id(); _id < _max; _id++) \
735  if(((n) = (mesh)->get_node(_id))->used) \
736  if(!(n)->type)
737 
738  #define for_all_edge_nodes(n, mesh) \
739  for (int _id = 0, _max = (mesh)->get_max_node_id(); _id < _max; _id++) \
740  if(((n) = (mesh)->get_node(_id))->used) \
741  if((n)->type)
742 
743  const int TOP_LEVEL_REF = 123456;
744  }
745 }
746 #endif