Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_seq_check.cc
Go to the documentation of this file.
1
21#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
29namespace rheolef {
30
31// TODO: quadrangle may also be convex
32// ----------------------------------------------------------------------------
33// compute measure
34// ----------------------------------------------------------------------------
35template <class T>
36static
37T
38meas_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}
45template <class T>
46static
47T
48meas_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}
56template <class T>
57static
58void
59meas_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}
74template <class T>
75static
76bool
77check_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;
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 }
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 }
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 }
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 }
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 }
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// ----------------------------------------------------------------------------
148template <class T>
149bool
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
161template <class T>
162bool
164{
165 // TODO
166 return true;
167}
168#endif // _RHEOLEF_HAVE_MPI
169// ----------------------------------------------------------------------------
170// instanciation in library
171// ----------------------------------------------------------------------------
172template class geo_rep<Float,sequential>;
173#ifdef _RHEOLEF_HAVE_MPI
174template class geo_rep<Float,distributed>;
175#endif // _RHEOLEF_HAVE_MPI
176
177} // namespace rheolef
see the geo_element page for the full documentation
base::const_iterator const_iterator
Definition geo.h:797
sequential mesh representation
Definition geo.h:778
static const variant_type H
static const variant_type q
static const variant_type e
static const variant_type p
static const variant_type T
static const variant_type P
static const variant_type t
#define error_macro(message)
Definition dis_macros.h:49
#define warning_macro(message)
Definition dis_macros.h:53
Expr1::float_type T
Definition field_expr.h:230
This file is part of Rheolef.
T mixt(const point_basic< T > &u, const point_basic< T > &v, const point_basic< T > &w)
Definition point.h:564
T vect2d(const point_basic< T > &v, const point_basic< T > &w)
Definition point.h:557
Definition sphere.icc:25