00001 #include "difcalc.h"
00002 #include "basefun.h"
00003 #include "ordering.h"
00004 #include <math.h>
00005 #include <stdlib.h>
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 void derivatives_1d (vector &dx,double &jac,vector &x,double xi)
00019 {
00020 long i;
00021 double d1,*nx;
00022
00023 nx = new double [x.n];
00024
00025 switch (x.n){
00026 case 2:{
00027 dx_bf_lin_1d (nx);
00028 break;
00029 }
00030 case 3:{
00031 dx_bf_quad_1d (nx,xi);
00032 break;
00033 }
00034 default:{
00035 fprintf (stderr,"\n\n wrong number of nodes on element in function");
00036 fprintf (stderr,"\n derivatives_1d (%s, line %d).\n",__FILE__,__LINE__);
00037 }
00038 }
00039
00040 d1=0.0;
00041 for (i=0;i<x.n;i++){
00042 d1+=nx[i]*x[i];
00043 }
00044
00045 jac = d1;
00046
00047 for (i=0;i<dx.n;i++){
00048 dx[i]=nx[i]/jac;
00049 }
00050
00051 delete [] nx;
00052 }
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 void jac_1d (double &jac,vector &x,double xi)
00064 {
00065 long i;
00066 double d1,*nx;
00067
00068 nx = new double [x.n];
00069
00070 switch (x.n){
00071 case 2:{
00072 dx_bf_lin_1d (nx);
00073 break;
00074 }
00075 case 3:{
00076 dx_bf_quad_1d (nx,xi);
00077 break;
00078 }
00079 default:{
00080 fprintf (stderr,"\n\n wrong number of nodes on element in function");
00081 fprintf (stderr,"\n jac_1d (%s, line %d).\n",__FILE__,__LINE__);
00082 }
00083 }
00084
00085 d1=0.0;
00086 for (i=0;i<x.n;i++){
00087 d1+=nx[i]*x[i];
00088 }
00089
00090 jac = d1;
00091
00092 delete [] nx;
00093 }
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106 void derivatives_2d (vector &dx,vector &dy,double &jac,
00107 vector &x,vector &y,double xi,double eta)
00108 {
00109 long i;
00110 double d1,d2,d3,d4,ddx,ddy,*nx,*ny;
00111
00112 nx = new double [x.n];
00113 ny = new double [x.n];
00114
00115 switch (x.n){
00116 case 3:{
00117 dx_bf_lin_3_2d (nx);
00118 dy_bf_lin_3_2d (ny);
00119 break;
00120 }
00121 case 4:{
00122 dx_bf_lin_4_2d (nx,eta);
00123 dy_bf_lin_4_2d (ny,xi);
00124 break;
00125 }
00126 case 6:{
00127 dx_bf_quad_3_2d (nx,xi,eta);
00128 dy_bf_quad_3_2d (ny,xi,eta);
00129 break;
00130 }
00131 case 8:{
00132 dx_bf_quad_4_2d (nx,xi,eta);
00133 dy_bf_quad_4_2d (ny,xi,eta);
00134 break;
00135 }
00136 default:{
00137 fprintf (stderr,"\n\n wrong number of nodes on 2D element in the function");
00138 fprintf (stderr,"\n derivatives_2d (%s, line %d).\n",__FILE__,__LINE__);
00139 }
00140 }
00141
00142 d1=0.0; d2=0.0; d3=0.0; d4=0.0;
00143 for (i=0;i<x.n;i++){
00144 d1+=nx[i]*x[i];
00145 d2+=ny[i]*x[i];
00146 d3+=nx[i]*y[i];
00147 d4+=ny[i]*y[i];
00148 }
00149
00150 jac = d1*d4 - d2*d3;
00151
00152 for (i=0;i<dx.n;i++){
00153 ddx=dx[i]; ddy=dy[i];
00154 dx[i]=(ddx*d4-ddy*d3)/jac;
00155 dy[i]=(ddy*d1-ddx*d2)/jac;
00156 }
00157 delete [] nx; delete [] ny;
00158 }
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171 void jac_2d (double &jac,vector &x,vector &y,double xi,double eta)
00172 {
00173 long i;
00174 double d1,d2,d3,d4,*nx,*ny;
00175
00176 nx = new double [x.n];
00177 ny = new double [x.n];
00178
00179 switch (x.n){
00180 case 3:{
00181 dx_bf_lin_3_2d (nx);
00182 dy_bf_lin_3_2d (ny);
00183 break;
00184 }
00185 case 4:{
00186 dx_bf_lin_4_2d (nx,eta);
00187 dy_bf_lin_4_2d (ny,xi);
00188 break;
00189 }
00190 case 6:{
00191 dx_bf_quad_3_2d (nx,xi,eta);
00192 dy_bf_quad_3_2d (ny,xi,eta);
00193 break;
00194 }
00195 case 8:{
00196 dx_bf_quad_4_2d (nx,xi,eta);
00197 dy_bf_quad_4_2d (ny,xi,eta);
00198 break;
00199 }
00200 default:{
00201 fprintf (stderr,"\n\n wrong number of nodes on 2D element in the function");
00202 fprintf (stderr,"\n jac_2d (%s, line %d).\n",__FILE__,__LINE__);
00203 }
00204 }
00205
00206 d1=0.0; d2=0.0; d3=0.0; d4=0.0;
00207 for (i=0;i<x.n;i++){
00208 d1+=nx[i]*x[i];
00209 d2+=ny[i]*x[i];
00210 d3+=nx[i]*y[i];
00211 d4+=ny[i]*y[i];
00212 }
00213
00214 jac = d1*d4 - d2*d3;
00215
00216 delete [] nx; delete [] ny;
00217 }
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232 void derivatives_3d (vector &dx,vector &dy,vector &dz,double &jac,
00233 vector &x,vector &y,vector &z,
00234 double xi,double eta,double zeta)
00235 {
00236 long i;
00237 double d1,d2,d3,d4,d5,d6,d7,d8,d9,ddx,ddy,ddz,*nx,*ny,*nz;
00238
00239 nx = new double [x.n];
00240 ny = new double [x.n];
00241 nz = new double [x.n];
00242
00243 switch (x.n){
00244 case 4:{
00245 dx_bf_lin_tet (nx);
00246 dy_bf_lin_tet (ny);
00247 dz_bf_lin_tet (nz);
00248 break;
00249 }
00250 case 6:{
00251 dx_bf_lin_wed_3d (nx,zeta);
00252 dy_bf_lin_wed_3d (ny,zeta);
00253 dz_bf_lin_wed_3d (nz,xi,eta);
00254 break;
00255 }
00256 case 8:{
00257 dx_bf_lin_hex_3d (nx,eta,zeta);
00258 dy_bf_lin_hex_3d (ny,xi,zeta);
00259 dz_bf_lin_hex_3d (nz,xi,eta);
00260 break;
00261 }
00262 case 10:{
00263 dx_bf_quad_tet(nx,xi,eta,zeta);
00264 dy_bf_quad_tet(ny,xi,eta,zeta);
00265 dz_bf_quad_tet(nz,xi,eta,zeta);
00266 break;
00267 }
00268 case 15:{
00269 dx_bf_quad_wed_3d (nx,xi,eta,zeta);
00270 dy_bf_quad_wed_3d (ny,xi,eta,zeta);
00271 dz_bf_quad_wed_3d (nz,xi,eta,zeta);
00272 break;
00273 }
00274 case 20:{
00275 dx_bf_quad_hex_3d (nx,xi,eta,zeta);
00276 dy_bf_quad_hex_3d (ny,xi,eta,zeta);
00277 dz_bf_quad_hex_3d (nz,xi,eta,zeta);
00278 break;
00279 }
00280 default:{
00281 fprintf (stderr,"\n\n wrong number of nodes on 3D element in the function");
00282 fprintf (stderr,"\n derivatives_3d (%s, line %d).\n",__FILE__,__LINE__);
00283 }
00284 }
00285
00286 d1=0.0; d2=0.0; d3=0.0; d4=0.0;
00287 d5=0.0; d6=0.0; d7=0.0; d8=0.0; d9=0.0;
00288 for (i=0;i<x.n;i++){
00289 d1+=nx[i]*x[i];
00290 d2+=ny[i]*x[i];
00291 d3+=nz[i]*x[i];
00292 d4+=nx[i]*y[i];
00293 d5+=ny[i]*y[i];
00294 d6+=nz[i]*y[i];
00295 d7+=nx[i]*z[i];
00296 d8+=ny[i]*z[i];
00297 d9+=nz[i]*z[i];
00298 }
00299
00300 jac = d1*d5*d9 + d4*d8*d3 + d7*d2*d6 - d7*d5*d3 - d8*d6*d1 - d9*d4*d2;
00301
00302 for (i=0;i<dx.n;i++){
00303 ddx=dx[i]; ddy=dy[i]; ddz=dz[i];
00304 dx[i]=(ddx*(d5*d9-d8*d6)+ddy*(d7*d6-d9*d4)+ddz*(d4*d8-d7*d5))/jac;
00305 dy[i]=(ddx*(d8*d3-d9*d2)+ddy*(d1*d9-d7*d3)+ddz*(d7*d2-d8*d1))/jac;
00306 dz[i]=(ddx*(d2*d6-d5*d3)+ddy*(d4*d3-d6*d1)+ddz*(d1*d5-d4*d2))/jac;
00307 }
00308
00309 delete [] nx; delete [] ny; delete [] nz;
00310 }
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323 void jac_3d (double &jac,vector &x,vector &y,vector &z,
00324 double xi,double eta,double zeta)
00325 {
00326 long i;
00327 double d1,d2,d3,d4,d5,d6,d7,d8,d9,*nx,*ny,*nz;
00328
00329 nx = new double [x.n];
00330 ny = new double [x.n];
00331 nz = new double [x.n];
00332
00333 switch (x.n){
00334 case 4:{
00335 dx_bf_lin_tet (nx);
00336 dy_bf_lin_tet (ny);
00337 dz_bf_lin_tet (nz);
00338 break;
00339 }
00340 case 6:{
00341 dx_bf_lin_wed_3d (nx,zeta);
00342 dy_bf_lin_wed_3d (ny,zeta);
00343 dz_bf_lin_wed_3d (nz,xi,eta);
00344 break;
00345 }
00346 case 8:{
00347 dx_bf_lin_hex_3d (nx,eta,zeta);
00348 dy_bf_lin_hex_3d (ny,xi,zeta);
00349 dz_bf_lin_hex_3d (nz,xi,eta);
00350 break;
00351 }
00352 case 10:{
00353 dx_bf_quad_tet (nx,xi,eta,zeta);
00354 dy_bf_quad_tet (ny,xi,eta,zeta);
00355 dz_bf_quad_tet (nz,xi,eta,zeta);
00356 break;
00357 }
00358 case 15:{
00359 dx_bf_quad_wed_3d (nx,xi,eta,zeta);
00360 dy_bf_quad_wed_3d (ny,xi,eta,zeta);
00361 dz_bf_quad_wed_3d (nz,xi,eta,zeta);
00362 break;
00363 }
00364 case 20:{
00365 dx_bf_quad_hex_3d (nx,xi,eta,zeta);
00366 dy_bf_quad_hex_3d (ny,xi,eta,zeta);
00367 dz_bf_quad_hex_3d (nz,xi,eta,zeta);
00368 break;
00369 }
00370 default:{
00371 fprintf (stderr,"\n\n wrong number of nodes on 3D element in the function");
00372 fprintf (stderr,"\n jac_3d (%s, line %d).\n",__FILE__,__LINE__);
00373 }
00374 }
00375
00376 d1=0.0; d2=0.0; d3=0.0; d4=0.0;
00377 d5=0.0; d6=0.0; d7=0.0; d8=0.0; d9=0.0;
00378 for (i=0;i<x.n;i++){
00379 d1+=nx[i]*x[i];
00380 d2+=ny[i]*x[i];
00381 d3+=nz[i]*x[i];
00382 d4+=nx[i]*y[i];
00383 d5+=ny[i]*y[i];
00384 d6+=nz[i]*y[i];
00385 d7+=nx[i]*z[i];
00386 d8+=ny[i]*z[i];
00387 d9+=nz[i]*z[i];
00388 }
00389
00390 jac = d1*d5*d9 + d4*d8*d3 + d7*d2*d6 - d7*d5*d3 - d8*d6*d1 - d9*d4*d2;
00391
00392 delete [] nx; delete [] ny; delete [] nz;
00393 }
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406 void jac1d2d (double &jac,vector &x,vector &y,double xi)
00407 {
00408 long i,n;
00409 double dx,dy,*nx;
00410
00411 n=x.n;
00412 nx = new double [n];
00413
00414 switch (n){
00415 case 2:{
00416 dx_bf_lin_1d (nx);
00417 break;
00418 }
00419 case 3:{
00420 dx_bf_quad_1d (nx,xi);
00421 break;
00422 }
00423 default:{
00424 fprintf (stderr,"\n\n wrong number of nodes on element in function");
00425 fprintf (stderr,"\n jac1d (%s, line %d).\n",__FILE__,__LINE__);
00426 }
00427 }
00428
00429 dx=0.0; dy=0.0;
00430 for (i=0;i<n;i++){
00431 dx+=nx[i]*x[i];
00432 dy+=nx[i]*y[i];
00433 }
00434
00435 jac = sqrt (dx*dx+dy*dy);
00436
00437 delete [] nx;
00438 }
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449 void jac1d3d (double &jac,vector &x,vector &y,vector &z,double xi)
00450 {
00451 long i,n;
00452 double dx,dy,dz,*nx;
00453
00454 n=x.n;
00455 nx = new double [n];
00456
00457 switch (n){
00458 case 2:{
00459 dx_bf_lin_1d (nx);
00460 break;
00461 }
00462 case 3:{
00463 dx_bf_quad_1d (nx,xi);
00464 break;
00465 }
00466 default:{
00467 fprintf (stderr,"\n\n wrong number of nodes on element in function");
00468 fprintf (stderr,"\n jac1d (%s, line %d).\n",__FILE__,__LINE__);
00469 }
00470 }
00471
00472 dx=0.0; dy=0.0; dz=0.0;
00473 for (i=0;i<n;i++){
00474 dx+=nx[i]*x[i];
00475 dy+=nx[i]*y[i];
00476 dz+=nx[i]*z[i];
00477 }
00478
00479 jac = sqrt (dx*dx+dy*dy+dz*dz);
00480
00481 delete [] nx;
00482 }
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493 void jac2d3d (double &jac,vector &x,vector &y,vector &z,double xi,double eta)
00494 {
00495 long i,n;
00496 double dxu,dxv,dyu,dyv,dzu,dzv,g11,g22,g12;
00497 double *nu,*nv;
00498
00499 n=x.n;
00500 nu = new double [n]; nv = new double [n];
00501
00502 switch (n){
00503 case 3:{
00504 dx_bf_lin_3_2d (nu);
00505 dy_bf_lin_3_2d (nv);
00506 break;
00507 }
00508 case 4:{
00509 dx_bf_lin_4_2d (nu,eta);
00510 dy_bf_lin_4_2d (nv,xi);
00511 break;
00512 }
00513 case 6:{
00514 dx_bf_quad_3_2d (nu,xi,eta);
00515 dy_bf_quad_3_2d (nv,xi,eta);
00516 break;
00517 }
00518 case 8:{
00519 dx_bf_quad_4_2d (nu,xi,eta);
00520 dy_bf_quad_4_2d (nv,xi,eta);
00521 break;
00522 }
00523 default:{
00524 fprintf (stderr,"\n\n wrong number of nodes on 2D element in the function");
00525 fprintf (stderr,"\n jac2d3d (%s, line %d).\n",__FILE__,__LINE__);
00526 }
00527 }
00528
00529 dxu=0.0; dyu=0.0; dzu=0.0; dxv=0.0; dyv=0.0; dzv=0.0;
00530 for (i=0;i<n;i++){
00531 dxu+=nu[i]*x[i];
00532 dyu+=nu[i]*y[i];
00533 dzu+=nu[i]*z[i];
00534 dxv+=nv[i]*x[i];
00535 dyv+=nv[i]*y[i];
00536 dzv+=nv[i]*z[i];
00537 }
00538
00539 g11=dxu*dxu+dyu*dyu+dzu*dzu;
00540 g22=dxv*dxv+dyv*dyv+dzv*dzv;
00541 g12=dxu*dxv+dyu*dyv+dzu*dzv;
00542
00543 jac = sqrt (g11*g22-g12*g12);
00544
00545 delete [] nu; delete [] nv;
00546 }
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559 void jac1d_2d (double &jac,vector &x,vector &y,double xi,long edid)
00560 {
00561 vector lx,ly;
00562
00563 if (x.n==3){
00564 allocv (2,lx); allocv (2,ly);
00565 if (edid==0){
00566 lx[0]=x[0]; lx[1]=x[1];
00567 ly[0]=y[0]; ly[1]=y[1];
00568 }
00569 if (edid==1){
00570 lx[0]=x[1]; lx[1]=x[2];
00571 ly[0]=y[1]; ly[1]=y[2];
00572 }
00573 if (edid==2){
00574 lx[0]=x[2]; lx[1]=x[0];
00575 ly[0]=y[2]; ly[1]=y[0];
00576 }
00577 }
00578
00579 if (x.n==4){
00580 allocv (2,lx); allocv (2,ly);
00581 if (edid==0){
00582 lx[0]=x[0]; lx[1]=x[1];
00583 ly[0]=y[0]; ly[1]=y[1];
00584 }
00585 if (edid==1){
00586 lx[0]=x[1]; lx[1]=x[2];
00587 ly[0]=y[1]; ly[1]=y[2];
00588 }
00589 if (edid==2){
00590 lx[0]=x[2]; lx[1]=x[3];
00591 ly[0]=y[2]; ly[1]=y[3];
00592 }
00593 if (edid==3){
00594 lx[0]=x[3]; lx[1]=x[0];
00595 ly[0]=y[3]; ly[1]=y[0];
00596 }
00597 }
00598
00599 if (x.n==6){
00600 allocv (3,lx); allocv (3,ly);
00601 if (edid==0){
00602 lx[0]=x[0]; lx[1]=x[1]; lx[2]=x[3];
00603 ly[0]=y[0]; ly[1]=y[1]; ly[2]=y[3];
00604 }
00605 if (edid==1){
00606 lx[0]=x[1]; lx[1]=x[2]; lx[2]=x[4];
00607 ly[0]=y[1]; ly[1]=y[2]; ly[2]=y[4];
00608 }
00609 if (edid==2){
00610 lx[0]=x[2]; lx[1]=x[0]; lx[2]=x[5];
00611 ly[0]=y[2]; ly[1]=y[0]; ly[2]=y[5];
00612 }
00613 }
00614
00615 if (x.n==8){
00616 allocv (3,lx); allocv (3,ly);
00617 if (edid==0){
00618 lx[0]=x[0]; lx[1]=x[1]; lx[2]=x[4];
00619 ly[0]=y[0]; ly[1]=y[1]; ly[2]=y[4];
00620 }
00621 if (edid==1){
00622 lx[0]=x[1]; lx[1]=x[2]; lx[2]=x[5];
00623 ly[0]=y[1]; ly[1]=y[2]; ly[2]=y[5];
00624 }
00625 if (edid==2){
00626 lx[0]=x[2]; lx[1]=x[3]; lx[2]=x[6];
00627 ly[0]=y[2]; ly[1]=y[3]; ly[2]=y[6];
00628 }
00629 if (edid==3){
00630 lx[0]=x[3]; lx[1]=x[0]; lx[2]=x[7];
00631 ly[0]=y[3]; ly[1]=y[0]; ly[2]=y[7];
00632 }
00633 }
00634
00635 jac1d2d (jac,lx,ly,xi);
00636
00637 destrv (lx); destrv (ly);
00638 }
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650 void jac2d_3d (double &jac,vector &x,vector &y,vector &z,double xi,double eta,long sid)
00651 {
00652 vector lx,ly,lz;
00653 ivector sn;
00654 long i;
00655
00656 if (x.n==4){
00657
00658 allocv (3,lx); allocv (3,ly); allocv (3,lz);
00659 allocv (3,sn);
00660 if (sid==0){
00661
00662
00663
00664
00665
00666 lx[0]=x[2]; lx[1]=x[1]; lx[2]=x[3];
00667 ly[0]=y[2]; ly[1]=y[1]; ly[2]=y[3];
00668 lz[0]=z[2]; lz[1]=z[1]; lz[2]=z[3];
00669 }
00670 if (sid==1){
00671 lx[0]=x[0]; lx[1]=x[2]; lx[2]=x[3];
00672 ly[0]=y[0]; ly[1]=y[2]; ly[2]=y[3];
00673 lz[0]=z[0]; lz[1]=z[2]; lz[2]=z[3];
00674 }
00675 if (sid==2){
00676 lx[0]=x[1]; lx[1]=x[0]; lx[2]=x[3];
00677 ly[0]=y[1]; ly[1]=y[0]; ly[2]=y[3];
00678 lz[0]=z[1]; lz[1]=z[0]; lz[2]=z[3];
00679 }
00680 if (sid==3){
00681 lx[0]=x[0]; lx[1]=x[1]; lx[2]=x[2];
00682 ly[0]=y[0]; ly[1]=y[1]; ly[2]=y[2];
00683 lz[0]=z[0]; lz[1]=z[1]; lz[2]=z[2];
00684 }
00685 }
00686
00687 if (x.n==10)
00688 {
00689
00690 allocv (6,lx); allocv (6,ly); allocv (6,lz);
00691 allocv (6,sn);
00692 quadtetrahedral_surfnod (sn.a,sid);
00693 for (i=0; i<6; i++)
00694 {
00695 lx[i] = x[sn[i]];
00696 ly[i] = y[sn[i]];
00697 lz[i] = z[sn[i]];
00698 }
00699 }
00700
00701
00702 if (x.n==8){
00703
00704 allocv (4,lx); allocv (4,ly); allocv (4,lz);
00705 allocv (4,sn);
00706 if (sid==0){
00707 lx[0]=x[0]; lx[1]=x[3]; lx[2]=x[7]; lx[3]=x[4];
00708 ly[0]=y[0]; ly[1]=y[3]; ly[2]=y[7]; ly[3]=y[4];
00709 lz[0]=z[0]; lz[1]=z[3]; lz[2]=z[7]; lz[3]=z[4];
00710 }
00711 if (sid==1){
00712 lx[0]=x[1]; lx[1]=x[0]; lx[2]=x[4]; lx[3]=x[5];
00713 ly[0]=y[1]; ly[1]=y[0]; ly[2]=y[4]; ly[3]=y[5];
00714 lz[0]=z[1]; lz[1]=z[0]; lz[2]=z[4]; lz[3]=z[5];
00715 }
00716 if (sid==2){
00717 lx[0]=x[2]; lx[1]=x[1]; lx[2]=x[5]; lx[3]=x[6];
00718 ly[0]=y[2]; ly[1]=y[1]; ly[2]=y[5]; ly[3]=y[6];
00719 lz[0]=z[2]; lz[1]=z[1]; lz[2]=z[5]; lz[3]=z[6];
00720 }
00721 if (sid==3){
00722 lx[0]=x[3]; lx[1]=x[2]; lx[2]=x[6]; lx[3]=x[7];
00723 ly[0]=y[3]; ly[1]=y[2]; ly[2]=y[6]; ly[3]=y[7];
00724 lz[0]=z[3]; lz[1]=z[2]; lz[2]=z[6]; lz[3]=z[7];
00725 }
00726 if (sid==4){
00727 lx[0]=x[0]; lx[1]=x[1]; lx[2]=x[2]; lx[3]=x[3];
00728 ly[0]=y[0]; ly[1]=y[1]; ly[2]=y[2]; ly[3]=y[3];
00729 lz[0]=z[0]; lz[1]=z[1]; lz[2]=z[2]; lz[3]=z[3];
00730 }
00731 if (sid==5){
00732 lx[0]=x[4]; lx[1]=x[7]; lx[2]=x[6]; lx[3]=x[5];
00733 ly[0]=y[4]; ly[1]=y[7]; ly[2]=y[6]; ly[3]=y[5];
00734 lz[0]=z[4]; lz[1]=z[7]; lz[2]=z[6]; lz[3]=z[5];
00735 }
00736 }
00737
00738 if (x.n==20){
00739
00740 allocv (8,lx); allocv (8,ly); allocv (8,lz);
00741 allocv (8,sn);
00742 if (sid==0){
00743 lx[0]=x[0]; lx[1]=x[3]; lx[2]=x[7]; lx[3]=x[4]; lx[4]=x[11]; lx[5]=x[15]; lx[6]=x[19]; lx[7]=x[12];
00744 ly[0]=y[0]; ly[1]=y[3]; ly[2]=y[7]; ly[3]=y[4]; ly[4]=y[11]; ly[5]=y[15]; ly[6]=y[19]; ly[7]=y[12];
00745 lz[0]=z[0]; lz[1]=z[3]; lz[2]=z[7]; lz[3]=z[4]; lz[4]=z[11]; lz[5]=z[15]; lz[6]=z[19]; lz[7]=z[12];
00746 }
00747 if (sid==1){
00748 lx[0]=x[1]; lx[1]=x[0]; lx[2]=x[4]; lx[3]=x[5]; lx[4]=x[8]; lx[5]=x[12]; lx[6]=x[16]; lx[7]=x[13];
00749 ly[0]=y[1]; ly[1]=y[0]; ly[2]=y[4]; ly[3]=y[5]; ly[4]=y[8]; ly[5]=y[12]; ly[6]=y[16]; ly[7]=y[13];
00750 lz[0]=z[1]; lz[1]=z[0]; lz[2]=z[4]; lz[3]=z[5]; lz[4]=z[8]; lz[5]=z[12]; lz[6]=z[16]; lz[7]=z[13];
00751 }
00752 if (sid==2){
00753 lx[0]=x[2]; lx[1]=x[1]; lx[2]=x[5]; lx[3]=x[6]; lx[4]=x[9]; lx[5]=x[13]; lx[6]=x[17]; lx[7]=x[14];
00754 ly[0]=y[2]; ly[1]=y[1]; ly[2]=y[5]; ly[3]=y[6]; ly[4]=y[9]; ly[5]=y[13]; ly[6]=y[17]; ly[7]=y[14];
00755 lz[0]=z[2]; lz[1]=z[1]; lz[2]=z[5]; lz[3]=z[6]; lz[4]=z[9]; lz[5]=z[13]; lz[6]=z[17]; lz[7]=z[14];
00756 }
00757 if (sid==3){
00758 lx[0]=x[3]; lx[1]=x[2]; lx[2]=x[6]; lx[3]=x[7]; lx[4]=x[10]; lx[5]=x[14]; lx[6]=x[18]; lx[7]=x[15];
00759 ly[0]=y[3]; ly[1]=y[2]; ly[2]=y[6]; ly[3]=y[7]; ly[4]=y[10]; ly[5]=y[14]; ly[6]=y[18]; ly[7]=y[15];
00760 lz[0]=z[3]; lz[1]=z[2]; lz[2]=z[6]; lz[3]=z[7]; lz[4]=z[10]; lz[5]=z[14]; lz[6]=z[18]; lz[7]=z[15];
00761 }
00762 if (sid==4){
00763 lx[0]=x[0]; lx[1]=x[1]; lx[2]=x[2]; lx[3]=x[3]; lx[4]=x[8]; lx[5]=x[9]; lx[6]=x[10]; lx[7]=x[11];
00764 ly[0]=y[0]; ly[1]=y[1]; ly[2]=y[2]; ly[3]=y[3]; ly[4]=y[8]; ly[5]=y[9]; ly[6]=y[10]; ly[7]=y[11];
00765 lz[0]=z[0]; lz[1]=z[1]; lz[2]=z[2]; lz[3]=z[3]; lz[4]=z[8]; lz[5]=z[9]; lz[6]=z[10]; lz[7]=z[11];
00766 }
00767 if (sid==5){
00768 lx[0]=x[4]; lx[1]=x[7]; lx[2]=x[6]; lx[3]=x[5]; lx[4]=x[19]; lx[5]=x[18]; lx[6]=x[17]; lx[7]=x[16];
00769 ly[0]=y[4]; ly[1]=y[7]; ly[2]=y[6]; ly[3]=y[5]; ly[4]=y[19]; ly[5]=y[18]; ly[6]=y[17]; ly[7]=y[16];
00770 lz[0]=z[4]; lz[1]=z[7]; lz[2]=z[6]; lz[3]=z[5]; lz[4]=z[19]; lz[5]=z[18]; lz[6]=z[17]; lz[7]=z[16];
00771 }
00772 }
00773
00774
00775 jac2d3d (jac,lx,ly,lz,xi,eta);
00776
00777 destrv (lx); destrv (ly); destrv (lz);
00778 destrv (sn);
00779 }
00780