1 /* This file contains the routines which generate hierarchic shape 2 functions for a any (right now hexahedral) element, using the 3 concept of Blend functions and entity level shapefunctions. 4 5 The shapefunction for a mesh entity Mj^dj is defined as 6 7 N = psi(Mj^di,Me^de)*phi(Mj^dj) 8 9 where psi(...) is the blending function and phi(..) is the entity 10 level function which takes care of the polynomial order . 11 */ 12 13 14 15 #include <math.h> 16 17 #include "topo_shapedefs.h" 18 19 int nshp =0; 20 double LP(int j, double x) 21 { 22 /* Bonnet's recursion formula for Legendre Polynomials */ 23 if( j==0) return 1; 24 if( j==1) return x; 25 26 double ply; 27 ply = ((2*j-1)*x*LP(j-1,x)-(j-1)*LP(j-2,x))/j; 28 return ply; 29 30 } 31 /***************************************************************************/ 32 double LPdrv(int j, double x) 33 { 34 if(j == 0) return 0; 35 if(j == 1) return 1; 36 37 double plydrv; 38 plydrv = (2*j -1)*LP(j-1,x)+LPdrv(j-2,x); 39 return plydrv; 40 } 41 42 /***************************************************************************/ 43 double phi(int p, double x) 44 { 45 double PHI; 46 47 PHI = LP(p,x)-LP(p-2,x); 48 PHI = PHI / (2*p-1); 49 /* PHI = PHI/sqrt(2*(2*p-1)); */ 50 return PHI; 51 } 52 53 double phiDrv(int p, double x) 54 { 55 double Phidrv; 56 Phidrv = LP(p-1,x); 57 /* Phidrv = sqrt(0.5*(2*p-1))*LP(p-1,x); */ 58 return Phidrv; 59 } 60 61 62 /*******************************************************************/ 63 /* Construction of Blending functions */ 64 65 /* Line element edge blend and its derivatives */ 66 67 double Line_eB(double xi1) 68 { 69 /* edge parameterization I defined in Saikat's thesis */ 70 return 0.5*(xi1*xi1 - 1); 71 } 72 73 double dLEBdxi1(double xi1){ return xi1; } 74 double dLEBdxi2(double xi1){return 0;} 75 double dLEBdxi3(double xi1){return 0;} 76 77 /* Quadrilateral edge blend and its derivatives */ 78 79 80 double Quad_eB(double xi1, double xi2, int sign) 81 { 82 /* Quad element edge blend, for edge along xi1 */ 83 return 0.25*(xi1*xi1 - 1)*(1+sign*xi2); 84 } 85 86 double dQEBdxi1(double xi1, double xi2, int sign) 87 { return 0.5*xi1*(1+sign*xi2); } 88 89 double dQEBdxi2(double xi1, double xi2, int sign) 90 { return 0.25*sign*(xi1*xi1 - 1); } 91 92 double dQEBdxi3(double xi1, double xi2, int sign) 93 { return 0; } 94 95 96 /* Quadrilateral face blend and its derivatives */ 97 98 double Quad_fB(double xi1, double xi2) 99 { 100 return 0.25*(xi1*xi1 - 1)*(xi2*xi2 - 1); 101 } 102 103 double dQFBdxi1(double xi1, double xi2) 104 { return 0.5*xi1*(xi2*xi2 - 1); } 105 106 double dQFBdxi2(double xi1, double xi2) 107 { return 0.5*xi2*(xi1*xi1 - 1); } 108 109 double dQFBdxi3(double xi1, double xi2) 110 { return 0.0 ; } 111 112 113 /* The hexahedral element edge blend and its derivatives */ 114 115 double Hex_eB(double xi[3], int sign2, int sign3) 116 { 117 /* For edge along xi1 */ 118 119 // return 0.125*(xi[0]*xi[0] - 1)*(1+sign2*xi[1])*(1+sign3*xi[2]); 120 return 0.25*(1+sign2*xi[1])*(1+sign3*xi[2]); 121 } 122 123 double dHEBdxi1(double xi[3], int sign2, int sign3) 124 { 125 /* For edge along xi1 */ 126 127 // return 0.25*xi[0]*(1+sign2*xi[1])*(1+sign3*xi[2]); 128 return 0; 129 } 130 131 double dHEBdxi2(double xi[3], int sign2, int sign3) 132 { 133 134 /* For edge along xi1 */ 135 136 // return 0.125*sign2*(xi[0]*xi[0] - 1)*(1+sign3*xi[2]); 137 138 return 0.25*sign2*(1+sign3*xi[2]); 139 } 140 141 double dHEBdxi3(double xi[3], int sign2, int sign3) 142 { 143 /* For edge along xi1 */ 144 145 return 0.25*sign3*(1+sign2*xi[1]); 146 } 147 148 149 /* The hexahedral face blend and its derivatives */ 150 151 double Hex_fB(double xi[3], int sign3) 152 { 153 /* for face perpendicular to xi3 */ 154 155 return 0.5*(1+sign3*xi[2]); 156 // return 0.125*(xi[0]*xi[0] - 1)*(xi[1]*xi[1] - 1)*(1+sign3*xi[2]); 157 } 158 159 double dHFBdxi1(double xi[3], int sign3) 160 //{ return 0.25*xi[0]*(xi[1]*xi[1] - 1)*(1+sign3*xi[2]);} 161 { return 0;} 162 163 double dHFBdxi2(double xi[3], int sign3) 164 //{ return 0.25*xi[1]*(xi[0]*xi[0] - 1)*(1+sign3*xi[2]);} 165 { return 0; } 166 167 double dHFBdxi3(double xi[3], int sign3) 168 //{ return 0.125*(xi[0]*xi[0] - 1)*(xi[1]*xi[1] - 1)*sign3; } 169 { return 0.5*sign3; } 170 171 172 // The Construction of higher-order blending functions for edge and face of pyramid */ 173 174 /* The pyramid element edge blend and its derivatives */ 175 double Pyr_eB(double xi[3], int sign[3], int k, int m, int along) 176 { 177 double psi; 178 double tmp0 =1/(1-xi[2]); // xi3->xi[2] 179 180 k--; m--; along--; 181 182 if (along==2) { // if actually along=3 183 // For edge along xi3, which bounds the triangular face 184 psi =0.125*(xi[2]*xi[2]-1)*(1+sign[k]*(2*xi[k]*tmp0))*(1+sign[m]*(2*xi[m]*tmp0)); 185 } else { // if actually along=k=1 or 2 186 // For edge along xik, k=1, 2, m=2, 1 which bounds the quadrilater face 187 double tmp1 =2*xi[k]*tmp0; 188 psi =0.125*(tmp1*tmp1-1)*(1+sign[m]*(2*xi[m]*tmp0))*(1-xi[2]); 189 } 190 return psi; 191 } 192 193 double dPeBdxi(double xi[3], int sign[3], int k, int m, int along, int byWhich) 194 { 195 double dPsidxi = 0; 196 double tmp0 =1/(1-xi[2]); // xi3->xi[2] 197 198 // byWhich was decreased by 1. When it comes in as j it is 0 initially 199 // then it becomes -1, wrong. Do not decrease byWhich. 200 k--; m--; along--; 201 if (along==2) { // for edges along xi3 202 if (byWhich==k) 203 dPsidxi =-0.25*(1+xi[2])*sign[k]*(1+sign[m]*(2*xi[m]*tmp0)); 204 else if (byWhich==m) 205 dPsidxi =-0.25*(1+xi[2])*sign[m]*(1+sign[k]*(2*xi[k]*tmp0)); 206 else if (byWhich==2) // actually byWhich =2 207 dPsidxi =0.25*(xi[2]*(1+sign[k]*2*xi[k]*tmp0)*(1+sign[m]*2*xi[m]*tmp0)- 208 (1+xi[2])*tmp0*(sign[k]*xi[k]*(1+sign[m]*2*xi[m]*tmp0)+ 209 sign[m]*xi[m]*(1+sign[k]*2*xi[k]*tmp0))); 210 } else { // for edges along xik with k=1 or 2 211 double tmp1 =2*xi[k]*tmp0; 212 if (byWhich==k) 213 dPsidxi =0.5*tmp1*(1+sign[m]*(2*xi[m]*tmp0)); 214 else if (byWhich==m) 215 dPsidxi =0.25*(tmp1*tmp1-1)*sign[m]; 216 else if (byWhich==2) // actually byWhich =2 217 dPsidxi =0.25*(tmp1*tmp1)*(1+sign[m]*(2*xi[m]*tmp0))-0.125*(tmp1*tmp1-1); 218 } 219 return dPsidxi; 220 } 221 222 /* The pyramid element face blend and its derivatives */ 223 /* Needs to be corrected as we did for the edge blend */ 224 double Pyr_fB (double xi[3], int sign[3], int k, int m, int faceType) 225 { 226 double tmp0 =1/(1-xi[1]); // xi2->xi[1] 227 double psi; 228 229 if (faceType==4) { 230 // for the quadrilater face 231 double tmp1 =2*xi[0]*tmp0; 232 double tmp2 =2*xi[2]*tmp0; 233 234 psi =0.125*(1-tmp1*tmp1)*(1-tmp2*tmp2)*(1-xi[1]); 235 } else { 236 // for the triangular faces with k=1, 3 and m=3, 1 237 k--; m--; 238 double tmp1 =2*xi[m]*tmp0; 239 240 psi =0.125*(1+sign[k]*2*xi[k]*tmp0)*(1-tmp1*tmp1)*(1-xi[1]*xi[1]); 241 } 242 return psi; 243 } 244 245 double dPfBdxi(double xi[3], int sign[3], int k, int m, int faceType, int byWhich) 246 { 247 double dPsidxi = 0; 248 double tmp0 =1/(1-xi[2]); // xi3->xi[2] 249 250 if (faceType==4) { 251 // for the quadrilater face 252 double tmp1 =2*xi[0]*tmp0; 253 double tmp2 =2*xi[2]*tmp0; 254 255 switch (byWhich) { 256 case 1: 257 dPsidxi =-0.5*tmp1*(1-tmp2*tmp2); 258 break; 259 case 2: 260 dPsidxi = -0.125* (1+tmp1*tmp1+tmp2*tmp2-3*tmp1*tmp1*tmp2*tmp2); 261 break; 262 case 3: 263 dPsidxi =-0.5*(1-tmp1*tmp1)*tmp2; 264 break; 265 } 266 } else { // for the triangular face with k=1,3 and m=3,1 267 k--; m--; byWhich--; 268 if (byWhich==k) { 269 double tmp1 =2*xi[m]*tmp0; 270 dPsidxi =0.25*sign[k]*(1-tmp1*tmp1)*(1+xi[1]); 271 } 272 else if (byWhich==m) 273 dPsidxi =-(1+sign[k]*2*xi[k]*tmp0)*(xi[m]*tmp0)*(1+xi[1]); 274 else if (byWhich==1) { // actually byWhich =2 275 double tmp1 =xi[k]*tmp0; 276 double tmp2 =xi[m]*tmp0; 277 dPsidxi =(1+xi[1])*(0.25*sign[k]*tmp1*(1-4*tmp2*tmp2)- 278 (1+2*sign[k]*tmp1)*tmp2*tmp2)- 279 0.25*(1+2*sign[k]*tmp1)*(1-4*tmp2*tmp2)*xi[1]; 280 } 281 } 282 return dPsidxi; 283 } 284 285 286 /* Add further Blending functions and their derivatives before this line */ 287 /*******************************************************************/ 288 289 /* Entity level functions phi(....) */ 290 291 292 int mesh_edge(double xi1, int gOrd[3], int p,double* entfn,double** edrv) 293 { 294 int nem = p-1; 295 // double leb = Line_eB(xi1); 296 // double dlebdxi1 = dLEBdxi1(xi1); 297 298 if(nem > 0){ 299 300 // entfn = new double [nem]; 301 // edrv = new double* [nem]; 302 303 for(int i =0; i< nem; i++) { 304 305 // edrv[i] = new double [3]; 306 307 // entfn[i] = phi(i+2, xi1)/leb; 308 // edrv[i][gOrd[0]] = (leb*phiDrv(i+2,xi1)-phi(i+2,xi1)*dlebdxi1)/leb*leb; 309 // edrv[i][gOrd[1]] = 0.0; 310 // edrv[i][gOrd[2]] = 0.0; 311 312 entfn[i] = phi(i+2, xi1); 313 edrv[i][gOrd[0]] = phiDrv(i+2,xi1); 314 edrv[i][gOrd[1]] = 0.0; 315 edrv[i][gOrd[2]] = 0.0; 316 317 } 318 } 319 320 return nem; 321 } 322 323 int quad_face(double xi[3], int gOrd[3], int p, double* entfn, double** edrv) 324 { 325 int nfm; 326 int a1, a2; 327 int mc=0; /* mode counter*/ 328 // double temp1, temp2; 329 330 double xi1 = xi[0]; 331 double xi2 = xi[1]; 332 333 // double qfb = Quad_fB(xi1,xi2); 334 // double dqfbdxi1 = dQFBdxi1(xi1,xi2); 335 // double dqfbdxi2 = dQFBdxi2(xi1,xi2); 336 337 if(p > 3){ 338 339 nfm = (p-2)*(p-3)/2; 340 // entfn = new double [nfm]; 341 // edrv = new double* [nfm]; 342 // for(int i=0; i< nfm; i++){ 343 // edrv[i]=new double [3]; 344 // } 345 346 for(int ip =3; ip <p+1; ip++){ /* for each p */ 347 348 for(a1 = 2; a1 < ip-1; a1++){ /* a1,a2 = 2....p-2 */ 349 for(a2 = 2; a2 < ip-1; a2++){ /* a1+a2 = p */ 350 if( a1+a1 == ip ){ 351 352 // entfn[mc] = phi(a1,xi1)*phi(a2,xi2)/qfb; 353 354 // temp1 = (phi(a2,xi2)/qfb*qfb); 355 // temp2 = (qfb*phiDrv(a1,xi1)- dqfbdxi1*phi(a1,xi1)); 356 // edrv[mc][gOrd[0]]= temp1*temp2; 357 358 // temp1 = (phi(a1,xi1)/qfb*qfb); 359 // temp2 = (qfb*phiDrv(a2,xi2)- dqfbdxi2*phi(a2,xi2)); 360 // edrv[mc][gOrd[1]]= temp1*temp2; 361 362 // edrv[mc++][gOrd[2]]=0.0; 363 364 entfn[mc] = phi(a1,xi1)*phi(a2,xi2); 365 edrv[mc][gOrd[0]] = phiDrv(a1, xi1)*phi(a2, xi2); 366 edrv[mc][gOrd[1]] = phi(a1, xi1)*phiDrv(a2, xi2); 367 edrv[mc++][gOrd[2]]=0.0; 368 369 } 370 } 371 } 372 } 373 } else nfm =0; 374 return nfm; 375 } 376 377 int hex_regn(double xi[3],int p, double* entfn, double** edrv) 378 { 379 int a1, a2, a3; 380 int nrm, mc=0; 381 382 double xi1 = xi[0]; 383 double xi2 = xi[1]; 384 double xi3 = xi[2]; 385 386 if( p > 5){ 387 388 nrm = (p-3)*(p-4)*(p-5)/6; 389 // entfn = new double [nrm]; 390 // edrv = new double* [nrm]; 391 // for(int i=0; i< nrm; i++){ 392 // edrv[i]=new double [3]; 393 // } 394 395 for(int ip =6; ip< p+1; ip++){ 396 397 for(a1=2; a1 < ip-3; a1++){ 398 for(a2=2; a2 < ip-3; a2++){ 399 for(a3=2; a3 < ip-3; a3++){ 400 if(a1+a2+a3 == ip){ 401 entfn[mc]=phi(a1,xi1)*phi(a2,xi2)*phi(a3,xi3); 402 403 edrv[mc][0]=phiDrv(a1,xi1)*phi(a2,xi2)*phi(a3,xi3); 404 405 edrv[mc][1]=phi(a1,xi1)*phiDrv(a2,xi2)*phi(a3,xi3); 406 407 edrv[mc++][2]=phi(a1,xi1)*phi(a2,xi2)*phiDrv(a3,xi3); 408 } 409 } 410 } 411 } 412 } 413 }else nrm =0; 414 return nrm; 415 } 416 417 418 419 /* hex hierarchic shape function */ 420 int HexShapeAndDrv(int p,double par[3],double N[],double dN[][3]) 421 { 422 int nshp = 0; 423 int tmp1[4]; 424 int a,b; 425 double EdgeBlend,dEBdxi,dEBdeta,dEBdzeta; 426 double arg[3]; 427 int arg2[3]; 428 double* entfn; 429 double** endrv; 430 int num_e_modes, num_f_modes, num_r_modes; 431 432 int** edge[12]; 433 int n[8][3]={{-1,-1,-1},{1,-1,-1},{1,1,-1},{-1,1,-1}, 434 {-1,-1,1},{1,-1,1},{1,1,1},{-1,1,1}}; 435 436 int face[6][4] = {{0,3,2,1},{0,1,5,4},{1,2,6,5}, 437 {0,4,7,3},{2,3,7,6},{4,5,6,7}}; 438 int vrt[4], z; 439 int normal = 0,sign; 440 double FaceBlend, dFBdxi, dFBdeta, dFBdzeta; 441 442 if(p<1) return nshp; 443 444 double xi = par[0]; 445 double eta = par[1]; 446 double zeta = par[2]; 447 448 double xim=1-xi; 449 double etam=1-eta; 450 double zetam=1-zeta; 451 452 double xip=1+xi; 453 double etap=1+eta; 454 double zetap=1+zeta; 455 456 /* Shape functions for the Nodes. 457 * There are eight nodal shapefunctions. These are same as the 458 * standard shape functions used in the eight-noded hexahedral 459 * elements 460 */ 461 462 N[0]= 0.125* xim * etam * zetam ; 463 N[1]= 0.125* xip * etam * zetam ; 464 N[2]= 0.125* xip * etap * zetam ; 465 N[3]= 0.125* xim * etap * zetam ; 466 N[4]= 0.125* xim * etam * zetap ; 467 N[5]= 0.125* xip * etam * zetap ; 468 N[6]= 0.125* xip * etap * zetap ; 469 N[7]= 0.125* xim * etap * zetap ; 470 471 /* Derivative of the above Shape Functions */ 472 473 dN[0][0]=-0.125*etam*zetam; 474 dN[0][1]=-0.125*xim*zetam; 475 dN[0][2]=-0.125*xim*etam; 476 477 dN[1][0]=0.125 * etam * zetam; 478 dN[1][1]=-0.125 * xip * zetam; 479 dN[1][2]=-0.125 * xip * etam; 480 481 dN[2][0]= 0.125 * etap * zetam; 482 dN[2][1] = 0.125 * xip * zetam; 483 dN[2][2] = -0.125 * xip * etap; 484 485 dN[3][0] = -0.125 * etap * zetam; 486 dN[3][1] = 0.125 * xim * zetam; 487 dN[3][2] =-0.125 * xim * etap; 488 489 dN[4][0] = -0.125 * etam * zetap; 490 dN[4][1] = -0.125 * xim * zetap; 491 dN[4][2] = 0.125 * xim * etam; 492 493 494 dN[5][0] = 0.125 * etam * zetap; 495 dN[5][1] =-0.125 * xip * zetap; 496 dN[5][2] = 0.125 * xip * etam; 497 498 dN[6][0] = 0.125 * etap * zetap; 499 dN[6][1] = 0.125 * xip * zetap; 500 dN[6][2] = 0.125 * xip * etap; 501 502 dN[7][0] = -0.125 * etap * zetap ; 503 dN[7][1] = 0.125 * xim * zetap; 504 dN[7][2] = 0.125 * xim * etap; 505 506 nshp = 8; 507 508 if( p > 1) { 509 510 /* Generate Shape Functions for Edge Modes. 511 * For a polynomial Order of p there will be 12*(p-1) 512 * edge modes for the entire element. 513 */ 514 515 /* 516 * edge order description; 517 */ 518 519 for(int y=0;y<12;y++) 520 { 521 edge[y]=new int * [2]; 522 } 523 524 edge[0][0]=n[0];edge[0][1]=n[1]; 525 edge[1][0]=n[1];edge[1][1]=n[2]; 526 edge[2][0]=n[2];edge[2][1]=n[3]; 527 edge[3][0]=n[3];edge[3][1]=n[0]; 528 edge[4][0]=n[4];edge[4][1]=n[5]; 529 edge[5][0]=n[5];edge[5][1]=n[6]; 530 edge[6][0]=n[6];edge[6][1]=n[7]; 531 edge[7][0]=n[7];edge[7][1]=n[4]; 532 edge[8][0]=n[0];edge[8][1]=n[4]; 533 edge[9][0]=n[1];edge[9][1]=n[5]; 534 edge[10][0]=n[2];edge[10][1]=n[6]; 535 edge[11][0]=n[3];edge[11][1]=n[7]; 536 537 538 int nem = p-1; 539 540 for(int e=0; e < 12; e++){ 541 542 for(z=0;z<3;z++){ 543 if(!(edge[e][0][z]==edge[e][1][z])) 544 tmp1[3]=z; 545 tmp1[z]=edge[e][1][z]; 546 } 547 548 if((arg2[0] = tmp1[3]) == 0) { arg2[1]=1;arg2[2]=2 ;} 549 else if(tmp1[3]==1) { arg2[1]=2;arg2[2]=0;} 550 else { arg2[1]=0;arg2[2]=1;} 551 552 /* arg[0]=par[arg2[0]]*tmp1[arg2[0]]; */ 553 arg[0]=par[arg2[0]]; 554 arg[1]=par[arg2[1]]; 555 arg[2]=par[arg2[2]]; 556 557 a = arg2[1]; 558 b = arg2[2]; 559 560 EdgeBlend = Hex_eB(arg,tmp1[a],tmp1[b]); 561 562 if( tmp1[3] == 0){ 563 dEBdxi = dHEBdxi1(arg,tmp1[a],tmp1[b]); 564 dEBdeta = dHEBdxi2(arg,tmp1[a],tmp1[b]); 565 dEBdzeta = dHEBdxi3(arg,tmp1[a],tmp1[b]); 566 } else if( tmp1[3] == 1 ){ 567 dEBdeta = dHEBdxi1(arg,tmp1[a],tmp1[b]); 568 dEBdzeta = dHEBdxi2(arg,tmp1[a],tmp1[b]); 569 dEBdxi = dHEBdxi3(arg,tmp1[a],tmp1[b]); 570 } else { 571 dEBdzeta = dHEBdxi1(arg,tmp1[a],tmp1[b]); 572 dEBdxi = dHEBdxi2(arg,tmp1[a],tmp1[b]); 573 dEBdeta = dHEBdxi3(arg,tmp1[a],tmp1[b]); 574 } 575 576 entfn = new double [nem]; 577 endrv = new double* [nem]; 578 for(z =0; z < nem ; z++){ 579 endrv[z]=new double [3]; 580 } 581 582 num_e_modes = mesh_edge(arg[0],arg2,p,entfn,endrv); 583 584 for(int nm=0; nm < num_e_modes; nm++){ 585 586 N[nshp]= EdgeBlend * entfn[nm]; 587 dN[nshp][0]= EdgeBlend * endrv[nm][0] + dEBdxi * entfn[nm]; 588 dN[nshp][1]= EdgeBlend * endrv[nm][1] + dEBdeta * entfn[nm]; 589 dN[nshp][2]= EdgeBlend * endrv[nm][2] + dEBdzeta * entfn[nm]; 590 591 nshp++; 592 } 593 594 delete [] entfn; 595 for(z=0; z< num_e_modes; z++){ 596 delete [] endrv[z]; } 597 delete [] endrv; 598 599 } 600 } 601 602 if( p > 3 ) { 603 604 /* Generating Shape functions for the face modes . 605 * For a Hexahedral Element there are 3*(p-2)*(p-3) shape 606 * functions associated with the face modes [ p>=4] 607 */ 608 609 int nfm = (p-2)*(p-3)/2; 610 611 for(int f=0; f< 6; f++) { 612 /* Identitfying the normal to the face */ 613 614 for(int u=0;u<4;u++){ 615 vrt[u]=face[f][u]; 616 } 617 for(int v=0;v<3;v++) 618 { 619 if((n[vrt[0]][v]==n[vrt[1]][v])&&(n[vrt[1]][v]==n[vrt[2]][v])) 620 normal=v; 621 sign=n[vrt[0]][v]; 622 } 623 624 625 arg2[2] = normal; 626 627 if( normal==0) { arg2[0]=1;arg2[1]=2 ;} 628 else if(normal==1) { arg2[0]=2;arg2[1]=0;} 629 else { arg2[0]=0;arg2[1]=1;} 630 631 arg[0] = par[arg2[0]]; 632 arg[1] = par[arg2[1]]; 633 arg[2] = par[arg2[2]]; 634 635 FaceBlend = Hex_fB(arg,sign); 636 637 if( normal == 0){ 638 dFBdxi = dHFBdxi3(arg,sign); 639 dFBdeta = dHFBdxi1(arg,sign); 640 dFBdzeta = dHFBdxi2(arg,sign); 641 } else if( normal == 1 ){ 642 dFBdeta = dHFBdxi3(arg,sign); 643 dFBdzeta = dHFBdxi1(arg,sign); 644 dFBdxi = dHFBdxi2(arg,sign); 645 } else { 646 dFBdzeta = dHFBdxi3(arg,sign); 647 dFBdxi = dHFBdxi1(arg,sign); 648 dFBdeta = dHFBdxi2(arg,sign); 649 } 650 651 entfn = new double [nfm]; 652 endrv = new double* [nfm]; 653 for(z=0; z < nfm ; z++){ 654 endrv[z]=new double [3]; 655 } 656 657 num_f_modes = quad_face(arg,arg2,p,entfn,endrv); 658 659 for(int nm =0; nm < num_f_modes ; nm++){ 660 661 N[nshp] = FaceBlend * entfn[nm]; 662 dN[nshp][0] = FaceBlend * endrv[nm][0] + dFBdxi * entfn[nm]; 663 dN[nshp][1] = FaceBlend * endrv[nm][1] + dFBdeta * entfn[nm]; 664 dN[nshp][2] = FaceBlend * endrv[nm][2] + dFBdzeta * entfn[nm]; 665 666 nshp++; 667 } 668 669 delete [] entfn; 670 for( z=0; z< num_f_modes; z++){ 671 delete [] endrv[z]; } 672 delete [] endrv; 673 } 674 } 675 676 if ( p > 5 ) { 677 678 int nrm = (p-3)*(p-4)*(p-5)/6; 679 680 entfn = new double [nrm]; 681 endrv = new double* [nrm]; 682 for(z =0; z < nrm ; z++){ 683 endrv[z]=new double [3]; 684 } 685 686 num_r_modes = hex_regn(par,p,entfn,endrv); 687 688 for(int nm =0; nm < num_r_modes ; nm++){ 689 N[nshp] = entfn[nm]; 690 for(int d=0; d < 3; d++){ 691 dN[nshp][d] = endrv[nm][d]; 692 } 693 nshp++; 694 } 695 delete [] entfn; 696 for(z=0; z< num_r_modes; z++){ 697 delete [] endrv[z]; } 698 delete [] endrv; 699 700 } 701 return nshp; 702 } 703 704 /* hierarchic wedge element shape function */ 705 706 int WedgeShapeAndDrv( int p, double Inputpar[3], double N[], double dN[][3] ) 707 { 708 int i,j; 709 double FaceBlend, FaceBlendDrv[4]; 710 double FaceEnt; //, FaceEntDrv[2][3]; 711 double par[4]; 712 // int sign; 713 // int temp[4]={0,0,0,0}; 714 double EdgeBlend[9], EdgeBlendDrv[9][3]; 715 // double arg[3]; 716 // double arg2[2]; 717 // int ** edge[9]; 718 double entfn[9]; 719 double endrv[9][3]; 720 // int num_e_modes; 721 // int n[6][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; 722 723 int nshp = 0; 724 par[0] = 1.0 - Inputpar[0]- Inputpar[1]; 725 par[1] = Inputpar[0]; 726 par[2] = Inputpar[1]; 727 par[3] = Inputpar[2]; 728 729 730 if(p<1) return nshp; 731 732 /* there are six nodal shape functions same as the standard 733 shape functions used in the six-noded wedge element */ 734 735 N[0]=0.5*par[0]*(1.0-par[3]); 736 N[1]=0.5*par[1]*(1.0-par[3]); 737 N[2]=0.5*par[2]*(1.0-par[3]); 738 N[3]=0.5*par[0]*(1.0+par[3]); 739 N[4]=0.5*par[1]*(1.0+par[3]); 740 N[5]=0.5*par[2]*(1.0+par[3]); 741 742 /* Derivative of the above Shape functions */ 743 dN[0][0]=-0.25*(1.0-par[3]); 744 dN[0][1]=-0.25*(1.0-par[3]); 745 dN[0][2]=-0.5*par[0]; 746 747 dN[1][0]=0.25*(1.0-par[3]); 748 dN[1][1]=0.0; 749 dN[1][2]=-0.5*par[1]; 750 751 dN[2][0]=0.0; 752 dN[2][1]=0.25*(1.0-par[3]); 753 dN[2][2]=-0.5*par[2]; 754 755 dN[3][0]=-0.25*(1.0+par[3]); 756 dN[3][1]=-0.25*(1.0+par[3]); 757 dN[3][2]=0.5*par[0]; 758 759 dN[4][0]=0.25*(1.0+par[3]); 760 dN[4][1]=0.0; 761 dN[4][2]=0.5*par[1]; 762 763 dN[5][0]=0.0; 764 dN[5][1]=0.25*(1.0+par[3]); 765 dN[5][2]=0.5*par[2]; 766 767 768 nshp = 6; 769 770 // if( p > 1 ){ 771 772 // /* calculate the blending shape functions and their drivitative */ 773 // EdgeBlend[0] = -par[0]*par[1]*(1.0-par[3]); 774 // EdgeBlendDrv[0][0] = -0.5*(par[0]-par[1])*(1.0-par[3]); 775 // EdgeBlendDrv[0][1] = 0.5*par[1]*(1.0-par[3]); 776 // EdgeBlendDrv[0][2] = par[0]*par[1]; 777 778 // EdgeBlend[1] = -par[1]*par[2]*(1.0-par[3]); 779 // EdgeBlendDrv[1][0] = -0.5*par[2]*(1.0-par[3]); 780 // EdgeBlendDrv[1][1] = -0.5*par[1]*(1.0-par[3]); 781 // EdgeBlendDrv[1][2] = par[1]*par[2]; 782 783 // EdgeBlend[2] = -par[0]*par[2]*(1.0-par[3]); 784 // EdgeBlendDrv[2][0] = 0.5*par[2]*(1.0-par[3]); 785 // EdgeBlendDrv[2][1] = -0.5*(par[0]-par[2])*(1.0-par[3]); 786 // EdgeBlendDrv[2][2] = par[0]*par[2]; 787 788 // EdgeBlend[3] = -par[0]*par[1]*(1.0+par[3]); 789 // EdgeBlendDrv[3][0] = -0.5*(par[0]-par[1])*(1.0+par[3]); 790 // EdgeBlendDrv[3][1] = 0.5*par[1]*(1.0+par[3]); 791 // EdgeBlendDrv[3][2] = -par[0]*par[1]; 792 793 // EdgeBlend[4] = -par[1]*par[2]*(1.0+par[3]); 794 // EdgeBlendDrv[4][0] = -0.5*par[2]*(1.0+par[3]); 795 // EdgeBlendDrv[4][1] = -0.5*par[1]*(1.0+par[3]); 796 // EdgeBlendDrv[4][2] = -par[1]*par[2]; 797 798 // EdgeBlend[5] = -par[0]*par[2]*(1.0+par[3]); 799 // EdgeBlendDrv[5][0] = 0.5*par[2]*(1.0+par[3]); 800 // EdgeBlendDrv[5][1] = -0.5*(par[0]-par[2])*(1.0+par[3]); 801 // EdgeBlendDrv[5][2] = -par[0]*par[2]; 802 803 // EdgeBlend[6] = 0.5*par[0]*(par[3]*par[3]-1.0); 804 // EdgeBlendDrv[6][0] = -0.25*(par[3]*par[3] - 1.0); 805 // EdgeBlendDrv[6][1] = -0.25*(par[3]*par[3] - 1.0); 806 // EdgeBlendDrv[6][2] = par[0]*par[3]; 807 808 // EdgeBlend[7] = 0.5*par[1]*(par[3]*par[3]-1.0); 809 // EdgeBlendDrv[7][0] = 0.25*(par[3]*par[3] - 1.0); 810 // EdgeBlendDrv[7][1] = 0.0; 811 // EdgeBlendDrv[7][2] = par[1]*par[3]; 812 813 // EdgeBlend[8] = 0.5*par[2]*(par[3]*par[3]-1.0); 814 // EdgeBlendDrv[8][0] = 0.0; 815 // EdgeBlendDrv[8][1] = 0.25*(par[3]*par[3] - 1.0); 816 // EdgeBlendDrv[8][2] = par[2]*par[3]; 817 818 // /* calculate the mesh entity shape function */ 819 820 // /* for the edge in the two triangle faces */ 821 // for(j=0; j<9; j++){ 822 // /* this is where I think the change in phi is for wedges */ 823 // entfn[j] = 1.0; 824 // /* entfn[j] = sqrt(1.5); */ 825 // for(i=0; i<3; i++) 826 // endrv[j][i] = 0.0; 827 // } 828 829 // // /* for the edge in the perpendicular to triangle faces */ 830 // // for(j=6; j<9; j++){ 831 // // entfn[j] = phi(2, par[3]); 832 // // for(i=0; i<3; i++) 833 // // endrv[j][i] = 0.0; 834 // // } 835 836 // for(j=0; j<9; j++){ 837 // N[nshp] = EdgeBlend[j]*entfn[j]; 838 // dN[nshp][0] = EdgeBlend[j]*endrv[j][0] + EdgeBlendDrv[j][0]*entfn[j]; 839 // dN[nshp][1] = EdgeBlend[j]*endrv[j][1] + EdgeBlendDrv[j][1]*entfn[j]; 840 // dN[nshp][2] = EdgeBlend[j]*endrv[j][2] + EdgeBlendDrv[j][2]*entfn[j]; 841 842 // ++nshp; 843 // } 844 845 // } 846 847 // /* generate the tri face shape fuction */ 848 // if( p > 2 ){ 849 // for(i=0;i<11;i++){ 850 // /* get the face blending functions */ 851 // if(par[3]>0){ 852 // FaceBlend = 0.5*par[0]*par[1]*par[2]*(1.0-par[3]); 853 854 // /* derivate the face blending functions */ 855 // FaceBlendDrv[0]= 0.5*par[1]*par[2]*(1.0-par[3]); 856 // FaceBlendDrv[1]= 0.5*par[0]*par[2]*(1.0-par[3]); 857 // FaceBlendDrv[2]= 0.5*par[0]*par[1]*(1.0-par[3]); 858 // FaceBlendDrv[3]=-0.5*par[0]*par[1]*par[2]; 859 // } 860 // else{ 861 // FaceBlend = 0.5*par[0]*par[1]*par[2]*(1.0+par[3]); 862 863 // /* derivate the face blending functions */ 864 // FaceBlendDrv[0]= 0.5*par[1]*par[2]*(1.0+par[3]); 865 // FaceBlendDrv[1]= 0.5*par[0]*par[2]*(1.0+par[3]); 866 // FaceBlendDrv[2]= 0.5*par[0]*par[1]*(1.0+par[3]); 867 // FaceBlendDrv[3]= 0.5*par[0]*par[1]*par[2]; 868 // } 869 870 // /*calculate the mesh entity function for cubic triangular face */ 871 // FaceEnt = 1.0; 872 873 // /*calculate the shape functions*/ 874 // N[nshp] = FaceBlend*FaceEnt; 875 // dN[nshp][0]=FaceBlendDrv[0]; 876 // dN[nshp][1]=FaceBlendDrv[1]; 877 // dN[nshp][2]=FaceBlendDrv[2]; 878 // dN[nshp][3]=FaceBlendDrv[3]; 879 880 // nshp++; 881 // } 882 // } 883 884 return nshp; 885 } 886 887 /* Pyramid hierarchic shape function up to order 3*/ 888 int PyrShapeAndDrv(int p,double par[3],double N[],double dN[][3]) 889 { 890 int i; 891 double EdgeBlend, EdgeBlendDrv[3]; 892 double EdgeEntity[2], EdgeEntityDrv[2][3]; 893 double FaceBlend, FaceBlendDrv[3]; 894 double FaceEntity[2], FaceEntityDrv[2][3]; 895 896 //dEBdxi,dEBdeta,dEBdzeta; 897 // double arg[3]; 898 // int arg2[3]; 899 // double* entfn; 900 // double** endrv; 901 // int num_e_modes, num_f_modes, num_r_modes; 902 903 int** edge[8]; // total 8 edges 904 int n[5][3]={{-1,-1, -1},{1,-1,-1},{1,1,-1},{-1,1, -1}, {0, 0, 1}}; // vertex coordinates 905 906 // int n[5][3]={{-1,-1, 1},{-1,-1,-1},{1,-1,-1},{1,-1, 1}, {0, 1, 0}}; // vertex coordinates 907 908 // int face[5][4] = {{0,3,2,1},{0,1,4, -1},{1,2,4, -1},{2,3,4,-1},{3,0,4,-1}}; // face vertices 909 // int face[5][4] = {{0,1,2,3},{1,0,4, -1},{0,3,4, -1},{3,2,4,-1},{2,1,4,-1}}; // face vertices 910 // int vrt[4]; 911 // int normal; 912 int sign[3]; 913 914 if(p<1) return nshp; 915 916 917 // when p=1, there are only nodal shapefunctions 918 double zeta = par[2]; 919 double xi = par[0]; 920 double eta = par[1]; 921 922 // double xim=1-xi; 923 // double etam=1-eta; 924 double zetam=1-zeta; 925 926 double zetamap=2.0/zetam; 927 928 // double xip=1+xi; 929 // double etap=1+eta; 930 double zetap=1+zeta; 931 932 double xipp=1+xi*zetamap; 933 double etapp=1+eta*zetamap; 934 935 double ximp=1-xi*zetamap; 936 double etamp=1-eta*zetamap; 937 938 // Shape functions for the Nodes. 939 // There are five nodal shapefunctions. 940 941 N[0]= 0.125* ximp * etamp * zetam ; 942 N[1]= 0.125* xipp * etamp * zetam ; 943 N[2]= 0.125* xipp * etapp * zetam ; 944 N[3]= 0.125* ximp * etapp * zetam ; 945 N[4]= 0.5* zetap; 946 947 // Derivative of the above Shape Functions 948 dN[0][0] =-0.25 * etamp; 949 dN[0][1] =-0.25 * ximp; 950 dN[0][2] = 0.125 * (xi*eta*zetamap*zetamap-1); 951 952 dN[1][0] = 0.25 * etamp; 953 dN[1][1] =-0.25 * xipp; 954 dN[1][2] =-0.125 * (xi*eta*zetamap*zetamap+1); 955 956 dN[2][0] = 0.25 * etapp; 957 dN[2][1] = 0.25 * xipp; 958 dN[2][2] = dN[0][2]; 959 960 dN[3][0] =-0.25 * etapp; 961 dN[3][1] = 0.25 * ximp; 962 dN[3][2] = dN[1][2]; 963 964 dN[4][0] = 0; 965 dN[4][1] = 0; 966 dN[4][2] = 0.5; 967 968 nshp = 5; 969 970 if( p>1) { 971 972 // Generate Shape Functions for Edges. 973 // For a polynomial Order of p there are 8*(p-1) 974 // edge modes for the entire element. 975 976 // allocate spaces for edge order description 977 for(i=0; i<8; i++) 978 edge[i]=new int * [2]; 979 980 // for the edges on the quadrilateral face 981 edge[0][0]=n[0]; edge[0][1]=n[1]; 982 edge[1][0]=n[1]; edge[1][1]=n[2]; 983 edge[2][0]=n[2]; edge[2][1]=n[3]; 984 edge[3][0]=n[3]; edge[3][1]=n[0]; 985 986 // for other edges on triangular faces 987 edge[4][0]=n[0]; edge[4][1]=n[4]; 988 edge[5][0]=n[1]; edge[5][1]=n[4]; 989 edge[6][0]=n[2]; edge[6][1]=n[4]; 990 edge[7][0]=n[3]; edge[7][1]=n[4]; 991 992 993 // calculate the edge blending functions and their derivatives 994 for(i=0; i < 8; i++) { 995 int k, m, along, j; 996 switch (i) { 997 // first four are on quadrilateral face 998 999 // P.N anil fix this 1000 case 0: 1001 k =1; m =2; sign[m-1] =-1; along =k; 1002 break; 1003 case 1: 1004 k =2; m =1; sign[m-1] =1; along =k; 1005 break; 1006 case 2: 1007 k =1; m =2; sign[m-1] =1; along =k; 1008 break; 1009 case 3: 1010 k =2; m =1; sign[m-1] =-1; along =k; 1011 break; 1012 // next four are on triangular faces 1013 case 4: 1014 k =1; m =2; sign[k-1] =-1; sign[m-1] =-1; along =3; 1015 break; 1016 case 5: 1017 k =2; m =1; sign[k-1] =-1; sign[m-1] =1; along =3; 1018 break; 1019 case 6: 1020 k =1; m =2; sign[k-1] =1; sign[m-1] =1; along =3; 1021 break; 1022 case 7: 1023 k =2; m =1; sign[k-1] =1; sign[m-1] =-1; along =3; 1024 } 1025 EdgeBlend =Pyr_eB (par, sign, k, m, along); 1026 for (j=0; j<3; j++) 1027 EdgeBlendDrv[j] =dPeBdxi(par, sign, k, m, along, j); 1028 1029 // calculate the mesh entity shape function 1030 1031 // calculate the edge entity function for p=2 1032 EdgeEntity[0] =1; 1033 EdgeEntityDrv[0][0] =EdgeEntityDrv[0][1] =EdgeEntityDrv[0][2] =0; 1034 1035 if (p>2) { 1036 // calculate the edge entity function for p=3 1037 if (i<4) { 1038 // for the edges of the quadrilateral face with parameterization I 1039 // equation (33) 1040 EdgeEntity[1] =par[k-1]; 1041 for (j=0; j<3; j++) 1042 EdgeEntityDrv[1][j] =(int)(k-1==j); 1043 } else { 1044 // for the edges of the triangular faces with parameterization II 1045 // First of all, we need to represent these edges use xi's 1046 // In our definition, u[0] points to the quadrilateral base 1047 // u[1] points to the top 1048 // then, we get u[0] =(1-xi[1])/2; u[1] =(1+xi[1])/2; 1049 // EdgeEntity[1] =u[1]-u[0] =xi[1]; 1050 EdgeEntity[1] =par[1]; 1051 EdgeEntityDrv[1][1] =1; EdgeEntityDrv[1][0]=EdgeEntityDrv[1][2]=0; 1052 } 1053 } 1054 1055 for(j=0; j < p-1; j++){ 1056 1057 N[nshp]= EdgeBlend * EdgeEntity[j]; 1058 dN[nshp][0]= EdgeBlend * EdgeEntityDrv[j][0] + EdgeBlendDrv[0] * EdgeEntity[j]; 1059 dN[nshp][1]= EdgeBlend * EdgeEntityDrv[j][1] + EdgeBlendDrv[1] * EdgeEntity[j]; 1060 dN[nshp][2]= EdgeBlend * EdgeEntityDrv[j][2] + EdgeBlendDrv[2] * EdgeEntity[j]; 1061 nshp++; 1062 } 1063 } 1064 1065 // calculate the shape function for triangular faces 1066 if (p==3) { 1067 for (i=1; i<5; i++) { 1068 // calculate the face blending functions 1069 int k, m, j; 1070 switch (i) { 1071 case 1: 1072 k =1; m =3; sign[k-1] =-1; 1073 break; 1074 case 2: 1075 k =3; m =1; sign[k-1] =-1; 1076 break; 1077 case 3: 1078 k =1; m =3; sign[k-1] =1; 1079 break; 1080 case 4: 1081 k =3; m =1; sign[k-1] =1; 1082 break; 1083 } 1084 FaceBlend =Pyr_fB (par, sign, k, m, 3); 1085 for (j=0; j<3; j++) 1086 FaceBlendDrv[j] =dPfBdxi(par, sign, k, m, 3, j); 1087 1088 // for p=3 1089 FaceEntity[0] =1; 1090 FaceEntityDrv[0][0] =FaceEntityDrv[0][1] =FaceEntityDrv[0][2] =0; 1091 1092 for(j=0; j < p-2; j++){ 1093 1094 N[nshp]= FaceBlend * EdgeEntity[j]; 1095 dN[nshp][0]= FaceBlend * FaceEntityDrv[j][0] + FaceBlendDrv[0] * FaceEntity[j]; 1096 dN[nshp][1]= FaceBlend * FaceEntityDrv[j][1] + FaceBlendDrv[1] * FaceEntity[j]; 1097 dN[nshp][2]= FaceBlend * FaceEntityDrv[j][2] + FaceBlendDrv[2] * FaceEntity[j]; 1098 nshp++; 1099 } 1100 } 1101 } 1102 } 1103 1104 return nshp; 1105 } 1106 1107 1108 1109 1110