Rheolef  7.1
an efficient C++ finite element environment
geo_seq_check.cc
Go to the documentation of this file.
1 #include "rheolef/geo.h"
22 #include "rheolef/geo_domain.h"
23 #include "rheolef/rheostream.h"
24 #include "rheolef/iorheo.h"
25 #include "rheolef/space_numbering.h"
26 
27 #include "geo_header.h"
28 
29 namespace rheolef {
30 
31 // TODO: quadrangle may also be convex
32 // ----------------------------------------------------------------------------
33 // compute measure
34 // ----------------------------------------------------------------------------
35 template <class T>
36 static
37 T
38 meas_t (
39  const point_basic<T>& a,
40  const point_basic<T>& b,
41  const point_basic<T>& c)
42 {
43  return vect2d(b-a, c-a)/2.0;
44 }
45 template <class T>
46 static
47 T
48 meas_T (
49  const point_basic<T>& a,
50  const point_basic<T>& b,
51  const point_basic<T>& c,
52  const point_basic<T>& d)
53 {
54  return mixt(b-a, c-a, d-a)/6;
55 }
56 template <class T>
57 static
58 void
59 meas_P (
60  const point_basic<T>& a0,
61  const point_basic<T>& a1,
62  const point_basic<T>& a2,
63  const point_basic<T>& a3,
64  const point_basic<T>& a4,
65  const point_basic<T>& a5,
66  T& meas1,
67  T& meas2,
68  T& meas3)
69 {
70  meas1 = meas_T (a0, a1, a2, a3);
71  meas2 = meas_T (a1, a2, a3, a4);
72  meas3 = meas_T (a2, a3, a4, a5);
73 }
74 template <class T>
75 static
76 bool
77 check_element (const geo_element& K, const disarray<point_basic<T>,sequential>& p, bool verbose = false)
78 {
79  std::cerr << std::setprecision(std::numeric_limits<T>::digits10);
80  bool orient = true;
81  switch (K.variant()) {
83  return true;
84  case reference_element::e: {
85  T meas = p[K[1]][0] - p[K[0]][0];
86  orient = (meas > 0);
87  if (verbose && !orient) {
88  warning_macro("meas(K)="<<meas<<": invalid negative orientation for element K.dis_ie="<<K.dis_ie());
89  }
90  break;
91  }
92  case reference_element::t: {
93  T meas = meas_t (p[K[0]], p[K[1]], p[K[2]]);
94  orient = (meas > 0);
95  if (verbose && !orient) {
96  warning_macro("meas(K)="<<meas<<": invalid negative orientation for element K.dis_ie="<<K.dis_ie());
97  }
98  break;
99  }
100  case reference_element::q: {
101  T meas1 = meas_t (p[K[0]], p[K[1]], p[K[2]]); // t=012
102  T meas2 = meas_t (p[K[0]], p[K[2]], p[K[3]]); // t=023
103  orient = (meas1 > 0 && meas2 > 0);
104  if (verbose && !orient) {
105  warning_macro("meas(K)="<<meas1<<"+"<<meas2<<": invalid negative orientation for element K.dis_ie="<<K.dis_ie());
106  }
107  break;
108  }
109  case reference_element::T: {
110  T meas = meas_T (p[K[0]], p[K[1]], p[K[2]], p[K[3]]);
111  orient = (meas > 0);
112  if (verbose && !orient) {
113  warning_macro("meas(K)="<<meas<<": invalid negative orientation for element K.dis_ie="<<K.dis_ie());
114  }
115  break;
116  }
117  case reference_element::P: {
118  T meas1, meas2, meas3;
119  meas_P (p[K[0]], p[K[1]], p[K[2]], p[K[3]], p[K[4]], p[K[5]], meas1, meas2, meas3);
120  orient = (meas1 > 0 && meas2 > 0 && meas3 > 0);
121  if (verbose && !orient) {
122  warning_macro("meas(K)="<<meas1<<"+"<<meas2<<"+"<<meas3
123  <<": invalid negative orientation for element K.dis_ie="<<K.dis_ie());
124  }
125  break;
126  }
127  case reference_element::H: {
128  T meas1, meas2, meas3, meas4, meas5, meas6;
129  meas_P (p[K[0]], p[K[1]], p[K[2]], p[K[4]], p[K[5]], p[K[6]], meas1, meas2, meas3);
130  meas_P (p[K[0]], p[K[2]], p[K[3]], p[K[4]], p[K[6]], p[K[7]], meas4, meas5, meas6);
131  orient = (meas1 > 0 && meas2 > 0 && meas3 > 0 && meas4 > 0 && meas5 > 0 && meas6 > 0);
132  if (verbose && !orient) {
133  warning_macro("meas(K)="<<meas1<<"+"<<meas2<<"+"<<meas3
134  <<"+"<<meas4<<"+"<<meas5<<"+"<<meas6
135  <<": invalid negative orientation for element K.dis_ie="<<K.dis_ie());
136  }
137  break;
138  }
139  default:
140  error_macro ("unexpected element variant `"<<K.variant()<<"'");
141  }
142  bool is_convex = true; // TODO: q,P,H
143  return (orient && is_convex);
144 }
145 // ----------------------------------------------------------------------------
146 // check element orientation
147 // ----------------------------------------------------------------------------
148 template <class T>
149 bool
151 {
152  if (base::_dimension != base::_gs._map_dimension) return true;
153  bool status = true;
154  for (const_iterator iter = begin(base::_gs._map_dimension), last = end(base::_gs._map_dimension); iter != last; iter++) {
155  const geo_element& K = *iter;
156  status = status && check_element (K, base::_node, verbose);
157  }
158  return status;
159 }
160 #ifdef _RHEOLEF_HAVE_MPI
161 template <class T>
162 bool
164 {
165  // TODO
166  return true;
167 }
168 #endif // _RHEOLEF_HAVE_MPI
169 // ----------------------------------------------------------------------------
170 // instanciation in library
171 // ----------------------------------------------------------------------------
172 template class geo_rep<Float,sequential>;
173 #ifdef _RHEOLEF_HAVE_MPI
174 template class geo_rep<Float,distributed>;
175 #endif // _RHEOLEF_HAVE_MPI
176 
177 } // namespace rheolef
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
rheolef::reference_element::H
static const variant_type H
Definition: reference_element.h:81
warning_macro
#define warning_macro(message)
Definition: dis_macros.h:53
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
geo_header.h
rheolef::reference_element::T
static const variant_type T
Definition: reference_element.h:79
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::geo_rep
sequential mesh representation
Definition: geo.h:778
p
Definition: sphere.icc:25
a
Definition: diffusion_isotropic.h:25
rheolef::vect2d
T vect2d(const point_basic< T > &v, const point_basic< T > &w)
Definition: point.h:558
rheolef::mixt
T mixt(const point_basic< T > &u, const point_basic< T > &v, const point_basic< T > &w)
Definition: point.h:565
rheolef::geo_iterator
geo iterator
Definition: geo.h:193
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::point_basic< T >
point_basic< T >
Definition: piola_fem.h:135
mkgeo_ball.verbose
verbose
Definition: mkgeo_ball.sh:133
rheolef::reference_element::p
static const variant_type p
Definition: reference_element.h:75
rheolef::reference_element::q
static const variant_type q
Definition: reference_element.h:78
rheolef::reference_element::P
static const variant_type P
Definition: reference_element.h:80
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
check
void check(Float p, const field &uh)
Definition: p_laplacian_post.cc:49
mkgeo_contraction.status
status
Definition: mkgeo_contraction.sh:290
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153
T
Expr1::float_type T
Definition: field_expr.h:261