Changeset 18078


Ignore:
Timestamp:
05/31/14 23:49:20 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: moved element_type to element rather than reference element

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

Legend:

Unmodified
Added
Removed

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

    r18068 r18078  
    2626        this->inputs     = NULL;
    2727        this->parameters = NULL;
     28        this->element_type_list=NULL;
    2829}/*}}}*/
    2930Element::~Element(){/*{{{*/
     31        xDelete<int>(element_type_list);
    3032        delete inputs;
    3133}

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

    r18062 r18078  
    4646                Matpar      *matpar;
    4747                Parameters  *parameters;
     48
     49                int* element_type_list;
     50                int  element_type;
    4851
    4952        public:

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

    r18068 r18078  
    2525/*}}}*/
    2626Penta::Penta(int penta_id, int penta_sid, int index, IoModel* iomodel,int nummodels)/*{{{*/
    27         :PentaRef(nummodels)
    28         ,ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
     27        :ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
    2928
    3029        int penta_elements_ids[2];
    … …  
    3534
    3635        /*id: */
    37         this->id=penta_id;
    38         this->sid=penta_sid;
     36        this->id  = penta_id;
     37        this->sid = penta_sid;
    3938
    4039        /*Build neighbors list*/
    … …  
    5756        this->matpar            = NULL;
    5857        this->verticalneighbors = NULL;
     58
     59        /*Only allocate pointer*/
     60        this->element_type_list=xNew<int>(nummodels);
    5961}
    6062/*}}}*/
    … …  
    869871
    870872        _assert_(nodes);
    871         int numnodes = this->NumberofNodes();
     873        int numnodes = this->NumberofNodes(this->element_type);
    872874
    873875        for(int i=0;i<numnodes;i++){
    … …  
    879881/*}}}*/
    880882int        Penta::GetNumberOfNodes(void){/*{{{*/
    881         return this->NumberofNodes();
     883        return this->NumberofNodes(this->element_type);
    882884}
    883885/*}}}*/
    … …  
    888890Node* Penta::GetNode(int node_number){/*{{{*/
    889891        _assert_(node_number>=0);
    890         _assert_(node_number<this->NumberofNodes());
     892        _assert_(node_number<this->NumberofNodes(this->element_type));
    891893        return this->nodes[node_number];
    892894}
    … …  
    11691171                /*Step3: Vertically integrate A COPY of the original*/
    11701172                if(original_input->ObjectEnum()==PentaInputEnum){
    1171                         if(((PentaInput*)original_input)->element_type==P0Enum){
     1173                        if(((PentaInput*)original_input)->interpolation_type==P0Enum){
    11721174                                original_input->GetInputValue(&p0top1_list[i]);
    11731175                                element_integrated_input= new  PentaInput(original_input->InstanceEnum(),p0top1_list,P1Enum);
    … …  
    13581360
    13591361        /*Fetch number of nodes for this finite element*/
    1360         int numnodes = this->NumberofNodes();
     1362        int numnodes = this->NumberofNodes(this->element_type);
    13611363
    13621364        /*Fetch dof list and allocate solution vector*/
    … …  
    16451647
    16461648        _assert_(gauss->Enum()==GaussPentaEnum);
    1647         this->GetNodalFunctions(basis,(GaussPenta*)gauss);
     1649        this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->element_type);
    16481650
    16491651}
    … …  
    16521654
    16531655        _assert_(gauss->Enum()==GaussPentaEnum);
    1654         this->GetNodalFunctionsP1(basis,(GaussPenta*)gauss);
     1656        this->GetNodalFunctions(basis,(GaussPenta*)gauss,P1Enum);
    16551657
    16561658}
    … …  
    16591661
    16601662        _assert_(gauss->Enum()==GaussPentaEnum);
    1661         this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss);
     1663        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->element_type);
    16621664
    16631665}
    … …  
    16661668
    16671669        _assert_(gauss->Enum()==GaussPentaEnum);
    1668         this->GetNodalFunctionsP1Derivatives(dbasis,xyz_list,(GaussPenta*)gauss);
     1670        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1Enum);
    16691671
    16701672}
    … …  
    16731675
    16741676        _assert_(gauss->Enum()==GaussPentaEnum);
    1675         this->GetNodalFunctionsMINIDerivatives(dbasis,xyz_list,(GaussPenta*)gauss);
     1677        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1bubbleEnum);
    16761678
    16771679}
    … …  
    21032105
    21042106        /*Recover element type*/
    2105         this->SetElementType(finiteelement_type,analysis_counter);
     2107        this->element_type_list[analysis_counter]=finiteelement_type;
    21062108
    21072109        /*Recover vertices ids needed to initialize inputs*/
    … …  
    24002402
    24012403        GaussPenta* gauss=new GaussPenta();
    2402         for(int iv=0;iv<this->NumberofNodes();iv++){
     2404        for(int iv=0;iv<this->NumberofNodes(this->element_type);iv++){
    24032405                gauss->GaussNode(this->element_type,iv);
    24042406                onbase->GetInputValue(&isonbase,gauss);
    … …  
    24382440/*}}}*/
    24392441void       Penta::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    2440         PentaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss);
     2442        PentaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
    24412443}
    24422444/*}}}*/
    … …  
    24762478/*}}}*/
    24772479int        Penta::VelocityInterpolation(void){/*{{{*/
    2478         return PentaRef::VelocityInterpolation();
     2480        return PentaRef::VelocityInterpolation(this->element_type);
    24792481}
    24802482/*}}}*/
    24812483int        Penta::PressureInterpolation(void){/*{{{*/
    2482         return PentaRef::PressureInterpolation();
     2484        return PentaRef::PressureInterpolation(this->element_type);
    24832485}
    24842486/*}}}*/

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

    r18075 r18078  
    2828/*Object constructors and destructor*/
    2929PentaRef::PentaRef(){/*{{{*/
    30         this->element_type_list=NULL;
    31 }
    32 /*}}}*/
    33 PentaRef::PentaRef(const int nummodels){/*{{{*/
    34 
    35         /*Only allocate pointer*/
    36         element_type_list=xNew<int>(nummodels);
    37 
    3830}
    3931/*}}}*/
    4032PentaRef::~PentaRef(){/*{{{*/
    41         xDelete<int>(element_type_list);
    42 }
    43 /*}}}*/
    44 
    45 /*Management*/
    46 void PentaRef::SetElementType(int type,int type_counter){/*{{{*/
    47 
    48         /*initialize element type*/
    49         this->element_type_list[type_counter]=type;
    5033}
    5134/*}}}*/
    … …  
    167150        /*Invert Jacobian matrix: */
    168151        Matrix3x3Invert(Jinv,&J[0][0]);
    169 }
    170 /*}}}*/
    171 void PentaRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss){/*{{{*/
    172         /*This routine returns the values of the nodal functions  at the gaussian point.*/
    173 
    174         _assert_(basis);
    175         GetNodalFunctions(basis,gauss,this->element_type);
    176 
    177152}
    178153/*}}}*/
    … …  
    321296}
    322297/*}}}*/
    323 void PentaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss){/*{{{*/
    324         GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,this->element_type);
    325 }
    326 /*}}}*/
    327298void PentaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
    328299
    … …  
    356327        /*Clean up*/
    357328        xDelete<IssmDouble>(dbasis_ref);
    358 }
    359 /*}}}*/
    360 void PentaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss){/*{{{*/
    361         GetNodalFunctionsDerivativesReference(dbasis,gauss,this->element_type);
    362329}
    363330/*}}}*/
    … …  
    784751}
    785752/*}}}*/
    786 void PentaRef::GetNodalFunctionsMINIDerivatives(IssmDouble* dbasismini,IssmDouble* xyz_list, Gauss* gauss){/*{{{*/
    787 
    788         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    789          * actual coordinate system): */
    790 
    791         IssmDouble    dbasismini_ref[3][NUMNODESP1b];
    792         IssmDouble    Jinv[3][3];
    793 
    794         /*Get derivative values with respect to parametric coordinate system: */
    795         GetNodalFunctionsMINIDerivativesReference(&dbasismini_ref[0][0], gauss);
    796 
    797         /*Get Jacobian invert: */
    798         GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
    799 
    800         /*Build dbasis:
    801          *
    802          * [dhi/dx]= Jinv'*[dhi/dr]
    803          * [dhi/dy]        [dhi/ds]
    804          * [dhi/dz]        [dhi/dzeta]
    805          */
    806 
    807         for(int i=0;i<NUMNODESP1b;i++){
    808                 *(dbasismini+NUMNODESP1b*0+i)=Jinv[0][0]*dbasismini_ref[0][i]+Jinv[0][1]*dbasismini_ref[1][i]+Jinv[0][2]*dbasismini_ref[2][i];
    809                 *(dbasismini+NUMNODESP1b*1+i)=Jinv[1][0]*dbasismini_ref[0][i]+Jinv[1][1]*dbasismini_ref[1][i]+Jinv[1][2]*dbasismini_ref[2][i];
    810                 *(dbasismini+NUMNODESP1b*2+i)=Jinv[2][0]*dbasismini_ref[0][i]+Jinv[2][1]*dbasismini_ref[1][i]+Jinv[2][2]*dbasismini_ref[2][i];
    811         }
    812 
    813 }
    814 /*}}}*/
    815 void PentaRef::GetNodalFunctionsMINIDerivativesReference(IssmDouble* dbasis,Gauss* gauss_in){/*{{{*/
    816         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    817          * natural coordinate system) at the gaussian point. */
    818 
    819         /*Cast gauss to GaussPenta*/
    820         _assert_(gauss_in->Enum()==GaussPentaEnum);
    821         GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
    822 
    823 
    824         IssmDouble zeta=gauss->coord4;
    825 
    826         /*Nodal function 1*/
    827         dbasis[NUMNODESP1b*0+0]=-0.5*(1.0-zeta)/2.0;
    828         dbasis[NUMNODESP1b*1+0]=-SQRT3/6.0*(1.0-zeta)/2.0;
    829         dbasis[NUMNODESP1b*2+0]=-0.5*gauss->coord1;
    830         /*Nodal function 2*/
    831         dbasis[NUMNODESP1b*0+1]=0.5*(1.0-zeta)/2.0;
    832         dbasis[NUMNODESP1b*1+1]=-SQRT3/6.0*(1.0-zeta)/2.0;
    833         dbasis[NUMNODESP1b*2+1]=-0.5*gauss->coord2;
    834         /*Nodal function 3*/
    835         dbasis[NUMNODESP1b*0+2]=0.;
    836         dbasis[NUMNODESP1b*1+2]=SQRT3/3.0*(1.0-zeta)/2.0;
    837         dbasis[NUMNODESP1b*2+2]=-0.5*gauss->coord3;
    838         /*Nodal function 4*/
    839         dbasis[NUMNODESP1b*0+3]=-0.5*(1.0+zeta)/2.0;
    840         dbasis[NUMNODESP1b*1+3]=-SQRT3/6.0*(1.0+zeta)/2.0;
    841         dbasis[NUMNODESP1b*2+3]=0.5*gauss->coord1;
    842         /*Nodal function 5*/
    843         dbasis[NUMNODESP1b*0+4]=0.5*(1.0+zeta)/2.0;
    844         dbasis[NUMNODESP1b*1+4]=-SQRT3/6.0*(1.0+zeta)/2.0;
    845         dbasis[NUMNODESP1b*2+4]=0.5*gauss->coord2;
    846         /*Nodal function 6*/
    847         dbasis[NUMNODESP1b*0+5]=0.;
    848         dbasis[NUMNODESP1b*1+5]=SQRT3/3.0*(1.0+zeta)/2.0;
    849         dbasis[NUMNODESP1b*2+5]=0.5*gauss->coord3;
    850         /*Nodal function 7*/
    851         dbasis[NUMNODESP1b*0+6]=27.*(1.+zeta)*(1.-zeta)*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3);
    852         dbasis[NUMNODESP1b*1+6]=27.*(1.+zeta)*(1.-zeta)*SQRT3*(-1./6.*gauss->coord2*gauss->coord3 - 1./6.*gauss->coord1*gauss->coord3 +1./3.*gauss->coord1*gauss->coord2);
    853         dbasis[NUMNODESP1b*2+6]=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);
    854 }
    855 /*}}}*/
    856 void PentaRef::GetNodalFunctionsP1(IssmDouble* basis, Gauss* gauss_in){/*{{{*/
    857         /*This routine returns the values of the nodal functions  at the gaussian point.*/
    858 
    859         /*Cast gauss to GaussPenta*/
    860         _assert_(gauss_in->Enum()==GaussPentaEnum);
    861         GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
    862 
    863         basis[0]=gauss->coord1*(1-gauss->coord4)/2.0;
    864         basis[1]=gauss->coord2*(1-gauss->coord4)/2.0;
    865         basis[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    866         basis[3]=gauss->coord1*(1+gauss->coord4)/2.0;
    867         basis[4]=gauss->coord2*(1+gauss->coord4)/2.0;
    868         basis[5]=gauss->coord3*(1+gauss->coord4)/2.0;
    869 
    870 }
    871 /*}}}*/
    872 void PentaRef::GetNodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss){/*{{{*/
    873 
    874         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    875          * actual coordinate system): */
    876         IssmDouble    dbasis_ref[NDOF3][NUMNODESP1];
    877         IssmDouble    Jinv[NDOF3][NDOF3];
    878 
    879         /*Get derivative values with respect to parametric coordinate system: */
    880         GetNodalFunctionsP1DerivativesReference(&dbasis_ref[0][0], gauss);
    881 
    882         /*Get Jacobian invert: */
    883         GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
    884 
    885         /*Build basis function derivatives:
    886          *
    887          * [dhi/dx]= Jinv*[dhi/dr]
    888          * [dhi/dy]       [dhi/ds]
    889          * [dhi/dz]       [dhi/dn]
    890          */
    891 
    892         for (int i=0;i<NUMNODESP1;i++){
    893                 *(dbasis+NUMNODESP1*0+i)=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i]+Jinv[0][2]*dbasis_ref[2][i];
    894                 *(dbasis+NUMNODESP1*1+i)=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i]+Jinv[1][2]*dbasis_ref[2][i];
    895                 *(dbasis+NUMNODESP1*2+i)=Jinv[2][0]*dbasis_ref[0][i]+Jinv[2][1]*dbasis_ref[1][i]+Jinv[2][2]*dbasis_ref[2][i];
    896         }
    897 
    898 }
    899 /*}}}*/
    900 void PentaRef::GetNodalFunctionsP1DerivativesReference(IssmDouble* dbasis,Gauss* gauss_in){/*{{{*/
    901 
    902         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    903          * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */
    904 
    905         /*Cast gauss to GaussPenta*/
    906         _assert_(gauss_in->Enum()==GaussPentaEnum);
    907         GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
    908 
    909         IssmDouble zeta=gauss->coord4;
    910 
    911         /*Nodal function 1*/
    912         dbasis[NUMNODESP1*0+0]=-0.5*(1.0-zeta)/2.0;
    913         dbasis[NUMNODESP1*1+0]=-SQRT3/6.0*(1.0-zeta)/2.0;
    914         dbasis[NUMNODESP1*2+0]=-0.5*gauss->coord1;
    915         /*Nodal function 2*/
    916         dbasis[NUMNODESP1*0+1]=0.5*(1.0-zeta)/2.0;
    917         dbasis[NUMNODESP1*1+1]=-SQRT3/6.0*(1.0-zeta)/2.0;
    918         dbasis[NUMNODESP1*2+1]=-0.5*gauss->coord2;
    919         /*Nodal function 3*/
    920         dbasis[NUMNODESP1*0+2]=0.;
    921         dbasis[NUMNODESP1*1+2]=SQRT3/3.0*(1.0-zeta)/2.0;
    922         dbasis[NUMNODESP1*2+2]=-0.5*gauss->coord3;
    923         /*Nodal function 4*/
    924         dbasis[NUMNODESP1*0+3]=-0.5*(1.0+zeta)/2.0;
    925         dbasis[NUMNODESP1*1+3]=-SQRT3/6.0*(1.0+zeta)/2.0;
    926         dbasis[NUMNODESP1*2+3]=0.5*gauss->coord1;
    927         /*Nodal function 5*/
    928         dbasis[NUMNODESP1*0+4]=0.5*(1.0+zeta)/2.0;
    929         dbasis[NUMNODESP1*1+4]=-SQRT3/6.0*(1.0+zeta)/2.0;
    930         dbasis[NUMNODESP1*2+4]=0.5*gauss->coord2;
    931         /*Nodal function 6*/
    932         dbasis[NUMNODESP1*0+5]=0.;
    933         dbasis[NUMNODESP1*1+5]=SQRT3/3.0*(1.0+zeta)/2.0;
    934         dbasis[NUMNODESP1*2+5]=0.5*gauss->coord3;
    935 }
    936 /*}}}*/
    937753void PentaRef::GetQuadJacobianDeterminant(IssmDouble* Jdet,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    938754        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    … …  
    960776}
    961777/*}}}*/
    962 void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss){/*{{{*/
    963 
    964         GetInputValue(pvalue,plist,gauss,this->element_type);
    965 
    966 }
    967 /*}}}*/
    968778void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss,int finiteelement){/*{{{*/
    969779
    … …  
    987797}
    988798/*}}}*/
    989 void PentaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss){/*{{{*/
     799void PentaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
    990800        /*From node values of parameter p (p_list[0], p_list[1], p_list[2],
    991801         * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
    … …  
    1004814
    1005815        /*Fetch number of nodes for this finite element*/
    1006         int numnodes = this->NumberofNodes();
     816        int numnodes = this->NumberofNodes(finiteelement);
    1007817
    1008818        /*Get nodal functions derivatives*/
    1009819        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
    1010         GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     820        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,finiteelement);
    1011821
    1012822        /*Calculate parameter for this Gauss point*/
    … …  
    1021831        p[2]=dpz;
    1022832
    1023 }
    1024 /*}}}*/
    1025 int  PentaRef::NumberofNodes(void){/*{{{*/
    1026 
    1027         return this->NumberofNodes(this->element_type);
    1028833}
    1029834/*}}}*/
    … …  
    1046851                case P2xP4Enum:             return NUMNODESP2xP4;
    1047852                case P1xP3Enum:             return NUMNODESP1xP3;
    1048                 default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     853                default:       _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    1049854        }
    1050855
    … …  
    1052857}
    1053858/*}}}*/
    1054 int  PentaRef::VelocityInterpolation(void){/*{{{*/
    1055 
    1056         switch(this->element_type){
     859int  PentaRef::VelocityInterpolation(int fe_stokes){/*{{{*/
     860
     861        switch(fe_stokes){
    1057862                case P1P1Enum:          return P1Enum;
    1058863                case P1P1GLSEnum:       return P1Enum;
    … …  
    1061866                case TaylorHoodEnum:    return P2Enum;
    1062867                case OneLayerP4zEnum:   return P2xP4Enum;
    1063                 default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     868                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    1064869        }
    1065870
    … …  
    1067872}
    1068873/*}}}*/
    1069 int  PentaRef::PressureInterpolation(void){/*{{{*/
    1070 
    1071         switch(this->element_type){
     874int  PentaRef::PressureInterpolation(int fe_stokes){/*{{{*/
     875
     876        switch(fe_stokes){
    1072877                case P1P1Enum:          return P1Enum;
    1073878                case P1P1GLSEnum:       return P1Enum;
    … …  
    1076881                case TaylorHoodEnum:    return P1Enum;
    1077882                case OneLayerP4zEnum:   return P1Enum;
    1078                 default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     883                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    1079884        }
    1080885
    … …  
    1082887}
    1083888/*}}}*/
    1084 int  PentaRef::TensorInterpolation(void){/*{{{*/
    1085 
    1086         switch(this->element_type){
     889int  PentaRef::TensorInterpolation(int fe_stokes){/*{{{*/
     890
     891        switch(fe_stokes){
    1087892                case XTaylorHoodEnum:    return P1DGEnum;
    1088                 default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     893                default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    1089894        }
    1090895
    … …  
    1158963                        break;
    1159964                default:
    1160                         _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     965                        _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    1161966        }
    1162967
    … …  
    12151020                        break;
    12161021                default:
    1217                         _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     1022                        _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    12181023        }
    12191024

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

    r18075 r18078  
    1111
    1212        public:
    13                 int* element_type_list; //P1CG, P1DG, MINI, P2...
    14                 int  element_type;
    15 
    1613                PentaRef();
    17                 PentaRef(const int nummodels);
    1814                ~PentaRef();
    1915
    20                 /*Management*/
    21                 void SetElementType(int type,int type_counter);
    22 
    2316                /*Numerics*/
    24                 void GetNodalFunctions(IssmDouble* basis, Gauss* gauss);
    2517                void GetNodalFunctions(IssmDouble* basis, Gauss* gauss,int finiteelement);
    26                 void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    2718                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss,int finiteelement);
    28                 void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss);
    2919                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss,int finiteelement);
    30                 void GetNodalFunctionsP1(IssmDouble* basis, Gauss* gauss);
    31                 void GetNodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss);
    32                 void GetNodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss);
    33                 void GetNodalFunctionsP1DerivativesReference(IssmDouble* dl1dl6,Gauss* gauss);
    34                 void GetNodalFunctionsMINIDerivativesReference(IssmDouble* dl1dl7,Gauss* gauss);
    3520                void GetQuadJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
    3621                void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss);
    … …  
    4025                void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,Gauss* gauss);
    4126                void GetLprimeFSSSA(IssmDouble* LprimeFSSSA, IssmDouble* xyz_list, Gauss* gauss);
    42                 void GetInputValue(IssmDouble* pvalue,IssmDouble* plist, Gauss* gauss);
    4327                void GetInputValue(IssmDouble* pvalue,IssmDouble* plist, Gauss* gauss,int finiteelement);
    44                 void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss);
     28                void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
    4529
    4630                void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement);
    4731                void SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement);
    48                 int  NumberofNodes(void);
    4932                int  NumberofNodes(int finiteelement);
    50                 int  VelocityInterpolation(void);
    51                 int  PressureInterpolation(void);
    52                 int  TensorInterpolation(void);
     33                int  VelocityInterpolation(int fe_stokes);
     34                int  PressureInterpolation(int fe_stokes);
     35                int  TensorInterpolation(int fe_stokes);
    5336};
    5437#endif

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

    r18062 r18078  
    2020/*Constructors/destructor/copy*/
    2121Seg::Seg(int seg_id, int seg_sid, int index, IoModel* iomodel,int nummodels)/*{{{*/
    22                 :SegRef(nummodels),ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
     22                :ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
    2323
    2424                        /*id: */
    … …  
    3737                        this->material = NULL;
    3838                        this->matpar   = NULL;
     39
     40                        /*Only allocate pointer*/
     41                        this->element_type_list=xNew<int>(nummodels);
    3942                }
    4043/*}}}*/
    … …  
    100103}/*}}}*/
    101104int        Seg::GetNumberOfNodes(void){/*{{{*/
    102         return this->NumberofNodes();
     105        return this->NumberofNodes(this->element_type);
    103106}
    104107/*}}}*/
    … …  
    177180
    178181        _assert_(gauss->Enum()==GaussSegEnum);
    179         this->GetNodalFunctions(basis,(GaussSeg*)gauss);
     182        this->GetNodalFunctions(basis,(GaussSeg*)gauss,this->element_type);
    180183
    181184}
    … …  
    184187
    185188        _assert_(gauss->Enum()==GaussSegEnum);
    186         this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussSeg*)gauss);
     189        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussSeg*)gauss,this->element_type);
    187190
    188191}

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

    r17962 r18078  
    2121/*Object constructors and destructor*/
    2222SegRef::SegRef(){/*{{{*/
    23         this->element_type_list=NULL;
    24 }
    25 /*}}}*/
    26 SegRef::SegRef(const int nummodels){/*{{{*/
    27 
    28         /*Only allocate pointer*/
    29         element_type_list=xNew<int>(nummodels);
    30 
    3123}
    3224/*}}}*/
    3325SegRef::~SegRef(){/*{{{*/
    34         xDelete<int>(element_type_list);
    35 }
    36 /*}}}*/
    37 
    38 /*Management*/
    39 void SegRef::SetElementType(int type,int type_counter){/*{{{*/
    40 
    41         /*initialize element type*/
    42         this->element_type_list[type_counter]=type;
    4326}
    4427/*}}}*/
    4528
    4629/*Reference Element numerics*/
    47 void SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss){/*{{{*/
    48         /*This routine returns the values of the nodal functions  at the gaussian point.*/
    49 
    50         _assert_(basis);
    51 
    52         GetNodalFunctions(basis,gauss,this->element_type);
    53 }
    54 /*}}}*/
    5530void SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss,int finiteelement){/*{{{*/
    5631        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    … …  
    7449                        _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    7550        }
    76 }
    77 /*}}}*/
    78 void SegRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss){/*{{{*/
    79 
    80         GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,this->element_type);
    81 
    8251}
    8352/*}}}*/
    … …  
    10776        /*Clean up*/
    10877        xDelete<IssmDouble>(dbasis_ref);
    109 
    110 }
    111 /*}}}*/
    112 void SegRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss){/*{{{*/
    113         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    114          * natural coordinate system) at the gaussian point. */
    115 
    116         GetNodalFunctionsDerivativesReference(dbasis,gauss,this->element_type);
    11778
    11879}
    … …  
    147108}
    148109/*}}}*/
    149 void SegRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss){/*{{{*/
     110void SegRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement){/*{{{*/
    150111
    151112        /*From node values of parameter p (plist[0],plist[1]), return parameter derivative value at gaussian
    … …  
    160121
    161122        /*Fetch number of nodes for this finite element*/
    162         int numnodes = this->NumberofNodes();
     123        int numnodes = this->NumberofNodes(finiteelement);
    163124
    164125        /*Get nodal functions derivatives*/
    165126        IssmDouble* dbasis=xNew<IssmDouble>(1*numnodes);
    166         GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     127        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,finiteelement);
    167128
    168129        /*Calculate parameter for this Gauss point*/
    … …  
    175136}
    176137/*}}}*/
    177 void SegRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss){/*{{{*/
    178 
    179         GetInputValue(p,plist,gauss,this->element_type);
    180 }
    181 /*}}}*/
    182138void SegRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss,int finiteelement){/*{{{*/
    183139
    … …  
    230186        /*Invert Jacobian matrix: */
    231187        *Jinv = 1./J;
    232 }
    233 /*}}}*/
    234 int  SegRef::NumberofNodes(void){/*{{{*/
    235 
    236         return this->NumberofNodes(this->element_type);
    237188}
    238189/*}}}*/
    … …  
    243194                case P1Enum:                return NUMNODESP1;
    244195                case P1DGEnum:              return NUMNODESP1;
    245                 default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     196                default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    246197        }
    247198

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

    r17926 r18078  
    1313
    1414        public:
    15                 int* element_type_list;
    16                 int  element_type;
    17 
    1815                SegRef();
    19                 SegRef(const int nummodels);
    2016                ~SegRef();
    2117
    22                 /*Management*/
    23                 void SetElementType(int type,int type_counter);
    2418                void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussSeg* gauss);
    2519                void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussSeg* gauss);
    2620                void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussSeg* gauss);
    27                 void GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss);
    2821                void GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss,int finiteelement);
    29                 void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss);
    3022                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement);
    31                 void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss);
    3223                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss,int finiteelement);
    33                 void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss);
    34                 void GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss);
     24                void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement);
    3525                void GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss,int finiteelement);
    36 
    37                 int  NumberofNodes(void);
    3826                int  NumberofNodes(int finiteelement);
    3927};

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

    r17962 r18078  
    2121/*Constructors/destructor/copy*/
    2222Tetra::Tetra(int seg_id, int seg_sid, int index, IoModel* iomodel,int nummodels)/*{{{*/
    23                 :TetraRef(nummodels),ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
     23                :ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
    2424
    2525                        /*id: */
    … …  
    3838                        this->material = NULL;
    3939                        this->matpar   = NULL;
     40
     41                        /*Only allocate pointer*/
     42                        this->element_type_list=xNew<int>(nummodels);
    4043                }
    4144/*}}}*/
    … …  
    202205
    203206        _assert_(nodes);
    204         int numnodes = this->NumberofNodes();
     207        int numnodes = this->NumberofNodes(this->element_type);
    205208
    206209        for(int i=0;i<numnodes;i++){
    … …  
    212215/*}}}*/
    213216int      Tetra::GetNumberOfNodes(void){/*{{{*/
    214         return this->NumberofNodes();
     217        return this->NumberofNodes(this->element_type);
    215218}
    216219/*}}}*/
    … …  
    400403
    401404        /*Fetch number of nodes for this finite element*/
    402         int numnodes = this->NumberofNodes();
     405        int numnodes = this->NumberofNodes(this->element_type);
    403406
    404407        /*Fetch dof list and allocate solution vector*/
    … …  
    502505
    503506        _assert_(gauss->Enum()==GaussTetraEnum);
    504         this->GetNodalFunctions(basis,(GaussTetra*)gauss);
     507        this->GetNodalFunctions(basis,(GaussTetra*)gauss,this->element_type);
    505508
    506509}
    … …  
    530533
    531534        _assert_(gauss->Enum()==GaussTetraEnum);
    532         this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTetra*)gauss);
     535        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTetra*)gauss,this->element_type);
    533536
    534537}
    … …  
    788791
    789792        /*Recover element type*/
    790         this->SetElementType(finiteelement_type,analysis_counter);
     793        this->element_type_list[analysis_counter]=finiteelement_type;
    791794
    792795        /*Recover vertices ids needed to initialize inputs*/
    … …  
    872875/*}}}*/
    873876int      Tetra::VelocityInterpolation(void){/*{{{*/
    874         return TetraRef::VelocityInterpolation();
     877        return TetraRef::VelocityInterpolation(this->element_type);
    875878}
    876879/*}}}*/
    … …  
    892895/*}}}*/
    893896int      Tetra::PressureInterpolation(void){/*{{{*/
    894         return TetraRef::PressureInterpolation();
     897        return TetraRef::PressureInterpolation(this->element_type);
    895898}
    896899/*}}}*/
    897900int      Tetra::TensorInterpolation(void){/*{{{*/
    898         return TetraRef::TensorInterpolation();
     901        return TetraRef::TensorInterpolation(this->element_type);
    899902}
    900903/*}}}*/

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

    r17962 r18078  
    44
    55/*Headers:*/
    6 /*{{{*//*{{{*/
     6/*{{{*/
    77#ifdef HAVE_CONFIG_H
    88#include <config.h>
    … …  
    2323/*Object constructors and destructor*/
    2424TetraRef::TetraRef(){/*{{{*/
    25         this->element_type_list=NULL;
    26 }
    27 /*}}}*/
    28 TetraRef::TetraRef(const int nummodels){/*{{{*/
    29 
    30         /*Only allocate pointer*/
    31         element_type_list=xNew<int>(nummodels);
    32 
    3325}
    3426/*}}}*/
    3527TetraRef::~TetraRef(){/*{{{*/
    36         xDelete<int>(element_type_list);
    37 }
    38 /*}}}*/
    39 
    40 /*Management*/
    41 void TetraRef::SetElementType(int type,int type_counter){/*{{{*/
    42 
    43         /*initialize element type*/
    44         this->element_type_list[type_counter]=type;
    4528}
    4629/*}}}*/
    4730
    4831/*Reference Element numerics*/
    49 void TetraRef::GetNodalFunctions(IssmDouble* basis,GaussTetra* gauss){/*{{{*/
    50         /*This routine returns the values of the nodal functions  at the gaussian point.*/
    51         _assert_(basis);
    52         GetNodalFunctions(basis,gauss,this->element_type);
    53 }
    54 /*}}}*/
    5532void TetraRef::GetNodalFunctions(IssmDouble* basis,GaussTetra* gauss,int finiteelement){/*{{{*/
    5633        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    … …  
    9673}
    9774/*}}}*/
    98 void TetraRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss){/*{{{*/
    99         GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,this->element_type);
    100 }
    101 /*}}}*/
    10275void TetraRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement){/*{{{*/
    10376
    … …  
    133106}
    134107/*}}}*/
    135 void TetraRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss){/*{{{*/
    136         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    137          * natural coordinate system) at the gaussian point. */
    138 
    139         GetNodalFunctionsDerivativesReference(dbasis,gauss,this->element_type);
    140 
    141 }
    142 /*}}}*/
    143108void TetraRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss,int finiteelement){/*{{{*/
    144109        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    … …  
    239204}
    240205/*}}}*/
    241 void TetraRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss){/*{{{*/
     206void TetraRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement){/*{{{*/
    242207        /*From node values of parameter p (p_list[0], p_list[1], p_list[2],
    243208         * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
    … …  
    256221
    257222        /*Fetch number of nodes for this finite element*/
    258         int numnodes = this->NumberofNodes();
     223        int numnodes = this->NumberofNodes(finiteelement);
    259224
    260225        /*Get nodal functions derivatives*/
    261226        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
    262         GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     227        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,finiteelement);
    263228
    264229        /*Calculate parameter for this Gauss point*/
    … …  
    272237        p[1]=dpy;
    273238        p[2]=dpz;
    274 }
    275 /*}}}*/
    276 void TetraRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTetra* gauss){/*{{{*/
    277 
    278         GetInputValue(p,plist,gauss,this->element_type);
    279239}
    280240/*}}}*/
    … …  
    374334        /*Invert Jacobian matrix: */
    375335        Matrix3x3Invert(Jinv,&J[0][0]);
    376 }
    377 /*}}}*/
    378 int  TetraRef::NumberofNodes(void){/*{{{*/
    379 
    380         return this->NumberofNodes(this->element_type);
    381336}
    382337/*}}}*/
    … …  
    395350                case MINIEnum:              return NUMNODESP1b+NUMNODESP1;
    396351                case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
    397                 default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     352                default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    398353        }
    399354
    … …  
    401356}
    402357/*}}}*/
    403 int  TetraRef::VelocityInterpolation(void){/*{{{*/
    404 
    405         switch(this->element_type){
     358int  TetraRef::VelocityInterpolation(int fe_stokes){/*{{{*/
     359
     360        switch(fe_stokes){
    406361                case P1P1Enum:          return P1Enum;
    407362                case P1P1GLSEnum:       return P1Enum;
    … …  
    409364                case MINIEnum:          return P1bubbleEnum;
    410365                case TaylorHoodEnum:    return P2Enum;
    411                 default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     366                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    412367        }
    413368
    … …  
    415370}
    416371/*}}}*/
    417 int TetraRef::PressureInterpolation(void){/*{{{*/
    418 
    419         switch(this->element_type){
     372int TetraRef::PressureInterpolation(int fe_stokes){/*{{{*/
     373
     374        switch(fe_stokes){
    420375                case P1P1Enum:          return P1Enum;
    421376                case P1P1GLSEnum:       return P1Enum;
    … …  
    423378                case MINIEnum:          return P1Enum;
    424379                case TaylorHoodEnum:    return P1Enum;
    425                 default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     380                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    426381        }
    427382
    428383        return -1;
    429384}/*}}}*/
    430 int  TetraRef::TensorInterpolation(void){/*{{{*/
     385int  TetraRef::TensorInterpolation(int fe_stokes){/*{{{*/
    431386        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    432387
    433         switch(this->element_type){
     388        switch(fe_stokes){
    434389                case XTaylorHoodEnum: return P1DGEnum;
    435                 default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
    436         }
    437 }
    438 /*}}}*/
     390                default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
     391        }
     392}
     393/*}}}*/

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

    r17926 r18078  
    1313
    1414        public:
    15                 int* element_type_list;
    16                 int  element_type;
    17 
    1815                TetraRef();
    19                 TetraRef(const int nummodels);
    2016                ~TetraRef();
    2117
    22                 /*Management*/
    23                 void SetElementType(int type,int type_counter);
    2418                void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTetra* gauss);
    2519                void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss);
    2620                void GetJacobianDeterminantFace(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss);
    2721                void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTetra* gauss);
    28                 void GetNodalFunctions(IssmDouble* basis,GaussTetra* gauss);
    2922                void GetNodalFunctions(IssmDouble* basis,GaussTetra* gauss,int finiteelement);
    30                 void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss);
    3123                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement);
    32                 void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss);
    3324                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss,int finiteelement);
    34                 void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss);
    35                 void GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTetra* gauss);
     25                void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement);
    3626                void GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTetra* gauss,int finiteelement);
    3727
    38                 int  NumberofNodes(void);
    3928                int  NumberofNodes(int finiteelement);
    40                 int  VelocityInterpolation(void);
    41                 int  PressureInterpolation(void);
    42                 int  TensorInterpolation(void);
     29                int  VelocityInterpolation(int fe_stokes);
     30                int  PressureInterpolation(int fe_stokes);
     31                int  TensorInterpolation(int fe_stokes);
    4332};
    4433#endif

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

    r18068 r18078  
    2525/*Constructors/destructor/copy*/
    2626Tria::Tria(int tria_id, int tria_sid, int index, IoModel* iomodel,int nummodels)/*{{{*/
    27         :TriaRef(nummodels),ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
     27        :ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
    2828
    2929                /*id: */
    … …  
    4242                this->material = NULL;
    4343                this->matpar   = NULL;
     44                this->element_type_list=xNew<int>(nummodels);
    4445}
    4546/*}}}*/
    … …  
    900901/*}}}*/
    901902int        Tria::GetNumberOfNodes(void){/*{{{*/
    902         return this->NumberofNodes();
     903        return this->NumberofNodes(this->element_type);
    903904}
    904905/*}}}*/
    … …  
    921922Node*      Tria::GetNode(int node_number){/*{{{*/
    922923        _assert_(node_number>=0);
    923         _assert_(node_number<this->NumberofNodes());
     924        _assert_(node_number<this->NumberofNodes(this->element_type));
    924925        return this->nodes[node_number];
    925926
    … …  
    10591060
    10601061        /*Fetch number of nodes for this finite element*/
    1061         int numnodes = this->NumberofNodes();
     1062        int numnodes = this->NumberofNodes(this->element_type);
    10621063
    10631064        /*Fetch dof list and allocate solution vector*/
    … …  
    11101111
    11111112                /*Get number of nodes and dof list: */
    1112                 numnodes = this->NumberofNodes();
     1113                numnodes = this->NumberofNodes(this->element_type);
    11131114                values   = xNew<IssmDouble>(numnodes);
    11141115                GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    … …  
    11241125
    11251126                /*Get number of nodes and dof list: */
    1126                 numnodes = this->NumberofNodes();
     1127                numnodes = this->NumberofNodes(this->element_type);
    11271128                values   = xNew<IssmDouble>(numnodes);
    11281129
    … …  
    14531454
    14541455        _assert_(gauss->Enum()==GaussTriaEnum);
    1455         this->GetNodalFunctions(basis,(GaussTria*)gauss);
     1456        this->GetNodalFunctions(basis,(GaussTria*)gauss,this->element_type);
    14561457
    14571458}
    … …  
    14671468
    14681469        _assert_(gauss->Enum()==GaussTriaEnum);
    1469         this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss);
     1470        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->element_type);
    14701471
    14711472}
    … …  
    15731574/*}}}*/
    15741575int        Tria::VelocityInterpolation(void){/*{{{*/
    1575         return TriaRef::VelocityInterpolation();
     1576        return TriaRef::VelocityInterpolation(this->element_type);
    15761577}
    15771578/*}}}*/
    15781579int        Tria::PressureInterpolation(void){/*{{{*/
    1579         return TriaRef::PressureInterpolation();
     1580        return TriaRef::PressureInterpolation(this->element_type);
    15801581}
    15811582/*}}}*/
    15821583int        Tria::TensorInterpolation(void){/*{{{*/
    1583         return TriaRef::TensorInterpolation();
     1584        return TriaRef::TensorInterpolation(this->element_type);
    15841585}
    15851586/*}}}*/
    … …  
    18901891
    18911892        /*Recover element type*/
    1892         this->SetElementType(finiteelement_type,analysis_counter);
     1893        this->element_type_list[analysis_counter]=finiteelement_type;
    18931894
    18941895        /*Recover nodes ids needed to initialize the node hook.*/
    … …  
    19881989
    19891990        GaussTria* gauss=new GaussTria();
    1990         for(int iv=0;iv<this->NumberofNodes();iv++){
     1991        for(int iv=0;iv<this->NumberofNodes(this->element_type);iv++){
    19911992                gauss->GaussNode(this->element_type,iv);
    19921993                onbase->GetInputValue(&isonbase,gauss);
    … …  
    20272028/*}}}*/
    20282029void       Tria::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    2029         TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss);
     2030        TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
    20302031}
    20312032/*}}}*/
    … …  
    28282829
    28292830        /*Fetch number of nodes for this finite element*/
    2830         int numnodes = this->NumberofNodes();
     2831        int numnodes = this->NumberofNodes(this->element_type);
    28312832
    28322833        /*Fetch dof list and allocate solution vector*/

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

    r18060 r18078  
    2323/*Object constructors and destructor*/
    2424TriaRef::TriaRef(){/*{{{*/
    25         this->element_type_list=NULL;
    26 }
    27 /*}}}*/
    28 TriaRef::TriaRef(const int nummodels){/*{{{*/
    29 
    30         /*Only allocate pointer*/
    31         element_type_list=xNew<int>(nummodels);
    32 
    3325}
    3426/*}}}*/
    3527TriaRef::~TriaRef(){/*{{{*/
    36         xDelete<int>(element_type_list);
    37 }
    38 /*}}}*/
    39 
    40 /*Management*/
    41 void TriaRef::SetElementType(int type,int type_counter){/*{{{*/
    42 
    43         /*initialize element type*/
    44         this->element_type_list[type_counter]=type;
    4528}
    4629/*}}}*/
    4730
    4831/*Reference Element numerics*/
    49 void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2){/*{{{*/
     32void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
    5033        /*Compute B  matrix. B=[phi1 phi2 -phi3 -phi4]
    5134         *
    … …  
    5639
    5740        /*Fetch number of nodes for this finite element*/
    58         int numnodes = this->NumberofNodes();
     41        int numnodes = this->NumberofNodes(finiteelement);
    5942
    6043        /*Get nodal functions*/
    6144        IssmDouble* basis=xNew<IssmDouble>(numnodes);
    62         GetNodalFunctions(basis,gauss);
     45        GetNodalFunctions(basis,gauss,finiteelement);
    6346
    6447        /*Build B for this segment*/
    … …  
    7255}
    7356/*}}}*/
    74 void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2){/*{{{*/
     57void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
    7558        /*Compute Bprime  matrix. Bprime=[phi1 phi2 phi3 phi4]
    7659         *
    … …  
    8164
    8265        /*Fetch number of nodes for this finite element*/
    83         int numnodes = this->NumberofNodes();
     66        int numnodes = this->NumberofNodes(finiteelement);
    8467
    8568        /*Get nodal functions*/
    8669        IssmDouble* basis=xNew<IssmDouble>(numnodes);
    87         GetNodalFunctions(basis,gauss);
     70        GetNodalFunctions(basis,gauss,finiteelement);
    8871
    8972        /*Build B'*/
    … …  
    152135        /*Invert Jacobian matrix: */
    153136        Matrix2x2Invert(Jinv,&J[0][0]);
    154 
    155 }
    156 /*}}}*/
    157 void TriaRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss){/*{{{*/
    158         /*This routine returns the values of the nodal functions  at the gaussian point.*/
    159 
    160         _assert_(basis);
    161         GetNodalFunctions(basis,gauss,this->element_type);
    162137
    163138}
    … …  
    204179}
    205180/*}}}*/
    206 void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss,int index1,int index2){/*{{{*/
     181void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss,int index1,int index2,int finiteelement){/*{{{*/
    207182        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    208183
    … …  
    211186
    212187        /*Fetch number of nodes for this finite element*/
    213         int numnodes = this->NumberofNodes();
     188        int numnodes = this->NumberofNodes(finiteelement);
    214189
    215190        /*Get nodal functions*/
    216191        IssmDouble* triabasis=xNew<IssmDouble>(numnodes);
    217         GetNodalFunctions(triabasis,gauss);
    218 
    219         switch(this->element_type){
     192        GetNodalFunctions(triabasis,gauss,finiteelement);
     193
     194        switch(finiteelement){
    220195                case P1Enum: case P1DGEnum:
    221196                        basis[0]=triabasis[index1];
    … …  
    236211                        return;
    237212                default:
    238                         _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     213                        _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    239214        }
    240215
    241216        /*Clean up*/
    242217        xDelete<IssmDouble>(triabasis);
    243 }
    244 /*}}}*/
    245 void TriaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss){/*{{{*/
    246 
    247         GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,this->element_type);
    248 
    249218}
    250219/*}}}*/
    … …  
    276245        /*Clean up*/
    277246        xDelete<IssmDouble>(dbasis_ref);
    278 
    279 }
    280 /*}}}*/
    281 void TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss){/*{{{*/
    282         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    283          * natural coordinate system) at the gaussian point. */
    284 
    285         GetNodalFunctionsDerivativesReference(dbasis,gauss,this->element_type);
    286247
    287248}
    … …  
    354315}
    355316/*}}}*/
    356 void TriaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss){/*{{{*/
     317void TriaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
    357318
    358319        /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian
    … …  
    369330
    370331        /*Fetch number of nodes for this finite element*/
    371         int numnodes = this->NumberofNodes();
     332        int numnodes = this->NumberofNodes(finiteelement);
    372333
    373334        /*Get nodal functions derivatives*/
    374335        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
    375         GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     336        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,finiteelement);
    376337
    377338        /*Calculate parameter for this Gauss point*/
    … …  
    386347}
    387348/*}}}*/
    388 void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss){/*{{{*/
    389 
    390         GetInputValue(p,plist,gauss,this->element_type);
    391 }
    392 /*}}}*/
    393349void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss,int finiteelement){/*{{{*/
    394350
    … …  
    409365        xDelete<IssmDouble>(basis);
    410366        *p = value;
    411 }
    412 /*}}}*/
    413 int  TriaRef::NumberofNodes(void){/*{{{*/
    414 
    415         return this->NumberofNodes(this->element_type);
    416367}
    417368/*}}}*/
    … …  
    431382                case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
    432383                case XTaylorHoodEnum:       return NUMNODESP2+NUMNODESP1;
    433                 default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     384                default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    434385        }
    435386
    … …  
    437388}
    438389/*}}}*/
    439 int  TriaRef::VelocityInterpolation(void){/*{{{*/
    440 
    441         switch(this->element_type){
     390int  TriaRef::VelocityInterpolation(int fe_stokes){/*{{{*/
     391
     392        switch(fe_stokes){
    442393                case P1P1Enum:          return P1Enum;
    443394                case P1P1GLSEnum:       return P1Enum;
    … …  
    446397                case TaylorHoodEnum:    return P2Enum;
    447398                case XTaylorHoodEnum:   return P2Enum;
    448                 default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     399                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    449400        }
    450401
    … …  
    452403}
    453404/*}}}*/
    454 int  TriaRef::PressureInterpolation(void){/*{{{*/
    455 
    456         switch(this->element_type){
     405int  TriaRef::PressureInterpolation(int fe_stokes){/*{{{*/
     406
     407        switch(fe_stokes){
    457408                case P1P1Enum:          return P1Enum;
    458409                case P1P1GLSEnum:       return P1Enum;
    … …  
    461412                case TaylorHoodEnum:    return P1Enum;
    462413                case XTaylorHoodEnum:   return P1Enum;
    463                 default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     414                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    464415        }
    465416
    … …  
    467418}
    468419/*}}}*/
    469 int  TriaRef::TensorInterpolation(void){/*{{{*/
     420int  TriaRef::TensorInterpolation(int fe_stokes){/*{{{*/
    470421        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    471422
    472         switch(this->element_type){
     423        switch(fe_stokes){
    473424                case XTaylorHoodEnum: return P1DGEnum;
    474                 default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     425                default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    475426        }
    476427}
    … …  
    527478                        break;
    528479                default:
    529                         _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     480                        _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    530481        }
    531482

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

    r17875 r18078  
    1212
    1313        public:
    14                 int* element_type_list; //P1CG, P1DG, MINI, P2...
    15                 int  element_type;
    16 
    1714                TriaRef();
    18                 TriaRef(const int nummodels);
    1915                ~TriaRef();
    20 
    21                 /*Management*/
    22                 void SetElementType(int type,int type_counter);
    2316
    2417                /*Numerics*/
    … …  
    2720                void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss);
    2821                void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,Gauss* gauss);
    29                 void GetNodalFunctions(IssmDouble* basis,Gauss* gauss);
    3022                void GetNodalFunctions(IssmDouble* basis,Gauss* gauss,int finiteelement);
    31                 void GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss, int index1,int index2);
    32                 void GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2);
    33                 void GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2);
    34                 void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss);
     23                void GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss, int index1,int index2,int finiteelement);
     24                void GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement);
     25                void GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement);
    3526                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
    36                 void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss);
    3727                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss,int finiteelement);
    38                 void GetInputValue(IssmDouble* pp, IssmDouble* plist, Gauss* gauss);
    3928                void GetInputValue(IssmDouble* pp, IssmDouble* plist, Gauss* gauss,int finiteelement);
    40                 void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss);
     29                void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
    4130
    4231                void NodeOnEdgeIndices(int* pnumindices,int** pindices,int index,int finiteelement);
    43                 int  NumberofNodes(void);
    4432                int  NumberofNodes(int finiteelement);
    45                 int  VelocityInterpolation(void);
    46                 int  PressureInterpolation(void);
    47                 int  TensorInterpolation(void);
     33                int  VelocityInterpolation(int fe_stokes);
     34                int  PressureInterpolation(int fe_stokes);
     35                int  TensorInterpolation(int fe_stokes);
    4836};
    4937#endif

  • issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp

    r18064 r18078  
    1717}
    1818/*}}}*/
    19 PentaInput::PentaInput(int in_enum_type,IssmDouble* in_values,int element_type_in)/*{{{*/
    20                 :PentaRef(1)
    21 {
    22 
    23         /*Set PentaRef*/
    24         this->SetElementType(element_type_in,0);
    25         this->element_type=element_type_in;
     19PentaInput::PentaInput(int in_enum_type,IssmDouble* in_values,int interpolation_type_in){/*{{{*/
    2620
    2721        /*Set Enum*/
    2822        enum_type=in_enum_type;
     23        this->interpolation_type=interpolation_type_in;
    2924
    3025        /*Set values*/
    31         this->values=xNew<IssmDouble>(this->NumberofNodes());
    32         for(int i=0;i<this->NumberofNodes();i++) values[i]=in_values[i];
     26        this->values=xNew<IssmDouble>(this->NumberofNodes(this->interpolation_type));
     27        for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) values[i]=in_values[i];
    3328}
    3429/*}}}*/
    … …  
    4641
    4742        _printf_(setw(15)<<"   PentaInput "<<setw(25)<<left<<EnumToStringx(this->enum_type)<<" [");
    48         for(int i=0;i<this->NumberofNodes();i++) _printf_(" "<<this->values[i]);
     43        for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) _printf_(" "<<this->values[i]);
    4944        _printf_("]\n");
    5045}
    … …  
    6055Object* PentaInput::copy() {/*{{{*/
    6156
    62         return new PentaInput(this->enum_type,this->values,this->element_type);
     57        return new PentaInput(this->enum_type,this->values,this->interpolation_type);
    6358
    6459}
    … …  
    7772        TriaInput* outinput=NULL;
    7873
    79         if(this->element_type==P0Enum){
     74        if(this->interpolation_type==P0Enum){
    8075                outinput=new TriaInput(this->enum_type,&this->values[0],P0Enum);
    8176        }
    … …  
    109104int  PentaInput::GetResultInterpolation(void){/*{{{*/
    110105
    111         if(this->element_type==P0Enum){
     106        if(this->interpolation_type==P0Enum){
    112107                return P0Enum;
    113108        }
    … …  
    118113int  PentaInput::GetResultNumberOfNodes(void){/*{{{*/
    119114
    120         return this->NumberofNodes();;
     115        return this->NumberofNodes(this->interpolation_type);;
    121116
    122117}
    … …  
    124119void PentaInput::ResultToPatch(IssmDouble* values,int nodesperelement,int sid){/*{{{*/
    125120
    126         int numnodes = this->NumberofNodes();
     121        int numnodes = this->NumberofNodes(this->interpolation_type);
    127122
    128123        /*Some checks*/
    … …  
    138133void PentaInput::GetInputValue(IssmDouble* pvalue){/*{{{*/
    139134
    140         if(this->element_type==P0Enum){
     135        if(this->interpolation_type==P0Enum){
    141136                pvalue=&values[0];
    142137        }
    … …  
    148143        /*Call PentaRef function*/
    149144        _assert_(gauss->Enum()==GaussPentaEnum);
    150         PentaRef::GetInputValue(pvalue,&values[0],(GaussPenta*)gauss);
     145        PentaRef::GetInputValue(pvalue,&values[0],(GaussPenta*)gauss,this->interpolation_type);
    151146
    152147}
    … …  
    156151        /*Call PentaRef function*/
    157152        _assert_(gauss->Enum()==GaussPentaEnum);
    158         PentaRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussPenta*)gauss);
     153        PentaRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussPenta*)gauss,this->interpolation_type);
    159154}
    160155/*}}}*/
    … …  
    165160void PentaInput::GetInputAverage(IssmDouble* pvalue){/*{{{*/
    166161
    167         int        numnodes  = this->NumberofNodes();
     162        int        numnodes  = this->NumberofNodes(this->interpolation_type);
    168163        IssmDouble numnodesd = reCast<int,IssmDouble>(numnodes);
    169164        IssmDouble value     = 0.;
    … …  
    179174void PentaInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/
    180175
    181         int        numnodes=this->NumberofNodes();
     176        int        numnodes=this->NumberofNodes(this->interpolation_type);
    182177        IssmDouble squaremin;
    183178
    … …  
    193188void PentaInput::ConstrainMin(IssmDouble minimum){/*{{{*/
    194189
    195         int numnodes = this->NumberofNodes();
     190        int numnodes = this->NumberofNodes(this->interpolation_type);
    196191        for(int i=0;i<numnodes;i++) if (values[i]<minimum) values[i]=minimum;
    197192}
    … …  
    201196        /*Output*/
    202197        IssmDouble norm=0.;
    203         int numnodes=this->NumberofNodes();
     198        int numnodes=this->NumberofNodes(this->interpolation_type);
    204199
    205200        for(int i=0;i<numnodes;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
    … …  
    209204IssmDouble PentaInput::Max(void){/*{{{*/
    210205
    211         int  numnodes=this->NumberofNodes();
     206        int  numnodes=this->NumberofNodes(this->interpolation_type);
    212207        IssmDouble max=values[0];
    213208
    … …  
    220215IssmDouble PentaInput::MaxAbs(void){/*{{{*/
    221216
    222         int  numnodes=this->NumberofNodes();
     217        int  numnodes=this->NumberofNodes(this->interpolation_type);
    223218        IssmDouble max=fabs(values[0]);
    224219
    … …  
    231226IssmDouble PentaInput::Min(void){/*{{{*/
    232227
    233         const int  numnodes=this->NumberofNodes();
     228        const int  numnodes=this->NumberofNodes(this->interpolation_type);
    234229        IssmDouble min=values[0];
    235230
    … …  
    242237IssmDouble PentaInput::MinAbs(void){/*{{{*/
    243238
    244         const int  numnodes=this->NumberofNodes();
     239        const int  numnodes=this->NumberofNodes(this->interpolation_type);
    245240        IssmDouble min=fabs(values[0]);
    246241
    … …  
    253248void PentaInput::Scale(IssmDouble scale_factor){/*{{{*/
    254249
    255         const int numnodes=this->NumberofNodes();
     250        const int numnodes=this->NumberofNodes(this->interpolation_type);
    256251        for(int i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
    257252}
    … …  
    259254void PentaInput::AXPY(Input* xinput,IssmDouble scalar){/*{{{*/
    260255
    261         const int numnodes=this->NumberofNodes();
     256        const int numnodes=this->NumberofNodes(this->interpolation_type);
    262257        PentaInput* xpentainput=NULL;
    263258
    … …  
    271266          _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
    272267        xpentainput=(PentaInput*)xinput;
    273         if(xpentainput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xpentainput->element_type));
     268        if(xpentainput->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xpentainput->interpolation_type));
    274269
    275270        /*Carry out the AXPY operation depending on type:*/
    … …  
    281276
    282277        int i;
    283         const int numnodes=this->NumberofNodes();
     278        const int numnodes=this->NumberofNodes(this->interpolation_type);
    284279
    285280        if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
    … …  
    290285void PentaInput::Extrude(void){/*{{{*/
    291286
    292         switch(this->element_type){
     287        switch(this->interpolation_type){
    293288                case P1Enum:
    294289                        for(int i=0;i<3;i++) this->values[3+i]=this->values[i];
    295290                        break;
    296291                default:
    297                         _error_("not supported yet for type "<<EnumToStringx(this->element_type));
     292                        _error_("not supported yet for type "<<EnumToStringx(this->interpolation_type));
    298293        }
    299294}
    … …  
    308303
    309304        /*vertically integrate depending on type (and use P1 interpolation from now on)*/
    310         switch(this->element_type){
     305        switch(this->interpolation_type){
    311306                case P1Enum:
    312307                case P1bubbleEnum:
    313308                case P2Enum:
    314309                          {
    315                                 this->element_type=P1Enum;
     310                                this->interpolation_type=P1Enum;
    316311                                GaussPenta *gauss=new GaussPenta();
    317312                                for(int iv=0;iv<3;iv++){
    … …  
    325320                          }
    326321                default:
    327                         _error_("not supported yet for type "<<EnumToStringx(this->element_type));
     322                        _error_("not supported yet for type "<<EnumToStringx(this->interpolation_type));
    328323        }
    329324}
    … …  
    336331        /*Intermediaries*/
    337332        PentaInput *xinputB  = NULL;
    338         const int   numnodes = this->NumberofNodes();
     333        const int   numnodes = this->NumberofNodes(this->interpolation_type);
    339334
    340335        /*Check that inputB is of the same type*/
    341336        if(inputB->ObjectEnum()!=PentaInputEnum)     _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
    342337        xinputB=(PentaInput*)inputB;
    343         if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
     338        if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
    344339
    345340        /*Allocate intermediary*/
    … …  
    353348
    354349        /*Create new Penta vertex input (copy of current input)*/
    355         outinput=new PentaInput(this->enum_type,AdotBvalues,this->element_type);
     350        outinput=new PentaInput(this->enum_type,AdotBvalues,this->interpolation_type);
    356351
    357352        /*Return output pointer*/
    … …  
    369364        int         i;
    370365        PentaInput  *xinputB   = NULL;
    371         const int   numnodes  = this->NumberofNodes();
     366        const int   numnodes  = this->NumberofNodes(this->interpolation_type);
    372367        IssmDouble *minvalues = xNew<IssmDouble>(numnodes);
    373368
    … …  
    375370        if(inputB->ObjectEnum()!=PentaInputEnum)       _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
    376371        xinputB=(PentaInput*)inputB;
    377         if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
     372        if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
    378373
    379374        /*Create point wise min*/
    … …  
    384379
    385380        /*Create new Penta vertex input (copy of current input)*/
    386         outinput=new PentaInput(this->enum_type,&minvalues[0],this->element_type);
     381        outinput=new PentaInput(this->enum_type,&minvalues[0],this->interpolation_type);
    387382
    388383        /*Return output pointer*/
    … …  
    399394        int         i;
    400395        PentaInput  *xinputB   = NULL;
    401         const int   numnodes  = this->NumberofNodes();
     396        const int   numnodes  = this->NumberofNodes(this->interpolation_type);
    402397        IssmDouble *maxvalues = xNew<IssmDouble>(numnodes);
    403398
    … …  
    405400        if(inputB->ObjectEnum()!=PentaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
    406401        xinputB=(PentaInput*)inputB;
    407         if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
     402        if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
    408403
    409404        /*Create point wise max*/
    … …  
    414409
    415410        /*Create new Penta vertex input (copy of current input)*/
    416         outinput=new PentaInput(this->enum_type,&maxvalues[0],this->element_type);
     411        outinput=new PentaInput(this->enum_type,&maxvalues[0],this->interpolation_type);
    417412
    418413        /*Return output pointer*/

  • issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h

    r17514 r18078  
    1717
    1818        public:
    19                 int        enum_type;
     19                int         enum_type;
     20                int         interpolation_type;
    2021                IssmDouble* values;
    2122

  • issm/trunk-jpl/src/c/classes/Inputs/SegInput.cpp

    r18064 r18078  
    1717}
    1818/*}}}*/
    19 SegInput::SegInput(int in_enum_type,IssmDouble* in_values,int element_type_in)/*{{{*/
    20         :SegRef(1)
    21 {
    22 
    23         /*Set SegRef*/
    24         this->SetElementType(element_type_in,0);
    25         this->element_type=element_type_in;
     19SegInput::SegInput(int in_enum_type,IssmDouble* in_values,int interpolation_type_in){/*{{{*/
    2620
    2721        /*Set Enum*/
    2822        enum_type=in_enum_type;
     23        this->interpolation_type=interpolation_type_in;
    2924
    3025        /*Set values*/
    31         this->values=xNew<IssmDouble>(this->NumberofNodes());
    32         for(int i=0;i<this->NumberofNodes();i++) values[i]=in_values[i];
     26        this->values=xNew<IssmDouble>(this->NumberofNodes(this->interpolation_type));
     27        for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) values[i]=in_values[i];
    3328}
    3429/*}}}*/
    … …  
    4641
    4742        _printf_(setw(15)<<"   SegInput "<<setw(25)<<left<<EnumToStringx(this->enum_type)<<" [");
    48         for(int i=0;i<this->NumberofNodes();i++) _printf_(" "<<this->values[i]);
     43        for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) _printf_(" "<<this->values[i]);
    4944        _printf_("]\n");
    5045}
    … …  
    6055Object* SegInput::copy() {/*{{{*/
    6156
    62         return new SegInput(this->enum_type,this->values,this->element_type);
     57        return new SegInput(this->enum_type,this->values,this->interpolation_type);
    6358
    6459}
    … …  
    7671void SegInput::GetInputAverage(IssmDouble* pvalue){/*{{{*/
    7772
    78         int        numnodes  = this->NumberofNodes();
     73        int        numnodes  = this->NumberofNodes(this->interpolation_type);
    7974        IssmDouble numnodesd = reCast<int,IssmDouble>(numnodes);
    8075        IssmDouble value     = 0.;
    … …  
    9085        /*Call SegRef function*/
    9186        _assert_(gauss->Enum()==GaussSegEnum);
    92         SegRef::GetInputValue(pvalue,&values[0],(GaussSeg*)gauss);
     87        SegRef::GetInputValue(pvalue,&values[0],(GaussSeg*)gauss,this->interpolation_type);
    9388
    9489}
    … …  
    9893        /*Call SegRef function*/
    9994        _assert_(gauss->Enum()==GaussSegEnum);
    100         SegRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussSeg*)gauss);
     95        SegRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussSeg*)gauss,this->interpolation_type);
    10196}
    10297/*}}}*/
    … …  
    107102IssmDouble SegInput::Min(void){/*{{{*/
    108103
    109         const int  numnodes=this->NumberofNodes();
     104        const int  numnodes=this->NumberofNodes(this->interpolation_type);
    110105        IssmDouble min=values[0];
    111106

  • issm/trunk-jpl/src/c/classes/Inputs/SegInput.h

    r17514 r18078  
    1818        public:
    1919                int         enum_type;
     20                int         interpolation_type;
    2021                IssmDouble* values;
    2122

  • issm/trunk-jpl/src/c/classes/Inputs/TetraInput.cpp

    r18064 r18078  
    1717}
    1818/*}}}*/
    19 TetraInput::TetraInput(int in_enum_type,IssmDouble* in_values,int element_type_in)/*{{{*/
    20         :TetraRef(1)
    21 {
    22 
    23         /*Set TetraRef*/
    24         this->SetElementType(element_type_in,0);
    25         this->element_type=element_type_in;
     19TetraInput::TetraInput(int in_enum_type,IssmDouble* in_values,int interpolation_type_in){/*{{{*/
    2620
    2721        /*Set Enum*/
    2822        enum_type=in_enum_type;
     23        this->interpolation_type=interpolation_type_in;
    2924
    3025        /*Set values*/
    31         this->values=xNew<IssmDouble>(this->NumberofNodes());
    32         for(int i=0;i<this->NumberofNodes();i++) values[i]=in_values[i];
     26        this->values=xNew<IssmDouble>(this->NumberofNodes(this->interpolation_type));
     27        for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) values[i]=in_values[i];
    3328}
    3429/*}}}*/
    … …  
    4641
    4742        _printf_(setw(15)<<"   TetraInput "<<setw(25)<<left<<EnumToStringx(this->enum_type)<<" [");
    48         for(int i=0;i<this->NumberofNodes();i++) _printf_(" "<<this->values[i]);
    49         _printf_("] ("<<EnumToStringx(this->element_type)<<")\n");
     43        for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) _printf_(" "<<this->values[i]);
     44        _printf_("] ("<<EnumToStringx(this->interpolation_type)<<")\n");
    5045}
    5146/*}}}*/
    … …  
    6055Object* TetraInput::copy() {/*{{{*/
    6156
    62         return new TetraInput(this->enum_type,this->values,this->element_type);
     57        return new TetraInput(this->enum_type,this->values,this->interpolation_type);
    6358
    6459}
    … …  
    7469int  TetraInput::GetResultInterpolation(void){/*{{{*/
    7570
    76         if(this->element_type==P0Enum){
     71        if(this->interpolation_type==P0Enum){
    7772                return P0Enum;
    7873        }
    … …  
    8378int  TetraInput::GetResultNumberOfNodes(void){/*{{{*/
    8479
    85         return this->NumberofNodes();
     80        return this->NumberofNodes(this->interpolation_type);
    8681
    8782}
    … …  
    8984void TetraInput::ResultToPatch(IssmDouble* values,int nodesperelement,int sid){/*{{{*/
    9085
    91         int numnodes = this->NumberofNodes();
     86        int numnodes = this->NumberofNodes(this->interpolation_type);
    9287
    9388        /*Some checks*/
    … …  
    105100        /*Call TetraRef function*/
    106101        _assert_(gauss->Enum()==GaussTetraEnum);
    107         TetraRef::GetInputValue(pvalue,&values[0],(GaussTetra*)gauss);
     102        TetraRef::GetInputValue(pvalue,&values[0],(GaussTetra*)gauss,this->interpolation_type);
    108103
    109104}
    … …  
    113108        /*Call TetraRef function*/
    114109        _assert_(gauss->Enum()==GaussTetraEnum);
    115         TetraRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussTetra*)gauss);
     110        TetraRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussTetra*)gauss,this->interpolation_type);
    116111}
    117112/*}}}*/
    … …  
    122117void TetraInput::GetInputAverage(IssmDouble* pvalue){/*{{{*/
    123118
    124         int        numnodes  = this->NumberofNodes();
     119        int        numnodes  = this->NumberofNodes(this->interpolation_type);
    125120        IssmDouble numnodesd = reCast<int,IssmDouble>(numnodes);
    126121        IssmDouble value     = 0.;
    … …  
    175170        TriaInput* outinput=NULL;
    176171
    177         if(this->element_type==P0Enum){
     172        if(this->interpolation_type==P0Enum){
    178173                outinput=new TriaInput(this->enum_type,&this->values[0],P0Enum);
    179174        }
    … …  
    204199void TetraInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/
    205200
    206         int        numnodes=this->NumberofNodes();
     201        int        numnodes=this->NumberofNodes(this->interpolation_type);
    207202        IssmDouble squaremin;
    208203
    … …  
    218213void TetraInput::ConstrainMin(IssmDouble minimum){/*{{{*/
    219214
    220         int numnodes = this->NumberofNodes();
     215        int numnodes = this->NumberofNodes(this->interpolation_type);
    221216        for(int i=0;i<numnodes;i++) if (values[i]<minimum) values[i]=minimum;
    222217}
    … …  
    226221        /*Output*/
    227222        IssmDouble norm=0.;
    228         int numnodes=this->NumberofNodes();
     223        int numnodes=this->NumberofNodes(this->interpolation_type);
    229224
    230225        for(int i=0;i<numnodes;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
    … …  
    234229IssmDouble TetraInput::Max(void){/*{{{*/
    235230
    236         int  numnodes=this->NumberofNodes();
     231        int  numnodes=this->NumberofNodes(this->interpolation_type);
    237232        IssmDouble max=values[0];
    238233
    … …  
    245240IssmDouble TetraInput::MaxAbs(void){/*{{{*/
    246241
    247         int  numnodes=this->NumberofNodes();
     242        int  numnodes=this->NumberofNodes(this->interpolation_type);
    248243        IssmDouble max=fabs(values[0]);
    249244
    … …  
    256251IssmDouble TetraInput::Min(void){/*{{{*/
    257252
    258         const int  numnodes=this->NumberofNodes();
     253        const int  numnodes=this->NumberofNodes(this->interpolation_type);
    259254        IssmDouble min=values[0];
    260255
    … …  
    267262IssmDouble TetraInput::MinAbs(void){/*{{{*/
    268263
    269         const int  numnodes=this->NumberofNodes();
     264        const int  numnodes=this->NumberofNodes(this->interpolation_type);
    270265        IssmDouble min=fabs(values[0]);
    271266
    … …  
    278273void TetraInput::Scale(IssmDouble scale_factor){/*{{{*/
    279274
    280         const int numnodes=this->NumberofNodes();
     275        const int numnodes=this->NumberofNodes(this->interpolation_type);
    281276        for(int i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
    282277}
    … …  
    284279void TetraInput::Set(IssmDouble setvalue){/*{{{*/
    285280
    286         const int numnodes=this->NumberofNodes();
     281        const int numnodes=this->NumberofNodes(this->interpolation_type);
    287282        for(int i=0;i<numnodes;i++)values[i]=setvalue;
    288283}
    … …  
    290285void TetraInput::AXPY(Input* xinput,IssmDouble scalar){/*{{{*/
    291286
    292         const int numnodes=this->NumberofNodes();
     287        const int numnodes=this->NumberofNodes(this->interpolation_type);
    293288        TetraInput*  xtriainput=NULL;
    294289
    … …  
    296291        if(xinput->ObjectEnum()!=TetraInputEnum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
    297292        xtriainput=(TetraInput*)xinput;
    298         if(xtriainput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
     293        if(xtriainput->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
    299294
    300295        /*Carry out the AXPY operation depending on type:*/
    … …  
    306301
    307302        int i;
    308         const int numnodes=this->NumberofNodes();
     303        const int numnodes=this->NumberofNodes(this->interpolation_type);
    309304
    310305        if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
    … …  
    325320        int         i;
    326321        TetraInput  *xinputB   = NULL;
    327         const int   numnodes  = this->NumberofNodes();
     322        const int   numnodes  = this->NumberofNodes(this->interpolation_type);
    328323        IssmDouble *minvalues = xNew<IssmDouble>(numnodes);
    329324
    … …  
    331326        if(inputB->ObjectEnum()!=TetraInputEnum)       _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
    332327        xinputB=(TetraInput*)inputB;
    333         if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
     328        if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
    334329
    335330        /*Create point wise min*/
    … …  
    340335
    341336        /*Create new Tetra vertex input (copy of current input)*/
    342         outinput=new TetraInput(this->enum_type,&minvalues[0],this->element_type);
     337        outinput=new TetraInput(this->enum_type,&minvalues[0],this->interpolation_type);
    343338
    344339        /*Return output pointer*/
    … …  
    356351        int         i;
    357352        TetraInput  *xinputB   = NULL;
    358         const int   numnodes  = this->NumberofNodes();
     353        const int   numnodes  = this->NumberofNodes(this->interpolation_type);
    359354        IssmDouble *maxvalues = xNew<IssmDouble>(numnodes);
    360355
    … …  
    362357        if(inputB->ObjectEnum()!=TetraInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
    363358        xinputB=(TetraInput*)inputB;
    364         if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
     359        if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
    365360
    366361        /*Create point wise max*/
    … …  
    371366
    372367        /*Create new Tetra vertex input (copy of current input)*/
    373         outinput=new TetraInput(this->enum_type,&maxvalues[0],this->element_type);
     368        outinput=new TetraInput(this->enum_type,&maxvalues[0],this->interpolation_type);
    374369
    375370        /*Return output pointer*/
    … …  
    386381        /*Intermediaries*/
    387382        TetraInput *xinputB  = NULL;
    388         const int   numnodes = this->NumberofNodes();
     383        const int   numnodes = this->NumberofNodes(this->interpolation_type);
    389384
    390385        /*Check that inputB is of the same type*/
    391386        if(inputB->ObjectEnum()!=TetraInputEnum)     _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
    392387        xinputB=(TetraInput*)inputB;
    393         if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
     388        if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
    394389
    395390        /*Allocate intermediary*/
    … …  
    403398
    404399        /*Create new Tetra vertex input (copy of current input)*/
    405         outinput=new TetraInput(this->enum_type,AdotBvalues,this->element_type);
     400        outinput=new TetraInput(this->enum_type,AdotBvalues,this->interpolation_type);
    406401
    407402        /*Return output pointer*/

  • issm/trunk-jpl/src/c/classes/Inputs/TetraInput.h

    r17514 r18078  
    1818        public:
    1919                int         enum_type;
     20                int         interpolation_type;
    2021                IssmDouble* values;
    2122

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp

    r18064 r18078  
    1717}
    1818/*}}}*/
    19 TriaInput::TriaInput(int in_enum_type,IssmDouble* in_values,int element_type_in)/*{{{*/
    20         :TriaRef(1)
    21 {
    22 
    23         /*Set TriaRef*/
    24         this->SetElementType(element_type_in,0);
    25         this->element_type=element_type_in;
     19TriaInput::TriaInput(int in_enum_type,IssmDouble* in_values,int interpolation_type_in){/*{{{*/
    2620
    2721        /*Set Enum*/
    2822        enum_type=in_enum_type;
     23        this->interpolation_type=interpolation_type_in;
    2924
    3025        /*Set values*/
    31         this->values=xNew<IssmDouble>(this->NumberofNodes());
    32         for(int i=0;i<this->NumberofNodes();i++) values[i]=in_values[i];
     26        this->values=xNew<IssmDouble>(this->NumberofNodes(this->interpolation_type));
     27        for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) values[i]=in_values[i];
    3328}
    3429/*}}}*/
    … …  
    4641
    4742        _printf_(setw(15)<<"   TriaInput "<<setw(25)<<left<<EnumToStringx(this->enum_type)<<" [");
    48         for(int i=0;i<this->NumberofNodes();i++) _printf_(" "<<this->values[i]);
    49         _printf_("] ("<<EnumToStringx(this->element_type)<<")\n");
     43        for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) _printf_(" "<<this->values[i]);
     44        _printf_("] ("<<EnumToStringx(this->interpolation_type)<<")\n");
    5045}
    5146/*}}}*/
    … …  
    6055Object* TriaInput::copy() {/*{{{*/
    6156
    62         return new TriaInput(this->enum_type,this->values,this->element_type);
     57        return new TriaInput(this->enum_type,this->values,this->interpolation_type);
    6358
    6459}
    … …  
    7873
    7974        /*Create new Tria input (copy of current input)*/
    80         outinput=new TriaInput(this->enum_type,&this->values[0],this->element_type);
     75        outinput=new TriaInput(this->enum_type,&this->values[0],this->interpolation_type);
    8176
    8277        /*Assign output*/
    … …  
    9085        SegInput* outinput=NULL;
    9186
    92         if(this->element_type==P0Enum){
     87        if(this->interpolation_type==P0Enum){
    9388                outinput=new SegInput(this->enum_type,&this->values[0],P0Enum);
    9489        }
    … …  
    112107int  TriaInput::GetResultInterpolation(void){/*{{{*/
    113108
    114         if(this->element_type==P0Enum){
     109        if(this->interpolation_type==P0Enum){
    115110                return P0Enum;
    116111        }
    … …  
    121116int  TriaInput::GetResultNumberOfNodes(void){/*{{{*/
    122117
    123         return this->NumberofNodes();
     118        return this->NumberofNodes(this->interpolation_type);
    124119
    125120}
    … …  
    127122void TriaInput::ResultToPatch(IssmDouble* values,int nodesperelement,int sid){/*{{{*/
    128123
    129         int numnodes = this->NumberofNodes();
     124        int numnodes = this->NumberofNodes(this->interpolation_type);
    130125
    131126        /*Some checks*/
    … …  
    143138        /*Call TriaRef function*/
    144139        _assert_(gauss->Enum()==GaussTriaEnum);
    145         TriaRef::GetInputValue(pvalue,&values[0],(GaussTria*)gauss);
     140        TriaRef::GetInputValue(pvalue,&values[0],(GaussTria*)gauss,this->interpolation_type);
    146141
    147142}
    … …  
    151146        /*Call TriaRef function*/
    152147        _assert_(gauss->Enum()==GaussTriaEnum);
    153         TriaRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussTria*)gauss);
     148        TriaRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussTria*)gauss,this->interpolation_type);
    154149}
    155150/*}}}*/
    … …  
    160155void TriaInput::GetInputAverage(IssmDouble* pvalue){/*{{{*/
    161156
    162         int        numnodes  = this->NumberofNodes();
     157        int        numnodes  = this->NumberofNodes(this->interpolation_type);
    163158        IssmDouble numnodesd = reCast<int,IssmDouble>(numnodes);
    164159        IssmDouble value     = 0.;
    … …  
    212207void TriaInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/
    213208
    214         int        numnodes=this->NumberofNodes();
     209        int        numnodes=this->NumberofNodes(this->interpolation_type);
    215210        IssmDouble squaremin;
    216211
    … …  
    226221void TriaInput::ConstrainMin(IssmDouble minimum){/*{{{*/
    227222
    228         int numnodes = this->NumberofNodes();
     223        int numnodes = this->NumberofNodes(this->interpolation_type);
    229224        for(int i=0;i<numnodes;i++) if (values[i]<minimum) values[i]=minimum;
    230225}
    … …  
    234229        /*Output*/
    235230        IssmDouble norm=0.;
    236         int numnodes=this->NumberofNodes();
     231        int numnodes=this->NumberofNodes(this->interpolation_type);
    237232
    238233        for(int i=0;i<numnodes;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
    … …  
    242237IssmDouble TriaInput::Max(void){/*{{{*/
    243238
    244         int  numnodes=this->NumberofNodes();
     239        int  numnodes=this->NumberofNodes(this->interpolation_type);
    245240        IssmDouble max=values[0];
    246241
    … …  
    253248IssmDouble TriaInput::MaxAbs(void){/*{{{*/
    254249
    255         int  numnodes=this->NumberofNodes();
     250        int  numnodes=this->NumberofNodes(this->interpolation_type);
    256251        IssmDouble max=fabs(values[0]);
    257252
    … …  
    264259IssmDouble TriaInput::Min(void){/*{{{*/
    265260
    266         const int  numnodes=this->NumberofNodes();
     261        const int  numnodes=this->NumberofNodes(this->interpolation_type);
    267262        IssmDouble min=values[0];
    268263
    … …  
    275270IssmDouble TriaInput::MinAbs(void){/*{{{*/
    276271
    277         const int  numnodes=this->NumberofNodes();
     272        const int  numnodes=this->NumberofNodes(this->interpolation_type);
    278273        IssmDouble min=fabs(values[0]);
    279274
    … …  
    286281void TriaInput::Scale(IssmDouble scale_factor){/*{{{*/
    287282
    288         const int numnodes=this->NumberofNodes();
     283        const int numnodes=this->NumberofNodes(this->interpolation_type);
    289284        for(int i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
    290285}
    … …  
    292287void TriaInput::Set(IssmDouble setvalue){/*{{{*/
    293288
    294         const int numnodes=this->NumberofNodes();
     289        const int numnodes=this->NumberofNodes(this->interpolation_type);
    295290        for(int i=0;i<numnodes;i++)values[i]=setvalue;
    296291}
    … …  
    298293void TriaInput::AXPY(Input* xinput,IssmDouble scalar){/*{{{*/
    299294
    300         const int numnodes=this->NumberofNodes();
     295        const int numnodes=this->NumberofNodes(this->interpolation_type);
    301296        TriaInput*  xtriainput=NULL;
    302297
    … …  
    304299        if(xinput->ObjectEnum()!=TriaInputEnum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
    305300        xtriainput=(TriaInput*)xinput;
    306         if(xtriainput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
     301        if(xtriainput->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
    307302
    308303        /*Carry out the AXPY operation depending on type:*/
    … …  
    314309
    315310        int i;
    316         const int numnodes=this->NumberofNodes();
     311        const int numnodes=this->NumberofNodes(this->interpolation_type);
    317312
    318313        if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
    … …  
    333328        int         i;
    334329        TriaInput  *xinputB   = NULL;
    335         const int   numnodes  = this->NumberofNodes();
     330        const int   numnodes  = this->NumberofNodes(this->interpolation_type);
    336331        IssmDouble *minvalues = xNew<IssmDouble>(numnodes);
    337332
    … …  
    339334        if(inputB->ObjectEnum()!=TriaInputEnum)       _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
    340335        xinputB=(TriaInput*)inputB;
    341         if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
     336        if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
    342337
    343338        /*Create point wise min*/
    … …  
    348343
    349344        /*Create new Tria vertex input (copy of current input)*/
    350         outinput=new TriaInput(this->enum_type,&minvalues[0],this->element_type);
     345        outinput=new TriaInput(this->enum_type,&minvalues[0],this->interpolation_type);
    351346
    352347        /*Return output pointer*/
    … …  
    364359        int         i;
    365360        TriaInput  *xinputB   = NULL;
    366         const int   numnodes  = this->NumberofNodes();
     361        const int   numnodes  = this->NumberofNodes(this->interpolation_type);
    367362        IssmDouble *maxvalues = xNew<IssmDouble>(numnodes);
    368363
    … …  
    370365        if(inputB->ObjectEnum()!=TriaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
    371366        xinputB=(TriaInput*)inputB;
    372         if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
     367        if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
    373368
    374369        /*Create point wise max*/
    … …  
    379374
    380375        /*Create new Tria vertex input (copy of current input)*/
    381         outinput=new TriaInput(this->enum_type,&maxvalues[0],this->element_type);
     376        outinput=new TriaInput(this->enum_type,&maxvalues[0],this->interpolation_type);
    382377
    383378        /*Return output pointer*/
    … …  
    394389        /*Intermediaries*/
    395390        TriaInput *xinputB  = NULL;
    396         const int   numnodes = this->NumberofNodes();
     391        const int   numnodes = this->NumberofNodes(this->interpolation_type);
    397392
    398393        /*Check that inputB is of the same type*/
    399394        if(inputB->ObjectEnum()!=TriaInputEnum)     _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
    400395        xinputB=(TriaInput*)inputB;
    401         if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
     396        if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
    402397
    403398        /*Allocate intermediary*/
    … …  
    411406
    412407        /*Create new Tria vertex input (copy of current input)*/
    413         outinput=new TriaInput(this->enum_type,AdotBvalues,this->element_type);
     408        outinput=new TriaInput(this->enum_type,AdotBvalues,this->interpolation_type);
    414409
    415410        /*Return output pointer*/

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h

    r17514 r18078  
    1818        public:
    1919                int         enum_type;
     20                int         interpolation_type;
    2021                IssmDouble* values;
    2122

  • issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp

    r18064 r18078  
    451451                gauss->GaussPoint(ig);
    452452
    453                 tria->GetSegmentBFlux(&B[0],gauss,index1,index2);
    454                 tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2);
     453                tria->GetSegmentBFlux(&B[0],gauss,index1,index2,tria->FiniteElement());
     454                tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2,tria->FiniteElement());
    455455
    456456                vxaverage_input->GetInputValue(&vx,gauss);
    … …  
    529529                gauss->GaussPoint(ig);
    530530
    531                 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
     531                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
    532532
    533533                vxaverage_input->GetInputValue(&vx,gauss);
    … …  
    597597                gauss->GaussPoint(ig);
    598598
    599                 tria->GetSegmentBFlux(&B[0],gauss,index1,index2);
    600                 tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2);
     599                tria->GetSegmentBFlux(&B[0],gauss,index1,index2,tria->FiniteElement());
     600                tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2,tria->FiniteElement());
    601601
    602602                vxaverage_input->GetInputValue(&vx,gauss);
    … …  
    674674                gauss->GaussPoint(ig);
    675675
    676                 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
     676                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
    677677
    678678                vxaverage_input->GetInputValue(&vx,gauss);
    … …  
    790790                gauss->GaussPoint(ig);
    791791
    792                 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
     792                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
    793793
    794794                vxaverage_input->GetInputValue(&vx,gauss);
    … …  
    876876                gauss->GaussPoint(ig);
    877877
    878                 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
     878                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
    879879
    880880                vxaverage_input->GetInputValue(&vx,gauss);
Note: See TracChangeset for help on using the changeset viewer.