Changeset 8653


Ignore:
Timestamp:
06/16/11 16:08:22 (14 years ago)
Author:
Mathieu Morlighem
Message:

moved some l1l2l3 to basis and dh1dh3 to dbasis

Location:
issm/trunk/src/c/objects/Elements
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r8649 r8653  
    17731773        double    surface_normal[3];
    17741774        double    Jdet2d,DL_scalar;
    1775         double    l1l6[NUMVERTICES];
     1775        double    basis[NUMVERTICES];
    17761776        GaussPenta *gauss=NULL;
    17771777        double Ke_g[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    … …  
    17921792
    17931793                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    1794                 GetNodalFunctionsP1(&l1l6[0], gauss);
     1794                GetNodalFunctionsP1(&basis[0], gauss);
    17951795
    17961796                /**********************Do not forget the sign**********************************/
    … …  
    17981798                /******************************************************************************/
    17991799
    1800                 TripleMultiply( l1l6,1,numdof,1,
     1800                TripleMultiply( basis,1,numdof,1,
    18011801                                        &DL_scalar,1,1,0,
    1802                                         l1l6,1,numdof,0,
     1802                                        basis,1,numdof,0,
    18031803                                        &Ke_g[0][0],0);
    18041804
    … …  
    18481848        double     Bprime_advec[3][numdof];
    18491849        double     L[numdof];
    1850         double     dh1dh6[3][6];
     1850        double     dbasis[3][6];
    18511851        double     D_scalar_conduct,D_scalar_advec;
    18521852        double     D_scalar_trans,D_scalar_artdiff;
    … …  
    19701970                else if(artdiff==2){
    19711971
    1972                         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
     1972                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    19731973
    19741974                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
    … …  
    19771977                                for(j=0;j<numdof;j++){
    19781978                                        Ke_gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*
    1979                                           ((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i])*((u-um)*dh1dh6[0][j]+(v-vm)*dh1dh6[1][j]+(w-wm)*dh1dh6[2][j]);
     1979                                          ((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i])*((u-um)*dbasis[0][j]+(v-vm)*dbasis[1][j]+(w-wm)*dbasis[2][j]);
    19801980                                }
    19811981                        }
    … …  
    19831983                                for(i=0;i<numdof;i++){
    19841984                                        for(j=0;j<numdof;j++){
    1985                                                 Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i]);
     1985                                                Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]);
    19861986                                        }
    19871987                                }
    … …  
    20132013        double    xyz_list[NUMVERTICES][3];
    20142014        double   xyz_list_tria[NUMVERTICES2D][3];
    2015         double    l1l6[NUMVERTICES];
     2015        double    basis[NUMVERTICES];
    20162016        double    D_scalar;
    20172017        double    Ke_gaussian[numdof][numdof]={0.0};
    … …  
    20392039               
    20402040                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    2041                 GetNodalFunctionsP1(&l1l6[0], gauss);
     2041                GetNodalFunctionsP1(&basis[0], gauss);
    20422042                               
    20432043                D_scalar=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(rho_ice*heatcapacity);
    20442044                if(dt) D_scalar=dt*D_scalar;
    20452045
    2046                 TripleMultiply(&l1l6[0],numdof,1,0,
     2046                TripleMultiply(&basis[0],numdof,1,0,
    20472047                                        &D_scalar,1,1,0,
    2048                                         &l1l6[0],1,numdof,0,
     2048                                        &basis[0],1,numdof,0,
    20492049                                        &Ke_gaussian[0][0],0);
    20502050
    … …  
    21382138        double     Bprime_advec[3][numdof];
    21392139        double     L[numdof];
    2140         double     dh1dh6[3][6];
     2140        double     dbasis[3][6];
    21412141        double     D_scalar_conduct,D_scalar_advec;
    21422142        double     D_scalar_trans,D_scalar_artdiff;
    … …  
    22532253                else if(artdiff==2){
    22542254
    2255                         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
     2255                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    22562256
    22572257                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
    … …  
    22602260                                for(j=0;j<numdof;j++){
    22612261                                        Ke_gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*
    2262                                           ((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i])*((u-um)*dh1dh6[0][j]+(v-vm)*dh1dh6[1][j]+(w-wm)*dh1dh6[2][j]);
     2262                                          ((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i])*((u-um)*dbasis[0][j]+(v-vm)*dbasis[1][j]+(w-wm)*dbasis[2][j]);
    22632263                                }
    22642264                        }
    … …  
    22662266                                for(i=0;i<numdof;i++){
    22672267                                        for(j=0;j<numdof;j++){
    2268                                                 Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i]);
     2268                                                Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]);
    22692269                                        }
    22702270                                }
    … …  
    22972297        double    xyz_list[NUMVERTICES][3];
    22982298        double   xyz_list_tria[NUMVERTICES2D][3];
    2299         double    l1l6[NUMVERTICES];
     2299        double    basis[NUMVERTICES];
    23002300        double    D_scalar;
    23012301        double    Ke_gaussian[numdof][numdof]={0.0};
    … …  
    23232323               
    23242324                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    2325                 GetNodalFunctionsP1(&l1l6[0], gauss);
     2325                GetNodalFunctionsP1(&basis[0], gauss);
    23262326                               
    23272327                D_scalar=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice);
    23282328                if(dt) D_scalar=dt*D_scalar;
    23292329
    2330                 TripleMultiply(&l1l6[0],numdof,1,0,
     2330                TripleMultiply(&basis[0],numdof,1,0,
    23312331                                        &D_scalar,1,1,0,
    2332                                         &l1l6[0],1,numdof,0,
     2332                                        &basis[0],1,numdof,0,
    23332333                                        &Ke_gaussian[0][0],0);
    23342334
    … …  
    25222522        double      dw[3];
    25232523        double      xyz_list[NUMVERTICES][3];
    2524         double      l1l6[6]; //for the six nodes of the penta
    2525         double      dh1dh6[3][6]; //for the six nodes of the penta
     2524        double      basis[6]; //for the six nodes of the penta
     2525        double      dbasis[3][6]; //for the six nodes of the penta
    25262526        GaussPenta *gauss=NULL;
    25272527
    … …  
    25462546
    25472547                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    2548                 GetNodalFunctionsP1(&l1l6[0], gauss);
    2549                 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
     2548                GetNodalFunctionsP1(&basis[0], gauss);
     2549                GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    25502550
    25512551                vzmacayeal_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    … …  
    25552555
    25562556                for(i=0;i<NUMVERTICES;i++){
    2557                         pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i];
    2558                         pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i];
    2559                         pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dh1dh6[0][i]+dw[1]*dh1dh6[1][i]+2*dw[2]*dh1dh6[2][i]);
    2560                         pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i];
     2557                        pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
     2558                        pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
     2559                        pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
     2560                        pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*basis[i];
    25612561                }
    25622562        }
    … …  
    25842584        double     xyz_list_tria[NUMVERTICES2D][3];
    25852585        double      xyz_list[NUMVERTICES][3];
    2586         double      l1l6[6]; //for the six nodes of the penta
     2586        double      basis[6]; //for the six nodes of the penta
    25872587        Tria*       tria=NULL;
    25882588        Friction*   friction=NULL;
    … …  
    26162616
    26172617                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    2618                 GetNodalFunctionsP1(l1l6, gauss);
     2618                GetNodalFunctionsP1(basis, gauss);
    26192619
    26202620                vzmacayeal_input->GetParameterValue(&w, gauss);
    … …  
    26272627
    26282628                for(i=0;i<NUMVERTICES2D;i++){
    2629                         pe->values[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*l1l6[i];
    2630                         pe->values[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*l1l6[i];
    2631                         pe->values[i*NDOF4+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*l1l6[i];
     2629                        pe->values[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
     2630                        pe->values[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
     2631                        pe->values[i*NDOF4+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
    26322632                }
    26332633        }
    … …  
    26672667        double      dw[3];
    26682668        double      xyz_list[NUMVERTICES][3];
    2669         double      l1l6[6]; //for the six nodes of the penta
    2670         double      dh1dh6[3][6]; //for the six nodes of the penta
     2669        double      basis[6]; //for the six nodes of the penta
     2670        double      dbasis[3][6]; //for the six nodes of the penta
    26712671        GaussPenta *gauss=NULL;
    26722672
    … …  
    26912691
    26922692                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    2693                 GetNodalFunctionsP1(&l1l6[0], gauss);
    2694                 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
     2693                GetNodalFunctionsP1(&basis[0], gauss);
     2694                GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    26952695
    26962696                vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    … …  
    27002700
    27012701                for(i=0;i<NUMVERTICES;i++){
    2702                         pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i];
    2703                         pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i];
    2704                         pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dh1dh6[0][i]+dw[1]*dh1dh6[1][i]+2*dw[2]*dh1dh6[2][i]);
    2705                         pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i];
     2702                        pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
     2703                        pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
     2704                        pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
     2705                        pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*basis[i];
    27062706                }
    27072707        }
    … …  
    27292729        double     xyz_list_tria[NUMVERTICES2D][3];
    27302730        double      xyz_list[NUMVERTICES][3];
    2731         double      l1l6[6]; //for the six nodes of the penta
     2731        double      basis[6]; //for the six nodes of the penta
    27322732        Tria*       tria=NULL;
    27332733        Friction*   friction=NULL;
    … …  
    27612761
    27622762                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    2763                 GetNodalFunctionsP1(l1l6, gauss);
     2763                GetNodalFunctionsP1(basis, gauss);
    27642764
    27652765                vzpattyn_input->GetParameterValue(&w, gauss);
    … …  
    27722772
    27732773                for(i=0;i<NUMVERTICES2D;i++){
    2774                         pe->values[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*l1l6[i];
    2775                         pe->values[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*l1l6[i];
    2776                         pe->values[i*NDOF4+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*l1l6[i];
     2774                        pe->values[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
     2775                        pe->values[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
     2776                        pe->values[i*NDOF4+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
    27772777                }
    27782778        }
    … …  
    29722972        double      driving_stress_baseline,thickness;
    29732973        double      xyz_list[NUMVERTICES][3];
    2974         double      l1l6[6];
     2974        double      basis[6];
    29752975        GaussPenta  *gauss=NULL;
    29762976
    … …  
    29902990
    29912991                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    2992                 GetNodalFunctionsP1(l1l6, gauss);
     2992                GetNodalFunctionsP1(basis, gauss);
    29932993
    29942994                thickness_input->GetParameterValue(&thickness, gauss);
    … …  
    29972997                driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG();
    29982998
    2999                 for(i=0;i<NUMVERTICES;i++) for(j=0;j<NDOF2;j++) pe->values[i*NDOF2+j]+= -driving_stress_baseline*slope[j]*Jdet*gauss->weight*l1l6[i];
     2999                for(i=0;i<NUMVERTICES;i++) for(j=0;j<NDOF2;j++) pe->values[i*NDOF2+j]+= -driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i];
    30003000        }
    30013001
    … …  
    31103110        double          bed_normal[3];
    31113111        double      dz[3];
    3112         double      l1l6[6]; //for the six nodes of the penta
     3112        double      basis[6]; //for the six nodes of the penta
    31133113        double      Jdet2d;
    31143114        GaussPenta  *gauss=NULL;
    … …  
    31393139
    31403140                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    3141                 GetNodalFunctionsP1(l1l6, gauss);
     3141                GetNodalFunctionsP1(basis, gauss);
    31423142
    31433143                BedNormal(&bed_normal[0],xyz_list_tria);
    … …  
    31553155                water_pressure=gravity*rho_water*bed;
    31563156
    3157                 for(i=0;i<NUMVERTICES;i++) for(j=0;j<3;j++) pe->values[i*NDOF4+j]+=(water_pressure+damper)*gauss->weight*Jdet2d*l1l6[i]*bed_normal[j];
     3157                for(i=0;i<NUMVERTICES;i++) for(j=0;j<3;j++) pe->values[i*NDOF4+j]+=(water_pressure+damper)*gauss->weight*Jdet2d*basis[i]*bed_normal[j];
    31583158        }
    31593159
    … …  
    32043204        double     dudx,dvdy,dwdz;
    32053205        double     du[3],dv[3],dw[3];
    3206         double     l1l6[6];
     3206        double     basis[6];
    32073207        GaussPenta *gauss=NULL;
    32083208
    … …  
    32273227
    32283228                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3229                 GetNodalFunctionsP1(l1l6, gauss);
     3229                GetNodalFunctionsP1(basis, gauss);
    32303230
    32313231                vx_input->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss);
    … …  
    32393239                dvdy=dv[1];
    32403240
    3241                 for (i=0;i<numdof;i++) pe->values[i] += (dudx+dvdy+dwdz)*Jdet*gauss->weight*l1l6[i];
     3241                for (i=0;i<numdof;i++) pe->values[i] += (dudx+dvdy+dwdz)*Jdet*gauss->weight*basis[i];
    32423242        }
    32433243
    … …  
    32623262        double     vx,vy,vz,dbdx,dbdy,basalmeltingvalue;
    32633263        double     slope[3];
    3264         double     l1l6[NUMVERTICES];
     3264        double     basis[NUMVERTICES];
    32653265        GaussPenta* gauss=NULL;
    32663266
    … …  
    33023302
    33033303                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    3304                 GetNodalFunctionsP1(&l1l6[0], gauss);
    3305 
    3306                 for(i=0;i<numdof;i++) pe->values[i]+=-Jdet2d*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*l1l6[i];
     3304                GetNodalFunctionsP1(&basis[0], gauss);
     3305
     3306                for(i=0;i<numdof;i++) pe->values[i]+=-Jdet2d*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*basis[i];
    33073307        }
    33083308
    … …  
    33473347        double xyz_list[NUMVERTICES][3];
    33483348        double L[numdof];
    3349         double dh1dh6[3][6];
     3349        double dbasis[3][6];
    33503350        double epsilon[6];
    33513351        GaussPenta *gauss=NULL;
    … …  
    33943394
    33953395                if(artdiff==2){
    3396                         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
     3396                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    33973397
    33983398                        vx_input->GetParameterValue(&u, gauss);
    … …  
    34023402                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
    34033403
    3404                         for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
     3404                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
    34053405                        if(dt){
    3406                                 for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_transient*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
     3406                                for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
    34073407                        }
    34083408                }
    … …  
    34283428        double     xyz_list[NUMVERTICES][3];
    34293429        double     xyz_list_tria[NUMVERTICES2D][3];
    3430         double     l1l6[NUMVERTICES];
     3430        double     basis[NUMVERTICES];
    34313431        GaussPenta* gauss=NULL;
    34323432
    … …  
    34553455
    34563456                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    3457                 GetNodalFunctionsP1(&l1l6[0], gauss);
     3457                GetNodalFunctionsP1(&basis[0], gauss);
    34583458
    34593459                pressure_input->GetParameterValue(&pressure,gauss);
    … …  
    34633463                if(dt) scalar_ocean=dt*scalar_ocean;
    34643464
    3465                 for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*l1l6[i];
     3465                for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*basis[i];
    34663466        }
    34673467
    … …  
    34863486        double     basalfriction,alpha2,vx,vy;
    34873487        double     scalar;
    3488         double     l1l6[NUMVERTICES];
     3488        double     basis[NUMVERTICES];
    34893489        Friction*  friction=NULL;
    34903490        GaussPenta* gauss=NULL;
    … …  
    35203520
    35213521                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    3522                 GetNodalFunctionsP1(&l1l6[0], gauss);
     3522                GetNodalFunctionsP1(&basis[0], gauss);
    35233523
    35243524                geothermalflux_input->GetParameterValue(&geothermalflux_value,gauss);
    … …  
    35313531                if(dt) scalar=dt*scalar;
    35323532
    3533                 for(i=0;i<numdof;i++) pe->values[i]+=scalar*l1l6[i];
     3533                for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    35343534        }
    35353535
    … …  
    36163616        double xyz_list[NUMVERTICES][3];
    36173617        double L[numdof];
    3618         double dh1dh6[3][6];
     3618        double dbasis[3][6];
    36193619        double epsilon[6];
    36203620        GaussPenta *gauss=NULL;
    … …  
    36633663
    36643664                if(artdiff==2){
    3665                         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
     3665                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    36663666
    36673667                        vx_input->GetParameterValue(&u, gauss);
    … …  
    36713671                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
    36723672
    3673                         for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
     3673                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
    36743674                        if(dt){
    3675                                 for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_transient*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
     3675                                for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
    36763676                        }
    36773677                }
    … …  
    36973697        double     xyz_list[NUMVERTICES][3];
    36983698        double     xyz_list_tria[NUMVERTICES2D][3];
    3699         double     l1l6[NUMVERTICES];
     3699        double     basis[NUMVERTICES];
    37003700        GaussPenta* gauss=NULL;
    37013701
    … …  
    37243724
    37253725                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    3726                 GetNodalFunctionsP1(&l1l6[0], gauss);
     3726                GetNodalFunctionsP1(&basis[0], gauss);
    37273727
    37283728                pressure_input->GetParameterValue(&pressure,gauss);
    … …  
    37323732                if(dt) scalar_ocean=dt*scalar_ocean;
    37333733
    3734                 for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*l1l6[i];
     3734                for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*basis[i];
    37353735        }
    37363736
    … …  
    37553755        double     basalfriction,alpha2,vx,vy;
    37563756        double     scalar;
    3757         double     l1l6[NUMVERTICES];
     3757        double     basis[NUMVERTICES];
    37583758        Friction*  friction=NULL;
    37593759        GaussPenta* gauss=NULL;
    … …  
    37893789
    37903790                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    3791                 GetNodalFunctionsP1(&l1l6[0], gauss);
     3791                GetNodalFunctionsP1(&basis[0], gauss);
    37923792
    37933793                geothermalflux_input->GetParameterValue(&geothermalflux_value,gauss);
    … …  
    38003800                if(dt) scalar=dt*scalar;
    38013801
    3802                 for(i=0;i<numdof;i++) pe->values[i]+=scalar*l1l6[i];
     3802                for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    38033803        }
    38043804

  • issm/trunk/src/c/objects/Elements/PentaRef.cpp

    r8417 r8653  
    7070         */
    7171
    72         double dh1dh6[3][NUMNODESP1];
    73 
    74         /*Get dh1dh6 in actual coordinate system: */
    75         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
     72        double dbasis[3][NUMNODESP1];
     73
     74        /*Get dbasis in actual coordinate system: */
     75        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
    7676
    7777        /*Build B: */
    7878        for (int i=0;i<NUMNODESP1;i++){
    79                 *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dh1dh6[0][i];
     79                *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i];
    8080                *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.0;
    8181
    8282                *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0.0;
    83                 *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dh1dh6[1][i];
    84 
    85                 *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dh1dh6[1][i];
    86                 *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];
    87 
    88         }
    89 
     83                *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
     84
     85                *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dbasis[1][i];
     86                *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][i];
     87        }
    9088}
    9189/*}}}*/
    … …  
    151149         */
    152150
    153         double dh1dh6[3][NUMNODESP1];
    154 
    155         /*Get dh1dh6 in actual coordinate system: */
    156         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
     151        double dbasis[3][NUMNODESP1];
     152
     153        /*Get dbasis in actual coordinate system: */
     154        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
    157155
    158156        /*Build B: */
    159157        for (int i=0;i<NUMNODESP1;i++){
    160                 *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dh1dh6[0][i];
     158                *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i];
    161159                *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.0;
    162160
    163161                *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0.0;
    164                 *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dh1dh6[1][i];
    165 
    166                 *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dh1dh6[1][i];
    167                 *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];
    168 
    169                 *(B+NDOF2*NUMNODESP1*3+NDOF2*i)=(float).5*dh1dh6[2][i];
     162                *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
     163
     164                *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dbasis[1][i];
     165                *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][i];
     166
     167                *(B+NDOF2*NUMNODESP1*3+NDOF2*i)=(float).5*dbasis[2][i];
    170168                *(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
    171169
    172170                *(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
    173                 *(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=(float).5*dh1dh6[2][i];
     171                *(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=(float).5*dbasis[2][i];
    174172        }
    175173
    … …  
    190188         * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
    191189         */
    192         double dh1dh6[3][NUMNODESP1];
    193 
    194         /*Get dh1dh6 in actual coordinate system: */
    195         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
     190        double dbasis[3][NUMNODESP1];
     191
     192        /*Get dbasis in actual coordinate system: */
     193        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss_coord);
    196194
    197195        /*Build BPrime: */
    198196        for (int i=0;i<NUMNODESP1;i++){
    199                 *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=2.0*dh1dh6[0][i];
    200                 *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=dh1dh6[1][i];
    201 
    202                 *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=dh1dh6[0][i];
    203                 *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2.0*dh1dh6[1][i];
    204 
    205                 *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=dh1dh6[1][i];
    206                 *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dh1dh6[0][i];
    207 
    208                 *(B+NDOF2*NUMNODESP1*3+NDOF2*i)=dh1dh6[2][i];
     197                *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=2.0*dbasis[0][i];
     198                *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=dbasis[1][i];
     199
     200                *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=dbasis[0][i];
     201                *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2.0*dbasis[1][i];
     202
     203                *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=dbasis[1][i];
     204                *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dbasis[0][i];
     205
     206                *(B+NDOF2*NUMNODESP1*3+NDOF2*i)=dbasis[2][i];
    209207                *(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
    210208
    211209                *(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
    212                 *(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=dh1dh6[2][i];
     210                *(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=dbasis[2][i];
    213211        }
    214212}
    … …  
    401399
    402400        /*Same thing in the actual coordinate system: */
    403         double dh1dh6[3][NUMNODESP1];
    404 
    405         /*Get dh1dh6 in actual coordinates system : */
    406         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
     401        double dbasis[3][NUMNODESP1];
     402
     403        /*Get dbasis in actual coordinates system : */
     404        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
    407405
    408406        /*Build B': */
    409407        for (int i=0;i<NUMNODESP1;i++){
    410                 *(B_artdiff+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i];
    411                 *(B_artdiff+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i];
     408                *(B_artdiff+NDOF1*NUMNODESP1*0+NDOF1*i)=dbasis[0][i];
     409                *(B_artdiff+NDOF1*NUMNODESP1*1+NDOF1*i)=dbasis[1][i];
    412410        }
    413411}

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r8651 r8653  
    15381538        double         xyz_list[NUMVERTICES][3];
    15391539        double         slope[2];
    1540         double         l1l2l3[3];
     1540        double         basis[3];
    15411541        double         pe_g_gaussian[numdof];
    15421542        GaussTria*     gauss=NULL;
    … …  
    15591559
    15601560                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    1561                 GetNodalFunctions(l1l2l3, gauss);
     1561                GetNodalFunctions(basis, gauss);
    15621562
    15631563                thickness_input->GetParameterValue(&thickness,gauss);
    … …  
    15741574                        for (i=0;i<NUMVERTICES;i++){
    15751575                                for (j=0;j<NDOF2;j++){
    1576                                         pe->values[i*NDOF2+j]+=(-driving_stress_baseline*slope[j]-plastic_stress)*Jdet*gauss->weight*l1l2l3[i];
     1576                                        pe->values[i*NDOF2+j]+=(-driving_stress_baseline*slope[j]-plastic_stress)*Jdet*gauss->weight*basis[i];
    15771577                                }
    15781578                        }
    … …  
    15811581                        for (i=0;i<NUMVERTICES;i++){
    15821582                                for (j=0;j<NDOF2;j++){
    1583                                         pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*l1l2l3[i];
     1583                                        pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i];
    15841584                                }
    15851585                        }
    … …  
    16051605        int         num_responses;
    16061606        double      xyz_list[NUMVERTICES][3];
    1607         double      l1l2l3[3];
     1607        double      basis[3];
    16081608        double      dbasis[NDOF2][NUMVERTICES];
    16091609        double      dH[2];
    … …  
    16281628
    16291629                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    1630                 GetNodalFunctions(l1l2l3, gauss);
     1630                GetNodalFunctions(basis, gauss);
    16311631                GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    16321632
    … …  
    16401640                        case ThicknessAbsMisfitEnum:
    16411641                                weights_input->GetParameterValue(&weight, gauss,resp);
    1642                                 for(i=0;i<numdof;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*l1l2l3[i];
     1642                                for(i=0;i<numdof;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i];
    16431643                                break;
    16441644                        case ThicknessAbsGradientEnum:
    … …  
    16731673        double     vx,vy,vxobs,vyobs,weight;
    16741674        double     xyz_list[NUMVERTICES][3];
    1675         double     l1l2l3[3];
     1675        double     basis[3];
    16761676        GaussTria* gauss=NULL;
    16771677
    … …  
    17121712                vxobs_input->GetParameterValue(&vxobs,gauss);
    17131713                vyobs_input->GetParameterValue(&vyobs,gauss);
    1714                 GetNodalFunctions(l1l2l3, gauss);
     1714                GetNodalFunctions(basis, gauss);
    17151715
    17161716                /*Loop over all requested responses*/
    … …  
    17331733                                                dux=vxobs-vx;
    17341734                                                duy=vyobs-vy;
    1735                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1736                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1735                                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
     1736                                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    17371737                                        }
    17381738                                        break;
    … …  
    17541754                                                dux=scalex*(vxobs-vx);
    17551755                                                duy=scaley*(vyobs-vy);
    1756                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1757                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1756                                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
     1757                                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    17581758                                        }
    17591759                                        break;
    … …  
    17761776                                                dux=scale*vx;
    17771777                                                duy=scale*vy;
    1778                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1779                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1778                                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
     1779                                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    17801780                                        }
    17811781                                        break;
    … …  
    17941794                                                dux=scale*(vxobs-vx);
    17951795                                                duy=scale*(vyobs-vy);
    1796                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1797                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1796                                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
     1797                                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    17981798                                        }
    17991799                                        break;
    … …  
    18111811                                                dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
    18121812                                                duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
    1813                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1814                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1813                                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
     1814                                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    18151815                                        }
    18161816                                        break;
    … …  
    18491849        double     vx,vy,vxobs,vyobs,weight;
    18501850        double     xyz_list[NUMVERTICES][3];
    1851         double     l1l2l3[3];
     1851        double     basis[3];
    18521852        GaussTria* gauss=NULL;
    18531853
    … …  
    18881888                vxobs_input->GetParameterValue(&vxobs,gauss);
    18891889                vyobs_input->GetParameterValue(&vyobs,gauss);
    1890                 GetNodalFunctions(l1l2l3, gauss);
     1890                GetNodalFunctions(basis, gauss);
    18911891
    18921892                /*Loop over all requested responses*/
    … …  
    19101910                                                dux=vxobs-vx;
    19111911                                                duy=vyobs-vy;
    1912                                                 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1913                                                 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1912                                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
     1913                                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    19141914                                        }
    19151915                                        break;
    … …  
    19311931                                                dux=scalex*(vxobs-vx);
    19321932                                                duy=scaley*(vyobs-vy);
    1933                                                 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1934                                                 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1933                                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
     1934                                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    19351935                                        }
    19361936                                        break;
    … …  
    19531953                                                dux=scale*vx;
    19541954                                                duy=scale*vy;
    1955                                                 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1956                                                 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1955                                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
     1956                                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    19571957                                        }
    19581958                                        break;
    … …  
    19711971                                                dux=scale*(vxobs-vx);
    19721972                                                duy=scale*(vyobs-vy);
    1973                                                 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1974                                                 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1973                                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
     1974                                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    19751975                                        }
    19761976                                        break;
    … …  
    19881988                                                dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
    19891989                                                duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
    1990                                                 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i];
    1991                                                 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i];
     1990                                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
     1991                                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    19921992                                        }
    19931993                                        break;
    … …  
    22292229        double     xyz_list[NUMVERTICES][3];
    22302230        double     slope[2];
    2231         double     l1l2l3[3];
     2231        double     basis[3];
    22322232        GaussTria* gauss=NULL;
    22332233
    … …  
    22532253
    22542254                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    2255                 GetNodalFunctions(l1l2l3, gauss);
     2255                GetNodalFunctions(basis, gauss);
    22562256
    22572257                slope_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    22582258
    22592259                if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==BedSlopeXAnalysisEnum)){
    2260                         for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[0]*l1l2l3[i];
     2260                        for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[0]*basis[i];
    22612261                }
    22622262                if ( (analysis_type==SurfaceSlopeYAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){
    2263                         for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[1]*l1l2l3[i];
     2263                        for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[1]*basis[i];
    22642264                }
    22652265        }
    … …  
    29492949        double     Jdet,weight;
    29502950        double     xyz_list[NUMVERTICES][3];
    2951         double     dh1dh3[NDOF2][NUMVERTICES];
     2951        double     dbasis[NDOF2][NUMVERTICES];
    29522952        double     dk[NDOF2];
    29532953        double     grade_g[NUMVERTICES]={0.0};
    … …  
    29672967
    29682968                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    2969                 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
     2969                GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    29702970                weights_input->GetParameterValue(&weight,gauss,weight_index);
    29712971
    … …  
    29742974
    29752975                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    2976                 for (i=0;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
     2976                for (i=0;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
    29772977        }
    29782978        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
    … …  
    30503050        double     bed,thickness,Neff,drag;
    30513051        double     xyz_list[NUMVERTICES][3];
    3052         double     dh1dh3[NDOF2][NUMVERTICES];
     3052        double     dbasis[NDOF2][NUMVERTICES];
    30533053        double     dk[NDOF2];
    30543054        double     grade_g[NUMVERTICES]={0.0};
    30553055        double     grade_g_gaussian[NUMVERTICES];
    3056         double     l1l2l3[3];
     3056        double     basis[3];
    30573057        double     epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    30583058        Friction*  friction=NULL;
    … …  
    30843084
    30853085                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    3086                 GetNodalFunctions(l1l2l3, gauss);
    3087                 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
     3086                GetNodalFunctions(basis, gauss);
     3087                GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    30883088
    30893089                /*Build alpha_complement_list: */
    … …  
    31003100                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    31013101                for (i=0;i<NUMVERTICES;i++){
    3102                         grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*l1l2l3[i];
     3102                        grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
    31033103                }
    31043104               
    … …  
    31283128        double     surface_normal[3],bed_normal[3];
    31293129        double     xyz_list[NUMVERTICES][3];
    3130         double     dh1dh3[NDOF2][NUMVERTICES];
     3130        double     dbasis[NDOF2][NUMVERTICES];
    31313131        double     dk[NDOF2];
    3132         double     l1l2l3[3];
     3132        double     basis[3];
    31333133        double     epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    31343134        double     grade_g[NUMVERTICES]={0.0};
    … …  
    31933193
    31943194                /* Get nodal functions value at gaussian point:*/
    3195                 GetNodalFunctions(l1l2l3, gauss);
     3195                GetNodalFunctions(basis, gauss);
    31963196
    31973197                /*Get nodal functions derivatives*/
    3198                 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
     3198                GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    31993199
    32003200                /*Get k derivative: dk/dx */
    … …  
    32083208                                                -mu    *(2*drag*alpha_complement*(vy - vz*bed_normal[1]*bed_normal[2]))
    32093209                                                -xi    *(2*drag*alpha_complement*(-vx*bed_normal[0]*bed_normal[2]-vy*bed_normal[1]*bed_normal[2]))
    3210                                                 )*Jdet*gauss->weight*l1l2l3[i];
     3210                                                )*Jdet*gauss->weight*basis[i];
    32113211                }
    32123212
    … …  
    32283228        double     Jdet,weight;
    32293229        double     xyz_list[NUMVERTICES][3];
    3230         double     dh1dh3[NDOF2][NUMVERTICES];
     3230        double     dbasis[NDOF2][NUMVERTICES];
    32313231        double     dk[NDOF2];
    32323232        double     grade_g[NUMVERTICES]={0.0};
    … …  
    32473247
    32483248                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    3249                 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
     3249                GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    32503250                weights_input->GetParameterValue(&weight,gauss,weight_index);
    32513251
    … …  
    32553255                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    32563256                for (i=0;i<NUMVERTICES;i++){
    3257                         grade_g[i]+=-weight*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
     3257                        grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
    32583258                        _assert_(!isnan(grade_g[i]));
    32593259                }
    … …  
    32893289        int        doflist1[NUMVERTICES];
    32903290        double     thickness,Jdet;
    3291         double     l1l2l3[3];
     3291        double     basis[3];
    32923292        double     dbasis[NDOF2][NUMVERTICES];
    32933293        double     Dlambda[2],dp[2];
    … …  
    33113311
    33123312                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    3313                 GetNodalFunctions(l1l2l3, gauss);
     3313                GetNodalFunctions(basis, gauss);
    33143314                GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    33153315               
    … …  
    33183318                thickness_input->GetParameterDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
    33193319
    3320                 for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[0]*Jdet*gauss->weight*l1l2l3[i];
     3320                for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[0]*Jdet*gauss->weight*basis[i];
    33213321        }
    33223322
    … …  
    33343334        int        doflist1[NUMVERTICES];
    33353335        double     thickness,Jdet;
    3336         double     l1l2l3[3];
     3336        double     basis[3];
    33373337        double     dbasis[NDOF2][NUMVERTICES];
    33383338        double     Dlambda[2],dp[2];
    … …  
    33563356
    33573357                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    3358                 GetNodalFunctions(l1l2l3, gauss);
     3358                GetNodalFunctions(basis, gauss);
    33593359                GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    33603360
    … …  
    33633363                thickness_input->GetParameterDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
    33643364
    3365                 for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[1]*Jdet*gauss->weight*l1l2l3[i];
     3365                for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[1]*Jdet*gauss->weight*basis[i];
    33663366        }
    33673367        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);

  • issm/trunk/src/c/objects/Elements/TriaRef.cpp

    r8303 r8653  
    7070
    7171        int i;
    72         double dh1dh3[NDOF2][NUMNODES];
     72        double dbasis[NDOF2][NUMNODES];
    7373
    7474        /*Get dh1dh2dh3 in actual coordinate system: */
    75         GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
     75        GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
    7676
    7777        /*Build B: */
    7878        for (i=0;i<NUMNODES;i++){
    79                 *(B+NDOF2*NUMNODES*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];
     79                *(B+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i]; //B[0][NDOF2*i]=dbasis[0][i];
    8080                *(B+NDOF2*NUMNODES*0+NDOF2*i+1)=0;
    8181                *(B+NDOF2*NUMNODES*1+NDOF2*i)=0;
    82                 *(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dh1dh3[1][i];
    83                 *(B+NDOF2*NUMNODES*2+NDOF2*i)=(float).5*dh1dh3[1][i];
    84                 *(B+NDOF2*NUMNODES*2+NDOF2*i+1)=(float).5*dh1dh3[0][i];
     82                *(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i];
     83                *(B+NDOF2*NUMNODES*2+NDOF2*i)=(float).5*dbasis[1][i];
     84                *(B+NDOF2*NUMNODES*2+NDOF2*i+1)=(float).5*dbasis[0][i];
    8585        }
    8686}
    … …  
    101101
    102102        /*Same thing in the actual coordinate system: */
    103         double dh1dh3[NDOF2][NUMNODES];
     103        double dbasis[NDOF2][NUMNODES];
    104104
    105105        /*Get dh1dh2dh3 in actual coordinates system : */
    106         GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
     106        GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
    107107
    108108        /*Build B': */
    109109        for (int i=0;i<NUMNODES;i++){
    110                 *(B+NDOF2*NUMNODES*0+NDOF2*i)=dh1dh3[0][i];
     110                *(B+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i];
    111111                *(B+NDOF2*NUMNODES*0+NDOF2*i+1)=0;
    112112                *(B+NDOF2*NUMNODES*1+NDOF2*i)=0;
    113                 *(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dh1dh3[1][i];
    114                 *(B+NDOF2*NUMNODES*2+NDOF2*i)=0.5*dh1dh3[1][i];
    115                 *(B+NDOF2*NUMNODES*2+NDOF2*i+1)=0.5*dh1dh3[0][i];
     113                *(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i];
     114                *(B+NDOF2*NUMNODES*2+NDOF2*i)=0.5*dbasis[1][i];
     115                *(B+NDOF2*NUMNODES*2+NDOF2*i+1)=0.5*dbasis[0][i];
    116116        }
    117117}
    … …  
    167167         */
    168168
    169         double l1l2l3[NUMNODES];
     169        double basis[NUMNODES];
    170170
    171171        /*Get dh1dh2dh3 in actual coordinate system: */
    172         GetNodalFunctions(&l1l2l3[0],gauss);
     172        GetNodalFunctions(&basis[0],gauss);
    173173
    174174        /*Build B_prog: */
    175175        for (int i=0;i<NUMNODES;i++){
    176                 *(B_prog+NDOF1*NUMNODES*0+NDOF1*i)=l1l2l3[i];
    177                 *(B_prog+NDOF1*NUMNODES*1+NDOF1*i)=l1l2l3[i];
     176                *(B_prog+NDOF1*NUMNODES*0+NDOF1*i)=basis[i];
     177                *(B_prog+NDOF1*NUMNODES*1+NDOF1*i)=basis[i];
    178178        }
    179179}
    … …  
    194194
    195195        /*Same thing in the actual coordinate system: */
    196         double dh1dh3[NDOF2][NUMNODES];
     196        double dbasis[NDOF2][NUMNODES];
    197197
    198198        /*Get dh1dh2dh3 in actual coordinates system : */
    199         GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
     199        GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
    200200
    201201        /*Build B': */
    202202        for (int i=0;i<NUMNODES;i++){
    203                 *(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=2*dh1dh3[0][i];
    204                 *(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=dh1dh3[1][i];
    205                 *(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=dh1dh3[0][i];
    206                 *(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=2*dh1dh3[1][i];
    207                 *(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dh1dh3[1][i];
    208                 *(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dh1dh3[0][i];
     203                *(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=2*dbasis[0][i];
     204                *(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=dbasis[1][i];
     205                *(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=dbasis[0][i];
     206                *(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=2*dbasis[1][i];
     207                *(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dbasis[1][i];
     208                *(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dbasis[0][i];
    209209        }
    210210}
    … …  
    226226
    227227        /*Same thing in the actual coordinate system: */
    228         double dh1dh3[NDOF2][NUMNODES];
     228        double dbasis[NDOF2][NUMNODES];
    229229
    230230        /*Get dh1dh2dh3 in actual coordinates system : */
    231         GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
     231        GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
    232232
    233233        /*Build Bprime: */
    234234        for (int i=0;i<NUMNODES;i++){
    235                 *(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=dh1dh3[0][i];
     235                *(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i];
    236236                *(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=0;
    237237                *(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=0;
    238                 *(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=dh1dh3[1][i];
    239                 *(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dh1dh3[1][i];
    240                 *(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dh1dh3[0][i];
    241                 *(Bprime+NDOF2*NUMNODES*3+NDOF2*i)=dh1dh3[0][i];
    242                 *(Bprime+NDOF2*NUMNODES*3+NDOF2*i+1)=dh1dh3[1][i];
     238                *(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i];
     239                *(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dbasis[1][i];
     240                *(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dbasis[0][i];
     241                *(Bprime+NDOF2*NUMNODES*3+NDOF2*i)=dbasis[0][i];
     242                *(Bprime+NDOF2*NUMNODES*3+NDOF2*i+1)=dbasis[1][i];
    243243        }
    244244}
    … …  
    257257
    258258        /*Same thing in the actual coordinate system: */
    259         double dh1dh3[NDOF2][NUMNODES];
     259        double dbasis[NDOF2][NUMNODES];
    260260
    261261        /*Get dh1dh2dh3 in actual coordinates system : */
    262         GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
     262        GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
    263263
    264264        /*Build B': */
    265265        for (int i=0;i<NUMNODES;i++){
    266                 *(Bprime_prog+NDOF1*NUMNODES*0+NDOF1*i)=dh1dh3[0][i];
    267                 *(Bprime_prog+NDOF1*NUMNODES*1+NDOF1*i)=dh1dh3[1][i];
     266                *(Bprime_prog+NDOF1*NUMNODES*0+NDOF1*i)=dbasis[0][i];
     267                *(Bprime_prog+NDOF1*NUMNODES*1+NDOF1*i)=dbasis[1][i];
    268268        }
    269269}
    … …  
    285285
    286286        int i;
    287         double l1l2l3[3];
    288 
    289         /*Get l1l2l3 in actual coordinate system: */
    290         GetNodalFunctions(l1l2l3,gauss);
     287        double basis[3];
     288
     289        /*Get basis in actual coordinate system: */
     290        GetNodalFunctions(basis,gauss);
    291291
    292292        /*Build L: */
    293293        if(numdof==1){
    294294                for (i=0;i<NUMNODES;i++){
    295                         L[i]=l1l2l3[i];
     295                        L[i]=basis[i];
    296296                }
    297297        }
    298298        else{
    299299                for (i=0;i<NUMNODES;i++){
    300                         *(L+numdof*NUMNODES*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i];
     300                        *(L+numdof*NUMNODES*0+numdof*i)=basis[i]; //L[0][NDOF2*i]=dbasis[0][i];
    301301                        *(L+numdof*NUMNODES*0+numdof*i+1)=0;
    302302                        *(L+numdof*NUMNODES*1+numdof*i)=0;
    303                         *(L+numdof*NUMNODES*1+numdof*i+1)=l1l2l3[i];
     303                        *(L+numdof*NUMNODES*1+numdof*i+1)=basis[i];
    304304                }
    305305        }
    … …  
    394394/*}}}*/
    395395/*FUNCTION TriaRef::GetNodalFunctions{{{1*/
    396 void TriaRef::GetNodalFunctions(double* l1l2l3,GaussTria* gauss){
     396void TriaRef::GetNodalFunctions(double* basis,GaussTria* gauss){
    397397        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    398398
    399         l1l2l3[0]=gauss->coord1;
    400         l1l2l3[1]=gauss->coord2;
    401         l1l2l3[2]=gauss->coord3;
     399        basis[0]=gauss->coord1;
     400        basis[1]=gauss->coord2;
     401        basis[2]=gauss->coord3;
    402402
    403403}
    404404/*}}}*/
    405405/*FUNCTION TriaRef::GetSegmentNodalFunctions{{{1*/
    406 void TriaRef::GetSegmentNodalFunctions(double* l1l2,GaussTria* gauss,int index1,int index2){
     406void TriaRef::GetSegmentNodalFunctions(double* basis,GaussTria* gauss,int index1,int index2){
    407407        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    408408
    … …  
    413413        _assert_(index1>=0 && index1<3);
    414414        _assert_(index2>=0 && index2<3);
    415         l1l2[0]=BasisFunctions[index1];
    416         l1l2[1]=BasisFunctions[index2];
     415        basis[0]=BasisFunctions[index1];
     416        basis[1]=BasisFunctions[index2];
    417417}
    418418/*}}}*/
    419419/*FUNCTION TriaRef::GetNodalFunctionsDerivatives{{{1*/
    420 void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, GaussTria* gauss){
     420void TriaRef::GetNodalFunctionsDerivatives(double* dbasis,double* xyz_list, GaussTria* gauss){
    421421
    422422        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    423423         * actual coordinate system): */
    424424        int       i;
    425         double    dh1dh3_ref[NDOF2][NUMNODES];
     425        double    dbasis_ref[NDOF2][NUMNODES];
    426426        double    Jinv[NDOF2][NDOF2];
    427427
    428428        /*Get derivative values with respect to parametric coordinate system: */
    429         GetNodalFunctionsDerivativesReference(&dh1dh3_ref[0][0], gauss);
     429        GetNodalFunctionsDerivativesReference(&dbasis_ref[0][0], gauss);
    430430
    431431        /*Get Jacobian invert: */
    432432        GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
    433433
    434         /*Build dh1dh3:
     434        /*Build dbasis:
    435435         *
    436436         * [dhi/dx]= Jinv*[dhi/dr]
    … …  
    438438         */
    439439        for (i=0;i<NUMNODES;i++){
    440                 dh1dh3[NUMNODES*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
    441                 dh1dh3[NUMNODES*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
     440                dbasis[NUMNODES*0+i]=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i];
     441                dbasis[NUMNODES*1+i]=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i];
    442442        }
    443443
    … …  
    467467
    468468        /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian
    469          * point specified by gauss_l1l2l3:
     469         * point specified by gauss_basis:
    470470         *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
    471471         *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
    … …  
    475475
    476476        /*Nodal Derivatives*/
    477         double dh1dh3[2][3]; //nodal derivative functions in actual coordinate system.
     477        double dbasis[2][3]; //nodal derivative functions in actual coordinate system.
    478478
    479479        /*Get dh1dh2dh3 in actual coordinate system: */
    480         GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss);
     480        GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list, gauss);
    481481
    482482        /*Assign values*/
    483         *(p+0)=plist[0]*dh1dh3[0][0]+plist[1]*dh1dh3[0][1]+plist[2]*dh1dh3[0][2];
    484         *(p+1)=plist[0]*dh1dh3[1][0]+plist[1]*dh1dh3[1][1]+plist[2]*dh1dh3[1][2];
     483        *(p+0)=plist[0]*dbasis[0][0]+plist[1]*dbasis[0][1]+plist[2]*dbasis[0][2];
     484        *(p+1)=plist[0]*dbasis[1][0]+plist[1]*dbasis[1][1]+plist[2]*dbasis[1][2];
    485485
    486486}
    … …  
    493493
    494494        /*nodal functions annd output: */
    495         double l1l2l3[3];
     495        double basis[3];
    496496
    497497        /*Get nodal functions*/
    498         GetNodalFunctions(l1l2l3, gauss);
     498        GetNodalFunctions(basis, gauss);
    499499
    500500        /*Get parameter*/
    501         *p=l1l2l3[0]*plist[0]+l1l2l3[1]*plist[1]+l1l2l3[2]*plist[2];
    502 }
    503 /*}}}*/
     501        *p=basis[0]*plist[0]+basis[1]*plist[1]+basis[2]*plist[2];
     502}
     503/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.