source: issm/oecreview/Archive/16554-17801/ISSM-17515-17516.diff@ 17802

Last change on this file since 17802 was 17802, checked in by Mathieu Morlighem, 11 years ago

Added archives
File size: 73.3 KB

  • ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp

     
    213213        IssmDouble *xyz_list = NULL;
    214214
    215215        /*Fetch number of nodes and dof for this finite element*/
    216         int vnumnodes = element->GetNumberOfNodesVelocity();
    217         int pnumnodes = element->GetNumberOfNodesPressure();
     216        int vnumnodes = element->NumberofNodesVelocity();
     217        int pnumnodes = element->NumberofNodesPressure();
    218218        int numdof    = vnumnodes*3 + pnumnodes;
    219219
    220220        /*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/
    … …  
    319319        }
    320320
    321321        /*Fetch number of nodes and dof for this finite element*/
    322         int vnumnodes = element->GetNumberOfNodesVelocity();
    323         int pnumnodes = element->GetNumberOfNodesPressure();
     322        int vnumnodes = element->NumberofNodesVelocity();
     323        int pnumnodes = element->NumberofNodesPressure();
    324324
    325325        /*Prepare coordinate system list*/
    326326        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    … …  
    938938        IssmDouble   FSreconditioning;
    939939
    940940        /*Fetch number of nodes and dof for this finite element*/
    941         int vnumnodes = element->GetNumberOfNodesVelocity();
    942         int pnumnodes = element->GetNumberOfNodesPressure();
     941        int vnumnodes = element->NumberofNodesVelocity();
     942        int pnumnodes = element->NumberofNodesPressure();
    943943        int vnumdof   = vnumnodes*3;
    944944        int pnumdof   = pnumnodes*1;
    945945

  • ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

     
    194194        IssmDouble *vertex_pairing=NULL;
    195195        IssmDouble *nodeonbed=NULL;
    196196        iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum);
    197         iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
     197        if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
    198198
    199199        for(int i=0;i<numvertex_pairing;i++){
    200200
    … …  
    204204                        _assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);
    205205
    206206                        /*Skip if one of the two is not on the bed*/
    207                         if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
     207                        if(iomodel->meshtype!=Mesh2DhorizontalEnum){
     208                                if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
     209                        }
    208210
    209211                        /*Get node ids*/
    210212                        penpair_ids[0]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+0]);

  • ../trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp

     
    3434        iomodel->FetchDataToInput(elements,SurfaceEnum);
    3535        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
    3636        iomodel->FetchDataToInput(elements,VxEnum);
    37         iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
     37        if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
    3838        if(iomodel->meshtype==Mesh3DEnum){
    3939                iomodel->FetchDataToInput(elements,VzEnum);
    4040                iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
    … …  
    6868        IssmDouble *vertex_pairing=NULL;
    6969        IssmDouble *nodeonsurface=NULL;
    7070        iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum);
    71         iomodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum);
     71        if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum);
    7272        for(int i=0;i<numvertex_pairing;i++){
    7373
    7474                if(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+0])-1]){
    … …  
    7777                        _assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);
    7878
    7979                        /*Skip if one of the two is not on the bed*/
    80                         if(!(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
     80                        if(iomodel->meshtype!=Mesh2DhorizontalEnum){
     81                                if(!(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
     82                        }
    8183
    8284                        /*Get node ids*/
    8385                        penpair_ids[0]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+0]);

  • ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

     
    27532753        }
    27542754
    27552755        /*Fetch number of nodes and dof for this finite element*/
    2756         int vnumnodes = element->GetNumberOfNodesVelocity();
    2757         int pnumnodes = element->GetNumberOfNodesPressure();
     2756        int vnumnodes = element->NumberofNodesVelocity();
     2757        int pnumnodes = element->NumberofNodesPressure();
    27582758
    27592759        /*Initialize output vector*/
    27602760        ElementVector* de = element->NewElementVector(FSvelocityEnum);
    … …  
    27832783        IssmDouble *xyz_list = NULL;
    27842784
    27852785        /*Fetch number of nodes and dof for this finite element*/
    2786         int vnumnodes = element->GetNumberOfNodesVelocity();
    2787         int pnumnodes = element->GetNumberOfNodesPressure();
     2786        int vnumnodes = element->NumberofNodesVelocity();
     2787        int pnumnodes = element->NumberofNodesPressure();
    27882788        int numdof    = vnumnodes*NDOF3 + pnumnodes*NDOF1;
    27892789
    27902790        /*Prepare coordinate system list*/
    … …  
    28852885        }
    28862886
    28872887        /*Fetch number of nodes and dof for this finite element*/
    2888         int vnumnodes = element->GetNumberOfNodesVelocity();
    2889         int pnumnodes = element->GetNumberOfNodesPressure();
     2888        int vnumnodes = element->NumberofNodesVelocity();
     2889        int pnumnodes = element->NumberofNodesPressure();
    28902890        int numdof    = vnumnodes*dim + pnumnodes;
    28912891        int bsize     = epssize + 2;
    28922892
    … …  
    29692969        }
    29702970
    29712971        /*Fetch number of nodes and dof for this finite element*/
    2972         int vnumnodes = element->GetNumberOfNodesVelocity();
    2973         int pnumnodes = element->GetNumberOfNodesPressure();
     2972        int vnumnodes = element->NumberofNodesVelocity();
     2973        int pnumnodes = element->NumberofNodesPressure();
    29742974        int numdof    = vnumnodes*dim + pnumnodes;
    29752975
    29762976        /*Initialize Element matrix and vectors*/
    … …  
    30643064        }
    30653065
    30663066        /*Fetch number of nodes and dof for this finite element*/
    3067         int vnumnodes = element->GetNumberOfNodesVelocity();
    3068         int pnumnodes = element->GetNumberOfNodesPressure();
     3067        int vnumnodes = element->NumberofNodesVelocity();
     3068        int pnumnodes = element->NumberofNodesPressure();
    30693069        int numdof    = vnumnodes*dim + pnumnodes;
    30703070
    30713071        /*Initialize Element matrix and vectors*/
    … …  
    31253125        }
    31263126
    31273127        /*Fetch number of nodes and dof for this finite element*/
    3128         int vnumnodes = element->GetNumberOfNodesVelocity();
    3129         int pnumnodes = element->GetNumberOfNodesPressure();
     3128        int vnumnodes = element->NumberofNodesVelocity();
     3129        int pnumnodes = element->NumberofNodesPressure();
    31303130
    31313131        /*Prepare coordinate system list*/
    31323132        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    … …  
    32083208        }
    32093209
    32103210        /*Fetch number of nodes and dof for this finite element*/
    3211         int vnumnodes = element->GetNumberOfNodesVelocity();
    3212         int pnumnodes = element->GetNumberOfNodesPressure();
     3211        int vnumnodes = element->NumberofNodesVelocity();
     3212        int pnumnodes = element->NumberofNodesPressure();
    32133213
    32143214        /*Prepare coordinate system list*/
    32153215        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    … …  
    32873287        }
    32883288
    32893289        /*Fetch number of nodes and dof for this finite element*/
    3290         int vnumnodes = element->GetNumberOfNodesVelocity();
    3291         int pnumnodes = element->GetNumberOfNodesPressure();
     3290        int vnumnodes = element->NumberofNodesVelocity();
     3291        int pnumnodes = element->NumberofNodesPressure();
    32923292
    32933293        /*Prepare coordinate system list*/
    32943294        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    … …  
    33643364        }
    33653365
    33663366        /*Fetch number of nodes and dof for this finite element*/
    3367         int vnumnodes   = element->GetNumberOfNodesVelocity();
    3368         int pnumnodes   = element->GetNumberOfNodesPressure();
     3367        int vnumnodes   = element->NumberofNodesVelocity();
     3368        int pnumnodes   = element->NumberofNodesPressure();
    33693369        int numvertices = element->GetNumberOfVertices();
    33703370
    33713371        /*Prepare coordinate system list*/
    … …  
    36583658         */
    36593659
    36603660        /*Fetch number of nodes for this finite element*/
    3661         int pnumnodes = element->GetNumberOfNodesPressure();
    3662         int vnumnodes = element->GetNumberOfNodesVelocity();
     3661        int pnumnodes = element->NumberofNodesPressure();
     3662        int vnumnodes = element->NumberofNodesVelocity();
    36633663        int pnumdof   = pnumnodes;
    36643664        int vnumdof   = vnumnodes*dim;
    36653665
    … …  
    37853785        }
    37863786
    37873787        /*Fetch number of nodes and dof for this finite element*/
    3788         int vnumnodes = element->GetNumberOfNodesVelocity();
    3789         int pnumnodes = element->GetNumberOfNodesPressure();
     3788        int vnumnodes = element->NumberofNodesVelocity();
     3789        int pnumnodes = element->NumberofNodesPressure();
    37903790        int vnumdof   = vnumnodes*dim;
    37913791        int pnumdof   = pnumnodes*1;
    37923792
    … …  
    43204320        element->GetInputValue(&approximation,ApproximationEnum);
    43214321        if(element->IsFloating() || !element->IsOnBed()) return NULL;
    43224322
    4323         int vnumnodes = element->GetNumberOfNodesVelocity();
    4324         int pnumnodes = element->GetNumberOfNodesPressure();
     4323        int vnumnodes = element->NumberofNodesVelocity();
     4324        int pnumnodes = element->NumberofNodesPressure();
    43254325        int numnodes  = 2*vnumnodes-1+pnumnodes;
    43264326
    43274327        /*Prepare node list*/
    … …  
    44384438        Element* pentabase=element->GetBasalElement();
    44394439        Element* basaltria=pentabase->SpawnBasalElement();
    44404440
    4441         int vnumnodes = element->GetNumberOfNodesVelocity();
    4442         int pnumnodes = element->GetNumberOfNodesPressure();
     4441        int vnumnodes = element->NumberofNodesVelocity();
     4442        int pnumnodes = element->NumberofNodesPressure();
    44434443        int numnodes  = 2*vnumnodes-1+pnumnodes;
    44444444
    44454445        /*Prepare node list*/
    … …  
    46074607        element->GetInputValue(&approximation,ApproximationEnum);
    46084608        if(approximation!=HOFSApproximationEnum) return NULL;
    46094609
    4610         int vnumnodes = element->GetNumberOfNodesVelocity();
    4611         int pnumnodes = element->GetNumberOfNodesPressure();
     4610        int vnumnodes = element->NumberofNodesVelocity();
     4611        int pnumnodes = element->NumberofNodesPressure();
    46124612        int numnodes  = vnumnodes+pnumnodes;
    46134613
    46144614        /*Prepare coordinate system list*/
    … …  
    46864686        /*Initialize Element vector and return if necessary*/
    46874687        element->GetInputValue(&approximation,ApproximationEnum);
    46884688        if(approximation!=HOFSApproximationEnum) return NULL;
    4689         int   vnumnodes = element->GetNumberOfNodesVelocity();
    4690         int   pnumnodes = element->GetNumberOfNodesPressure();
     4689        int   vnumnodes = element->NumberofNodesVelocity();
     4690        int   pnumnodes = element->NumberofNodesPressure();
    46914691
    46924692        /*Prepare coordinate system list*/
    46934693        int*   cs_list   = xNew<int>(vnumnodes+pnumnodes);
    … …  
    47704770        if(!element->IsOnBed() || element->IsFloating()) return NULL;
    47714771        element->GetInputValue(&approximation,ApproximationEnum);
    47724772        if(approximation!=SSAFSApproximationEnum) return NULL;
    4773         int vnumnodes = element->GetNumberOfNodesVelocity();
    4774         int pnumnodes = element->GetNumberOfNodesPressure();
     4773        int vnumnodes = element->NumberofNodesVelocity();
     4774        int pnumnodes = element->NumberofNodesPressure();
    47754775
    47764776        /*Prepare coordinate system list*/
    47774777        int* cs_list     = xNew<int>(vnumnodes+pnumnodes);
    … …  
    48464846        /*Initialize Element vector and return if necessary*/
    48474847        element->GetInputValue(&approximation,ApproximationEnum);
    48484848        if(approximation!=SSAFSApproximationEnum) return NULL;
    4849         int vnumnodes = element->GetNumberOfNodesVelocity();
    4850         int pnumnodes = element->GetNumberOfNodesPressure();
     4849        int vnumnodes = element->NumberofNodesVelocity();
     4850        int pnumnodes = element->NumberofNodesPressure();
    48514851
    48524852        /*Prepare coordinate system list*/
    48534853        int* cs_list = xNew<int>(vnumnodes+pnumnodes);

  • ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp

     
    7272}/*}}}*/
    7373void DamageEvolutionAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
    7474
    75         iomodel->FetchData(1,MeshVertexonbedEnum);
     75        if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(1,MeshVertexonbedEnum);
    7676        ::CreateNodes(nodes,iomodel,DamageEvolutionAnalysisEnum,P1Enum);
    7777        iomodel->DeleteData(1,MeshVertexonbedEnum);
    7878}/*}}}*/
    … …  
    8686}/*}}}*/
    8787void DamageEvolutionAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
    8888
    89         /*create penalties for nodes: no node can have a damage > 1*/
    90         iomodel->FetchData(1,DamageSpcdamageEnum);
    91         CreateSingleNodeToElementConnectivity(iomodel);
     89        /*Nothing for now*/
    9290
    93         for(int i=0;i<iomodel->numberofvertices;i++){
    94 
    95                 /*keep only this partition's nodes:*/
    96                 if(iomodel->my_vertices[i]){
    97                         if (xIsNan<IssmDouble>(iomodel->Data(DamageSpcdamageEnum)[i])){ //No penalty applied on spc nodes!
    98                                 loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,DamageEvolutionAnalysisEnum));
    99                         }
    100                 }
    101         }
    102         iomodel->DeleteData(1,DamageSpcdamageEnum);
    103 
    10491}/*}}}*/
    10592
    10693/*Finite Element Analysis*/

  • ../trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp

     
    5959        IssmDouble *vertex_pairing=NULL;
    6060        IssmDouble *nodeonbed=NULL;
    6161        iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum);
    62         iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
     62        if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
    6363        for(int i=0;i<numvertex_pairing;i++){
    6464
    6565                if(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+0])-1]){
    … …  
    6868                        _assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);
    6969
    7070                        /*Skip if one of the two is not on the bed*/
    71                         if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
     71                        if(iomodel->meshtype!=Mesh2DhorizontalEnum){
     72                                if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
     73                        }
    7274
    7375                        /*Get node ids*/
    7476                        penpair_ids[0]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+0]);

  • ../trunk-jpl/src/c/solutionsequences/solutionsequence_damage_nonlinear.cpp

     
    1 /*
    2  * \brief: solutionsequence_damage_nonlinear.cpp: core of the damage solution
    3  */
    4 
    5 #include "../toolkits/toolkits.h"
    6 #include "../classes/classes.h"
    7 #include "../shared/shared.h"
    8 #include "../modules/modules.h"
    9 
    10 void solutionsequence_damage_nonlinear(FemModel* femmodel){
    11 
    12         /*solution : */
    13         Vector<IssmDouble>* Dg=NULL;
    14         Vector<IssmDouble>* Df=NULL;
    15         Vector<IssmDouble>* Df_old=NULL;
    16         Vector<IssmDouble>* ys=NULL;
    17 
    18         /*intermediary: */
    19         Matrix<IssmDouble>* Kff=NULL;
    20         Matrix<IssmDouble>* Kfs=NULL;
    21         Vector<IssmDouble>* pf=NULL;
    22         Vector<IssmDouble>* df=NULL;
    23 
    24         bool converged;
    25         int constraints_converged;
    26         int num_unstable_constraints;
    27         int count;
    28         int damage_penalty_threshold;
    29         int damage_maxiter;
    30 
    31         /*parameters:*/
    32         int  configuration_type;
    33 
    34         /*Recover parameters: */
    35         femmodel->parameters->FindParam(&damage_penalty_threshold,DamagePenaltyThresholdEnum);
    36         femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    37         femmodel->parameters->FindParam(&damage_maxiter,DamageMaxiterEnum);
    38 
    39         count=1;
    40         converged=false;
    41 
    42         InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum);
    43         InputUpdateFromConstantx(femmodel,false,ConvergedEnum);
    44         femmodel->UpdateConstraintsx();
    45 
    46         for(;;){
    47 
    48                 delete Df_old; Df_old=Df;
    49                 SystemMatricesx(&Kff, &Kfs, &pf,&df, NULL,femmodel);
    50                 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
    51                 Reduceloadx(pf, Kfs, ys); delete Kfs;
    52                 Solverx(&Df, Kff, pf,Df_old, df, femmodel->parameters);
    53                 delete Kff;delete pf;delete Dg; delete df;
    54                 Mergesolutionfromftogx(&Dg, Df,ys,femmodel->nodes,femmodel->parameters); delete ys;
    55                 InputUpdateFromSolutionx(femmodel,Dg);
    56 
    57                 ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
    58 
    59                 if (!converged){
    60                         if(VerboseConvergence()) _printf0_("   #unstable constraints = " << num_unstable_constraints << "\n");
    61                         if (num_unstable_constraints <= damage_penalty_threshold)converged=true;
    62                         if (count>=damage_maxiter){
    63                                 converged=true;
    64                                 _printf0_("   maximum number of iterations (" << damage_maxiter << ") exceeded\n");
    65                         }
    66                 }
    67                 count++;
    68 
    69                 InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
    70 
    71                 if(converged)break;
    72         }
    73 
    74         InputUpdateFromSolutionx(femmodel,Dg);
    75 
    76         /*Free ressources: */
    77         delete Dg;
    78         delete Df;
    79         delete Df_old;
    80 }

  • ../trunk-jpl/src/c/solutionsequences/solutionsequences.h

     
    1212#include "../shared/Numerics/types.h"
    1313
    1414void solutionsequence_thermal_nonlinear(FemModel* femmodel);
    15 void solutionsequence_damage_nonlinear(FemModel* femmodel);
    1615void solutionsequence_hydro_nonlinear(FemModel* femmodel);
    1716void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads);
    1817void solutionsequence_newton(FemModel* femmodel);

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

     
    242242                case HydrologyDCInefficientAnalysisEnum:
    243243                        Ke=PenaltyCreateKMatrixHydrologyDCInefficient(kmax);
    244244                        break;
    245                 case DamageEvolutionAnalysisEnum:
    246                         Ke=PenaltyCreateKMatrixDamageEvolution(kmax);
    247                         break;
    248245                default:
    249246                        _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    250247        }
    … …  
    276273                case HydrologyDCInefficientAnalysisEnum:
    277274                        pe=PenaltyCreatePVectorHydrologyDCInefficient(kmax);
    278275                        break;
    279                 case DamageEvolutionAnalysisEnum:
    280                         pe=PenaltyCreatePVectorDamageEvolution(kmax);
    281                         break;
    282276                default:
    283277                        _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    284278        }
    … …  
    416410                ConstraintActivateHydrologyDCInefficient(punstable);
    417411                return;
    418412        }
    419         else if (analysis_type==DamageEvolutionAnalysisEnum){
    420                 ConstraintActivateDamageEvolution(punstable);
    421                 return;
    422         }
    423 
    424413        else{
    425414                _error_("analysis: " << EnumToStringx(analysis_type) << " not supported yet");
    426415        }
    … …  
    709698        return pe;
    710699}
    711700/*}}}*/
    712 /*FUNCTION Pengrid::ConstraintActivateDamageEvolution {{{*/
    713 void  Pengrid::ConstraintActivateDamageEvolution(int* punstable){
    714 
    715         //   The penalty is stable if it doesn't change during to successive iterations.   
    716         IssmDouble max_damage;
    717         IssmDouble damage;
    718         int        new_active;
    719         int        unstable=0;
    720         int        penalty_lock;
    721 
    722         /*check that pengrid is not a clone (penalty to be added only once)*/
    723         if (node->IsClone()){
    724                 unstable=0;
    725                 *punstable=unstable;
    726                 return;
    727         }
    728 
    729         //First recover damage  using the element: */
    730         element->GetInputValue(&damage,node,DamageDEnum);
    731 
    732         //Recover our data:
    733         parameters->FindParam(&penalty_lock,DamagePenaltyLockEnum);
    734         parameters->FindParam(&max_damage,DamageMaxDamageEnum);
    735        
    736         //Figure out if damage>max_damage, in which case, this penalty needs to be activated.
    737         //Would need to do the same for damage<0 if penalties are used.  For now, ConstraintStatex
    738         //is not called in solutionsequence_damage_nonlinear, so no penalties are applied.
    739 
    740         if (damage>max_damage){
    741                 new_active=1;
    742         }
    743         else{
    744                 new_active=0;
    745         }
    746 
    747         //Figure out stability of this penalty
    748         if (active==new_active){
    749                 unstable=0;
    750         }
    751         else{
    752                 unstable=1;
    753                 if(penalty_lock)zigzag_counter++;
    754         }
    755 
    756         /*If penalty keeps zigzagging more than penalty_lock times: */
    757         if(penalty_lock){
    758                 if(zigzag_counter>penalty_lock){
    759                         unstable=0;
    760                         active=1;
    761                 }
    762         }
    763 
    764         //Set penalty flag
    765         active=new_active;
    766 
    767         //*Assign output pointers:*/
    768         *punstable=unstable;
    769 }
    770 /*}}}*/
    771 /*FUNCTION Pengrid::PenaltyCreateKMatrixDamageEvolution {{{*/
    772 ElementMatrix* Pengrid::PenaltyCreateKMatrixDamageEvolution(IssmDouble kmax){
    773 
    774         IssmDouble    penalty_factor;
    775 
    776         /*Initialize Element matrix and return if necessary*/
    777         if(!this->active) return NULL;
    778         ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters);
    779 
    780         /*recover parameters: */
    781         parameters->FindParam(&penalty_factor,DamagePenaltyFactorEnum);
    782 
    783         Ke->values[0]=kmax*pow(10.,penalty_factor);
    784 
    785         /*Clean up and return*/
    786         return Ke;
    787 }
    788 /*}}}*/
    789 /*FUNCTION Pengrid::PenaltyCreatePVectorDamageEvolution {{{*/
    790 ElementVector* Pengrid::PenaltyCreatePVectorDamageEvolution(IssmDouble kmax){
    791 
    792         IssmDouble penalty_factor;
    793         IssmDouble max_damage;
    794 
    795         /*Initialize Element matrix and return if necessary*/
    796         if(!this->active) return NULL;
    797         ElementVector* pe=new ElementVector(&node,1,this->parameters);
    798 
    799         //Recover our data:
    800         parameters->FindParam(&penalty_factor,DamagePenaltyFactorEnum);
    801         parameters->FindParam(&max_damage,DamageMaxDamageEnum);
    802 
    803         //right hand side penalizes to max_damage
    804         pe->values[0]=kmax*pow(10.,penalty_factor)*max_damage;
    805 
    806         /*Clean up and return*/
    807         return pe;
    808 }
    809 /*}}}*/
    810701/*FUNCTION Pengrid::CreatePVectorHydrologyDCInefficient {{{*/
    811702ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){
    812703

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

     
    8686                ElementVector* PenaltyCreatePVectorThermal(IssmDouble kmax);
    8787                ElementVector* PenaltyCreatePVectorMelting(IssmDouble kmax);
    8888                void           ConstraintActivateThermal(int* punstable);
    89                 ElementMatrix* PenaltyCreateKMatrixDamageEvolution(IssmDouble kmax);
    90                 ElementVector* PenaltyCreatePVectorDamageEvolution(IssmDouble kmax);
    91                 void           ConstraintActivateDamageEvolution(int* punstable);
    9289                ElementMatrix* PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax);
    9390                ElementVector* PenaltyCreatePVectorHydrologyDCInefficient(IssmDouble kmax);
    9491                void           ConstraintActivateHydrologyDCInefficient(int* punstable);

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

     
    3737class Element: public Object,public Update{
    3838
    3939        public:
     40                int          id;
     41                int          sid;
    4042                Inputs      *inputs;
    4143                Node       **nodes;
    4244                Vertex     **vertices;
    … …  
    5456                /* bool       AllActive(void); */
    5557                /* bool       AnyActive(void); */
    5658                void       CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array);
     59                void       Echo();
     60                void       DeepEcho();
    5761                void       DeleteMaterials(void);
    5862                IssmDouble Divergence(void);
     63                void       ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure);
     64                void       EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
     65                IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
     66                IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);
    5967                void       FindParam(bool* pvalue,int paramenum);
    6068                void       FindParam(int* pvalue,int paramenum);
    6169                void       FindParam(IssmDouble* pvalue,int paramenum);
    … …  
    8088                void       GetVerticesCoordinates(IssmDouble** xyz_list);
    8189                void       GetVerticesSidList(int* sidlist);
    8290                void       GetVerticesConnectivityList(int* connectivitylist);
     91                bool       HasNodeOnBed();
     92                bool       HasNodeOnSurface();
     93                int        Id();
     94                int        Sid();
    8395                void       InputChangeName(int enum_type,int enum_type_old);
    8496                void       InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code);
    8597                void       InputUpdateFromConstant(IssmDouble constant, int name);
    8698                void       InputUpdateFromConstant(int constant, int name);
    8799                void       InputUpdateFromConstant(bool constant, int name);
     100                bool       IsIceInElement();
    88101                bool         IsInput(int name);
    89102                bool       IsFloating();
    90103                ElementVector*  NewElementVector(int approximation_enum=NoneApproximationEnum);
    91104                ElementMatrix*  NewElementMatrix(int approximation_enum=NoneApproximationEnum);
    92105                ElementMatrix*  NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum);
     106                IssmDouble PureIceEnthalpy(IssmDouble pressure);
    93107                void       ResultInterpolation(int* pinterpolation,int*nodesperelement,int output_enum);
    94108                void       ResultToVector(Vector<IssmDouble>* vector,int output_enum);
    95109                void       ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum);
     110                void       SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
    96111                void       StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    97112                void       StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input);
    98113                void       StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    … …  
    134149                virtual IssmDouble CharacteristicLength(void)=0;
    135150                virtual void       Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
    136151                virtual void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0;
    137                 virtual void       SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;
    138152                virtual void   ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
    139                 virtual void   ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure)=0;
    140                 virtual void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0;
    141                 virtual IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure)=0;
    142                 virtual IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure)=0;
     153
    143154                virtual int    FiniteElement(void)=0;
    144155                virtual IssmDouble MinEdgeLength(IssmDouble* xyz_list)=0;
    145156                virtual void   NodalFunctions(IssmDouble* basis,Gauss* gauss)=0;
    … …  
    153164                virtual void   NormalSection(IssmDouble* normal,IssmDouble* xyz_list)=0;
    154165                virtual void   NormalTop(IssmDouble* normal,IssmDouble* xyz_list)=0;
    155166                virtual void   NormalBase(IssmDouble* normal,IssmDouble* xyz_list)=0;
    156                 virtual IssmDouble PureIceEnthalpy(IssmDouble pressure)=0;
    157167
    158168                virtual Element* GetUpperElement(void)=0;
    159169                virtual Element* GetLowerElement(void)=0;
    … …  
    168178                virtual void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int solutionenum)=0;
    169179                virtual int    GetNodeIndex(Node* node)=0;
    170180                virtual int    GetNumberOfNodes(void)=0;
    171                 virtual int    GetNumberOfNodesVelocity(void)=0;
    172                 virtual int    GetNumberOfNodesPressure(void)=0;
    173181                virtual int    GetNumberOfVertices(void)=0;
    174182
    175                 virtual int    Id()=0;
    176                 virtual int    Sid()=0;
    177183                virtual bool   IsNodeOnShelfFromFlags(IssmDouble* flags)=0;
    178184                virtual bool   IsOnBed()=0;
    179185                virtual bool   IsOnSurface()=0;
    180                 virtual bool   IsIceInElement()=0;
    181186                virtual void   GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating)=0;
    182187                virtual IssmDouble GetGroundedPortion(IssmDouble* xyz_list)=0;
    183188                virtual void   GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;

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

     
    9494/*}}}*/
    9595
    9696/*Other*/
    97 /*FUNCTION Tria::SetwiseNodeConnectivity{{{*/
    98 void Tria::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){
    99 
    100         /*Intermediaries*/
    101         const int numnodes = this->NumberofNodes();
    102 
    103         /*Output */
    104         int d_nz = 0;
    105         int o_nz = 0;
    106 
    107         /*Loop over all nodes*/
    108         for(int i=0;i<numnodes;i++){
    109 
    110                 if(!flags[this->nodes[i]->Lid()]){
    111 
    112                         /*flag current node so that no other element processes it*/
    113                         flags[this->nodes[i]->Lid()]=true;
    114 
    115                         int counter=0;
    116                         while(flagsindices[counter]>=0) counter++;
    117                         flagsindices[counter]=this->nodes[i]->Lid();
    118 
    119                         /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
    120                         switch(set2_enum){
    121                                 case FsetEnum:
    122                                         if(nodes[i]->indexing.fsize){
    123                                                 if(this->nodes[i]->IsClone())
    124                                                  o_nz += 1;
    125                                                 else
    126                                                  d_nz += 1;
    127                                         }
    128                                         break;
    129                                 case GsetEnum:
    130                                         if(nodes[i]->indexing.gsize){
    131                                                 if(this->nodes[i]->IsClone())
    132                                                  o_nz += 1;
    133                                                 else
    134                                                  d_nz += 1;
    135                                         }
    136                                         break;
    137                                 case SsetEnum:
    138                                         if(nodes[i]->indexing.ssize){
    139                                                 if(this->nodes[i]->IsClone())
    140                                                  o_nz += 1;
    141                                                 else
    142                                                  d_nz += 1;
    143                                         }
    144                                         break;
    145                                 default: _error_("not supported");
    146                         }
    147                 }
    148         }
    149 
    150         /*Assign output pointers: */
    151         *pd_nz=d_nz;
    152         *po_nz=o_nz;
    153 }
    154 /*}}}*/
    15597/*FUNCTION Tria::AddBasalInput{{{*/
    15698void  Tria::AddBasalInput(int input_enum,IssmDouble* values, int interpolation_enum){
    15799
    … …  
    291233
    292234}
    293235/*}}}*/
    294 /*FUNCTION Tria::DeepEcho{{{*/
    295 void Tria::DeepEcho(void){
    296 
    297         _printf_("Tria:\n");
    298         _printf_("   id: " << id << "\n");
    299         if(nodes){
    300                 nodes[0]->DeepEcho();
    301                 nodes[1]->DeepEcho();
    302                 nodes[2]->DeepEcho();
    303         }
    304         else _printf_("nodes = NULL\n");
    305 
    306         if (material) material->DeepEcho();
    307         else _printf_("material = NULL\n");
    308 
    309         if (matpar) matpar->DeepEcho();
    310         else _printf_("matpar = NULL\n");
    311 
    312         _printf_("   parameters\n");
    313         if (parameters) parameters->DeepEcho();
    314         else _printf_("parameters = NULL\n");
    315 
    316         _printf_("   inputs\n");
    317         if (inputs) inputs->DeepEcho();
    318         else _printf_("inputs=NULL\n");
    319 
    320         return;
    321 }
    322 /*}}}*/
    323236/*FUNCTION Tria::Delta18oParameterization{{{*/
    324237void  Tria::Delta18oParameterization(void){
    325238
    … …  
    415328        *hz=0.;
    416329}
    417330/*}}}*/
    418 /*FUNCTION Tria::Echo{{{*/
    419 void Tria::Echo(void){
    420         _printf_("Tria:\n");
    421         _printf_("   id: " << id << "\n");
    422         if(nodes){
    423                 nodes[0]->Echo();
    424                 nodes[1]->Echo();
    425                 nodes[2]->Echo();
    426         }
    427         else _printf_("nodes = NULL\n");
    428 
    429         if (material) material->Echo();
    430         else _printf_("material = NULL\n");
    431 
    432         if (matpar) matpar->Echo();
    433         else _printf_("matpar = NULL\n");
    434 
    435         _printf_("   parameters\n");
    436         if (parameters) parameters->Echo();
    437         else _printf_("parameters = NULL\n");
    438 
    439         _printf_("   inputs\n");
    440         if (inputs) inputs->Echo();
    441         else _printf_("inputs=NULL\n");
    442 }
    443 /*}}}*/
    444331/*FUNCTION Tria::FiniteElement{{{*/
    445332int Tria::FiniteElement(void){
    446333        return this->element_type;
    … …  
    921808        return this->NumberofNodes();
    922809}
    923810/*}}}*/
    924 /*FUNCTION Tria::GetNumberOfNodesPressure;{{{*/
    925 int Tria::GetNumberOfNodesPressure(void){
    926         return this->NumberofNodesPressure();
    927 }
    928 /*}}}*/
    929 /*FUNCTION Tria::GetNumberOfNodesVelocity;{{{*/
    930 int Tria::GetNumberOfNodesVelocity(void){
    931         return this->NumberofNodesVelocity();
    932 }
    933 /*}}}*/
    934811/*FUNCTION Tria::GetNumberOfVertices;{{{*/
    935812int Tria::GetNumberOfVertices(void){
    936813        return NUMVERTICES;
    … …  
    1009886        return y;
    1010887}
    1011888/*}}}*/
    1012 /*FUNCTION Tria::Id {{{*/
    1013 int    Tria::Id(){
    1014 
    1015         return id;
    1016 
    1017 }
    1018 /*}}}*/
    1019 /*FUNCTION Tria::Sid {{{*/
    1020 int    Tria::Sid(){
    1021 
    1022         return sid;
    1023 
    1024 }
    1025 /*}}}*/
    1026889/*FUNCTION Tria::InputDepthAverageAtBase {{{*/
    1027890void  Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){
    1028891
    … …  
    13501213        }
    13511214}
    13521215/*}}}*/
    1353 /*FUNCTION Tria::HasNodeOnBed {{{*/
    1354 bool Tria::HasNodeOnBed(){
    1355 
    1356         IssmDouble values[NUMVERTICES];
    1357         IssmDouble sum;
    1358 
    1359         /*Retrieve all inputs and parameters*/
    1360         GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
    1361         sum = values[0]+values[1]+values[2];
    1362 
    1363         return sum>0.;
    1364 }
    1365 /*}}}*/
    13661216/*FUNCTION Tria::HasEdgeOnSurface {{{*/
    13671217bool Tria::HasEdgeOnSurface(){
    13681218
    … …  
    13851235        }
    13861236}
    13871237/*}}}*/
    1388 /*FUNCTION Tria::HasNodeOnSurface {{{*/
    1389 bool Tria::HasNodeOnSurface(){
    1390 
    1391         IssmDouble values[NUMVERTICES];
    1392         IssmDouble sum;
    1393 
    1394         /*Retrieve all inputs and parameters*/
    1395         GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
    1396         sum = values[0]+values[1]+values[2];
    1397 
    1398         return sum>0.;
    1399 }
    1400 /*}}}*/
    14011238/*FUNCTION Tria::EdgeOnBedIndices{{{*/
    14021239void Tria::EdgeOnBedIndices(int* pindex1,int* pindex2){
    14031240
    … …  
    15411378        return new GaussTria(indices[0],indices[1],order);
    15421379}
    15431380/*}}}*/
    1544 /*FUNCTION Tria::IsIceInElement {{{*/
    1545 bool   Tria::IsIceInElement(){
    1546 
    1547         /*Get levelset*/
    1548         IssmDouble ls[NUMVERTICES];
    1549         GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
    1550 
    1551         /*If the level set on one of the nodes is <0, ice is present in this element*/
    1552         if(ls[0]<0. || ls[1]<0. || ls[2]<0.)
    1553                 return true;
    1554         else
    1555                 return false;
    1556 }
    1557 /*}}}*/
    15581381/*FUNCTION Tria::NodalFunctions{{{*/
    15591382void Tria::NodalFunctions(IssmDouble* basis, Gauss* gauss){
    15601383

  • ../trunk-jpl/src/c/classes/Elements/Tria.h

     
    3131
    3232        public:
    3333
    34                 int id;
    35                 int sid;
    36 
    3734                /*Tria constructors, destructors {{{*/
    3835                Tria(){};
    3936                Tria(int tria_id,int tria_sid,int i, IoModel* iomodel,int nummodels);
    4037                ~Tria();
    4138                /*}}}*/
    4239                /*Object virtual functions definitions:{{{ */
    43                 void    Echo();
    44                 void    DeepEcho();
    45                 int     Id();
    4640                int     ObjectEnum();
    4741                Object *copy();
    4842                /*}}}*/
    … …  
    6256                void        ComputeSurfaceNormalVelocity();
    6357                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
    6458                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
    65                 void        SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
    6659                void        Delta18oParameterization(void);
    6760                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
    68                 void        ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");};
    69                 void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
    7061                int         FiniteElement(void);
    7162                Element*    GetUpperElement(void){_error_("not implemented yet");};
    7263                Element*    GetLowerElement(void){_error_("not implemented yet");};
    … …  
    7667                IssmDouble  GetGroundedPortion(IssmDouble* xyz_list);
    7768                int         GetNodeIndex(Node* node);
    7869                int         GetNumberOfNodes(void);
    79                 int         GetNumberOfNodesPressure(void);
    80                 int         GetNumberOfNodesVelocity(void);
    8170                int         GetNumberOfVertices(void);
    82                 int         Sid();
    8371                bool        IsOnBed();
    8472                bool        IsOnSurface();
    8573                bool        HasEdgeOnBed();
    86                 bool        HasNodeOnBed();
    8774                bool        HasEdgeOnSurface();
    88                 bool        HasNodeOnSurface();
    8975                void        EdgeOnSurfaceIndices(int* pindex1,int* pindex);
    9076                void        EdgeOnBedIndices(int* pindex1,int* pindex);
    9177                int         EdgeOnBedIndex();
    … …  
    9379                bool        IsNodeOnShelfFromFlags(IssmDouble* flags);
    9480                int         NumberofNodesVelocity(void);
    9581                int         NumberofNodesPressure(void);
    96                 bool        IsIceInElement();
    9782                void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
    9883                void        GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum);
    9984                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
    … …  
    11095           Element*    SpawnBasalElement(void);
    11196                Element*    SpawnTopElement(void);
    11297                int         VelocityInterpolation();
    113                 IssmDouble  PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};
    11498                int         PressureInterpolation();
    11599                IssmDouble  SurfaceArea(void);
    116100                void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
    … …  
    191175                /*Tria specific routines:{{{*/
    192176                void           AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum);
    193177                void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
    194                 IssmDouble     EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
    195                 IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
    196178                IssmDouble     GetArea(void);
    197179                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
    198180                int            GetElementType(void);

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

     
    370370        this->inputs->Configure(parameters);
    371371}
    372372/*}}}*/
    373 /*FUNCTION Penta::DeepEcho{{{*/
    374 void Penta::DeepEcho(void){
    375 
    376         _printf_("Penta:\n");
    377         _printf_("   id: " << id << "\n");
    378         nodes[0]->DeepEcho();
    379         nodes[1]->DeepEcho();
    380         nodes[2]->DeepEcho();
    381         nodes[3]->DeepEcho();
    382         nodes[4]->DeepEcho();
    383         nodes[5]->DeepEcho();
    384         material->DeepEcho();
    385         matpar->DeepEcho();
    386         _printf_("   neighbor ids: " << verticalneighbors[0]->Id() << "-" << verticalneighbors[1]->Id() << "\n");
    387         _printf_("   parameters\n");
    388         parameters->DeepEcho();
    389         _printf_("   inputs\n");
    390         inputs->DeepEcho();
    391 }
    392 /*}}}*/
    393373/*FUNCTION Penta::Delta18oParameterization{{{*/
    394374void  Penta::Delta18oParameterization(void){
    395375        /*Are we on the base? If not, return*/
    … …  
    464444        delete gauss;
    465445}
    466446/*}}}*/
    467 /*FUNCTION Penta::Echo{{{*/
    468 
    469 void Penta::Echo(void){
    470         this->DeepEcho();
    471 }
    472 /*}}}*/
    473 /*FUNCTION Penta::ThermalToEnthalpy{{{*/
    474 void Penta::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){
    475         matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);
    476 }
    477 /*}}}*/
    478 /*FUNCTION Penta::EnthalpyToThermal{{{*/
    479 void Penta::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){
    480         matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure);
    481 }
    482 /*}}}*/
    483 /*FUNCTION Penta::EnthalpyDiffusionParameter{{{*/
    484 IssmDouble Penta::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){
    485         return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
    486 }
    487 /*}}}*/
    488 /*FUNCTION Penta::EnthalpyDiffusionParameterVolume{{{*/
    489 IssmDouble Penta::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){
    490         return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure);
    491 }
    492 /*}}}*/
    493447/*FUNCTION Penta::FiniteElement{{{*/
    494448int Penta::FiniteElement(void){
    495449        return this->element_type;
    … …  
    847801        return this->NumberofNodes();
    848802}
    849803/*}}}*/
    850 /*FUNCTION Penta::GetNumberOfNodesPressure;{{{*/
    851 int Penta::GetNumberOfNodesPressure(void){
    852         return this->NumberofNodesPressure();
    853 }
    854 /*}}}*/
    855 /*FUNCTION Penta::GetNumberOfNodesVelocity;{{{*/
    856 int Penta::GetNumberOfNodesVelocity(void){
    857         return this->NumberofNodesVelocity();
    858 }
    859 /*}}}*/
    860804/*FUNCTION Penta::GetNumberOfVertices;{{{*/
    861805int Penta::GetNumberOfVertices(void){
    862806        return NUMVERTICES;
    … …  
    11611105        *pxyz_zero= xyz_zero;
    11621106}
    11631107/*}}}*/
    1164 /*FUNCTION Penta::Sid {{{*/
    1165 int    Penta::Sid(){
    1166 
    1167         return sid;
    1168 
    1169 }
    1170 /*}}}*/
    1171 /*FUNCTION Penta::Id {{{*/
    1172 int    Penta::Id(void){
    1173         return id;
    1174 }
    1175 /*}}}*/
    11761108/*FUNCTION Penta::InputDepthAverageAtBase{{{*/
    11771109void  Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){
    11781110
    … …  
    16421574
    16431575}
    16441576/*}}}*/
    1645 /*FUNCTION Penta::IsIceInElement {{{*/
    1646 bool   Penta::IsIceInElement(){
    1647 
    1648         /*Get levelset*/
    1649         IssmDouble ls[NUMVERTICES];
    1650         GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
    1651 
    1652         /*If the level set one one of the nodes is <0, ice is present in this element*/
    1653         if(ls[0]<0. || ls[1]<0. || ls[2]<0.)
    1654                 return true;
    1655         else
    1656                 return false;
    1657 }
    1658 /*}}}*/
    16591577/*FUNCTION Penta::MinEdgeLength{{{*/
    16601578IssmDouble Penta::MinEdgeLength(IssmDouble* xyz_list){
    16611579        /*Return the minimum lenght of the nine egdes of the penta*/
    … …  
    18121730/*FUNCTION Penta::NormalBase {{{*/
    18131731void Penta::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){
    18141732
    1815         int i;
    18161733        IssmDouble v13[3],v23[3];
    18171734        IssmDouble normal[3];
    18181735        IssmDouble normal_norm;
    18191736
    1820         for (i=0;i<3;i++){
     1737        for(int i=0;i<3;i++){
    18211738                v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i];
    18221739                v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i];
    18231740        }
    … …  
    18251742        normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
    18261743        normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
    18271744        normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
    1828         normal_norm=sqrt( pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2) );
     1745        normal_norm=sqrt(normal[0]*normal[0]+ normal[1]*normal[1]+ normal[2]*normal[2]);
    18291746
    18301747        /*Bed normal is opposite to surface normal*/
    18311748        bed_normal[0]=-normal[0]/normal_norm;
    … …  
    19181835        delete gauss;
    19191836}
    19201837/*}}}*/
    1921 /*FUNCTION Penta::PureIceEnthalpy{{{*/
    1922 IssmDouble Penta::PureIceEnthalpy(IssmDouble pressure){
    1923 
    1924         return this->matpar->PureIceEnthalpy(pressure);
    1925 }
    1926 /*}}}*/
    19271838/*FUNCTION Penta::ReduceMatrices{{{*/
    19281839void Penta::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){
    19291840
    … …  
    20711982        this->element_type=element_type_in;
    20721983}
    20731984/*}}}*/
    2074 /*FUNCTION Penta::SetwiseNodeConnectivity{{{*/
    2075 void Penta::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){
    2076 
    2077         /*Intermediaries*/
    2078         const int numnodes = this->NumberofNodes();
    2079 
    2080         /*Output */
    2081         int d_nz = 0;
    2082         int o_nz = 0;
    2083 
    2084         /*Loop over all nodes*/
    2085         for(int i=0;i<numnodes;i++){
    2086 
    2087                 if(!flags[this->nodes[i]->Lid()]){
    2088 
    2089                         /*flag current node so that no other element processes it*/
    2090                         flags[this->nodes[i]->Lid()]=true;
    2091 
    2092                         int counter=0;
    2093                         while(flagsindices[counter]>=0) counter++;
    2094                         flagsindices[counter]=this->nodes[i]->Lid();
    2095 
    2096                         /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
    2097                         switch(set2_enum){
    2098                                 case FsetEnum:
    2099                                         if(nodes[i]->indexing.fsize){
    2100                                                 if(this->nodes[i]->IsClone())
    2101                                                  o_nz += 1;
    2102                                                 else
    2103                                                  d_nz += 1;
    2104                                         }
    2105                                         break;
    2106                                 case GsetEnum:
    2107                                         if(nodes[i]->indexing.gsize){
    2108                                                 if(this->nodes[i]->IsClone())
    2109                                                  o_nz += 1;
    2110                                                 else
    2111                                                  d_nz += 1;
    2112                                         }
    2113                                         break;
    2114                                 case SsetEnum:
    2115                                         if(nodes[i]->indexing.ssize){
    2116                                                 if(this->nodes[i]->IsClone())
    2117                                                  o_nz += 1;
    2118                                                 else
    2119                                                  d_nz += 1;
    2120                                         }
    2121                                         break;
    2122                                 default: _error_("not supported");
    2123                         }
    2124                 }
    2125         }
    2126 
    2127         /*Special case: 2d/3d coupling, the node of this element might be connected
    2128          *to the basal element*/
    2129         int analysis_type,approximation,numlayers;
    2130         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    2131         if(analysis_type==StressbalanceAnalysisEnum){
    2132                 inputs->GetInputValue(&approximation,ApproximationEnum);
    2133                 if(approximation==SSAHOApproximationEnum || approximation==SSAFSApproximationEnum){
    2134                         parameters->FindParam(&numlayers,MeshNumberoflayersEnum);
    2135                         o_nz += numlayers*3;
    2136                         d_nz += numlayers*3;
    2137                 }
    2138         }
    2139 
    2140         /*Assign output pointers: */
    2141         *pd_nz=d_nz;
    2142         *po_nz=o_nz;
    2143 }
    2144 /*}}}*/
    21451985/*FUNCTION Penta::SpawnTria {{{*/
    21461986Tria*  Penta::SpawnTria(int index1,int index2,int index3){
    21471987

  • ../trunk-jpl/src/c/classes/Elements/Penta.h

     
    3131
    3232        public:
    3333
    34                 int id;
    35                 int sid;
    36 
    3734                Penta      **verticalneighbors;           // 2 neighbors: first one under, second one above
    3835
    3936                /*Penta constructors and destructor: {{{*/
    … …  
    4340                /*}}}*/
    4441                /*Object virtual functions definitions: {{{*/
    4542                Object *copy();
    46                 void    DeepEcho();
    47                 void    Echo();
    4843                int     ObjectEnum();
    49                 int     Id();
    5044                /*}}}*/
    5145                /*Update virtual functions definitions: {{{*/
    5246                void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
    … …  
    6660                void   ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
    6761                int    FiniteElement(void);
    6862                void   SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
    69                 void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
    7063                void   Delta18oParameterization(void);
    71                 void   ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure);
    72                 void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
    7364                Penta* GetUpperPenta(void);
    7465                Penta* GetLowerPenta(void);
    7566                Penta* GetSurfacePenta(void);
    … …  
    8273                IssmDouble GetGroundedPortion(IssmDouble* xyz_list);
    8374                int    GetNodeIndex(Node* node);
    8475                int    GetNumberOfNodes(void);
    85                 int    GetNumberOfNodesPressure(void);
    86                 int    GetNumberOfNodesVelocity(void);
    8776                int    GetNumberOfVertices(void);
    8877                void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
    8978                IssmDouble GetXcoord(Gauss* gauss);
    … …  
    9382                void   GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
    9483                void   GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
    9584
    96                 int    Sid();
    9785                void   InputDepthAverageAtBase(int enum_type,int average_enum_type);
    9886                void   InputDuplicate(int original_enum,int new_enum);
    9987                void   InputScale(int enum_type,IssmDouble scale_factor);
    10088                int    NumberofNodesVelocity(void);
    10189                int    NumberofNodesPressure(void);
    10290                int    VelocityInterpolation();
    103                 IssmDouble PureIceEnthalpy(IssmDouble pressure);
    10491                int    PressureInterpolation();
    10592                bool   IsZeroLevelset(int levelset_enum);
    10693                bool   IsIcefront(void);
    … …  
    188175                void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
    189176                void             NormalTop(IssmDouble* bed_normal, IssmDouble* xyz_list);
    190177                ElementMatrix* CreateBasalMassMatrix(void);
    191                 IssmDouble     EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
    192                 IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);
    193178                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
    194179
    195180                void             GetVertexPidList(int* doflist);
    … …  
    207192                void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
    208193                void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
    209194                void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
    210                 bool           IsIceInElement(void);
    211195                Gauss*         NewGauss(void);
    212196                Gauss*         NewGauss(int order);
    213197                Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};

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

     
    6464        return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
    6565}
    6666/*}}}*/
    67 /*FUNCTION Seg::Echo{{{*/
    68 void Seg::Echo(void){
    69         _printf_("Seg:\n");
    70         _printf_("   id: " << id << "\n");
    71         if(nodes){
    72                 nodes[0]->Echo();
    73                 nodes[1]->Echo();
    74         }
    75         else _printf_("nodes = NULL\n");
    76 
    77         if (material) material->Echo();
    78         else _printf_("material = NULL\n");
    79 
    80         if (matpar) matpar->Echo();
    81         else _printf_("matpar = NULL\n");
    82 
    83         _printf_("   parameters\n");
    84         if (parameters) parameters->Echo();
    85         else _printf_("parameters = NULL\n");
    86 
    87         _printf_("   inputs\n");
    88         if (inputs) inputs->Echo();
    89         else _printf_("inputs=NULL\n");
    90 }
    91 /*}}}*/
    9267/*FUNCTION Seg::FiniteElement{{{*/
    9368int Seg::FiniteElement(void){
    9469        return this->element_type;
    9570}
    9671/*}}}*/
    97 /*FUNCTION Seg::DeepEcho{{{*/
    98 void Seg::DeepEcho(void){
    99 
    100         _printf_("Seg:\n");
    101         _printf_("   id: " << id << "\n");
    102         if(nodes){
    103                 nodes[0]->DeepEcho();
    104                 nodes[1]->DeepEcho();
    105         }
    106         else _printf_("nodes = NULL\n");
    107 
    108         if (material) material->DeepEcho();
    109         else _printf_("material = NULL\n");
    110 
    111         if (matpar) matpar->DeepEcho();
    112         else _printf_("matpar = NULL\n");
    113 
    114         _printf_("   parameters\n");
    115         if (parameters) parameters->DeepEcho();
    116         else _printf_("parameters = NULL\n");
    117 
    118         _printf_("   inputs\n");
    119         if (inputs) inputs->DeepEcho();
    120         else _printf_("inputs=NULL\n");
    121 
    122         return;
    123 }
    124 /*}}}*/
    12572/*FUNCTION Seg::ObjectEnum{{{*/
    12673int Seg::ObjectEnum(void){
    12774
    … …  
    12976
    13077}
    13178/*}}}*/
    132 /*FUNCTION Seg::Id {{{*/
    133 int    Seg::Id(){
    13479
    135         return id;
    136 
    137 }
    138 /*}}}*/
    139 
    14080void  Seg::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
    14181       
    14282        /* Intermediaries */
    … …  
    185125        *pxyz_list = xyz_list;
    186126
    187127}/*}}}*/
    188 /*FUNCTION Seg::IsIceInElement {{{*/
    189 bool   Seg::IsIceInElement(){
    190 
    191         /*Get levelset*/
    192         IssmDouble ls[NUMVERTICES];
    193         GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
    194 
    195         /*If the level set on one of the nodes is <0, ice is present in this element*/
    196         if(ls[0]<0. || ls[1]<0.)
    197          return true;
    198         else
    199          return false;
    200 }
    201 /*}}}*/
    202128bool Seg::IsIcefront(void){/*{{{*/
    203129
    204130        bool isicefront;

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

     
    5151}
    5252/*}}}*/
    5353
    54 /*FUNCTION Tetra::Echo{{{*/
    55 void Tetra::Echo(void){
    56         _printf_("Tetra:\n");
    57         _printf_("   id: " << id << "\n");
    58         if(nodes){
    59                 nodes[0]->Echo();
    60                 nodes[1]->Echo();
    61                 nodes[2]->Echo();
    62                 nodes[3]->Echo();
    63         }
    64         else _printf_("nodes = NULL\n");
    65 
    66         if (material) material->Echo();
    67         else _printf_("material = NULL\n");
    68 
    69         if (matpar) matpar->Echo();
    70         else _printf_("matpar = NULL\n");
    71 
    72         _printf_("   parameters\n");
    73         if (parameters) parameters->Echo();
    74         else _printf_("parameters = NULL\n");
    75 
    76         _printf_("   inputs\n");
    77         if (inputs) inputs->Echo();
    78         else _printf_("inputs=NULL\n");
    79 }
    80 /*}}}*/
    8154/*FUNCTION Tetra::FiniteElement{{{*/
    8255int Tetra::FiniteElement(void){
    8356        return this->element_type;
    84 }
    85 /*}}}*/
    86 /*FUNCTION Tetra::DeepEcho{{{*/
    87 void Tetra::DeepEcho(void){
    88 
    89         _printf_("Tetra:\n");
    90         _printf_("   id: " << id << "\n");
    91         if(nodes){
    92                 nodes[0]->DeepEcho();
    93                 nodes[1]->DeepEcho();
    94                 nodes[2]->DeepEcho();
    95                 nodes[3]->DeepEcho();
    96         }
    97         else _printf_("nodes = NULL\n");
    98 
    99         if (material) material->DeepEcho();
    100         else _printf_("material = NULL\n");
    101 
    102         if (matpar) matpar->DeepEcho();
    103         else _printf_("matpar = NULL\n");
    104 
    105         _printf_("   parameters\n");
    106         if (parameters) parameters->DeepEcho();
    107         else _printf_("parameters = NULL\n");
    108 
    109         _printf_("   inputs\n");
    110         if (inputs) inputs->DeepEcho();
    111         else _printf_("inputs=NULL\n");
    112 
    113         return;
    114 }
    115 /*}}}*/
     57} /*}}}*/
    11658/*FUNCTION Tetra::ObjectEnum{{{*/
    11759int Tetra::ObjectEnum(void){
    11860
    11961        return TetraEnum;
    12062
    121 }
    122 /*}}}*/
    123 /*FUNCTION Tetra::Id {{{*/
    124 int Tetra::Id(){
    125         return id;
    126 }
    127 /*}}}*/
     63}/*}}}*/
    12864
    12965/*FUNCTION Tetra::AddInput{{{*/
    13066void  Tetra::AddInput(int input_enum,IssmDouble* values, int interpolation_enum){
    … …  
    239175        return NUMVERTICES;
    240176}
    241177/*}}}*/
    242 /*FUNCTION Tetra::GetNumberOfNodesPressure         THIS ONE (and corresponding TetraRef function){{{*/
    243 int Tetra::GetNumberOfNodesPressure(void){
    244         return this->NumberofNodesPressure();
    245 }
    246 /*}}}*/
    247 /*FUNCTION Tetra::GetNumberOfNodesVelocity;{{{*/
    248 int Tetra::GetNumberOfNodesVelocity(void){
    249         return this->NumberofNodesVelocity();
    250 }
    251 /*}}}*/
    252178/*FUNCTION Tetra::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){{{*/
    253179void Tetra::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){
    254180
    … …  
    340266        }
    341267}
    342268/*}}}*/
    343 /*FUNCTION Tetra::HasNodeOnBed           THIS ONE{{{*/
    344 bool Tetra::HasNodeOnBed(){
    345 
    346         IssmDouble values[NUMVERTICES];
    347         IssmDouble sum;
    348 
    349         /*Retrieve all inputs and parameters*/
    350         GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
    351         sum = values[0]+values[1]+values[2]+values[3];
    352 
    353         return sum>0.;
    354 }
    355 /*}}}*/
    356269/*FUNCTION Tetra::InputUpdateFromIoModel {{{*/
    357270void Tetra::InputUpdateFromIoModel(int index,IoModel* iomodel){
    358271
    … …  
    484397        xDelete<int>(doflist);
    485398}
    486399/*}}}*/
    487 /*FUNCTION Tetra::IsIceInElement    THIS ONE{{{*/
    488 bool   Tetra::IsIceInElement(){
    489 
    490         /*Get levelset*/
    491         IssmDouble ls[NUMVERTICES];
    492         GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
    493 
    494         /*If the level set on one of the nodes is <0, ice is present in this element*/
    495         for(int i=0;i<NUMVERTICES;i++){
    496                 if(ls[i]<0.) return true;
    497         }
    498 
    499         return false;
    500 }
    501 /*}}}*/
    502400/*FUNCTION Tetra::IsOnBed {{{*/
    503401bool Tetra::IsOnBed(){
    504402        return HasFaceOnBed();
    … …  
    650548        for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
    651549}
    652550/*}}}*/
     551/*FUNCTION Tetra::NormalBase (THIS ONE){{{*/
     552void Tetra::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){
     553
     554        IssmDouble v13[3],v23[3];
     555        IssmDouble normal[3];
     556        IssmDouble normal_norm;
     557
     558        for(int i=0;i<3;i++){
     559                v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i];
     560                v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i];
     561        }
     562
     563        normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
     564        normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
     565        normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
     566        normal_norm=sqrt(normal[0]*normal[0]+ normal[1]*normal[1]+ normal[2]*normal[2]);
     567
     568        /*Bed normal is opposite to surface normal*/
     569        bed_normal[0]=-normal[0]/normal_norm;
     570        bed_normal[1]=-normal[1]/normal_norm;
     571        bed_normal[2]=-normal[2]/normal_norm;
     572}
     573/*}}}*/
     574/*FUNCTION Tetra::NormalTop (THIS ONE){{{*/
     575void Tetra::NormalTop(IssmDouble* top_normal,IssmDouble* xyz_list){
     576
     577        IssmDouble v13[3],v23[3];
     578        IssmDouble normal[3];
     579        IssmDouble normal_norm;
     580
     581        for(int i=0;i<3;i++){
     582                v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i];
     583                v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i];
     584        }
     585
     586        normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
     587        normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
     588        normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
     589        normal_norm=sqrt(normal[0]*normal[0]+ normal[1]*normal[1]+ normal[2]*normal[2]);
     590
     591        top_normal[0]=normal[0]/normal_norm;
     592        top_normal[1]=normal[1]/normal_norm;
     593        top_normal[2]=normal[2]/normal_norm;
     594}
     595/*}}}*/
    653596/*FUNCTION Tetra::NumberofNodesPressure{{{*/
    654597int Tetra::NumberofNodesPressure(void){
    655598        return TetraRef::NumberofNodesPressure();
    … …  
    665608
    666609        if(pe){
    667610                if(this->element_type==MINIcondensedEnum){
    668                         _error_("Not implemented");
     611                        int indices[3]={12,13,14};
     612                        pe->StaticCondensation(Ke,3,&indices[0]);
    669613                }
    670614                else if(this->element_type==P1bubblecondensedEnum){
    671615                        int size   = nodes[4]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
    … …  
    680624
    681625        if(Ke){
    682626                if(this->element_type==MINIcondensedEnum){
    683                         _error_("Not implemented");
     627                        int indices[3]={12,13,14};
     628                        Ke->StaticCondensation(3,&indices[0]);
    684629                }
    685630                else if(this->element_type==P1bubblecondensedEnum){
    686631                        int size   = nodes[4]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
    … …  
    770715
    771716}
    772717/*}}}*/
    773 /*FUNCTION Tetra::SetwiseNodeConnectivity   THIS ONE (take from Penta.cpp){{{*/
    774 void Tetra::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){
    775 
    776         /*Intermediaries*/
    777         const int numnodes = this->NumberofNodes();
    778 
    779         /*Output */
    780         int d_nz = 0;
    781         int o_nz = 0;
    782 
    783         /*Loop over all nodes*/
    784         for(int i=0;i<numnodes;i++){
    785 
    786                 if(!flags[this->nodes[i]->Lid()]){
    787 
    788                         /*flag current node so that no other element processes it*/
    789                         flags[this->nodes[i]->Lid()]=true;
    790 
    791                         int counter=0;
    792                         while(flagsindices[counter]>=0) counter++;
    793                         flagsindices[counter]=this->nodes[i]->Lid();
    794 
    795                         /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
    796                         switch(set2_enum){
    797                                 case FsetEnum:
    798                                         if(nodes[i]->indexing.fsize){
    799                                                 if(this->nodes[i]->IsClone())
    800                                                  o_nz += 1;
    801                                                 else
    802                                                  d_nz += 1;
    803                                         }
    804                                         break;
    805                                 case GsetEnum:
    806                                         if(nodes[i]->indexing.gsize){
    807                                                 if(this->nodes[i]->IsClone())
    808                                                  o_nz += 1;
    809                                                 else
    810                                                  d_nz += 1;
    811                                         }
    812                                         break;
    813                                 case SsetEnum:
    814                                         if(nodes[i]->indexing.ssize){
    815                                                 if(this->nodes[i]->IsClone())
    816                                                  o_nz += 1;
    817                                                 else
    818                                                  d_nz += 1;
    819                                         }
    820                                         break;
    821                                 default: _error_("not supported");
    822                         }
    823                 }
    824         }
    825 
    826         /*Assign output pointers: */
    827         *pd_nz=d_nz;
    828         *po_nz=o_nz;
    829 }
    830 /*}}}*/
    831 /*FUNCTION Tetra::Sid  THIS ONE{{{*/
    832 int    Tetra::Sid(){
    833 
    834         return sid;
    835 
    836 }
    837 /*}}}*/
    838718/*FUNCTION Tetra::SpawnBasalElement{{{*/
    839719Element*  Tetra::SpawnBasalElement(void){
    840720

  • ../trunk-jpl/src/c/classes/Elements/Seg.h

     
    2929
    3030        public:
    3131
    32                 int id;
    33                 int sid;
    34 
    3532                /*Seg constructors, destructors {{{*/
    3633                Seg(){};
    3734                Seg(int seg_id,int seg_sid,int i, IoModel* iomodel,int nummodels);
    3835                ~Seg();
    3936                /*}}}*/
    4037                /*Object virtual functions definitions:{{{ */
    41                 void    Echo();
    42                 void    DeepEcho();
    43                 int     Id();
    4438                int     ObjectEnum();
    4539                Object *copy();
    4640                /*}}}*/
    … …  
    6357                void        ComputeStressTensor(){_error_("not implemented yet");};
    6458                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
    6559                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
    66                 void        SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){_error_("not implemented yet");};
    6760                void        Delta18oParameterization(void){_error_("not implemented yet");};
    6861                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
    69                 void        ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");};
    70                 void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
    71                 IssmDouble  EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
    72                 IssmDouble  EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
    7362                int         FiniteElement(void);
    7463                Element*    GetUpperElement(void){_error_("not implemented yet");};
    7564                Element*    GetLowerElement(void){_error_("not implemented yet");};
    … …  
    7766                Element*    GetBasalElement(void){_error_("not implemented yet");};
    7867                int         GetNodeIndex(Node* node){_error_("not implemented yet");};
    7968                int         GetNumberOfNodes(void);
    80                 int         GetNumberOfNodesVelocity(void){_error_("not implemented yet");};
    81                 int         GetNumberOfNodesPressure(void){_error_("not implemented yet");};
    8269                int         GetNumberOfVertices(void);
    8370                void        GetVerticesCoordinates(IssmDouble** pxyz_list);
    8471                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list){_error_("not implemented yet");};
    8572                void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");};
    86                 int         Sid(){_error_("not implemented yet");};
    8773                bool        IsOnBed(){_error_("not implemented yet");};
    8874                bool        IsOnSurface(){_error_("not implemented yet");};
    8975                bool        IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");};
    … …  
    10187                void        NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    10288                void        NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    10389                void        NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    104                 bool        IsIceInElement();
    10590                void        NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
    10691                void        NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
    10792                void        NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
    … …  
    11095           Element*    SpawnBasalElement(void){_error_("not implemented yet");};
    11196                Element*    SpawnTopElement(void){_error_("not implemented yet");};
    11297                IssmDouble  StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
    113                 IssmDouble  PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};
    11498                int         PressureInterpolation(void){_error_("not implemented yet");};
    11599                void        ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){_error_("not implemented yet");};
    116100                void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};

  • ../trunk-jpl/src/c/classes/Elements/Tetra.h

     
    2929
    3030        public:
    3131
    32                 int id;
    33                 int sid;
    34 
    3532                /*Tetra constructors, destructors {{{*/
    3633                Tetra(){};
    3734                Tetra(int seg_id,int seg_sid,int i, IoModel* iomodel,int nummodels);
    3835                ~Tetra();
    3936                /*}}}*/
    4037                /*Object virtual functions definitions:{{{ */
    41                 void    Echo();
    42                 void    DeepEcho();
    43                 int     Id();
    4438                int     ObjectEnum();
    4539                Object *copy();
    4640                /*}}}*/
    … …  
    6357                void        ComputeStressTensor(){_error_("not implemented yet");};
    6458                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
    6559                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
    66                 void        SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
    6760                void        Delta18oParameterization(void){_error_("not implemented yet");};
    6861                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
    69                 void        ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");};
    70                 void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
    71                 IssmDouble  EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
    72                 IssmDouble  EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
    7362                void        FaceOnFrontIndices(int* pindex1,int* pindex2,int* pindex3);
    7463                void        FaceOnBedIndices(int* pindex1,int* pindex2,int* pindex3);
    7564                void        FaceOnSurfaceIndices(int* pindex1,int* pindex2,int* pindex3);
    … …  
    8069                Element*    GetBasalElement(void){_error_("not implemented yet");};
    8170                int         GetNodeIndex(Node* node){_error_("not implemented yet");};
    8271                int         GetNumberOfNodes(void);
    83                 int         GetNumberOfNodesVelocity(void);
    84                 int         GetNumberOfNodesPressure(void);
    8572                int         GetNumberOfVertices(void);
    8673                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
    8774                void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
    8875                bool        HasFaceOnBed();
    8976                bool        HasFaceOnSurface();
    90                 bool        HasNodeOnBed();
    91                 int         Sid();
    9277                bool        IsOnBed();
    9378                bool        IsOnSurface();
    9479                bool        IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");};
    … …  
    10691                void        NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    10792                void        NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    10893                void        NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    109                 bool        IsIceInElement();
    11094                void        NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
    111                 void        NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
    112                 void        NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
     95                void        NormalTop(IssmDouble* normal,IssmDouble* xyz_list);
     96                void        NormalBase(IssmDouble* normal,IssmDouble* xyz_list);
    11397                int         NumberofNodesVelocity(void);
    11498                int         NumberofNodesPressure(void);
    11599           Element*    SpawnBasalElement(void);
    116100                Element*    SpawnTopElement(void);
    117101                Tria*       SpawnTria(int index1,int index2,int index3);
    118102                IssmDouble  StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
    119                 IssmDouble  PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};
    120103                int         PressureInterpolation(void);
    121104                void        ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){_error_("not implemented yet");};
    122105                void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    … …  
    227210                int    UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){_error_("not implemented yet");};
    228211                /*}}}*/
    229212};
    230 #endif  /* _SEG_H */
     213#endif  /* _TETRA_H_*/

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

     
    1616
    1717/*Constructors/destructor/copy*/
    1818Element::Element(){/*{{{*/
     19        this->id  = -1;
     20        this->sid = -1;
    1921        this->inputs     = NULL;
    2022        this->nodes      = NULL;
    2123        this->vertices   = NULL;
    … …  
    129131void Element::DeleteMaterials(void){/*{{{*/
    130132        delete this->material;
    131133}/*}}}*/
     134void Element::DeepEcho(void){/*{{{*/
     135
     136        _printf_(EnumToStringx(this->ObjectEnum())<<" element:\n");
     137        _printf_("   id : "<<this->id <<"\n");
     138        _printf_("   sid: "<<this->sid<<"\n");
     139        if(vertices){
     140                int numvertices = this->GetNumberOfVertices();
     141                for(int i=0;i<numvertices;i++) vertices[i]->Echo();
     142        }
     143        else _printf_("vertices = NULL\n");
     144
     145        if(nodes){
     146                int numnodes = this->GetNumberOfNodes();
     147                for(int i=0;i<numnodes;i++) nodes[i]->DeepEcho();
     148        }
     149        else _printf_("nodes = NULL\n");
     150
     151        if (material) material->DeepEcho();
     152        else _printf_("material = NULL\n");
     153
     154        if (matpar) matpar->DeepEcho();
     155        else _printf_("matpar = NULL\n");
     156
     157        _printf_("   parameters\n");
     158        if (parameters) parameters->DeepEcho();
     159        else _printf_("parameters = NULL\n");
     160
     161        _printf_("   inputs\n");
     162        if (inputs) inputs->DeepEcho();
     163        else _printf_("inputs=NULL\n");
     164
     165        return;
     166}
     167/*}}}*/
     168void Element::Echo(void){/*{{{*/
     169        _printf_(EnumToStringx(this->ObjectEnum())<<" element:\n");
     170        _printf_("   id : "<<this->id <<"\n");
     171        _printf_("   sid: "<<this->sid<<"\n");
     172        if(vertices){
     173                int numvertices = this->GetNumberOfVertices();
     174                for(int i=0;i<numvertices;i++) vertices[i]->Echo();
     175        }
     176        else _printf_("vertices = NULL\n");
     177
     178        if(nodes){
     179                int numnodes = this->GetNumberOfNodes();
     180                for(int i=0;i<numnodes;i++) nodes[i]->Echo();
     181        }
     182        else _printf_("nodes = NULL\n");
     183
     184        if (material) material->Echo();
     185        else _printf_("material = NULL\n");
     186
     187        if (matpar) matpar->Echo();
     188        else _printf_("matpar = NULL\n");
     189
     190        _printf_("   parameters\n");
     191        if (parameters) parameters->Echo();
     192        else _printf_("parameters = NULL\n");
     193
     194        _printf_("   inputs\n");
     195        if (inputs) inputs->Echo();
     196        else _printf_("inputs=NULL\n");
     197}
     198/*}}}*/
    132199IssmDouble Element::Divergence(void){/*{{{*/
    133200        /*Compute element divergence*/
    134201
    … …  
    161228        delete gauss;
    162229        return divergence;
    163230}/*}}}*/
     231void Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
     232        matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);
     233}/*}}}*/
     234void Element::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
     235        matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure);
     236}/*}}}*/
     237IssmDouble Element::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
     238        return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
     239}/*}}}*/
     240IssmDouble Element::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/
     241        return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure);
     242}/*}}}*/
    164243void Element::FindParam(bool* pvalue,int paramenum){/*{{{*/
    165244        this->parameters->FindParam(pvalue,paramenum);
    166245}/*}}}*/
    … …  
    199278void Element::GetDofListVelocity(int** pdoflist,int setenum){/*{{{*/
    200279
    201280        /*Fetch number of nodes and dof for this finite element*/
    202         int numnodes = this->GetNumberOfNodesVelocity();
     281        int numnodes = this->NumberofNodesVelocity();
    203282
    204283        /*First, figure out size of doflist and create it: */
    205284        int numberofdofs=0;
    … …  
    222301void Element::GetDofListPressure(int** pdoflist,int setenum){/*{{{*/
    223302
    224303        /*Fetch number of nodes and dof for this finite element*/
    225         int vnumnodes = this->GetNumberOfNodesVelocity();
    226         int pnumnodes = this->GetNumberOfNodesPressure();
     304        int vnumnodes = this->NumberofNodesVelocity();
     305        int pnumnodes = this->NumberofNodesPressure();
    227306
    228307        /*First, figure out size of doflist and create it: */
    229308        int numberofdofs=0;
    … …  
    387466
    388467        _assert_(pvalue);
    389468
    390         int    numnodes = this->GetNumberOfNodesVelocity();
     469        int    numnodes = this->NumberofNodesVelocity();
    391470        Input *input    = this->GetInput(enumtype);
    392471        if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
    393472
    … …  
    469548        for(int i=0;i<numvertices;i++) connectivity[i]=this->vertices[i]->Connectivity();
    470549}
    471550/*}}}*/
     551bool Element::HasNodeOnBed(){/*{{{*/
     552        return (this->inputs->Max(MeshVertexonbedEnum)>0.);
     553}/*}}}*/
     554bool Element::HasNodeOnSurface(){/*{{{*/
     555        return (this->inputs->Max(MeshVertexonsurfaceEnum)>0.);
     556}/*}}}*/
     557int  Element::Id(){/*{{{*/
     558
     559        return this->id;
     560
     561}
     562/*}}}*/
    472563void Element::InputChangeName(int original_enum,int new_enum){/*{{{*/
    473564        this->inputs->ChangeEnum(original_enum,new_enum);
    474565}
    … …  
    586677
    587678        return shelf;
    588679}/*}}}*/
     680bool Element::IsIceInElement(){/*{{{*/
     681        return (this->inputs->Min(MaskIceLevelsetEnum)<0.);
     682}
     683/*}}}*/
    589684bool Element::IsInput(int name){/*{{{*/
    590685        if (
    591686                                name==ThicknessEnum ||
    … …  
    674769        return new ElementMatrix(nodes,number_nodes,this->parameters,approximation_enum);
    675770}
    676771/*}}}*/
     772IssmDouble Element::PureIceEnthalpy(IssmDouble pressure){/*{{{*/
     773        return this->matpar->PureIceEnthalpy(pressure);
     774}/*}}}*/
    677775void Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int output_enum){/*{{{*/
    678776
    679777        Input* input=this->inputs->GetInput(output_enum);
    … …  
    772870        input->ResultToPatch(values,nodesperelement,this->Sid());
    773871
    774872} /*}}}*/
     873void Element::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/
     874
     875        /*Intermediaries*/
     876        const int numnodes = this->GetNumberOfNodes();
     877
     878        /*Output */
     879        int d_nz = 0;
     880        int o_nz = 0;
     881
     882        /*Loop over all nodes*/
     883        for(int i=0;i<numnodes;i++){
     884
     885                if(!flags[this->nodes[i]->Lid()]){
     886
     887                        /*flag current node so that no other element processes it*/
     888                        flags[this->nodes[i]->Lid()]=true;
     889
     890                        int counter=0;
     891                        while(flagsindices[counter]>=0) counter++;
     892                        flagsindices[counter]=this->nodes[i]->Lid();
     893
     894                        /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
     895                        switch(set2_enum){
     896                                case FsetEnum:
     897                                        if(nodes[i]->indexing.fsize){
     898                                                if(this->nodes[i]->IsClone())
     899                                                 o_nz += 1;
     900                                                else
     901                                                 d_nz += 1;
     902                                        }
     903                                        break;
     904                                case GsetEnum:
     905                                        if(nodes[i]->indexing.gsize){
     906                                                if(this->nodes[i]->IsClone())
     907                                                 o_nz += 1;
     908                                                else
     909                                                 d_nz += 1;
     910                                        }
     911                                        break;
     912                                case SsetEnum:
     913                                        if(nodes[i]->indexing.ssize){
     914                                                if(this->nodes[i]->IsClone())
     915                                                 o_nz += 1;
     916                                                else
     917                                                 d_nz += 1;
     918                                        }
     919                                        break;
     920                                default: _error_("not supported");
     921                        }
     922                }
     923        }
     924
     925        /*Special case: 2d/3d coupling, the node of this element might be connected
     926         *to the basal element*/
     927        int analysis_type,approximation,numlayers;
     928        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     929        if(analysis_type==StressbalanceAnalysisEnum){
     930                inputs->GetInputValue(&approximation,ApproximationEnum);
     931                if(approximation==SSAHOApproximationEnum || approximation==SSAFSApproximationEnum){
     932                        parameters->FindParam(&numlayers,MeshNumberoflayersEnum);
     933                        o_nz += numlayers*3;
     934                        d_nz += numlayers*3;
     935                }
     936        }
     937
     938        /*Assign output pointers: */
     939        *pd_nz=d_nz;
     940        *po_nz=o_nz;
     941}
     942/*}}}*/
     943int  Element::Sid(){/*{{{*/
     944
     945        return this->sid;
     946
     947}
     948/*}}}*/
    775949IssmDouble Element::TMeltingPoint(IssmDouble pressure){/*{{{*/
    776950        _assert_(matpar);
    777951        return this->matpar->TMeltingPoint(pressure);
Note: See TracBrowser for help on using the repository browser.