Changeset 18926


Ignore:
Timestamp:
12/04/14 09:21:03 (10 years ago)
Author:
seroussi
Message:

CHG: reordering loads

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

Legend:

Unmodified
Added
Removed

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

    r18805 r18926  
    4040}
    4141/*}}}*/
    42 void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
    43 
    44         switch(this->law){
    45                 case 1:
    46                         GetAlpha2Viscous(palpha2,gauss);
    47                         break;
    48                 case 2:
    49                         GetAlpha2Weertman(palpha2,gauss);
    50                         break;
    51                 case 3:
    52                         GetAlpha2Hydro(palpha2,gauss);
    53                         break;
    54                 case 4:
    55                         GetAlpha2Temp(palpha2,gauss);
    56                         break;
    57                 case 5:
    58                         GetAlpha2WaterLayer(palpha2,gauss);
    59                         break;
    60                 case 6:
    61                         GetAlpha2WeertmanTemp(palpha2,gauss);
    62                         break;
    63           default:
    64                         _error_("not supported");
    65         }
    66 
    67 }/*}}}*/
    68 void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
    69 
    70         /*This routine calculates the basal friction coefficient
    71           alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
     42void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
     43
     44        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.
     45         * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
     46         * alpha_complement= Neff ^r * vel ^s*/
     47
     48        if(this->law!=1)_error_("not supported");
    7249
    7350        /*diverse: */
    7451        IssmDouble  r,s;
    75         IssmDouble  drag_p, drag_q;
     52        IssmDouble  vx,vy,vz,vmag;
     53        IssmDouble  drag_p,drag_q;
    7654        IssmDouble  Neff;
    77         IssmDouble  thickness,bed;
    78         IssmDouble  vx,vy,vz,vmag;
    7955        IssmDouble  drag_coefficient;
    80         IssmDouble  alpha2;
     56        IssmDouble  bed,thickness;
     57        IssmDouble  alpha_complement;
    8158
    8259        /*Recover parameters: */
    … …  
    9875        if(Neff<0)Neff=0;
    9976
     77        //We need the velocity magnitude to evaluate the basal stress:
    10078        switch(dim){
    10179                case 1:
    … …  
    11997
    12098        /*Check to prevent dividing by zero if vmag==0*/
    121         if(vmag==0. && (s-1.)<0.) alpha2=0.;
    122         else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
    123         _assert_(!xIsNan<IssmDouble>(alpha2));
    124 
    125         /*Assign output pointers:*/
    126         *palpha2=alpha2;
    127 }/*}}}*/
    128 void Friction::GetAlpha2Weertman(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
    129 
    130         /*This routine calculates the basal friction coefficient alpha2= C^-1/m |v|^(1/m-1) */
    131 
    132         /*diverse: */
    133         IssmDouble  C,m;
    134         IssmDouble  vx,vy,vz,vmag;
    135         IssmDouble  alpha2;
    136 
    137         /*Recover parameters: */
    138         element->GetInputValue(&C,gauss,FrictionCEnum);
    139         element->GetInputValue(&m,FrictionMEnum);
    140 
    141         switch(dim){
    142                 case 1:
    143                         element->GetInputValue(&vx,gauss,VxEnum);
    144                         vmag=sqrt(vx*vx);
    145                         break;
    146                 case 2:
    147                         element->GetInputValue(&vx,gauss,VxEnum);
    148                         element->GetInputValue(&vy,gauss,VyEnum);
    149                         vmag=sqrt(vx*vx+vy*vy);
    150                         break;
    151                 case 3:
    152                         element->GetInputValue(&vx,gauss,VxEnum);
    153                         element->GetInputValue(&vy,gauss,VyEnum);
    154                         element->GetInputValue(&vz,gauss,VzEnum);
    155                         vmag=sqrt(vx*vx+vy*vy+vz*vz);
    156                         break;
    157                 default:
    158                         _error_("not supported");
    159         }
    160 
    161         /*Check to prevent dividing by zero if vmag==0*/
    162         if(vmag==0. && (1./m-1.)<0.) alpha2=0.;
    163         else alpha2=pow(C,-1./m)*pow(vmag,(1./m-1.));
    164         _assert_(!xIsNan<IssmDouble>(alpha2));
    165 
    166         /*Assign output pointers:*/
    167         *palpha2=alpha2;
     99        if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
     100        else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));_assert_(!xIsNan<IssmDouble>(alpha_complement));
     101
     102        /*Assign output pointers:*/
     103        *palpha_complement=alpha_complement;
     104}
     105/*}}}*/
     106void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
     107
     108        switch(this->law){
     109                case 1:
     110                        GetAlpha2Viscous(palpha2,gauss);
     111                        break;
     112                case 2:
     113                        GetAlpha2Weertman(palpha2,gauss);
     114                        break;
     115                case 3:
     116                        GetAlpha2Hydro(palpha2,gauss);
     117                        break;
     118                case 4:
     119                        GetAlpha2Temp(palpha2,gauss);
     120                        break;
     121                case 5:
     122                        GetAlpha2WaterLayer(palpha2,gauss);
     123                        break;
     124                case 6:
     125                        GetAlpha2WeertmanTemp(palpha2,gauss);
     126                        break;
     127          default:
     128                        _error_("not supported");
     129        }
     130
    168131}/*}}}*/
    169132void Friction::GetAlpha2Hydro(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
    … …  
    256219        element->parameters->FindParam(&gamma,FrictionGammaEnum);
    257220        alpha2 = alpha2 / exp((T-Tpmp)/gamma);
     221
     222        /*Assign output pointers:*/
     223        *palpha2=alpha2;
     224}/*}}}*/
     225void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
     226
     227        /*This routine calculates the basal friction coefficient
     228          alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
     229
     230        /*diverse: */
     231        IssmDouble  r,s;
     232        IssmDouble  drag_p, drag_q;
     233        IssmDouble  Neff;
     234        IssmDouble  thickness,bed;
     235        IssmDouble  vx,vy,vz,vmag;
     236        IssmDouble  drag_coefficient;
     237        IssmDouble  alpha2;
     238
     239        /*Recover parameters: */
     240        element->GetInputValue(&drag_p,FrictionPEnum);
     241        element->GetInputValue(&drag_q,FrictionQEnum);
     242        element->GetInputValue(&thickness, gauss,ThicknessEnum);
     243        element->GetInputValue(&bed, gauss,BaseEnum);
     244        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
     245        IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     246        IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
     247        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
     248
     249        //compute r and q coefficients: */
     250        r=drag_q/drag_p;
     251        s=1./drag_p;
     252
     253        //From bed and thickness, compute effective pressure when drag is viscous:
     254        Neff=gravity*(rho_ice*thickness+rho_water*bed);
     255        if(Neff<0)Neff=0;
     256
     257        switch(dim){
     258                case 1:
     259                        element->GetInputValue(&vx,gauss,VxEnum);
     260                        vmag=sqrt(vx*vx);
     261                        break;
     262                case 2:
     263                        element->GetInputValue(&vx,gauss,VxEnum);
     264                        element->GetInputValue(&vy,gauss,VyEnum);
     265                        vmag=sqrt(vx*vx+vy*vy);
     266                        break;
     267                case 3:
     268                        element->GetInputValue(&vx,gauss,VxEnum);
     269                        element->GetInputValue(&vy,gauss,VyEnum);
     270                        element->GetInputValue(&vz,gauss,VzEnum);
     271                        vmag=sqrt(vx*vx+vy*vy+vz*vz);
     272                        break;
     273                default:
     274                        _error_("not supported");
     275        }
     276
     277        /*Check to prevent dividing by zero if vmag==0*/
     278        if(vmag==0. && (s-1.)<0.) alpha2=0.;
     279        else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
     280        _assert_(!xIsNan<IssmDouble>(alpha2));
    258281
    259282        /*Assign output pointers:*/
    … …  
    322345        *palpha2=alpha2;
    323346}/*}}}*/
     347void Friction::GetAlpha2Weertman(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
     348
     349        /*This routine calculates the basal friction coefficient alpha2= C^-1/m |v|^(1/m-1) */
     350
     351        /*diverse: */
     352        IssmDouble  C,m;
     353        IssmDouble  vx,vy,vz,vmag;
     354        IssmDouble  alpha2;
     355
     356        /*Recover parameters: */
     357        element->GetInputValue(&C,gauss,FrictionCEnum);
     358        element->GetInputValue(&m,FrictionMEnum);
     359
     360        switch(dim){
     361                case 1:
     362                        element->GetInputValue(&vx,gauss,VxEnum);
     363                        vmag=sqrt(vx*vx);
     364                        break;
     365                case 2:
     366                        element->GetInputValue(&vx,gauss,VxEnum);
     367                        element->GetInputValue(&vy,gauss,VyEnum);
     368                        vmag=sqrt(vx*vx+vy*vy);
     369                        break;
     370                case 3:
     371                        element->GetInputValue(&vx,gauss,VxEnum);
     372                        element->GetInputValue(&vy,gauss,VyEnum);
     373                        element->GetInputValue(&vz,gauss,VzEnum);
     374                        vmag=sqrt(vx*vx+vy*vy+vz*vz);
     375                        break;
     376                default:
     377                        _error_("not supported");
     378        }
     379
     380        /*Check to prevent dividing by zero if vmag==0*/
     381        if(vmag==0. && (1./m-1.)<0.) alpha2=0.;
     382        else alpha2=pow(C,-1./m)*pow(vmag,(1./m-1.));
     383        _assert_(!xIsNan<IssmDouble>(alpha2));
     384
     385        /*Assign output pointers:*/
     386        *palpha2=alpha2;
     387}/*}}}*/
    324388void Friction::GetAlpha2WeertmanTemp(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
    325389        /*Here, we want to parameterize the friction as a function of temperature
    … …  
    349413        *palpha2=alpha2;
    350414}/*}}}*/
    351 void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
    352 
    353         /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.
    354          * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
    355          * alpha_complement= Neff ^r * vel ^s*/
    356 
    357         if(this->law!=1)_error_("not supported");
    358 
    359         /*diverse: */
    360         IssmDouble  r,s;
    361         IssmDouble  vx,vy,vz,vmag;
    362         IssmDouble  drag_p,drag_q;
    363         IssmDouble  Neff;
    364         IssmDouble  drag_coefficient;
    365         IssmDouble  bed,thickness;
    366         IssmDouble  alpha_complement;
    367 
    368         /*Recover parameters: */
    369         element->GetInputValue(&drag_p,FrictionPEnum);
    370         element->GetInputValue(&drag_q,FrictionQEnum);
    371         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    372         element->GetInputValue(&bed, gauss,BaseEnum);
    373         element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    374         IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    375         IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
    376         IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    377 
    378         //compute r and q coefficients: */
    379         r=drag_q/drag_p;
    380         s=1./drag_p;
    381 
    382         //From bed and thickness, compute effective pressure when drag is viscous:
    383         Neff=gravity*(rho_ice*thickness+rho_water*bed);
    384         if(Neff<0)Neff=0;
    385 
    386         //We need the velocity magnitude to evaluate the basal stress:
    387         switch(dim){
    388                 case 1:
    389                         element->GetInputValue(&vx,gauss,VxEnum);
    390                         vmag=sqrt(vx*vx);
    391                         break;
    392                 case 2:
    393                         element->GetInputValue(&vx,gauss,VxEnum);
    394                         element->GetInputValue(&vy,gauss,VyEnum);
    395                         vmag=sqrt(vx*vx+vy*vy);
    396                         break;
    397                 case 3:
    398                         element->GetInputValue(&vx,gauss,VxEnum);
    399                         element->GetInputValue(&vy,gauss,VyEnum);
    400                         element->GetInputValue(&vz,gauss,VzEnum);
    401                         vmag=sqrt(vx*vx+vy*vy+vz*vz);
    402                         break;
    403                 default:
    404                         _error_("not supported");
    405         }
    406 
    407         /*Check to prevent dividing by zero if vmag==0*/
    408         if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
    409         else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));_assert_(!xIsNan<IssmDouble>(alpha_complement));
    410 
    411         /*Assign output pointers:*/
    412         *palpha_complement=alpha_complement;
    413 }
    414 /*}}}*/

  • issm/trunk-jpl/src/c/classes/Loads/Friction.h

    r18804 r18926  
    2929
    3030                void  Echo(void);
     31                void  GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss);
    3132                void  GetAlpha2(IssmDouble* palpha2,Gauss* gauss);
    32                 void  GetAlpha2Viscous(IssmDouble* palpha2,Gauss* gauss);
    33                 void  GetAlpha2Weertman(IssmDouble* palpha2,Gauss* gauss);
    3433                void  GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss);
    3534                void  GetAlpha2Temp(IssmDouble* palpha2,Gauss* gauss);
     35                void  GetAlpha2Viscous(IssmDouble* palpha2,Gauss* gauss);
    3636                void  GetAlpha2WaterLayer(IssmDouble* palpha2,Gauss* gauss);
     37                void  GetAlpha2Weertman(IssmDouble* palpha2,Gauss* gauss);
    3738                void  GetAlpha2WeertmanTemp(IssmDouble* palpha2,Gauss* gauss);
    38                 void  GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss);
    3939};
    4040

  • issm/trunk-jpl/src/c/classes/Loads/Load.h

    r18237 r18926  
    2626                virtual       ~Load(){};
    2727                virtual void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
    28                 virtual void  ResetHooks()=0;
    29                 virtual bool  IsPenalty(void)=0;
    30                 virtual int   GetNumberOfNodes(void)=0;
    31                 virtual void  GetNodesSidList(int* sidlist)=0;
    32                 virtual void  GetNodesLidList(int* lidlist)=0;
    33                 virtual void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
     28                virtual void  CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0;
    3429                virtual void  CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs)=0;
    3530                virtual void  CreatePVector(Vector<IssmDouble>* pf)=0;
    36                 virtual void  CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0;
     31                virtual void  GetNodesLidList(int* lidlist)=0;
     32                virtual void  GetNodesSidList(int* sidlist)=0;
     33                virtual int   GetNumberOfNodes(void)=0;
     34                virtual bool  InAnalysis(int analysis_type)=0;
     35                virtual bool  IsPenalty(void)=0;
    3736                virtual void  PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax)=0;
    3837                virtual void  PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs, IssmDouble kmax)=0;
    3938                virtual void  PenaltyCreatePVector(Vector<IssmDouble>* pf, IssmDouble kmax)=0;
    40                 virtual bool  InAnalysis(int analysis_type)=0;
     39                virtual void  ResetHooks()=0;
     40                virtual void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
    4141                virtual void  SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;
    4242};

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

    r18521 r18926  
    4848}
    4949/*}}}*/
    50 void Loads::ResetHooks(){/*{{{*/
    51 
    52         vector<Object*>::iterator object;
    53         Load* load=NULL;
    54 
    55         for ( object=objects.begin() ; object < objects.end(); object++ ){
    56 
    57                 load=xDynamicCast<Load*>((*object));
    58                 load->ResetHooks();
    59 
    60         }
    61 
    62 }
    63 /*}}}*/
    6450bool Loads::IsPenalty(int analysis_type){/*{{{*/
    6551
    … …  
    8672}
    8773/*}}}*/
    88 int Loads::MaxNumNodes(int analysis_type){/*{{{*/
     74int  Loads::MaxNumNodes(int analysis_type){/*{{{*/
    8975
    9076        int max=0;
    … …  
    10995}
    11096/*}}}*/
    111 int Loads::NumberOfLoads(void){/*{{{*/
     97int  Loads::NumberOfLoads(void){/*{{{*/
    11298
    11399        int localloads;
    … …  
    124110}
    125111/*}}}*/
    126 int Loads::NumberOfLoads(int analysis_type){/*{{{*/
     112int  Loads::NumberOfLoads(int analysis_type){/*{{{*/
    127113
    128114        int localloads = 0;
    … …  
    145131}
    146132/*}}}*/
    147 int Loads::Size(void){/*{{{*/
     133void Loads::ResetHooks(){/*{{{*/
    148134
    149         return this->DataSet::Size();
    150 }
    151 /*}}}*/
    152 int Loads::Size(int analysis_type){/*{{{*/
     135        vector<Object*>::iterator object;
     136        Load* load=NULL;
    153137
    154         int localloads = 0;
     138        for ( object=objects.begin() ; object < objects.end(); object++ ){
    155139
    156         /*Get number of local loads*/
    157         for(int i=0;i<this->Size();i++){
     140                load=xDynamicCast<Load*>((*object));
     141                load->ResetHooks();
    158142
    159                 Load* load=xDynamicCast<Load*>(this->GetObjectByOffset(i));
    160 
    161                 /*Check that this load corresponds to our analysis currently being carried out: */
    162                 if (load->InAnalysis(analysis_type)) localloads++;
    163143        }
    164144
    165         return localloads;
    166145}
    167146/*}}}*/
    … …  
    180159}
    181160/*}}}*/
     161int  Loads::Size(void){/*{{{*/
     162
     163        return this->DataSet::Size();
     164}
     165/*}}}*/
     166int  Loads::Size(int analysis_type){/*{{{*/
     167
     168        int localloads = 0;
     169
     170        /*Get number of local loads*/
     171        for(int i=0;i<this->Size();i++){
     172
     173                Load* load=xDynamicCast<Load*>(this->GetObjectByOffset(i));
     174
     175                /*Check that this load corresponds to our analysis currently being carried out: */
     176                if (load->InAnalysis(analysis_type)) localloads++;
     177        }
     178
     179        return localloads;
     180}
     181/*}}}*/

  • issm/trunk-jpl/src/c/classes/Loads/Loads.h

    r18237 r18926  
    2424                /*numerics*/
    2525                void  Configure(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters);
    26                 void  ResetHooks();
    2726                bool  IsPenalty(int analysis);
    2827                int   MaxNumNodes(int analysis);
    2928                int   NumberOfLoads(void);
    3029                int   NumberOfLoads(int analysis);
     30                void  ResetHooks();
    3131                void  SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters);
    3232                int   Size(int analysis);

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

    r18237 r18926  
    131131
    132132/*Object virtual functions definitions:*/
    133 void Numericalflux::Echo(void){/*{{{*/
    134         _printf_("Numericalflux:\n");
    135         _printf_("   id: " << id << "\n");
    136         _printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
    137         _printf_("   flux_type: " << this->flux_type<< "\n");
    138         hnodes->Echo();
    139         hvertices->Echo();
    140         helement->Echo();
    141         _printf_("   parameters: " << parameters << "\n");
    142 }
    143 /*}}}*/
    144 void Numericalflux::DeepEcho(void){/*{{{*/
     133Object* Numericalflux::copy() {/*{{{*/
     134
     135        Numericalflux* numericalflux=NULL;
     136
     137        numericalflux=new Numericalflux();
     138
     139        /*copy fields: */
     140        numericalflux->id=this->id;
     141        numericalflux->analysis_type=this->analysis_type;
     142        numericalflux->flux_type=this->flux_type;
     143
     144        /*point parameters: */
     145        numericalflux->parameters=this->parameters;
     146
     147        /*now deal with hooks and objects: */
     148        numericalflux->hnodes    = (Hook*)this->hnodes->copy();
     149        numericalflux->hvertices = (Hook*)this->hvertices->copy();
     150        numericalflux->helement  = (Hook*)this->helement->copy();
     151
     152        /*corresponding fields*/
     153        numericalflux->nodes    = (Node**)numericalflux->hnodes->deliverp();
     154        numericalflux->vertices = (Vertex**)numericalflux->hvertices->deliverp();
     155        numericalflux->element  = (Element*)numericalflux->helement->delivers();
     156
     157        return numericalflux;
     158}
     159/*}}}*/
     160void    Numericalflux::DeepEcho(void){/*{{{*/
    145161
    146162        _printf_("Numericalflux:\n");
    … …  
    158174}               
    159175/*}}}*/
    160 int Numericalflux::Id(void){/*{{{*/
     176void    Numericalflux::Echo(void){/*{{{*/
     177        _printf_("Numericalflux:\n");
     178        _printf_("   id: " << id << "\n");
     179        _printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
     180        _printf_("   flux_type: " << this->flux_type<< "\n");
     181        hnodes->Echo();
     182        hvertices->Echo();
     183        helement->Echo();
     184        _printf_("   parameters: " << parameters << "\n");
     185}
     186/*}}}*/
     187int     Numericalflux::Id(void){/*{{{*/
    161188        return id;
    162189}
    163190/*}}}*/
    164 int Numericalflux::ObjectEnum(void){/*{{{*/
     191int     Numericalflux::ObjectEnum(void){/*{{{*/
    165192
    166193        return NumericalfluxEnum;
    167194
    168 }
    169 /*}}}*/
    170 Object* Numericalflux::copy() {/*{{{*/
    171 
    172         Numericalflux* numericalflux=NULL;
    173 
    174         numericalflux=new Numericalflux();
    175 
    176         /*copy fields: */
    177         numericalflux->id=this->id;
    178         numericalflux->analysis_type=this->analysis_type;
    179         numericalflux->flux_type=this->flux_type;
    180 
    181         /*point parameters: */
    182         numericalflux->parameters=this->parameters;
    183 
    184         /*now deal with hooks and objects: */
    185         numericalflux->hnodes    = (Hook*)this->hnodes->copy();
    186         numericalflux->hvertices = (Hook*)this->hvertices->copy();
    187         numericalflux->helement  = (Hook*)this->helement->copy();
    188 
    189         /*corresponding fields*/
    190         numericalflux->nodes    = (Node**)numericalflux->hnodes->deliverp();
    191         numericalflux->vertices = (Vertex**)numericalflux->hvertices->deliverp();
    192         numericalflux->element  = (Element*)numericalflux->helement->delivers();
    193 
    194         return numericalflux;
    195195}
    196196/*}}}*/
    … …  
    212212        /*point parameters to real dataset: */
    213213        this->parameters=parametersin;
    214 }
    215 /*}}}*/
    216 void  Numericalflux::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
    217 
    218 }
    219 /*}}}*/
    220 void  Numericalflux::ResetHooks(){/*{{{*/
    221 
    222         this->nodes=NULL;
    223         this->vertices=NULL;
    224         this->element=NULL;
    225         this->parameters=NULL;
    226 
    227         /*Get Element type*/
    228         this->hnodes->reset();
    229         this->hvertices->reset();
    230         this->helement->reset();
    231 
    232214}
    233215/*}}}*/
    … …  
    291273}
    292274/*}}}*/
     275void  Numericalflux::GetNodesLidList(int* lidlist){/*{{{*/
     276
     277        _assert_(lidlist);
     278        _assert_(nodes);
     279
     280        switch(this->flux_type){
     281                case InternalEnum:
     282                        for(int i=0;i<NUMNODES_INTERNAL;i++) lidlist[i]=nodes[i]->Lid();
     283                        return;
     284                case BoundaryEnum:
     285                        for(int i=0;i<NUMNODES_BOUNDARY;i++) lidlist[i]=nodes[i]->Lid();
     286                        return;
     287                default:
     288                        _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
     289        }
     290}
     291/*}}}*/
    293292void  Numericalflux::GetNodesSidList(int* sidlist){/*{{{*/
    294293
    … …  
    308307}
    309308/*}}}*/
    310 void  Numericalflux::GetNodesLidList(int* lidlist){/*{{{*/
    311 
    312         _assert_(lidlist);
    313         _assert_(nodes);
    314 
    315         switch(this->flux_type){
    316                 case InternalEnum:
    317                         for(int i=0;i<NUMNODES_INTERNAL;i++) lidlist[i]=nodes[i]->Lid();
    318                         return;
    319                 case BoundaryEnum:
    320                         for(int i=0;i<NUMNODES_BOUNDARY;i++) lidlist[i]=nodes[i]->Lid();
    321                         return;
    322                 default:
    323                         _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
    324         }
    325 }
    326 /*}}}*/
    327309int   Numericalflux::GetNumberOfNodes(void){/*{{{*/
    328310
    … …  
    338320}
    339321/*}}}*/
    340 bool  Numericalflux::IsPenalty(void){/*{{{*/
    341         return false;
    342 }
    343 /*}}}*/
    344 void  Numericalflux::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){/*{{{*/
    345 
    346         /*No stiffness loads applied, do nothing: */
    347         return;
    348 
    349 }
    350 /*}}}*/
    351 void  Numericalflux::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){/*{{{*/
    352 
    353         /*No penalty loads applied, do nothing: */
    354         return;
    355 
    356 }
    357 /*}}}*/
    358322bool  Numericalflux::InAnalysis(int in_analysis_type){/*{{{*/
    359323        if (in_analysis_type==this->analysis_type) return true;
    360324        else return false;
     325}
     326/*}}}*/
     327bool  Numericalflux::IsPenalty(void){/*{{{*/
     328        return false;
     329}
     330/*}}}*/
     331void  Numericalflux::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){/*{{{*/
     332
     333        /*No stiffness loads applied, do nothing: */
     334        return;
     335
     336}
     337/*}}}*/
     338void  Numericalflux::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){/*{{{*/
     339
     340        /*No penalty loads applied, do nothing: */
     341        return;
     342
     343}
     344/*}}}*/
     345void  Numericalflux::ResetHooks(){/*{{{*/
     346
     347        this->nodes=NULL;
     348        this->vertices=NULL;
     349        this->element=NULL;
     350        this->parameters=NULL;
     351
     352        /*Get Element type*/
     353        this->hnodes->reset();
     354        this->hvertices->reset();
     355        this->helement->reset();
     356
     357}
     358/*}}}*/
     359void  Numericalflux::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
     360
    361361}
    362362/*}}}*/
    … …  
    417417
    418418/*Numericalflux management*/
     419ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethickness(void){/*{{{*/
     420
     421        switch(this->flux_type){
     422                case InternalEnum:
     423                        return CreateKMatrixAdjointBalancethicknessInternal();
     424                case BoundaryEnum:
     425                        return CreateKMatrixAdjointBalancethicknessBoundary();
     426                default:
     427                        _error_("type not supported yet");
     428        }
     429}
     430/*}}}*/
     431ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessBoundary(void){/*{{{*/
     432
     433        ElementMatrix* Ke=CreateKMatrixBalancethicknessBoundary();
     434        if(Ke) Ke->Transpose();
     435        return Ke;
     436}
     437/*}}}*/
     438ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessInternal(void){/*{{{*/
     439
     440        ElementMatrix* Ke=CreateKMatrixBalancethicknessInternal();
     441        if (Ke) Ke->Transpose();
     442        return Ke;
     443}
     444/*}}}*/
     445ElementMatrix* Numericalflux::CreateKMatrixBalancethickness(void){/*{{{*/
     446
     447        switch(this->flux_type){
     448                case InternalEnum:
     449                        return CreateKMatrixBalancethicknessInternal();
     450                case BoundaryEnum:
     451                        return CreateKMatrixBalancethicknessBoundary();
     452                default:
     453                        _error_("type not supported yet");
     454        }
     455}
     456/*}}}*/
     457ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessBoundary(void){/*{{{*/
     458
     459        /* constants*/
     460        const int numdof=NDOF1*NUMNODES_BOUNDARY;
     461
     462        /* Intermediaries*/
     463        int        i,j,ig,index1,index2;
     464        IssmDouble     DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN;
     465        IssmDouble     xyz_list[NUMVERTICES][3];
     466        IssmDouble     normal[2];
     467        IssmDouble     L[numdof];
     468        IssmDouble     Ke_g[numdof][numdof];
     469        GaussTria *gauss;
     470
     471        /*Initialize Element matrix and return if necessary*/
     472        ElementMatrix* Ke = NULL;
     473        Tria*  tria=(Tria*)element;
     474        if(!tria->IsIceInElement()) return NULL;
     475
     476        /*Retrieve all inputs and parameters*/
     477        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     478        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
     479        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
     480        GetNormal(&normal[0],xyz_list);
     481
     482        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
     483        index1=tria->GetNodeIndex(nodes[0]);
     484        index2=tria->GetNodeIndex(nodes[1]);
     485
     486        gauss=new GaussTria();
     487        gauss->GaussEdgeCenter(index1,index2);
     488        vxaverage_input->GetInputValue(&mean_vx,gauss);
     489        vyaverage_input->GetInputValue(&mean_vy,gauss);
     490        delete gauss;
     491
     492        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
     493        if (UdotN<=0){
     494                return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
     495        }
     496        else{
     497                Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters);
     498        }
     499
     500        /* Start  looping on the number of gaussian points: */
     501        gauss=new GaussTria(index1,index2,2);
     502        for(ig=gauss->begin();ig<gauss->end();ig++){
     503
     504                gauss->GaussPoint(ig);
     505
     506                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
     507
     508                vxaverage_input->GetInputValue(&vx,gauss);
     509                vyaverage_input->GetInputValue(&vy,gauss);
     510                UdotN=vx*normal[0]+vy*normal[1];
     511                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
     512                DL=gauss->weight*Jdet*UdotN;
     513
     514                TripleMultiply(&L[0],1,numdof,1,
     515                                        &DL,1,1,0,
     516                                        &L[0],1,numdof,0,
     517                                        &Ke_g[0][0],0);
     518
     519                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
     520        }
     521
     522        /*Clean up and return*/
     523        delete gauss;
     524        return Ke;
     525}
     526/*}}}*/
     527ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessInternal(void){/*{{{*/
     528
     529        /* constants*/
     530        const int numdof=NDOF1*NUMNODES_INTERNAL;
     531
     532        /* Intermediaries*/
     533        int        i,j,ig,index1,index2;
     534        IssmDouble     DL1,DL2,Jdet,vx,vy,UdotN;
     535        IssmDouble     xyz_list[NUMVERTICES][3];
     536        IssmDouble     normal[2];
     537        IssmDouble     B[numdof];
     538        IssmDouble     Bprime[numdof];
     539        IssmDouble     Ke_g1[numdof][numdof];
     540        IssmDouble     Ke_g2[numdof][numdof];
     541        GaussTria *gauss;
     542
     543        /*Initialize Element matrix and return if necessary*/
     544        Tria*  tria=(Tria*)element;
     545        if(!tria->IsIceInElement()) return NULL;
     546        ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES_INTERNAL,this->parameters);
     547
     548        /*Retrieve all inputs and parameters*/
     549        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     550        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
     551        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
     552        GetNormal(&normal[0],xyz_list);
     553
     554        /* Start  looping on the number of gaussian points: */
     555        index1=tria->GetNodeIndex(nodes[0]);
     556        index2=tria->GetNodeIndex(nodes[1]);
     557        gauss=new GaussTria(index1,index2,2);
     558        for(ig=gauss->begin();ig<gauss->end();ig++){
     559
     560                gauss->GaussPoint(ig);
     561
     562                tria->GetSegmentBFlux(&B[0],gauss,index1,index2,tria->FiniteElement());
     563                tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2,tria->FiniteElement());
     564
     565                vxaverage_input->GetInputValue(&vx,gauss);
     566                vyaverage_input->GetInputValue(&vy,gauss);
     567                UdotN=vx*normal[0]+vy*normal[1];
     568                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
     569                DL1=gauss->weight*Jdet*UdotN/2;
     570                DL2=gauss->weight*Jdet*fabs(UdotN)/2;
     571
     572                TripleMultiply(&B[0],1,numdof,1,
     573                                        &DL1,1,1,0,
     574                                        &Bprime[0],1,numdof,0,
     575                                        &Ke_g1[0][0],0);
     576                TripleMultiply(&B[0],1,numdof,1,
     577                                        &DL2,1,1,0,
     578                                        &B[0],1,numdof,0,
     579                                        &Ke_g2[0][0],0);
     580
     581                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g1[i][j];
     582                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g2[i][j];
     583        }
     584
     585        /*Clean up and return*/
     586        delete gauss;
     587        return Ke;
     588}
     589/*}}}*/
    419590ElementMatrix* Numericalflux::CreateKMatrixMasstransport(void){/*{{{*/
    420591
    … …  
    427598                        _error_("type not supported yet");
    428599        }
     600}
     601/*}}}*/
     602ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/
     603
     604        /* constants*/
     605        const int numdof=NDOF1*NUMNODES_BOUNDARY;
     606
     607        /* Intermediaries*/
     608        int        i,j,ig,index1,index2;
     609        IssmDouble     DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN;
     610        IssmDouble     xyz_list[NUMVERTICES][3];
     611        IssmDouble     normal[2];
     612        IssmDouble     L[numdof];
     613        IssmDouble     Ke_g[numdof][numdof];
     614        GaussTria *gauss;
     615
     616        /*Initialize Element matrix and return if necessary*/
     617        ElementMatrix* Ke = NULL;
     618        Tria*  tria=(Tria*)element;
     619        if(!tria->IsIceInElement()) return NULL;
     620
     621        /*Retrieve all inputs and parameters*/
     622        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     623        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     624        Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
     625        Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
     626        GetNormal(&normal[0],xyz_list);
     627
     628        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
     629        index1=tria->GetNodeIndex(nodes[0]);
     630        index2=tria->GetNodeIndex(nodes[1]);
     631
     632        gauss=new GaussTria();
     633        gauss->GaussEdgeCenter(index1,index2);
     634        vxaverage_input->GetInputValue(&mean_vx,gauss);
     635        vyaverage_input->GetInputValue(&mean_vy,gauss);
     636        delete gauss;
     637
     638        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
     639        if (UdotN<=0){
     640                return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
     641        }
     642        else{
     643                Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters);
     644        }
     645
     646        /* Start  looping on the number of gaussian points: */
     647        gauss=new GaussTria(index1,index2,2);
     648        for(ig=gauss->begin();ig<gauss->end();ig++){
     649
     650                gauss->GaussPoint(ig);
     651
     652                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
     653
     654                vxaverage_input->GetInputValue(&vx,gauss);
     655                vyaverage_input->GetInputValue(&vy,gauss);
     656                UdotN=vx*normal[0]+vy*normal[1];
     657                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
     658                DL=gauss->weight*Jdet*dt*UdotN;
     659
     660                TripleMultiply(&L[0],1,numdof,1,
     661                                        &DL,1,1,0,
     662                                        &L[0],1,numdof,0,
     663                                        &Ke_g[0][0],0);
     664
     665                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
     666        }
     667
     668        /*Clean up and return*/
     669        delete gauss;
     670        return Ke;
    429671}
    430672/*}}}*/
    … …  
    493735}
    494736/*}}}*/
    495 ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/
    496 
    497         /* constants*/
    498         const int numdof=NDOF1*NUMNODES_BOUNDARY;
    499 
    500         /* Intermediaries*/
    501         int        i,j,ig,index1,index2;
    502         IssmDouble     DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN;
    503         IssmDouble     xyz_list[NUMVERTICES][3];
    504         IssmDouble     normal[2];
    505         IssmDouble     L[numdof];
    506         IssmDouble     Ke_g[numdof][numdof];
    507         GaussTria *gauss;
    508 
    509         /*Initialize Element matrix and return if necessary*/
    510         ElementMatrix* Ke = NULL;
    511         Tria*  tria=(Tria*)element;
    512         if(!tria->IsIceInElement()) return NULL;
    513 
    514         /*Retrieve all inputs and parameters*/
    515         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    516         parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    517         Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
    518         Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
    519         GetNormal(&normal[0],xyz_list);
    520 
    521         /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    522         index1=tria->GetNodeIndex(nodes[0]);
    523         index2=tria->GetNodeIndex(nodes[1]);
    524 
    525         gauss=new GaussTria();
    526         gauss->GaussEdgeCenter(index1,index2);
    527         vxaverage_input->GetInputValue(&mean_vx,gauss);
    528         vyaverage_input->GetInputValue(&mean_vy,gauss);
    529         delete gauss;
    530 
    531         UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    532         if (UdotN<=0){
    533                 return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
    534         }
    535         else{
    536                 Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters);
    537         }
    538 
    539         /* Start  looping on the number of gaussian points: */
    540         gauss=new GaussTria(index1,index2,2);
    541         for(ig=gauss->begin();ig<gauss->end();ig++){
    542 
    543                 gauss->GaussPoint(ig);
    544 
    545                 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
    546 
    547                 vxaverage_input->GetInputValue(&vx,gauss);
    548                 vyaverage_input->GetInputValue(&vy,gauss);
    549                 UdotN=vx*normal[0]+vy*normal[1];
    550                 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
    551                 DL=gauss->weight*Jdet*dt*UdotN;
    552 
    553                 TripleMultiply(&L[0],1,numdof,1,
    554                                         &DL,1,1,0,
    555                                         &L[0],1,numdof,0,
    556                                         &Ke_g[0][0],0);
    557 
    558                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
    559         }
    560 
    561         /*Clean up and return*/
    562         delete gauss;
    563         return Ke;
    564 }
    565 /*}}}*/
    566 ElementMatrix* Numericalflux::CreateKMatrixBalancethickness(void){/*{{{*/
    567 
    568         switch(this->flux_type){
    569                 case InternalEnum:
    570                         return CreateKMatrixBalancethicknessInternal();
    571                 case BoundaryEnum:
    572                         return CreateKMatrixBalancethicknessBoundary();
    573                 default:
    574                         _error_("type not supported yet");
    575         }
    576 }
    577 /*}}}*/
    578 ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessInternal(void){/*{{{*/
    579 
    580         /* constants*/
    581         const int numdof=NDOF1*NUMNODES_INTERNAL;
    582 
    583         /* Intermediaries*/
    584         int        i,j,ig,index1,index2;
    585         IssmDouble     DL1,DL2,Jdet,vx,vy,UdotN;
    586         IssmDouble     xyz_list[NUMVERTICES][3];
    587         IssmDouble     normal[2];
    588         IssmDouble     B[numdof];
    589         IssmDouble     Bprime[numdof];
    590         IssmDouble     Ke_g1[numdof][numdof];
    591         IssmDouble     Ke_g2[numdof][numdof];
    592         GaussTria *gauss;
    593 
    594         /*Initialize Element matrix and return if necessary*/
    595         Tria*  tria=(Tria*)element;
    596         if(!tria->IsIceInElement()) return NULL;
    597         ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES_INTERNAL,this->parameters);
    598 
    599         /*Retrieve all inputs and parameters*/
    600         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    601         Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
    602         Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
    603         GetNormal(&normal[0],xyz_list);
    604 
    605         /* Start  looping on the number of gaussian points: */
    606         index1=tria->GetNodeIndex(nodes[0]);
    607         index2=tria->GetNodeIndex(nodes[1]);
    608         gauss=new GaussTria(index1,index2,2);
    609         for(ig=gauss->begin();ig<gauss->end();ig++){
    610 
    611                 gauss->GaussPoint(ig);
    612 
    613                 tria->GetSegmentBFlux(&B[0],gauss,index1,index2,tria->FiniteElement());
    614                 tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2,tria->FiniteElement());
    615 
    616                 vxaverage_input->GetInputValue(&vx,gauss);
    617                 vyaverage_input->GetInputValue(&vy,gauss);
    618                 UdotN=vx*normal[0]+vy*normal[1];
    619                 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
    620                 DL1=gauss->weight*Jdet*UdotN/2;
    621                 DL2=gauss->weight*Jdet*fabs(UdotN)/2;
    622 
    623                 TripleMultiply(&B[0],1,numdof,1,
    624                                         &DL1,1,1,0,
    625                                         &Bprime[0],1,numdof,0,
    626                                         &Ke_g1[0][0],0);
    627                 TripleMultiply(&B[0],1,numdof,1,
    628                                         &DL2,1,1,0,
    629                                         &B[0],1,numdof,0,
    630                                         &Ke_g2[0][0],0);
    631 
    632                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g1[i][j];
    633                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g2[i][j];
    634         }
    635 
    636         /*Clean up and return*/
    637         delete gauss;
    638         return Ke;
    639 }
    640 /*}}}*/
    641 ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessBoundary(void){/*{{{*/
    642 
    643         /* constants*/
    644         const int numdof=NDOF1*NUMNODES_BOUNDARY;
    645 
    646         /* Intermediaries*/
    647         int        i,j,ig,index1,index2;
    648         IssmDouble     DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN;
    649         IssmDouble     xyz_list[NUMVERTICES][3];
    650         IssmDouble     normal[2];
    651         IssmDouble     L[numdof];
    652         IssmDouble     Ke_g[numdof][numdof];
    653         GaussTria *gauss;
    654 
    655         /*Initialize Element matrix and return if necessary*/
    656         ElementMatrix* Ke = NULL;
    657         Tria*  tria=(Tria*)element;
    658         if(!tria->IsIceInElement()) return NULL;
    659 
    660         /*Retrieve all inputs and parameters*/
    661         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    662         Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
    663         Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
    664         GetNormal(&normal[0],xyz_list);
    665 
    666         /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    667         index1=tria->GetNodeIndex(nodes[0]);
    668         index2=tria->GetNodeIndex(nodes[1]);
    669 
    670         gauss=new GaussTria();
    671         gauss->GaussEdgeCenter(index1,index2);
    672         vxaverage_input->GetInputValue(&mean_vx,gauss);
    673         vyaverage_input->GetInputValue(&mean_vy,gauss);
    674         delete gauss;
    675 
    676         UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    677         if (UdotN<=0){
    678                 return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
    679         }
    680         else{
    681                 Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters);
    682         }
    683 
    684         /* Start  looping on the number of gaussian points: */
    685         gauss=new GaussTria(index1,index2,2);
    686         for(ig=gauss->begin();ig<gauss->end();ig++){
    687 
    688                 gauss->GaussPoint(ig);
    689 
    690                 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
    691 
    692                 vxaverage_input->GetInputValue(&vx,gauss);
    693                 vyaverage_input->GetInputValue(&vy,gauss);
    694                 UdotN=vx*normal[0]+vy*normal[1];
    695                 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
    696                 DL=gauss->weight*Jdet*UdotN;
    697 
    698                 TripleMultiply(&L[0],1,numdof,1,
    699                                         &DL,1,1,0,
    700                                         &L[0],1,numdof,0,
    701                                         &Ke_g[0][0],0);
    702 
    703                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
    704         }
    705 
    706         /*Clean up and return*/
    707         delete gauss;
    708         return Ke;
    709 }
    710 /*}}}*/
    711 ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethickness(void){/*{{{*/
    712 
    713         switch(this->flux_type){
    714                 case InternalEnum:
    715                         return CreateKMatrixAdjointBalancethicknessInternal();
    716                 case BoundaryEnum:
    717                         return CreateKMatrixAdjointBalancethicknessBoundary();
    718                 default:
    719                         _error_("type not supported yet");
    720         }
    721 }
    722 /*}}}*/
    723 ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessInternal(void){/*{{{*/
    724 
    725         ElementMatrix* Ke=CreateKMatrixBalancethicknessInternal();
    726         if (Ke) Ke->Transpose();
    727         return Ke;
    728 }
    729 /*}}}*/
    730 ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessBoundary(void){/*{{{*/
    731 
    732         ElementMatrix* Ke=CreateKMatrixBalancethicknessBoundary();
    733         if(Ke) Ke->Transpose();
    734         return Ke;
    735 }
    736 /*}}}*/
    737 ElementVector* Numericalflux::CreatePVectorMasstransport(void){/*{{{*/
    738 
    739         switch(this->flux_type){
    740                 case InternalEnum:
    741                         return CreatePVectorMasstransportInternal();
    742                 case BoundaryEnum:
    743                         return CreatePVectorMasstransportBoundary();
    744                 default:
    745                         _error_("type not supported yet");
    746         }
    747 }
    748 /*}}}*/
    749 ElementVector* Numericalflux::CreatePVectorMasstransportInternal(void){/*{{{*/
    750 
    751         /*Nothing added to PVector*/
     737ElementVector* Numericalflux::CreatePVectorAdjointBalancethickness(void){/*{{{*/
     738
     739        /*No PVector for the Adjoint*/
    752740        return NULL;
    753 
    754 }
    755 /*}}}*/
    756 ElementVector* Numericalflux::CreatePVectorMasstransportBoundary(void){/*{{{*/
    757 
    758         /* constants*/
    759         const int numdof=NDOF1*NUMNODES_BOUNDARY;
    760 
    761         /* Intermediaries*/
    762         int        i,ig,index1,index2;
    763         IssmDouble     DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN,thickness;
    764         IssmDouble     xyz_list[NUMVERTICES][3];
    765         IssmDouble     normal[2];
    766         IssmDouble     L[numdof];
    767         GaussTria *gauss;
    768 
    769         /*Initialize Load Vector and return if necessary*/
    770         ElementVector* pe = NULL;
    771         Tria*  tria=(Tria*)element;
    772         if(!tria->IsIceInElement()) return NULL;
    773 
    774         /*Retrieve all inputs and parameters*/
    775         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    776         parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    777         Input* vxaverage_input   =tria->inputs->GetInput(VxEnum);                     _assert_(vxaverage_input);
    778         Input* vyaverage_input   =tria->inputs->GetInput(VyEnum);                     _assert_(vyaverage_input);
    779         Input* spcthickness_input=tria->inputs->GetInput(MasstransportSpcthicknessEnum); _assert_(spcthickness_input);
    780         GetNormal(&normal[0],xyz_list);
    781 
    782         /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    783         index1=tria->GetNodeIndex(nodes[0]);
    784         index2=tria->GetNodeIndex(nodes[1]);
    785 
    786         gauss=new GaussTria();
    787         gauss->GaussEdgeCenter(index1,index2);
    788         vxaverage_input->GetInputValue(&mean_vx,gauss);
    789         vyaverage_input->GetInputValue(&mean_vy,gauss);
    790         delete gauss;
    791 
    792         UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    793         if (UdotN>0){
    794                 return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
    795         }
    796         else{
    797                 pe=new ElementVector(nodes,NUMNODES_BOUNDARY,this->parameters);
    798         }
    799 
    800         /* Start  looping on the number of gaussian points: */
    801         gauss=new GaussTria(index1,index2,2);
    802         for(ig=gauss->begin();ig<gauss->end();ig++){
    803 
    804                 gauss->GaussPoint(ig);
    805 
    806                 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
    807 
    808                 vxaverage_input->GetInputValue(&vx,gauss);
    809                 vyaverage_input->GetInputValue(&vy,gauss);
    810                 spcthickness_input->GetInputValue(&thickness,gauss);
    811                 if(xIsNan<IssmDouble>(thickness)) _error_("Cannot weakly apply constraint because NaN was provided");
    812 
    813                 UdotN=vx*normal[0]+vy*normal[1];
    814                 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
    815                 DL= - gauss->weight*Jdet*dt*UdotN*thickness;
    816 
    817                 for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
    818         }
    819 
    820         /*Clean up and return*/
    821         delete gauss;
    822         return pe;
    823741}
    824742/*}}}*/
    … …  
    833751                        _error_("type not supported yet");
    834752        }
    835 }
    836 /*}}}*/
    837 ElementVector* Numericalflux::CreatePVectorBalancethicknessInternal(void){/*{{{*/
    838 
    839         /*Nothing added to PVector*/
    840         return NULL;
    841 
    842753}
    843754/*}}}*/
    … …  
    908819}
    909820/*}}}*/
    910 ElementVector* Numericalflux::CreatePVectorAdjointBalancethickness(void){/*{{{*/
    911 
    912         /*No PVector for the Adjoint*/
     821ElementVector* Numericalflux::CreatePVectorBalancethicknessInternal(void){/*{{{*/
     822
     823        /*Nothing added to PVector*/
    913824        return NULL;
    914 }
    915 /*}}}*/
    916 void Numericalflux:: GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){/*{{{*/
     825
     826}
     827/*}}}*/
     828ElementVector* Numericalflux::CreatePVectorMasstransport(void){/*{{{*/
     829
     830        switch(this->flux_type){
     831                case InternalEnum:
     832                        return CreatePVectorMasstransportInternal();
     833                case BoundaryEnum:
     834                        return CreatePVectorMasstransportBoundary();
     835                default:
     836                        _error_("type not supported yet");
     837        }
     838}
     839/*}}}*/
     840ElementVector* Numericalflux::CreatePVectorMasstransportBoundary(void){/*{{{*/
     841
     842        /* constants*/
     843        const int numdof=NDOF1*NUMNODES_BOUNDARY;
     844
     845        /* Intermediaries*/
     846        int        i,ig,index1,index2;
     847        IssmDouble     DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN,thickness;
     848        IssmDouble     xyz_list[NUMVERTICES][3];
     849        IssmDouble     normal[2];
     850        IssmDouble     L[numdof];
     851        GaussTria *gauss;
     852
     853        /*Initialize Load Vector and return if necessary*/
     854        ElementVector* pe = NULL;
     855        Tria*  tria=(Tria*)element;
     856        if(!tria->IsIceInElement()) return NULL;
     857
     858        /*Retrieve all inputs and parameters*/
     859        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     860        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     861        Input* vxaverage_input   =tria->inputs->GetInput(VxEnum);                     _assert_(vxaverage_input);
     862        Input* vyaverage_input   =tria->inputs->GetInput(VyEnum);                     _assert_(vyaverage_input);
     863        Input* spcthickness_input=tria->inputs->GetInput(MasstransportSpcthicknessEnum); _assert_(spcthickness_input);
     864        GetNormal(&normal[0],xyz_list);
     865
     866        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
     867        index1=tria->GetNodeIndex(nodes[0]);
     868        index2=tria->GetNodeIndex(nodes[1]);
     869
     870        gauss=new GaussTria();
     871        gauss->GaussEdgeCenter(index1,index2);
     872        vxaverage_input->GetInputValue(&mean_vx,gauss);
     873        vyaverage_input->GetInputValue(&mean_vy,gauss);
     874        delete gauss;
     875
     876        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
     877        if (UdotN>0){
     878                return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
     879        }
     880        else{
     881                pe=new ElementVector(nodes,NUMNODES_BOUNDARY,this->parameters);
     882        }
     883
     884        /* Start  looping on the number of gaussian points: */
     885        gauss=new GaussTria(index1,index2,2);
     886        for(ig=gauss->begin();ig<gauss->end();ig++){
     887
     888                gauss->GaussPoint(ig);
     889
     890                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
     891
     892                vxaverage_input->GetInputValue(&vx,gauss);
     893                vyaverage_input->GetInputValue(&vy,gauss);
     894                spcthickness_input->GetInputValue(&thickness,gauss);
     895                if(xIsNan<IssmDouble>(thickness)) _error_("Cannot weakly apply constraint because NaN was provided");
     896
     897                UdotN=vx*normal[0]+vy*normal[1];
     898                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
     899                DL= - gauss->weight*Jdet*dt*UdotN*thickness;
     900
     901                for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
     902        }
     903
     904        /*Clean up and return*/
     905        delete gauss;
     906        return pe;
     907}
     908/*}}}*/
     909ElementVector* Numericalflux::CreatePVectorMasstransportInternal(void){/*{{{*/
     910
     911        /*Nothing added to PVector*/
     912        return NULL;
     913
     914}
     915/*}}}*/
     916void           Numericalflux::GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){/*{{{*/
    917917
    918918        /*Build unit outward pointing vector*/

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

    r18237 r18926  
    4040                /*}}}*/
    4141                /*Object virtual functions definitions:{{{ */
     42                Object *copy();
     43                void    DeepEcho();
    4244                void    Echo();
    43                 void    DeepEcho();
    4445                int     Id();
    4546                int     ObjectEnum();
    46                 Object *copy();
    4747                /*}}}*/
    4848                /*Update virtual functions resolution: {{{*/
    49                 void InputUpdateFromVector(IssmDouble* vector, int name, int type){/*Do nothing*/}
    50                 void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*Do nothing*/}
    51                 void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*Do nothing*/}
    5249                void InputUpdateFromConstant(IssmDouble constant, int name){/*Do nothing*/};
    5350                void InputUpdateFromConstant(int constant, int name){/*Do nothing*/};
    5451                void InputUpdateFromConstant(bool constant, int name){_error_("Not implemented yet!");}
    5552                void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
     53                void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*Do nothing*/}
     54                void InputUpdateFromVector(IssmDouble* vector, int name, int type){/*Do nothing*/}
     55                void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*Do nothing*/}
    5656                /*}}}*/
    5757                /*Load virtual functions definitions: {{{*/
    5858                void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    59                 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    60                 void ResetHooks();
     59                void CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("Not implemented yet");};
    6160                void CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
    6261                void CreatePVector(Vector<IssmDouble>* pf);
     62                void GetNodesLidList(int* lidlist);
    6363                void GetNodesSidList(int* sidlist);
    64                 void GetNodesLidList(int* lidlist);
    6564                int  GetNumberOfNodes(void);
    66                 void CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("Not implemented yet");};
     65                bool InAnalysis(int analysis_type);
    6766                bool IsPenalty(void);
    6867                void PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){_error_("Not implemented yet");};
    6968                void PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* kfs, IssmDouble kmax);
    7069                void PenaltyCreatePVector(Vector<IssmDouble>* pf, IssmDouble kmax);
     70                void ResetHooks();
    7171                void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
    72                 bool InAnalysis(int analysis_type);
     72                void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    7373                /*}}}*/
    7474                /*Numericalflux management:{{{*/
     75                ElementMatrix* CreateKMatrixAdjointBalancethickness(void);
     76                ElementMatrix* CreateKMatrixAdjointBalancethicknessBoundary(void);
     77                ElementMatrix* CreateKMatrixAdjointBalancethicknessInternal(void);
     78                ElementMatrix* CreateKMatrixBalancethickness(void);
     79                ElementMatrix* CreateKMatrixBalancethicknessBoundary(void);
     80                ElementMatrix* CreateKMatrixBalancethicknessInternal(void);
     81                ElementMatrix* CreateKMatrixMasstransport(void);
     82                ElementMatrix* CreateKMatrixMasstransportBoundary(void);
     83                ElementMatrix* CreateKMatrixMasstransportInternal(void);
     84                ElementVector* CreatePVectorAdjointBalancethickness(void);
     85                ElementVector* CreatePVectorBalancethickness(void);
     86                ElementVector* CreatePVectorBalancethicknessBoundary(void);
     87                ElementVector* CreatePVectorBalancethicknessInternal(void);
     88                ElementVector* CreatePVectorMasstransport(void);
     89                ElementVector* CreatePVectorMasstransportBoundary(void);
     90                ElementVector* CreatePVectorMasstransportInternal(void);
    7591                void           GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
    76                 ElementMatrix* CreateKMatrixMasstransport(void);
    77                 ElementMatrix* CreateKMatrixMasstransportInternal(void);
    78                 ElementMatrix* CreateKMatrixMasstransportBoundary(void);
    79                 ElementMatrix* CreateKMatrixBalancethickness(void);
    80                 ElementMatrix* CreateKMatrixBalancethicknessInternal(void);
    81                 ElementMatrix* CreateKMatrixBalancethicknessBoundary(void);
    82                 ElementMatrix* CreateKMatrixAdjointBalancethickness(void);
    83                 ElementMatrix* CreateKMatrixAdjointBalancethicknessInternal(void);
    84                 ElementMatrix* CreateKMatrixAdjointBalancethicknessBoundary(void);
    85                 ElementVector* CreatePVectorMasstransport(void);
    86                 ElementVector* CreatePVectorMasstransportInternal(void);
    87                 ElementVector* CreatePVectorMasstransportBoundary(void);
    88                 ElementVector* CreatePVectorBalancethickness(void);
    89                 ElementVector* CreatePVectorBalancethicknessInternal(void);
    90                 ElementVector* CreatePVectorBalancethicknessBoundary(void);
    91                 ElementVector* CreatePVectorAdjointBalancethickness(void);
    9292                /*}}}*/
    9393

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

    r18237 r18926  
    8181
    8282/*Object virtual functions definitions:*/
    83 void Pengrid::Echo(void){/*{{{*/
    84         this->DeepEcho();
    85 }
    86 /*}}}*/
    87 void Pengrid::DeepEcho(void){/*{{{*/
     83Object* Pengrid::copy() {/*{{{*/
     84
     85        Pengrid* pengrid=NULL;
     86
     87        pengrid=new Pengrid();
     88
     89        /*copy fields: */
     90        pengrid->id=this->id;
     91        pengrid->analysis_type=this->analysis_type;
     92
     93        /*point parameters: */
     94        pengrid->parameters=this->parameters;
     95
     96        /*now deal with hooks and objects: */
     97        pengrid->hnode=(Hook*)this->hnode->copy();
     98        pengrid->hmatpar=(Hook*)this->hmatpar->copy();
     99        pengrid->helement=(Hook*)this->helement->copy();
     100
     101        /*corresponding fields*/
     102        pengrid->node  =(Node*)pengrid->hnode->delivers();
     103        pengrid->matpar =(Matpar*)pengrid->hmatpar->delivers();
     104        pengrid->element=(Element*)pengrid->helement->delivers();
     105
     106        //let's not forget internals
     107        pengrid->active=this->active=0;
     108        pengrid->zigzag_counter=this->zigzag_counter=0;
     109
     110        return pengrid;
     111
     112}
     113/*}}}*/
     114void    Pengrid::DeepEcho(void){/*{{{*/
    88115
    89116        _printf_("Pengrid:\n");
    … …  
    99126}
    100127/*}}}*/
    101 int    Pengrid::Id(void){ return id; }/*{{{*/
    102 /*}}}*/
    103 int Pengrid::ObjectEnum(void){/*{{{*/
     128void    Pengrid::Echo(void){/*{{{*/
     129        this->DeepEcho();
     130}
     131/*}}}*/
     132int     Pengrid::Id(void){ return id; }/*{{{*/
     133/*}}}*/
     134int     Pengrid::ObjectEnum(void){/*{{{*/
    104135
    105136        return PengridEnum;
    106 }
    107 /*}}}*/
    108 Object* Pengrid::copy() {/*{{{*/
    109 
    110         Pengrid* pengrid=NULL;
    111 
    112         pengrid=new Pengrid();
    113 
    114         /*copy fields: */
    115         pengrid->id=this->id;
    116         pengrid->analysis_type=this->analysis_type;
    117 
    118         /*point parameters: */
    119         pengrid->parameters=this->parameters;
    120 
    121         /*now deal with hooks and objects: */
    122         pengrid->hnode=(Hook*)this->hnode->copy();
    123         pengrid->hmatpar=(Hook*)this->hmatpar->copy();
    124         pengrid->helement=(Hook*)this->helement->copy();
    125 
    126         /*corresponding fields*/
    127         pengrid->node  =(Node*)pengrid->hnode->delivers();
    128         pengrid->matpar =(Matpar*)pengrid->hmatpar->delivers();
    129         pengrid->element=(Element*)pengrid->helement->delivers();
    130 
    131         //let's not forget internals
    132         pengrid->active=this->active=0;
    133         pengrid->zigzag_counter=this->zigzag_counter=0;
    134 
    135         return pengrid;
    136 
    137137}
    138138/*}}}*/
    … …  
    154154        /*point parameters to real dataset: */
    155155        this->parameters=parametersin;
    156 }
    157 /*}}}*/
    158 void  Pengrid::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
    159 
    160 }
    161 /*}}}*/
    162 void  Pengrid::ResetHooks(){/*{{{*/
    163 
    164         this->node=NULL;
    165         this->element=NULL;
    166         this->matpar=NULL;
    167         this->parameters=NULL;
    168 
    169         /*Get Element type*/
    170         this->hnode->reset();
    171         this->helement->reset();
    172         this->hmatpar->reset();
    173 
    174156}
    175157/*}}}*/
    … …  
    203185}
    204186/*}}}*/
     187void  Pengrid::GetNodesLidList(int* lidlist){/*{{{*/
     188
     189        _assert_(lidlist);
     190        _assert_(node);
     191
     192        lidlist[0]=node->Lid();
     193}
     194/*}}}*/
    205195void  Pengrid::GetNodesSidList(int* sidlist){/*{{{*/
    206196
    … …  
    211201}
    212202/*}}}*/
    213 void  Pengrid::GetNodesLidList(int* lidlist){/*{{{*/
    214 
    215         _assert_(lidlist);
    216         _assert_(node);
    217 
    218         lidlist[0]=node->Lid();
    219 }
    220 /*}}}*/
    221203int   Pengrid::GetNumberOfNodes(void){/*{{{*/
    222204
    223205        return NUMVERTICES;
     206}
     207/*}}}*/
     208bool  Pengrid::InAnalysis(int in_analysis_type){/*{{{*/
     209        if (in_analysis_type==this->analysis_type)return true;
     210        else return false;
     211}
     212/*}}}*/
     213bool  Pengrid::IsPenalty(void){/*{{{*/
     214        return true;
    224215}
    225216/*}}}*/
    … …  
    282273}
    283274/*}}}*/
    284 bool  Pengrid::InAnalysis(int in_analysis_type){/*{{{*/
    285         if (in_analysis_type==this->analysis_type)return true;
    286         else return false;
    287 }
    288 /*}}}*/
    289 bool  Pengrid::IsPenalty(void){/*{{{*/
    290         return true;
     275void  Pengrid::ResetHooks(){/*{{{*/
     276
     277        this->node=NULL;
     278        this->element=NULL;
     279        this->matpar=NULL;
     280        this->parameters=NULL;
     281
     282        /*Get Element type*/
     283        this->hnode->reset();
     284        this->helement->reset();
     285        this->hmatpar->reset();
     286
     287}
     288/*}}}*/
     289void  Pengrid::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
     290
    291291}
    292292/*}}}*/
    … …  
    343343
    344344/*Update virtual functions definitions:*/
    345 void  Pengrid::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/
    346         /*Nothing updated yet*/
    347 }
    348 /*}}}*/
    349 void  Pengrid::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
    350         /*Nothing updated yet*/
    351 }
    352 /*}}}*/
    353 void  Pengrid::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/
    354         /*Nothing updated yet*/
    355 }
    356 /*}}}*/
    357345void  Pengrid::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/
    358346        /*Nothing*/
    … …  
    374362}
    375363/*}}}*/
     364void  Pengrid::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
     365        /*Nothing updated yet*/
     366}
     367/*}}}*/
     368void  Pengrid::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/
     369        /*Nothing updated yet*/
     370}
     371/*}}}*/
     372void  Pengrid::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/
     373        /*Nothing updated yet*/
     374}
     375/*}}}*/
    376376
    377377/*Pengrid management:*/
    378 void  Pengrid::ConstraintActivate(int* punstable){/*{{{*/
     378void           Pengrid::ConstraintActivate(int* punstable){/*{{{*/
    379379
    380380        int analysis_type;
    … …  
    404404}
    405405/*}}}*/
    406 void  Pengrid::ConstraintActivateThermal(int* punstable){/*{{{*/
    407 
    408         //   The penalty is stable if it doesn't change during to successive iterations.   
    409         IssmDouble pressure;
    410         IssmDouble temperature;
    411         IssmDouble t_pmp;
    412         int        new_active;
    413         int        unstable=0;
    414         int        penalty_lock;
    415 
    416         /*recover pointers: */
    417         Penta* penta=(Penta*)element;
    418 
    419         /*check that pengrid is not a clone (penalty to be added only once)*/
    420         if (node->IsClone()){
    421                 unstable=0;
    422                 *punstable=unstable;
    423                 return;
    424         }
    425 
    426         //First recover pressure and temperature values, using the element: */
    427         penta->GetInputValue(&pressure,node,PressureEnum);
    428         penta->GetInputValue(&temperature,node,TemperaturePicardEnum);
    429 
    430         //Recover our data:
    431         parameters->FindParam(&penalty_lock,ThermalPenaltyLockEnum);
    432 
    433         //Compute pressure melting point
    434         t_pmp=matpar->TMeltingPoint(pressure);
    435 
    436         //Figure out if temperature is over melting_point, in which case, this penalty needs to be activated.
    437 
    438         if (temperature>t_pmp){
    439                 new_active=1;
    440         }
    441         else{
    442                 new_active=0;
    443         }
    444 
    445         //Figure out stability of this penalty
    446         if (active==new_active){
    447                 unstable=0;
    448         }
    449         else{
    450                 unstable=1;
    451                 if(penalty_lock)zigzag_counter++;
    452         }
    453 
    454         /*If penalty keeps zigzagging more than 5 times: */
    455         if(penalty_lock){
    456                 if(zigzag_counter>penalty_lock){
    457                         unstable=0;
    458                         active=1;
    459                 }
    460         }
    461 
    462         //Set penalty flag
    463         active=new_active;
    464 
    465         //*Assign output pointers:*/
    466         *punstable=unstable;
    467 }
    468 /*}}}*/
    469 ElementMatrix* Pengrid::PenaltyCreateKMatrixMelting(IssmDouble kmax){/*{{{*/
    470 
    471         IssmDouble pressure,temperature,t_pmp;
    472         IssmDouble penalty_factor;
    473 
    474         Penta* penta=(Penta*)element;
    475 
    476         /*check that pengrid is not a clone (penalty to be added only once)*/
    477         if (node->IsClone()) return NULL;
    478         ElementMatrix* Ke=new ElementMatrix(&node,1,this->parameters);
    479 
    480         /*Retrieve all parameters*/
    481         penta->GetInputValue(&pressure,node,PressureEnum);
    482         penta->GetInputValue(&temperature,node,TemperatureEnum);
    483         parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
    484 
    485         /*Compute pressure melting point*/
    486         t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
    487 
    488         /*Add penalty load*/
    489         if (temperature<t_pmp){ //If T<Tpmp, there must be no melting. Therefore, melting should be  constrained to 0 when T<Tpmp, instead of using spcs, use penalties
    490                 Ke->values[0]=kmax*pow(10.,penalty_factor);
    491         }
    492 
    493         /*Clean up and return*/
    494         return Ke;
    495 }
    496 /*}}}*/
    497 ElementMatrix* Pengrid::PenaltyCreateKMatrixThermal(IssmDouble kmax){/*{{{*/
    498 
    499         IssmDouble    penalty_factor;
    500 
    501         /*Initialize Element matrix and return if necessary*/
    502         if(!this->active) return NULL;
    503         ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters);
    504 
    505         /*recover parameters: */
    506         parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
    507 
    508         Ke->values[0]=kmax*pow(10.,penalty_factor);
    509 
    510         /*Clean up and return*/
    511         return Ke;
    512 }
    513 /*}}}*/
    514 ElementVector* Pengrid::PenaltyCreatePVectorMelting(IssmDouble kmax){/*{{{*/
    515 
    516         IssmDouble pressure;
    517         IssmDouble temperature;
    518         IssmDouble melting_offset;
    519         IssmDouble t_pmp;
    520         IssmDouble dt,penalty_factor;
    521 
    522         /*recover pointers: */
    523         Penta* penta=(Penta*)element;
    524 
    525         /*check that pengrid is not a clone (penalty to be added only once)*/
    526         if (node->IsClone()) return NULL;
    527         ElementVector* pe=new ElementVector(&node,NUMVERTICES,this->parameters);
    528 
    529         /*Retrieve all parameters*/
    530         penta->GetInputValue(&pressure,node,PressureEnum);
    531         penta->GetInputValue(&temperature,node,TemperatureEnum);
    532         parameters->FindParam(&melting_offset,MeltingOffsetEnum);
    533         parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    534         parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
    535 
    536         /*Compute pressure melting point*/
    537         t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
    538 
    539         /*Add penalty load
    540           This time, the penalty must have the same value as the one used for the thermal computation
    541           so that the corresponding melting can be computed correctly
    542           In the thermal computation, we used kmax=melting_offset, and the same penalty_factor*/
    543         if (temperature<t_pmp){ //%no melting
    544                 pe->values[0]=0;
    545         }
    546         else{
    547                 if (reCast<bool>(dt)) pe->values[0]=melting_offset*pow(10.,penalty_factor)*(temperature-t_pmp)/dt;
    548                 else    pe->values[0]=melting_offset*pow(10.,penalty_factor)*(temperature-t_pmp);
    549         }
    550 
    551         /*Clean up and return*/
    552         return pe;
    553 }
    554 /*}}}*/
    555 ElementVector* Pengrid::PenaltyCreatePVectorThermal(IssmDouble kmax){/*{{{*/
    556 
    557         IssmDouble pressure;
    558         IssmDouble t_pmp;
    559         IssmDouble penalty_factor;
    560 
    561         Penta* penta=(Penta*)element;
    562 
    563         /*Initialize Element matrix and return if necessary*/
    564         if(!this->active) return NULL;
    565         ElementVector* pe=new ElementVector(&node,1,this->parameters);
    566 
    567         /*Retrieve all parameters*/
    568         penta->GetInputValue(&pressure,node,PressureEnum);
    569         parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
    570 
    571         /*Compute pressure melting point*/
    572         t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
    573 
    574         pe->values[0]=kmax*pow(10.,penalty_factor)*t_pmp;
    575 
    576         /*Clean up and return*/
    577         return pe;
    578 }
    579 /*}}}*/
    580 void  Pengrid::ConstraintActivateHydrologyDCInefficient(int* punstable){/*{{{*/
     406void           Pengrid::ConstraintActivateHydrologyDCInefficient(int* punstable){/*{{{*/
    581407
    582408        //   The penalty is stable if it doesn't change during two consecutive iterations.   
    … …  
    638464}
    639465/*}}}*/
     466void           Pengrid::ConstraintActivateThermal(int* punstable){/*{{{*/
     467
     468        //   The penalty is stable if it doesn't change during to successive iterations.   
     469        IssmDouble pressure;
     470        IssmDouble temperature;
     471        IssmDouble t_pmp;
     472        int        new_active;
     473        int        unstable=0;
     474        int        penalty_lock;
     475
     476        /*recover pointers: */
     477        Penta* penta=(Penta*)element;
     478
     479        /*check that pengrid is not a clone (penalty to be added only once)*/
     480        if (node->IsClone()){
     481                unstable=0;
     482                *punstable=unstable;
     483                return;
     484        }
     485
     486        //First recover pressure and temperature values, using the element: */
     487        penta->GetInputValue(&pressure,node,PressureEnum);
     488        penta->GetInputValue(&temperature,node,TemperaturePicardEnum);
     489
     490        //Recover our data:
     491        parameters->FindParam(&penalty_lock,ThermalPenaltyLockEnum);
     492
     493        //Compute pressure melting point
     494        t_pmp=matpar->TMeltingPoint(pressure);
     495
     496        //Figure out if temperature is over melting_point, in which case, this penalty needs to be activated.
     497
     498        if (temperature>t_pmp){
     499                new_active=1;
     500        }
     501        else{
     502                new_active=0;
     503        }
     504
     505        //Figure out stability of this penalty
     506        if (active==new_active){
     507                unstable=0;
     508        }
     509        else{
     510                unstable=1;
     511                if(penalty_lock)zigzag_counter++;
     512        }
     513
     514        /*If penalty keeps zigzagging more than 5 times: */
     515        if(penalty_lock){
     516                if(zigzag_counter>penalty_lock){
     517                        unstable=0;
     518                        active=1;
     519                }
     520        }
     521
     522        //Set penalty flag
     523        active=new_active;
     524
     525        //*Assign output pointers:*/
     526        *punstable=unstable;
     527}
     528/*}}}*/
     529ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){/*{{{*/
     530
     531        IssmDouble moulin_load,dt;
     532
     533        /*Initialize Element matrix*/
     534        ElementVector* pe=new ElementVector(&node,1,this->parameters);
     535
     536        this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
     537        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     538
     539        if(dt!=0.0) pe->values[0]=moulin_load*dt;
     540
     541        /*Clean up and return*/
     542        return pe;
     543 }
     544/*}}}*/
    640545ElementMatrix* Pengrid::PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax){/*{{{*/
    641546        IssmDouble    penalty_factor;
    … …  
    647552        if(!this->active) return NULL;
    648553        ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters);
     554
     555        Ke->values[0]=kmax*pow(10.,penalty_factor);
     556
     557        /*Clean up and return*/
     558        return Ke;
     559}
     560/*}}}*/
     561ElementMatrix* Pengrid::PenaltyCreateKMatrixMelting(IssmDouble kmax){/*{{{*/
     562
     563        IssmDouble pressure,temperature,t_pmp;
     564        IssmDouble penalty_factor;
     565
     566        Penta* penta=(Penta*)element;
     567
     568        /*check that pengrid is not a clone (penalty to be added only once)*/
     569        if (node->IsClone()) return NULL;
     570        ElementMatrix* Ke=new ElementMatrix(&node,1,this->parameters);
     571
     572        /*Retrieve all parameters*/
     573        penta->GetInputValue(&pressure,node,PressureEnum);
     574        penta->GetInputValue(&temperature,node,TemperatureEnum);
     575        parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
     576
     577        /*Compute pressure melting point*/
     578        t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
     579
     580        /*Add penalty load*/
     581        if (temperature<t_pmp){ //If T<Tpmp, there must be no melting. Therefore, melting should be  constrained to 0 when T<Tpmp, instead of using spcs, use penalties
     582                Ke->values[0]=kmax*pow(10.,penalty_factor);
     583        }
     584
     585        /*Clean up and return*/
     586        return Ke;
     587}
     588/*}}}*/
     589ElementMatrix* Pengrid::PenaltyCreateKMatrixThermal(IssmDouble kmax){/*{{{*/
     590
     591        IssmDouble    penalty_factor;
     592
     593        /*Initialize Element matrix and return if necessary*/
     594        if(!this->active) return NULL;
     595        ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters);
     596
     597        /*recover parameters: */
     598        parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
    649599
    650600        Ke->values[0]=kmax*pow(10.,penalty_factor);
    … …  
    678628}
    679629/*}}}*/
    680 ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){/*{{{*/
    681 
    682         IssmDouble moulin_load,dt;
    683 
    684         /*Initialize Element matrix*/
    685         ElementVector* pe=new ElementVector(&node,1,this->parameters);
    686 
    687         this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
     630ElementVector* Pengrid::PenaltyCreatePVectorMelting(IssmDouble kmax){/*{{{*/
     631
     632        IssmDouble pressure;
     633        IssmDouble temperature;
     634        IssmDouble melting_offset;
     635        IssmDouble t_pmp;
     636        IssmDouble dt,penalty_factor;
     637
     638        /*recover pointers: */
     639        Penta* penta=(Penta*)element;
     640
     641        /*check that pengrid is not a clone (penalty to be added only once)*/
     642        if (node->IsClone()) return NULL;
     643        ElementVector* pe=new ElementVector(&node,NUMVERTICES,this->parameters);
     644
     645        /*Retrieve all parameters*/
     646        penta->GetInputValue(&pressure,node,PressureEnum);
     647        penta->GetInputValue(&temperature,node,TemperatureEnum);
     648        parameters->FindParam(&melting_offset,MeltingOffsetEnum);
    688649        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    689 
    690         if(dt!=0.0) pe->values[0]=moulin_load*dt;
     650        parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
     651
     652        /*Compute pressure melting point*/
     653        t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
     654
     655        /*Add penalty load
     656          This time, the penalty must have the same value as the one used for the thermal computation
     657          so that the corresponding melting can be computed correctly
     658          In the thermal computation, we used kmax=melting_offset, and the same penalty_factor*/
     659        if (temperature<t_pmp){ //%no melting
     660                pe->values[0]=0;
     661        }
     662        else{
     663                if (reCast<bool>(dt)) pe->values[0]=melting_offset*pow(10.,penalty_factor)*(temperature-t_pmp)/dt;
     664                else    pe->values[0]=melting_offset*pow(10.,penalty_factor)*(temperature-t_pmp);
     665        }
    691666
    692667        /*Clean up and return*/
    693668        return pe;
    694  }
    695 /*}}}*/
    696 void  Pengrid::ResetConstraint(void){/*{{{*/
     669}
     670/*}}}*/
     671ElementVector* Pengrid::PenaltyCreatePVectorThermal(IssmDouble kmax){/*{{{*/
     672
     673        IssmDouble pressure;
     674        IssmDouble t_pmp;
     675        IssmDouble penalty_factor;
     676
     677        Penta* penta=(Penta*)element;
     678
     679        /*Initialize Element matrix and return if necessary*/
     680        if(!this->active) return NULL;
     681        ElementVector* pe=new ElementVector(&node,1,this->parameters);
     682
     683        /*Retrieve all parameters*/
     684        penta->GetInputValue(&pressure,node,PressureEnum);
     685        parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
     686
     687        /*Compute pressure melting point*/
     688        t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
     689
     690        pe->values[0]=kmax*pow(10.,penalty_factor)*t_pmp;
     691
     692        /*Clean up and return*/
     693        return pe;
     694}
     695/*}}}*/
     696void           Pengrid::ResetConstraint(void){/*{{{*/
    697697        active         = 0;
    698698        zigzag_counter = 0;

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

    r18237 r18926  
    4545
    4646/*Object virtual functions definitions:*/
    47 void Penpair::Echo(void){/*{{{*/
     47Object* Penpair::copy() {/*{{{*/
     48
     49        Penpair* penpair=NULL;
     50
     51        penpair=new Penpair();
     52
     53        /*copy fields: */
     54        penpair->id=this->id;
     55        penpair->analysis_type=this->analysis_type;
     56
     57        /*now deal with hooks and objects: */
     58        penpair->hnodes=(Hook*)this->hnodes->copy();
     59        penpair->nodes =(Node**)penpair->hnodes->deliverp();
     60
     61        /*point parameters: */
     62        penpair->parameters=this->parameters;
     63
     64        return penpair;
     65
     66}
     67/*}}}*/
     68void    Penpair::DeepEcho(void){/*{{{*/
     69
     70        _printf_("Penpair:\n");
     71        _printf_("   id: " << id << "\n");
     72        _printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
     73        hnodes->DeepEcho();
     74
     75        return;
     76}               
     77/*}}}*/
     78void    Penpair::Echo(void){/*{{{*/
    4879
    4980        _printf_("Penpair:\n");
    … …  
    5586}
    5687/*}}}*/
    57 void Penpair::DeepEcho(void){/*{{{*/
    58 
    59         _printf_("Penpair:\n");
    60         _printf_("   id: " << id << "\n");
    61         _printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
    62         hnodes->DeepEcho();
    63 
    64         return;
    65 }               
    66 /*}}}*/
    67 int  Penpair::Id(void){ return id; }/*{{{*/
    68 /*}}}*/
    69 int  Penpair::ObjectEnum(void){/*{{{*/
     88int     Penpair::Id(void){ return id; }/*{{{*/
     89/*}}}*/
     90int     Penpair::ObjectEnum(void){/*{{{*/
    7091
    7192        return PenpairEnum;
    72 }
    73 /*}}}*/
    74 Object* Penpair::copy() {/*{{{*/
    75 
    76         Penpair* penpair=NULL;
    77 
    78         penpair=new Penpair();
    79 
    80         /*copy fields: */
    81         penpair->id=this->id;
    82         penpair->analysis_type=this->analysis_type;
    83 
    84         /*now deal with hooks and objects: */
    85         penpair->hnodes=(Hook*)this->hnodes->copy();
    86         penpair->nodes =(Node**)penpair->hnodes->deliverp();
    87 
    88         /*point parameters: */
    89         penpair->parameters=this->parameters;
    90 
    91         return penpair;
    92 
    9393}
    9494/*}}}*/
    … …  
    109109}
    110110/*}}}*/
    111 void  Penpair::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
    112 
    113 }
    114 /*}}}*/
    115 void  Penpair::ResetHooks(){/*{{{*/
    116 
    117         this->nodes=NULL;
    118         this->parameters=NULL;
    119 
    120         /*Get Element type*/
    121         this->hnodes->reset();
    122 
     111void  Penpair::CreateJacobianMatrix(Matrix<IssmDouble>* Jff){/*{{{*/
     112        this->CreateKMatrix(Jff,NULL);
    123113}
    124114/*}}}*/
    … …  
    137127}
    138128/*}}}*/
    139 void  Penpair::CreateJacobianMatrix(Matrix<IssmDouble>* Jff){/*{{{*/
    140         this->CreateKMatrix(Jff,NULL);
    141 }
    142 /*}}}*/
    143129void  Penpair::GetNodesSidList(int* sidlist){/*{{{*/
    144130
    … …  
    162148}
    163149/*}}}*/
     150bool  Penpair::InAnalysis(int in_analysis_type){/*{{{*/
     151        if (in_analysis_type==this->analysis_type)return true;
     152        else return false;
     153}
     154/*}}}*/
    164155bool  Penpair::IsPenalty(void){/*{{{*/
    165156        return true;
     157}
     158/*}}}*/
     159void  Penpair::PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){/*{{{*/
     160        this->PenaltyCreateKMatrix(Jff,NULL,kmax);
    166161}
    167162/*}}}*/
    … …  
    196191}
    197192/*}}}*/
    198 void  Penpair::PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){/*{{{*/
    199         this->PenaltyCreateKMatrix(Jff,NULL,kmax);
    200 }
    201 /*}}}*/
    202 bool  Penpair::InAnalysis(int in_analysis_type){/*{{{*/
    203         if (in_analysis_type==this->analysis_type)return true;
    204         else return false;
     193void  Penpair::ResetHooks(){/*{{{*/
     194
     195        this->nodes=NULL;
     196        this->parameters=NULL;
     197
     198        /*Get Element type*/
     199        this->hnodes->reset();
     200
     201}
     202/*}}}*/
     203void  Penpair::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
     204
    205205}
    206206/*}}}*/
    … …  
    279279
    280280/*Penpair management:*/
     281ElementMatrix* Penpair::PenaltyCreateKMatrixMasstransport(IssmDouble kmax){/*{{{*/
     282
     283        const int numdof=NUMVERTICES*NDOF1;
     284        IssmDouble penalty_factor;
     285
     286        /*Initialize Element vector and return if necessary*/
     287        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);
     288
     289        /*recover parameters: */
     290        parameters->FindParam(&penalty_factor,MasstransportPenaltyFactorEnum);
     291
     292        //Create elementary matrix: add penalty to
     293        Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_factor);
     294        Ke->values[0*numdof+1]=-kmax*pow(10.,penalty_factor);
     295        Ke->values[1*numdof+0]=-kmax*pow(10.,penalty_factor);
     296        Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_factor);
     297
     298        /*Clean up and return*/
     299        return Ke;
     300}
     301/*}}}*/
     302ElementMatrix* Penpair::PenaltyCreateKMatrixStressbalanceFS(IssmDouble kmax){/*{{{*/
     303
     304        const int  numdof=NUMVERTICES*NDOF3;
     305        IssmDouble penalty_offset;
     306
     307        /*Initialize Element vector and return if necessary*/
     308        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);
     309
     310        /*recover parameters: */
     311        parameters->FindParam(&penalty_offset,StressbalancePenaltyFactorEnum);
     312
     313        //Create elementary matrix: add penalty to
     314        Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_offset);
     315        Ke->values[0*numdof+3]=-kmax*pow(10.,penalty_offset);
     316        Ke->values[3*numdof+0]=-kmax*pow(10.,penalty_offset);
     317        Ke->values[3*numdof+3]=+kmax*pow(10.,penalty_offset);
     318
     319        Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_offset);
     320        Ke->values[1*numdof+4]=-kmax*pow(10.,penalty_offset);
     321        Ke->values[4*numdof+1]=-kmax*pow(10.,penalty_offset);
     322        Ke->values[4*numdof+4]=+kmax*pow(10.,penalty_offset);
     323
     324        Ke->values[2*numdof+2]=+kmax*pow(10.,penalty_offset);
     325        Ke->values[2*numdof+5]=-kmax*pow(10.,penalty_offset);
     326        Ke->values[5*numdof+2]=-kmax*pow(10.,penalty_offset);
     327        Ke->values[5*numdof+5]=+kmax*pow(10.,penalty_offset);
     328
     329        /*Clean up and return*/
     330        return Ke;
     331}
     332/*}}}*/
    281333ElementMatrix* Penpair::PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax){/*{{{*/
    282334
    … …  
    338390}
    339391/*}}}*/
    340 ElementMatrix* Penpair::PenaltyCreateKMatrixStressbalanceFS(IssmDouble kmax){/*{{{*/
    341 
    342         const int  numdof=NUMVERTICES*NDOF3;
    343         IssmDouble penalty_offset;
    344 
    345         /*Initialize Element vector and return if necessary*/
    346         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);
    347 
    348         /*recover parameters: */
    349         parameters->FindParam(&penalty_offset,StressbalancePenaltyFactorEnum);
    350 
    351         //Create elementary matrix: add penalty to
    352         Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_offset);
    353         Ke->values[0*numdof+3]=-kmax*pow(10.,penalty_offset);
    354         Ke->values[3*numdof+0]=-kmax*pow(10.,penalty_offset);
    355         Ke->values[3*numdof+3]=+kmax*pow(10.,penalty_offset);
    356 
    357         Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_offset);
    358         Ke->values[1*numdof+4]=-kmax*pow(10.,penalty_offset);
    359         Ke->values[4*numdof+1]=-kmax*pow(10.,penalty_offset);
    360         Ke->values[4*numdof+4]=+kmax*pow(10.,penalty_offset);
    361 
    362         Ke->values[2*numdof+2]=+kmax*pow(10.,penalty_offset);
    363         Ke->values[2*numdof+5]=-kmax*pow(10.,penalty_offset);
    364         Ke->values[5*numdof+2]=-kmax*pow(10.,penalty_offset);
    365         Ke->values[5*numdof+5]=+kmax*pow(10.,penalty_offset);
    366 
    367         /*Clean up and return*/
    368         return Ke;
    369 }
    370 /*}}}*/
    371 ElementMatrix* Penpair::PenaltyCreateKMatrixMasstransport(IssmDouble kmax){/*{{{*/
    372 
    373         const int numdof=NUMVERTICES*NDOF1;
    374         IssmDouble penalty_factor;
    375 
    376         /*Initialize Element vector and return if necessary*/
    377         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);
    378 
    379         /*recover parameters: */
    380         parameters->FindParam(&penalty_factor,MasstransportPenaltyFactorEnum);
    381 
    382         //Create elementary matrix: add penalty to
    383         Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_factor);
    384         Ke->values[0*numdof+1]=-kmax*pow(10.,penalty_factor);
    385         Ke->values[1*numdof+0]=-kmax*pow(10.,penalty_factor);
    386         Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_factor);
    387 
    388         /*Clean up and return*/
    389         return Ke;
    390 }
    391 /*}}}*/

  • issm/trunk-jpl/src/c/classes/Loads/Penpair.h

    r18237 r18926  
    3131                /*}}}*/
    3232                /*Object virtual functions definitions:{{{ */
    33                 void  Echo();
    34                 void  DeepEcho();
    35                 int   Id();
    36                 int   ObjectEnum();
    3733                Object* copy();
     34                void     DeepEcho();
     35                void     Echo();
     36                int      Id();
     37                int      ObjectEnum();
    3838                /*}}}*/
    3939                /*Update virtual functions resolution: {{{*/
    40                 void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
    41                 void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrow, int ncols,int name, int type){_error_("Not implemented yet!");}
    42                 void  InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("Not implemented yet!");}
    4340                void  InputUpdateFromConstant(IssmDouble constant, int name);
    4441                void  InputUpdateFromConstant(int constant, int name);
    4542                void  InputUpdateFromConstant(bool constant, int name);
    4643                void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
     44                void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrow, int ncols,int name, int type){_error_("Not implemented yet!");}
     45                void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
     46                void  InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("Not implemented yet!");}
    4747                /*}}}*/
    4848                        /*Load virtual functions definitions: {{{*/
    4949                void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
     50                void  CreateJacobianMatrix(Matrix<IssmDouble>* Jff);
     51                void  CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
     52                void  CreatePVector(Vector<IssmDouble>* pf);
     53                void  GetNodesLidList(int* lidlist);
     54                void  GetNodesSidList(int* sidlist);
     55                int   GetNumberOfNodes(void);
     56                bool  InAnalysis(int analysis_type);
     57                bool  IsPenalty(void);
     58                void  PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax);
     59                void  PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff,Matrix<IssmDouble>* Kfs,IssmDouble kmax);
     60                void  PenaltyCreatePVector(Vector<IssmDouble>* pf, IssmDouble kmax);
    5061                void  ResetHooks();
    5162                void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    52                 void  CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
    53                 void  CreatePVector(Vector<IssmDouble>* pf);
    54                 void  CreateJacobianMatrix(Matrix<IssmDouble>* Jff);
    55                 void  GetNodesSidList(int* sidlist);
    56                 void  GetNodesLidList(int* lidlist);
    57                 int   GetNumberOfNodes(void);
    58                 bool  IsPenalty(void);
    59                 void  PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff,Matrix<IssmDouble>* Kfs,IssmDouble kmax);
    60                 void  PenaltyCreatePVector(Vector<IssmDouble>* pf, IssmDouble kmax);
    61                 void  PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax);
    6263                void  SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
    63                 bool  InAnalysis(int analysis_type);
    6464                /*}}}*/
    6565                        /*Penpair management: {{{*/
     66                ElementMatrix* PenaltyCreateKMatrixMasstransport(IssmDouble kmax);
     67                ElementMatrix* PenaltyCreateKMatrixStressbalanceFS(IssmDouble kmax);
    6668                ElementMatrix* PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax);
    6769                ElementMatrix* PenaltyCreateKMatrixStressbalanceSSAHO(IssmDouble kmax);
    68                 ElementMatrix* PenaltyCreateKMatrixStressbalanceFS(IssmDouble kmax);
    69                 ElementMatrix* PenaltyCreateKMatrixMasstransport(IssmDouble kmax);
    7070                /*}}}*/
    7171};

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

    r18237 r18926  
    108108
    109109/*Object virtual functions definitions:*/
    110 void Riftfront::Echo(void){/*{{{*/
     110Object* Riftfront::copy() {/*{{{*/
     111
     112        Riftfront* riftfront=NULL;
     113
     114        riftfront=new Riftfront();
     115
     116        /*copy fields: */
     117        riftfront->id=this->id;
     118        riftfront->analysis_type=this->analysis_type;
     119        riftfront->type=this->type;
     120        riftfront->fill=this->fill;
     121        riftfront->friction=this->friction;
     122        riftfront->fractionincrement=this->fractionincrement;
     123        riftfront->shelf=this->shelf;
     124
     125        /*point parameters: */
     126        riftfront->parameters=this->parameters;
     127
     128        /*now deal with hooks and objects: */
     129        riftfront->hnodes=(Hook*)this->hnodes->copy();
     130        riftfront->helements=(Hook*)this->helements->copy();
     131        riftfront->hmatpar=(Hook*)this->hmatpar->copy();
     132
     133        /*corresponding fields*/
     134        riftfront->nodes   =(Node**)riftfront->hnodes->deliverp();
     135        riftfront->elements=(Element**)riftfront->helements->deliverp();
     136        riftfront->matpar  =(Matpar*)riftfront->hmatpar->delivers();
     137
     138        /*internal data: */
     139        riftfront->penalty_lock=this->penalty_lock;
     140        riftfront->active=this->active;
     141        riftfront->frozen=this->frozen;
     142        riftfront->state=this->state;
     143        riftfront->counter=this->counter;
     144        riftfront->prestable=this->prestable;
     145        riftfront->material_converged=this->material_converged;
     146        riftfront->normal[0]=this->normal[0];
     147        riftfront->normal[1]=this->normal[1];
     148        riftfront->length=this->length;
     149        riftfront->fraction=this->fraction;
     150
     151        return riftfront;
     152
     153}
     154/*}}}*/
     155void    Riftfront::DeepEcho(void){/*{{{*/
     156
     157        _printf_("Riftfront:\n");
     158        _printf_("   id: " << id << "\n");
     159        _printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
     160        hnodes->DeepEcho();
     161        helements->DeepEcho();
     162        hmatpar->DeepEcho();
     163        _printf_("   parameters\n");
     164        if(parameters)parameters->DeepEcho();
     165}
     166/*}}}*/
     167void    Riftfront::Echo(void){/*{{{*/
    111168
    112169        _printf_("Riftfront:\n");
    … …  
    134191}
    135192/*}}}*/
    136 void Riftfront::DeepEcho(void){/*{{{*/
    137 
    138         _printf_("Riftfront:\n");
    139         _printf_("   id: " << id << "\n");
    140         _printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
    141         hnodes->DeepEcho();
    142         helements->DeepEcho();
    143         hmatpar->DeepEcho();
    144         _printf_("   parameters\n");
    145         if(parameters)parameters->DeepEcho();
    146 }
    147 /*}}}*/
    148 int    Riftfront::Id(void){ return id; }/*{{{*/
    149 /*}}}*/
    150 int Riftfront::ObjectEnum(void){/*{{{*/
     193int     Riftfront::Id(void){ return id; }/*{{{*/
     194/*}}}*/
     195int     Riftfront::ObjectEnum(void){/*{{{*/
    151196
    152197        return RiftfrontEnum;
    153 
    154 }
    155 /*}}}*/
    156 Object* Riftfront::copy() {/*{{{*/
    157 
    158         Riftfront* riftfront=NULL;
    159 
    160         riftfront=new Riftfront();
    161 
    162         /*copy fields: */
    163         riftfront->id=this->id;
    164         riftfront->analysis_type=this->analysis_type;
    165         riftfront->type=this->type;
    166         riftfront->fill=this->fill;
    167         riftfront->friction=this->friction;
    168         riftfront->fractionincrement=this->fractionincrement;
    169         riftfront->shelf=this->shelf;
    170 
    171         /*point parameters: */
    172         riftfront->parameters=this->parameters;
    173 
    174         /*now deal with hooks and objects: */
    175         riftfront->hnodes=(Hook*)this->hnodes->copy();
    176         riftfront->helements=(Hook*)this->helements->copy();
    177         riftfront->hmatpar=(Hook*)this->hmatpar->copy();
    178 
    179         /*corresponding fields*/
    180         riftfront->nodes   =(Node**)riftfront->hnodes->deliverp();
    181         riftfront->elements=(Element**)riftfront->helements->deliverp();
    182         riftfront->matpar  =(Matpar*)riftfront->hmatpar->delivers();
    183 
    184         /*internal data: */
    185         riftfront->penalty_lock=this->penalty_lock;
    186         riftfront->active=this->active;
    187         riftfront->frozen=this->frozen;
    188         riftfront->state=this->state;
    189         riftfront->counter=this->counter;
    190         riftfront->prestable=this->prestable;
    191         riftfront->material_converged=this->material_converged;
    192         riftfront->normal[0]=this->normal[0];
    193         riftfront->normal[1]=this->normal[1];
    194         riftfront->length=this->length;
    195         riftfront->fraction=this->fraction;
    196 
    197         return riftfront;
    198198
    199199}
    … …  
    234234}
    235235/*}}}*/
    236 void  Riftfront::ResetHooks(){/*{{{*/
    237 
    238         this->nodes=NULL;
    239         this->elements=NULL;
    240         this->matpar=NULL;
    241         this->parameters=NULL;
    242 
    243         /*Get Element type*/
    244         this->hnodes->reset();
    245         this->helements->reset();
    246         this->hmatpar->reset();
    247 
     236void  Riftfront::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
     237        /*do nothing: */
     238        return;
     239}
     240/*}}}*/
     241void  Riftfront::CreatePVector(Vector<IssmDouble>* pf){/*{{{*/
     242        /*do nothing: */
     243        return;
     244}
     245/*}}}*/
     246void  Riftfront::GetNodesLidList(int* lidlist){/*{{{*/
     247
     248        _assert_(lidlist);
     249        _assert_(nodes);
     250
     251        for(int i=0;i<NUMVERTICES;i++) lidlist[i]=nodes[i]->Lid();
     252}
     253/*}}}*/
     254void  Riftfront::GetNodesSidList(int* sidlist){/*{{{*/
     255
     256        _assert_(sidlist);
     257        _assert_(nodes);
     258
     259        for(int i=0;i<NUMVERTICES;i++) sidlist[i]=nodes[i]->Sid();
     260}
     261/*}}}*/
     262int   Riftfront::GetNumberOfNodes(void){/*{{{*/
     263
     264        return NUMVERTICES;
     265}
     266/*}}}*/
     267bool  Riftfront::InAnalysis(int in_analysis_type){/*{{{*/
     268        if (in_analysis_type==this->analysis_type) return true;
     269        else return false;
    248270}
    249271/*}}}*/
    250272bool  Riftfront::IsPenalty(void){/*{{{*/
    251273        return true;
    252 }
    253 /*}}}*/
    254 void  Riftfront::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
    255 
    256274}
    257275/*}}}*/
    … …  
    306324}
    307325/*}}}*/
    308 void  Riftfront::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
    309         /*do nothing: */
    310         return;
    311 }
    312 /*}}}*/
    313 void  Riftfront::CreatePVector(Vector<IssmDouble>* pf){/*{{{*/
    314         /*do nothing: */
    315         return;
    316 }
    317 /*}}}*/
    318 void  Riftfront::GetNodesSidList(int* sidlist){/*{{{*/
    319 
    320         _assert_(sidlist);
    321         _assert_(nodes);
    322 
    323         for(int i=0;i<NUMVERTICES;i++) sidlist[i]=nodes[i]->Sid();
    324 }
    325 /*}}}*/
    326 void  Riftfront::GetNodesLidList(int* lidlist){/*{{{*/
    327 
    328         _assert_(lidlist);
    329         _assert_(nodes);
    330 
    331         for(int i=0;i<NUMVERTICES;i++) lidlist[i]=nodes[i]->Lid();
    332 }
    333 /*}}}*/
    334 int   Riftfront::GetNumberOfNodes(void){/*{{{*/
    335 
    336         return NUMVERTICES;
    337 }
    338 /*}}}*/
    339 bool  Riftfront::InAnalysis(int in_analysis_type){/*{{{*/
    340         if (in_analysis_type==this->analysis_type) return true;
    341         else return false;
     326void  Riftfront::ResetHooks(){/*{{{*/
     327
     328        this->nodes=NULL;
     329        this->elements=NULL;
     330        this->matpar=NULL;
     331        this->parameters=NULL;
     332
     333        /*Get Element type*/
     334        this->hnodes->reset();
     335        this->helements->reset();
     336        this->hmatpar->reset();
     337
     338}
     339/*}}}*/
     340void  Riftfront::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
     341
    342342}
    343343/*}}}*/

  • issm/trunk-jpl/src/c/classes/Loads/Riftfront.h

    r18237 r18926  
    5757                /*}}}*/
    5858                /*Object virtual functions definitions:{{{ */
    59                 void  Echo();
    60                 void  DeepEcho();
    61                 int   Id();
    62                 int   ObjectEnum();
    63                 Object* copy();
     59                Object*  copy();
     60                void     DeepEcho();
     61                void     Echo();
     62                int      Id();
     63                int      ObjectEnum();
    6464                /*}}}*/
    6565                /*Update virtual functions resolution: {{{*/
    66                 void    InputUpdateFromVector(IssmDouble* vector, int name, int type);
    67                 void    InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows,int ncols, int name, int type){_error_("Not implemented yet!");}
    68                 void    InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("Not implemented yet!");}
    6966                void    InputUpdateFromConstant(IssmDouble constant, int name);
    7067                void    InputUpdateFromConstant(int constant, int name){_error_("Not implemented yet!");}
    7168                void    InputUpdateFromConstant(bool constant, int name);
    7269                void    InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
     70                void    InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows,int ncols, int name, int type){_error_("Not implemented yet!");}
     71                void    InputUpdateFromVector(IssmDouble* vector, int name, int type);
     72                void    InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("Not implemented yet!");}
    7373                /*}}}*/
    7474                /*Load virtual functions definitions: {{{*/
    7575                void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    76                 void  ResetHooks();
    77                 void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
     76                void  CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("Not implemented yet");};
    7877                void  CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
    7978                void  CreatePVector(Vector<IssmDouble>* pf);
    80                 void  CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("Not implemented yet");};
     79                void  GetNodesLidList(int* lidlist);
    8180                void  GetNodesSidList(int* sidlist);
    82                 void  GetNodesLidList(int* lidlist);
    8381                int   GetNumberOfNodes(void);
     82                bool  InAnalysis(int analysis_type);
    8483                bool  IsPenalty(void);
    8584                void  PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){_error_("Not implemented yet");};
    8685                void  PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* kfs, IssmDouble kmax);
    8786                void  PenaltyCreatePVector(Vector<IssmDouble>* pf, IssmDouble kmax);
     87                void  ResetHooks();
     88                void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    8889                void  SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
    89                 bool  InAnalysis(int analysis_type);
    9090                /*}}}*/
    9191                /*Riftfront specific routines: {{{*/
     92                int            Constrain(int* punstable);
     93                void           FreezeConstraints(void);
     94                bool           IsFrozen(void);
    9295                ElementMatrix* PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax);
    9396                ElementVector* PenaltyCreatePVectorStressbalanceHoriz(IssmDouble kmax);
    94                 int   Constrain(int* punstable);
    95                 void  FreezeConstraints(void);
    96                 bool  IsFrozen(void);
    9797                /*}}}*/
    9898};
Note: See TracChangeset for help on using the changeset viewer.