Changeset 15314


Ignore:
Timestamp:
06/22/13 09:20:42 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: cosmetics

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r15313 r15314  
    6363         */
    6464
    65         int i;
    66         IssmDouble dbasis[NDOF2][NUMNODESP1];
    67 
    68         /*Get dh1dh2dh3 in actual coordinate system: */
    69         GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
     65        /*Fetch number of nodes for this finite element*/
     66        int numnodes = this->NumberofNodes();
     67
     68        /*Get nodal functions*/
     69        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     70        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    7071
    7172        /*Build B: */
    72         for (i=0;i<NUMNODESP1;i++){
    73                 B[NDOF1*NUMNODESP1*0+NDOF1*i]=dbasis[0][i];
    74                 B[NDOF1*NUMNODESP1*1+NDOF1*i]=dbasis[1][i];
    75         }
     73        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];
     76        }
     77
     78        /*Clean-up*/
     79        xDelete<IssmDouble>(dbasis);
    7680}
    7781/*}}}*/
    … …  
    132136        /*Build B': */
    133137        for (int i=0;i<NUMNODESP1;i++){
    134                 *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i];
    135                 *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0;
    136                 *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0;
    137                 *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
    138                 *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=0.5*dbasis[1][i];
    139                 *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=0.5*dbasis[0][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];
    140144        }
    141145}
    … …  
    150154         */
    151155
    152         IssmDouble l1l3[NUMNODESP1];
    153 
    154         GetNodalFunctions(&l1l3[0],gauss);
    155 
    156         B[0] = +l1l3[index1];
    157         B[1] = +l1l3[index2];
    158         B[2] = -l1l3[index1];
    159         B[3] = -l1l3[index2];
     156        IssmDouble basis[NUMNODESP1];
     157        GetNodalFunctions(&basis[0],gauss);
     158
     159        B[0] = +basis[index1];
     160        B[1] = +basis[index2];
     161        B[2] = -basis[index1];
     162        B[3] = -basis[index2];
    160163}
    161164/*}}}*/
    … …  
    169172         */
    170173
    171         IssmDouble l1l3[NUMNODESP1];
    172 
    173         GetNodalFunctions(&l1l3[0],gauss);
    174 
    175         Bprime[0] = l1l3[index1];
    176         Bprime[1] = l1l3[index2];
    177         Bprime[2] = l1l3[index1];
    178         Bprime[3] = l1l3[index2];
     174        IssmDouble basis[NUMNODESP1];
     175        GetNodalFunctions(&basis[0],gauss);
     176
     177        Bprime[0] = basis[index1];
     178        Bprime[1] = basis[index2];
     179        Bprime[2] = basis[index1];
     180        Bprime[3] = basis[index2];
    179181}
    180182/*}}}*/
    181183/*FUNCTION TriaRef::GetBPrognostic{{{*/
    182 void TriaRef::GetBPrognostic(IssmDouble* B_prog, IssmDouble* xyz_list, GaussTria* gauss){
     184void TriaRef::GetBPrognostic(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
    183185        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    184186         * For node i, Bi can be expressed in the actual coordinate system
    … …  
    196198        GetNodalFunctions(&basis[0],gauss);
    197199
    198         /*Build B_prog: */
     200        /*Build B: */
    199201        for (int i=0;i<NUMNODESP1;i++){
    200                 *(B_prog+NDOF1*NUMNODESP1*0+NDOF1*i)=basis[i];
    201                 *(B_prog+NDOF1*NUMNODESP1*1+NDOF1*i)=basis[i];
     202                B[NUMNODESP1*0+i]=basis[i];
     203                B[NUMNODESP1*1+i]=basis[i];
    202204        }
    203205}
    … …  
    220222        int numnodes = this->NumberofNodes();
    221223
    222         /*Get nodal functions*/
     224        /*Get nodal functions derivatives*/
    223225        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
    224226        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    … …  
    260262
    261263        /*Build Bprime: */
    262         for (int i=0;i<NUMNODESP1;i++){
    263                 *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i];
    264                 *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0;
    265                 *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i)=0;
    266                 *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
    267                 *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i)=dbasis[1][i];
    268                 *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dbasis[0][i];
    269                 *(Bprime+NDOF2*NUMNODESP1*3+NDOF2*i)=dbasis[0][i];
    270                 *(Bprime+NDOF2*NUMNODESP1*3+NDOF2*i+1)=dbasis[1][i];
     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];
    271273        }
    272274}
    273275/*}}}*/
    274276/*FUNCTION TriaRef::GetBprimePrognostic{{{*/
    275 void TriaRef::GetBprimePrognostic(IssmDouble* Bprime_prog, IssmDouble* xyz_list, GaussTria* gauss){
     277void TriaRef::GetBprimePrognostic(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss){
    276278        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    277279         * For node i, Bi' can be expressed in the actual coordinate system
    … …  
    291293
    292294        /*Build B': */
    293         for (int i=0;i<NUMNODESP1;i++){
    294                 *(Bprime_prog+NDOF1*NUMNODESP1*0+NDOF1*i)=dbasis[0][i];
    295                 *(Bprime_prog+NDOF1*NUMNODESP1*1+NDOF1*i)=dbasis[1][i];
    296         }
    297 }
    298 /*}}}*/
    299 /*FUNCTION TriaRef::GetL{{{*/
     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        }
     299}
     300/*}}}*/
     301/*FUNCTION TriaRef::GetBMacAyealFriction{{{*/
    300302void TriaRef::GetBMacAyealFriction(IssmDouble* B, IssmDouble* xyz_list,GaussTria* gauss){
    301303        /*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2.
    … …  
    317319        /*Build L: */
    318320        for (i=0;i<NUMNODESP1;i++){
    319                 *(B+2*NUMNODESP1*0+2*i)=basis[i]; //L[0][NDOF2*i]=dbasis[0][i];
    320                 *(B+2*NUMNODESP1*0+2*i+1)=0;
    321                 *(B+2*NUMNODESP1*1+2*i)=0;
    322                 *(B+2*NUMNODESP1*1+2*i+1)=basis[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];
    323325        }
    324326}
    … …  
    328330        /*The Jacobian is constant over the element, discard the gaussian points.
    329331         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    330         IssmDouble x1,y1,x2,y2,x3,y3;
    331 
    332         x1=*(xyz_list+NUMNODESP1*0+0);
    333         y1=*(xyz_list+NUMNODESP1*0+1);
    334         x2=*(xyz_list+NUMNODESP1*1+0);
    335         y2=*(xyz_list+NUMNODESP1*1+1);
    336         x3=*(xyz_list+NUMNODESP1*2+0);
    337         y3=*(xyz_list+NUMNODESP1*2+1);
    338 
    339         *(J+NDOF2*0+0)=0.5*(x2-x1);
    340         *(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);
    341         *(J+NDOF2*0+1)=0.5*(y2-y1);
    342         *(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
     332
     333        IssmDouble x1 = xyz_list[3*0+0];
     334        IssmDouble y1 = xyz_list[3*0+1];
     335        IssmDouble x2 = xyz_list[3*1+0];
     336        IssmDouble y2 = xyz_list[3*1+1];
     337        IssmDouble x3 = xyz_list[3*2+0];
     338        IssmDouble y3 = xyz_list[3*2+1];
     339
     340        J[2*0+0] = 0.5*(x2-x1);
     341        J[2*1+0] = SQRT3/6.0*(2*x3-x1-x2);
     342        J[2*0+1] = 0.5*(y2-y1);
     343        J[2*1+1] = SQRT3/6.0*(2*y3-y1-y2);
    343344}
    344345/*}}}*/
    … …  
    347348        /*The Jacobian determinant is constant over the element, discard the gaussian points.
    348349         * J is assumed to have been allocated*/
    349         IssmDouble x1,y1,x2,y2;
    350 
    351         x1=*(xyz_list+3*0+0);
    352         y1=*(xyz_list+3*0+1);
    353         x2=*(xyz_list+3*1+0);
    354         y2=*(xyz_list+3*1+1);
    355 
    356         *Jdet=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.));
     350
     351        IssmDouble x1 = xyz_list[3*0+0];
     352        IssmDouble y1 = xyz_list[3*0+1];
     353        IssmDouble x2 = xyz_list[3*1+0];
     354        IssmDouble y2 = xyz_list[3*1+1];
     355
     356        *Jdet = .5*sqrt(pow(x2-x1,2) + pow(y2-y1,2));
    357357        if(*Jdet<0) _error_("negative jacobian determinant!");
    358358
    … …  
    391391void TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss){
    392392        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     393
     394        _assert_(basis);
    393395
    394396        switch(this->element_type){
    … …  
    483485        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    484486         * natural coordinate system) at the gaussian point. */
     487
     488        _assert_(dbasis && gauss);
    485489
    486490        switch(this->element_type){
Note: See TracChangeset for help on using the changeset viewer.