1 /*! \file evm.h 2 * \ingroup gm 3 */ 4 5 /****************************************************************************/ 6 /* */ 7 /* File: evm.h */ 8 /* */ 9 /* Purpose: elementary vector manipulations, header for evm.c */ 10 /* */ 11 /* Author: Klaus Johannsen */ 12 /* Institut fuer Computeranwendungen */ 13 /* Universitaet Stuttgart */ 14 /* Pfaffenwaldring 27 */ 15 /* 70569 Stuttgart */ 16 /* internet: ug@ica3.uni-stuttgart.de */ 17 /* */ 18 /* History: 8.12.94 begin, ug3-version */ 19 /* */ 20 /* Remarks: */ 21 /* */ 22 /****************************************************************************/ 23 24 25 /****************************************************************************/ 26 /* */ 27 /* auto include mechanism and other include files */ 28 /* */ 29 /****************************************************************************/ 30 31 #ifndef __EVM__ 32 #define __EVM__ 33 34 #include <cmath> 35 36 #include <dune/uggrid/low/debug.h> 37 #include <dune/uggrid/low/namespace.h> 38 #include <dune/uggrid/low/ugtypes.h> 39 40 #include "gm.h" 41 42 START_UGDIM_NAMESPACE 43 44 /****************************************************************************/ 45 /* */ 46 /* defines in the following order */ 47 /* */ 48 /* compile time constants defining static data size (i.e. arrays) */ 49 /* other constants */ 50 /* macros */ 51 /* */ 52 /****************************************************************************/ 53 54 /****************************************************************************/ 55 /* */ 56 /* macros */ 57 /* */ 58 /****************************************************************************/ 59 60 61 /* space dimension indices */ 62 #define _X_ 0 63 #define _Y_ 1 64 #define _Z_ 2 65 66 /* misc macros */ 67 #define SQRT(a) sqrt((double)(a)) 68 #define POW(a,b) pow((double)(a),(double)(b)) 69 #define ISNaN(x) (!((x)-(x)==0) || ((x)-(x)!=0)) 70 71 /* macros for coord points */ 72 #define COPY_SC_TO_SH(p1,p2) (p2).x=(short)((p1).x);(p2).y=(short)((p1).y) 73 #define CP_SUBTRACT(A,B,C) {(C).x = (A).x - (B).x;\ 74 (C).y = (A).y - (B).y;} 75 #define CP_LIMCOMB(a,A,b,B,C) {(C).x = (DOUBLE)(a)*(A).x + (DOUBLE)(b)*(B).x;\ 76 (C).y = (DOUBLE)(a)*(A).y + (DOUBLE)(b)*(B).y;} 77 #define CP_SCALARPRODUCT(A,B,c) (c) = (A).x*(B).x + (A).y*(B).y; 78 #define CP_EUKLIDNORM(A,b) (b) = (DOUBLE)sqrt((double)((A).x*(A).x+(A).y*(A).y)); 79 80 /* macros for 1D vector operations */ 81 #define V1_COPY(A,C) {(C)[0] = (A)[0];} 82 83 /* macros for 2D vector operations */ 84 #define V2_LINCOMB(a,A,b,B,C) {(C)[0] = (a)*(A)[0] + (b)*(B)[0];\ 85 (C)[1] = (a)*(A)[1] + (b)*(B)[1];} 86 #define V2_COPY(A,C) {(C)[0] = (A)[0];\ 87 (C)[1] = (A)[1];} 88 #define V2_SUBTRACT(A,B,C) {(C)[0] = (A)[0] - (B)[0];\ 89 (C)[1] = (A)[1] - (B)[1];} 90 #define V2_ADD(A,B,C) {(C)[0] = (A)[0] + (B)[0];\ 91 (C)[1] = (A)[1] + (B)[1];} 92 #define V2_SCALE(c,C) {(C)[0] = (c)*(C)[0];\ 93 (C)[1] = (c)*(C)[1];} 94 #define V2_VECTOR_PRODUCT(A,B,c) (c) = (A)[0]*(B)[1] - (A)[1]*(B)[0]; 95 #define V2_ISEQUAL(A,B) ((std::abs((A)[0]-(B)[0])<SMALL_C)&&(std::abs((A)[1]-(B)[1])<SMALL_C)) 96 #define V2_EUKLIDNORM(A,b) (b) = sqrt((double)((A)[0]*(A)[0]+(A)[1]*(A)[1])); 97 #define V2_EUKLIDNORM_OF_DIFF(A,B,b) (b) = sqrt((double)(((A)[0]-(B)[0])*((A)[0]-(B)[0])+((A)[1]-(B)[1])*((A)[1]-(B)[1]))); 98 #define V2_CLEAR(A) {(A)[0] = 0.0; (A)[1]= 0.0;} 99 #define V2_SCALAR_PRODUCT(A,B,c) (c) = (A)[0]*(B)[0]+(A)[1]*(B)[1]; 100 #define V2_SCAL_PROD(A,B) ((A)[0]*(B)[0]+(A)[1]*(B)[1]) 101 102 103 /* macros for 2D matrix-vector operations */ 104 #define M2_TIMES_V2(M,A,B) {(B)[0] = (M)[0]*(A)[0] + (M)[2]*(A)[1];\ 105 (B)[1] = (M)[1]*(A)[0] + (M)[3]*(A)[1];} 106 #define MM2_TIMES_V2(M,A,B) {(B)[0] = (M)[0][0]*(A)[0] + (M)[0][1]*(A)[1];\ 107 (B)[1] = (M)[1][0]*(A)[0] + (M)[1][1]*(A)[1];} 108 109 #define MT2_TIMES_V2(M,A,B) {(B)[0] = (M)[0][0]*(A)[0] + (M)[1][0]*(A)[1];\ 110 (B)[1] = (M)[0][1]*(A)[0] + (M)[1][1]*(A)[1];} 111 112 #define MD2_TIMES_V2(M,A,B) {(B)[0] = (M)[0]*(A)[0];\ 113 (B)[1] = (M)[1]*(A)[1];} 114 115 /* macros for matrix operations */ 116 #define M2_DET(M) ((M)[0]*(M)[3] - (M)[1]*(M)[2]) 117 118 #define M2_INVERT(M,IM,det) \ 119 { DOUBLE invdet; \ 120 det = (M)[0][0]*(M)[1][1]-(M)[1][0]*(M)[0][1]; \ 121 if (std::abs((det))<SMALL_D*SMALL_D) det= 0.; \ 122 else { \ 123 invdet = 1.0 / (det); \ 124 (IM)[0][0] = (M)[1][1]*invdet; \ 125 (IM)[1][0] = -(M)[1][0]*invdet; \ 126 (IM)[0][1] = -(M)[0][1]*invdet; \ 127 (IM)[1][1] = (M)[0][0]*invdet;}} 128 129 /* macros for vector operations */ 130 #define V3_LINCOMB(a,A,b,B,C) {(C)[0] = (a)*(A)[0] + (b)*(B)[0];\ 131 (C)[1] = (a)*(A)[1] + (b)*(B)[1];\ 132 (C)[2] = (a)*(A)[2] + (b)*(B)[2];} 133 #define V3_COPY(A,C) {(C)[0] = (A)[0];\ 134 (C)[1] = (A)[1];\ 135 (C)[2] = (A)[2];} 136 #define V3_SUBTRACT(A,B,C) {(C)[0] = (A)[0] - (B)[0];\ 137 (C)[1] = (A)[1] - (B)[1];\ 138 (C)[2] = (A)[2] - (B)[2];} 139 #define V3_ADD(A,B,C) {(C)[0] = (A)[0] + (B)[0];\ 140 (C)[1] = (A)[1] + (B)[1];\ 141 (C)[2] = (A)[2] + (B)[2];} 142 #define V3_SCALE(c,C) {(C)[0] = (c)*(C)[0];\ 143 (C)[1] = (c)*(C)[1];\ 144 (C)[2] = (c)*(C)[2];} 145 #define V3_VECTOR_PRODUCT(A,B,C) {(C)[0] = (A)[1]*(B)[2] - (A)[2]*(B)[1];\ 146 (C)[1] = (A)[2]*(B)[0] - (A)[0]*(B)[2];\ 147 (C)[2] = (A)[0]*(B)[1] - (A)[1]*(B)[0];} 148 #define V3_EUKLIDNORM(A,b) (b) = (sqrt((double)((A)[0]*(A)[0]+(A)[1]*(A)[1]+(A)[2]*(A)[2]))); 149 #define V3_ISEQUAL(A,B) ((std::abs((A)[0]-(B)[0])<SMALL_C)&&(std::abs((A)[1]-(B)[1])<SMALL_C)&&(std::abs((A)[2]-(B)[2])<SMALL_C)) 150 #define V3_EUKLIDNORM_OF_DIFF(A,B,b) (b) = (sqrt((double)(((A)[0]-(B)[0])*((A)[0]-(B)[0])+((A)[1]-(B)[1])*((A)[1]-(B)[1])+((A)[2]-(B)[2])*((A)[2]-(B)[2])))); 151 #define V3_CLEAR(A) {(A)[0] = 0.0; (A)[1]= 0.0; (A)[2] = 0.0;} 152 #define V3_SCALAR_PRODUCT(A,B,c) (c) = ((A)[0]*(B)[0]+(A)[1]*(B)[1]+(A)[2]*(B)[2]); 153 #define V3_SCAL_PROD(A,B) ((A)[0]*(B)[0]+(A)[1]*(B)[1]+(A)[2]*(B)[2]) 154 155 /* macros for matrix-vector operations */ 156 #define M3_TIMES_V3(M,A,B) {(B)[0] = (M)[0]*(A)[0] + (M)[3]*(A)[1] + (M)[6]*(A)[2];\ 157 (B)[1] = (M)[1]*(A)[0] + (M)[4]*(A)[1] + (M)[7]*(A)[2];\ 158 (B)[2] = (M)[2]*(A)[0] + (M)[5]*(A)[1] + (M)[8]*(A)[2];} 159 160 #define MT3_TIMES_V3(M,A,B) {(B)[0] = (M)[0][0]*(A)[0] + (M)[1][0]*(A)[1] + (M)[2][0]*(A)[2];\ 161 (B)[1] = (M)[0][1]*(A)[0] + (M)[1][1]*(A)[1] + (M)[2][1]*(A)[2];\ 162 (B)[2] = (M)[0][2]*(A)[0] + (M)[1][2]*(A)[1] + (M)[2][2]*(A)[2];} 163 164 /* macros for matrix operations */ 165 #define M3_DET(M) ((M)[0]*(M)[4]*(M)[8] + (M)[1]*(M)[5]*(M)[6] + (M)[2]*(M)[3]*(M)[7] \ 166 -(M)[2]*(M)[4]*(M)[6] - (M)[0]*(M)[5]*(M)[7] - (M)[1]*(M)[3]*(M)[8]) 167 168 #define M3_INVERT(M,IM,det) \ 169 { DOUBLE invdet; \ 170 (det) = (M)[0][0]*(M)[1][1]*(M)[2][2] \ 171 + (M)[0][1]*(M)[1][2]*(M)[2][0] \ 172 + (M)[0][2]*(M)[1][0]*(M)[2][1] \ 173 - (M)[0][2]*(M)[1][1]*(M)[2][0] \ 174 - (M)[0][0]*(M)[1][2]*(M)[2][1] \ 175 - (M)[0][1]*(M)[1][0]*(M)[2][2]; \ 176 if (std::abs((det))<SMALL_D*SMALL_D) \ 177 return (1); \ 178 invdet = 1.0 / (det); \ 179 (IM)[0][0] = ( (M)[1][1]*(M)[2][2] - (M)[1][2]*(M)[2][1]) * invdet; \ 180 (IM)[0][1] = (-(M)[0][1]*(M)[2][2] + (M)[0][2]*(M)[2][1]) * invdet; \ 181 (IM)[0][2] = ( (M)[0][1]*(M)[1][2] - (M)[0][2]*(M)[1][1]) * invdet; \ 182 (IM)[1][0] = (-(M)[1][0]*(M)[2][2] + (M)[1][2]*(M)[2][0]) * invdet; \ 183 (IM)[1][1] = ( (M)[0][0]*(M)[2][2] - (M)[0][2]*(M)[2][0]) * invdet; \ 184 (IM)[1][2] = (-(M)[0][0]*(M)[1][2] + (M)[0][2]*(M)[1][0]) * invdet; \ 185 (IM)[2][0] = ( (M)[1][0]*(M)[2][1] - (M)[1][1]*(M)[2][0]) * invdet; \ 186 (IM)[2][1] = (-(M)[0][0]*(M)[2][1] + (M)[0][1]*(M)[2][0]) * invdet; \ 187 (IM)[2][2] = ( (M)[0][0]*(M)[1][1] - (M)[0][1]*(M)[1][0]) * invdet;} 188 189 /* macros for exact solver (EX) */ 190 #define EX_MAT(m,b,i,j) ((m)[2*(b)*(i) + (j)]) 191 192 /****************************************************************************/ 193 /* */ 194 /* typedef of DIM-routines */ 195 /* */ 196 /****************************************************************************/ 197 198 #ifdef __TWODIM__ 199 200 #define V_BDIM_COPY(A,C) V1_COPY(A,C) 201 202 #define V_DIM_LINCOMB(a,A,b,B,C) V2_LINCOMB(a,A,b,B,C) 203 #define V_DIM_COPY(A,C) V2_COPY(A,C) 204 #define V_DIM_SUBTRACT(A,B,C) V2_SUBTRACT(A,B,C) 205 #define V_DIM_ADD(A,B,C) V2_ADD(A,B,C) 206 #define V_DIM_SCALE(c,C) V2_SCALE(c,C) 207 #define V_DIM_ISEQUAL(A,B) V2_ISEQUAL(A,B) 208 #define V_DIM_EUKLIDNORM(A,b) V2_EUKLIDNORM(A,b) 209 #define V_DIM_EUKLIDNORM_OF_DIFF(A,B,b) V2_EUKLIDNORM_OF_DIFF(A,B,b) 210 #define V_DIM_CLEAR(A) V2_CLEAR(A) 211 #define V_DIM_SCALAR_PRODUCT(A,B,c) V2_SCALAR_PRODUCT(A,B,c) 212 #define V_DIM_SCAL_PROD(A,B) V2_SCAL_PROD(A,B) 213 214 #define MT_TIMES_V_DIM(M,A,B) MT2_TIMES_V2(M,A,B) 215 #define M_DIM_INVERT(M,IM,det) M2_INVERT(M,IM,det) 216 217 #endif 218 219 #ifdef __THREEDIM__ 220 221 #define V_BDIM_COPY(A,C) V2_COPY(A,C) 222 223 #define V_DIM_LINCOMB(a,A,b,B,C) V3_LINCOMB(a,A,b,B,C) 224 #define V_DIM_COPY(A,C) V3_COPY(A,C) 225 #define V_DIM_SUBTRACT(A,B,C) V3_SUBTRACT(A,B,C) 226 #define V_DIM_ADD(A,B,C) V3_ADD(A,B,C) 227 #define V_DIM_SCALE(c,C) V3_SCALE(c,C) 228 #define V_DIM_ISEQUAL(A,B) V3_ISEQUAL(A,B) 229 #define V_DIM_EUKLIDNORM(A,b) V3_EUKLIDNORM(A,b) 230 #define V_DIM_EUKLIDNORM_OF_DIFF(A,B,b) V3_EUKLIDNORM_OF_DIFF(A,B,b) 231 #define V_DIM_CLEAR(A) V3_CLEAR(A) 232 #define V_DIM_SCALAR_PRODUCT(A,B,c) V3_SCALAR_PRODUCT(A,B,c) 233 #define V_DIM_SCAL_PROD(A,B) V3_SCAL_PROD(A,B) 234 235 #define MT_TIMES_V_DIM(M,A,B) MT3_TIMES_V3(M,A,B) 236 #define M_DIM_INVERT(M,IM,det) M3_INVERT(M,IM,det) 237 238 #endif 239 240 /****************************************************************************/ 241 /* */ 242 /* typedef of 2d points for screen coordinates */ 243 /* */ 244 /****************************************************************************/ 245 246 struct coord_point 247 { 248 DOUBLE x; 249 DOUBLE y; 250 }; 251 252 typedef struct coord_point COORD_POINT; 253 254 /****************************************************************************/ 255 /* */ 256 /* definition of exported global variables */ 257 /* */ 258 /****************************************************************************/ 259 260 extern const DOUBLE unit_vec[DIM][DIM]; 261 262 /****************************************************************************/ 263 /* */ 264 /* function declarations */ 265 /* */ 266 /****************************************************************************/ 267 268 /* general routines */ 269 bool PointInPolygon (const COORD_POINT *Points, INT n, COORD_POINT Point); 270 271 /* 2D routines */ 272 DOUBLE vp (const DOUBLE x1, const DOUBLE y1, const DOUBLE x2, const DOUBLE y2); 273 INT V2_Normalize (DOUBLE *a); 274 DOUBLE c_tarea (const DOUBLE *x0, const DOUBLE *x1, const DOUBLE *x2); 275 DOUBLE c_qarea (const DOUBLE *x0, const DOUBLE *x1, const DOUBLE *x2, const DOUBLE *x3); 276 DOUBLE V_te (const DOUBLE *x0, const DOUBLE *x1, 277 const DOUBLE *x2, const DOUBLE *x3); 278 DOUBLE V_py (const DOUBLE *x0, const DOUBLE *x1, const DOUBLE *x2, 279 const DOUBLE *x3, const DOUBLE *x4); 280 DOUBLE V_pr (const DOUBLE *x0, const DOUBLE *x1, const DOUBLE *x2, 281 const DOUBLE *x3, const DOUBLE *x4, const DOUBLE *x5); 282 DOUBLE V_he (const DOUBLE *x0, const DOUBLE *x1, const DOUBLE *x2, const DOUBLE *x3, 283 const DOUBLE *x4, const DOUBLE *x5, const DOUBLE *x6, const DOUBLE *x7); 284 285 286 /* 3D routines */ 287 INT M3_Invert (DOUBLE *Inverse, const DOUBLE *Matrix); 288 INT V3_Normalize (DOUBLE *a); 289 INT V3_Project (const DOUBLE *a, const DOUBLE *b, DOUBLE *r); 290 291 292 /* volume calculations*/ 293 DOUBLE GeneralElementVolume (INT tag, DOUBLE *x_co[]); 294 DOUBLE ElementVolume (const ELEMENT *elem); 295 296 END_UGDIM_NAMESPACE 297 298 #endif 299