Changeset 15323


Ignore:
Timestamp:
06/22/13 13:25:33 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: extended dynamic GetB for other solutions + cosmetics

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r15314 r15323  
    5656         * For node i, Bi can be expressed in the actual coordinate system
    5757         * by:
    58          *       Bi=[ dh/dx ]
    59          *          [ dh/dy ]
    60          * where h is the interpolation function for node i.
     58         *       Bi=[ dN/dx ]
     59         *          [ dN/dy ]
     60         * where N is the interpolation function for node i.
    6161         *
    6262         * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
    … …  
    6666        int numnodes = this->NumberofNodes();
    6767
    68         /*Get nodal functions*/
     68        /*Get nodal functions derivatives*/
    6969        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
    7070        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    … …  
    9696        int numnodes = this->NumberofNodes();
    9797
    98         /*Get nodal functions*/
     98        /*Get nodal functions derivatives*/
    9999        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
    100100        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    … …  
    120120         * For node i, Bi can be expressed in the actual coordinate system
    121121         * by:
    122          *       Bi=[   dh/dx         0     ]
    123          *          [       0       dh/dy   ]
    124          *          [  1/2*dh/dy  1/2*dh/dx ]
    125          * where h is the interpolation function for node i.
    126          *
    127          * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
    128          */
    129 
    130         /*Same thing in the actual coordinate system: */
    131         IssmDouble dbasis[NDOF2][NUMNODESP1];
    132 
    133         /*Get dh1dh2dh3 in actual coordinates system : */
    134         GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
     122         *       Bi=[   dN/dx         0     ]
     123         *          [       0       dN/dy   ]
     124         *          [  1/2*dN/dy  1/2*dN/dx ]
     125         * where N is the interpolation function for node i.
     126         *
     127         * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
     128         */
     129
     130        /*Fetch number of nodes for this finite element*/
     131        int numnodes = this->NumberofNodes();
     132
     133        /*Get nodal functions derivatives*/
     134        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     135        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    135136
    136137        /*Build B': */
    137         for (int i=0;i<NUMNODESP1;i++){
    138                 B[NDOF2*NUMNODESP1*0+NDOF2*i+0]=dbasis[0][i];
    139                 B[NDOF2*NUMNODESP1*0+NDOF2*i+1]=0.;
    140                 B[NDOF2*NUMNODESP1*1+NDOF2*i+0]=0.;
    141                 B[NDOF2*NUMNODESP1*1+NDOF2*i+1]=dbasis[1][i];
    142                 B[NDOF2*NUMNODESP1*2+NDOF2*i+0]=0.5*dbasis[1][i];
    143                 B[NDOF2*NUMNODESP1*2+NDOF2*i+1]=0.5*dbasis[0][i];
    144         }
     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];
     145        }
     146
     147        /*Clean-up*/
     148        xDelete<IssmDouble>(dbasis);
    145149}
    146150/*}}}*/
    … …  
    154158         */
    155159
    156         IssmDouble basis[NUMNODESP1];
    157         GetNodalFunctions(&basis[0],gauss);
    158 
     160        /*Fetch number of nodes for this finite element*/
     161        int numnodes = this->NumberofNodes();
     162
     163        /*Get nodal functions*/
     164        IssmDouble* basis=xNew<IssmDouble>(numnodes);
     165        GetNodalFunctions(basis,gauss);
     166
     167        /*Build B for this segment*/
    159168        B[0] = +basis[index1];
    160169        B[1] = +basis[index2];
    161170        B[2] = -basis[index1];
    162171        B[3] = -basis[index2];
     172
     173        /*Clean-up*/
     174        xDelete<IssmDouble>(basis);
    163175}
    164176/*}}}*/
    … …  
    172184         */
    173185
    174         IssmDouble basis[NUMNODESP1];
    175         GetNodalFunctions(&basis[0],gauss);
    176 
     186        /*Fetch number of nodes for this finite element*/
     187        int numnodes = this->NumberofNodes();
     188
     189        /*Get nodal functions*/
     190        IssmDouble* basis=xNew<IssmDouble>(numnodes);
     191        GetNodalFunctions(basis,gauss);
     192
     193        /*Build B'*/
    177194        Bprime[0] = basis[index1];
    178195        Bprime[1] = basis[index2];
    179196        Bprime[2] = basis[index1];
    180197        Bprime[3] = basis[index2];
     198
     199        /*Clean-up*/
     200        xDelete<IssmDouble>(basis);
    181201}
    182202/*}}}*/
    … …  
    186206         * For node i, Bi can be expressed in the actual coordinate system
    187207         * by:
    188          *       Bi=[ h ]
    189          *          [ h ]
    190          * where h is the interpolation function for node i.
    191          *
    192          * We assume B_prog has been allocated already, of size: 2x(NDOF1*NUMNODESP1)
    193          */
    194 
    195         IssmDouble basis[NUMNODESP1];
    196 
    197         /*Get dh1dh2dh3 in actual coordinate system: */
    198         GetNodalFunctions(&basis[0],gauss);
     208         *       Bi=[ N ]
     209         *          [ N ]
     210         * where N is the interpolation function for node i.
     211         *
     212         * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
     213         */
     214
     215        /*Fetch number of nodes for this finite element*/
     216        int numnodes = this->NumberofNodes();
     217
     218        /*Get nodal functions*/
     219        IssmDouble* basis=xNew<IssmDouble>(numnodes);
     220        GetNodalFunctions(basis,gauss);
    199221
    200222        /*Build B: */
    201         for (int i=0;i<NUMNODESP1;i++){
    202                 B[NUMNODESP1*0+i]=basis[i];
    203                 B[NUMNODESP1*1+i]=basis[i];
    204         }
     223        for (int i=0;i<numnodes;i++){
     224                B[numnodes*0+i]=basis[i];
     225                B[numnodes*1+i]=basis[i];
     226        }
     227
     228        /*Clean-up*/
     229        xDelete<IssmDouble>(basis);
    205230}
    206231/*}}}*/
    … …  
    227252
    228253        /*Build B': */
    229         for (int i=0;i<NUMNODESP1;i++){
    230                 Bprime[NDOF2*NUMNODESP1*0+NDOF2*i+0] = 2.*dbasis[0*numnodes+i];
    231                 Bprime[NDOF2*NUMNODESP1*0+NDOF2*i+1] =    dbasis[1*numnodes+i];
    232                 Bprime[NDOF2*NUMNODESP1*1+NDOF2*i+0] =    dbasis[0*numnodes+i];
    233                 Bprime[NDOF2*NUMNODESP1*1+NDOF2*i+1] = 2.*dbasis[1*numnodes+i];
    234                 Bprime[NDOF2*NUMNODESP1*2+NDOF2*i+0] =    dbasis[1*numnodes+i];
    235                 Bprime[NDOF2*NUMNODESP1*2+NDOF2*i+1] =    dbasis[0*numnodes+i];
     254        for (int i=0;i<numnodes;i++){
     255                Bprime[NDOF2*numnodes*0+NDOF2*i+0] = 2.*dbasis[0*numnodes+i];
     256                Bprime[NDOF2*numnodes*0+NDOF2*i+1] =    dbasis[1*numnodes+i];
     257                Bprime[NDOF2*numnodes*1+NDOF2*i+0] =    dbasis[0*numnodes+i];
     258                Bprime[NDOF2*numnodes*1+NDOF2*i+1] = 2.*dbasis[1*numnodes+i];
     259                Bprime[NDOF2*numnodes*2+NDOF2*i+0] =    dbasis[1*numnodes+i];
     260                Bprime[NDOF2*numnodes*2+NDOF2*i+1] =    dbasis[0*numnodes+i];
    236261        }
    237262
    … …  
    242267/*FUNCTION TriaRef::GetBprimeMacAyealStokes {{{*/
    243268void TriaRef::GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss){
    244 
    245269        /*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2.
    246270         * For node i, Bprimei can be expressed in the actual coordinate system
    247271         * by:
    248          *       Bprimei=[  dh/dx    0   ]
    249          *               [    0    dh/dy ]
    250          *               [  dh/dy  dh/dx ]
    251          *               [  dh/dx  dh/dy ]
    252          * where h is the interpolation function for node i.
    253          *
    254          * We assume Bprime has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
    255          */
    256 
    257         /*Same thing in the actual coordinate system: */
    258         IssmDouble dbasis[NDOF2][NUMNODESP1];
    259 
    260         /*Get dh1dh2dh3 in actual coordinates system : */
    261         GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
     272         *       Bprimei=[  dN/dx    0   ]
     273         *               [    0    dN/dy ]
     274         *               [  dN/dy  dN/dx ]
     275         N               [  dN/dx  dN/dy ]
     276         * where N is the interpolation function for node i.
     277         *
     278         * We assume Bprime has been allocated already, of size: 3x(NDOF2*numnodes)
     279         */
     280
     281        /*Fetch number of nodes for this finite element*/
     282        int numnodes = this->NumberofNodes();
     283
     284        /*Get nodal functions*/
     285        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     286        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    262287
    263288        /*Build Bprime: */
    264         for(int i=0;i<NUMNODESP1;i++){
    265                 Bprime[NDOF2*NUMNODESP1*0+NDOF2*i+0]=dbasis[0][i];
    266                 Bprime[NDOF2*NUMNODESP1*0+NDOF2*i+1]=0.;
    267                 Bprime[NDOF2*NUMNODESP1*1+NDOF2*i+0]=0.;
    268                 Bprime[NDOF2*NUMNODESP1*1+NDOF2*i+1]=dbasis[1][i];
    269                 Bprime[NDOF2*NUMNODESP1*2+NDOF2*i+0]=dbasis[1][i];
    270                 Bprime[NDOF2*NUMNODESP1*2+NDOF2*i+1]=dbasis[0][i];
    271                 Bprime[NDOF2*NUMNODESP1*3+NDOF2*i+0]=dbasis[0][i];
    272                 Bprime[NDOF2*NUMNODESP1*3+NDOF2*i+1]=dbasis[1][i];
    273         }
     289        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];
     298        }
     299
     300        /*Clean-up*/
     301        xDelete<IssmDouble>(dbasis);
    274302}
    275303/*}}}*/
    … …  
    279307         * For node i, Bi' can be expressed in the actual coordinate system
    280308         * by:
    281          *       Bi_prime=[ dh/dx ]
    282          *                [ dh/dy ]
    283          * where h is the interpolation function for node i.
    284          *
    285          * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
    286          */
    287 
    288         /*Same thing in the actual coordinate system: */
    289         IssmDouble dbasis[NDOF2][NUMNODESP1];
    290 
    291         /*Get dh1dh2dh3 in actual coordinates system : */
    292         GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
     309         *       Bi_prime=[ dN/dx ]
     310         *                [ dN/dy ]
     311         * where N is the interpolation function for node i.
     312         *
     313         * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
     314         */
     315
     316        /*Fetch number of nodes for this finite element*/
     317        int numnodes = this->NumberofNodes();
     318
     319        /*Get nodal functions derivatives*/
     320        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     321        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    293322
    294323        /*Build B': */
    295         for(int i=0;i<NUMNODESP1;i++){
    296                 Bprime[NUMNODESP1*0+i]=dbasis[0][i];
    297                 Bprime[NUMNODESP1*1+i]=dbasis[1][i];
    298         }
     324        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];
     327        }
     328
     329        /*Clean-up*/
     330        xDelete<IssmDouble>(dbasis);
    299331}
    300332/*}}}*/
    … …  
    308340         * where N is the interpolation function for node i.
    309341         *
    310          * We assume B has been allocated already, of size: 2 x (numdof*NUMNODESP1)
    311          */
    312 
    313         int i;
    314         IssmDouble basis[3];
    315 
    316         /*Get basis in actual coordinate system: */
     342         * We assume B has been allocated already, of size: 2 x (numdof*numnodes)
     343         */
     344
     345        /*Fetch number of nodes for this finite element*/
     346        int numnodes = this->NumberofNodes();
     347
     348        /*Get nodal functions derivatives*/
     349        IssmDouble* basis=xNew<IssmDouble>(numnodes);
    317350        GetNodalFunctions(basis,gauss);
    318351
    319352        /*Build L: */
    320         for (i=0;i<NUMNODESP1;i++){
    321                 B[2*NUMNODESP1*0+2*i+0]=basis[i]; //L[0][NDOF2*i]=dbasis[0][i];
    322                 B[2*NUMNODESP1*0+2*i+1]=0.;
    323                 B[2*NUMNODESP1*1+2*i+0]=0.;
    324                 B[2*NUMNODESP1*1+2*i+1]=basis[i];
    325         }
     353        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];
     358        }
     359
     360        /*Clean-up*/
     361        xDelete<IssmDouble>(basis);
    326362}
    327363/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.