Changeset 5745


Ignore:
Timestamp:
09/10/10 10:53:39 (15 years ago)
Author:
Mathieu Morlighem
Message:

Some cleaning

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

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/Element.h

    r5743 r5745  
    3232                virtual void   GetSolutionFromInputs(Vec solution)=0;
    3333                virtual int    GetNodeIndex(Node* node)=0;
    34                 virtual bool   GetShelf()=0;
    35                 virtual bool   GetOnBed()=0;
     34                virtual bool   IsOnShelf()=0;
     35                virtual bool   IsOnBed()=0;
    3636                virtual void   GetParameterListOnVertices(double* pvalue,int enumtype)=0;
    3737                virtual void   GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue)=0;

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

    r5743 r5745  
    897897        ISSMERROR("Node provided not found among element nodes");
    898898
    899 }
    900 /*}}}*/
    901 /*FUNCTION Penta::GetOnBed {{{1*/
    902 bool Penta::GetOnBed(){
    903        
    904         bool onbed;
    905 
    906         inputs->GetParameterValue(&onbed,ElementOnBedEnum);
    907 
    908         return onbed;
    909 }
    910 /*}}}*/
    911 /*FUNCTION Penta::GetShelf {{{1*/
    912 bool   Penta::GetShelf(){
    913         bool onshelf;
    914 
    915         /*retrieve inputs :*/
    916         inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum);
    917 
    918         return onshelf;
    919899}
    920900/*}}}*/
    … …  
    20011981
    20021982        /*If shelf: hydrostatic equilibrium*/
    2003         if (this->GetShelf()){
     1983        if (this->IsOnShelf()){
    20041984
    20051985                /*recover material parameters: */
    … …  
    54835463}
    54845464/*}}}*/
     5465/*FUNCTION Penta::IsOnShelf {{{1*/
     5466bool   Penta::IsOnShelf(){
     5467
     5468        bool onshelf;
     5469        inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum);
     5470        return onshelf;
     5471
     5472}
     5473/*}}}*/
    54855474/*FUNCTION Penta::IsOnSurface{{{1*/
    54865475bool Penta::IsOnSurface(void){

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

    r5743 r5745  
    7979                void   DeleteResults(void);
    8080                int    GetNodeIndex(Node* node);
    81                 bool   GetOnBed();
    82                 bool   GetShelf();
    8381                void   GetSolutionFromInputs(Vec solution);
    8482                void   GetVectorFromInputs(Vec vector,int NameEnum);
    … …  
    179177                bool      IsOnSurface(void);
    180178                bool      IsOnBed(void);
     179                bool    IsOnShelf(void);
    181180                void      ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp);
    182181                void      ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp);

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

    r5743 r5745  
    1919#include "../../include/include.h"
    2020/*}}}*/
     21
     22/*Element macros*/
     23#define NUMNODESP1    6
     24#define NUMNODESP1_2d 3
     25#define NUMNODESMINI  7
     26#define NDOF1 1
     27#define NDOF2 2
     28#define NDOF3 3
     29#define NDOF4 4
    2130
    2231/*Object constructors and destructor*/
    … …  
    6271         * where h is the interpolation function for grid i.
    6372         *
    64          * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
    65          */
    66 
    67         int i;
    68         const int numgrids=6;
    69         const int NDOF3=3;
    70         const int NDOF2=2;
    71 
    72         double dh1dh6[NDOF3][numgrids];
     73         * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
     74         */
     75
     76        double dh1dh6[3][NUMNODESP1];
    7377
    7478        /*Get dh1dh6 in actual coordinate system: */
    … …  
    7680
    7781        /*Build B: */
    78         for (i=0;i<numgrids;i++){
    79                 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i];
    80                 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;
    81 
    82                 *(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
    83                 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];
    84 
    85                 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i];
    86                 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];
     82        for (int i=0;i<NUMNODESP1;i++){
     83                *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dh1dh6[0][i];
     84                *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.0;
     85
     86                *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0.0;
     87                *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dh1dh6[1][i];
     88
     89                *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dh1dh6[1][i];
     90                *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];
    8791
    8892        }
    … …  
    102106         * where h is the interpolation function for grid i.
    103107         *
    104          * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
    105          */
    106 
    107         int i;
    108         const int numgrids=6;
    109         const int NDOF3=3;
    110         const int NDOF2=2;
    111 
    112         double dh1dh6[NDOF3][numgrids];
     108         * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
     109         */
     110
     111        double dh1dh6[3][NUMNODESP1];
    113112
    114113        /*Get dh1dh6 in actual coordinate system: */
    … …  
    116115
    117116        /*Build B: */
    118         for (i=0;i<numgrids;i++){
    119                 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i];
    120                 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;
    121 
    122                 *(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
    123                 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];
    124 
    125                 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i];
    126                 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];
    127 
    128                 *(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh6[2][i];
    129                 *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
    130 
    131                 *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
    132                 *(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh6[2][i];
     117        for (int i=0;i<NUMNODESP1;i++){
     118                *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dh1dh6[0][i];
     119                *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.0;
     120
     121                *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0.0;
     122                *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dh1dh6[1][i];
     123
     124                *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dh1dh6[1][i];
     125                *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];
     126
     127                *(B+NDOF2*NUMNODESP1*3+NDOF2*i)=(float).5*dh1dh6[2][i];
     128                *(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
     129
     130                *(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
     131                *(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=(float).5*dh1dh6[2][i];
    133132        }
    134133
    … …  
    147146         * where h is the interpolation function for grid i.
    148147         *
    149          * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
    150          */
    151 
    152         int i;
    153         const int NDOF3=3;
    154         const int NDOF2=2;
    155         const int numgrids=6;
    156 
    157         double dh1dh6[NDOF3][numgrids];
     148         * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
     149         */
     150        double dh1dh6[3][NUMNODESP1];
    158151
    159152        /*Get dh1dh6 in actual coordinate system: */
    … …  
    161154
    162155        /*Build BPrime: */
    163         for (i=0;i<numgrids;i++){
    164                 *(B+NDOF2*numgrids*0+NDOF2*i)=2.0*dh1dh6[0][i];
    165                 *(B+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh6[1][i];
    166 
    167                 *(B+NDOF2*numgrids*1+NDOF2*i)=dh1dh6[0][i];
    168                 *(B+NDOF2*numgrids*1+NDOF2*i+1)=2.0*dh1dh6[1][i];
    169 
    170                 *(B+NDOF2*numgrids*2+NDOF2*i)=dh1dh6[1][i];
    171                 *(B+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh6[0][i];
    172 
    173                 *(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh6[2][i];
    174                 *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
    175 
    176                 *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
    177                 *(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh6[2][i];
     156        for (int i=0;i<NUMNODESP1;i++){
     157                *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=2.0*dh1dh6[0][i];
     158                *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=dh1dh6[1][i];
     159
     160                *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=dh1dh6[0][i];
     161                *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2.0*dh1dh6[1][i];
     162
     163                *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=dh1dh6[1][i];
     164                *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dh1dh6[0][i];
     165
     166                *(B+NDOF2*NUMNODESP1*3+NDOF2*i)=dh1dh6[2][i];
     167                *(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
     168
     169                *(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
     170                *(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=dh1dh6[2][i];
    178171        }
    179172}
    … …  
    182175void PentaRef::GetBStokes(double* B, double* xyz_list, GaussPenta* gauss){
    183176
    184         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID.
     177        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4.
    185178         * For grid i, Bi can be expressed in the actual coordinate system
    186179         * by:          Bi=[ dh/dx          0              0       0  ]
    … …  
    197190
    198191        int i;
    199         const int calculationdof=3;
    200         const int numgrids=6;
    201         int DOFPERGRID=4;
    202 
    203         double dh1dh7[calculationdof][numgrids+1];
    204         double l1l6[numgrids];
     192
     193        double dh1dh7[3][NUMNODESMINI];
     194        double l1l6[NUMNODESP1];
    205195
    206196
    … …  
    210200
    211201        /*Build B: */
    212         for (i=0;i<numgrids+1;i++){
    213                 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B[0][DOFPERGRID*i]=dh1dh6[0][i];
    214                 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
    215                 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
    216                 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
    217                 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
    218                 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
    219                 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
    220                 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
    221                 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
    222                 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7[1][i];
    223                 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7[0][i];
    224                 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
    225                 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7[2][i];
    226                 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
    227                 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7[0][i];
    228                 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
    229                 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7[2][i];
    230                 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7[1][i];
    231                 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=0;
    232                 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=0;
    233                 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=0;
    234                 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7[0][i];
    235                 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7[1][i];
    236                 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7[2][i];
    237         }
    238 
    239         for (i=0;i<numgrids;i++){ //last column not for the bubble function
    240                 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;
    241                 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;
    242                 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;
    243                 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;
    244                 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;
    245                 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;
    246                 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=l1l6[i];
    247                 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=0;
     202        for (i=0;i<NUMNODESMINI;i++){
     203                *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh7[0][i]; //B[0][NDOF4*i]=dh1dh6[0][i];
     204                *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0;
     205                *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
     206                *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0;
     207                *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh7[1][i];
     208                *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
     209                *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0;
     210                *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0;
     211                *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=dh1dh7[2][i];
     212                *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=(float).5*dh1dh7[1][i];
     213                *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=(float).5*dh1dh7[0][i];
     214                *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;
     215                *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i)=(float).5*dh1dh7[2][i];
     216                *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1)=0;
     217                *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2)=(float).5*dh1dh7[0][i];
     218                *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i)=0;
     219                *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1)=(float).5*dh1dh7[2][i];
     220                *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2)=(float).5*dh1dh7[1][i];
     221                *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i)=0;
     222                *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1)=0;
     223                *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2)=0;
     224                *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i)=dh1dh7[0][i];
     225                *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1)=dh1dh7[1][i];
     226                *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2)=dh1dh7[2][i];
     227        }
     228
     229        for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
     230                *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
     231                *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
     232                *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
     233                *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0;
     234                *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0;
     235                *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0;
     236                *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=l1l6[i];
     237                *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=0;
    248238        }
    249239
    … …  
    269259
    270260        int i;
    271         const int calculationdof=3;
    272         const int numgrids=6;
    273         int DOFPERGRID=4;
    274 
    275         double dh1dh7[calculationdof][numgrids+1];
    276         double l1l6[numgrids];
     261        double dh1dh7[3][NUMNODESMINI];
     262        double l1l6[NUMNODESP1];
    277263
    278264        /*Get dh1dh7 in actual coordinate system: */
    279265        GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
    280 
    281266        GetNodalFunctionsP1(l1l6, gauss);
    282267
    283268        /*B_primeuild B_prime: */
    284         for (i=0;i<numgrids+1;i++){
    285                 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B_prime[0][DOFPERGRID*i]=dh1dh6[0][i];
    286                 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
    287                 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
    288                 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
    289                 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
    290                 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
    291                 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
    292                 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
    293                 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
    294                 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=dh1dh7[1][i];
    295                 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=dh1dh7[0][i];
    296                 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
    297                 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=dh1dh7[2][i];
    298                 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
    299                 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=dh1dh7[0][i];
    300                 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
    301                 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=dh1dh7[2][i];
    302                 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=dh1dh7[1][i];
    303                 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=dh1dh7[0][i];
    304                 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=dh1dh7[1][i];
    305                 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=dh1dh7[2][i];
    306                 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=0;
    307                 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=0;
    308                 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=0;
    309         }
    310 
    311         for (i=0;i<numgrids;i++){ //last column not for the bubble function
    312                 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;
    313                 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;
    314                 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;
    315                 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;
    316                 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;
    317                 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;
    318                 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=0;
    319                 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=l1l6[i];
     269        for (i=0;i<NUMNODESMINI;i++){
     270                *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh7[0][i]; //B_prime[0][NDOF4*i]=dh1dh6[0][i];
     271                *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0;
     272                *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
     273                *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0;
     274                *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh7[1][i];
     275                *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
     276                *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0;
     277                *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0;
     278                *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=dh1dh7[2][i];
     279                *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=dh1dh7[1][i];
     280                *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=dh1dh7[0][i];
     281                *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;
     282                *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i)=dh1dh7[2][i];
     283                *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1)=0;
     284                *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2)=dh1dh7[0][i];
     285                *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i)=0;
     286                *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1)=dh1dh7[2][i];
     287                *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2)=dh1dh7[1][i];
     288                *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i)=dh1dh7[0][i];
     289                *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1)=dh1dh7[1][i];
     290                *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2)=dh1dh7[2][i];
     291                *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i)=0;
     292                *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1)=0;
     293                *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2)=0;
     294        }
     295
     296        for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
     297                *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
     298                *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
     299                *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
     300                *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0;
     301                *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0;
     302                *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0;
     303                *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=0;
     304                *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=l1l6[i];
    320305        }
    321306
    … …  
    324309/*FUNCTION PentaRef::GetBArtdiff {{{1*/
    325310void PentaRef::GetBArtdiff(double* B_artdiff, double* xyz_list, GaussPenta* gauss){
    326         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
     311        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
    327312         * For grid i, Bi' can be expressed in the actual coordinate system
    328313         * by:
    … …  
    331316         * where h is the interpolation function for grid i.
    332317         *
    333          * We assume B has been allocated already, of size: 2x(DOFPERGRID*numgrids)
    334          */
    335 
    336         int i;
    337         const int calculationdof=3;
    338         const int numgrids=6;
    339         int DOFPERGRID=1;
     318         * We assume B has been allocated already, of size: 2x(NDOF1*NUMNODESP1)
     319         */
    340320
    341321        /*Same thing in the actual coordinate system: */
    342         double dh1dh6[calculationdof][numgrids];
     322        double dh1dh6[3][NUMNODESP1];
    343323
    344324        /*Get dh1dh6 in actual coordinates system : */
    … …  
    346326
    347327        /*Build B': */
    348         for (i=0;i<numgrids;i++){
    349                 *(B_artdiff+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];
    350                 *(B_artdiff+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];
     328        for (int i=0;i<NUMNODESP1;i++){
     329                *(B_artdiff+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i];
     330                *(B_artdiff+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i];
    351331        }
    352332}
    … …  
    354334/*FUNCTION PentaRef::GetBAdvec{{{1*/
    355335void PentaRef::GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss){
    356         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
     336        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
    357337         * For grid i, Bi' can be expressed in the actual coordinate system
    358338         * by:
    … …  
    362342         * where h is the interpolation function for grid i.
    363343         *
    364          * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
    365          */
    366 
    367         int i;
    368         int calculationdof=3;
    369         int numgrids=6;
    370         int DOFPERGRID=1;
     344         * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
     345         */
    371346
    372347        /*Same thing in the actual coordinate system: */
    … …  
    377352
    378353        /*Build B': */
    379         for (i=0;i<numgrids;i++){
    380                 *(B_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=l1l6[i];
    381                 *(B_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=l1l6[i];
    382                 *(B_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=l1l6[i];
     354        for (int i=0;i<NUMNODESP1;i++){
     355                *(B_advec+NDOF1*NUMNODESP1*0+NDOF1*i)=l1l6[i];
     356                *(B_advec+NDOF1*NUMNODESP1*1+NDOF1*i)=l1l6[i];
     357                *(B_advec+NDOF1*NUMNODESP1*2+NDOF1*i)=l1l6[i];
    383358        }
    384359}
    … …  
    386361/*FUNCTION PentaRef::GetBConduct{{{1*/
    387362void PentaRef::GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss){
    388         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
     363        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
    389364         * For grid i, Bi' can be expressed in the actual coordinate system
    390365         * by:
    … …  
    394369         * where h is the interpolation function for grid i.
    395370         *
    396          * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
    397          */
    398 
    399         int i;
    400         const int calculationdof=3;
    401         const int numgrids=6;
    402         int DOFPERGRID=1;
     371         * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
     372         */
    403373
    404374        /*Same thing in the actual coordinate system: */
    405         double dh1dh6[calculationdof][numgrids];
     375        double dh1dh6[3][NUMNODESP1];
    406376
    407377        /*Get dh1dh2dh3 in actual coordinates system : */
    … …  
    409379
    410380        /*Build B': */
    411         for (i=0;i<numgrids;i++){
    412                 *(B_conduct+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];
    413                 *(B_conduct+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];
    414                 *(B_conduct+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i];
     381        for (int i=0;i<NUMNODESP1;i++){
     382                *(B_conduct+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i];
     383                *(B_conduct+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i];
     384                *(B_conduct+NDOF1*NUMNODESP1*2+NDOF1*i)=dh1dh6[2][i];
    415385        }
    416386}
    … …  
    422392
    423393        int i;
    424         const int NDOF3=3;
    425         const int numgrids=6;
    426         double dh1dh6[NDOF3][numgrids];
     394        double dh1dh6[3][NUMNODESP1];
    427395
    428396        /*Get dh1dh6 in actual coordinate system: */
    … …  
    430398
    431399        /*Build B: */
    432         for (i=0;i<numgrids;i++){
     400        for (i=0;i<NUMNODESP1;i++){
    433401                B[i]=dh1dh6[2][i]; 
    434402        }
    … …  
    438406/*FUNCTION PentaRef::GetBprimeAdvec{{{1*/
    439407void PentaRef::GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss){
    440         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
     408        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
    441409         * For grid i, Bi' can be expressed in the actual coordinate system
    442410         * by:
    … …  
    446414         * where h is the interpolation function for grid i.
    447415         *
    448          * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
    449          */
    450 
    451         int i;
    452         const int calculationdof=3;
    453         const int numgrids=6;
    454         int DOFPERGRID=1;
     416         * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
     417         */
    455418
    456419        /*Same thing in the actual coordinate system: */
    457         double dh1dh6[calculationdof][numgrids];
     420        double dh1dh6[3][NUMNODESP1];
    458421
    459422        /*Get dh1dh2dh3 in actual coordinates system : */
    … …  
    461424
    462425        /*Build B': */
    463         for (i=0;i<numgrids;i++){
    464                 *(Bprime_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];
    465                 *(Bprime_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];
    466                 *(Bprime_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i];
     426        for (int i=0;i<NUMNODESP1;i++){
     427                *(Bprime_advec+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i];
     428                *(Bprime_advec+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i];
     429                *(Bprime_advec+NDOF1*NUMNODESP1*2+NDOF1*i)=dh1dh6[2][i];
    467430        }
    468431}
    … …  
    500463
    501464        int i;
    502         const int numgrids2d=3;
    503465        int num_dof=4;
    504466
    505         double l1l2l3[numgrids2d];
     467        double l1l2l3[NUMNODESP1_2d];
    506468
    507469
    … …  
    513475        /*Build LStokes: */
    514476        for (i=0;i<3;i++){
    515                 *(LStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
    516                 *(LStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
    517                 *(LStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
    518                 *(LStokes+num_dof*numgrids2d*0+num_dof*i+3)=0;
    519                 *(LStokes+num_dof*numgrids2d*1+num_dof*i)=0;
    520                 *(LStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i];
    521                 *(LStokes+num_dof*numgrids2d*1+num_dof*i+2)=0;
    522                 *(LStokes+num_dof*numgrids2d*1+num_dof*i+3)=0;
    523                 *(LStokes+num_dof*numgrids2d*2+num_dof*i)=0;
    524                 *(LStokes+num_dof*numgrids2d*2+num_dof*i+1)=0;
    525                 *(LStokes+num_dof*numgrids2d*2+num_dof*i+2)=l1l2l3[i];
    526                 *(LStokes+num_dof*numgrids2d*2+num_dof*i+3)=0;
    527                 *(LStokes+num_dof*numgrids2d*3+num_dof*i)=0;
    528                 *(LStokes+num_dof*numgrids2d*3+num_dof*i+1)=0;
    529                 *(LStokes+num_dof*numgrids2d*3+num_dof*i+2)=l1l2l3[i];
    530                 *(LStokes+num_dof*numgrids2d*3+num_dof*i+3)=0;
    531                 *(LStokes+num_dof*numgrids2d*4+num_dof*i)=l1l2l3[i];
    532                 *(LStokes+num_dof*numgrids2d*4+num_dof*i+1)=0;
    533                 *(LStokes+num_dof*numgrids2d*4+num_dof*i+2)=0;
    534                 *(LStokes+num_dof*numgrids2d*4+num_dof*i+3)=0;
    535                 *(LStokes+num_dof*numgrids2d*5+num_dof*i)=0;
    536                 *(LStokes+num_dof*numgrids2d*5+num_dof*i+1)=l1l2l3[i];
    537                 *(LStokes+num_dof*numgrids2d*5+num_dof*i+2)=0;
    538                 *(LStokes+num_dof*numgrids2d*5+num_dof*i+3)=0;
    539                 *(LStokes+num_dof*numgrids2d*6+num_dof*i)=l1l2l3[i];
    540                 *(LStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;
    541                 *(LStokes+num_dof*numgrids2d*6+num_dof*i+2)=0;
    542                 *(LStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;
    543                 *(LStokes+num_dof*numgrids2d*7+num_dof*i)=0;
    544                 *(LStokes+num_dof*numgrids2d*7+num_dof*i+1)=l1l2l3[i];
    545                 *(LStokes+num_dof*numgrids2d*7+num_dof*i+2)=0;
    546                 *(LStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;
    547                 *(LStokes+num_dof*numgrids2d*8+num_dof*i)=0;
    548                 *(LStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;
    549                 *(LStokes+num_dof*numgrids2d*8+num_dof*i+2)=l1l2l3[i];
    550                 *(LStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;
    551                 *(LStokes+num_dof*numgrids2d*9+num_dof*i)=0;
    552                 *(LStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;
    553                 *(LStokes+num_dof*numgrids2d*9+num_dof*i+2)=l1l2l3[i];
    554                 *(LStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;
    555                 *(LStokes+num_dof*numgrids2d*10+num_dof*i)=0;
    556                 *(LStokes+num_dof*numgrids2d*10+num_dof*i+1)=0;
    557                 *(LStokes+num_dof*numgrids2d*10+num_dof*i+2)=l1l2l3[i];
    558                 *(LStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;
    559                 *(LStokes+num_dof*numgrids2d*11+num_dof*i)=l1l2l3[i];
    560                 *(LStokes+num_dof*numgrids2d*11+num_dof*i+1)=0;
    561                 *(LStokes+num_dof*numgrids2d*11+num_dof*i+2)=0;
    562                 *(LStokes+num_dof*numgrids2d*11+num_dof*i+3)=0;
    563                 *(LStokes+num_dof*numgrids2d*12+num_dof*i)=0;
    564                 *(LStokes+num_dof*numgrids2d*12+num_dof*i+1)=l1l2l3[i];
    565                 *(LStokes+num_dof*numgrids2d*12+num_dof*i+2)=0;
    566                 *(LStokes+num_dof*numgrids2d*12+num_dof*i+3)=0;
    567                 *(LStokes+num_dof*numgrids2d*13+num_dof*i)=0;
    568                 *(LStokes+num_dof*numgrids2d*13+num_dof*i+1)=0;
    569                 *(LStokes+num_dof*numgrids2d*13+num_dof*i+2)=l1l2l3[i];
    570                 *(LStokes+num_dof*numgrids2d*13+num_dof*i+3)=0;
     477                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
     478                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
     479                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
     480                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
     481                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
     482                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
     483                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
     484                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
     485                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0;
     486                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
     487                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i];
     488                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
     489                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
     490                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0;
     491                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i];
     492                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
     493                *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=l1l2l3[i];
     494                *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
     495                *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=0;
     496                *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0;
     497                *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
     498                *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=l1l2l3[i];
     499                *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=0;
     500                *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0;
     501                *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=l1l2l3[i];
     502                *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
     503                *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=0;
     504                *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=0;
     505                *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
     506                *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=l1l2l3[i];
     507                *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=0;
     508                *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=0;
     509                *(LStokes+num_dof*NUMNODESP1_2d*8+num_dof*i)=0;
     510                *(LStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+1)=0;
     511                *(LStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+2)=l1l2l3[i];
     512                *(LStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+3)=0;
     513                *(LStokes+num_dof*NUMNODESP1_2d*9+num_dof*i)=0;
     514                *(LStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+1)=0;
     515                *(LStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+2)=l1l2l3[i];
     516                *(LStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+3)=0;
     517                *(LStokes+num_dof*NUMNODESP1_2d*10+num_dof*i)=0;
     518                *(LStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+1)=0;
     519                *(LStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+2)=l1l2l3[i];
     520                *(LStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+3)=0;
     521                *(LStokes+num_dof*NUMNODESP1_2d*11+num_dof*i)=l1l2l3[i];
     522                *(LStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+1)=0;
     523                *(LStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+2)=0;
     524                *(LStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+3)=0;
     525                *(LStokes+num_dof*NUMNODESP1_2d*12+num_dof*i)=0;
     526                *(LStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+1)=l1l2l3[i];
     527                *(LStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+2)=0;
     528                *(LStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+3)=0;
     529                *(LStokes+num_dof*NUMNODESP1_2d*13+num_dof*i)=0;
     530                *(LStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+1)=0;
     531                *(LStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+2)=l1l2l3[i];
     532                *(LStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+3)=0;
    571533
    572534        }
    … …  
    597559         */
    598560        int i;
    599         const int numgrids2d=3;
    600561        int num_dof=4;
    601562
    602         double l1l2l3[numgrids2d];
    603         double dh1dh6[3][6];
     563        double l1l2l3[NUMNODESP1_2d];
     564        double dh1dh6[3][NUMNODESP1];
    604565
    605566        /*Get l1l2l3 in actual coordinate system: */
    … …  
    612573        /*Build LprimeStokes: */
    613574        for (i=0;i<3;i++){
    614                 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
    615                 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
    616                 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
    617                 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+3)=0;
    618                 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i)=0;
    619                 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i];
    620                 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+2)=0;
    621                 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+3)=0;
    622                 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i)=l1l2l3[i];
    623                 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+1)=0;
    624                 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+2)=0;
    625                 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+3)=0;
    626                 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i)=0;
    627                 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+1)=l1l2l3[i];
    628                 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+2)=0;
    629                 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+3)=0;
    630                 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i)=0;
    631                 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+1)=0;
    632                 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+2)=l1l2l3[i];
    633                 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+3)=0;
    634                 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i)=0;
    635                 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+1)=0;
    636                 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+2)=l1l2l3[i];
    637                 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+3)=0;
    638                 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i)=0;
    639                 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;
    640                 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+2)=dh1dh6[2][i];
    641                 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;
    642                 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i)=0;
    643                 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+1)=0;
    644                 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+2)=dh1dh6[2][i];
    645                 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;
    646                 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i)=0;
    647                 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;
    648                 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+2)=dh1dh6[2][i];
    649                 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;
    650                 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i)=dh1dh6[2][i];
    651                 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;
    652                 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+2)=dh1dh6[0][i];
    653                 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;
    654                 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i)=0;
    655                 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+1)=dh1dh6[2][i];
    656                 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+2)=dh1dh6[1][i];
    657                 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;
    658                 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i)=0;
    659                 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+1)=0;
    660                 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+2)=0;
    661                 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+3)=l1l2l3[i];
    662                 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i)=0;
    663                 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+1)=0;
    664                 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+2)=0;
    665                 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+3)=l1l2l3[i];
    666                 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i)=0;
    667                 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+1)=0;
    668                 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+2)=0;
    669                 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+3)=l1l2l3[i];
    670 
     575                *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
     576                *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
     577                *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
     578                *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
     579                *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
     580                *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
     581                *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
     582                *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
     583                *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
     584                *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
     585                *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=0;
     586                *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
     587                *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
     588                *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
     589                *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=0;
     590                *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
     591                *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=0;
     592                *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
     593                *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=l1l2l3[i];
     594                *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0;
     595                *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
     596                *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=0;
     597                *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=l1l2l3[i];
     598                *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0;
     599                *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=0;
     600                *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
     601                *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=dh1dh6[2][i];
     602                *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=0;
     603                *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
     604                *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=0;
     605                *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=dh1dh6[2][i];
     606                *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=0;
     607                *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i)=0;
     608                *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+1)=0;
     609                *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+2)=dh1dh6[2][i];
     610                *(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+3)=0;
     611                *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i)=dh1dh6[2][i];
     612                *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+1)=0;
     613                *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+2)=dh1dh6[0][i];
     614                *(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+3)=0;
     615                *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i)=0;
     616                *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+1)=dh1dh6[2][i];
     617                *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+2)=dh1dh6[1][i];
     618                *(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+3)=0;
     619                *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i)=0;
     620                *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+1)=0;
     621                *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+2)=0;
     622                *(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+3)=l1l2l3[i];
     623                *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i)=0;
     624                *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+1)=0;
     625                *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+2)=0;
     626                *(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+3)=l1l2l3[i];
     627                *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i)=0;
     628                *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+1)=0;
     629                *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+2)=0;
     630                *(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+3)=l1l2l3[i];
    671631        }
    672632}
    … …  
    675635void PentaRef::GetJacobian(double* J, double* xyz_list,GaussPenta* gauss){
    676636
    677         const int NDOF3=3;
    678637        int i,j;
    679638
    … …  
    822781
    823782        int       i;
    824         const int numgrids = 7;
    825         double    dh1dh7_ref[3][numgrids];
     783        double    dh1dh7_ref[3][NUMNODESMINI];
    826784        double    Jinv[3][3];
    827785
    … …  
    839797         */
    840798
    841         for (i=0;i<numgrids;i++){
    842                 *(dh1dh7+numgrids*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];
    843                 *(dh1dh7+numgrids*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];
    844                 *(dh1dh7+numgrids*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i];
     799        for (i=0;i<NUMNODESMINI;i++){
     800                *(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];
     801                *(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];
     802                *(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];
    845803        }
    846804
    … …  
    852810        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    853811         * natural coordinate system) at the gaussian point. */
    854 
    855         int    numgrids=7; //six plus bubble grids
    856 
    857812        double r=gauss->coord2-gauss->coord1;
    858813        double s=-3.0/SQRT3*(gauss->coord1+gauss->coord2-2.0/3.0);
    … …  
    860815
    861816        /*First nodal function: */
    862         *(dl1dl7+numgrids*0+0)=-0.5*(1.0-zeta)/2.0;
    863         *(dl1dl7+numgrids*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0;
    864         *(dl1dl7+numgrids*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
     817        *(dl1dl7+NUMNODESMINI*0+0)=-0.5*(1.0-zeta)/2.0;
     818        *(dl1dl7+NUMNODESMINI*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0;
     819        *(dl1dl7+NUMNODESMINI*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
    865820
    866821        /*Second nodal function: */
    867         *(dl1dl7+numgrids*0+1)=0.5*(1.0-zeta)/2.0;
    868         *(dl1dl7+numgrids*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0;
    869         *(dl1dl7+numgrids*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
     822        *(dl1dl7+NUMNODESMINI*0+1)=0.5*(1.0-zeta)/2.0;
     823        *(dl1dl7+NUMNODESMINI*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0;
     824        *(dl1dl7+NUMNODESMINI*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
    870825
    871826        /*Third nodal function: */
    872         *(dl1dl7+numgrids*0+2)=0;
    873         *(dl1dl7+numgrids*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;
    874         *(dl1dl7+numgrids*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);
     827        *(dl1dl7+NUMNODESMINI*0+2)=0;
     828        *(dl1dl7+NUMNODESMINI*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;
     829        *(dl1dl7+NUMNODESMINI*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);
    875830
    876831        /*Fourth nodal function: */
    877         *(dl1dl7+numgrids*0+3)=-0.5*(1.0+zeta)/2.0;
    878         *(dl1dl7+numgrids*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0;
    879         *(dl1dl7+numgrids*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
     832        *(dl1dl7+NUMNODESMINI*0+3)=-0.5*(1.0+zeta)/2.0;
     833        *(dl1dl7+NUMNODESMINI*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0;
     834        *(dl1dl7+NUMNODESMINI*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
    880835
    881836        /*Fith nodal function: */
    882         *(dl1dl7+numgrids*0+4)=0.5*(1.0+zeta)/2.0;
    883         *(dl1dl7+numgrids*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0;
    884         *(dl1dl7+numgrids*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
     837        *(dl1dl7+NUMNODESMINI*0+4)=0.5*(1.0+zeta)/2.0;
     838        *(dl1dl7+NUMNODESMINI*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0;
     839        *(dl1dl7+NUMNODESMINI*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
    885840
    886841        /*Sixth nodal function: */
    887         *(dl1dl7+numgrids*0+5)=0;
    888         *(dl1dl7+numgrids*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;
    889         *(dl1dl7+numgrids*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD);
     842        *(dl1dl7+NUMNODESMINI*0+5)=0;
     843        *(dl1dl7+NUMNODESMINI*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;
     844        *(dl1dl7+NUMNODESMINI*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD);
    890845
    891846        /*Seventh nodal function: */
    892         *(dl1dl7+numgrids*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);
    893         *(dl1dl7+numgrids*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));
    894         *(dl1dl7+numgrids*2+6)=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);
     847        *(dl1dl7+NUMNODESMINI*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);
     848        *(dl1dl7+NUMNODESMINI*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));
     849        *(dl1dl7+NUMNODESMINI*2+6)=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);
    895850
    896851}
    … …  
    914869        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    915870         * actual coordinate system): */
    916         int       i;
    917         const int NDOF3    = 3;
    918         const int numgrids = 6;
    919         double    dh1dh6_ref[NDOF3][numgrids];
     871        double    dh1dh6_ref[NDOF3][NUMNODESP1];
    920872        double    Jinv[NDOF3][NDOF3];
    921873
    … …  
    933885         */
    934886
    935         for (i=0;i<numgrids;i++){
    936                 *(dh1dh6+numgrids*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i];
    937                 *(dh1dh6+numgrids*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i];
    938                 *(dh1dh6+numgrids*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i];
     887        for (int i=0;i<NUMNODESP1;i++){
     888                *(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];
     889                *(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];
     890                *(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];
    939891        }
    940892
    … …  
    947899         * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */
    948900
    949         const int numgrids=6;
    950901        double A1,A2,A3,z;
    951902
    … …  
    956907
    957908        /*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
    958         *(dl1dl6+numgrids*0+0)=-0.5*(1.0-z)/2.0;
    959         *(dl1dl6+numgrids*1+0)=-0.5/SQRT3*(1.0-z)/2.0;
    960         *(dl1dl6+numgrids*2+0)=-0.5*A1;
     909        *(dl1dl6+NUMNODESP1*0+0)=-0.5*(1.0-z)/2.0;
     910        *(dl1dl6+NUMNODESP1*1+0)=-0.5/SQRT3*(1.0-z)/2.0;
     911        *(dl1dl6+NUMNODESP1*2+0)=-0.5*A1;
    961912
    962913        /*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
    963         *(dl1dl6+numgrids*0+1)=0.5*(1.0-z)/2.0;
    964         *(dl1dl6+numgrids*1+1)=-0.5/SQRT3*(1.0-z)/2.0;
    965         *(dl1dl6+numgrids*2+1)=-0.5*A2;
     914        *(dl1dl6+NUMNODESP1*0+1)=0.5*(1.0-z)/2.0;
     915        *(dl1dl6+NUMNODESP1*1+1)=-0.5/SQRT3*(1.0-z)/2.0;
     916        *(dl1dl6+NUMNODESP1*2+1)=-0.5*A2;
    966917
    967918        /*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/
    968         *(dl1dl6+numgrids*0+2)=0.0;
    969         *(dl1dl6+numgrids*1+2)=1.0/SQRT3*(1.0-z)/2.0;
    970         *(dl1dl6+numgrids*2+2)=-0.5*A3;
     919        *(dl1dl6+NUMNODESP1*0+2)=0.0;
     920        *(dl1dl6+NUMNODESP1*1+2)=1.0/SQRT3*(1.0-z)/2.0;
     921        *(dl1dl6+NUMNODESP1*2+2)=-0.5*A3;
    971922
    972923        /*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/
    973         *(dl1dl6+numgrids*0+3)=-0.5*(1.0+z)/2.0;
    974         *(dl1dl6+numgrids*1+3)=-0.5/SQRT3*(1.0+z)/2.0;
    975         *(dl1dl6+numgrids*2+3)=0.5*A1;
     924        *(dl1dl6+NUMNODESP1*0+3)=-0.5*(1.0+z)/2.0;
     925        *(dl1dl6+NUMNODESP1*1+3)=-0.5/SQRT3*(1.0+z)/2.0;
     926        *(dl1dl6+NUMNODESP1*2+3)=0.5*A1;
    976927
    977928        /*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/
    978         *(dl1dl6+numgrids*0+4)=0.5*(1.0+z)/2.0;
    979         *(dl1dl6+numgrids*1+4)=-0.5/SQRT3*(1.0+z)/2.0;
    980         *(dl1dl6+numgrids*2+4)=0.5*A2;
     929        *(dl1dl6+NUMNODESP1*0+4)=0.5*(1.0+z)/2.0;
     930        *(dl1dl6+NUMNODESP1*1+4)=-0.5/SQRT3*(1.0+z)/2.0;
     931        *(dl1dl6+NUMNODESP1*2+4)=0.5*A2;
    981932
    982933        /*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/
    983         *(dl1dl6+numgrids*0+5)=0.0;
    984         *(dl1dl6+numgrids*1+5)=1.0/SQRT3*(1.0+z)/2.0;
    985         *(dl1dl6+numgrids*2+5)=0.5*A3;
     934        *(dl1dl6+NUMNODESP1*0+5)=0.0;
     935        *(dl1dl6+NUMNODESP1*1+5)=1.0/SQRT3*(1.0+z)/2.0;
     936        *(dl1dl6+NUMNODESP1*2+5)=0.5*A3;
    986937}
    987938/*}}}*/
    … …  
    1010961         *   p is a vector of size 3x1 already allocated.
    1011962         */
    1012         double dh1dh6[3][6];
     963        double dh1dh6[3][NUMNODESP1];
    1013964
    1014965        /*Get nodal funnctions derivatives in actual coordinate system: */

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

    r5743 r5745  
    489489
    490490                                        /*build new bed and surface: */
    491                                         if (this->GetShelf()){
     491                                        if (this->IsOnShelf()){
    492492                                                /*hydrostatic equilibrium: */
    493493                                                double rho_ice,rho_water,di;
    … …  
    819819}
    820820/*}}}*/
    821 /*FUNCTION Tria::GetOnBed {{{1*/
    822 bool Tria::GetOnBed(){
     821/*FUNCTION Tria::IsOnBed {{{1*/
     822bool Tria::IsOnBed(){
    823823       
    824824        bool onbed;
    … …  
    828828}
    829829/*}}}*/
    830 /*FUNCTION Tria::GetShelf {{{1*/
    831 bool   Tria::GetShelf(){
     830/*FUNCTION Tria::IsOnShelf {{{1*/
     831bool   Tria::IsOnShelf(){
    832832        /*inputs: */
    833833        bool shelf;
    … …  
    24212421
    24222422        /*If shelf: hydrostatic equilibrium*/
    2423         if (this->GetShelf()){
     2423        if (this->IsOnShelf()){
    24242424
    24252425                /*recover material parameters: */

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

    r5743 r5745  
    7575                void   CreatePVector(Vec pg);
    7676                int    GetNodeIndex(Node* node);
    77                 bool   GetOnBed();
    78                 bool   GetShelf();
     77                bool   IsOnBed();
     78                bool   IsOnShelf();
    7979                void   GetSolutionFromInputs(Vec solution);
    8080                void   GetVectorFromInputs(Vec vector,int NameEnum);

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

    r5743 r5745  
    2020/*}}}*/
    2121
     22/*Element macros*/
     23#define NUMNODES 3
     24#define NDOF1 1
     25#define NDOF2 2
     26#define NDOF3 3
     27#define NDOF4 4
     28
    2229/*Object constructors and destructor*/
    2330/*FUNCTION TriaRef::TriaRef(){{{1*/
    … …  
    2734/*}}}*/
    2835/*FUNCTION TriaRef::TriaRef(int* types,int nummodels){{{1*/
     36
    2937TriaRef::TriaRef(const int nummodels){
    3038
    … …  
    6270         * where h is the interpolation function for grid i.
    6371         *
    64          * We assume B has been allocated already, of size: 3x(NDOF2*numgrids)
     72         * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES)
    6573         */
    6674
    6775        int i;
    68         const int NDOF2=2;
    69         const int numgrids=3;
    70 
    71         double dh1dh3[NDOF2][numgrids];
    72 
     76        double dh1dh3[NDOF2][NUMNODES];
    7377
    7478        /*Get dh1dh2dh3 in actual coordinate system: */
    … …  
    7680
    7781        /*Build B: */
    78         for (i=0;i<numgrids;i++){
    79                 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];
    80                 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0;
    81                 *(B+NDOF2*numgrids*1+NDOF2*i)=0;
    82                 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh3[1][i];
    83                 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i];
    84                 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i];
     82        for (i=0;i<NUMNODES;i++){
     83                *(B+NDOF2*NUMNODES*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];
     84                *(B+NDOF2*NUMNODES*0+NDOF2*i+1)=0;
     85                *(B+NDOF2*NUMNODES*1+NDOF2*i)=0;
     86                *(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dh1dh3[1][i];
     87                *(B+NDOF2*NUMNODES*2+NDOF2*i)=(float).5*dh1dh3[1][i];
     88                *(B+NDOF2*NUMNODES*2+NDOF2*i+1)=(float).5*dh1dh3[0][i];
    8589        }
    8690}
    … …  
    9599         */
    96100
    97         const int numgrids=3;
    98         double l1l3[numgrids];
     101        double l1l3[NUMNODES];
    99102
    100103        GetNodalFunctions(&l1l3[0],gauss);
    … …  
    115118         */
    116119
    117         const int numgrids=3;
    118         double l1l3[numgrids];
     120        double l1l3[NUMNODES];
    119121
    120122        GetNodalFunctions(&l1l3[0],gauss);
    … …  
    135137         * where h is the interpolation function for grid i.
    136138         *
    137          * We assume B_prog has been allocated already, of size: 2x(NDOF1*numgrids)
    138          */
    139 
    140         int i;
    141         const int NDOF1=1;
    142         const int numgrids=3;
    143 
    144         double l1l2l3[numgrids];
    145 
     139         * We assume B_prog has been allocated already, of size: 2x(NDOF1*NUMNODES)
     140         */
     141
     142        double l1l2l3[NUMNODES];
    146143
    147144        /*Get dh1dh2dh3 in actual coordinate system: */
    … …  
    149146
    150147        /*Build B_prog: */
    151         for (i=0;i<numgrids;i++){
    152                 *(B_prog+NDOF1*numgrids*0+NDOF1*i)=l1l2l3[i];
    153                 *(B_prog+NDOF1*numgrids*1+NDOF1*i)=l1l2l3[i];
     148        for (int i=0;i<NUMNODES;i++){
     149                *(B_prog+NDOF1*NUMNODES*0+NDOF1*i)=l1l2l3[i];
     150                *(B_prog+NDOF1*NUMNODES*1+NDOF1*i)=l1l2l3[i];
    154151        }
    155152}
    … …  
    166163         * where h is the interpolation function for grid i.
    167164         *
    168          * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
    169          */
    170 
    171         int i;
    172         const int NDOF2=2;
    173         const int numgrids=3;
     165         * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODES)
     166         */
    174167
    175168        /*Same thing in the actual coordinate system: */
    176         double dh1dh3[NDOF2][numgrids];
     169        double dh1dh3[NDOF2][NUMNODES];
    177170
    178171        /*Get dh1dh2dh3 in actual coordinates system : */
    … …  
    180173
    181174        /*Build B': */
    182         for (i=0;i<numgrids;i++){
    183                 *(Bprime+NDOF2*numgrids*0+NDOF2*i)=2*dh1dh3[0][i];
    184                 *(Bprime+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh3[1][i];
    185                 *(Bprime+NDOF2*numgrids*1+NDOF2*i)=dh1dh3[0][i];
    186                 *(Bprime+NDOF2*numgrids*1+NDOF2*i+1)=2*dh1dh3[1][i];
    187                 *(Bprime+NDOF2*numgrids*2+NDOF2*i)=dh1dh3[1][i];
    188                 *(Bprime+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh3[0][i];
     175        for (int i=0;i<NUMNODES;i++){
     176                *(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=2*dh1dh3[0][i];
     177                *(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=dh1dh3[1][i];
     178                *(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=dh1dh3[0][i];
     179                *(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=2*dh1dh3[1][i];
     180                *(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dh1dh3[1][i];
     181                *(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dh1dh3[0][i];
    189182        }
    190183}
    … …  
    199192         * where h is the interpolation function for grid i.
    200193         *
    201          * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
    202          */
    203 
    204         int i;
    205         const int NDOF1=1;
    206         const int NDOF2=2;
    207         const int numgrids=3;
     194         * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODES)
     195         */
    208196
    209197        /*Same thing in the actual coordinate system: */
    210         double dh1dh3[NDOF2][numgrids];
     198        double dh1dh3[NDOF2][NUMNODES];
    211199
    212200        /*Get dh1dh2dh3 in actual coordinates system : */
    … …  
    214202
    215203        /*Build B': */
    216         for (i=0;i<numgrids;i++){
    217                 *(Bprime_prog+NDOF1*numgrids*0+NDOF1*i)=dh1dh3[0][i];
    218                 *(Bprime_prog+NDOF1*numgrids*1+NDOF1*i)=dh1dh3[1][i];
     204        for (int i=0;i<NUMNODES;i++){
     205                *(Bprime_prog+NDOF1*NUMNODES*0+NDOF1*i)=dh1dh3[0][i];
     206                *(Bprime_prog+NDOF1*NUMNODES*1+NDOF1*i)=dh1dh3[1][i];
    219207        }
    220208}
    … …  
    232220         * where h is the interpolation function for grid i.
    233221         *
    234          * We assume L has been allocated already, of size: numgrids (numdof=1), or numdofx(numdof*numgrids) (numdof=2)
     222         * We assume L has been allocated already, of size: NUMNODES (numdof=1), or numdofx(numdof*NUMNODES) (numdof=2)
    235223         */
    236224
    237225        int i;
    238         const int NDOF2=2;
    239         const int numgrids=3;
    240226        double l1l2l3[3];
    241227
    … …  
    245231        /*Build L: */
    246232        if(numdof==1){
    247                 for (i=0;i<numgrids;i++){
     233                for (i=0;i<NUMNODES;i++){
    248234                        L[i]=l1l2l3[i];
    249235                }
    250236        }
    251237        else{
    252                 for (i=0;i<numgrids;i++){
    253                         *(L+numdof*numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i];
    254                         *(L+numdof*numgrids*0+numdof*i+1)=0;
    255                         *(L+numdof*numgrids*1+numdof*i)=0;
    256                         *(L+numdof*numgrids*1+numdof*i+1)=l1l2l3[i];
     238                for (i=0;i<NUMNODES;i++){
     239                        *(L+numdof*NUMNODES*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i];
     240                        *(L+numdof*NUMNODES*0+numdof*i+1)=0;
     241                        *(L+numdof*NUMNODES*1+numdof*i)=0;
     242                        *(L+numdof*NUMNODES*1+numdof*i+1)=l1l2l3[i];
    257243                }
    258244        }
    … …  
    263249        /*The Jacobian is constant over the element, discard the gaussian points.
    264250         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    265 
    266         const int NDOF2=2;
    267         const int numgrids=3;
    268251        double x1,y1,x2,y2,x3,y3;
    269252
    270         x1=*(xyz_list+numgrids*0+0);
    271         y1=*(xyz_list+numgrids*0+1);
    272         x2=*(xyz_list+numgrids*1+0);
    273         y2=*(xyz_list+numgrids*1+1);
    274         x3=*(xyz_list+numgrids*2+0);
    275         y3=*(xyz_list+numgrids*2+1);
     253        x1=*(xyz_list+NUMNODES*0+0);
     254        y1=*(xyz_list+NUMNODES*0+1);
     255        x2=*(xyz_list+NUMNODES*1+0);
     256        y2=*(xyz_list+NUMNODES*1+1);
     257        x3=*(xyz_list+NUMNODES*2+0);
     258        y3=*(xyz_list+NUMNODES*2+1);
    276259
    277260
    … …  
    388371         * actual coordinate system): */
    389372        int       i;
    390         const int NDOF2    = 2;
    391         const int numgrids = 3;
    392         double    dh1dh3_ref[NDOF2][numgrids];
     373        double    dh1dh3_ref[NDOF2][NUMNODES];
    393374        double    Jinv[NDOF2][NDOF2];
    394375
    … …  
    404385         * [dhi/dy]       [dhi/ds]
    405386         */
    406         for (i=0;i<numgrids;i++){
    407                 dh1dh3[numgrids*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
    408                 dh1dh3[numgrids*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
     387        for (i=0;i<NUMNODES;i++){
     388                dh1dh3[NUMNODES*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
     389                dh1dh3[NUMNODES*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
    409390        }
    410391
    … …  
    413394/*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference{{{1*/
    414395void TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss){
    415 
    416396        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    417397         * natural coordinate system) at the gaussian point. */
    418398
    419         const int NDOF2=2;
    420         const int numgrids=3;
    421 
    422399        /*First nodal function: */
    423         *(dl1dl3+numgrids*0+0)=-0.5;
    424         *(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3);
     400        *(dl1dl3+NUMNODES*0+0)=-0.5;
     401        *(dl1dl3+NUMNODES*1+0)=-1.0/(2.0*SQRT3);
    425402
    426403        /*Second nodal function: */
    427         *(dl1dl3+numgrids*0+1)=0.5;
    428         *(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3);
     404        *(dl1dl3+NUMNODES*0+1)=0.5;
     405        *(dl1dl3+NUMNODES*1+1)=-1.0/(2.0*SQRT3);
    429406
    430407        /*Third nodal function: */
    431         *(dl1dl3+numgrids*0+2)=0;
    432         *(dl1dl3+numgrids*1+2)=1.0/SQRT3;
     408        *(dl1dl3+NUMNODES*0+2)=0;
     409        *(dl1dl3+NUMNODES*1+2)=1.0/SQRT3;
    433410
    434411}
Note: See TracChangeset for help on using the changeset viewer.