Changeset 15470


Ignore:
Timestamp:
07/09/13 09:13:41 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: cosmetics

Location:
issm/trunk-jpl/src/c/classes/Elements
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp

    r15458 r15470  
    7171
    7272        /*Build B: */
    73         for (int i=0;i<NUMNODESP1;i++){
    74                 *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i];
    75                 *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.0;
    76 
    77                 *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0.0;
    78                 *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
    79 
    80                 *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dbasis[1][i];
    81                 *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][i];
     73        for(int i=0;i<NUMNODESP1;i++){
     74                B[NDOF2*NUMNODESP1*0+NDOF2*i+0] = dbasis[0][i];
     75                B[NDOF2*NUMNODESP1*0+NDOF2*i+1] = 0.;
     76
     77                B[NDOF2*NUMNODESP1*1+NDOF2*i+0] = 0.;
     78                B[NDOF2*NUMNODESP1*1+NDOF2*i+1] = dbasis[1][i];
     79
     80                B[NDOF2*NUMNODESP1*2+NDOF2*i+0] = .5*dbasis[1][i];
     81                B[NDOF2*NUMNODESP1*2+NDOF2*i+1] = .5*dbasis[0][i];
    8282        }
    8383}
    … …  
    9797         */
    9898
    99         int    i;
    100         IssmDouble dh1dh7[3][NUMNODESMINI];
    101         IssmDouble l1l6[NUMNODESP1];
    102 
    103         /*Get dh1dh6 in actual coordinate system: */
    104         GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
    105         GetNodalFunctionsP1(l1l6, gauss);
     99        int i;
     100        IssmDouble dbasismini[3][NUMNODESMINI];
     101        IssmDouble basis[NUMNODESP1];
     102
     103        /*Get dbasis in actual coordinate system: */
     104        GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
     105        GetNodalFunctionsP1(basis,gauss);
    106106
    107107        /*Build B: */
    108         for (i=0;i<NUMNODESMINI;i++){
    109                 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh7[0][i]; //B[0][NDOF4*i]=dh1dh6[0][i];
    110                 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0;
    111                 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
    112                 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0;
    113                 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh7[1][i];
    114                 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
    115                 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0.5*dh1dh7[1][i];
    116                 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0.5*dh1dh7[0][i];
    117                 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=0;
    118                 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=0;
    119                 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=0;
    120                 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;
    121         }
    122 
    123         for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
    124                 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
    125                 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
    126                 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
    127                 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=l1l6[i];
     108        for(i=0;i<NUMNODESMINI;i++){
     109                B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dbasismini[0][i];
     110                B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = 0.;
     111                B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
     112                B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = 0.;
     113                B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dbasismini[1][i];
     114                B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
     115                B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = 0.5*dbasismini[1][i];
     116                B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = 0.5*dbasismini[0][i];
     117                B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = 0.;
     118                B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = 0.;
     119                B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = 0.;
     120                B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2] = 0.;
     121        }
     122
     123        for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
     124                B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0;
     125                B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0;
     126                B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0;
     127                B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3] = basis[i];
    128128        }
    129129}
    … …  
    151151        /*Build B: */
    152152        for (int i=0;i<NUMNODESP1;i++){
    153                 *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i];
    154                 *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.0;
    155 
    156                 *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0.0;
    157                 *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
    158 
    159                 *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dbasis[1][i];
    160                 *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][i];
    161 
    162                 *(B+NDOF2*NUMNODESP1*3+NDOF2*i)=(float).5*dbasis[2][i];
    163                 *(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
    164 
    165                 *(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
    166                 *(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=(float).5*dbasis[2][i];
     153                B[NDOF2*NUMNODESP1*0+NDOF2*i+0] = dbasis[0][i];
     154                B[NDOF2*NUMNODESP1*0+NDOF2*i+1] = 0.;
     155
     156                B[NDOF2*NUMNODESP1*1+NDOF2*i+0] = 0.;
     157                B[NDOF2*NUMNODESP1*1+NDOF2*i+1] = dbasis[1][i];
     158
     159                B[NDOF2*NUMNODESP1*2+NDOF2*i+0] = .5*dbasis[1][i];
     160                B[NDOF2*NUMNODESP1*2+NDOF2*i+1] = .5*dbasis[0][i];
     161
     162                B[NDOF2*NUMNODESP1*3+NDOF2*i+0] = .5*dbasis[2][i];
     163                B[NDOF2*NUMNODESP1*3+NDOF2*i+1] = 0.;
     164
     165                B[NDOF2*NUMNODESP1*4+NDOF2*i+0] = 0.;
     166                B[NDOF2*NUMNODESP1*4+NDOF2*i+1] = .5*dbasis[2][i];
    167167        }
    168168
    … …  
    189189
    190190        /*Build BPrime: */
    191         for (int i=0;i<NUMNODESP1;i++){
    192                 *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=2.0*dbasis[0][i];
    193                 *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=dbasis[1][i];
    194 
    195                 *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=dbasis[0][i];
    196                 *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2.0*dbasis[1][i];
    197 
    198                 *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=dbasis[1][i];
    199                 *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dbasis[0][i];
    200 
    201                 *(B+NDOF2*NUMNODESP1*3+NDOF2*i)=dbasis[2][i];
    202                 *(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
    203 
    204                 *(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
    205                 *(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=dbasis[2][i];
     191        for(int i=0;i<NUMNODESP1;i++){
     192                B[NDOF2*NUMNODESP1*0+NDOF2*i+0]=2.*dbasis[0][i];
     193                B[NDOF2*NUMNODESP1*0+NDOF2*i+1]=dbasis[1][i];
     194
     195                B[NDOF2*NUMNODESP1*1+NDOF2*i+0]=dbasis[0][i];
     196                B[NDOF2*NUMNODESP1*1+NDOF2*i+1]=2.*dbasis[1][i];
     197
     198                B[NDOF2*NUMNODESP1*2+NDOF2*i+0]=dbasis[1][i];
     199                B[NDOF2*NUMNODESP1*2+NDOF2*i+1]=dbasis[0][i];
     200
     201                B[NDOF2*NUMNODESP1*3+NDOF2*i+0]=dbasis[2][i];
     202                B[NDOF2*NUMNODESP1*3+NDOF2*i+1]=0.;
     203
     204                B[NDOF2*NUMNODESP1*4+NDOF2*i+0]=0.;
     205                B[NDOF2*NUMNODESP1*4+NDOF2*i+1]=dbasis[2][i];
    206206        }
    207207}
    … …  
    221221
    222222        int    i;
    223         IssmDouble dh1dh7[3][NUMNODESMINI];
    224 
    225         /*Get dh1dh6 in actual coordinate system: */
    226         GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
     223        IssmDouble dbasismini[3][NUMNODESMINI];
     224
     225        /*Get dbasis in actual coordinate system: */
     226        GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
    227227
    228228        /*Build Bprime: */
    229         for (i=0;i<NUMNODESMINI;i++){
    230                 *(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=2*dh1dh7[0][i]; //Bprime[0][NDOF4*i]=dh1dh6[0][i];
    231                 *(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=dh1dh7[1][i];
    232                 *(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
    233                 *(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=dh1dh7[0][i];
    234                 *(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=2*dh1dh7[1][i];
    235                 *(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
    236                 *(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=dh1dh7[1][i];
    237                 *(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=dh1dh7[0][i];
    238                 *(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=0;
    239         }
    240 
    241         for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
    242                 *(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
    243                 *(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
    244                 *(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
     229        for(i=0;i<NUMNODESMINI;i++){
     230                Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = 2.*dbasismini[0][i];
     231                Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = dbasismini[1][i];
     232                Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
     233                Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = dbasismini[0][i];
     234                Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = 2.*dbasismini[1][i];
     235                Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
     236                Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = dbasismini[1][i];
     237                Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = dbasismini[0][i];
     238                Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = 0.;
     239        }
     240
     241        for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
     242                Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0.;
     243                Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0.;
     244                Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0.;
    245245        }
    246246
    … …  
    266266        int i;
    267267
    268         IssmDouble dh1dh7[3][NUMNODESMINI];
    269         IssmDouble l1l6[NUMNODESP1];
    270 
    271         /*Get dh1dh7 in actual coordinate system: */
    272         GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
    273         GetNodalFunctionsP1(l1l6, gauss);
     268        IssmDouble dbasismini[3][NUMNODESMINI];
     269        IssmDouble basis[NUMNODESP1];
     270
     271        /*Get dbasismini in actual coordinate system: */
     272        GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
     273        GetNodalFunctionsP1(basis, gauss);
    274274
    275275        /*Build B: */
    276         for (i=0;i<NUMNODESMINI;i++){
    277                 B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dh1dh7[0][i+0]; //B[0][NDOF4*i+0] = dh1dh6[0][i+0];
     276        for(i=0;i<NUMNODESMINI;i++){
     277                B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dbasismini[0][i+0]; //B[0][NDOF4*i+0] = dbasis[0][i+0];
    278278                B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = 0.;
    279279                B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
    280280                B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = 0.;
    281                 B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dh1dh7[1][i+0];
     281                B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dbasismini[1][i+0];
    282282                B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
    283283                B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = 0.;
    284284                B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = 0.;
    285                 B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = dh1dh7[2][i+0];
    286                 B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = (float).5*dh1dh7[1][i+0];
    287                 B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = (float).5*dh1dh7[0][i+0];
     285                B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = dbasismini[2][i+0];
     286                B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = .5*dbasismini[1][i+0];
     287                B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = .5*dbasismini[0][i+0];
    288288                B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2] = 0.;
    289                 B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+0] = (float).5*dh1dh7[2][i+0];
     289                B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+0] = .5*dbasismini[2][i+0];
    290290                B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1] = 0.;
    291                 B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2] = (float).5*dh1dh7[0][i+0];
     291                B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2] = .5*dbasismini[0][i+0];
    292292                B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+0] = 0.;
    293                 B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1] = (float).5*dh1dh7[2][i+0];
    294                 B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2] = (float).5*dh1dh7[1][i+0];
     293                B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1] = .5*dbasismini[2][i+0];
     294                B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2] = .5*dbasismini[1][i+0];
    295295                B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+0] = 0.;
    296296                B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1] = 0.;
    297297                B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2] = 0.;
    298                 B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+0] = dh1dh7[0][i+0];
    299                 B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1] = dh1dh7[1][i+0];
    300                 B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2] = dh1dh7[2][i+0];
    301         }
    302 
    303         for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
    304                 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
    305                 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
    306                 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
    307                 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0;
    308                 *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0;
    309                 *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0;
    310                 *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=l1l6[i];
    311                 *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=0;
     298                B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+0] = dbasismini[0][i+0];
     299                B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1] = dbasismini[1][i+0];
     300                B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2] = dbasismini[2][i+0];
     301        }
     302
     303        for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
     304                B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0.;
     305                B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0.;
     306                B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0.;
     307                B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3] = 0.;
     308                B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3] = 0.;
     309                B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3] = 0.;
     310                B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3] = basis[i];
     311                B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3] = 0.;
    312312        }
    313313
    … …  
    332332
    333333        int i;
    334 
    335         IssmDouble dh1dh6[3][NUMNODESP1];
    336         IssmDouble l1l6[NUMNODESP1];
    337 
    338         /*Get dh1dh7 in actual coordinate system: */
    339         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
    340         GetNodalFunctionsP1(&l1l6[0], gauss);
     334        IssmDouble dbasis[3][NUMNODESP1];
     335        IssmDouble basis[NUMNODESP1];
     336
     337        /*Get dbasismini in actual coordinate system: */
     338        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
     339        GetNodalFunctionsP1(&basis[0], gauss);
    341340
    342341        /*Build B: */
    343         for (i=0;i<NUMNODESP1;i++){
    344                 *(B+(NDOF4*NUMNODESP1)*0+NDOF4*i)=dh1dh6[0][i]; //B[0][NDOF4*i]=dh1dh6[0][i];
    345                 *(B+(NDOF4*NUMNODESP1)*0+NDOF4*i+1)=0.;
    346                 *(B+(NDOF4*NUMNODESP1)*0+NDOF4*i+2)=0.;
    347                 *(B+(NDOF4*NUMNODESP1)*1+NDOF4*i)=0.;
    348                 *(B+(NDOF4*NUMNODESP1)*1+NDOF4*i+1)=dh1dh6[1][i];
    349                 *(B+(NDOF4*NUMNODESP1)*1+NDOF4*i+2)=0.;
    350                 *(B+(NDOF4*NUMNODESP1)*2+NDOF4*i)=0.;
    351                 *(B+(NDOF4*NUMNODESP1)*2+NDOF4*i+1)=0.;
    352                 *(B+(NDOF4*NUMNODESP1)*2+NDOF4*i+2)=dh1dh6[2][i];
    353                 *(B+(NDOF4*NUMNODESP1)*3+NDOF4*i)=.5*dh1dh6[1][i];
    354                 *(B+(NDOF4*NUMNODESP1)*3+NDOF4*i+1)=.5*dh1dh6[0][i];
    355                 *(B+(NDOF4*NUMNODESP1)*3+NDOF4*i+2)=0.;
    356                 *(B+(NDOF4*NUMNODESP1)*4+NDOF4*i)=.5*dh1dh6[2][i];
    357                 *(B+(NDOF4*NUMNODESP1)*4+NDOF4*i+1)=0.;
    358                 *(B+(NDOF4*NUMNODESP1)*4+NDOF4*i+2)=.5*dh1dh6[0][i];
    359                 *(B+(NDOF4*NUMNODESP1)*5+NDOF4*i)=0.;
    360                 *(B+(NDOF4*NUMNODESP1)*5+NDOF4*i+1)=.5*dh1dh6[2][i];
    361                 *(B+(NDOF4*NUMNODESP1)*5+NDOF4*i+2)=.5*dh1dh6[1][i];
    362                 *(B+(NDOF4*NUMNODESP1)*6+NDOF4*i)=0.;
    363                 *(B+(NDOF4*NUMNODESP1)*6+NDOF4*i+1)=0.;
    364                 *(B+(NDOF4*NUMNODESP1)*6+NDOF4*i+2)=0.;
    365                 *(B+(NDOF4*NUMNODESP1)*7+NDOF4*i)=dh1dh6[0][i];
    366                 *(B+(NDOF4*NUMNODESP1)*7+NDOF4*i+1)=dh1dh6[1][i];
    367                 *(B+(NDOF4*NUMNODESP1)*7+NDOF4*i+2)=dh1dh6[2][i];
    368                 _assert_((NDOF4*NUMNODESP1)*7+NDOF4*i+2<8*24);
    369         }
    370 
    371         for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
    372                 B[(NDOF4*NUMNODESP1)*0+NDOF4*i+3]=0.;
    373                 B[(NDOF4*NUMNODESP1)*1+NDOF4*i+3]=0.;
    374                 B[(NDOF4*NUMNODESP1)*2+NDOF4*i+3]=0.;
    375                 B[(NDOF4*NUMNODESP1)*3+NDOF4*i+3]=0.;
    376                 B[(NDOF4*NUMNODESP1)*4+NDOF4*i+3]=0.;
    377                 B[(NDOF4*NUMNODESP1)*5+NDOF4*i+3]=0.;
    378                 B[(NDOF4*NUMNODESP1)*6+NDOF4*i+3]=l1l6[i];
    379                 B[(NDOF4*NUMNODESP1)*7+NDOF4*i+3]=0.;
    380                 _assert_(((NDOF4*NUMNODESP1)*7+NDOF4*i+3)<8*24);
     342        for(i=0;i<NUMNODESP1;i++){
     343                B[(NDOF4*NUMNODESP1)*0+NDOF4*i+0] = dbasis[0][i];
     344                B[(NDOF4*NUMNODESP1)*0+NDOF4*i+1] = 0.;
     345                B[(NDOF4*NUMNODESP1)*0+NDOF4*i+2] = 0.;
     346                B[(NDOF4*NUMNODESP1)*1+NDOF4*i+0] = 0.;
     347                B[(NDOF4*NUMNODESP1)*1+NDOF4*i+1] = dbasis[1][i];
     348                B[(NDOF4*NUMNODESP1)*1+NDOF4*i+2] = 0.;
     349                B[(NDOF4*NUMNODESP1)*2+NDOF4*i+0] = 0.;
     350                B[(NDOF4*NUMNODESP1)*2+NDOF4*i+1] = 0.;
     351                B[(NDOF4*NUMNODESP1)*2+NDOF4*i+2] = dbasis[2][i];
     352                B[(NDOF4*NUMNODESP1)*3+NDOF4*i+0] = .5*dbasis[1][i];
     353                B[(NDOF4*NUMNODESP1)*3+NDOF4*i+1] = .5*dbasis[0][i];
     354                B[(NDOF4*NUMNODESP1)*3+NDOF4*i+2] = 0.;
     355                B[(NDOF4*NUMNODESP1)*4+NDOF4*i+0] = .5*dbasis[2][i];
     356                B[(NDOF4*NUMNODESP1)*4+NDOF4*i+1] = 0.;
     357                B[(NDOF4*NUMNODESP1)*4+NDOF4*i+2] = .5*dbasis[0][i];
     358                B[(NDOF4*NUMNODESP1)*5+NDOF4*i+0] = 0.;
     359                B[(NDOF4*NUMNODESP1)*5+NDOF4*i+1] = .5*dbasis[2][i];
     360                B[(NDOF4*NUMNODESP1)*5+NDOF4*i+2] = .5*dbasis[1][i];
     361                B[(NDOF4*NUMNODESP1)*6+NDOF4*i+0] = 0.;
     362                B[(NDOF4*NUMNODESP1)*6+NDOF4*i+1] = 0.;
     363                B[(NDOF4*NUMNODESP1)*6+NDOF4*i+2] = 0.;
     364                B[(NDOF4*NUMNODESP1)*7+NDOF4*i+0] = dbasis[0][i];
     365                B[(NDOF4*NUMNODESP1)*7+NDOF4*i+1] = dbasis[1][i];
     366                B[(NDOF4*NUMNODESP1)*7+NDOF4*i+2] = dbasis[2][i];
     367        }
     368
     369        for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
     370                B[(NDOF4*NUMNODESP1)*0+NDOF4*i+3] = 0.;
     371                B[(NDOF4*NUMNODESP1)*1+NDOF4*i+3] = 0.;
     372                B[(NDOF4*NUMNODESP1)*2+NDOF4*i+3] = 0.;
     373                B[(NDOF4*NUMNODESP1)*3+NDOF4*i+3] = 0.;
     374                B[(NDOF4*NUMNODESP1)*4+NDOF4*i+3] = 0.;
     375                B[(NDOF4*NUMNODESP1)*5+NDOF4*i+3] = 0.;
     376                B[(NDOF4*NUMNODESP1)*6+NDOF4*i+3] = basis[i];
     377                B[(NDOF4*NUMNODESP1)*7+NDOF4*i+3] = 0.;
    381378        }
    382379
    … …  
    402399
    403400        int i;
    404         IssmDouble dh1dh7[3][NUMNODESMINI];
    405         IssmDouble l1l6[NUMNODESP1];
    406 
    407         /*Get dh1dh7 in actual coordinate system: */
    408         GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
    409         GetNodalFunctionsP1(l1l6, gauss);
     401        IssmDouble dbasismini[3][NUMNODESMINI];
     402        IssmDouble basis[NUMNODESP1];
     403
     404        /*Get dbasismini in actual coordinate system: */
     405        GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
     406        GetNodalFunctionsP1(basis, gauss);
    410407
    411408        /*B_primeuild B_prime: */
    412         for (i=0;i<NUMNODESMINI;i++){
    413                 *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh7[0][i]; //B_prime[0][NDOF4*i]=dh1dh6[0][i];
    414                 *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0;
    415                 *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
    416                 *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0;
    417                 *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh7[1][i];
    418                 *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
    419                 *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0;
    420                 *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0;
    421                 *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=dh1dh7[2][i];
    422                 *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=dh1dh7[1][i];
    423                 *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=dh1dh7[0][i];
    424                 *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;
    425                 *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i)=dh1dh7[2][i];
    426                 *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1)=0;
    427                 *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2)=dh1dh7[0][i];
    428                 *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i)=0;
    429                 *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1)=dh1dh7[2][i];
    430                 *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2)=dh1dh7[1][i];
    431                 *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i)=dh1dh7[0][i];
    432                 *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1)=dh1dh7[1][i];
    433                 *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2)=dh1dh7[2][i];
    434                 *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i)=0;
    435                 *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1)=0;
    436                 *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2)=0;
    437         }
    438 
    439         for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
    440                 *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
    441                 *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
    442                 *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
    443                 *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0;
    444                 *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0;
    445                 *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0;
    446                 *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=0;
    447                 *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=l1l6[i];
     409        for(i=0;i<NUMNODESMINI;i++){
     410                B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dbasismini[0][i];
     411                B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = 0.;
     412                B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
     413                B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = 0.;
     414                B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dbasismini[1][i];
     415                B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
     416                B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = 0.;
     417                B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = 0.;
     418                B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = dbasismini[2][i];
     419                B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = dbasismini[1][i];
     420                B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = dbasismini[0][i];
     421                B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2] = 0.;
     422                B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+0] = dbasismini[2][i];
     423                B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1] = 0.;
     424                B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2] = dbasismini[0][i];
     425                B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+0] = 0.;
     426                B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1] = dbasismini[2][i];
     427                B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2] = dbasismini[1][i];
     428                B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+0] = dbasismini[0][i];
     429                B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1] = dbasismini[1][i];
     430                B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2] = dbasismini[2][i];
     431                B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+0] = 0.;
     432                B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1] = 0.;
     433                B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2] = 0.;
     434        }
     435
     436        for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
     437                B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0.;
     438                B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0.;
     439                B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0.;
     440                B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3] = 0.;
     441                B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3] = 0.;
     442                B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3] = 0.;
     443                B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3] = 0.;
     444                B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3] = basis[i];
    448445        }
    449446
    … …  
    469466
    470467        int i;
    471         IssmDouble dh1dh6[3][NUMNODESP1];
    472         IssmDouble l1l6[NUMNODESP1];
    473 
    474         /*Get dh1dh7 in actual coordinate system: */
    475         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
    476         GetNodalFunctionsP1(l1l6, gauss);
     468        IssmDouble dbasis[3][NUMNODESP1];
     469        IssmDouble basis[NUMNODESP1];
     470
     471        /*Get dbasismini in actual coordinate system: */
     472        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
     473        GetNodalFunctionsP1(basis, gauss);
    477474
    478475        /*B_primeuild B_prime: */
    479         for (i=0;i<NUMNODESP1;i++){
    480                 *(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i)=dh1dh6[0][i]; //B_prime[0][NDOF4*i]=dh1dh6[0][i];
    481                 *(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i+1)=0.;
    482                 *(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i+2)=0.;
    483                 *(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i)=0.;
    484                 *(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i+1)=dh1dh6[1][i];
    485                 *(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i+2)=0.;
    486                 *(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i)=0.;
    487                 *(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i+1)=0.;
    488                 *(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i+2)=dh1dh6[2][i];
    489                 *(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i)=dh1dh6[1][i];
    490                 *(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i+1)=dh1dh6[0][i];
    491                 *(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i+2)=0.;
    492                 *(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i)=dh1dh6[2][i];
    493                 *(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i+1)=0.;
    494                 *(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i+2)=dh1dh6[0][i];
    495                 *(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i)=0.;
    496                 *(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i+1)=dh1dh6[2][i];
    497                 *(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i+2)=dh1dh6[1][i];
    498                 *(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i)=dh1dh6[0][i];
    499                 *(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i+1)=dh1dh6[1][i];
    500                 *(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i+2)=dh1dh6[2][i];
    501                 *(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i)=0.;
    502                 *(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i+1)=0.;
    503                 *(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i+2)=0.;
    504                 _assert_((NDOF4*NUMNODESP1)*7+NDOF4*i+2<8*24);
    505         }
    506 
    507         for (i=0;i<NUMNODESP1;i++){ //last column
    508                 *(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i+3)=0.;
    509                 *(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i+3)=0.;
    510                 *(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i+3)=0.;
    511                 *(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i+3)=0.;
    512                 *(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i+3)=0.;
    513                 *(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i+3)=0.;
    514                 *(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i+3)=0.;
    515                 *(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i+3)=l1l6[i];
    516                 _assert_((NDOF4*NUMNODESP1)*7+NDOF4*i+2<8*24);
     476        for(i=0;i<NUMNODESP1;i++){
     477                B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+0] = dbasis[0][i];
     478                B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+1] = 0.;
     479                B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+2] = 0.;
     480                B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+0] = 0.;
     481                B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+1] = dbasis[1][i];
     482                B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+2] = 0.;
     483                B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+0] = 0.;
     484                B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+1] = 0.;
     485                B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+2] = dbasis[2][i];
     486                B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+0] = dbasis[1][i];
     487                B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+1] = dbasis[0][i];
     488                B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+2] = 0.;
     489                B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+0] = dbasis[2][i];
     490                B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+1] = 0.;
     491                B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+2] = dbasis[0][i];
     492                B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+0] = 0.;
     493                B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+1] = dbasis[2][i];
     494                B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+2] = dbasis[1][i];
     495                B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+0] = dbasis[0][i];
     496                B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+1] = dbasis[1][i];
     497                B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+2] = dbasis[2][i];
     498                B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+0] = 0.;
     499                B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+1] = 0.;
     500                B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+2] = 0.;
     501        }
     502
     503        for(i=0;i<NUMNODESP1;i++){ //last column
     504                B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+3] = 0.;
     505                B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+3] = 0.;
     506                B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+3] = 0.;
     507                B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+3] = 0.;
     508                B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+3] = 0.;
     509                B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+3] = 0.;
     510                B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+3] = 0.;
     511                B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+3] = basis[i];
    517512        }
    518513
    … …  
    533528
    534529        /*Same thing in the actual coordinate system: */
    535         IssmDouble l1l6[6];
     530        IssmDouble basis[6];
    536531
    537532        /*Get dh1dh2dh3 in actual coordinates system : */
    538         GetNodalFunctionsP1(l1l6, gauss);
     533        GetNodalFunctionsP1(basis, gauss);
    539534
    540535        /*Build B': */
    541         for (int i=0;i<NUMNODESP1;i++){
    542                 *(B_advec+NDOF1*NUMNODESP1*0+NDOF1*i)=l1l6[i];
    543                 *(B_advec+NDOF1*NUMNODESP1*1+NDOF1*i)=l1l6[i];
    544                 *(B_advec+NDOF1*NUMNODESP1*2+NDOF1*i)=l1l6[i];
     536        for(int i=0;i<NUMNODESP1;i++){
     537                B_advec[NDOF1*NUMNODESP1*0+NDOF1*i] = basis[i];
     538                B_advec[NDOF1*NUMNODESP1*1+NDOF1*i] = basis[i];
     539                B_advec[NDOF1*NUMNODESP1*2+NDOF1*i] = basis[i];
    545540        }
    546541}
    … …  
    560555
    561556        /*Same thing in the actual coordinate system: */
    562         IssmDouble dh1dh6[3][NUMNODESP1];
     557        IssmDouble dbasis[3][NUMNODESP1];
    563558
    564559        /*Get dh1dh2dh3 in actual coordinates system : */
    565         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
     560        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
    566561
    567562        /*Build B': */
    568         for (int i=0;i<NUMNODESP1;i++){
    569                 *(B_conduct+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i];
    570                 *(B_conduct+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i];
    571                 *(B_conduct+NDOF1*NUMNODESP1*2+NDOF1*i)=dh1dh6[2][i];
     563        for(int i=0;i<NUMNODESP1;i++){
     564                B_conduct[NDOF1*NUMNODESP1*0+NDOF1*i] = dbasis[0][i];
     565                B_conduct[NDOF1*NUMNODESP1*1+NDOF1*i] = dbasis[1][i];
     566                B_conduct[NDOF1*NUMNODESP1*2+NDOF1*i] = dbasis[2][i];
    572567        }
    573568}
    … …  
    578573                where hi is the interpolation function for node i.*/
    579574
    580         int i;
    581         IssmDouble dh1dh6[3][NUMNODESP1];
    582 
    583         /*Get dh1dh6 in actual coordinate system: */
    584         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
     575        /*Get dbasis in actual coordinate system: */
     576        IssmDouble dbasis[3][NUMNODESP1];
     577        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
    585578
    586579        /*Build B: */
    587         for (i=0;i<NUMNODESP1;i++){
    588                 B[i]=dh1dh6[2][i]; 
     580        for(int i=0;i<NUMNODESP1;i++){
     581                B[i] = dbasis[2][i]; 
    589582        }
    590583
    … …  
    604597         */
    605598
    606         /*Same thing in the actual coordinate system: */
    607         IssmDouble dh1dh6[3][NUMNODESP1];
    608 
    609         /*Get dh1dh2dh3 in actual coordinates system : */
    610         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
     599        /*Get nodal function derivatives in actual coordinates system : */
     600        IssmDouble dbasis[3][NUMNODESP1];
     601        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
    611602
    612603        /*Build B': */
    613         for (int i=0;i<NUMNODESP1;i++){
    614                 *(Bprime_advec+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i];
    615                 *(Bprime_advec+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i];
    616                 *(Bprime_advec+NDOF1*NUMNODESP1*2+NDOF1*i)=dh1dh6[2][i];
     604        for(int i=0;i<NUMNODESP1;i++){
     605                Bprime_advec[NDOF1*NUMNODESP1*0+NDOF1*i] = dbasis[0][i];
     606                Bprime_advec[NDOF1*NUMNODESP1*1+NDOF1*i] = dbasis[1][i];
     607                Bprime_advec[NDOF1*NUMNODESP1*2+NDOF1*i] = dbasis[2][i];
    617608        }
    618609}
    … …  
    620611/*FUNCTION PentaRef::GetBprimeVert{{{*/
    621612void PentaRef::GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
    622         /* Compute Bprime  matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for node i*/
    623613
    624614        GetNodalFunctionsP1(B, gauss);
    … …  
    638628         **/
    639629
     630        /*Get basis in actual coordinate system: */
    640631        IssmDouble basis[6];
    641 
    642         /*Get l1l6 in actual coordinate system: */
    643632        GetNodalFunctionsP1(&basis[0],gauss);
    644633
    645634        for(int i=0;i<NUMNODESP1;i++){
    646                 B[2*NUMNODESP1*0+2*i+0]=basis[i];
    647                 B[2*NUMNODESP1*0+2*i+1]=0.;
    648                 B[2*NUMNODESP1*1+2*i+0]=0.;
    649                 B[2*NUMNODESP1*1+2*i+1]=basis[i];
     635                B[2*NUMNODESP1*0+2*i+0] = basis[i];
     636                B[2*NUMNODESP1*0+2*i+1] = 0.;
     637                B[2*NUMNODESP1*1+2*i+0] = 0.;
     638                B[2*NUMNODESP1*1+2*i+1] = basis[i];
    650639        }
    651640}
    … …  
    653642/*FUNCTION PentaRef::GetLStokes{{{*/
    654643void PentaRef::GetLStokes(IssmDouble* LStokes, GaussPenta* gauss){
    655         /*
    656          * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
     644        /* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    657645         * For node i, Li can be expressed in the actual coordinate system
    658646         * by:
    … …  
    673661
    674662        /*Build LStokes: */
    675         for (int i=0;i<3;i++){
    676                 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+0)=l1l2l3[i];
    677                 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0.;
    678                 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0.;
    679                 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0.;
    680 
    681                 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+0)=0.;
    682                 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
    683                 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0.;
    684                 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0.;
     663        for(int i=0;i<3;i++){
     664                LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
     665                LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
     666                LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
     667                LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
     668
     669                LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
     670                LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
     671                LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
     672                LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
    685673        }
    686674}
    … …  
    688676/*FUNCTION PentaRef::GetLprimeStokes {{{*/
    689677void PentaRef::GetLprimeStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
    690 
    691         /*
    692          * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
     678        /* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
    693679         * For node i, Lpi can be expressed in the actual coordinate system
    694680         * by:
    … …  
    728714
    729715        IssmDouble l1l2l3[NUMNODESP1_2d];
    730         IssmDouble dh1dh6[3][NUMNODESP1];
     716        IssmDouble dbasis[3][NUMNODESP1];
    731717
    732718        /*Get l1l2l3 in actual coordinate system: */
    … …  
    735721        l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    736722
    737         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
     723        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
    738724
    739725        /*Build LprimeStokes: */
    740         for (i=0;i<3;i++){
    741                 *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
    742                 *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
    743                 *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
    744                 *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
    745                 *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
    746                 *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
    747                 *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
    748                 *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
    749                 *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
    750                 *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
    751                 *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=0;
    752                 *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
    753                 *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
    754                 *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
    755                 *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=0;
    756                 *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
    757                 *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=0;
    758                 *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
    759                 *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=l1l2l3[i];
    760                 *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0;
    761                 *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
    762                 *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=0;
    763                 *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=l1l2l3[i];
    764                 *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0;
    765                 *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=0;
    766                 *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
    767                 *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=dh1dh6[2][i];
    768                 *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=0;
    769                 *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
    770                 *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=0;
    771                 *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=dh1dh6[2][i];
    772                 *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=0;
    773                 *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i)=0;
    774                 *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+1)=0;
    775                 *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+2)=dh1dh6[2][i];
    776                 *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+3)=0;
    777                 *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i)=dh1dh6[2][i];
    778                 *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+1)=0;
    779                 *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+2)=dh1dh6[0][i];
    780                 *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+3)=0;
    781                 *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i)=0;
    782                 *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+1)=dh1dh6[2][i];
    783                 *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+2)=dh1dh6[1][i];
    784                 *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+3)=0;
    785                 *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i)=0;
    786                 *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+1)=0;
    787                 *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+2)=0;
    788                 *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+3)=l1l2l3[i];
    789                 *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i)=0;
    790                 *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+1)=0;
    791                 *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+2)=0;
    792                 *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+3)=l1l2l3[i];
    793                 *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i)=0;
    794                 *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+1)=0;
    795                 *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+2)=0;
    796                 *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+3)=l1l2l3[i];
     726        for(int i=0;i<3;i++){
     727                LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0]  = l1l2l3[i];
     728                LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1]  = 0.;
     729                LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2]  = 0.;
     730                LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3]  = 0.;
     731                LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0]  = 0.;
     732                LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1]  = l1l2l3[i];
     733                LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2]  = 0.;
     734                LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3]  = 0.;
     735                LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0]  = l1l2l3[i];
     736                LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1]  = 0.;
     737                LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2]  = 0.;
     738                LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3]  = 0.;
     739                LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0]  = 0.;
     740                LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1]  = l1l2l3[i];
     741                LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2]  = 0.;
     742                LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3]  = 0.;
     743                LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0]  = 0.;
     744                LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1]  = 0.;
     745                LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+2]  = l1l2l3[i];
     746                LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+3]  = 0.;
     747                LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0]  = 0.;
     748                LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1]  = 0.;
     749                LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+2]  = l1l2l3[i];
     750                LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+3]  = 0.;
     751                LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0]  = 0.;
     752                LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1]  = 0.;
     753                LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+2]  = dbasis[2][i];
     754                LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+3]  = 0.;
     755                LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0]  = 0.;
     756                LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1]  = 0.;
     757                LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+2]  = dbasis[2][i];
     758                LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+3]  = 0.;
     759                LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+0]  = 0.;
     760                LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+1]  = 0.;
     761                LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+2]  = dbasis[2][i];
     762                LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+3]  = 0.;
     763                LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+0]  = dbasis[2][i];
     764                LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+1]  = 0.;
     765                LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+2]  = dbasis[0][i];
     766                LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+3]  = 0.;
     767                LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+0] = 0.;
     768                LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+1] = dbasis[2][i];
     769                LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+2] = dbasis[1][i];
     770                LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+3] = 0.;
     771                LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+0] = 0.;
     772                LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+1] = 0.;
     773                LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+2] = 0.;
     774                LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+3] = l1l2l3[i];
     775                LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+0] = 0.;
     776                LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+1] = 0.;
     777                LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+2] = 0.;
     778                LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+3] = l1l2l3[i];
     779                LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+0] = 0.;
     780                LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+1] = 0;
     781                LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+2] = 0;
     782                LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+3] = l1l2l3[i];
    797783        }
    798784}
    … …  
    815801         */
    816802
    817         int i;
    818803        int num_dof=2;
    819 
    820804        IssmDouble l1l2l3[NUMNODESP1_2d];
    821805
    … …  
    826810
    827811        /*Build LStokes: */
    828         for (i=0;i<3;i++){
    829                 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
    830                 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
    831                 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
    832                 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
    833                 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
    834                 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
    835                 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
    836                 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
    837                 *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=l1l2l3[i];
    838                 *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
    839                 *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
    840                 *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=l1l2l3[i];
    841                 *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=l1l2l3[i];
    842                 *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
    843                 *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
    844                 *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=l1l2l3[i];
    845 
     812        for(int i=0;i<3;i++){
     813                LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
     814                LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0;
     815                LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0;
     816                LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
     817                LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = l1l2l3[i];
     818                LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0;
     819                LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0;
     820                LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = l1l2l3[i];
     821                LStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = l1l2l3[i];
     822                LStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0;
     823                LStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0;
     824                LStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = l1l2l3[i];
     825                LStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = l1l2l3[i];
     826                LStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0;
     827                LStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0;
     828                LStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = l1l2l3[i];
    846829        }
    847830}
    … …  
    849832/*FUNCTION PentaRef::GetLprimeMacAyealStokes {{{*/
    850833void PentaRef::GetLprimeMacAyealStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
    851 
    852         /*
    853          * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
     834        /* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
    854835         * For node i, Lpi can be expressed in the actual coordinate system
    855836         * by:
    … …  
    864845         * where h is the interpolation function for node i.
    865846         */
    866         int i;
    867847        int num_dof=4;
    868 
    869848        IssmDouble l1l2l3[NUMNODESP1_2d];
    870         IssmDouble dh1dh6[3][NUMNODESP1];
     849        IssmDouble dbasis[3][NUMNODESP1];
    871850
    872851        /*Get l1l2l3 in actual coordinate system: */
    … …  
    875854        l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    876855
    877         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
     856        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
    878857
    879858        /*Build LprimeStokes: */
    880         for (i=0;i<3;i++){
    881                 *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
    882                 *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
    883                 *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
    884                 *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
    885                 *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
    886                 *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
    887                 *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
    888                 *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
    889                 *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0;
    890                 *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
    891                 *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i];
    892                 *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
    893                 *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
    894                 *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0;
    895                 *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i];
    896                 *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
    897                 *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=0;
    898                 *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
    899                 *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=dh1dh6[2][i];
    900                 *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0;
    901                 *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
    902                 *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=0;
    903                 *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=dh1dh6[2][i];
    904                 *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0;
    905                 *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=0;
    906                 *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
    907                 *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=0;
    908                 *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=l1l2l3[i];
    909                 *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
    910                 *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=0;
    911                 *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=0;
    912                 *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=l1l2l3[i];
     859        for(int i=0;i<3;i++){
     860                LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
     861                LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
     862                LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
     863                LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
     864                LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
     865                LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
     866                LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
     867                LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
     868                LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
     869                LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
     870                LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = l1l2l3[i];
     871                LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
     872                LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
     873                LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
     874                LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = l1l2l3[i];
     875                LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
     876                LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = 0.;
     877                LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0.;
     878                LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+2] = dbasis[2][i];
     879                LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+3] = 0.;
     880                LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0.;
     881                LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = 0.;
     882                LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+2] = dbasis[2][i];
     883                LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+3] = 0.;
     884                LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = 0.;
     885                LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0.;
     886                LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+2] = 0.;
     887                LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+3] = l1l2l3[i];
     888                LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0.;
     889                LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = 0.;
     890                LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+2] = 0.;
     891                LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+3] = l1l2l3[i];
    913892        }
    914893}
    … …  
    916895/*FUNCTION PentaRef::GetLStokesMacAyeal {{{*/
    917896void PentaRef::GetLStokesMacAyeal(IssmDouble* LStokes, GaussPenta* gauss){
    918         /*
    919          * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
     897        /* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    920898         * For node i, Li can be expressed in the actual coordinate system
    921899         * by:
    … …  
    927905         */
    928906
    929         int i;
    930907        int num_dof=4;
    931 
    932908        IssmDouble l1l2l3[NUMNODESP1_2d];
    933909
    … …  
    938914
    939915        /*Build LStokes: */
    940         for (i=0;i<3;i++){
    941                 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
    942                 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
    943                 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
    944                 *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
    945                 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
    946                 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
    947                 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
    948                 *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
    949                 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0;
    950                 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
    951                 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i];
    952                 *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
    953                 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
    954                 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0;
    955                 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i];
    956                 *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
    957 
     916        for(int i=0;i<3;i++){
     917                LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
     918                LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
     919                LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
     920                LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
     921                LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
     922                LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
     923                LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
     924                LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
     925                LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
     926                LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
     927                LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = l1l2l3[i];
     928                LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
     929                LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
     930                LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
     931                LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = l1l2l3[i];
     932                LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
    958933        }
    959934}
    … …  
    961936/*FUNCTION PentaRef::GetLprimeStokesMacAyeal {{{*/
    962937void PentaRef::GetLprimeStokesMacAyeal(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
    963 
    964         /*
    965          * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
     938        /* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
    966939         * For node i, Lpi can be expressed in the actual coordinate system
    967940         * by:
    … …  
    972945         * where h is the interpolation function for node i.
    973946         */
    974         int i;
    975947        int num_dof=2;
    976 
    977948        IssmDouble l1l2l3[NUMNODESP1_2d];
    978         IssmDouble dh1dh6[3][NUMNODESP1];
     949        IssmDouble dbasis[3][NUMNODESP1];
    979950
    980951        /*Get l1l2l3 in actual coordinate system: */
    … …  
    982953        l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
    983954        l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    984 
    985         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
     955        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
    986956
    987957        /*Build LprimeStokes: */
    988         for (i=0;i<3;i++){
    989                 *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
    990                 *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
    991                 *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
    992                 *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
    993                 *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
    994                 *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
    995                 *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
    996                 *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
     958        for(int i=0;i<3;i++){
     959                LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
     960                LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
     961                LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
     962                LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
     963                LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = l1l2l3[i];
     964                LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
     965                LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
     966                LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = l1l2l3[i];
    997967        }
    998968}
    … …  
    1000970/*FUNCTION PentaRef::GetJacobian {{{*/
    1001971void PentaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussPenta* gauss){
    1002 
    1003972        /*The Jacobian is constant over the element, discard the gaussian points.
    1004973         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    1005974
    1006         IssmDouble A1,A2,A3; //area coordinates
    1007         IssmDouble xi,eta,zi; //parametric coordinates
    1008 
     975        IssmDouble A1,A2,A3;  // area coordinates
     976        IssmDouble xi,eta,zi; // parametric coordinates
    1009977        IssmDouble x1,x2,x3,x4,x5,x6;
    1010978        IssmDouble y1,y2,y3,y4,y5,y6;
    … …  
    1012980
    1013981        /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
    1014         A1=gauss->coord1;
    1015         A2=gauss->coord2;
    1016         A3=gauss->coord3;
    1017 
     982        A1  = gauss->coord1;
     983        A2  = gauss->coord2;
     984        A3  = gauss->coord3;
    1018985        xi  = A2-A1;
    1019986        eta = SQRT3*A3;
    1020987        zi  = gauss->coord4;
    1021988
    1022         x1=*(xyz_list+3*0+0);
    1023         x2=*(xyz_list+3*1+0);
    1024         x3=*(xyz_list+3*2+0);
    1025         x4=*(xyz_list+3*3+0);
    1026         x5=*(xyz_list+3*4+0);
    1027         x6=*(xyz_list+3*5+0);
    1028 
    1029         y1=*(xyz_list+3*0+1);
    1030         y2=*(xyz_list+3*1+1);
    1031         y3=*(xyz_list+3*2+1);
    1032         y4=*(xyz_list+3*3+1);
    1033         y5=*(xyz_list+3*4+1);
    1034         y6=*(xyz_list+3*5+1);
    1035 
    1036         z1=*(xyz_list+3*0+2);
    1037         z2=*(xyz_list+3*1+2);
    1038         z3=*(xyz_list+3*2+2);
    1039         z4=*(xyz_list+3*3+2);
    1040         z5=*(xyz_list+3*4+2);
    1041         z6=*(xyz_list+3*5+2);
    1042 
    1043         *(J+NDOF3*0+0)=0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
    1044         *(J+NDOF3*1+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);
    1045         *(J+NDOF3*2+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);
    1046 
    1047         *(J+NDOF3*0+1)=0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
    1048         *(J+NDOF3*1+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);
    1049         *(J+NDOF3*2+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);
    1050 
    1051         *(J+NDOF3*0+2)=0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
    1052         *(J+NDOF3*1+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);
    1053         *(J+NDOF3*2+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);
    1054 
     989        x1=xyz_list[3*0+0];
     990        x2=xyz_list[3*1+0];
     991        x3=xyz_list[3*2+0];
     992        x4=xyz_list[3*3+0];
     993        x5=xyz_list[3*4+0];
     994        x6=xyz_list[3*5+0];
     995
     996        y1=xyz_list[3*0+1];
     997        y2=xyz_list[3*1+1];
     998        y3=xyz_list[3*2+1];
     999        y4=xyz_list[3*3+1];
     1000        y5=xyz_list[3*4+1];
     1001        y6=xyz_list[3*5+1];
     1002
     1003        z1=xyz_list[3*0+2];
     1004        z2=xyz_list[3*1+2];
     1005        z3=xyz_list[3*2+2];
     1006        z4=xyz_list[3*3+2];
     1007        z5=xyz_list[3*4+2];
     1008        z6=xyz_list[3*5+2];
     1009
     1010        J[NDOF3*0+0] = 0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
     1011        J[NDOF3*1+0] = SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);
     1012        J[NDOF3*2+0] = SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);
     1013
     1014        J[NDOF3*0+1] = 0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
     1015        J[NDOF3*1+1] = SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);
     1016        J[NDOF3*2+1] = SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);
     1017
     1018        J[NDOF3*0+2] = 0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
     1019        J[NDOF3*1+2] = SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);
     1020        J[NDOF3*2+2] = SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);
    10551021}
    10561022/*}}}*/
    … …  
    10751041         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    10761042
    1077         IssmDouble x1,x2,x3,y1,y2,y3,z1,z2,z3;
    1078 
    1079         x1=*(xyz_list+3*0+0);
    1080         y1=*(xyz_list+3*0+1);
    1081         z1=*(xyz_list+3*0+2);
    1082         x2=*(xyz_list+3*1+0);
    1083         y2=*(xyz_list+3*1+1);
    1084         z2=*(xyz_list+3*1+2);
    1085         x3=*(xyz_list+3*2+0);
    1086         y3=*(xyz_list+3*2+1);
    1087         z3=*(xyz_list+3*2+2);
     1043        IssmDouble x1=xyz_list[3*0+0];
     1044        IssmDouble y1=xyz_list[3*0+1];
     1045        IssmDouble z1=xyz_list[3*0+2];
     1046        IssmDouble x2=xyz_list[3*1+0];
     1047        IssmDouble y2=xyz_list[3*1+1];
     1048        IssmDouble z2=xyz_list[3*1+2];
     1049        IssmDouble x3=xyz_list[3*2+0];
     1050        IssmDouble y3=xyz_list[3*2+1];
     1051        IssmDouble z3=xyz_list[3*2+2];
    10881052
    10891053        /*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
    1090         *Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5);
     1054        *Jdet=SQRT3/6.*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2),0.5);
    10911055        if(*Jdet<0) _error_("negative jacobian determinant!");
    10921056}
    … …  
    10971061         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    10981062
    1099         IssmDouble x1,x2,y1,y2,z1,z2;
    1100 
    1101         x1=*(xyz_list+3*0+0);
    1102         y1=*(xyz_list+3*0+1);
    1103         z1=*(xyz_list+3*0+2);
    1104         x2=*(xyz_list+3*1+0);
    1105         y2=*(xyz_list+3*1+1);
    1106         z2=*(xyz_list+3*1+2);
    1107 
    1108         *Jdet=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.) + pow(z2-z1,2.));
     1063        IssmDouble x1=xyz_list[3*0+0];
     1064        IssmDouble y1=xyz_list[3*0+1];
     1065        IssmDouble z1=xyz_list[3*0+2];
     1066        IssmDouble x2=xyz_list[3*1+0];
     1067        IssmDouble y2=xyz_list[3*1+1];
     1068        IssmDouble z2=xyz_list[3*1+2];
     1069
     1070        *Jdet=.5*sqrt(pow(x2-x1,2) + pow(y2-y1,2) + pow(z2-z1,2));
    11091071        if(*Jdet<0) _error_("negative jacobian determinant!");
    11101072
    … …  
    12751237/*}}}*/
    12761238/*FUNCTION PentaRef::GetNodalFunctionsMINIDerivatives{{{*/
    1277 void PentaRef::GetNodalFunctionsMINIDerivatives(IssmDouble* dh1dh7,IssmDouble* xyz_list, GaussPenta* gauss){
     1239void PentaRef::GetNodalFunctionsMINIDerivatives(IssmDouble* dbasismini,IssmDouble* xyz_list, GaussPenta* gauss){
    12781240
    12791241        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    12801242         * actual coordinate system): */
    12811243
    1282         IssmDouble    dh1dh7_ref[3][NUMNODESMINI];
     1244        IssmDouble    dbasismini_ref[3][NUMNODESMINI];
    12831245        IssmDouble    Jinv[3][3];
    12841246
    12851247        /*Get derivative values with respect to parametric coordinate system: */
    1286         GetNodalFunctionsMINIDerivativesReference(&dh1dh7_ref[0][0], gauss);
     1248        GetNodalFunctionsMINIDerivativesReference(&dbasismini_ref[0][0], gauss);
    12871249
    12881250        /*Get Jacobian invert: */
    12891251        GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
    12901252
    1291         /*Build dh1dh6:
     1253        /*Build dbasis:
    12921254         *
    12931255         * [dhi/dx]= Jinv'*[dhi/dr]
    … …  
    12971259
    12981260        for(int i=0;i<NUMNODESMINI;i++){
    1299                 *(dh1dh7+NUMNODESMINI*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];
    1300                 *(dh1dh7+NUMNODESMINI*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];
    1301                 *(dh1dh7+NUMNODESMINI*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i];
     1261                *(dbasismini+NUMNODESMINI*0+i)=Jinv[0][0]*dbasismini_ref[0][i]+Jinv[0][1]*dbasismini_ref[1][i]+Jinv[0][2]*dbasismini_ref[2][i];
     1262                *(dbasismini+NUMNODESMINI*1+i)=Jinv[1][0]*dbasismini_ref[0][i]+Jinv[1][1]*dbasismini_ref[1][i]+Jinv[1][2]*dbasismini_ref[2][i];
     1263                *(dbasismini+NUMNODESMINI*2+i)=Jinv[2][0]*dbasismini_ref[0][i]+Jinv[2][1]*dbasismini_ref[1][i]+Jinv[2][2]*dbasismini_ref[2][i];
    13021264        }
    13031265
    … …  
    13421304/*}}}*/
    13431305/*FUNCTION PentaRef::GetNodalFunctionsP1 {{{*/
    1344 void PentaRef::GetNodalFunctionsP1(IssmDouble* l1l6, GaussPenta* gauss){
     1306void PentaRef::GetNodalFunctionsP1(IssmDouble* basis, GaussPenta* gauss){
    13451307        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    13461308
    1347         l1l6[0]=gauss->coord1*(1-gauss->coord4)/2.0;
    1348         l1l6[1]=gauss->coord2*(1-gauss->coord4)/2.0;
    1349         l1l6[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    1350         l1l6[3]=gauss->coord1*(1+gauss->coord4)/2.0;
    1351         l1l6[4]=gauss->coord2*(1+gauss->coord4)/2.0;
    1352         l1l6[5]=gauss->coord3*(1+gauss->coord4)/2.0;
     1309        basis[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     1310        basis[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     1311        basis[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     1312        basis[3]=gauss->coord1*(1+gauss->coord4)/2.0;
     1313        basis[4]=gauss->coord2*(1+gauss->coord4)/2.0;
     1314        basis[5]=gauss->coord3*(1+gauss->coord4)/2.0;
    13531315
    13541316}
    13551317/*}}}*/
    13561318/*FUNCTION PentaRef::GetNodalFunctionsP1Derivatives {{{*/
    1357 void PentaRef::GetNodalFunctionsP1Derivatives(IssmDouble* dh1dh6,IssmDouble* xyz_list, GaussPenta* gauss){
     1319void PentaRef::GetNodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){
    13581320
    13591321        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    13601322         * actual coordinate system): */
    1361         IssmDouble    dh1dh6_ref[NDOF3][NUMNODESP1];
     1323        IssmDouble    dbasis_ref[NDOF3][NUMNODESP1];
    13621324        IssmDouble    Jinv[NDOF3][NDOF3];
    13631325
    13641326        /*Get derivative values with respect to parametric coordinate system: */
    1365         GetNodalFunctionsP1DerivativesReference(&dh1dh6_ref[0][0], gauss);
     1327        GetNodalFunctionsP1DerivativesReference(&dbasis_ref[0][0], gauss);
    13661328
    13671329        /*Get Jacobian invert: */
    … …  
    13761338
    13771339        for (int i=0;i<NUMNODESP1;i++){
    1378                 *(dh1dh6+NUMNODESP1*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i];
    1379                 *(dh1dh6+NUMNODESP1*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i];
    1380                 *(dh1dh6+NUMNODESP1*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i];
     1340                *(dbasis+NUMNODESP1*0+i)=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i]+Jinv[0][2]*dbasis_ref[2][i];
     1341                *(dbasis+NUMNODESP1*1+i)=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i]+Jinv[1][2]*dbasis_ref[2][i];
     1342                *(dbasis+NUMNODESP1*2+i)=Jinv[2][0]*dbasis_ref[0][i]+Jinv[2][1]*dbasis_ref[1][i]+Jinv[2][2]*dbasis_ref[2][i];
    13811343        }
    13821344
    … …  
    14681430
    14691431        /*intermediary*/
    1470         IssmDouble l1l6[6];
     1432        IssmDouble basis[6];
    14711433
    14721434        /*nodal functions: */
    1473         GetNodalFunctionsP1(&l1l6[0],gauss);
     1435        GetNodalFunctionsP1(&basis[0],gauss);
    14741436
    14751437        /*Assign output pointers:*/
    1476         *pvalue=l1l6[0]*plist[0]+l1l6[1]*plist[1]+l1l6[2]*plist[2]+l1l6[3]*plist[3]+l1l6[4]*plist[4]+l1l6[5]*plist[5];
     1438        *pvalue=basis[0]*plist[0]+basis[1]*plist[1]+basis[2]*plist[2]+basis[3]*plist[3]+basis[4]*plist[4]+basis[5]*plist[5];
    14771439
    14781440}
    … …  
    14801442/*FUNCTION PentaRef::GetInputDerivativeValue{{{*/
    14811443void PentaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussPenta* gauss){
    1482         /*From node values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_coord:
     1444        /*From node values of parameter p (p_list[0], p_list[1], p_list[2],
     1445         * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
     1446         * gaussian point specified by gauss_coord:
    14831447         *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
    14841448         *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
    … …  
    14871451         *   p is a vector of size 3x1 already allocated.
    14881452         */
    1489         IssmDouble dh1dh6[3][NUMNODESP1];
     1453        IssmDouble dbasis[3][NUMNODESP1];
    14901454
    14911455        /*Get nodal funnctions derivatives in actual coordinate system: */
    1492         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
     1456        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
    14931457
    14941458        /*Assign output*/
    1495         p[0]=plist[0]*dh1dh6[0][0]+plist[1]*dh1dh6[0][1]+plist[2]*dh1dh6[0][2]+plist[3]*dh1dh6[0][3]+plist[4]*dh1dh6[0][4]+plist[5]*dh1dh6[0][5];
    1496         p[1]=plist[0]*dh1dh6[1][0]+plist[1]*dh1dh6[1][1]+plist[2]*dh1dh6[1][2]+plist[3]*dh1dh6[1][3]+plist[4]*dh1dh6[1][4]+plist[5]*dh1dh6[1][5];
    1497         p[2]=plist[0]*dh1dh6[2][0]+plist[1]*dh1dh6[2][1]+plist[2]*dh1dh6[2][2]+plist[3]*dh1dh6[2][3]+plist[4]*dh1dh6[2][4]+plist[5]*dh1dh6[2][5];
     1459        p[0]=plist[0]*dbasis[0][0]+plist[1]*dbasis[0][1]+plist[2]*dbasis[0][2]+plist[3]*dbasis[0][3]+plist[4]*dbasis[0][4]+plist[5]*dbasis[0][5];
     1460        p[1]=plist[0]*dbasis[1][0]+plist[1]*dbasis[1][1]+plist[2]*dbasis[1][2]+plist[3]*dbasis[1][3]+plist[4]*dbasis[1][4]+plist[5]*dbasis[1][5];
     1461        p[2]=plist[0]*dbasis[2][0]+plist[1]*dbasis[2][1]+plist[2]*dbasis[2][2]+plist[3]*dbasis[2][3]+plist[4]*dbasis[2][4]+plist[5]*dbasis[2][5];
    14981462
    14991463}

  • issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp

    r15446 r15470  
    7272        /*Build B: */
    7373        for(int i=0;i<numnodes;i++){
    74                 B[numnodes*0+i]=dbasis[0*numnodes+i];
    75                 B[numnodes*1+i]=dbasis[1*numnodes+i];
     74                B[numnodes*0+i] = dbasis[0*numnodes+i];
     75                B[numnodes*1+i] = dbasis[1*numnodes+i];
    7676        }
    7777
    … …  
    102102        /*Build B: */
    103103        for(int i=0;i<numnodes;i++){
    104                 B[NDOF2*numnodes*0+NDOF2*i+0]=dbasis[0*numnodes+i];
    105                 B[NDOF2*numnodes*0+NDOF2*i+1]=0.;
    106                 B[NDOF2*numnodes*1+NDOF2*i+0]=0.;
    107                 B[NDOF2*numnodes*1+NDOF2*i+1]=dbasis[1*numnodes+i];
    108                 B[NDOF2*numnodes*2+NDOF2*i+0]=.5*dbasis[1*numnodes+i];
    109                 B[NDOF2*numnodes*2+NDOF2*i+1]=.5*dbasis[0*numnodes+i];
     104                B[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];
     105                B[NDOF2*numnodes*0+NDOF2*i+1] = 0.;
     106                B[NDOF2*numnodes*1+NDOF2*i+0] = 0.;
     107                B[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];
     108                B[NDOF2*numnodes*2+NDOF2*i+0] = .5*dbasis[1*numnodes+i];
     109                B[NDOF2*numnodes*2+NDOF2*i+1] = .5*dbasis[0*numnodes+i];
    110110        }
    111111
    … …  
    136136
    137137        /*Build B': */
    138         for (int i=0;i<numnodes;i++){
    139                 B[NDOF2*numnodes*0+NDOF2*i+0]=dbasis[0*numnodes+i];
    140                 B[NDOF2*numnodes*0+NDOF2*i+1]=0.;
    141                 B[NDOF2*numnodes*1+NDOF2*i+0]=0.;
    142                 B[NDOF2*numnodes*1+NDOF2*i+1]=dbasis[1*numnodes+i];
    143                 B[NDOF2*numnodes*2+NDOF2*i+0]=0.5*dbasis[1*numnodes+i];
    144                 B[NDOF2*numnodes*2+NDOF2*i+1]=0.5*dbasis[0*numnodes+i];
     138        for(int i=0;i<numnodes;i++){
     139                B[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];
     140                B[NDOF2*numnodes*0+NDOF2*i+1] = 0.;
     141                B[NDOF2*numnodes*1+NDOF2*i+0] = 0.;
     142                B[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];
     143                B[NDOF2*numnodes*2+NDOF2*i+0] = 0.5*dbasis[1*numnodes+i];
     144                B[NDOF2*numnodes*2+NDOF2*i+1] = 0.5*dbasis[0*numnodes+i];
    145145        }
    146146
    … …  
    221221
    222222        /*Build B: */
    223         for (int i=0;i<numnodes;i++){
    224                 B[numnodes*0+i]=basis[i];
    225                 B[numnodes*1+i]=basis[i];
     223        for(int i=0;i<numnodes;i++){
     224                B[numnodes*0+i] = basis[i];
     225                B[numnodes*1+i] = basis[i];
    226226        }
    227227
    … …  
    252252
    253253        /*Build B': */
    254         for (int i=0;i<numnodes;i++){
     254        for(int i=0;i<numnodes;i++){
    255255                Bprime[NDOF2*numnodes*0+NDOF2*i+0] = 2.*dbasis[0*numnodes+i];
    256256                Bprime[NDOF2*numnodes*0+NDOF2*i+1] =    dbasis[1*numnodes+i];
    … …  
    288288        /*Build Bprime: */
    289289        for(int i=0;i<numnodes;i++){
    290                 Bprime[NDOF2*numnodes*0+NDOF2*i+0]=dbasis[0*numnodes+i];
    291                 Bprime[NDOF2*numnodes*0+NDOF2*i+1]=0.;
    292                 Bprime[NDOF2*numnodes*1+NDOF2*i+0]=0.;
    293                 Bprime[NDOF2*numnodes*1+NDOF2*i+1]=dbasis[1*numnodes+i];
    294                 Bprime[NDOF2*numnodes*2+NDOF2*i+0]=dbasis[1*numnodes+i];
    295                 Bprime[NDOF2*numnodes*2+NDOF2*i+1]=dbasis[0*numnodes+i];
    296                 Bprime[NDOF2*numnodes*3+NDOF2*i+0]=dbasis[0*numnodes+i];
    297                 Bprime[NDOF2*numnodes*3+NDOF2*i+1]=dbasis[1*numnodes+i];
     290                Bprime[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];
     291                Bprime[NDOF2*numnodes*0+NDOF2*i+1] = 0.;
     292                Bprime[NDOF2*numnodes*1+NDOF2*i+0] = 0.;
     293                Bprime[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];
     294                Bprime[NDOF2*numnodes*2+NDOF2*i+0] = dbasis[1*numnodes+i];
     295                Bprime[NDOF2*numnodes*2+NDOF2*i+1] = dbasis[0*numnodes+i];
     296                Bprime[NDOF2*numnodes*3+NDOF2*i+0] = dbasis[0*numnodes+i];
     297                Bprime[NDOF2*numnodes*3+NDOF2*i+1] = dbasis[1*numnodes+i];
    298298        }
    299299
    … …  
    323323        /*Build B': */
    324324        for(int i=0;i<numnodes;i++){
    325                 Bprime[numnodes*0+i]=dbasis[0*numnodes+i];
    326                 Bprime[numnodes*1+i]=dbasis[1*numnodes+i];
     325                Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
     326                Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
    327327        }
    328328
    … …  
    352352        /*Build L: */
    353353        for(int i=0;i<numnodes;i++){
    354                 B[2*numnodes*0+2*i+0]=basis[i];
    355                 B[2*numnodes*0+2*i+1]=0.;
    356                 B[2*numnodes*1+2*i+0]=0.;
    357                 B[2*numnodes*1+2*i+1]=basis[i];
     354                B[2*numnodes*0+2*i+0] = basis[i];
     355                B[2*numnodes*0+2*i+1] = 0.;
     356                B[2*numnodes*1+2*i+0] = 0.;
     357                B[2*numnodes*1+2*i+1] = basis[i];
    358358        }
    359359
    … …  
    508508         */
    509509        for(int i=0;i<numnodes;i++){
    510                 dbasis[numnodes*0+i]=Jinv[0][0]*dbasis_ref[0*numnodes+i]+Jinv[0][1]*dbasis_ref[1*numnodes+i];
    511                 dbasis[numnodes*1+i]=Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i];
     510                dbasis[numnodes*0+i] = Jinv[0][0]*dbasis_ref[0*numnodes+i]+Jinv[0][1]*dbasis_ref[1*numnodes+i];
     511                dbasis[numnodes*1+i] = Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i];
    512512        }
    513513
    … …  
    527527                case P1Enum: case P1DGEnum:
    528528                        /*Nodal function 1*/
    529                         dbasis[NUMNODESP1*0+0]=-0.5;
    530                         dbasis[NUMNODESP1*1+0]=-1.0/(2.0*SQRT3);
     529                        dbasis[NUMNODESP1*0+0] = -0.5;
     530                        dbasis[NUMNODESP1*1+0] = -1.0/(2.0*SQRT3);
    531531                        /*Nodal function 2*/
    532                         dbasis[NUMNODESP1*0+1]=0.5;
    533                         dbasis[NUMNODESP1*1+1]=-1.0/(2.0*SQRT3);
     532                        dbasis[NUMNODESP1*0+1] = 0.5;
     533                        dbasis[NUMNODESP1*1+1] = -1.0/(2.0*SQRT3);
    534534                        /*Nodal function 3*/
    535                         dbasis[NUMNODESP1*0+2]=0;
    536                         dbasis[NUMNODESP1*1+2]=1.0/SQRT3;
     535                        dbasis[NUMNODESP1*0+2] = 0;
     536                        dbasis[NUMNODESP1*1+2] = 1.0/SQRT3;
    537537                        return;
    538538                case P2Enum:
    539539                        /*Nodal function 1*/
    540                         dbasis[NUMNODESP2*0+0]=-2.*gauss->coord1 + 0.5;
    541                         dbasis[NUMNODESP2*1+0]=-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.;
     540                        dbasis[NUMNODESP2*0+0] = -2.*gauss->coord1 + 0.5;
     541                        dbasis[NUMNODESP2*1+0] = -2.*SQRT3/3.*gauss->coord1 + SQRT3/6.;
    542542                        /*Nodal function 2*/
    543                         dbasis[NUMNODESP2*0+1]=+2.*gauss->coord2 - 0.5;
    544                         dbasis[NUMNODESP2*1+1]=-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.;
     543                        dbasis[NUMNODESP2*0+1] = +2.*gauss->coord2 - 0.5;
     544                        dbasis[NUMNODESP2*1+1] = -2.*SQRT3/3.*gauss->coord2 + SQRT3/6.;
    545545                        /*Nodal function 3*/
    546                         dbasis[NUMNODESP2*0+2]=0.;
    547                         dbasis[NUMNODESP2*1+2]=+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.;
     546                        dbasis[NUMNODESP2*0+2] = 0.;
     547                        dbasis[NUMNODESP2*1+2] = +4.*SQRT3/3.*gauss->coord3 - SQRT3/3.;
    548548                        /*Nodal function 4*/
    549                         dbasis[NUMNODESP2*0+3]=+2.*gauss->coord3;
    550                         dbasis[NUMNODESP2*1+3]=+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3;
     549                        dbasis[NUMNODESP2*0+3] = +2.*gauss->coord3;
     550                        dbasis[NUMNODESP2*1+3] = +4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3;
    551551                        /*Nodal function 5*/
    552                         dbasis[NUMNODESP2*0+4]=-2.*gauss->coord3;
    553                         dbasis[NUMNODESP2*1+4]=+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3;
     552                        dbasis[NUMNODESP2*0+4] = -2.*gauss->coord3;
     553                        dbasis[NUMNODESP2*1+4] = +4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3;
    554554                        /*Nodal function 6*/
    555                         dbasis[NUMNODESP2*0+5]=2.*(gauss->coord1-gauss->coord2);
    556                         dbasis[NUMNODESP2*1+5]=-2.*SQRT3/3.*(gauss->coord1+gauss->coord2);
     555                        dbasis[NUMNODESP2*0+5] = 2.*(gauss->coord1-gauss->coord2);
     556                        dbasis[NUMNODESP2*1+5] = -2.*SQRT3/3.*(gauss->coord1+gauss->coord2);
    557557                        return;
    558558                default:
    … …  
    590590        /*Assign values*/
    591591        xDelete<IssmDouble>(dbasis);
    592         *(p+0)=dpx;
    593         *(p+1)=dpy;
     592        p[0]=dpx;
     593        p[1]=dpy;
    594594
    595595}
Note: See TracChangeset for help on using the changeset viewer.