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;
LP(int j,double x)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 /***************************************************************************/
LPdrv(int j,double x)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 /***************************************************************************/
phi(int p,double x)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
phiDrv(int p,double x)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
Line_eB(double xi1)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
dLEBdxi1(double xi1)73 double dLEBdxi1(double xi1){ return xi1; }
dLEBdxi2(double xi1)74 double dLEBdxi2(double xi1){return 0;}
dLEBdxi3(double xi1)75 double dLEBdxi3(double xi1){return 0;}
76
77 /* Quadrilateral edge blend and its derivatives */
78
79
Quad_eB(double xi1,double xi2,int sign)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
dQEBdxi1(double xi1,double xi2,int sign)86 double dQEBdxi1(double xi1, double xi2, int sign)
87 { return 0.5*xi1*(1+sign*xi2); }
88
dQEBdxi2(double xi1,double xi2,int sign)89 double dQEBdxi2(double xi1, double xi2, int sign)
90 { return 0.25*sign*(xi1*xi1 - 1); }
91
dQEBdxi3(double xi1,double xi2,int sign)92 double dQEBdxi3(double xi1, double xi2, int sign)
93 { return 0; }
94
95
96 /* Quadrilateral face blend and its derivatives */
97
Quad_fB(double xi1,double xi2)98 double Quad_fB(double xi1, double xi2)
99 {
100 return 0.25*(xi1*xi1 - 1)*(xi2*xi2 - 1);
101 }
102
dQFBdxi1(double xi1,double xi2)103 double dQFBdxi1(double xi1, double xi2)
104 { return 0.5*xi1*(xi2*xi2 - 1); }
105
dQFBdxi2(double xi1,double xi2)106 double dQFBdxi2(double xi1, double xi2)
107 { return 0.5*xi2*(xi1*xi1 - 1); }
108
dQFBdxi3(double xi1,double xi2)109 double dQFBdxi3(double xi1, double xi2)
110 { return 0.0 ; }
111
112
113 /* The hexahedral element edge blend and its derivatives */
114
Hex_eB(double xi[3],int sign2,int sign3)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
dHEBdxi1(double xi[3],int sign2,int sign3)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
dHEBdxi2(double xi[3],int sign2,int sign3)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
dHEBdxi3(double xi[3],int sign2,int sign3)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
Hex_fB(double xi[3],int sign3)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
dHFBdxi1(double xi[3],int sign3)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
dHFBdxi2(double xi[3],int sign3)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
dHFBdxi3(double xi[3],int sign3)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 */
Pyr_eB(double xi[3],int sign[3],int k,int m,int along)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
dPeBdxi(double xi[3],int sign[3],int k,int m,int along,int byWhich)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 */
Pyr_fB(double xi[3],int sign[3],int k,int m,int faceType)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
dPfBdxi(double xi[3],int sign[3],int k,int m,int faceType,int byWhich)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
mesh_edge(double xi1,int gOrd[3],int p,double * entfn,double ** edrv)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
quad_face(double xi[3],int gOrd[3],int p,double * entfn,double ** edrv)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
hex_regn(double xi[3],int p,double * entfn,double ** edrv)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 */
HexShapeAndDrv(int p,double par[3],double N[],double dN[][3])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
WedgeShapeAndDrv(int p,double Inputpar[3],double N[],double dN[][3])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*/
PyrShapeAndDrv(int p,double par[3],double N[],double dN[][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