Changeset 14767


Ignore:
Timestamp:
04/26/13 10:07:43 (12 years ago)
Author:
Mathieu Morlighem
Message:

NEW: added some functions for P2 Finite elements

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

Legend:

Unmodified
Added
Removed

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

    r14591 r14767  
    2121
    2222/*Element macros*/
    23 #define NUMNODES 3
     23#define NUMNODESP1 3
     24#define NUMNODESP2 6
    2425
    2526/*Object constructors and destructor*/
    … …  
    6566         * where h is the interpolation function for node i.
    6667         *
    67          * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES)
     68         * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
    6869         */
    6970
    7071        int i;
    71         IssmDouble dbasis[NDOF2][NUMNODES];
     72        IssmDouble dbasis[NDOF2][NUMNODESP1];
    7273
    7374        /*Get dh1dh2dh3 in actual coordinate system: */
    … …  
    7576
    7677        /*Build B: */
    77         for (i=0;i<NUMNODES;i++){
    78                 B[NDOF1*NUMNODES*0+NDOF1*i]=dbasis[0][i];
    79                 B[NDOF1*NUMNODES*1+NDOF1*i]=dbasis[1][i];
     78        for (i=0;i<NUMNODESP1;i++){
     79                B[NDOF1*NUMNODESP1*0+NDOF1*i]=dbasis[0][i];
     80                B[NDOF1*NUMNODESP1*1+NDOF1*i]=dbasis[1][i];
    8081        }
    8182}
    … …  
    9192         * where h is the interpolation function for node i.
    9293         *
    93          * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES)
     94         * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
    9495         */
    9596
    9697        int i;
    97         IssmDouble dbasis[NDOF2][NUMNODES];
     98        IssmDouble dbasis[NDOF2][NUMNODESP1];
    9899
    99100        /*Get dh1dh2dh3 in actual coordinate system: */
    … …  
    101102
    102103        /*Build B: */
    103         for (i=0;i<NUMNODES;i++){
    104                 *(B+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i]; //B[0][NDOF2*i]=dbasis[0][i];
    105                 *(B+NDOF2*NUMNODES*0+NDOF2*i+1)=0;
    106                 *(B+NDOF2*NUMNODES*1+NDOF2*i)=0;
    107                 *(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i];
    108                 *(B+NDOF2*NUMNODES*2+NDOF2*i)=(float).5*dbasis[1][i];
    109                 *(B+NDOF2*NUMNODES*2+NDOF2*i+1)=(float).5*dbasis[0][i];
     104        for (i=0;i<NUMNODESP1;i++){
     105                *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i]; //B[0][NDOF2*i]=dbasis[0][i];
     106                *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0;
     107                *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0;
     108                *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
     109                *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dbasis[1][i];
     110                *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][i];
    110111        }
    111112}
    … …  
    122123         * where h is the interpolation function for node i.
    123124         *
    124          * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES)
     125         * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
    125126         */
    126127
    127128        /*Same thing in the actual coordinate system: */
    128         IssmDouble dbasis[NDOF2][NUMNODES];
     129        IssmDouble dbasis[NDOF2][NUMNODESP1];
    129130
    130131        /*Get dh1dh2dh3 in actual coordinates system : */
    … …  
    132133
    133134        /*Build B': */
    134         for (int i=0;i<NUMNODES;i++){
    135                 *(B+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i];
    136                 *(B+NDOF2*NUMNODES*0+NDOF2*i+1)=0;
    137                 *(B+NDOF2*NUMNODES*1+NDOF2*i)=0;
    138                 *(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i];
    139                 *(B+NDOF2*NUMNODES*2+NDOF2*i)=0.5*dbasis[1][i];
    140                 *(B+NDOF2*NUMNODES*2+NDOF2*i+1)=0.5*dbasis[0][i];
     135        for (int i=0;i<NUMNODESP1;i++){
     136                *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i];
     137                *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0;
     138                *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0;
     139                *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
     140                *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=0.5*dbasis[1][i];
     141                *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=0.5*dbasis[0][i];
    141142        }
    142143}
    … …  
    151152         */
    152153
    153         IssmDouble l1l3[NUMNODES];
     154        IssmDouble l1l3[NUMNODESP1];
    154155
    155156        GetNodalFunctions(&l1l3[0],gauss);
    … …  
    170171         */
    171172
    172         IssmDouble l1l3[NUMNODES];
     173        IssmDouble l1l3[NUMNODESP1];
    173174
    174175        GetNodalFunctions(&l1l3[0],gauss);
    … …  
    189190         * where h is the interpolation function for node i.
    190191         *
    191          * We assume B_prog has been allocated already, of size: 2x(NDOF1*NUMNODES)
    192          */
    193 
    194         IssmDouble basis[NUMNODES];
     192         * We assume B_prog has been allocated already, of size: 2x(NDOF1*NUMNODESP1)
     193         */
     194
     195        IssmDouble basis[NUMNODESP1];
    195196
    196197        /*Get dh1dh2dh3 in actual coordinate system: */
    … …  
    198199
    199200        /*Build B_prog: */
    200         for (int i=0;i<NUMNODES;i++){
    201                 *(B_prog+NDOF1*NUMNODES*0+NDOF1*i)=basis[i];
    202                 *(B_prog+NDOF1*NUMNODES*1+NDOF1*i)=basis[i];
     201        for (int i=0;i<NUMNODESP1;i++){
     202                *(B_prog+NDOF1*NUMNODESP1*0+NDOF1*i)=basis[i];
     203                *(B_prog+NDOF1*NUMNODESP1*1+NDOF1*i)=basis[i];
    203204        }
    204205}
    … …  
    215216         * where h is the interpolation function for node i.
    216217         *
    217          * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODES)
     218         * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
    218219         */
    219220
    220221        /*Same thing in the actual coordinate system: */
    221         IssmDouble dbasis[NDOF2][NUMNODES];
     222        IssmDouble dbasis[NDOF2][NUMNODESP1];
    222223
    223224        /*Get dh1dh2dh3 in actual coordinates system : */
    … …  
    225226
    226227        /*Build B': */
    227         for (int i=0;i<NUMNODES;i++){
    228                 *(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=2*dbasis[0][i];
    229                 *(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=dbasis[1][i];
    230                 *(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=dbasis[0][i];
    231                 *(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=2*dbasis[1][i];
    232                 *(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dbasis[1][i];
    233                 *(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dbasis[0][i];
     228        for (int i=0;i<NUMNODESP1;i++){
     229                *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i)=2*dbasis[0][i];
     230                *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i+1)=dbasis[1][i];
     231                *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i)=dbasis[0][i];
     232                *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2*dbasis[1][i];
     233                *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i)=dbasis[1][i];
     234                *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dbasis[0][i];
    234235        }
    235236}
    … …  
    247248         * where h is the interpolation function for node i.
    248249         *
    249          * We assume Bprime has been allocated already, of size: 3x(NDOF2*NUMNODES)
     250         * We assume Bprime has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
    250251         */
    251252
    252253        /*Same thing in the actual coordinate system: */
    253         IssmDouble dbasis[NDOF2][NUMNODES];
     254        IssmDouble dbasis[NDOF2][NUMNODESP1];
    254255
    255256        /*Get dh1dh2dh3 in actual coordinates system : */
    … …  
    257258
    258259        /*Build Bprime: */
    259         for (int i=0;i<NUMNODES;i++){
    260                 *(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i];
    261                 *(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=0;
    262                 *(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=0;
    263                 *(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i];
    264                 *(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dbasis[1][i];
    265                 *(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dbasis[0][i];
    266                 *(Bprime+NDOF2*NUMNODES*3+NDOF2*i)=dbasis[0][i];
    267                 *(Bprime+NDOF2*NUMNODES*3+NDOF2*i+1)=dbasis[1][i];
     260        for (int i=0;i<NUMNODESP1;i++){
     261                *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i];
     262                *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0;
     263                *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i)=0;
     264                *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
     265                *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i)=dbasis[1][i];
     266                *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dbasis[0][i];
     267                *(Bprime+NDOF2*NUMNODESP1*3+NDOF2*i)=dbasis[0][i];
     268                *(Bprime+NDOF2*NUMNODESP1*3+NDOF2*i+1)=dbasis[1][i];
    268269        }
    269270}
    … …  
    278279         * where h is the interpolation function for node i.
    279280         *
    280          * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODES)
     281         * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
    281282         */
    282283
    283284        /*Same thing in the actual coordinate system: */
    284         IssmDouble dbasis[NDOF2][NUMNODES];
     285        IssmDouble dbasis[NDOF2][NUMNODESP1];
    285286
    286287        /*Get dh1dh2dh3 in actual coordinates system : */
    … …  
    288289
    289290        /*Build B': */
    290         for (int i=0;i<NUMNODES;i++){
    291                 *(Bprime_prog+NDOF1*NUMNODES*0+NDOF1*i)=dbasis[0][i];
    292                 *(Bprime_prog+NDOF1*NUMNODES*1+NDOF1*i)=dbasis[1][i];
     291        for (int i=0;i<NUMNODESP1;i++){
     292                *(Bprime_prog+NDOF1*NUMNODESP1*0+NDOF1*i)=dbasis[0][i];
     293                *(Bprime_prog+NDOF1*NUMNODESP1*1+NDOF1*i)=dbasis[1][i];
    293294        }
    294295}
    … …  
    306307         * where h is the interpolation function for node i.
    307308         *
    308          * We assume L has been allocated already, of size: NUMNODES (numdof=1), or numdofx(numdof*NUMNODES) (numdof=2)
     309         * We assume L has been allocated already, of size: NUMNODESP1 (numdof=1), or numdofx(numdof*NUMNODESP1) (numdof=2)
    309310         */
    310311
    … …  
    317318        /*Build L: */
    318319        if(numdof==1){
    319                 for (i=0;i<NUMNODES;i++){
     320                for (i=0;i<NUMNODESP1;i++){
    320321                        L[i]=basis[i];
    321322                }
    322323        }
    323324        else{
    324                 for (i=0;i<NUMNODES;i++){
    325                         *(L+numdof*NUMNODES*0+numdof*i)=basis[i]; //L[0][NDOF2*i]=dbasis[0][i];
    326                         *(L+numdof*NUMNODES*0+numdof*i+1)=0;
    327                         *(L+numdof*NUMNODES*1+numdof*i)=0;
    328                         *(L+numdof*NUMNODES*1+numdof*i+1)=basis[i];
     325                for (i=0;i<NUMNODESP1;i++){
     326                        *(L+numdof*NUMNODESP1*0+numdof*i)=basis[i]; //L[0][NDOF2*i]=dbasis[0][i];
     327                        *(L+numdof*NUMNODESP1*0+numdof*i+1)=0;
     328                        *(L+numdof*NUMNODESP1*1+numdof*i)=0;
     329                        *(L+numdof*NUMNODESP1*1+numdof*i+1)=basis[i];
    329330                }
    330331        }
    … …  
    337338        IssmDouble x1,y1,x2,y2,x3,y3;
    338339
    339         x1=*(xyz_list+NUMNODES*0+0);
    340         y1=*(xyz_list+NUMNODES*0+1);
    341         x2=*(xyz_list+NUMNODES*1+0);
    342         y2=*(xyz_list+NUMNODES*1+1);
    343         x3=*(xyz_list+NUMNODES*2+0);
    344         y3=*(xyz_list+NUMNODES*2+1);
     340        x1=*(xyz_list+NUMNODESP1*0+0);
     341        y1=*(xyz_list+NUMNODESP1*0+1);
     342        x2=*(xyz_list+NUMNODESP1*1+0);
     343        y2=*(xyz_list+NUMNODESP1*1+1);
     344        x3=*(xyz_list+NUMNODESP1*2+0);
     345        y3=*(xyz_list+NUMNODESP1*2+1);
    345346
    346347        *(J+NDOF2*0+0)=0.5*(x2-x1);
    … …  
    405406}
    406407/*}}}*/
     408/*FUNCTION TriaRef::GetNodalFunctionsP2{{{*/
     409void TriaRef::GetNodalFunctionsP2(IssmDouble** pbasis,GaussTria* gauss){
     410        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     411
     412        /*Allocate*/
     413        IssmDouble* basis =xNew<IssmDouble>(NUMNODESP2);
     414
     415        basis[0]=gauss->coord1*(2.*gauss->coord1-1.);
     416        basis[1]=gauss->coord2*(2.*gauss->coord2-1.);
     417        basis[2]=gauss->coord3*(2.*gauss->coord3-1.);
     418        basis[3]=4.*gauss->coord3*gauss->coord2;
     419        basis[4]=4.*gauss->coord3*gauss->coord1;
     420        basis[5]=4.*gauss->coord1*gauss->coord2;
     421
     422        /*Assign output pointer*/
     423        *pbasis = basis;
     424
     425}
     426/*}}}*/
    407427/*FUNCTION TriaRef::GetSegmentNodalFunctions{{{*/
    408428void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss,int index1,int index2){
    … …  
    425445         * actual coordinate system): */
    426446        int       i;
    427         IssmDouble    dbasis_ref[NDOF2][NUMNODES];
     447        IssmDouble    dbasis_ref[NDOF2][NUMNODESP1];
    428448        IssmDouble    Jinv[NDOF2][NDOF2];
    429449
    … …  
    439459         * [dhi/dy]       [dhi/ds]
    440460         */
    441         for (i=0;i<NUMNODES;i++){
    442                 dbasis[NUMNODES*0+i]=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i];
    443                 dbasis[NUMNODES*1+i]=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i];
     461        for (i=0;i<NUMNODESP1;i++){
     462                dbasis[NUMNODESP1*0+i]=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i];
     463                dbasis[NUMNODESP1*1+i]=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i];
    444464        }
    445465
    … …  
    452472
    453473        /*First nodal function: */
    454         *(dl1dl3+NUMNODES*0+0)=-0.5;
    455         *(dl1dl3+NUMNODES*1+0)=-1.0/(2.0*SQRT3);
     474        *(dl1dl3+NUMNODESP1*0+0)=-0.5;
     475        *(dl1dl3+NUMNODESP1*1+0)=-1.0/(2.0*SQRT3);
    456476
    457477        /*Second nodal function: */
    458         *(dl1dl3+NUMNODES*0+1)=0.5;
    459         *(dl1dl3+NUMNODES*1+1)=-1.0/(2.0*SQRT3);
     478        *(dl1dl3+NUMNODESP1*0+1)=0.5;
     479        *(dl1dl3+NUMNODESP1*1+1)=-1.0/(2.0*SQRT3);
    460480
    461481        /*Third nodal function: */
    462         *(dl1dl3+NUMNODES*0+2)=0;
    463         *(dl1dl3+NUMNODES*1+2)=1.0/SQRT3;
    464 
     482        *(dl1dl3+NUMNODESP1*0+2)=0;
     483        *(dl1dl3+NUMNODESP1*1+2)=1.0/SQRT3;
     484
     485}
     486/*}}}*/
     487/*FUNCTION TriaRef::GetNodalFunctionsDerivativesP2Reference{{{*/
     488void TriaRef::GetNodalFunctionsDerivativesP2Reference(IssmDouble** pdbasis,GaussTria* gauss){
     489        /*This routine returns the values of the nodal functions derivatives  (with respect to the
     490         * natural coordinate system) at the gaussian point. */
     491
     492        /*Allocate*/
     493        IssmDouble* dbasis =xNew<IssmDouble>(NUMNODESP2*2);
     494
     495        /*Nodal function 1*/
     496        dbasis[NUMNODESP2*0+0]=-2.*gauss->coord1 + 0.5;
     497        dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.;
     498        /*Nodal function 2*/
     499        dbasis[NUMNODESP2*0+0]=+2.*gauss->coord2 + 0.5;
     500        dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.;
     501        /*Nodal function 3*/
     502        dbasis[NUMNODESP2*0+0]=0.;
     503        dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.;
     504        /*Nodal function 4*/
     505        dbasis[NUMNODESP2*0+0]=+2.*gauss->coord3;
     506        dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3;
     507        /*Nodal function 5*/
     508        dbasis[NUMNODESP2*0+0]=-2.*gauss->coord3;
     509        dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3;
     510        /*Nodal function 6*/
     511        dbasis[NUMNODESP2*0+0]=2.*(gauss->coord1-gauss->coord2);
     512        dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*(gauss->coord1+gauss->coord2);
     513
     514        /*Assign output pointer*/
     515        *pdbasis = dbasis;
    465516}
    466517/*}}}*/

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

    r14591 r14767  
    3535                void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss);
    3636                void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,GaussTria* gauss);
    37                 void GetNodalFunctions(IssmDouble* l1l2l3,GaussTria* gauss);
    38                 void GetSegmentNodalFunctions(IssmDouble* l1l2l3,GaussTria* gauss, int index1,int index2);
     37                void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss);
     38                void GetNodalFunctionsP2(IssmDouble** pbasis,GaussTria* gauss);
     39                void GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss, int index1,int index2);
    3940                void GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2);
    4041                void GetSegmentBprimeFlux(IssmDouble* Bprime,GaussTria* gauss, int index1,int index2);
    41                 void GetNodalFunctionsDerivatives(IssmDouble* l1l2l3,IssmDouble* xyz_list, GaussTria* gauss);
     42                void GetNodalFunctionsDerivatives(IssmDouble* basis,IssmDouble* xyz_list, GaussTria* gauss);
    4243                void GetNodalFunctionsDerivativesReference(IssmDouble* dl1dl3,GaussTria* gauss);
     44                void GetNodalFunctionsDerivativesP2Reference(IssmDouble** pdbasis,GaussTria* gauss);
    4345                void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss);
    4446                void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss);
Note: See TracChangeset for help on using the changeset viewer.