Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
verdict_defines.hpp
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Module: $RCSfile: verdict_defines.hpp,v $
4 
5  Copyright (c) 2006 Sandia Corporation.
6  All rights reserved.
7  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
8 
9  This software is distributed WITHOUT ANY WARRANTY; without even
10  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11  PURPOSE. See the above copyright notice for more information.
12 
13 =========================================================================*/
14 
15 /*
16  * verdict_defines.cpp contains common definitions
17  *
18  * This file is part of VERDICT
19  *
20  */
21 
22 #ifndef VERDICT_DEFINES
23 #define VERDICT_DEFINES
24 
25 #include <cmath>
26 #include "v_vector.h"
27 #include "VerdictVector.hpp"
28 
29 enum VerdictBoolean
30 {
31  VERDICT_FALSE = 0,
32  VERDICT_TRUE = 1
33 };
34 
35 #define VERDICT_MIN( a, b ) ( ( a ) < ( b ) ? ( a ) : ( b ) )
36 #define VERDICT_MAX( a, b ) ( ( a ) > ( b ) ? ( a ) : ( b ) )
37 
38 inline double determinant( double a, double b, double c, double d )
39 {
40  return ( ( a ) * ( d ) - ( b ) * ( c ) );
41 }
42 
43 inline double determinant( const VerdictVector& v1, const VerdictVector& v2, const VerdictVector& v3 )
44 {
45  return v1 % ( v2 * v3 );
46 }
47 
48 #define jacobian_matrix( a, b, c, d, e, f, g ) \
49  double jac_mat_tmp; \
50  jac_mat_tmp = sqrt( a ); \
51  if( jac_mat_tmp == 0 ) \
52  { \
53  ( d ) = 0; \
54  ( e ) = 0; \
55  ( f ) = 0; \
56  ( g ) = 0; \
57  } \
58  else \
59  { \
60  ( d ) = jac_mat_tmp; \
61  ( e ) = 0; \
62  ( f ) = ( b ) / jac_mat_tmp; \
63  ( g ) = ( c ) / jac_mat_tmp; \
64  }
65 
66 // this assumes that detmw != 0
67 #define form_t( m11, m21, m12, m22, mw11, mw21, mw12, mw22, detmw, xm11, xm21, xm12, xm22 ) \
68  xm11 = ( ( m11 ) * ( mw22 ) - ( m12 ) * ( mw21 ) ) / ( detmw ); \
69  ( xm21 ) = ( ( m21 ) * ( mw22 ) - ( m22 ) * ( mw21 ) ) / ( detmw ); \
70  ( xm12 ) = ( ( m12 ) * ( mw11 ) - ( m11 ) * ( mw12 ) ) / ( detmw ); \
71  ( xm22 ) = ( ( m22 ) * ( mw11 ) - ( m21 ) * ( mw12 ) ) / ( detmw );
72 
73 extern double verdictSqrt2;
74 
75 inline double normalize_jacobian( double jacobi,
76  VerdictVector& v1,
77  VerdictVector& v2,
78  VerdictVector& v3,
79  int tet_flag = 0 )
80 {
81  double return_value = 0.0;
82 
83  if( jacobi != 0.0 )
84  {
85 
86  double l1, l2, l3, length_product;
87  // Note: there may be numerical problems if one is a lot shorter
88  // than the others this way. But scaling each vector before the
89  // triple product would involve 3 square roots instead of just
90  // one.
91  l1 = v1.length_squared();
92  l2 = v2.length_squared();
93  l3 = v3.length_squared();
94  length_product = sqrt( l1 * l2 * l3 );
95 
96  // if some numerical scaling problem, or just plain roundoff,
97  // then push back into range [-1,1].
98  if( length_product < fabs( jacobi ) )
99  {
100  length_product = fabs( jacobi );
101  }
102 
103  if( tet_flag == 1 )
104  return_value = verdictSqrt2 * jacobi / length_product;
105  else
106  return_value = jacobi / length_product;
107  }
108  return return_value;
109 }
110 
111 inline double norm_squared( double m11, double m21, double m12, double m22 )
112 {
113  return m11 * m11 + m21 * m21 + m12 * m12 + m22 * m22;
114 }
115 
116 #define metric_matrix( m11, m21, m12, m22, gm11, gm12, gm22 ) \
117  gm11 = ( m11 ) * ( m11 ) + ( m21 ) * ( m21 ); \
118  ( gm12 ) = ( m11 ) * ( m12 ) + ( m21 ) * ( m22 ); \
119  ( gm22 ) = ( m12 ) * ( m12 ) + ( m22 ) * ( m22 );
120 
121 inline int skew_matrix( double gm11,
122  double gm12,
123  double gm22,
124  double det,
125  double& qm11,
126  double& qm21,
127  double& qm12,
128  double& qm22 )
129 {
130  double tmp = sqrt( gm11 * gm22 );
131  if( tmp == 0 )
132  {
133  return false;
134  }
135 
136  qm11 = 1;
137  qm21 = 0;
138  qm12 = gm12 / tmp;
139  qm22 = det / tmp;
140  return true;
141 }
142 
143 inline void inverse( const VerdictVector& x1,
144  const VerdictVector& x2,
145  const VerdictVector& x3,
146  VerdictVector& u1,
147  VerdictVector& u2,
148  VerdictVector& u3 )
149 {
150  double detx = determinant( x1, x2, x3 );
151  VerdictVector rx1, rx2, rx3;
152 
153  rx1.set( x1.x(), x2.x(), x3.x() );
154  rx2.set( x1.y(), x2.y(), x3.y() );
155  rx3.set( x1.z(), x2.z(), x3.z() );
156 
157  u1 = rx2 * rx3;
158  u2 = rx3 * rx1;
159  u3 = rx1 * rx2;
160 
161  u1 /= detx;
162  u2 /= detx;
163  u3 /= detx;
164 }
165 
166 /*
167 inline void form_T(double a1[3],
168  double a2[3],
169  double a3[3],
170  double w1[3],
171  double w2[3],
172  double w3[3],
173  double t1[3],
174  double t2[3],
175  double t3[3] )
176 {
177 
178  double x1[3], x2[3], x3[3];
179  double ra1[3], ra2[3], ra3[3];
180 
181  x1[0] = a1[0];
182  x1[1] = a2[0];
183  x1[2] = a3[0];
184 
185  x2[0] = a1[1];
186  x2[1] = a2[1];
187  x2[2] = a3[1];
188 
189  x3[0] = a1[2];
190  x3[1] = a2[2];
191  x3[2] = a3[2];
192 
193  inverse(w1,w2,w3,x1,x2,x3);
194 
195  t1[0] = dot_product(ra1, x1);
196  t1[1] = dot_product(ra1, x2);
197  t1[2] = dot_product(ra1, x3);
198 
199  t2[0] = dot_product(ra2, x1);
200  t2[0] = dot_product(ra2, x2);
201  t2[0] = dot_product(ra2, x3);
202 
203  t3[0] = dot_product(ra3, x1);
204  t3[0] = dot_product(ra3, x2);
205  t3[0] = dot_product(ra3, x3);
206 
207 }
208 */
209 
210 inline void form_Q( const VerdictVector& v1,
211  const VerdictVector& v2,
212  const VerdictVector& v3,
213  VerdictVector& q1,
214  VerdictVector& q2,
215  VerdictVector& q3 )
216 {
217 
218  double g11, g12, g13, g22, g23, g33;
219 
220  g11 = v1 % v1;
221  g12 = v1 % v2;
222  g13 = v1 % v3;
223  g22 = v2 % v2;
224  g23 = v2 % v3;
225  g33 = v3 % v3;
226 
227  double rtg11 = sqrt( g11 );
228  double rtg22 = sqrt( g22 );
229  double rtg33 = sqrt( g33 );
230  VerdictVector temp1;
231 
232  temp1 = v1 * v2;
233 
234  double cross = sqrt( temp1 % temp1 );
235 
236  double q11, q21, q31;
237  double q12, q22, q32;
238  double q13, q23, q33;
239 
240  q11 = 1;
241  q21 = 0;
242  q31 = 0;
243 
244  q12 = g12 / rtg11 / rtg22;
245  q22 = cross / rtg11 / rtg22;
246  q32 = 0;
247 
248  q13 = g13 / rtg11 / rtg33;
249  q23 = ( g11 * g23 - g12 * g13 ) / rtg11 / rtg33 / cross;
250  temp1 = v2 * v3;
251  q33 = ( v1 % temp1 ) / rtg33 / cross;
252 
253  q1.set( q11, q21, q31 );
254  q2.set( q12, q22, q32 );
255  q3.set( q13, q23, q33 );
256 }
257 
258 inline void product( VerdictVector& a1,
259  VerdictVector& a2,
260  VerdictVector& a3,
261  VerdictVector& b1,
262  VerdictVector& b2,
263  VerdictVector& b3,
264  VerdictVector& c1,
265  VerdictVector& c2,
266  VerdictVector& c3 )
267 {
268 
269  VerdictVector x1, x2, x3;
270 
271  x1.set( a1.x(), a2.x(), a3.x() );
272  x2.set( a1.y(), a2.y(), a3.y() );
273  x3.set( a1.z(), a2.z(), a3.z() );
274 
275  c1.set( x1 % b1, x2 % b1, x3 % b1 );
276  c2.set( x1 % b2, x2 % b2, x3 % b2 );
277  c3.set( x1 % b3, x2 % b3, x3 % b3 );
278 }
279 
280 inline double norm_squared( VerdictVector& x1, VerdictVector& x2, VerdictVector& x3 )
281 
282 {
283  return ( x1 % x1 ) + ( x2 % x2 ) + ( x3 % x3 );
284 }
285 
286 inline double skew_x( VerdictVector& q1,
287  VerdictVector& q2,
288  VerdictVector& q3,
289  VerdictVector& qw1,
290  VerdictVector& qw2,
291  VerdictVector& qw3 )
292 {
293  double normsq1, normsq2, kappa;
294  VerdictVector u1, u2, u3;
295  VerdictVector x1, x2, x3;
296 
297  inverse( qw1, qw2, qw3, u1, u2, u3 );
298  product( q1, q2, q3, u1, u2, u3, x1, x2, x3 );
299  inverse( x1, x2, x3, u1, u2, u3 );
300  normsq1 = norm_squared( x1, x2, x3 );
301  normsq2 = norm_squared( u1, u2, u3 );
302  kappa = sqrt( normsq1 * normsq2 );
303 
304  double skew = 0;
305  if( kappa > VERDICT_DBL_MIN ) skew = 3 / kappa;
306 
307  return skew;
308 }
309 
310 #endif