1 /*! \file shapes.h 2 * \ingroup gm 3 */ 4 5 6 /****************************************************************************/ 7 /* */ 8 /* File: shapes.h */ 9 /* */ 10 /* Purpose: header file for shape functions */ 11 /* */ 12 /* Author: Klaus Johannsen */ 13 /* Interdisziplinaeres Zentrum fuer Wissenschaftliches Rechnen */ 14 /* Universitaet Heidelberg */ 15 /* Im Neuenheimer Feld 368 */ 16 /* 6900 Heidelberg */ 17 /* internet: ug@ica3.uni-stuttgart.de */ 18 /* */ 19 /* History: 28.11.95 begin, ug version 3.1 */ 20 /* */ 21 /* Remarks: */ 22 /* */ 23 /****************************************************************************/ 24 25 26 /****************************************************************************/ 27 /* */ 28 /* auto include mechanism and other include files */ 29 /* */ 30 /****************************************************************************/ 31 32 #ifndef __SHAPES__ 33 #define __SHAPES__ 34 35 #include <cmath> 36 37 #include "gm.h" 38 #include "evm.h" 39 40 #include <dune/uggrid/low/namespace.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 #ifdef __TWODIM__ 55 56 #define CORNER_COORDINATES_TRIANGLE(e,n,x) \ 57 {(n)=3; \ 58 (x)[0]=CVECT(MYVERTEX(CORNER((e),0))); \ 59 (x)[1]=CVECT(MYVERTEX(CORNER((e),1))); \ 60 (x)[2]=CVECT(MYVERTEX(CORNER((e),2)));} 61 62 #define CORNER_COORDINATES_QUADRILATERAL(e,n,x) \ 63 {(n)=4; \ 64 (x)[0]=CVECT(MYVERTEX(CORNER((e),0))); \ 65 (x)[1]=CVECT(MYVERTEX(CORNER((e),1))); \ 66 (x)[2]=CVECT(MYVERTEX(CORNER((e),2))); \ 67 (x)[3]=CVECT(MYVERTEX(CORNER((e),3)));} 68 69 #define CORNER_COORDINATES(e,n,x) \ 70 {if (TAG((e))==TRIANGLE) \ 71 CORNER_COORDINATES_TRIANGLE((e),(n),(x)) \ 72 else CORNER_COORDINATES_QUADRILATERAL((e),(n),(x))} 73 74 #endif /* __TWODIM__ */ 75 76 77 #define LOCAL_TO_GLOBAL_TRIANGLE(x,local,global) \ 78 {(global)[0] = (1.0-(local)[0]-(local)[1])*(x)[0][0] \ 79 +(local)[0]*(x)[1][0] + (local)[1]*(x)[2][0]; \ 80 (global)[1] = (1.0-(local)[0]-(local)[1])*(x)[0][1] \ 81 +(local)[0]*(x)[1][1] + (local)[1]*(x)[2][1]; } 82 83 #define LOCAL_TO_GLOBAL_QUADRILATERAL(x,local,global) \ 84 {(global)[0] = (1.0-(local)[0])*(1.0-(local)[1])*(x)[0][0] \ 85 + (local)[0]*(1.0-(local)[1])*(x)[1][0] \ 86 + (local)[0]*(local)[1]*(x)[2][0] \ 87 + (1.0-(local)[0])*(local)[1]*(x)[3][0]; \ 88 (global)[1] = (1.0-(local)[0])*(1.0-(local)[1])*(x)[0][1] \ 89 + (local)[0]*(1.0-(local)[1])*(x)[1][1] \ 90 + (local)[0]*(local)[1]*(x)[2][1] \ 91 + (1.0-(local)[0])*(local)[1]*(x)[3][1];} 92 93 #define LOCAL_TO_GLOBAL_2D(n,x,local,global) \ 94 {if ((n) == 3) LOCAL_TO_GLOBAL_TRIANGLE((x),(local),(global)) \ 95 else if ((n) == 4) LOCAL_TO_GLOBAL_QUADRILATERAL((x),(local),(global))} 96 97 #define AREA_OF_TRIANGLE(x,area) {DOUBLE detJ; DOUBLE_VECTOR a,b; \ 98 V2_SUBTRACT((x)[1],(x)[0],a); \ 99 V2_SUBTRACT((x)[2],(x)[0],b); \ 100 V2_VECTOR_PRODUCT(a,b,detJ); \ 101 (area) = std::abs(detJ) * 0.5;} 102 103 #define AREA_OF_QUADRILATERAL(x,area) {DOUBLE detJ,ar; DOUBLE_VECTOR a,b; \ 104 V2_SUBTRACT((x)[1],(x)[0],a); \ 105 V2_SUBTRACT((x)[2],(x)[0],b); \ 106 V2_VECTOR_PRODUCT(a,b,detJ); \ 107 ar = std::abs(detJ) * 0.5; \ 108 V2_SUBTRACT((x)[3],(x)[0],a); \ 109 V2_VECTOR_PRODUCT(a,b,detJ); \ 110 (area) = std::abs(detJ) * 0.5 + ar;} 111 112 #define AREA_OF_ELEMENT_2D(n,x,area) \ 113 {if ((n) == 3) AREA_OF_TRIANGLE((x),(area)) \ 114 else if ((n) == 4) AREA_OF_QUADRILATERAL((x),(area)) } 115 116 #define TRANSFORMATION_OF_TRIANGLE(x,M) \ 117 { V2_SUBTRACT((x)[1],(x)[0],(M)[0]); \ 118 V2_SUBTRACT((x)[2],(x)[0],(M)[1]); } 119 120 #define TRANSFORMATION_OF_QUADRILATERAL(x,local,M) \ 121 { DOUBLE a; \ 122 a = 1.0 - (local)[1]; \ 123 (M)[0][0] = a*((x)[1][0]-x[0][0])+(local)[1]*((x)[2][0]-(x)[3][0]); \ 124 (M)[0][1] = a*((x)[1][1]-x[0][1])+(local)[1]*((x)[2][1]-(x)[3][1]); \ 125 a = 1.0 - (local)[0]; \ 126 (M)[1][0] = a*((x)[3][0]-x[0][0])+(local)[0]*((x)[2][0]-(x)[1][0]); \ 127 (M)[1][1] = a*((x)[3][1]-x[0][1])+(local)[0]*((x)[2][1]-(x)[1][1]); } 128 129 #define TRANSFORMATION_2D(n,x,local,M) \ 130 {if ((n) == 3) {TRANSFORMATION_OF_TRIANGLE((x),(M));} \ 131 else TRANSFORMATION_OF_QUADRILATERAL((x),(local),(M)); } 132 133 134 #ifdef __THREEDIM__ 135 #define CORNER_COORDINATES_TETRAHEDRON(e,n,x) \ 136 {(n) = 4; \ 137 (x)[0]=CVECT(MYVERTEX(CORNER((e),0))); \ 138 (x)[1]=CVECT(MYVERTEX(CORNER((e),1))); \ 139 (x)[2]=CVECT(MYVERTEX(CORNER((e),2))); \ 140 (x)[3]=CVECT(MYVERTEX(CORNER((e),3)));} 141 142 #define CORNER_COORDINATES_PYRAMID(e,n,x) \ 143 {(n) = 5; \ 144 (x)[0]=CVECT(MYVERTEX(CORNER((e),0))); \ 145 (x)[1]=CVECT(MYVERTEX(CORNER((e),1))); \ 146 (x)[2]=CVECT(MYVERTEX(CORNER((e),2))); \ 147 (x)[3]=CVECT(MYVERTEX(CORNER((e),3))); \ 148 (x)[4]=CVECT(MYVERTEX(CORNER((e),4)));} 149 150 #define CORNER_COORDINATES_PRISM(e,n,x) \ 151 {(n) = 6; \ 152 (x)[0]=CVECT(MYVERTEX(CORNER((e),0))); \ 153 (x)[1]=CVECT(MYVERTEX(CORNER((e),1))); \ 154 (x)[2]=CVECT(MYVERTEX(CORNER((e),2))); \ 155 (x)[3]=CVECT(MYVERTEX(CORNER((e),3))); \ 156 (x)[4]=CVECT(MYVERTEX(CORNER((e),4))); \ 157 (x)[5]=CVECT(MYVERTEX(CORNER((e),5)));} 158 159 #define CORNER_COORDINATES_HEXAHEDRON(e,n,x) \ 160 {(n) = 8; \ 161 (x)[0]=CVECT(MYVERTEX(CORNER((e),0))); \ 162 (x)[1]=CVECT(MYVERTEX(CORNER((e),1))); \ 163 (x)[2]=CVECT(MYVERTEX(CORNER((e),2))); \ 164 (x)[3]=CVECT(MYVERTEX(CORNER((e),3))); \ 165 (x)[4]=CVECT(MYVERTEX(CORNER((e),4))); \ 166 (x)[5]=CVECT(MYVERTEX(CORNER((e),5))); \ 167 (x)[6]=CVECT(MYVERTEX(CORNER((e),6))); \ 168 (x)[7]=CVECT(MYVERTEX(CORNER((e),7)));} 169 170 #define CORNER_COORDINATES(e,n,x) \ 171 {if (TAG((e))==TETRAHEDRON) CORNER_COORDINATES_TETRAHEDRON((e),(n),(x))\ 172 else if (TAG((e))==PYRAMID) CORNER_COORDINATES_PYRAMID((e),(n),(x)) \ 173 else if (TAG((e))==PRISM) CORNER_COORDINATES_PRISM((e),(n),(x)) \ 174 else CORNER_COORDINATES_HEXAHEDRON((e),(n),(x))} 175 176 #endif /* __THREEDIM__ */ 177 178 179 #define LOCAL_TO_GLOBAL_TETRAHEDRON(x,local,global) \ 180 {(global)[0] = (1.0-(local)[0]-(local)[1]-(local)[2])*(x)[0][0] \ 181 +(local)[0]*(x)[1][0] + (local)[1]*(x)[2][0] + (local)[2]*(x)[3][0]; \ 182 (global)[1] = (1.0-(local)[0]-(local)[1]-(local)[2])*(x)[0][1] \ 183 +(local)[0]*(x)[1][1] + (local)[1]*(x)[2][1] + (local)[2]*(x)[3][1]; \ 184 (global)[2] = (1.0-(local)[0]-(local)[1]-(local)[2])*(x)[0][2] \ 185 +(local)[0]*(x)[1][2] + (local)[1]*(x)[2][2] + (local)[2]*(x)[3][2]; } 186 187 #define LOCAL_TO_GLOBAL_PYRAMID(x,local,global) \ 188 {DOUBLE a,b,a0,a1,a2,a3; \ 189 a = 1.0 - (local)[0]; \ 190 b = 1.0 - (local)[1]; \ 191 if ((local)[0] > (local)[1]) { \ 192 a0 = a * b - (local)[2] * b; \ 193 a1 = (local)[0] * b - (local)[2]*(local)[1]; \ 194 a2 = (local)[0] * (local)[1] + (local)[2]*(local)[1]; \ 195 a3 = a * (local)[1] - (local)[2] * (local)[1]; \ 196 (global)[0] = \ 197 a0*(x)[0][0]+a1*(x)[1][0]+a2*(x)[2][0]+a3*(x)[3][0]+(local)[2]*(x)[4][0]; \ 198 (global)[1] = \ 199 a0*(x)[0][1]+a1*(x)[1][1]+a2*(x)[2][1]+a3*(x)[3][1]+(local)[2]*(x)[4][1]; \ 200 (global)[2] = \ 201 a0*(x)[0][2]+a1*(x)[1][2]+a2*(x)[2][2]+a3*(x)[3][2]+(local)[2]*(x)[4][2];}\ 202 else { \ 203 a0 = a * b - (local)[2] * a; \ 204 a1 = (local)[0] * b - (local)[2]*(local)[0]; \ 205 a2 = (local)[0] * (local)[1] + (local)[2]*(local)[0]; \ 206 a3 = a * (local)[1] - (local)[2] * (local)[0]; \ 207 (global)[0] = \ 208 a0*(x)[0][0]+a1*(x)[1][0]+a2*(x)[2][0]+a3*(x)[3][0]+(local)[2]*(x)[4][0]; \ 209 (global)[1] = \ 210 a0*(x)[0][1]+a1*(x)[1][1]+a2*(x)[2][1]+a3*(x)[3][1]+(local)[2]*(x)[4][1]; \ 211 (global)[2] = \ 212 a0*(x)[0][2]+a1*(x)[1][2]+a2*(x)[2][2]+a3*(x)[3][2]+(local)[2]*(x)[4][2];}} 213 214 #define LOCAL_TO_GLOBAL_PRISM(x,local,global) \ 215 {DOUBLE a,b,a0,a1,a2,a3,a4,a5; \ 216 a = 1.0 - (local)[0] - (local)[1]; \ 217 b = 1.0 - (local)[2]; \ 218 a0 = a * b; \ 219 a1 = (local)[0] * b; \ 220 a2 = (local)[1] * b; \ 221 a3 = a * (local)[2]; \ 222 a4 = (local)[0] * (local)[2]; \ 223 a5 = (local)[1] * (local)[2]; \ 224 (global)[0] = \ 225 a0*(x)[0][0]+a1*(x)[1][0]+a2*(x)[2][0]+a3*(x)[3][0]+ \ 226 a4*(x)[4][0]+a5*(x)[5][0]; \ 227 (global)[1] = \ 228 a0*(x)[0][1]+a1*(x)[1][1]+a2*(x)[2][1]+a3*(x)[3][1]+ \ 229 a4*(x)[4][1]+a5*(x)[5][1]; \ 230 (global)[2] = \ 231 a0*(x)[0][2]+a1*(x)[1][2]+a2*(x)[2][2]+a3*(x)[3][2]+ \ 232 a4*(x)[4][2]+a5*(x)[5][2]; } 233 234 #define LOCAL_TO_GLOBAL_HEXAHEDRON(x,local,global) \ 235 {DOUBLE a,b,c,a0,a1,a2,a3,a4,a5,a6,a7; \ 236 a = 1.0 - (local)[0]; \ 237 b = 1.0 - (local)[1]; \ 238 c = 1.0 - (local)[2]; \ 239 a0 = a * b * c; \ 240 a1 = (local)[0] * b * c; \ 241 a2 = (local)[0] * (local)[1] * c; \ 242 a3 = a * (local)[1] * c; \ 243 a4 = a * b * (local)[2]; \ 244 a5 = (local)[0] * b * (local)[2]; \ 245 a6 = (local)[0] * (local)[1] * (local)[2]; \ 246 a7 = a * (local)[1] * (local)[2]; \ 247 (global)[0] = \ 248 a0*(x)[0][0]+a1*(x)[1][0]+a2*(x)[2][0]+a3*(x)[3][0]+ \ 249 a4*(x)[4][0]+a5*(x)[5][0]+a6*(x)[6][0]+a7*(x)[7][0]; \ 250 (global)[1] = \ 251 a0*(x)[0][1]+a1*(x)[1][1]+a2*(x)[2][1]+a3*(x)[3][1]+ \ 252 a4*(x)[4][1]+a5*(x)[5][1]+a6*(x)[6][1]+a7*(x)[7][1]; \ 253 (global)[2] = \ 254 a0*(x)[0][2]+a1*(x)[1][2]+a2*(x)[2][2]+a3*(x)[3][2]+ \ 255 a4*(x)[4][2]+a5*(x)[5][2]+a6*(x)[6][2]+a7*(x)[7][2]; } 256 257 #define LOCAL_TO_GLOBAL_3D(n,x,local,global) \ 258 {if ((n) == 4) LOCAL_TO_GLOBAL_TETRAHEDRON((x),(local),(global)) \ 259 else if ((n) == 5) LOCAL_TO_GLOBAL_PYRAMID((x),(local),(global)) \ 260 else if ((n) == 6) LOCAL_TO_GLOBAL_PRISM((x),(local),(global)) \ 261 else if ((n) == 8) LOCAL_TO_GLOBAL_HEXAHEDRON((x),(local),(global))} 262 263 #define AREA_OF_TETRAHEDRON(x,area) {DOUBLE detJ; DOUBLE_VECTOR a,b,c; \ 264 V3_SUBTRACT((x)[1],(x)[0],a); \ 265 V3_SUBTRACT((x)[2],(x)[0],b); \ 266 V3_VECTOR_PRODUCT(a,b,c); \ 267 V3_SUBTRACT((x)[3],(x)[0],a); \ 268 V3_SCALAR_PRODUCT(a,c,detJ); \ 269 (area) = std::abs(detJ)/6.0;} 270 271 #define AREA_OF_PYRAMID(x,area) {DOUBLE detJ,ar; DOUBLE_VECTOR a,b,c,d;\ 272 V3_SUBTRACT((x)[1],(x)[0],a); \ 273 V3_SUBTRACT((x)[2],(x)[0],b); \ 274 V3_VECTOR_PRODUCT(a,b,c); \ 275 V3_SUBTRACT((x)[4],(x)[0],d); \ 276 V3_SCALAR_PRODUCT(c,d,detJ); \ 277 ar = std::abs(detJ)/6.0; \ 278 V3_SUBTRACT((x)[3],(x)[0],a); \ 279 V3_VECTOR_PRODUCT(a,b,c); \ 280 V3_SCALAR_PRODUCT(c,d,detJ); \ 281 (area) = std::abs(detJ)/6.0 + ar;} 282 283 #define AREA_OF_PRISM(x,area) {DOUBLE detJ,ar; DOUBLE_VECTOR a,b,c; \ 284 V3_SUBTRACT((x)[1],(x)[0],a); \ 285 V3_SUBTRACT((x)[2],(x)[0],b); \ 286 V3_VECTOR_PRODUCT(a,b,c); \ 287 V3_SUBTRACT((x)[3],(x)[0],a); \ 288 V3_SCALAR_PRODUCT(a,c,detJ); \ 289 ar = std::abs(detJ)/6.0; \ 290 V3_SUBTRACT((x)[2],(x)[1],a); \ 291 V3_SUBTRACT((x)[3],(x)[1],b); \ 292 V3_VECTOR_PRODUCT(a,b,c); \ 293 V3_SUBTRACT((x)[4],(x)[1],a); \ 294 V3_SCALAR_PRODUCT(a,c,detJ); \ 295 ar += std::abs(detJ)/6.0; \ 296 V3_SUBTRACT((x)[2],(x)[5],a); \ 297 V3_SUBTRACT((x)[3],(x)[5],b); \ 298 V3_VECTOR_PRODUCT(a,b,c); \ 299 V3_SUBTRACT((x)[4],(x)[5],a); \ 300 V3_SCALAR_PRODUCT(a,c,detJ); \ 301 (area) = std::abs(detJ)/6.0 + ar;} 302 303 #define AREA_OF_HEXAHEDRON(x,area) {DOUBLE detJ,ar; DOUBLE_VECTOR a,b,c; \ 304 V3_SUBTRACT((x)[1],(x)[0],a); \ 305 V3_SUBTRACT((x)[2],(x)[0],b); \ 306 V3_VECTOR_PRODUCT(a,b,c); \ 307 V3_SUBTRACT((x)[5],(x)[0],a); \ 308 V3_SCALAR_PRODUCT(a,c,detJ); \ 309 ar = std::abs(detJ)/6.0; \ 310 V3_SUBTRACT((x)[2],(x)[0],a); \ 311 V3_SUBTRACT((x)[5],(x)[0],b); \ 312 V3_VECTOR_PRODUCT(a,b,c); \ 313 V3_SUBTRACT((x)[6],(x)[0],a); \ 314 V3_SCALAR_PRODUCT(a,c,detJ); \ 315 ar += std::abs(detJ)/6.0; \ 316 V3_SUBTRACT((x)[4],(x)[0],a); \ 317 V3_SUBTRACT((x)[5],(x)[0],b); \ 318 V3_VECTOR_PRODUCT(a,b,c); \ 319 V3_SUBTRACT((x)[6],(x)[0],a); \ 320 V3_SCALAR_PRODUCT(a,c,detJ); \ 321 ar += std::abs(detJ)/6.0; \ 322 V3_SUBTRACT((x)[2],(x)[0],a); \ 323 V3_SUBTRACT((x)[3],(x)[0],b); \ 324 V3_VECTOR_PRODUCT(a,b,c); \ 325 V3_SUBTRACT((x)[6],(x)[0],a); \ 326 V3_SCALAR_PRODUCT(a,c,detJ); \ 327 ar += std::abs(detJ)/6.0; \ 328 V3_SUBTRACT((x)[3],(x)[0],a); \ 329 V3_SUBTRACT((x)[4],(x)[0],b); \ 330 V3_VECTOR_PRODUCT(a,b,c); \ 331 V3_SUBTRACT((x)[6],(x)[0],a); \ 332 V3_SCALAR_PRODUCT(a,c,detJ); \ 333 ar += std::abs(detJ)/6.0; \ 334 V3_SUBTRACT((x)[3],(x)[7],a); \ 335 V3_SUBTRACT((x)[4],(x)[7],b); \ 336 V3_VECTOR_PRODUCT(a,b,c); \ 337 V3_SUBTRACT((x)[6],(x)[7],a); \ 338 V3_SCALAR_PRODUCT(a,c,detJ); \ 339 (area) = std::abs(detJ)/6.0 + ar;} 340 341 #define AREA_OF_ELEMENT_3D(n,x,area) \ 342 {if ((n) == 4) {AREA_OF_TETRAHEDRON((x),(area));} \ 343 else if ((n) == 5) {AREA_OF_PYRAMID((x),(area));} \ 344 else if ((n) == 6) {AREA_OF_PRISM((x),(area));} \ 345 else if ((n) == 8) {AREA_OF_HEXAHEDRON((x),(area));}} 346 347 #define TRANSFORMATION_OF_TETRAHEDRON(x,M) \ 348 { V3_SUBTRACT((x)[1],(x)[0],(M)[0]); \ 349 V3_SUBTRACT((x)[2],(x)[0],(M)[1]); \ 350 V3_SUBTRACT((x)[3],(x)[0],(M)[2]);} 351 352 #define TRANSFORMATION_OF_PYRAMID(x,local,M) \ 353 { DOUBLE a,b,c; \ 354 a = (x)[0][0]-(x)[1][0]+(x)[2][0]-(x)[3][0]; \ 355 b = (x)[0][1]-(x)[1][1]+(x)[2][1]-(x)[3][1]; \ 356 c = (x)[0][2]-(x)[1][2]+(x)[2][2]-(x)[3][2]; \ 357 if ((local)[0] > (local)[1]) { \ 358 (M)[0][0] = (x)[1][0]-(x)[0][0]+(local)[1]*a; \ 359 (M)[0][1] = (x)[1][1]-(x)[0][1]+(local)[1]*b; \ 360 (M)[0][2] = (x)[1][2]-(x)[0][2]+(local)[1]*c; \ 361 (M)[1][0] = (x)[3][0]-(x)[0][0]+((local)[0]+(local)[2])*a; \ 362 (M)[1][1] = (x)[3][1]-(x)[0][1]+((local)[0]+(local)[2])*b; \ 363 (M)[1][2] = (x)[3][2]-(x)[0][2]+((local)[0]+(local)[2])*c; \ 364 (M)[2][0] = (x)[4][0]-(x)[0][0]+(local)[1]*a; \ 365 (M)[2][1] = (x)[4][1]-(x)[0][1]+(local)[1]*b; \ 366 (M)[2][2] = (x)[4][2]-(x)[0][2]+(local)[1]*c;} \ 367 else { \ 368 (M)[0][0] = (x)[1][0]-(x)[0][0]+((local)[1]+(local)[2])*a; \ 369 (M)[0][1] = (x)[1][1]-(x)[0][1]+((local)[1]+(local)[2])*b; \ 370 (M)[0][2] = (x)[1][2]-(x)[0][2]+((local)[1]+(local)[2])*c; \ 371 (M)[1][0] = (x)[3][0]-(x)[0][0]+(local)[0]*a; \ 372 (M)[1][1] = (x)[3][1]-(x)[0][1]+(local)[0]*b; \ 373 (M)[1][2] = (x)[3][2]-(x)[0][2]+(local)[0]*c; \ 374 (M)[2][0] = (x)[4][0]-(x)[0][0]+(local)[0]*a; \ 375 (M)[2][1] = (x)[4][1]-(x)[0][1]+(local)[0]*b; \ 376 (M)[2][2] = (x)[4][2]-(x)[0][2]+(local)[0]*c;} } 377 378 #define TRANSFORMATION_OF_PRISM(x,local,M) \ 379 { DOUBLE a0,a1,a2,b0,b1,b2; \ 380 a0 = (x)[0][0]-(x)[1][0]-(x)[3][0]+(x)[4][0]; \ 381 a1 = (x)[0][1]-(x)[1][1]-(x)[3][1]+(x)[4][1]; \ 382 a2 = (x)[0][2]-(x)[1][2]-(x)[3][2]+(x)[4][2]; \ 383 b0 = (x)[0][0]-(x)[2][0]-(x)[3][0]+(x)[5][0]; \ 384 b1 = (x)[0][1]-(x)[2][1]-(x)[3][1]+(x)[5][1]; \ 385 b2 = (x)[0][2]-(x)[2][2]-(x)[3][2]+(x)[5][2]; \ 386 (M)[0][0] = (x)[1][0]-(x)[0][0]+(local)[2]*a0; \ 387 (M)[0][1] = (x)[1][1]-(x)[0][1]+(local)[2]*a1; \ 388 (M)[0][2] = (x)[1][2]-(x)[0][2]+(local)[2]*a2; \ 389 (M)[1][0] = (x)[2][0]-(x)[0][0]+(local)[2]*b0; \ 390 (M)[1][1] = (x)[2][1]-(x)[0][1]+(local)[2]*b1; \ 391 (M)[1][2] = (x)[2][2]-(x)[0][2]+(local)[2]*b2; \ 392 (M)[2][0] = (x)[3][0]-(x)[0][0]+(local)[0]*a0+(local)[1]*b0; \ 393 (M)[2][1] = (x)[3][1]-(x)[0][1]+(local)[0]*a1+(local)[1]*b1; \ 394 (M)[2][2] = (x)[3][2]-(x)[0][2]+(local)[0]*a2+(local)[1]*b2;} 395 396 #define TRANSFORMATION_OF_HEXAHEDRON(x,local,M) \ 397 { DOUBLE a,b,c,a0,a1,a2,a3; \ 398 a = 1.0 - (local)[0]; \ 399 b = 1.0 - (local)[1]; \ 400 c = 1.0 - (local)[2]; \ 401 a0 = b * c; \ 402 a1 = (local)[1] * c; \ 403 a2 = (local)[1] * (local)[2]; \ 404 a3 = b * (local)[2]; \ 405 (M)[0][0] = a0*((x)[1][0]-x[0][0])+a1*((x)[2][0]-(x)[3][0]) \ 406 + a2*((x)[6][0]-x[7][0])+a3*((x)[5][0]-(x)[4][0]); \ 407 (M)[0][1] = a0*((x)[1][1]-x[0][1])+a1*((x)[2][1]-(x)[3][1]) \ 408 + a2*((x)[6][1]-x[7][1])+a3*((x)[5][1]-(x)[4][1]); \ 409 (M)[0][2] = a0*((x)[1][2]-x[0][2])+a1*((x)[2][2]-(x)[3][2]) \ 410 + a2*((x)[6][2]-x[7][2])+a3*((x)[5][2]-(x)[4][2]); \ 411 a0 = a * c; \ 412 a1 = (local)[0] * c; \ 413 a2 = (local)[0] * (local)[2]; \ 414 a3 = a * (local)[2]; \ 415 (M)[1][0] = a0*((x)[3][0]-x[0][0])+a1*((x)[2][0]-(x)[1][0]) \ 416 + a2*((x)[6][0]-x[5][0])+a3*((x)[7][0]-(x)[4][0]); \ 417 (M)[1][1] = a0*((x)[3][1]-x[0][1])+a1*((x)[2][1]-(x)[1][1]) \ 418 + a2*((x)[6][1]-x[5][1])+a3*((x)[7][1]-(x)[4][1]); \ 419 (M)[1][2] = a0*((x)[3][2]-x[0][2])+a1*((x)[2][2]-(x)[1][2]) \ 420 + a2*((x)[6][2]-x[5][2])+a3*((x)[7][2]-(x)[4][2]); \ 421 a0 = a * b; \ 422 a1 = (local)[0] * b; \ 423 a2 = (local)[0] * (local)[1]; \ 424 a3 = a * (local)[1]; \ 425 (M)[2][0] = a0*((x)[4][0]-x[0][0])+a1*((x)[5][0]-(x)[1][0]) \ 426 + a2*((x)[6][0]-x[2][0])+a3*((x)[7][0]-(x)[3][0]); \ 427 (M)[2][1] = a0*((x)[4][1]-x[0][1])+a1*((x)[5][1]-(x)[1][1]) \ 428 + a2*((x)[6][1]-x[2][1])+a3*((x)[7][1]-(x)[3][1]); \ 429 (M)[2][2] = a0*((x)[4][2]-x[0][2])+a1*((x)[5][2]-(x)[1][2]) \ 430 + a2*((x)[6][2]-x[2][2])+a3*((x)[7][2]-(x)[3][2]); } 431 432 #define TRANSFORMATION_3D(n,x,local,M) \ 433 {if ((n) == 4) {TRANSFORMATION_OF_TETRAHEDRON((x),(M));} \ 434 else if ((n) == 5) {TRANSFORMATION_OF_PYRAMID((x),(local),(M));} \ 435 else if ((n) == 6) {TRANSFORMATION_OF_PRISM((x),(local),(M));} \ 436 else TRANSFORMATION_OF_HEXAHEDRON((x),(local),(M));} 437 438 439 #ifdef __TWODIM__ 440 #define LOCAL_TO_GLOBAL(n,x,local,global) LOCAL_TO_GLOBAL_2D(n,x,local,global) 441 #define AREA_OF_ELEMENT(n,x,area) AREA_OF_ELEMENT_2D(n,x,area) 442 #define TRANSFORMATION(n,x,local,M) TRANSFORMATION_2D(n,x,local,M) 443 #endif /* __TWODIM__ */ 444 445 #ifdef __THREEDIM__ 446 #define TRANSFORMATION_BND(n,x,local,M) TRANSFORMATION_2D(n,x,local,M) 447 #define LOCAL_TO_GLOBAL_BND(n,x,local,global) LOCAL_TO_GLOBAL_2D(n,x,local,global) 448 #define AREA_OF_ELEMENT_BND(n,x,area) AREA_OF_ELEMENT_2D(n,x,area) 449 #define AREA_OF_REF_BND(n,area) AREA_OF_REF_2D(n,area) 450 451 #define LOCAL_TO_GLOBAL(n,x,local,global) LOCAL_TO_GLOBAL_3D(n,x,local,global) 452 #define AREA_OF_ELEMENT(n,x,area) AREA_OF_ELEMENT_3D(n,x,area) 453 #define TRANSFORMATION(n,x,local,M) TRANSFORMATION_3D(n,x,local,M) 454 #endif /* __THREEDIM__ */ 455 456 #define INVERSE_TRANSFORMATION(n,x,local,Jinv,Jdet) \ 457 { DOUBLE_VECTOR J[DIM]; \ 458 TRANSFORMATION((n),(x),(local),J); \ 459 M_DIM_INVERT(J,(Jinv),(Jdet)); } 460 461 /****************************************************************************/ 462 /* */ 463 /* data structures exported by the corresponding source file */ 464 /* */ 465 /****************************************************************************/ 466 467 /****************************************************************************/ 468 /* */ 469 /* definition of exported global variables */ 470 /* */ 471 /****************************************************************************/ 472 473 /****************************************************************************/ 474 /* */ 475 /* function declarations */ 476 /* */ 477 /****************************************************************************/ 478 479 DOUBLE *LMP (INT n); 480 INT UG_GlobalToLocal (INT n, const DOUBLE **Corners, const DOUBLE *EvalPoint, DOUBLE *LocalCoord); 481 482 #ifdef __THREEDIM__ 483 DOUBLE N (const INT i, const DOUBLE *LocalCoord); 484 INT TetraSideNormals (ELEMENT *theElement, DOUBLE **theCorners, DOUBLE_VECTOR theNormals[MAX_SIDES_OF_ELEM]); 485 INT TetMaxSideAngle (ELEMENT *theElement, const DOUBLE **theCorners, DOUBLE *MaxAngle); 486 INT TetAngleAndLength (ELEMENT *theElement, const DOUBLE **theCorners, DOUBLE *Angle, DOUBLE *Length); 487 #endif 488 489 END_UGDIM_NAMESPACE 490 491 #endif 492