Changeset 18925


Ignore:
Timestamp:
12/04/14 08:48:31 (10 years ago)
Author:
seroussi
Message:

CHG: finished reordering elements

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

Legend:

Unmodified
Added
Removed

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

    r18243 r18925  
    7575/*}}}*/
    7676
     77void ElementHook::InitHookNeighbors(int* element_ids){/*{{{*/
     78        this->hneighbors=new Hook(element_ids,2);
     79}
     80/*}}}*/
    7781void ElementHook::SetHookNodes(int* node_ids,int numnodes,int analysis_counter){/*{{{*/
    7882        if(this->hnodes) this->hnodes[analysis_counter]= new Hook(node_ids,numnodes);
    7983}
    8084/*}}}*/
    81 void ElementHook::InitHookNeighbors(int* element_ids){/*{{{*/
    82         this->hneighbors=new Hook(element_ids,2);
     85void ElementHook::SpawnSegHook(ElementHook* triahook,int index1,int index2){/*{{{*/
     86
     87        triahook->numanalyses=this->numanalyses;
     88
     89        int indices[2];
     90        indices[0]=index1;
     91        indices[1]=index2;
     92
     93        /*Spawn nodes hook*/
     94        triahook->hnodes=new Hook*[this->numanalyses];
     95        for(int i=0;i<this->numanalyses;i++){
     96                /*Do not do anything if Hook is empty*/
     97                if (!this->hnodes[i] || this->hnodes[i]->GetNum()==0){
     98                        triahook->hnodes[i]=NULL;
     99                }
     100                else{
     101                        triahook->hnodes[i]=this->hnodes[i]->Spawn(indices,2);
     102                }
     103        }
     104
     105        /*do not spawn hmaterial. material will be taken care of by Tria*/
     106        triahook->hmaterial=NULL;
     107        triahook->hvertices=(Hook*)this->hvertices->Spawn(indices,2);
     108        triahook->hmatpar=(Hook*)this->hmatpar->copy();
    83109}
    84110/*}}}*/
    … …  
    111137}
    112138/*}}}*/
    113 void ElementHook::SpawnSegHook(ElementHook* triahook,int index1,int index2){/*{{{*/
    114 
    115         triahook->numanalyses=this->numanalyses;
    116 
    117         int indices[2];
    118         indices[0]=index1;
    119         indices[1]=index2;
    120 
    121         /*Spawn nodes hook*/
    122         triahook->hnodes=new Hook*[this->numanalyses];
    123         for(int i=0;i<this->numanalyses;i++){
    124                 /*Do not do anything if Hook is empty*/
    125                 if (!this->hnodes[i] || this->hnodes[i]->GetNum()==0){
    126                         triahook->hnodes[i]=NULL;
    127                 }
    128                 else{
    129                         triahook->hnodes[i]=this->hnodes[i]->Spawn(indices,2);
    130                 }
    131         }
    132 
    133         /*do not spawn hmaterial. material will be taken care of by Tria*/
    134         triahook->hmaterial=NULL;
    135         triahook->hvertices=(Hook*)this->hvertices->Spawn(indices,2);
    136         triahook->hmatpar=(Hook*)this->hmatpar->copy();
    137 }
    138 /*}}}*/

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

    r17514 r18925  
    2424                ~ElementHook();
    2525
     26                void InitHookNeighbors(int* element_ids);               //3d only
    2627                void SetHookNodes(int* node_ids,int numnodes,int analysis_counter);
     28                void SpawnSegHook(ElementHook* triahook,int ndex1,int index2); //2d only
    2729                void SpawnTriaHook(ElementHook* triahook,int index1,int index2,int index3); //3d only
    28                 void SpawnSegHook(ElementHook* triahook,int ndex1,int index2); //2d only
    29                 void InitHookNeighbors(int* element_ids);               //3d only
    3030};
    3131

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

    r18523 r18925  
    3737
    3838/*Reference Element numerics*/
     39void PentaRef::BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){/*{{{*/
     40
     41        /*Output*/
     42        int  numindices;
     43        int* indices = NULL;
     44
     45        switch(finiteelement){
     46                case P1Enum: case P1DGEnum:
     47                        numindices = 3;
     48                        indices    = xNew<int>(numindices);
     49                        indices[0] = 0;
     50                        indices[1] = 1;
     51                        indices[2] = 2;
     52                        break;
     53                case P1bubbleEnum: case P1bubblecondensedEnum:
     54                        numindices = 3;
     55                        indices    = xNew<int>(numindices);
     56                        indices[0] = 0;
     57                        indices[1] = 1;
     58                        indices[2] = 2;
     59                        break;
     60                case P2xP1Enum:
     61                        numindices = 6;
     62                        indices    = xNew<int>(numindices);
     63                        indices[0] = 0;
     64                        indices[1] = 1;
     65                        indices[2] = 2;
     66                        indices[3] = 6;
     67                        indices[4] = 7;
     68                        indices[5] = 8;
     69                        break;
     70                case P1xP2Enum:
     71                        numindices = 3;
     72                        indices    = xNew<int>(numindices);
     73                        indices[0] = 0;
     74                        indices[1] = 1;
     75                        indices[2] = 2;
     76                        break;
     77                case P1xP3Enum:
     78                        numindices = 3;
     79                        indices    = xNew<int>(numindices);
     80                        indices[0] = 0;
     81                        indices[1] = 1;
     82                        indices[2] = 2;
     83                        break;
     84                case P2Enum:
     85                        numindices = 6;
     86                        indices    = xNew<int>(numindices);
     87                        indices[0] = 0;
     88                        indices[1] = 1;
     89                        indices[2] = 2;
     90                        indices[3] = 9;
     91                        indices[4] = 10;
     92                        indices[5] = 11;
     93                        break;
     94                case P2bubbleEnum:
     95                        numindices = 6;
     96                        indices    = xNew<int>(numindices);
     97                        indices[0] = 0;
     98                        indices[1] = 1;
     99                        indices[2] = 2;
     100                        indices[3] = 9;
     101                        indices[4] = 10;
     102                        indices[5] = 11;
     103                        break;
     104                case P2xP4Enum:
     105                        numindices = 6;
     106                        indices    = xNew<int>(numindices);
     107                        indices[0] = 0;
     108                        indices[1] = 1;
     109                        indices[2] = 2;
     110                        indices[3] = 9;
     111                        indices[4] = 10;
     112                        indices[5] = 11;
     113                        break;
     114                default:
     115                        _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
     116        }
     117
     118        /*Assign output pointer*/
     119        *pnumindices = numindices;
     120        *pindices    = indices;
     121}
     122/*}}}*/
     123void PentaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
     124        /*From node values of parameter p (p_list[0], p_list[1], p_list[2],
     125         * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
     126         * gaussian point specified by gauss_coord:
     127         *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
     128         *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
     129         *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
     130         *
     131         *   p is a vector of size 3x1 already allocated.
     132         *
     133         * WARNING: For a significant gain in performance, it is better to use
     134         * static memory allocation instead of dynamic.
     135         */
     136
     137        /*Allocate derivatives of basis functions*/
     138        IssmDouble  dbasis[3*NUMNODESMAX];
     139
     140        /*Fetch number of nodes for this finite element*/
     141        int numnodes = this->NumberofNodes(finiteelement);
     142        _assert_(numnodes<=NUMNODESMAX);
     143
     144        /*Get basis functions derivatives at this point*/
     145        GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
     146
     147        /*Calculate parameter for this Gauss point*/
     148        IssmDouble dpx=0.;
     149        IssmDouble dpy=0.;
     150        IssmDouble dpz=0.;
     151        for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
     152        for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
     153        for(int i=0;i<numnodes;i++) dpz += dbasis[2*numnodes+i]*plist[i];
     154
     155        /*Assign values*/
     156        p[0]=dpx;
     157        p[1]=dpy;
     158        p[2]=dpz;
     159}
     160/*}}}*/
     161void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss,int finiteelement){/*{{{*/
     162        /* WARNING: For a significant gain in performance, it is better to use
     163         * static memory allocation instead of dynamic.*/
     164
     165        /*Allocate basis functions*/
     166        IssmDouble  basis[NUMNODESMAX];
     167
     168        /*Fetch number of nodes for this finite element*/
     169        int numnodes = this->NumberofNodes(finiteelement);
     170        _assert_(numnodes<=NUMNODESMAX);
     171
     172        /*Get basis functions at this point*/
     173        GetNodalFunctions(&basis[0],gauss,finiteelement);
     174
     175        /*Calculate parameter for this Gauss point*/
     176        IssmDouble value =0.;
     177        for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
     178
     179        /*Assign output pointer*/
     180        *pvalue = value;
     181}
     182/*}}}*/
    39183void PentaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
    40184        /*The Jacobian is constant over the element, discard the gaussian points.
    … …  
    103247        /*Get Determinant*/
    104248        Matrix3x3Determinant(Jdet,&J[0][0]);
    105         if(*Jdet<0) _error_("negative jacobian determinant!");
    106 
    107 }
    108 /*}}}*/
    109 void PentaRef::GetTriaJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    110         /*The Jacobian determinant is constant over the element, discard the gaussian points.
    111          * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    112 
    113         IssmDouble x1=xyz_list[3*0+0];
    114         IssmDouble y1=xyz_list[3*0+1];
    115         IssmDouble z1=xyz_list[3*0+2];
    116         IssmDouble x2=xyz_list[3*1+0];
    117         IssmDouble y2=xyz_list[3*1+1];
    118         IssmDouble z2=xyz_list[3*1+2];
    119         IssmDouble x3=xyz_list[3*2+0];
    120         IssmDouble y3=xyz_list[3*2+1];
    121         IssmDouble z3=xyz_list[3*2+2];
    122 
    123         /*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
    124         *Jdet=SQRT3/6.*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2),0.5);
    125         if(*Jdet<0) _error_("negative jacobian determinant!");
    126 }
    127 /*}}}*/
    128 void PentaRef::GetSegmentJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    129         /*The Jacobian determinant is constant over the element, discard the gaussian points.
    130          * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    131 
    132         IssmDouble x1=xyz_list[3*0+0];
    133         IssmDouble y1=xyz_list[3*0+1];
    134         IssmDouble z1=xyz_list[3*0+2];
    135         IssmDouble x2=xyz_list[3*1+0];
    136         IssmDouble y2=xyz_list[3*1+1];
    137         IssmDouble z2=xyz_list[3*1+2];
    138 
    139         *Jdet=.5*sqrt(pow(x2-x1,2) + pow(y2-y1,2) + pow(z2-z1,2));
    140249        if(*Jdet<0) _error_("negative jacobian determinant!");
    141250
    … …  
    886995}
    887996/*}}}*/
    888 void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss,int finiteelement){/*{{{*/
    889         /* WARNING: For a significant gain in performance, it is better to use
    890          * static memory allocation instead of dynamic.*/
    891 
    892         /*Allocate basis functions*/
    893         IssmDouble  basis[NUMNODESMAX];
    894 
    895         /*Fetch number of nodes for this finite element*/
    896         int numnodes = this->NumberofNodes(finiteelement);
    897         _assert_(numnodes<=NUMNODESMAX);
    898 
    899         /*Get basis functions at this point*/
    900         GetNodalFunctions(&basis[0],gauss,finiteelement);
    901 
    902         /*Calculate parameter for this Gauss point*/
    903         IssmDouble value =0.;
    904         for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
    905 
    906         /*Assign output pointer*/
    907         *pvalue = value;
    908 }
    909 /*}}}*/
    910 void PentaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
    911         /*From node values of parameter p (p_list[0], p_list[1], p_list[2],
    912          * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
    913          * gaussian point specified by gauss_coord:
    914          *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
    915          *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
    916          *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
    917          *
    918          *   p is a vector of size 3x1 already allocated.
    919          *
    920          * WARNING: For a significant gain in performance, it is better to use
    921          * static memory allocation instead of dynamic.
    922          */
    923 
    924         /*Allocate derivatives of basis functions*/
    925         IssmDouble  dbasis[3*NUMNODESMAX];
    926 
    927         /*Fetch number of nodes for this finite element*/
    928         int numnodes = this->NumberofNodes(finiteelement);
    929         _assert_(numnodes<=NUMNODESMAX);
    930 
    931         /*Get basis functions derivatives at this point*/
    932         GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
    933 
    934         /*Calculate parameter for this Gauss point*/
    935         IssmDouble dpx=0.;
    936         IssmDouble dpy=0.;
    937         IssmDouble dpz=0.;
    938         for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
    939         for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
    940         for(int i=0;i<numnodes;i++) dpz += dbasis[2*numnodes+i]*plist[i];
    941 
    942         /*Assign values*/
    943         p[0]=dpx;
    944         p[1]=dpy;
    945         p[2]=dpz;
     997void PentaRef::GetSegmentJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     998        /*The Jacobian determinant is constant over the element, discard the gaussian points.
     999         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
     1000
     1001        IssmDouble x1=xyz_list[3*0+0];
     1002        IssmDouble y1=xyz_list[3*0+1];
     1003        IssmDouble z1=xyz_list[3*0+2];
     1004        IssmDouble x2=xyz_list[3*1+0];
     1005        IssmDouble y2=xyz_list[3*1+1];
     1006        IssmDouble z2=xyz_list[3*1+2];
     1007
     1008        *Jdet=.5*sqrt(pow(x2-x1,2) + pow(y2-y1,2) + pow(z2-z1,2));
     1009        if(*Jdet<0) _error_("negative jacobian determinant!");
     1010
     1011}
     1012/*}}}*/
     1013void PentaRef::GetTriaJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1014        /*The Jacobian determinant is constant over the element, discard the gaussian points.
     1015         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
     1016
     1017        IssmDouble x1=xyz_list[3*0+0];
     1018        IssmDouble y1=xyz_list[3*0+1];
     1019        IssmDouble z1=xyz_list[3*0+2];
     1020        IssmDouble x2=xyz_list[3*1+0];
     1021        IssmDouble y2=xyz_list[3*1+1];
     1022        IssmDouble z2=xyz_list[3*1+2];
     1023        IssmDouble x3=xyz_list[3*2+0];
     1024        IssmDouble y3=xyz_list[3*2+1];
     1025        IssmDouble z3=xyz_list[3*2+2];
     1026
     1027        /*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
     1028        *Jdet=SQRT3/6.*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2),0.5);
     1029        if(*Jdet<0) _error_("negative jacobian determinant!");
    9461030}
    9471031/*}}}*/
    … …  
    9771061}
    9781062/*}}}*/
     1063int  PentaRef::PressureInterpolation(int fe_stokes){/*{{{*/
     1064
     1065        switch(fe_stokes){
     1066                case P1P1Enum:              return P1Enum;
     1067                case P1P1GLSEnum:           return P1Enum;
     1068                case MINIcondensedEnum:     return P1Enum;
     1069                case MINIEnum:              return P1Enum;
     1070                case TaylorHoodEnum:        return P1Enum;
     1071                case LATaylorHoodEnum:      return NoneEnum;
     1072                case OneLayerP4zEnum:       return P1Enum;
     1073                case CrouzeixRaviartEnum:   return P1DGEnum;
     1074                case LACrouzeixRaviartEnum: return NoneEnum;
     1075                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
     1076        }
     1077
     1078        return -1;
     1079}
     1080/*}}}*/
     1081void PentaRef::SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement){/*{{{*/
     1082
     1083        /*Output*/
     1084        int  numindices;
     1085        int* indices = NULL;
     1086
     1087        switch(finiteelement){
     1088                case P1Enum: case P1DGEnum:
     1089                        numindices = 3;
     1090                        indices    = xNew<int>(numindices);
     1091                        indices[0] = 3;
     1092                        indices[1] = 4;
     1093                        indices[2] = 5;
     1094                        break;
     1095                case P1bubbleEnum: case P1bubblecondensedEnum:
     1096                        numindices = 3;
     1097                        indices    = xNew<int>(numindices);
     1098                        indices[0] = 3;
     1099                        indices[1] = 4;
     1100                        indices[2] = 5;
     1101                        break;
     1102                case P2xP1Enum:
     1103                        numindices = 6;
     1104                        indices    = xNew<int>(numindices);
     1105                        indices[0] = 3;
     1106                        indices[1] = 4;
     1107                        indices[2] = 5;
     1108                        indices[3] = 9;
     1109                        indices[4] = 10;
     1110                        indices[5] = 11;
     1111                        break;
     1112                case P1xP2Enum:
     1113                        numindices = 3;
     1114                        indices    = xNew<int>(numindices);
     1115                        indices[0] = 3;
     1116                        indices[1] = 4;
     1117                        indices[2] = 5;
     1118                        return;
     1119                case P2Enum:
     1120                        numindices = 6;
     1121                        indices    = xNew<int>(numindices);
     1122                        indices[0] = 3;
     1123                        indices[1] = 4;
     1124                        indices[2] = 5;
     1125                        indices[3] = 12;
     1126                        indices[4] = 13;
     1127                        indices[5] = 14;
     1128                        break;
     1129                default:
     1130                        _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
     1131        }
     1132
     1133        /*Assign output pointer*/
     1134        *pnumindices = numindices;
     1135        *pindices    = indices;
     1136}
     1137/*}}}*/
     1138int  PentaRef::TensorInterpolation(int fe_stokes){/*{{{*/
     1139
     1140        switch(fe_stokes){
     1141                case XTaylorHoodEnum:    return P1DGEnum;
     1142                default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
     1143        }
     1144
     1145        return -1;
     1146}
     1147/*}}}*/
    9791148int  PentaRef::VelocityInterpolation(int fe_stokes){/*{{{*/
    9801149
    … …  
    9951164}
    9961165/*}}}*/
    997 int  PentaRef::PressureInterpolation(int fe_stokes){/*{{{*/
    998 
    999         switch(fe_stokes){
    1000                 case P1P1Enum:              return P1Enum;
    1001                 case P1P1GLSEnum:           return P1Enum;
    1002                 case MINIcondensedEnum:     return P1Enum;
    1003                 case MINIEnum:              return P1Enum;
    1004                 case TaylorHoodEnum:        return P1Enum;
    1005                 case LATaylorHoodEnum:      return NoneEnum;
    1006                 case OneLayerP4zEnum:       return P1Enum;
    1007                 case CrouzeixRaviartEnum:   return P1DGEnum;
    1008                 case LACrouzeixRaviartEnum: return NoneEnum;
    1009                 default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    1010         }
    1011 
    1012         return -1;
    1013 }
    1014 /*}}}*/
    1015 int  PentaRef::TensorInterpolation(int fe_stokes){/*{{{*/
    1016 
    1017         switch(fe_stokes){
    1018                 case XTaylorHoodEnum:    return P1DGEnum;
    1019                 default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    1020         }
    1021 
    1022         return -1;
    1023 }
    1024 /*}}}*/
    1025 void PentaRef::BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){/*{{{*/
    1026 
    1027         /*Output*/
    1028         int  numindices;
    1029         int* indices = NULL;
    1030 
    1031         switch(finiteelement){
    1032                 case P1Enum: case P1DGEnum:
    1033                         numindices = 3;
    1034                         indices    = xNew<int>(numindices);
    1035                         indices[0] = 0;
    1036                         indices[1] = 1;
    1037                         indices[2] = 2;
    1038                         break;
    1039                 case P1bubbleEnum: case P1bubblecondensedEnum:
    1040                         numindices = 3;
    1041                         indices    = xNew<int>(numindices);
    1042                         indices[0] = 0;
    1043                         indices[1] = 1;
    1044                         indices[2] = 2;
    1045                         break;
    1046                 case P2xP1Enum:
    1047                         numindices = 6;
    1048                         indices    = xNew<int>(numindices);
    1049                         indices[0] = 0;
    1050                         indices[1] = 1;
    1051                         indices[2] = 2;
    1052                         indices[3] = 6;
    1053                         indices[4] = 7;
    1054                         indices[5] = 8;
    1055                         break;
    1056                 case P1xP2Enum:
    1057                         numindices = 3;
    1058                         indices    = xNew<int>(numindices);
    1059                         indices[0] = 0;
    1060                         indices[1] = 1;
    1061                         indices[2] = 2;
    1062                         break;
    1063                 case P1xP3Enum:
    1064                         numindices = 3;
    1065                         indices    = xNew<int>(numindices);
    1066                         indices[0] = 0;
    1067                         indices[1] = 1;
    1068                         indices[2] = 2;
    1069                         break;
    1070                 case P2Enum:
    1071                         numindices = 6;
    1072                         indices    = xNew<int>(numindices);
    1073                         indices[0] = 0;
    1074                         indices[1] = 1;
    1075                         indices[2] = 2;
    1076                         indices[3] = 9;
    1077                         indices[4] = 10;
    1078                         indices[5] = 11;
    1079                         break;
    1080                 case P2bubbleEnum:
    1081                         numindices = 6;
    1082                         indices    = xNew<int>(numindices);
    1083                         indices[0] = 0;
    1084                         indices[1] = 1;
    1085                         indices[2] = 2;
    1086                         indices[3] = 9;
    1087                         indices[4] = 10;
    1088                         indices[5] = 11;
    1089                         break;
    1090                 case P2xP4Enum:
    1091                         numindices = 6;
    1092                         indices    = xNew<int>(numindices);
    1093                         indices[0] = 0;
    1094                         indices[1] = 1;
    1095                         indices[2] = 2;
    1096                         indices[3] = 9;
    1097                         indices[4] = 10;
    1098                         indices[5] = 11;
    1099                         break;
    1100                 default:
    1101                         _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    1102         }
    1103 
    1104         /*Assign output pointer*/
    1105         *pnumindices = numindices;
    1106         *pindices    = indices;
    1107 }
    1108 /*}}}*/
    1109 void PentaRef::SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement){/*{{{*/
    1110 
    1111         /*Output*/
    1112         int  numindices;
    1113         int* indices = NULL;
    1114 
    1115         switch(finiteelement){
    1116                 case P1Enum: case P1DGEnum:
    1117                         numindices = 3;
    1118                         indices    = xNew<int>(numindices);
    1119                         indices[0] = 3;
    1120                         indices[1] = 4;
    1121                         indices[2] = 5;
    1122                         break;
    1123                 case P1bubbleEnum: case P1bubblecondensedEnum:
    1124                         numindices = 3;
    1125                         indices    = xNew<int>(numindices);
    1126                         indices[0] = 3;
    1127                         indices[1] = 4;
    1128                         indices[2] = 5;
    1129                         break;
    1130                 case P2xP1Enum:
    1131                         numindices = 6;
    1132                         indices    = xNew<int>(numindices);
    1133                         indices[0] = 3;
    1134                         indices[1] = 4;
    1135                         indices[2] = 5;
    1136                         indices[3] = 9;
    1137                         indices[4] = 10;
    1138                         indices[5] = 11;
    1139                         break;
    1140                 case P1xP2Enum:
    1141                         numindices = 3;
    1142                         indices    = xNew<int>(numindices);
    1143                         indices[0] = 3;
    1144                         indices[1] = 4;
    1145                         indices[2] = 5;
    1146                         return;
    1147                 case P2Enum:
    1148                         numindices = 6;
    1149                         indices    = xNew<int>(numindices);
    1150                         indices[0] = 3;
    1151                         indices[1] = 4;
    1152                         indices[2] = 5;
    1153                         indices[3] = 12;
    1154                         indices[4] = 13;
    1155                         indices[5] = 14;
    1156                         break;
    1157                 default:
    1158                         _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    1159         }
    1160 
    1161         /*Assign output pointer*/
    1162         *pnumindices = numindices;
    1163         *pindices    = indices;
    1164 }
    1165 /*}}}*/

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

    r18078 r18925  
    1515
    1616                /*Numerics*/
     17                void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement);
     18                void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
     19                void GetInputValue(IssmDouble* pvalue,IssmDouble* plist, Gauss* gauss,int finiteelement);
     20                void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss);
     21                void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
     22                void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,Gauss* gauss);
     23                void GetLprimeFSSSA(IssmDouble* LprimeFSSSA, IssmDouble* xyz_list, Gauss* gauss);
    1724                void GetNodalFunctions(IssmDouble* basis, Gauss* gauss,int finiteelement);
    1825                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss,int finiteelement);
    1926                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss,int finiteelement);
    2027                void GetQuadJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
    21                 void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss);
    22                 void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
     28                void GetSegmentJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
    2329                void GetTriaJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
    24                 void GetSegmentJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
    25                 void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,Gauss* gauss);
    26                 void GetLprimeFSSSA(IssmDouble* LprimeFSSSA, IssmDouble* xyz_list, Gauss* gauss);
    27                 void GetInputValue(IssmDouble* pvalue,IssmDouble* plist, Gauss* gauss,int finiteelement);
    28                 void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
    29 
    30                 void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement);
     30                int  NumberofNodes(int finiteelement);
     31                int  PressureInterpolation(int fe_stokes);
    3132                void SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement);
    32                 int  NumberofNodes(int finiteelement);
     33                int  TensorInterpolation(int fe_stokes);
    3334                int  VelocityInterpolation(int fe_stokes);
    34                 int  PressureInterpolation(int fe_stokes);
    35                 int  TensorInterpolation(int fe_stokes);
    3635};
    3736#endif

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

    r18523 r18925  
    2929
    3030/*Reference Element numerics*/
     31void SegRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement){/*{{{*/
     32
     33        /*From node values of parameter p (plist[0],plist[1]), return parameter derivative value at gaussian
     34         * point specified by gauss_basis:
     35         *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx
     36         *
     37         * p is a vector already allocated.
     38         *
     39         * WARNING: For a significant gain in performance, it is better to use
     40         * static memory allocation instead of dynamic.
     41         */
     42
     43        /*Allocate derivatives of basis functions*/
     44        IssmDouble  dbasis[1*NUMNODESMAX];
     45
     46        /*Fetch number of nodes for this finite element*/
     47        int numnodes = this->NumberofNodes(finiteelement);
     48        _assert_(numnodes<=NUMNODESMAX);
     49
     50        /*Get basis functions derivatives at this point*/
     51        GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
     52
     53        /*Calculate parameter for this Gauss point*/
     54        IssmDouble dpx=0.;
     55        for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
     56
     57        /*Assign values*/
     58        p[0]=dpx;
     59}
     60/*}}}*/
     61void SegRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss,int finiteelement){/*{{{*/
     62        /* WARNING: For a significant gain in performance, it is better to use
     63         * static memory allocation instead of dynamic.*/
     64
     65        /*Allocate basis functions*/
     66        IssmDouble  basis[NUMNODESMAX];
     67
     68        /*Fetch number of nodes for this finite element*/
     69        int numnodes = this->NumberofNodes(finiteelement);
     70        _assert_(numnodes<=NUMNODESMAX);
     71
     72        /*Get basis functions at this point*/
     73        GetNodalFunctions(&basis[0],gauss,finiteelement);
     74
     75        /*Calculate parameter for this Gauss point*/
     76        IssmDouble value =0.;
     77        for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
     78
     79        /*Assign output pointer*/
     80        *p = value;
     81}
     82/*}}}*/
     83void SegRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussSeg* gauss){/*{{{*/
     84        /*The Jacobian is constant over the element, discard the gaussian points.
     85         * J is assumed to have been allocated of size 1*/
     86
     87        IssmDouble x1=xyz_list[3*0+0];
     88        IssmDouble x2=xyz_list[3*1+0];
     89
     90        *J=.5*fabs(x2-x1);
     91}
     92/*}}}*/
     93void SegRef::GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussSeg* gauss){/*{{{*/
     94        /*The Jacobian determinant is constant over the element, discard the gaussian points.
     95         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
     96
     97        /*Call Jacobian routine to get the jacobian:*/
     98        GetJacobian(Jdet, xyz_list, gauss);
     99        if(*Jdet<0) _error_("negative jacobian determinant!");
     100
     101}
     102/*}}}*/
     103void SegRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussSeg* gauss){/*{{{*/
     104
     105        /*Jacobian*/
     106        IssmDouble J;
     107
     108        /*Call Jacobian routine to get the jacobian:*/
     109        GetJacobian(&J, xyz_list, gauss);
     110
     111        /*Invert Jacobian matrix: */
     112        *Jinv = 1./J;
     113}
     114/*}}}*/
    31115void SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss,int finiteelement){/*{{{*/
    32116        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    … …  
    109193}
    110194/*}}}*/
    111 void SegRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement){/*{{{*/
    112 
    113         /*From node values of parameter p (plist[0],plist[1]), return parameter derivative value at gaussian
    114          * point specified by gauss_basis:
    115          *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx
    116          *
    117          * p is a vector already allocated.
    118          *
    119          * WARNING: For a significant gain in performance, it is better to use
    120          * static memory allocation instead of dynamic.
    121          */
    122 
    123         /*Allocate derivatives of basis functions*/
    124         IssmDouble  dbasis[1*NUMNODESMAX];
    125 
    126         /*Fetch number of nodes for this finite element*/
    127         int numnodes = this->NumberofNodes(finiteelement);
    128         _assert_(numnodes<=NUMNODESMAX);
    129 
    130         /*Get basis functions derivatives at this point*/
    131         GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
    132 
    133         /*Calculate parameter for this Gauss point*/
    134         IssmDouble dpx=0.;
    135         for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
    136 
    137         /*Assign values*/
    138         p[0]=dpx;
    139 }
    140 /*}}}*/
    141 void SegRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss,int finiteelement){/*{{{*/
    142         /* WARNING: For a significant gain in performance, it is better to use
    143          * static memory allocation instead of dynamic.*/
    144 
    145         /*Allocate basis functions*/
    146         IssmDouble  basis[NUMNODESMAX];
    147 
    148         /*Fetch number of nodes for this finite element*/
    149         int numnodes = this->NumberofNodes(finiteelement);
    150         _assert_(numnodes<=NUMNODESMAX);
    151 
    152         /*Get basis functions at this point*/
    153         GetNodalFunctions(&basis[0],gauss,finiteelement);
    154 
    155         /*Calculate parameter for this Gauss point*/
    156         IssmDouble value =0.;
    157         for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
    158 
    159         /*Assign output pointer*/
    160         *p = value;
    161 }
    162 /*}}}*/
    163 void SegRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussSeg* gauss){/*{{{*/
    164         /*The Jacobian is constant over the element, discard the gaussian points.
    165          * J is assumed to have been allocated of size 1*/
    166 
    167         IssmDouble x1=xyz_list[3*0+0];
    168         IssmDouble x2=xyz_list[3*1+0];
    169 
    170         *J=.5*fabs(x2-x1);
    171 }
    172 /*}}}*/
    173 void SegRef::GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussSeg* gauss){/*{{{*/
    174         /*The Jacobian determinant is constant over the element, discard the gaussian points.
    175          * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    176 
    177         /*Call Jacobian routine to get the jacobian:*/
    178         GetJacobian(Jdet, xyz_list, gauss);
    179         if(*Jdet<0) _error_("negative jacobian determinant!");
    180 
    181 }
    182 /*}}}*/
    183 void SegRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussSeg* gauss){/*{{{*/
    184 
    185         /*Jacobian*/
    186         IssmDouble J;
    187 
    188         /*Call Jacobian routine to get the jacobian:*/
    189         GetJacobian(&J, xyz_list, gauss);
    190 
    191         /*Invert Jacobian matrix: */
    192         *Jinv = 1./J;
    193 }
    194 /*}}}*/
    195195int  SegRef::NumberofNodes(int finiteelement){/*{{{*/
    196196

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

    r18078 r18925  
    1616                ~SegRef();
    1717
     18                void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement);
     19                void GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss,int finiteelement);
    1820                void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussSeg* gauss);
    1921                void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussSeg* gauss);
    … …  
    2224                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement);
    2325                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss,int finiteelement);
    24                 void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement);
    25                 void GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss,int finiteelement);
    2626                int  NumberofNodes(int finiteelement);
    2727};

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

    r18523 r18925  
    3131
    3232/*Reference Element numerics*/
     33void TetraRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement){/*{{{*/
     34        /*From node values of parameter p (p_list[0], p_list[1], p_list[2],
     35         * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
     36         * gaussian point specified by gauss_coord:
     37         *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
     38         *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
     39         *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
     40         *
     41         *   p is a vector of size 3x1 already allocated.
     42         *
     43         * WARNING: For a significant gain in performance, it is better to use
     44         * static memory allocation instead of dynamic.
     45         */
     46
     47        /*Allocate derivatives of basis functions*/
     48        IssmDouble  dbasis[3*NUMNODESMAX];
     49
     50        /*Fetch number of nodes for this finite element*/
     51        int numnodes = this->NumberofNodes(finiteelement);
     52        _assert_(numnodes<=NUMNODESMAX);
     53
     54        /*Get basis functions derivatives at this point*/
     55        GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
     56
     57        /*Calculate parameter for this Gauss point*/
     58        IssmDouble dpx=0.;
     59        IssmDouble dpy=0.;
     60        IssmDouble dpz=0.;
     61        for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
     62        for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
     63        for(int i=0;i<numnodes;i++) dpz += dbasis[2*numnodes+i]*plist[i];
     64
     65        /*Assign values*/
     66        p[0]=dpx;
     67        p[1]=dpy;
     68        p[2]=dpz;
     69}
     70/*}}}*/
     71void TetraRef::GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss,int finiteelement){/*{{{*/
     72        /* WARNING: For a significant gain in performance, it is better to use
     73         * static memory allocation instead of dynamic.*/
     74
     75        /*Allocate basis functions*/
     76        IssmDouble  basis[NUMNODESMAX];
     77
     78        /*Fetch number of nodes for this finite element*/
     79        int numnodes = this->NumberofNodes(finiteelement);
     80        _assert_(numnodes<=NUMNODESMAX);
     81
     82        /*Get basis functions at this point*/
     83        GetNodalFunctions(&basis[0],gauss,finiteelement);
     84
     85        /*Calculate parameter for this Gauss point*/
     86        IssmDouble value =0.;
     87        for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
     88
     89        /*Assign output pointer*/
     90        *p = value;
     91}
     92/*}}}*/
     93void TetraRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
     94        /*The Jacobian is constant over the element, discard the gaussian points.
     95         * J is assumed to have been allocated of size 1*/
     96
     97        IssmDouble x1=xyz_list[3*0+0];
     98        IssmDouble x2=xyz_list[3*1+0];
     99        IssmDouble x3=xyz_list[3*2+0];
     100        IssmDouble x4=xyz_list[3*3+0];
     101
     102        IssmDouble y1=xyz_list[3*0+1];
     103        IssmDouble y2=xyz_list[3*1+1];
     104        IssmDouble y3=xyz_list[3*2+1];
     105        IssmDouble y4=xyz_list[3*3+1];
     106
     107        IssmDouble z1=xyz_list[3*0+2];
     108        IssmDouble z2=xyz_list[3*1+2];
     109        IssmDouble z3=xyz_list[3*2+2];
     110        IssmDouble z4=xyz_list[3*3+2];
     111
     112        J[NDOF3*0+0] = x2-x1;
     113        J[NDOF3*0+1] = y2-y1;
     114        J[NDOF3*0+2] = z2-z1;
     115
     116        J[NDOF3*1+0] = x3-x1;
     117        J[NDOF3*1+1] = y3-y1;
     118        J[NDOF3*1+2] = z3-z1;
     119
     120        J[NDOF3*2+0] = x4-x1;
     121        J[NDOF3*2+1] = y4-y1;
     122        J[NDOF3*2+2] = z4-z1;
     123}
     124/*}}}*/
     125void TetraRef::GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
     126        /*The Jacobian determinant is constant over the element, discard the gaussian points.
     127         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
     128        IssmDouble J[3][3];
     129
     130        /*Call Jacobian routine to get the jacobian:*/
     131        GetJacobian(&J[0][0],xyz_list, gauss);
     132
     133        /*Get Determinant*/
     134        Matrix3x3Determinant(Jdet,&J[0][0]);
     135        if(*Jdet<0) _error_("negative jacobian determinant!");
     136
     137}
     138/*}}}*/
     139void TetraRef::GetJacobianDeterminantFace(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
     140        /*The Jacobian determinant is constant over the element, discard the gaussian points.
     141         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
     142
     143        IssmDouble x1=xyz_list[3*0+0];
     144        IssmDouble y1=xyz_list[3*0+1];
     145        IssmDouble z1=xyz_list[3*0+2];
     146        IssmDouble x2=xyz_list[3*1+0];
     147        IssmDouble y2=xyz_list[3*1+1];
     148        IssmDouble z2=xyz_list[3*1+2];
     149        IssmDouble x3=xyz_list[3*2+0];
     150        IssmDouble y3=xyz_list[3*2+1];
     151        IssmDouble z3=xyz_list[3*2+2];
     152
     153        /*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
     154        *Jdet=SQRT3/6.*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2),0.5);
     155        if(*Jdet<0) _error_("negative jacobian determinant!");
     156}
     157/*}}}*/
     158void TetraRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
     159
     160        /*Jacobian*/
     161        IssmDouble J[3][3];
     162
     163        /*Call Jacobian routine to get the jacobian:*/
     164        GetJacobian(&J[0][0], xyz_list, gauss);
     165
     166        /*Invert Jacobian matrix: */
     167        Matrix3x3Invert(Jinv,&J[0][0]);
     168}
     169/*}}}*/
    33170void TetraRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss_in,int finiteelement){/*{{{*/
    34171        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    … …  
    207344        }
    208345
    209 }
    210 /*}}}*/
    211 void TetraRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement){/*{{{*/
    212         /*From node values of parameter p (p_list[0], p_list[1], p_list[2],
    213          * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
    214          * gaussian point specified by gauss_coord:
    215          *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
    216          *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
    217          *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
    218          *
    219          *   p is a vector of size 3x1 already allocated.
    220          *
    221          * WARNING: For a significant gain in performance, it is better to use
    222          * static memory allocation instead of dynamic.
    223          */
    224 
    225         /*Allocate derivatives of basis functions*/
    226         IssmDouble  dbasis[3*NUMNODESMAX];
    227 
    228         /*Fetch number of nodes for this finite element*/
    229         int numnodes = this->NumberofNodes(finiteelement);
    230         _assert_(numnodes<=NUMNODESMAX);
    231 
    232         /*Get basis functions derivatives at this point*/
    233         GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
    234 
    235         /*Calculate parameter for this Gauss point*/
    236         IssmDouble dpx=0.;
    237         IssmDouble dpy=0.;
    238         IssmDouble dpz=0.;
    239         for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
    240         for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
    241         for(int i=0;i<numnodes;i++) dpz += dbasis[2*numnodes+i]*plist[i];
    242 
    243         /*Assign values*/
    244         p[0]=dpx;
    245         p[1]=dpy;
    246         p[2]=dpz;
    247 }
    248 /*}}}*/
    249 void TetraRef::GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss,int finiteelement){/*{{{*/
    250         /* WARNING: For a significant gain in performance, it is better to use
    251          * static memory allocation instead of dynamic.*/
    252 
    253         /*Allocate basis functions*/
    254         IssmDouble  basis[NUMNODESMAX];
    255 
    256         /*Fetch number of nodes for this finite element*/
    257         int numnodes = this->NumberofNodes(finiteelement);
    258         _assert_(numnodes<=NUMNODESMAX);
    259 
    260         /*Get basis functions at this point*/
    261         GetNodalFunctions(&basis[0],gauss,finiteelement);
    262 
    263         /*Calculate parameter for this Gauss point*/
    264         IssmDouble value =0.;
    265         for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
    266 
    267         /*Assign output pointer*/
    268         *p = value;
    269 }
    270 /*}}}*/
    271 void TetraRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
    272         /*The Jacobian is constant over the element, discard the gaussian points.
    273          * J is assumed to have been allocated of size 1*/
    274 
    275         IssmDouble x1=xyz_list[3*0+0];
    276         IssmDouble x2=xyz_list[3*1+0];
    277         IssmDouble x3=xyz_list[3*2+0];
    278         IssmDouble x4=xyz_list[3*3+0];
    279 
    280         IssmDouble y1=xyz_list[3*0+1];
    281         IssmDouble y2=xyz_list[3*1+1];
    282         IssmDouble y3=xyz_list[3*2+1];
    283         IssmDouble y4=xyz_list[3*3+1];
    284 
    285         IssmDouble z1=xyz_list[3*0+2];
    286         IssmDouble z2=xyz_list[3*1+2];
    287         IssmDouble z3=xyz_list[3*2+2];
    288         IssmDouble z4=xyz_list[3*3+2];
    289 
    290         J[NDOF3*0+0] = x2-x1;
    291         J[NDOF3*0+1] = y2-y1;
    292         J[NDOF3*0+2] = z2-z1;
    293 
    294         J[NDOF3*1+0] = x3-x1;
    295         J[NDOF3*1+1] = y3-y1;
    296         J[NDOF3*1+2] = z3-z1;
    297 
    298         J[NDOF3*2+0] = x4-x1;
    299         J[NDOF3*2+1] = y4-y1;
    300         J[NDOF3*2+2] = z4-z1;
    301 }
    302 /*}}}*/
    303 void TetraRef::GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
    304         /*The Jacobian determinant is constant over the element, discard the gaussian points.
    305          * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    306         IssmDouble J[3][3];
    307 
    308         /*Call Jacobian routine to get the jacobian:*/
    309         GetJacobian(&J[0][0],xyz_list, gauss);
    310 
    311         /*Get Determinant*/
    312         Matrix3x3Determinant(Jdet,&J[0][0]);
    313         if(*Jdet<0) _error_("negative jacobian determinant!");
    314 
    315 }
    316 /*}}}*/
    317 void TetraRef::GetJacobianDeterminantFace(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
    318         /*The Jacobian determinant is constant over the element, discard the gaussian points.
    319          * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    320 
    321         IssmDouble x1=xyz_list[3*0+0];
    322         IssmDouble y1=xyz_list[3*0+1];
    323         IssmDouble z1=xyz_list[3*0+2];
    324         IssmDouble x2=xyz_list[3*1+0];
    325         IssmDouble y2=xyz_list[3*1+1];
    326         IssmDouble z2=xyz_list[3*1+2];
    327         IssmDouble x3=xyz_list[3*2+0];
    328         IssmDouble y3=xyz_list[3*2+1];
    329         IssmDouble z3=xyz_list[3*2+2];
    330 
    331         /*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
    332         *Jdet=SQRT3/6.*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2),0.5);
    333         if(*Jdet<0) _error_("negative jacobian determinant!");
    334 }
    335 /*}}}*/
    336 void TetraRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
    337 
    338         /*Jacobian*/
    339         IssmDouble J[3][3];
    340 
    341         /*Call Jacobian routine to get the jacobian:*/
    342         GetJacobian(&J[0][0], xyz_list, gauss);
    343 
    344         /*Invert Jacobian matrix: */
    345         Matrix3x3Invert(Jinv,&J[0][0]);
    346346}
    347347/*}}}*/
    … …  
    368368}
    369369/*}}}*/
     370int  TetraRef::PressureInterpolation(int fe_stokes){/*{{{*/
     371
     372        switch(fe_stokes){
     373                case P1P1Enum:          return P1Enum;
     374                case P1P1GLSEnum:       return P1Enum;
     375                case MINIcondensedEnum: return P1Enum;
     376                case MINIEnum:          return P1Enum;
     377                case TaylorHoodEnum:    return P1Enum;
     378                case LATaylorHoodEnum:  return NoneEnum;
     379                case XTaylorHoodEnum:   return P1Enum;
     380                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
     381        }
     382
     383        return -1;
     384}/*}}}*/
     385int  TetraRef::TensorInterpolation(int fe_stokes){/*{{{*/
     386        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     387
     388        switch(fe_stokes){
     389                case XTaylorHoodEnum: return P1DGEnum;
     390                default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
     391        }
     392}
     393/*}}}*/
    370394int  TetraRef::VelocityInterpolation(int fe_stokes){/*{{{*/
    371395
    … …  
    384408}
    385409/*}}}*/
    386 int  TetraRef::PressureInterpolation(int fe_stokes){/*{{{*/
    387 
    388         switch(fe_stokes){
    389                 case P1P1Enum:          return P1Enum;
    390                 case P1P1GLSEnum:       return P1Enum;
    391                 case MINIcondensedEnum: return P1Enum;
    392                 case MINIEnum:          return P1Enum;
    393                 case TaylorHoodEnum:    return P1Enum;
    394                 case LATaylorHoodEnum:  return NoneEnum;
    395                 case XTaylorHoodEnum:   return P1Enum;
    396                 default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    397         }
    398 
    399         return -1;
    400 }/*}}}*/
    401 int  TetraRef::TensorInterpolation(int fe_stokes){/*{{{*/
    402         /*This routine returns the values of the nodal functions  at the gaussian point.*/
    403 
    404         switch(fe_stokes){
    405                 case XTaylorHoodEnum: return P1DGEnum;
    406                 default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    407         }
    408 }
    409 /*}}}*/

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

    r18139 r18925  
    1616                ~TetraRef();
    1717
     18                void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement);
     19                void GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss,int finiteelement);
    1820                void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTetra* gauss);
    1921                void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss);
    … …  
    2325                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement);
    2426                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss,int finiteelement);
    25                 void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement);
    26                 void GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss,int finiteelement);
    27 
    2827                int  NumberofNodes(int finiteelement);
    29                 int  VelocityInterpolation(int fe_stokes);
    3028                int  PressureInterpolation(int fe_stokes);
    3129                int  TensorInterpolation(int fe_stokes);
     30                int  VelocityInterpolation(int fe_stokes);
    3231};
    3332#endif

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

    r18523 r18925  
    3232
    3333/*Reference Element numerics*/
    34 void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
    35         /*Compute B  matrix. B=[phi1 phi2 -phi3 -phi4]
     34void TriaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
     35        /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian
     36         * point specified by gauss_basis:
     37         *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
     38         *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
    3639         *
    37          * and phi1=phi3 phi2=phi4
     40         * p is a vector already allocated.
    3841         *
    39          * We assume B has been allocated already, of size: 1x4
     42         * WARNING: For a significant gain in performance, it is better to use
     43         * static memory allocation instead of dynamic.
    4044         */
     45
     46        /*Allocate derivatives of basis functions*/
     47        IssmDouble  dbasis[2*NUMNODESMAX];
    4148
    4249        /*Fetch number of nodes for this finite element*/
    4350        int numnodes = this->NumberofNodes(finiteelement);
    44 
    45         /*Get nodal functions*/
    46         IssmDouble* basis=xNew<IssmDouble>(numnodes);
    47         GetNodalFunctions(basis,gauss,finiteelement);
    48 
    49         /*Build B for this segment*/
    50         B[0] = +basis[index1];
    51         B[1] = +basis[index2];
    52         B[2] = -basis[index1];
    53         B[3] = -basis[index2];
    54 
    55         /*Clean-up*/
    56         xDelete<IssmDouble>(basis);
    57 }
    58 /*}}}*/
    59 void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
    60         /*Compute Bprime  matrix. Bprime=[phi1 phi2 phi3 phi4]
    61          *
    62          * and phi1=phi3 phi2=phi4
    63          *
    64          * We assume Bprime has been allocated already, of size: 1x4
    65          */
     51        _assert_(numnodes<=NUMNODESMAX);
     52
     53        /*Get basis functions derivatives at this point*/
     54        GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
     55
     56        /*Calculate parameter for this Gauss point*/
     57        IssmDouble dpx=0.;
     58        IssmDouble dpy=0.;
     59        for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
     60        for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
     61
     62        /*Assign values*/
     63        p[0]=dpx;
     64        p[1]=dpy;
     65
     66}
     67/*}}}*/
     68void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss,int finiteelement){/*{{{*/
     69        /* WARNING: For a significant gain in performance, it is better to use
     70         * static memory allocation instead of dynamic.*/
     71
     72        /*Allocate basis functions*/
     73        IssmDouble  basis[NUMNODESMAX];
    6674
    6775        /*Fetch number of nodes for this finite element*/
    6876        int numnodes = this->NumberofNodes(finiteelement);
    69 
    70         /*Get nodal functions*/
    71         IssmDouble* basis=xNew<IssmDouble>(numnodes);
    72         GetNodalFunctions(basis,gauss,finiteelement);
    73 
    74         /*Build B'*/
    75         Bprime[0] = basis[index1];
    76         Bprime[1] = basis[index2];
    77         Bprime[2] = basis[index1];
    78         Bprime[3] = basis[index2];
    79 
    80         /*Clean-up*/
    81         xDelete<IssmDouble>(basis);
     77        _assert_(numnodes<=NUMNODESMAX);
     78
     79        /*Get basis functions at this point*/
     80        GetNodalFunctions(&basis[0],gauss,finiteelement);
     81
     82        /*Calculate parameter for this Gauss point*/
     83        IssmDouble value =0.;
     84        for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
     85
     86        /*Assign output pointer*/
     87        *p = value;
    8288}
    8389/*}}}*/
    … …  
    97103        J[2*0+1] = 0.5*(y2-y1);
    98104        J[2*1+1] = SQRT3/6.0*(2*y3-y1-y2);
    99 }
    100 /*}}}*/
    101 void TriaRef::GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    102         /*The Jacobian determinant is constant over the element, discard the gaussian points.
    103          * J is assumed to have been allocated*/
    104 
    105         IssmDouble x1 = xyz_list[3*0+0];
    106         IssmDouble y1 = xyz_list[3*0+1];
    107         IssmDouble x2 = xyz_list[3*1+0];
    108         IssmDouble y2 = xyz_list[3*1+1];
    109 
    110         *Jdet = .5*sqrt(pow(x2-x1,2) + pow(y2-y1,2));
    111         if(*Jdet<0) _error_("negative jacobian determinant!");
    112 
    113105}
    114106/*}}}*/
    … …  
    193185                        _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    194186        }
    195 }
    196 /*}}}*/
    197 void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss,int index1,int index2,int finiteelement){/*{{{*/
    198         /*This routine returns the values of the nodal functions  at the gaussian point.*/
    199 
    200         _assert_(index1>=0 && index1<3);
    201         _assert_(index2>=0 && index2<3);
    202 
    203         /*Fetch number of nodes for this finite element*/
    204         int numnodes = this->NumberofNodes(finiteelement);
    205 
    206         /*Get nodal functions*/
    207         IssmDouble* triabasis=xNew<IssmDouble>(numnodes);
    208         GetNodalFunctions(triabasis,gauss,finiteelement);
    209 
    210         switch(finiteelement){
    211                 case P1Enum: case P1DGEnum:
    212                         basis[0]=triabasis[index1];
    213                         basis[1]=triabasis[index2];
    214                         xDelete<IssmDouble>(triabasis);
    215                         return;
    216                 case P1bubbleEnum: case P1bubblecondensedEnum:
    217                         basis[0]=triabasis[index1];
    218                         basis[1]=triabasis[index2];
    219                         xDelete<IssmDouble>(triabasis);
    220                         return;
    221                 case P2Enum:
    222                         _assert_(index2<index1);
    223                         basis[0]=triabasis[index1];
    224                         basis[1]=triabasis[index2];
    225                         basis[2]=triabasis[3+index2-1];
    226                         xDelete<IssmDouble>(triabasis);
    227                         return;
    228                 default:
    229                         _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    230         }
    231 
    232         /*Clean up*/
    233         xDelete<IssmDouble>(triabasis);
    234187}
    235188/*}}}*/
    … …  
    354307}
    355308/*}}}*/
    356 void TriaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
    357         /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian
    358          * point specified by gauss_basis:
    359          *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
    360          *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
     309void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
     310        /*Compute B  matrix. B=[phi1 phi2 -phi3 -phi4]
    361311         *
    362          * p is a vector already allocated.
     312         * and phi1=phi3 phi2=phi4
    363313         *
    364          * WARNING: For a significant gain in performance, it is better to use
    365          * static memory allocation instead of dynamic.
     314         * We assume B has been allocated already, of size: 1x4
    366315         */
    367 
    368         /*Allocate derivatives of basis functions*/
    369         IssmDouble  dbasis[2*NUMNODESMAX];
    370316
    371317        /*Fetch number of nodes for this finite element*/
    372318        int numnodes = this->NumberofNodes(finiteelement);
    373         _assert_(numnodes<=NUMNODESMAX);
    374 
    375         /*Get basis functions derivatives at this point*/
    376         GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
    377 
    378         /*Calculate parameter for this Gauss point*/
    379         IssmDouble dpx=0.;
    380         IssmDouble dpy=0.;
    381         for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
    382         for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
    383 
    384         /*Assign values*/
    385         p[0]=dpx;
    386         p[1]=dpy;
    387 
    388 }
    389 /*}}}*/
    390 void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss,int finiteelement){/*{{{*/
    391         /* WARNING: For a significant gain in performance, it is better to use
    392          * static memory allocation instead of dynamic.*/
    393 
    394         /*Allocate basis functions*/
    395         IssmDouble  basis[NUMNODESMAX];
     319
     320        /*Get nodal functions*/
     321        IssmDouble* basis=xNew<IssmDouble>(numnodes);
     322        GetNodalFunctions(basis,gauss,finiteelement);
     323
     324        /*Build B for this segment*/
     325        B[0] = +basis[index1];
     326        B[1] = +basis[index2];
     327        B[2] = -basis[index1];
     328        B[3] = -basis[index2];
     329
     330        /*Clean-up*/
     331        xDelete<IssmDouble>(basis);
     332}
     333/*}}}*/
     334void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
     335        /*Compute Bprime  matrix. Bprime=[phi1 phi2 phi3 phi4]
     336         *
     337         * and phi1=phi3 phi2=phi4
     338         *
     339         * We assume Bprime has been allocated already, of size: 1x4
     340         */
    396341
    397342        /*Fetch number of nodes for this finite element*/
    398343        int numnodes = this->NumberofNodes(finiteelement);
    399         _assert_(numnodes<=NUMNODESMAX);
    400 
    401         /*Get basis functions at this point*/
    402         GetNodalFunctions(&basis[0],gauss,finiteelement);
    403 
    404         /*Calculate parameter for this Gauss point*/
    405         IssmDouble value =0.;
    406         for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
    407 
    408         /*Assign output pointer*/
    409         *p = value;
    410 }
    411 /*}}}*/
    412 int  TriaRef::NumberofNodes(int finiteelement){/*{{{*/
     344
     345        /*Get nodal functions*/
     346        IssmDouble* basis=xNew<IssmDouble>(numnodes);
     347        GetNodalFunctions(basis,gauss,finiteelement);
     348
     349        /*Build B'*/
     350        Bprime[0] = basis[index1];
     351        Bprime[1] = basis[index2];
     352        Bprime[2] = basis[index1];
     353        Bprime[3] = basis[index2];
     354
     355        /*Clean-up*/
     356        xDelete<IssmDouble>(basis);
     357}
     358/*}}}*/
     359void TriaRef::GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     360        /*The Jacobian determinant is constant over the element, discard the gaussian points.
     361         * J is assumed to have been allocated*/
     362
     363        IssmDouble x1 = xyz_list[3*0+0];
     364        IssmDouble y1 = xyz_list[3*0+1];
     365        IssmDouble x2 = xyz_list[3*1+0];
     366        IssmDouble y2 = xyz_list[3*1+1];
     367
     368        *Jdet = .5*sqrt(pow(x2-x1,2) + pow(y2-y1,2));
     369        if(*Jdet<0) _error_("negative jacobian determinant!");
     370
     371}
     372/*}}}*/
     373void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss,int index1,int index2,int finiteelement){/*{{{*/
     374        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     375
     376        _assert_(index1>=0 && index1<3);
     377        _assert_(index2>=0 && index2<3);
     378
     379        /*Fetch number of nodes for this finite element*/
     380        int numnodes = this->NumberofNodes(finiteelement);
     381
     382        /*Get nodal functions*/
     383        IssmDouble* triabasis=xNew<IssmDouble>(numnodes);
     384        GetNodalFunctions(triabasis,gauss,finiteelement);
    413385
    414386        switch(finiteelement){
    415                 case NoneEnum:                return 0;
    416                 case P0Enum:                  return NUMNODESP0;
    417                 case P1Enum:                  return NUMNODESP1;
    418                 case P1DGEnum:                return NUMNODESP1;
    419                 case P1bubbleEnum:            return NUMNODESP1b;
    420                 case P1bubblecondensedEnum:   return NUMNODESP1b;
    421                 case P2Enum:                  return NUMNODESP2;
    422                 case P2bubbleEnum:            return NUMNODESP2b;
    423                 case P2bubblecondensedEnum:   return NUMNODESP2b;
    424                 case P1P1Enum:                return NUMNODESP1*2;
    425                 case P1P1GLSEnum:             return NUMNODESP1*2;
    426                 case MINIcondensedEnum:       return NUMNODESP1b+NUMNODESP1;
    427                 case MINIEnum:                return NUMNODESP1b+NUMNODESP1;
    428                 case TaylorHoodEnum:          return NUMNODESP2+NUMNODESP1;
    429                 case LATaylorHoodEnum:        return NUMNODESP2;
    430                 case XTaylorHoodEnum:         return NUMNODESP2+NUMNODESP1;
    431                 case CrouzeixRaviartEnum:     return NUMNODESP2b+NUMNODESP1;
    432                 case LACrouzeixRaviartEnum:   return NUMNODESP2b;
    433                 default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    434         }
    435 
    436         return -1;
    437 }
    438 /*}}}*/
    439 int  TriaRef::VelocityInterpolation(int fe_stokes){/*{{{*/
    440 
    441         switch(fe_stokes){
    442                 case P1P1Enum:              return P1Enum;
    443                 case P1P1GLSEnum:           return P1Enum;
    444                 case MINIcondensedEnum:     return P1bubbleEnum;
    445                 case MINIEnum:              return P1bubbleEnum;
    446                 case TaylorHoodEnum:        return P2Enum;
    447                 case LATaylorHoodEnum:      return P2Enum;
    448                 case XTaylorHoodEnum:       return P2Enum;
    449                 case CrouzeixRaviartEnum:   return P2bubbleEnum;
    450                 case LACrouzeixRaviartEnum: return P2bubbleEnum;
    451                 default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    452         }
    453 
    454         return -1;
    455 }
    456 /*}}}*/
    457 int  TriaRef::PressureInterpolation(int fe_stokes){/*{{{*/
    458 
    459         switch(fe_stokes){
    460                 case P1P1Enum:              return P1Enum;
    461                 case P1P1GLSEnum:           return P1Enum;
    462                 case MINIcondensedEnum:     return P1Enum;
    463                 case MINIEnum:              return P1Enum;
    464                 case TaylorHoodEnum:        return P1Enum;
    465                 case LATaylorHoodEnum:      return NoneEnum;
    466                 case XTaylorHoodEnum:       return P1Enum;
    467                 case CrouzeixRaviartEnum:   return P1DGEnum;
    468                 case LACrouzeixRaviartEnum: return NoneEnum;
    469                 default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    470         }
    471 
    472         return -1;
    473 }
    474 /*}}}*/
    475 int  TriaRef::TensorInterpolation(int fe_stokes){/*{{{*/
    476         /*This routine returns the values of the nodal functions  at the gaussian point.*/
    477 
    478         switch(fe_stokes){
    479                 case XTaylorHoodEnum: return P1DGEnum;
    480                 default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    481         }
     387                case P1Enum: case P1DGEnum:
     388                        basis[0]=triabasis[index1];
     389                        basis[1]=triabasis[index2];
     390                        xDelete<IssmDouble>(triabasis);
     391                        return;
     392                case P1bubbleEnum: case P1bubblecondensedEnum:
     393                        basis[0]=triabasis[index1];
     394                        basis[1]=triabasis[index2];
     395                        xDelete<IssmDouble>(triabasis);
     396                        return;
     397                case P2Enum:
     398                        _assert_(index2<index1);
     399                        basis[0]=triabasis[index1];
     400                        basis[1]=triabasis[index2];
     401                        basis[2]=triabasis[3+index2-1];
     402                        xDelete<IssmDouble>(triabasis);
     403                        return;
     404                default:
     405                        _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
     406        }
     407
     408        /*Clean up*/
     409        xDelete<IssmDouble>(triabasis);
    482410}
    483411/*}}}*/
    … …  
    541469}
    542470/*}}}*/
     471int  TriaRef::NumberofNodes(int finiteelement){/*{{{*/
     472
     473        switch(finiteelement){
     474                case NoneEnum:                return 0;
     475                case P0Enum:                  return NUMNODESP0;
     476                case P1Enum:                  return NUMNODESP1;
     477                case P1DGEnum:                return NUMNODESP1;
     478                case P1bubbleEnum:            return NUMNODESP1b;
     479                case P1bubblecondensedEnum:   return NUMNODESP1b;
     480                case P2Enum:                  return NUMNODESP2;
     481                case P2bubbleEnum:            return NUMNODESP2b;
     482                case P2bubblecondensedEnum:   return NUMNODESP2b;
     483                case P1P1Enum:                return NUMNODESP1*2;
     484                case P1P1GLSEnum:             return NUMNODESP1*2;
     485                case MINIcondensedEnum:       return NUMNODESP1b+NUMNODESP1;
     486                case MINIEnum:                return NUMNODESP1b+NUMNODESP1;
     487                case TaylorHoodEnum:          return NUMNODESP2+NUMNODESP1;
     488                case LATaylorHoodEnum:        return NUMNODESP2;
     489                case XTaylorHoodEnum:         return NUMNODESP2+NUMNODESP1;
     490                case CrouzeixRaviartEnum:     return NUMNODESP2b+NUMNODESP1;
     491                case LACrouzeixRaviartEnum:   return NUMNODESP2b;
     492                default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
     493        }
     494
     495        return -1;
     496}
     497/*}}}*/
     498int  TriaRef::PressureInterpolation(int fe_stokes){/*{{{*/
     499
     500        switch(fe_stokes){
     501                case P1P1Enum:              return P1Enum;
     502                case P1P1GLSEnum:           return P1Enum;
     503                case MINIcondensedEnum:     return P1Enum;
     504                case MINIEnum:              return P1Enum;
     505                case TaylorHoodEnum:        return P1Enum;
     506                case LATaylorHoodEnum:      return NoneEnum;
     507                case XTaylorHoodEnum:       return P1Enum;
     508                case CrouzeixRaviartEnum:   return P1DGEnum;
     509                case LACrouzeixRaviartEnum: return NoneEnum;
     510                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
     511        }
     512
     513        return -1;
     514}
     515/*}}}*/
     516int  TriaRef::TensorInterpolation(int fe_stokes){/*{{{*/
     517        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     518
     519        switch(fe_stokes){
     520                case XTaylorHoodEnum: return P1DGEnum;
     521                default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
     522        }
     523}
     524/*}}}*/
     525int  TriaRef::VelocityInterpolation(int fe_stokes){/*{{{*/
     526
     527        switch(fe_stokes){
     528                case P1P1Enum:              return P1Enum;
     529                case P1P1GLSEnum:           return P1Enum;
     530                case MINIcondensedEnum:     return P1bubbleEnum;
     531                case MINIEnum:              return P1bubbleEnum;
     532                case TaylorHoodEnum:        return P2Enum;
     533                case LATaylorHoodEnum:      return P2Enum;
     534                case XTaylorHoodEnum:       return P2Enum;
     535                case CrouzeixRaviartEnum:   return P2bubbleEnum;
     536                case LACrouzeixRaviartEnum: return P2bubbleEnum;
     537                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
     538        }
     539
     540        return -1;
     541}
     542/*}}}*/

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

    r18078 r18925  
    1616
    1717                /*Numerics*/
     18                void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
     19                void GetInputValue(IssmDouble* pp, IssmDouble* plist, Gauss* gauss,int finiteelement);
    1820                void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss);
    19                 void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss);
    2021                void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss);
    2122                void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,Gauss* gauss);
    2223                void GetNodalFunctions(IssmDouble* basis,Gauss* gauss,int finiteelement);
    23                 void GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss, int index1,int index2,int finiteelement);
     24                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
     25                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss,int finiteelement);
    2426                void GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement);
    2527                void GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement);
    26                 void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
    27                 void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss,int finiteelement);
    28                 void GetInputValue(IssmDouble* pp, IssmDouble* plist, Gauss* gauss,int finiteelement);
    29                 void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
    30 
     28                void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss);
     29                void GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss, int index1,int index2,int finiteelement);
    3130                void NodeOnEdgeIndices(int* pnumindices,int** pindices,int index,int finiteelement);
    3231                int  NumberofNodes(int finiteelement);
    33                 int  VelocityInterpolation(int fe_stokes);
    3432                int  PressureInterpolation(int fe_stokes);
    3533                int  TensorInterpolation(int fe_stokes);
     34                int  VelocityInterpolation(int fe_stokes);
    3635};
    3736#endif
Note: See TracChangeset for help on using the changeset viewer.