Changeset 16539


Ignore:
Timestamp:
10/24/13 10:12:44 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: moving model processor stuff to analyses

Location:
issm/trunk-jpl/src/c
Files:
21 deleted
58 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r16534 r16539  
    376376#}}}
    377377#Masstransport sources  {{{
    378 masstransport_sources = ./modules/ModelProcessorx/Masstransport/UpdateElementsMasstransport.cpp\
    379                                                                 ./modules/ModelProcessorx/Masstransport/CreateNodesMasstransport.cpp\
    380                                                                 ./modules/ModelProcessorx/Masstransport/CreateConstraintsMasstransport.cpp\
    381                                                                 ./modules/ModelProcessorx/Masstransport/CreateLoadsMasstransport.cpp\
    382                                                                 ./modules/ModelProcessorx/FreeSurfaceTop/UpdateElementsFreeSurfaceTop.cpp\
    383                                                                 ./modules/ModelProcessorx/FreeSurfaceTop/CreateNodesFreeSurfaceTop.cpp\
    384                                                                 ./modules/ModelProcessorx/FreeSurfaceTop/CreateConstraintsFreeSurfaceTop.cpp\
    385                                                                 ./modules/ModelProcessorx/FreeSurfaceTop/CreateLoadsFreeSurfaceTop.cpp\
    386                                                                 ./modules/ModelProcessorx/FreeSurfaceBase/UpdateElementsFreeSurfaceBase.cpp\
    387                                                                 ./modules/ModelProcessorx/FreeSurfaceBase/CreateNodesFreeSurfaceBase.cpp\
    388                                                                 ./modules/ModelProcessorx/FreeSurfaceBase/CreateConstraintsFreeSurfaceBase.cpp\
    389                                                                 ./modules/ModelProcessorx/FreeSurfaceBase/CreateLoadsFreeSurfaceBase.cpp\
    390                                                                 ./modules/ModelProcessorx/ExtrudeFromBase/UpdateElementsExtrudeFromBase.cpp\
    391                                                                 ./modules/ModelProcessorx/ExtrudeFromBase/CreateNodesExtrudeFromBase.cpp \
    392                                                                 ./modules/ModelProcessorx/ExtrudeFromBase/CreateConstraintsExtrudeFromBase.cpp\
    393                                                                 ./modules/ModelProcessorx/ExtrudeFromBase/CreateLoadsExtrudeFromBase.cpp\
    394                                                                 ./modules/ModelProcessorx/ExtrudeFromTop/UpdateElementsExtrudeFromTop.cpp\
    395                                                                 ./modules/ModelProcessorx/ExtrudeFromTop/CreateNodesExtrudeFromTop.cpp \
    396                                                                 ./modules/ModelProcessorx/ExtrudeFromTop/CreateConstraintsExtrudeFromTop.cpp\
    397                                                                 ./modules/ModelProcessorx/ExtrudeFromTop/CreateLoadsExtrudeFromTop.cpp\
    398                                                                 ./analyses/ExtrudeFromBaseAnalysis.h\
     378masstransport_sources = ./analyses/ExtrudeFromBaseAnalysis.h\
    399379                                                                ./analyses/ExtrudeFromBaseAnalysis.cpp\
    400380                                                                ./analyses/ExtrudeFromTopAnalysis.h\
    … …  
    412392#}}}
    413393#Thermal sources  {{{
    414 thermal_sources = ./modules/ModelProcessorx/Thermal/UpdateElementsThermal.cpp\
    415                                            ./modules/ModelProcessorx/Thermal/CreateNodesThermal.cpp\
    416                                            ./modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp\
    417                                            ./modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp\
    418                                            ./modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp\
    419                                            ./modules/ModelProcessorx/Enthalpy/CreateNodesEnthalpy.cpp\
    420                                            ./modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp\
    421                                            ./modules/ModelProcessorx/Enthalpy/CreateLoadsEnthalpy.cpp\
    422                                            ./modules/ModelProcessorx/Melting/UpdateElementsMelting.cpp\
    423                                            ./modules/ModelProcessorx/Melting/CreateNodesMelting.cpp\
    424                                            ./modules/ModelProcessorx/Melting/CreateConstraintsMelting.cpp\
    425                                            ./modules/ModelProcessorx/Melting/CreateLoadsMelting.cpp\
    426                                            ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.h\
     394thermal_sources = ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.h\
    427395                                           ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.cpp\
    428396                                                ./analyses/ThermalAnalysis.h\
    … …  
    493461#}}}
    494462#Hydrology sources  {{{
    495 hydrology_sources  = ./modules/ModelProcessorx/HydrologyShreve/UpdateElementsHydrologyShreve.cpp\
    496                                               ./modules/ModelProcessorx/HydrologyShreve/CreateNodesHydrologyShreve.cpp\
    497                                               ./modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp\
    498                                               ./modules/ModelProcessorx/HydrologyShreve/CreateLoadsHydrologyShreve.cpp \
    499                                                         ./modules/ModelProcessorx/HydrologyShreve/CreateParametersHydrologyShreve.cpp \
    500                                                         ./modules/ModelProcessorx/HydrologyDCInefficient/UpdateElementsHydrologyDCInefficient.cpp\
    501                                                         ./modules/ModelProcessorx/HydrologyDCInefficient/CreateNodesHydrologyDCInefficient.cpp\
    502                                                         ./modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp\
    503                                                         ./modules/ModelProcessorx/HydrologyDCInefficient/CreateLoadsHydrologyDCInefficient.cpp \
    504                                                         ./modules/ModelProcessorx/HydrologyDCInefficient/CreateParametersHydrologyDCInefficient.cpp \
    505                                                         ./modules/ModelProcessorx/HydrologyDCEfficient/UpdateElementsHydrologyDCEfficient.cpp\
    506                                                         ./modules/ModelProcessorx/HydrologyDCEfficient/CreateNodesHydrologyDCEfficient.cpp\
    507                                                         ./modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp\
    508                                                         ./modules/ModelProcessorx/HydrologyDCEfficient/CreateLoadsHydrologyDCEfficient.cpp \
    509                                                         ./modules/ModelProcessorx/HydrologyDCEfficient/CreateParametersHydrologyDCEfficient.cpp \
    510                                                         ./analyses/HydrologyDCEfficientAnalysis.h\
     463hydrology_sources  = ./analyses/HydrologyDCEfficientAnalysis.h\
    511464                                                        ./analyses/HydrologyDCEfficientAnalysis.cpp\
    512465                                                        ./analyses/HydrologyDCInefficientAnalysis.h\
    … …  
    518471#}}}
    519472#Stressbalance sources  {{{
    520 stressbalance_sources = ./modules/ModelProcessorx/Stressbalance/UpdateElementsStressbalance.cpp\
    521                                               ./modules/ModelProcessorx/Stressbalance/CreateNodesStressbalance.cpp \
    522                                               ./modules/ModelProcessorx/Stressbalance/CreateConstraintsStressbalance.cpp \
    523                                               ./modules/ModelProcessorx/Stressbalance/CreateLoadsStressbalance.cpp\
    524                                                         ./modules/ModelProcessorx/Stressbalance/CreateParametersStressbalance.cpp\
    525                                               ./modules/ModelProcessorx/StressbalanceVertical/UpdateElementsStressbalanceVertical.cpp\
    526                                               ./modules/ModelProcessorx/StressbalanceVertical/CreateNodesStressbalanceVertical.cpp \
    527                                               ./modules/ModelProcessorx/StressbalanceVertical/CreateConstraintsStressbalanceVertical.cpp \
    528                                               ./modules/ModelProcessorx/StressbalanceVertical/CreateLoadsStressbalanceVertical.cpp\
    529                                               ./modules/ModelProcessorx/StressbalanceSIA/UpdateElementsStressbalanceSIA.cpp\
    530                                               ./modules/ModelProcessorx/StressbalanceSIA/CreateNodesStressbalanceSIA.cpp \
    531                                               ./modules/ModelProcessorx/StressbalanceSIA/CreateConstraintsStressbalanceSIA.cpp \
    532                                                         ./modules/ModelProcessorx/StressbalanceSIA/CreateLoadsStressbalanceSIA.cpp \
    533                                                         ./analyses/StressbalanceAnalysis.h\
     473stressbalance_sources = ./analyses/StressbalanceAnalysis.h\
    534474                                                        ./analyses/StressbalanceAnalysis.cpp\
    535475                                                        ./analyses/StressbalanceSIAAnalysis.h\
    … …  
    541481#}}}
    542482#Balanced sources  {{{
    543 balanced_sources = ./modules/ModelProcessorx/Balancethickness/UpdateElementsBalancethickness.cpp\
    544                                             ./modules/ModelProcessorx/Balancethickness/CreateNodesBalancethickness.cpp\
    545                                             ./modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp\
    546                                                  ./modules/ModelProcessorx/Balancethickness/CreateLoadsBalancethickness.cpp\
    547                                                  ./modules/ModelProcessorx/Balancevelocity/UpdateElementsBalancevelocity.cpp\
    548                                                  ./modules/ModelProcessorx/Balancevelocity/CreateNodesBalancevelocity.cpp\
    549                                                  ./modules/ModelProcessorx/Balancevelocity/CreateConstraintsBalancevelocity.cpp\
    550                                                  ./modules/ModelProcessorx/Balancevelocity/CreateLoadsBalancevelocity.cpp\
    551                                                  ./analyses/BalancevelocityAnalysis.h\
     483balanced_sources = ./analyses/BalancevelocityAnalysis.h\
    552484                                                 ./analyses/BalancevelocityAnalysis.cpp\
    553485                                                 ./analyses/SmoothedSurfaceSlopeXAnalysis.h\
    … …  
    564496#}}}
    565497#Slope sources  {{{
    566 slope_sources =  ./modules/ModelProcessorx/L2ProjectionBase/UpdateElementsL2ProjectionBase.cpp\
    567                                           ./modules/ModelProcessorx/L2ProjectionBase/CreateNodesL2ProjectionBase.cpp \
    568                                           ./modules/ModelProcessorx/L2ProjectionBase/CreateConstraintsL2ProjectionBase.cpp\
    569                                           ./modules/ModelProcessorx/L2ProjectionBase/CreateLoadsL2ProjectionBase.cpp\
    570                                           ./modules/ModelProcessorx/L2ProjectionTop/UpdateElementsL2ProjectionTop.cpp\
    571                                           ./modules/ModelProcessorx/L2ProjectionTop/CreateNodesL2ProjectionTop.cpp \
    572                                           ./modules/ModelProcessorx/L2ProjectionTop/CreateConstraintsL2ProjectionTop.cpp\
    573                                           ./modules/ModelProcessorx/L2ProjectionTop/CreateLoadsL2ProjectionTop.cpp\
    574                                           ./analyses/L2ProjectionBaseAnalysis.h\
     498slope_sources =  ./analyses/L2ProjectionBaseAnalysis.h\
    575499                                          ./analyses/L2ProjectionBaseAnalysis.cpp\
    576500                                          ./cores/surfaceslope_core.cpp\
    … …  
    578502#}}}
    579503#MeshDeformation sources  {{{
    580 meshdeformation_sources =  ./modules/ModelProcessorx/MeshDeformation/UpdateElementsMeshDeformation.cpp\
    581                                           ./modules/ModelProcessorx/MeshDeformation/CreateNodesMeshDeformation.cpp \
    582                                           ./modules/ModelProcessorx/MeshDeformation/CreateConstraintsMeshDeformation.cpp\
    583                                           ./modules/ModelProcessorx/MeshDeformation/CreateLoadsMeshDeformation.cpp\
    584                                           ./analyses/MeshdeformationAnalysis.h\
     504meshdeformation_sources = ./analyses/MeshdeformationAnalysis.h\
    585505                                          ./analyses/MeshdeformationAnalysis.cpp\
    586506                                          ./cores/meshdeformation_core.cpp
    … …  
    590510                                        ./analyses/GiaAnalysis.h\
    591511                                        ./analyses/GiaAnalysis.cpp\
    592                                         ./modules/ModelProcessorx/Gia/UpdateElementsGia.cpp\
    593                                         ./modules/ModelProcessorx/Gia/CreateNodesGia.cpp \
    594                                         ./modules/ModelProcessorx/Gia/CreateConstraintsGia.cpp\
    595                                         ./modules/ModelProcessorx/Gia/CreateLoadsGia.cpp\
    596512                                        ./modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp\
    597513                                        ./modules/GiaDeflectionCorex/GiaDeflectionCorex.h\
    … …  
    610526                                                ./analyses/DamageEvolutionAnalysis.h\
    611527                                                ./analyses/DamageEvolutionAnalysis.cpp\
    612                                                 ./modules/ModelProcessorx/Damage/UpdateElementsDamage.cpp\
    613                                                 ./modules/ModelProcessorx/Damage/CreateNodesDamage.cpp \
    614                                                 ./modules/ModelProcessorx/Damage/CreateConstraintsDamage.cpp\
    615                                                 ./modules/ModelProcessorx/Damage/CreateParametersDamage.cpp\
    616                                                 ./modules/ModelProcessorx/Damage/CreateLoadsDamage.cpp\
    617528                                                ./solutionsequences/solutionsequence_damage_nonlinear.cpp
    618529

  • issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processor*/
    8 int AdjointBalancethicknessAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  AdjointBalancethicknessAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void AdjointBalancethicknessAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12        _error_("not implemented yet");
     13}/*}}}*/
     14void AdjointBalancethicknessAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     15        _error_("not implemented yet");
     16}/*}}}*/
     17void AdjointBalancethicknessAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     18        _error_("not implemented yet");
     19}/*}}}*/
     20void AdjointBalancethicknessAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     21        _error_("not implemented yet");
     22}/*}}}*/
     23void AdjointBalancethicknessAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     24        _error_("not implemented yet");
     25}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

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

    r16534 r16539  
    66
    77/*Model processing*/
    8 int AdjointHorizAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  AdjointHorizAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        _error_("not implemented");
    1010}/*}}}*/
     11void AdjointHorizAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12           _error_("not implemented yet");
     13}/*}}}*/
     14void AdjointHorizAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     15           _error_("not implemented yet");
     16}/*}}}*/
     17void AdjointHorizAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     18           _error_("not implemented yet");
     19}/*}}}*/
     20void AdjointHorizAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     21           _error_("not implemented yet");
     22}/*}}}*/
     23void AdjointHorizAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     24           _error_("not implemented yet");
     25}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/Analysis.h

    r16534 r16539  
    66#define _ANALYSIS_H_
    77
     8class Parameters;
     9class IoModel;
     10class Elements;
     11class Nodes;
     12class Constraints;
     13class Loads;
     14
    815class Analysis{
    916
    1017        public:
    1118
    12                 virtual         ~Analysis(){};
    13                 virtual int DofsPerNode(int** doflist,int meshtype,int approximation)=0;
     19                virtual      ~Analysis(){};
     20                virtual int  DofsPerNode(int** doflist,int meshtype,int approximation)=0;
     21                virtual void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum)=0;
     22                virtual void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type)=0;
     23                virtual void CreateNodes(Nodes** pnodes,IoModel* iomodel)=0;
     24                virtual void CreateConstraints(Constraints** pconstraints,IoModel* iomodel)=0;
     25                virtual void CreateLoads(Loads** ploads, IoModel* iomodel)=0;
    1426};
    1527#endif

  • issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int BalancethicknessAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  BalancethicknessAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void BalancethicknessAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12}/*}}}*/
     13void BalancethicknessAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     14
     15        int    stabilization,finiteelement;
     16
     17        /*Fetch data needed: */
     18        iomodel->Constant(&stabilization,BalancethicknessStabilizationEnum);
     19
     20        /*Finite element type*/
     21        finiteelement = P1Enum;
     22        if(stabilization==3){
     23                finiteelement = P1DGEnum;
     24        }
     25
     26        /*Update elements: */
     27        int counter=0;
     28        for(int i=0;i<iomodel->numberofelements;i++){
     29                if(iomodel->my_elements[i]){
     30                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     31                        element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement);
     32                        counter++;
     33                }
     34        }
     35
     36        iomodel->FetchDataToInput(elements,ThicknessEnum);
     37        iomodel->FetchDataToInput(elements,SurfaceEnum);
     38        iomodel->FetchDataToInput(elements,BedEnum);
     39        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     40        iomodel->FetchDataToInput(elements,VxEnum);
     41        iomodel->FetchDataToInput(elements,VyEnum);
     42        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
     43        iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
     44        iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum);
     45
     46        if(iomodel->meshtype==Mesh3DEnum){
     47                iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     48                iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     49        }
     50}/*}}}*/
     51void BalancethicknessAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     52
     53        int  stabilization;
     54        iomodel->Constant(&stabilization,BalancethicknessStabilizationEnum);
     55
     56        /*Check in 3d*/
     57        if(stabilization==3 && iomodel->meshtype==Mesh3DEnum) _error_("DG 3d not implemented yet");
     58
     59        /*First fetch data: */
     60        if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     61        if(stabilization!=3){
     62                ::CreateNodes(pnodes,iomodel,BalancethicknessAnalysisEnum,P1Enum);
     63        }
     64        else{
     65                ::CreateNodes(pnodes,iomodel,BalancethicknessAnalysisEnum,P1DGEnum);
     66        }
     67        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     68}/*}}}*/
     69void BalancethicknessAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     70
     71        /*Fetch parameters: */
     72        int    stabilization;   
     73        iomodel->Constant(&stabilization,BalancethicknessStabilizationEnum);
     74
     75        /*Recover pointer: */
     76        Constraints* constraints=*pconstraints;
     77
     78        /*Do not add constraints in DG*/
     79        if(stabilization!=3){
     80                IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,BalancethicknessAnalysisEnum,P1Enum);
     81        }
     82
     83        /*Assign output pointer: */
     84        *pconstraints=constraints;
     85}/*}}}*/
     86void BalancethicknessAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     87
     88        /*Intermediary*/
     89        int element;
     90        int stabilization;
     91
     92        /*Fetch parameters: */
     93        iomodel->Constant(&stabilization,BalancethicknessStabilizationEnum);
     94
     95        /*Recover pointer: */
     96        Loads* loads=*ploads;
     97
     98        /*Loads only in DG*/
     99        if (stabilization==3){
     100
     101                /*Get faces and elements*/
     102                CreateFaces(iomodel);
     103                iomodel->FetchData(1,ThicknessEnum);
     104
     105                /*First load data:*/
     106                for(int i=0;i<iomodel->numberoffaces;i++){
     107
     108                        /*Get left and right elements*/
     109                        element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
     110
     111                        /*Now, if this element is not in the partition, pass: */
     112                        if(!iomodel->my_elements[element]) continue;
     113
     114                        /* Add load */
     115                        loads->AddObject(new Numericalflux(iomodel->loadcounter+i+1,i,i,iomodel,BalancethicknessAnalysisEnum));
     116                }
     117
     118                /*Free data: */
     119                iomodel->DeleteData(1,ThicknessEnum);
     120        }
     121
     122        /*Assign output pointer: */
     123        *ploads=loads;
     124}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/BalancethicknessSoftAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int BalancethicknessSoftAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  BalancethicknessSoftAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        _error_("not implemented");
    1010}/*}}}*/
     11void BalancethicknessSoftAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12           _error_("not implemented yet");
     13}/*}}}*/
     14void BalancethicknessSoftAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     15           _error_("not implemented yet");
     16}/*}}}*/
     17void BalancethicknessSoftAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     18           _error_("not implemented yet");
     19}/*}}}*/
     20void BalancethicknessSoftAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     21           _error_("not implemented yet");
     22}/*}}}*/
     23void BalancethicknessSoftAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     24           _error_("not implemented yet");
     25}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/BalancethicknessSoftAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int BalancevelocityAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  BalancevelocityAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void BalancevelocityAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12}/*}}}*/
     13void BalancevelocityAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     14
     15        /*Update elements: */
     16        int counter=0;
     17        for(int i=0;i<iomodel->numberofelements;i++){
     18                if(iomodel->my_elements[i]){
     19                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     20                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     21                        counter++;
     22                }
     23        }
     24
     25        iomodel->FetchDataToInput(elements,ThicknessEnum);
     26        iomodel->FetchDataToInput(elements,SurfaceEnum);
     27        iomodel->FetchDataToInput(elements,BedEnum);
     28        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     29        iomodel->FetchDataToInput(elements,VxEnum);
     30        iomodel->FetchDataToInput(elements,VyEnum);
     31        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
     32        iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
     33        iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum);
     34
     35        if(iomodel->meshtype==Mesh3DEnum){
     36                iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     37                iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     38        }
     39}/*}}}*/
     40void BalancevelocityAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     41
     42        /*Check in 3d*/
     43        if(iomodel->meshtype==Mesh3DEnum) _error_("DG 3d not implemented yet");
     44
     45        /*First fetch data: */
     46        iomodel->FetchData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationVertexEquationEnum);
     47        ::CreateNodes(pnodes,iomodel,BalancevelocityAnalysisEnum,P1Enum);
     48        iomodel->DeleteData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationVertexEquationEnum);
     49}/*}}}*/
     50void BalancevelocityAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     51
     52        /*Recover pointer: */
     53        Constraints* constraints=*pconstraints;
     54
     55        /*No constraints for now*/
     56        //IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,BalancevelocityAnalysisEnum,P1Enum);
     57
     58        /*Assign output pointer: */
     59        *pconstraints=constraints;
     60}/*}}}*/
     61void BalancevelocityAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     62
     63        /*No loads*/
     64}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

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

    r16534 r16539  
    66
    77/*Model processing*/
    8 int DamageEvolutionAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  DamageEvolutionAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void DamageEvolutionAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12}/*}}}*/
     13void DamageEvolutionAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     14
     15        /*Update elements: */
     16        int counter=0;
     17        for(int i=0;i<iomodel->numberofelements;i++){
     18                if(iomodel->my_elements[i]){
     19                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     20                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     21                        counter++;
     22                }
     23        }
     24
     25        /*What input do I need to run my damage evolution model?*/
     26        iomodel->FetchDataToInput(elements,VxEnum);
     27        iomodel->FetchDataToInput(elements,VyEnum);
     28        iomodel->FetchDataToInput(elements,VzEnum);
     29        iomodel->FetchDataToInput(elements,DamageDEnum);
     30        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     31        iomodel->FetchDataToInput(elements,PressureEnum);
     32
     33        bool dakota_analysis;
     34        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
     35        if(dakota_analysis){
     36                elements->InputDuplicate(DamageDEnum, QmuDamageDEnum);
     37        }
     38}/*}}}*/
     39void DamageEvolutionAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     40
     41        iomodel->FetchData(1,MeshVertexonbedEnum);
     42        ::CreateNodes(pnodes,iomodel,DamageEvolutionAnalysisEnum,P1Enum);
     43        iomodel->DeleteData(1,MeshVertexonbedEnum);
     44}/*}}}*/
     45void DamageEvolutionAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     46
     47        /*Fetch parameters: */
     48        int stabilization;
     49        iomodel->Constant(&stabilization,DamageStabilizationEnum);
     50
     51        /*Recover pointer: */
     52        Constraints* constraints=*pconstraints;
     53
     54        IoModelToConstraintsx(constraints,iomodel,DamageSpcdamageEnum,DamageEvolutionAnalysisEnum,P1Enum);
     55
     56        /*Assign output pointer: */
     57        *pconstraints=constraints;
     58}/*}}}*/
     59void DamageEvolutionAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     60
     61        /*Recover pointer: */
     62        Loads* loads=*ploads;
     63
     64        /*create penalties for nodes: no node can have a damage > 1*/
     65        iomodel->FetchData(1,DamageSpcdamageEnum);
     66        CreateSingleNodeToElementConnectivity(iomodel);
     67
     68        for(int i=0;i<iomodel->numberofvertices;i++){
     69
     70                /*keep only this partition's nodes:*/
     71                if(iomodel->my_vertices[i]){
     72                        if (xIsNan<IssmDouble>(iomodel->Data(DamageSpcdamageEnum)[i])){ //No penalty applied on spc nodes!
     73                                loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,DamageEvolutionAnalysisEnum));
     74                        }
     75                }
     76        }
     77        iomodel->DeleteData(1,DamageSpcdamageEnum);
     78
     79        /*Assign output pointer: */
     80        *ploads=loads;
     81}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int EnthalpyAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  EnthalpyAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void EnthalpyAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12}/*}}}*/
     13void EnthalpyAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     14
     15        bool dakota_analysis;
     16        bool isenthalpy;
     17
     18        /*Now, is the model 3d? otherwise, do nothing: */
     19        if(iomodel->meshtype==Mesh2DhorizontalEnum)return;
     20
     21        /*Is enthalpy requested?*/
     22        iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum);
     23        if(!isenthalpy) return;
     24
     25        /*Fetch data needed: */
     26        iomodel->FetchData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);
     27
     28        /*Update elements: */
     29        int counter=0;
     30        for(int i=0;i<iomodel->numberofelements;i++){
     31                if(iomodel->my_elements[i]){
     32                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     33                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     34                        counter++;
     35                }
     36        }
     37
     38        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
     39
     40        iomodel->FetchDataToInput(elements,ThicknessEnum);
     41        iomodel->FetchDataToInput(elements,SurfaceEnum);
     42        iomodel->FetchDataToInput(elements,BedEnum);
     43        iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
     44        iomodel->FetchDataToInput(elements,FrictionPEnum);
     45        iomodel->FetchDataToInput(elements,FrictionQEnum);
     46        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     47        iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
     48        iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     49        iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     50        iomodel->FetchDataToInput(elements,FlowequationElementEquationEnum);
     51        iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
     52        iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
     53        iomodel->FetchDataToInput(elements,PressureEnum);
     54        iomodel->FetchDataToInput(elements,TemperatureEnum);
     55        iomodel->FetchDataToInput(elements,WaterfractionEnum);
     56        iomodel->FetchDataToInput(elements,EnthalpyEnum);
     57        iomodel->FetchDataToInput(elements,BasalforcingsGeothermalfluxEnum);
     58        iomodel->FetchDataToInput(elements,WatercolumnEnum);
     59        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
     60        iomodel->FetchDataToInput(elements,VxEnum);
     61        iomodel->FetchDataToInput(elements,VyEnum);
     62        iomodel->FetchDataToInput(elements,VzEnum);
     63        InputUpdateFromConstantx(elements,0.,VxMeshEnum);
     64        InputUpdateFromConstantx(elements,0.,VyMeshEnum);
     65        InputUpdateFromConstantx(elements,0.,VzMeshEnum);
     66        if(dakota_analysis){
     67                elements->InputDuplicate(TemperatureEnum,QmuTemperatureEnum);
     68                elements->InputDuplicate(BasalforcingsMeltingRateEnum,QmuMeltingEnum);
     69                elements->InputDuplicate(VxMeshEnum,QmuVxMeshEnum);
     70                elements->InputDuplicate(VxMeshEnum,QmuVyMeshEnum);
     71                elements->InputDuplicate(VxMeshEnum,QmuVzMeshEnum);
     72        }
     73
     74        /*Free data: */
     75        iomodel->DeleteData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);
     76}/*}}}*/
     77void EnthalpyAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     78
     79        if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     80        ::CreateNodes(pnodes,iomodel,EnthalpyAnalysisEnum,P1Enum);
     81        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     82}/*}}}*/
     83void EnthalpyAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     84
     85        /*Intermediary*/
     86        int        count;
     87        int        M,N;
     88        bool       spcpresent = false;
     89        IssmDouble heatcapacity;
     90        IssmDouble referencetemperature;
     91
     92        /*Output*/
     93        IssmDouble *spcvector  = NULL;
     94        IssmDouble* times=NULL;
     95        IssmDouble* values=NULL;
     96
     97        /*Fetch parameters: */
     98        iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum);
     99        iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum);
     100
     101        /*Recover pointer: */
     102        Constraints* constraints=*pconstraints;
     103
     104        /*return if 2d mesh*/
     105        if(iomodel->meshtype==Mesh2DhorizontalEnum) return;
     106
     107        /*Fetch data: */
     108        iomodel->FetchData(&spcvector,&M,&N,ThermalSpctemperatureEnum);
     109
     110        //FIX ME: SHOULD USE IOMODELCREATECONSTRAINTS
     111        /*Transient or static?:*/
     112        if(M==iomodel->numberofvertices){
     113                /*static: just create Constraints objects*/
     114                count=0;
     115
     116                for(int i=0;i<iomodel->numberofvertices;i++){
     117                        /*keep only this partition's nodes:*/
     118                        if((iomodel->my_vertices[i])){
     119
     120                                if (!xIsNan<IssmDouble>(spcvector[i])){
     121
     122                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,heatcapacity*(spcvector[i]-referencetemperature),EnthalpyAnalysisEnum));
     123                                        count++;
     124
     125                                }
     126                        }
     127                }
     128        }
     129        else if (M==(iomodel->numberofvertices+1)){
     130                /*transient: create transient SpcTransient objects. Same logic, except we need to retrieve
     131                 * various times and values to initialize an SpcTransient object: */
     132                count=0;
     133
     134                /*figure out times: */
     135                times=xNew<IssmDouble>(N);
     136                for(int j=0;j<N;j++){
     137                        times[j]=spcvector[(M-1)*N+j];
     138                }
     139
     140                /*Create constraints from x,y,z: */
     141                for(int i=0;i<iomodel->numberofvertices;i++){
     142
     143                        /*keep only this partition's nodes:*/
     144                        if((iomodel->my_vertices[i])){
     145
     146                                /*figure out times and values: */
     147                                values=xNew<IssmDouble>(N);
     148                                spcpresent=false;
     149                                for(int j=0;j<N;j++){
     150                                        values[j]=heatcapacity*(spcvector[i*N+j]-referencetemperature);
     151                                        if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
     152                                }
     153
     154                                if(spcpresent){
     155                                        constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,N,times,values,EnthalpyAnalysisEnum));
     156                                        count++;
     157                                }
     158                                xDelete<IssmDouble>(values);
     159                        }
     160                }
     161        }
     162        else{
     163                _error_("Size of field " << EnumToStringx(ThermalSpctemperatureEnum) << " not supported");
     164        }
     165
     166        /*Free ressources:*/
     167        iomodel->DeleteData(spcvector,ThermalSpctemperatureEnum);
     168        xDelete<IssmDouble>(times);
     169        xDelete<IssmDouble>(values);
     170
     171        /*Assign output pointer: */
     172        *pconstraints=constraints;
     173}/*}}}*/
     174void EnthalpyAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     175
     176        /*No loads */
     177}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/EnumToAnalysis.h

    r16534 r16539  
    77
    88#endif
     9                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     10                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     11                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     12                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     13                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     14                void CreateLoads(Loads** ploads, IoModel* iomodel);

  • issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int ExtrudeFromBaseAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  ExtrudeFromBaseAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void ExtrudeFromBaseAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12}/*}}}*/
     13void ExtrudeFromBaseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     14
     15        int counter=0;
     16        for(int i=0;i<iomodel->numberofelements;i++){
     17                if(iomodel->my_elements[i]){
     18                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     19                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     20                        counter++;
     21                }
     22        }
     23
     24        if(iomodel->meshtype==Mesh2DverticalEnum){
     25                iomodel->FetchDataToInput(elements,MeshVertexonbedEnum);
     26        }
     27}/*}}}*/
     28void ExtrudeFromBaseAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     29
     30        ::CreateNodes(pnodes,iomodel,ExtrudeFromBaseAnalysisEnum,P1Enum);
     31
     32}/*}}}*/
     33void ExtrudeFromBaseAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     34}/*}}}*/
     35void ExtrudeFromBaseAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     36}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int ExtrudeFromTopAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  ExtrudeFromTopAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void ExtrudeFromTopAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12}/*}}}*/
     13void ExtrudeFromTopAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     14
     15        int counter=0;
     16        for(int i=0;i<iomodel->numberofelements;i++){
     17                if(iomodel->my_elements[i]){
     18                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     19                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     20                        counter++;
     21                }
     22        }
     23
     24        if(iomodel->meshtype==Mesh2DverticalEnum){
     25                iomodel->FetchDataToInput(elements,MeshVertexonbedEnum);
     26        }
     27}/*}}}*/
     28void ExtrudeFromTopAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     29
     30        ::CreateNodes(pnodes,iomodel,ExtrudeFromBaseAnalysisEnum,P1Enum);
     31
     32}/*}}}*/
     33void ExtrudeFromTopAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     34}/*}}}*/
     35void ExtrudeFromTopAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     36}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

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

    r16534 r16539  
    66
    77/*Model processing*/
    8 int FreeSurfaceBaseAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  FreeSurfaceBaseAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void FreeSurfaceBaseAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12}/*}}}*/
     13void FreeSurfaceBaseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     14
     15        /*Now, is the model 3d? otherwise, do nothing: */
     16        if (iomodel->meshtype==Mesh2DhorizontalEnum)return;
     17
     18        /*Finite element type*/
     19        int finiteelement = P1Enum;
     20
     21        /*Update elements: */
     22        int counter=0;
     23        for(int i=0;i<iomodel->numberofelements;i++){
     24                if(iomodel->my_elements[i]){
     25                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     26                        element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement);
     27                        counter++;
     28                }
     29        }
     30
     31        iomodel->FetchDataToInput(elements,SurfaceEnum);
     32        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     33        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
     34        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateCorrectionEnum,0.);
     35        iomodel->FetchDataToInput(elements,VxEnum);
     36        iomodel->FetchDataToInput(elements,VyEnum);
     37        if(iomodel->meshtype==Mesh3DEnum){
     38                iomodel->FetchDataToInput(elements,VzEnum);
     39                iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     40                iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     41        }
     42}/*}}}*/
     43void FreeSurfaceBaseAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     44
     45        if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     46        ::CreateNodes(pnodes,iomodel,FreeSurfaceBaseAnalysisEnum,P1Enum);
     47        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     48}/*}}}*/
     49void FreeSurfaceBaseAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     50}/*}}}*/
     51void FreeSurfaceBaseAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     52
     53        /*Intermediaries*/
     54        int penpair_ids[2];
     55        int count=0;
     56        int numvertex_pairing;
     57
     58        /*Recover pointer: */
     59        Loads* loads=*ploads;
     60
     61        /*Create Penpair for vertex_pairing: */
     62        IssmDouble *vertex_pairing=NULL;
     63        IssmDouble *nodeonbed=NULL;
     64        iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum);
     65        iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
     66        for(int i=0;i<numvertex_pairing;i++){
     67
     68                if(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+0])-1]){
     69
     70                        /*In debugging mode, check that the second node is in the same cpu*/
     71                        _assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);
     72
     73                        /*Skip if one of the two is not on the bed*/
     74                        if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
     75
     76                        /*Get node ids*/
     77                        penpair_ids[0]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+0]);
     78                        penpair_ids[1]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+1]);
     79
     80                        /*Create Load*/
     81                        loads->AddObject(new Penpair(
     82                                                        iomodel->loadcounter+count+1,
     83                                                        &penpair_ids[0],
     84                                                        FreeSurfaceBaseAnalysisEnum));
     85                        count++;
     86                }
     87        }
     88
     89        /*free ressources: */
     90        iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum);
     91        iomodel->DeleteData(nodeonbed,MeshVertexonbedEnum);
     92
     93        /*Assign output pointer: */
     94        *ploads=loads;
     95}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

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

    r16534 r16539  
    66
    77/*Model processing*/
    8 int FreeSurfaceTopAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  FreeSurfaceTopAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void FreeSurfaceTopAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12}/*}}}*/
     13void FreeSurfaceTopAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     14
     15        /*Now, is the model 3d? otherwise, do nothing: */
     16        if (iomodel->meshtype==Mesh2DhorizontalEnum)return;
     17
     18        int finiteelement = P1Enum;
     19
     20        /*Update elements: */
     21        int counter=0;
     22        for(int i=0;i<iomodel->numberofelements;i++){
     23                if(iomodel->my_elements[i]){
     24                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     25                        element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement);
     26                        counter++;
     27                }
     28        }
     29
     30        iomodel->FetchDataToInput(elements,SurfaceEnum);
     31        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     32        iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum,0.);
     33        iomodel->FetchDataToInput(elements,VxEnum);
     34        iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
     35        if(iomodel->meshtype==Mesh3DEnum){
     36                iomodel->FetchDataToInput(elements,VzEnum);
     37                iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     38                iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     39        }
     40}/*}}}*/
     41void FreeSurfaceTopAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     42
     43        if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     44        ::CreateNodes(pnodes,iomodel,FreeSurfaceTopAnalysisEnum,P1Enum);
     45        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     46}/*}}}*/
     47void FreeSurfaceTopAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     48}/*}}}*/
     49void FreeSurfaceTopAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     50
     51        /*Intermediaries*/
     52        int penpair_ids[2];
     53        int count=0;
     54        int numvertex_pairing;
     55
     56        /*Recover pointer: */
     57        Loads* loads=*ploads;
     58
     59        /*Create Penpair for vertex_pairing: */
     60        IssmDouble *vertex_pairing=NULL;
     61        IssmDouble *nodeonsurface=NULL;
     62        iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum);
     63        iomodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum);
     64        for(int i=0;i<numvertex_pairing;i++){
     65
     66                if(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+0])-1]){
     67
     68                        /*In debugging mode, check that the second node is in the same cpu*/
     69                        _assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);
     70
     71                        /*Skip if one of the two is not on the bed*/
     72                        if(!(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
     73
     74                        /*Get node ids*/
     75                        penpair_ids[0]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+0]);
     76                        penpair_ids[1]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+1]);
     77
     78                        /*Create Load*/
     79                        loads->AddObject(new Penpair(
     80                                                        iomodel->loadcounter+count+1,
     81                                                        &penpair_ids[0],
     82                                                        FreeSurfaceTopAnalysisEnum));
     83                        count++;
     84                }
     85        }
     86
     87        /*free ressources: */
     88        iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum);
     89        iomodel->DeleteData(nodeonsurface,MeshVertexonsurfaceEnum);
     90
     91        /*Assign output pointer: */
     92        *ploads=loads;
     93}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/GiaAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int GiaAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  GiaAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void GiaAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12}/*}}}*/
     13void GiaAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     14
     15        /*Update elements: */
     16        int counter=0;
     17        for(int i=0;i<iomodel->numberofelements;i++){
     18                if(iomodel->my_elements[i]){
     19                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     20                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     21                        counter++;
     22                }
     23        }
     24
     25        iomodel->FetchDataToInput(elements,ThicknessEnum);
     26        iomodel->FetchDataToInput(elements,GiaMantleViscosityEnum);
     27        iomodel->FetchDataToInput(elements,GiaLithosphereThicknessEnum);
     28}/*}}}*/
     29void GiaAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     30        ::CreateNodes(pnodes,iomodel,GiaAnalysisEnum,P1Enum);
     31}/*}}}*/
     32void GiaAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     33        /*No constraints*/
     34}/*}}}*/
     35void GiaAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     36        /*No loads*/
     37}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/GiaAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int HydrologyDCEfficientAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  HydrologyDCEfficientAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void HydrologyDCEfficientAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12
     13        Parameters *parameters = NULL;
     14        int         hydrology_model;
     15        bool        isefficientlayer;
     16
     17        /*Get parameters: */
     18        parameters=*pparameters;
     19
     20        /*retrieve some parameters: */
     21        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
     22
     23        /*Now, do we really want DC?*/
     24        if(hydrology_model!=HydrologydcEnum){
     25                *pparameters=parameters;
     26                return;
     27        }
     28
     29        /*Do we want an efficient layer*/
     30        iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     31        if(!isefficientlayer){
     32                *pparameters=parameters;
     33                return;
     34        }
     35
     36        /*Nothing for now*/
     37
     38        /*Assign output pointer: */
     39        *pparameters=parameters;
     40}/*}}}*/
     41void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     42
     43        bool   isefficientlayer;
     44        int    hydrology_model;
     45
     46        /*Now, do we really want DC?*/
     47        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
     48        if(hydrology_model!=HydrologydcEnum) return;
     49
     50        /*Do we want an efficient layer*/
     51        iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     52        if(!isefficientlayer) return;
     53
     54        /*Update elements: */
     55        int counter=0;
     56        for(int i=0;i<iomodel->numberofelements;i++){
     57                if(iomodel->my_elements[i]){
     58                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     59                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     60                        counter++;
     61                }
     62        }
     63
     64        iomodel->FetchDataToInput(elements,ThicknessEnum);
     65        iomodel->FetchDataToInput(elements,SurfaceEnum);
     66        iomodel->FetchDataToInput(elements,BedEnum);
     67        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     68        iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     69        iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     70        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
     71        iomodel->FetchDataToInput(elements,EplHeadEnum);
     72}/*}}}*/
     73void HydrologyDCEfficientAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     74
     75        /*Now, do we really want DC?*/
     76        int  hydrology_model;
     77        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
     78        if(hydrology_model!=HydrologydcEnum) return;
     79
     80        /*Do we want an efficient layer*/
     81        bool isefficientlayer;
     82        iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     83        if(!isefficientlayer) return;
     84
     85        iomodel->FetchData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationVertexEquationEnum);
     86        ::CreateNodes(pnodes,iomodel,HydrologyDCEfficientAnalysisEnum,P1Enum);
     87        iomodel->DeleteData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationVertexEquationEnum);
     88}/*}}}*/
     89void HydrologyDCEfficientAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     90
     91        /*Recover pointer: */
     92        Constraints* constraints=*pconstraints;
     93
     94        /*Do we really want DC?*/
     95        int  hydrology_model;
     96        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
     97        if(hydrology_model!=HydrologydcEnum) return;
     98
     99        /*Do we want an efficient layer*/
     100        bool isefficientlayer;
     101        iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     102        if(!isefficientlayer) return;
     103
     104        IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum,P1Enum);
     105
     106        /*Assign output pointer: */
     107        *pconstraints=constraints;
     108}/*}}}*/
     109void HydrologyDCEfficientAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     110
     111        /*Recover pointer: */
     112        Loads* loads=*ploads;
     113
     114        /*Do we really want DC?*/
     115        int  hydrology_model;
     116        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
     117        if(hydrology_model!=HydrologydcEnum) return;
     118
     119        /*Do we want an efficient layer*/
     120        bool isefficientlayer;
     121        iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     122        if(!isefficientlayer) return;
     123
     124        /*Nothing for now*/
     125
     126        /*Assign output pointer: */
     127        *ploads=loads;
     128}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int HydrologyDCInefficientAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  HydrologyDCInefficientAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void HydrologyDCInefficientAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12
     13        Parameters *parameters = NULL;
     14        int         hydrology_model;
     15        int         sedimentlimit_flag;
     16        int         transfer_flag;
     17        bool        isefficientlayer;
     18        IssmDouble  sedimentlimit;
     19        IssmDouble  penalty_factor;
     20        IssmDouble  leakagefactor;
     21        IssmDouble  rel_tol;
     22
     23        /*Get parameters: */
     24        parameters=*pparameters;
     25
     26        /*retrieve some parameters: */
     27        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
     28
     29        /*Now, do we really want DC?*/
     30        if(hydrology_model!=HydrologydcEnum){
     31                *pparameters=parameters;
     32                return;
     33        }
     34
     35        iomodel->FetchData(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     36        iomodel->FetchData(&sedimentlimit_flag,HydrologydcSedimentlimitFlagEnum);
     37        iomodel->FetchData(&transfer_flag,HydrologydcTransferFlagEnum);
     38        iomodel->FetchData(&penalty_factor,HydrologydcPenaltyFactorEnum);
     39        iomodel->FetchData(&rel_tol,HydrologydcRelTolEnum);
     40
     41        if(sedimentlimit_flag==1){
     42                iomodel->FetchData(&sedimentlimit,HydrologydcSedimentlimitEnum);
     43                parameters->AddObject(new DoubleParam(HydrologydcSedimentlimitEnum,sedimentlimit));
     44        }
     45
     46        if(transfer_flag==1){
     47                iomodel->FetchData(&leakagefactor,HydrologydcLeakageFactorEnum);
     48                parameters->AddObject(new DoubleParam(HydrologydcLeakageFactorEnum,leakagefactor));
     49        }
     50
     51        parameters->AddObject(new DoubleParam(HydrologydcPenaltyFactorEnum,penalty_factor));
     52        parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
     53        parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer));
     54        parameters->AddObject(new IntParam(HydrologydcSedimentlimitFlagEnum,sedimentlimit_flag));
     55        parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag));
     56        parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol));
     57
     58        /*Assign output pointer: */
     59        *pparameters=parameters;
     60}/*}}}*/
     61void HydrologyDCInefficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     62
     63        bool   isefficientlayer;
     64        int    hydrology_model;
     65
     66        /*Fetch data needed: */
     67        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
     68
     69        /*Now, do we really want DC?*/
     70        if(hydrology_model!=HydrologydcEnum) return;
     71
     72        /*Fetch data needed: */
     73        iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     74
     75        /*Update elements: */
     76        int counter=0;
     77        for(int i=0;i<iomodel->numberofelements;i++){
     78                if(iomodel->my_elements[i]){
     79                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     80                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     81                        counter++;
     82                }
     83        }
     84
     85        iomodel->FetchDataToInput(elements,ThicknessEnum);
     86        iomodel->FetchDataToInput(elements,SurfaceEnum);
     87        iomodel->FetchDataToInput(elements,BedEnum);
     88        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     89        iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     90        iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     91        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
     92        iomodel->FetchDataToInput(elements,SedimentHeadEnum);
     93        if(isefficientlayer)iomodel->FetchDataToInput(elements,HydrologydcMaskEplactiveEnum);
     94}/*}}}*/
     95void HydrologyDCInefficientAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     96
     97        /*Fetch parameters: */
     98        int  hydrology_model;
     99        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
     100
     101        /*Now, do we really want DC?*/
     102        if(hydrology_model!=HydrologydcEnum) return;
     103
     104        if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     105        ::CreateNodes(pnodes,iomodel,HydrologyDCInefficientAnalysisEnum,P1Enum);
     106        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     107}/*}}}*/
     108void HydrologyDCInefficientAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     109
     110        /*Recover pointer: */
     111        Constraints* constraints=*pconstraints;
     112
     113        /*retrieve some parameters: */
     114        int hydrology_model;
     115        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
     116        if(hydrology_model!=HydrologydcEnum) return;
     117
     118        IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyDCInefficientAnalysisEnum,P1Enum);
     119
     120        /*Assign output pointer: */
     121        *pconstraints=constraints;
     122}/*}}}*/
     123void HydrologyDCInefficientAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     124
     125        /*Recover pointer: */
     126        Loads* loads=*ploads;
     127
     128        /*Fetch parameters: */
     129        int hydrology_model;
     130        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
     131        if(hydrology_model!=HydrologydcEnum) return;
     132
     133        iomodel->FetchData(1,MeshVertexonbedEnum);
     134
     135        //create penalties for nodes: no node can have a temperature over the melting point
     136        CreateSingleNodeToElementConnectivity(iomodel);
     137        for(int i=0;i<iomodel->numberofvertices;i++){
     138                if (iomodel->meshtype==Mesh3DEnum){
     139                        /*keep only this partition's nodes:*/
     140                        if(iomodel->my_vertices[i]){
     141                                loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum));
     142                        }
     143                }
     144                else if(reCast<int>(iomodel->Data(MeshVertexonbedEnum)[i])){
     145                        if(iomodel->my_vertices[i]){
     146                                loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum));
     147                        }       
     148                }
     149        }
     150        /*Assign output pointer: */
     151        *ploads=loads;
     152        iomodel->DeleteData(1,MeshVertexonbedEnum);
     153}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int HydrologyShreveAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  HydrologyShreveAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void HydrologyShreveAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12
     13        Parameters *parameters = NULL;
     14        int         hydrology_model;
     15
     16        /*Get parameters: */
     17        parameters=*pparameters;
     18
     19        /*retrieve some parameters: */
     20        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
     21
     22        /*Now, do we really want Shreve?*/
     23        if(hydrology_model!=HydrologyshreveEnum){
     24                *pparameters=parameters;
     25                return;
     26        }
     27
     28        parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
     29        parameters->AddObject(iomodel->CopyConstantObject(HydrologyshreveStabilizationEnum));
     30
     31        /*Assign output pointer: */
     32        *pparameters=parameters;
     33}/*}}}*/
     34void HydrologyShreveAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     35
     36        /*Fetch data needed: */
     37        int    hydrology_model;
     38        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
     39
     40        /*Now, do we really want Shreve?*/
     41        if(hydrology_model!=HydrologyshreveEnum) return;
     42
     43        /*Update elements: */
     44        int counter=0;
     45        for(int i=0;i<iomodel->numberofelements;i++){
     46                if(iomodel->my_elements[i]){
     47                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     48                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     49                        counter++;
     50                }
     51        }
     52
     53        iomodel->FetchDataToInput(elements,ThicknessEnum);
     54        iomodel->FetchDataToInput(elements,SurfaceEnum);
     55        iomodel->FetchDataToInput(elements,BedEnum);
     56        iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     57        iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     58        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     59        iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
     60        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
     61        iomodel->FetchDataToInput(elements,WatercolumnEnum);
     62
     63        elements->InputDuplicate(WatercolumnEnum,WaterColumnOldEnum);
     64}/*}}}*/
     65void HydrologyShreveAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     66
     67        /*Fetch parameters: */
     68        int  hydrology_model;
     69        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
     70
     71        /*Now, do we really want Shreve?*/
     72        if(hydrology_model!=HydrologyshreveEnum) return;
     73
     74        if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     75        ::CreateNodes(pnodes,iomodel,HydrologyShreveAnalysisEnum,P1Enum);
     76        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     77}/*}}}*/
     78void HydrologyShreveAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     79
     80        /*Recover pointer: */
     81        Constraints* constraints=*pconstraints;
     82
     83        /*retrieve some parameters: */
     84        int          hydrology_model;
     85        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
     86
     87        if(hydrology_model!=HydrologyshreveEnum) return;
     88
     89        IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum,P1Enum);
     90
     91        /*Assign output pointer: */
     92        *pconstraints=constraints;
     93}/*}}}*/
     94void HydrologyShreveAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     95        /*No loads*/
     96}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int L2ProjectionBaseAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  L2ProjectionBaseAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void L2ProjectionBaseAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12}/*}}}*/
     13void L2ProjectionBaseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     14
     15        /*Update elements: */
     16        int counter=0;
     17        for(int i=0;i<iomodel->numberofelements;i++){
     18                if(iomodel->my_elements[i]){
     19                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     20                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     21                        counter++;
     22                }
     23        }
     24
     25        iomodel->FetchDataToInput(elements,SurfaceEnum);
     26        iomodel->FetchDataToInput(elements,BedEnum);
     27        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     28        if(iomodel->meshtype==Mesh3DEnum){
     29                iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     30                iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     31        }
     32        if(iomodel->meshtype==Mesh2DverticalEnum){
     33                iomodel->FetchDataToInput(elements,MeshVertexonbedEnum);
     34        }
     35}/*}}}*/
     36void L2ProjectionBaseAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     37
     38        if(iomodel->meshtype==Mesh3DEnum){
     39                iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     40        }
     41        else if(iomodel->meshtype==Mesh2DverticalEnum){
     42                iomodel->FetchData(1,MeshVertexonbedEnum);
     43        }
     44        ::CreateNodes(pnodes,iomodel,L2ProjectionBaseAnalysisEnum,P1Enum);
     45        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     46}/*}}}*/
     47void L2ProjectionBaseAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     48
     49        /*No constraints*/
     50}/*}}}*/
     51void L2ProjectionBaseAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     52
     53        /*No loads*/
     54}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

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

    r16534 r16539  
    66
    77/*Model processing*/
    8 int MasstransportAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  MasstransportAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void MasstransportAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12
     13        /*Get parameters: */
     14        Parameters *parameters=*pparameters;
     15
     16        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum));
     17
     18        /*Assign output pointer: */
     19        *pparameters = parameters;
     20}/*}}}*/
     21void MasstransportAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     22
     23        int    stabilization,finiteelement;
     24        bool   dakota_analysis;
     25        bool   issmbgradients;
     26        bool   ispdd;
     27        bool   isdelta18o;
     28        bool   isgroundingline;
     29
     30        /*Fetch data needed: */
     31        iomodel->Constant(&stabilization,MasstransportStabilizationEnum);
     32        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
     33        iomodel->Constant(&ispdd,SurfaceforcingsIspddEnum);
     34        iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
     35        iomodel->Constant(&issmbgradients,SurfaceforcingsIssmbgradientsEnum);
     36        iomodel->Constant(&isgroundingline,TransientIsgroundinglineEnum);
     37
     38        /*Finite element type*/
     39        finiteelement = P1Enum;
     40        if(stabilization==3){
     41                finiteelement = P1DGEnum;
     42        }
     43
     44        /*Update elements: */
     45        int counter=0;
     46        for(int i=0;i<iomodel->numberofelements;i++){
     47                if(iomodel->my_elements[i]){
     48                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     49                        element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement);
     50                        counter++;
     51                }
     52        }
     53
     54        iomodel->FetchDataToInput(elements,ThicknessEnum);
     55        iomodel->FetchDataToInput(elements,SurfaceEnum);
     56        iomodel->FetchDataToInput(elements,BedEnum);
     57        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     58        iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
     59        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
     60        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateCorrectionEnum,0.);
     61        iomodel->FetchDataToInput(elements,VxEnum);
     62        iomodel->FetchDataToInput(elements,VyEnum);
     63
     64        if(stabilization==3){
     65                iomodel->FetchDataToInput(elements,MasstransportSpcthicknessEnum); //for DG, we need the spc in the element
     66        }
     67
     68        if(dakota_analysis){
     69                elements->InputDuplicate(BedEnum,QmuBedEnum);
     70                elements->InputDuplicate(ThicknessEnum,QmuThicknessEnum);
     71                elements->InputDuplicate(SurfaceEnum,QmuSurfaceEnum);
     72                elements->InputDuplicate(MaskIceLevelsetEnum,QmuMaskIceLevelsetEnum);
     73                if(isgroundingline) elements->InputDuplicate(MaskGroundediceLevelsetEnum,QmuMaskGroundediceLevelsetEnum);
     74        }
     75
     76        if(iomodel->meshtype==Mesh3DEnum){
     77                iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     78                iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     79        }
     80
     81        if(issmbgradients){
     82                iomodel->FetchDataToInput(elements,SurfaceforcingsHrefEnum);
     83                iomodel->FetchDataToInput(elements,SurfaceforcingsSmbrefEnum);
     84                iomodel->FetchDataToInput(elements,SurfaceforcingsBPosEnum);
     85                iomodel->FetchDataToInput(elements,SurfaceforcingsBNegEnum);
     86        }
     87        if(ispdd){
     88                iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum);
     89                if(isdelta18o){
     90                        iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesLgmEnum);
     91                        iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesPresentdayEnum);
     92                        iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsPresentdayEnum);
     93                }
     94                else{
     95                        iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum);
     96                        iomodel->FetchDataToInput(elements,SurfaceforcingsMonthlytemperaturesEnum);
     97                }
     98        }
     99        if(~ispdd && ~issmbgradients){
     100                iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum,0.);
     101        }
     102}/*}}}*/
     103void MasstransportAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     104
     105        /*Fetch parameters: */
     106        int  stabilization;
     107        iomodel->Constant(&stabilization,MasstransportStabilizationEnum);
     108
     109        /*Check in 3d*/
     110        if(stabilization==3 && iomodel->meshtype==Mesh3DEnum) _error_("DG 3d not implemented yet");
     111
     112        /*Create Nodes either DG or CG depending on stabilization*/
     113        if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     114        if(stabilization!=3){
     115                ::CreateNodes(pnodes,iomodel,MasstransportAnalysisEnum,P1Enum);
     116        }
     117        else{
     118                ::CreateNodes(pnodes,iomodel,MasstransportAnalysisEnum,P1DGEnum);
     119        }
     120        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     121}/*}}}*/
     122void MasstransportAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     123
     124        /*Fetch parameters: */
     125        int stabilization;
     126        iomodel->Constant(&stabilization,MasstransportStabilizationEnum);
     127
     128        /*Recover pointer: */
     129        Constraints* constraints=*pconstraints;
     130
     131        /*Do not add constraints in DG, they are weakly imposed*/
     132        if(stabilization!=3){
     133                IoModelToConstraintsx(constraints,iomodel,MasstransportSpcthicknessEnum,MasstransportAnalysisEnum,P1Enum);
     134        }
     135
     136        /*Assign output pointer: */
     137        *pconstraints=constraints;
     138}/*}}}*/
     139void MasstransportAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     140
     141        /*Intermediaries*/
     142        int element;
     143        int penpair_ids[2];
     144        int count=0;
     145        int stabilization;
     146        int numvertex_pairing;
     147
     148        /*Fetch parameters: */
     149        iomodel->Constant(&stabilization,MasstransportStabilizationEnum);
     150
     151        /*Recover pointer: */
     152        Loads* loads=*ploads;
     153
     154        /*Loads only in DG*/
     155        if (stabilization==3){
     156
     157                /*Get faces and elements*/
     158                CreateFaces(iomodel);
     159                iomodel->FetchData(1,ThicknessEnum);
     160
     161                /*First load data:*/
     162                for(int i=0;i<iomodel->numberoffaces;i++){
     163
     164                        /*Get left and right elements*/
     165                        element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
     166
     167                        /*Now, if this element is not in the partition, pass: */
     168                        if(!iomodel->my_elements[element]) continue;
     169
     170                        /* Add load */
     171                        loads->AddObject(new Numericalflux(iomodel->loadcounter+i+1,i,i,iomodel,MasstransportAnalysisEnum));
     172                }
     173
     174                /*Free data: */
     175                iomodel->DeleteData(1,ThicknessEnum);
     176        }
     177
     178        /*Create Penpair for vertex_pairing: */
     179        IssmDouble *vertex_pairing=NULL;
     180        IssmDouble *nodeonbed=NULL;
     181        iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum);
     182        iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
     183
     184        for(int i=0;i<numvertex_pairing;i++){
     185
     186                if(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+0])-1]){
     187
     188                        /*In debugging mode, check that the second node is in the same cpu*/
     189                        _assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);
     190
     191                        /*Skip if one of the two is not on the bed*/
     192                        if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
     193
     194                        /*Get node ids*/
     195                        penpair_ids[0]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+0]);
     196                        penpair_ids[1]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+1]);
     197
     198                        /*Create Load*/
     199                        loads->AddObject(new Penpair(
     200                                                        iomodel->loadcounter+count+1,
     201                                                        &penpair_ids[0],
     202                                                        MasstransportAnalysisEnum));
     203                        count++;
     204                }
     205        }
     206
     207        /*free ressources: */
     208        iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum);
     209        iomodel->DeleteData(nodeonbed,MeshVertexonbedEnum);
     210
     211        /*Assign output pointer: */
     212        *ploads=loads;
     213}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int MeltingAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  MeltingAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void MeltingAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12}/*}}}*/
     13void MeltingAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     14
     15        /*Now, is the model 3d? otherwise, do nothing: */
     16        if(iomodel->meshtype==Mesh2DhorizontalEnum)return;
     17
     18        /*Update elements: */
     19        int counter=0;
     20        for(int i=0;i<iomodel->numberofelements;i++){
     21                if(iomodel->my_elements[i]){
     22                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     23                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     24                        counter++;
     25                }
     26        }
     27
     28        /*Create inputs: */
     29        iomodel->FetchDataToInput(elements,ThicknessEnum);
     30        iomodel->FetchDataToInput(elements,SurfaceEnum);
     31        iomodel->FetchDataToInput(elements,BedEnum);
     32        iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
     33        iomodel->FetchDataToInput(elements,FrictionPEnum);
     34        iomodel->FetchDataToInput(elements,FrictionQEnum);
     35        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     36        iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     37        iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     38        iomodel->FetchDataToInput(elements,FlowequationElementEquationEnum);
     39        iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
     40        iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
     41        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
     42        iomodel->FetchDataToInput(elements,PressureEnum);
     43}/*}}}*/
     44void MeltingAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     45
     46        if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     47        ::CreateNodes(pnodes,iomodel,MeltingAnalysisEnum,P1Enum);
     48        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     49}/*}}}*/
     50void MeltingAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     51        /*No Constraints*/
     52}/*}}}*/
     53void MeltingAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     54
     55        /*if 2d: Error*/
     56        if(iomodel->meshtype==Mesh2DhorizontalEnum) _error_("2d meshes not supported yet");
     57
     58        /*Recover pointer: */
     59        Loads* loads=*ploads;
     60
     61        //create penalties for nodes: no node can have a temperature over the melting point
     62        iomodel->FetchData(1,MeshVertexonbedEnum);
     63        CreateSingleNodeToElementConnectivity(iomodel);
     64
     65        for(int i=0;i<iomodel->numberofvertices;i++){
     66                if(iomodel->my_vertices[i]){
     67                        if (reCast<int>(iomodel->Data(MeshVertexonbedEnum)[i])){
     68                                loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,MeltingAnalysisEnum));
     69                        }
     70                }
     71        }
     72        iomodel->DeleteData(1,MeshVertexonbedEnum);
     73
     74        /*Assign output pointer: */
     75        *ploads=loads;
     76}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/MeltingAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/MeshdeformationAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int MeshdeformationAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  MeshdeformationAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        _error_("not implemented");
    1010}/*}}}*/
     11void MeshdeformationAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12           _error_("not implemented yet");
     13}/*}}}*/
     14void MeshdeformationAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     15           _error_("not implemented yet");
     16}/*}}}*/
     17void MeshdeformationAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     18           _error_("not implemented yet");
     19}/*}}}*/
     20void MeshdeformationAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     21           _error_("not implemented yet");
     22}/*}}}*/
     23void MeshdeformationAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     24           _error_("not implemented yet");
     25}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/MeshdeformationAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int SmoothedSurfaceSlopeXAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  SmoothedSurfaceSlopeXAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void SmoothedSurfaceSlopeXAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12}/*}}}*/
     13void SmoothedSurfaceSlopeXAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     14
     15        /*Update elements: */
     16        int counter=0;
     17        for(int i=0;i<iomodel->numberofelements;i++){
     18                if(iomodel->my_elements[i]){
     19                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     20                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     21                        counter++;
     22                }
     23        }
     24
     25        iomodel->FetchDataToInput(elements,SurfaceEnum);
     26        iomodel->FetchDataToInput(elements,BedEnum);
     27        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     28        if(iomodel->meshtype==Mesh3DEnum){
     29                iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     30                iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     31        }
     32        if(iomodel->meshtype==Mesh2DverticalEnum){
     33                iomodel->FetchDataToInput(elements,MeshVertexonbedEnum);
     34        }
     35}/*}}}*/
     36void SmoothedSurfaceSlopeXAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     37        iomodel->FetchData(1,MeshVertexonbedEnum);
     38        ::CreateNodes(pnodes,iomodel,SmoothedSurfaceSlopeXAnalysisEnum,P1Enum);
     39        iomodel->DeleteData(1,MeshVertexonbedEnum);
     40}/*}}}*/
     41void SmoothedSurfaceSlopeXAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     42}/*}}}*/
     43void SmoothedSurfaceSlopeXAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     44}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int SmoothedSurfaceSlopeYAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  SmoothedSurfaceSlopeYAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void SmoothedSurfaceSlopeYAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12}/*}}}*/
     13void SmoothedSurfaceSlopeYAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     14
     15        /*Update elements: */
     16        int counter=0;
     17        for(int i=0;i<iomodel->numberofelements;i++){
     18                if(iomodel->my_elements[i]){
     19                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     20                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     21                        counter++;
     22                }
     23        }
     24
     25        iomodel->FetchDataToInput(elements,SurfaceEnum);
     26        iomodel->FetchDataToInput(elements,BedEnum);
     27        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     28        if(iomodel->meshtype==Mesh3DEnum){
     29                iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     30                iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     31        }
     32        if(iomodel->meshtype==Mesh2DverticalEnum){
     33                iomodel->FetchDataToInput(elements,MeshVertexonbedEnum);
     34        }
     35}/*}}}*/
     36void SmoothedSurfaceSlopeYAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     37        iomodel->FetchData(1,MeshVertexonbedEnum);
     38        ::CreateNodes(pnodes,iomodel,SmoothedSurfaceSlopeYAnalysisEnum,P1Enum);
     39        iomodel->DeleteData(1,MeshVertexonbedEnum);
     40}/*}}}*/
     41void SmoothedSurfaceSlopeYAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     42}/*}}}*/
     43void SmoothedSurfaceSlopeYAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     44}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

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

    r16534 r16539  
    66
    77/*Model processing*/
    8 int StressbalanceAnalysis::DofsPerNode(int** pdoftype,int meshtype,int approximation){/*{{{*/
     8int  StressbalanceAnalysis::DofsPerNode(int** pdoftype,int meshtype,int approximation){/*{{{*/
    99
    1010        /*output*/
    … …  
    6666        return numdofs;
    6767}/*}}}*/
     68void StressbalanceAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     69
     70        /*Intermediaries*/
     71        int         numoutputs;
     72        char**      requestedoutputs = NULL;
     73
     74        /*Get parameters: */
     75        Parameters *parameters=*pparameters;
     76
     77        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSIAEnum));
     78        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSSAEnum));
     79        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsL1L2Enum));
     80        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsHOEnum));
     81        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum));
     82        parameters->AddObject(iomodel->CopyConstantObject(FlowequationFeFSEnum));
     83        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum));
     84        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum));
     85        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum));
     86        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceIsnewtonEnum));
     87        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum));
     88        parameters->AddObject(iomodel->CopyConstantObject(StressbalancePenaltyFactorEnum));
     89        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum));
     90        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceFSreconditioningEnum));
     91        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceShelfDampeningEnum));
     92        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceViscosityOvershootEnum));
     93
     94        /*Requested outputs*/
     95        iomodel->FetchData(&requestedoutputs,&numoutputs,StressbalanceRequestedOutputsEnum);
     96        parameters->AddObject(new IntParam(StressbalanceNumRequestedOutputsEnum,numoutputs));
     97        if(numoutputs)parameters->AddObject(new StringArrayParam(StressbalanceRequestedOutputsEnum,requestedoutputs,numoutputs));
     98        iomodel->DeleteData(&requestedoutputs,numoutputs,StressbalanceRequestedOutputsEnum);
     99
     100        /*Assign output pointer: */
     101        *pparameters = parameters;
     102}/*}}}*/
     103void StressbalanceAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     104
     105        /*Intermediaries*/
     106        int    materials_type,finiteelement;
     107        int    approximation;
     108        int*   finiteelement_list=NULL;
     109        bool   isSSA,isL1L2,isHO,isFS,iscoupling;
     110        bool   control_analysis;
     111        bool   dakota_analysis;
     112
     113        /*Fetch constants needed: */
     114        iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
     115        iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
     116        iomodel->Constant(&isHO,FlowequationIsHOEnum);
     117        iomodel->Constant(&isFS,FlowequationIsFSEnum);
     118        iomodel->Constant(&control_analysis,InversionIscontrolEnum);
     119        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
     120        iomodel->Constant(&materials_type,MaterialsEnum);
     121
     122        /*return if no processing required*/
     123        if(!isSSA & !isL1L2 & !isHO & !isFS) return;
     124
     125        /*Fetch data needed and allocate vectors: */
     126        iomodel->FetchData(1,FlowequationElementEquationEnum);
     127        finiteelement_list=xNewZeroInit<int>(iomodel->numberofelements);
     128
     129        /*Do we have coupling*/
     130        if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
     131         iscoupling = true;
     132        else
     133         iscoupling = false;
     134
     135        /*Get finite element type*/
     136        if(!iscoupling){
     137                if(isSSA)       iomodel->Constant(&finiteelement,FlowequationFeSSAEnum);
     138                else if(isL1L2) finiteelement = P1Enum;
     139                else if(isHO)   iomodel->Constant(&finiteelement,FlowequationFeHOEnum);
     140                else if(isFS)   iomodel->Constant(&finiteelement,FlowequationFeFSEnum);
     141                for(int i=0;i<iomodel->numberofelements;i++){
     142                        finiteelement_list[i]=finiteelement;
     143                }
     144        }
     145        else{
     146                if(isFS){
     147                        for(int i=0;i<iomodel->numberofelements;i++){
     148                                approximation=reCast<int>(iomodel->Data(FlowequationElementEquationEnum)[i]);
     149                                if(approximation==FSApproximationEnum || approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
     150                                        finiteelement_list[i]=MINIcondensedEnum;
     151                                }
     152                                else{
     153                                        finiteelement_list[i]=P1Enum;
     154                                }
     155                        }
     156                }
     157                else{
     158                        finiteelement = P1Enum;
     159                        for(int i=0;i<iomodel->numberofelements;i++){
     160                                finiteelement_list[i]=finiteelement;
     161                        }
     162                }
     163        }
     164
     165        /*Update elements: */
     166        int counter=0;
     167        for(int i=0;i<iomodel->numberofelements;i++){
     168                if(iomodel->my_elements[i]){
     169                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     170                        element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement_list[i]);
     171                        counter++;
     172                }
     173        }
     174
     175        /*Create inputs: */
     176        iomodel->FetchDataToInput(elements,ThicknessEnum);
     177        iomodel->FetchDataToInput(elements,SurfaceEnum);
     178        iomodel->FetchDataToInput(elements,BedEnum);
     179        iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
     180        iomodel->FetchDataToInput(elements,FrictionPEnum);
     181        iomodel->FetchDataToInput(elements,FrictionQEnum);
     182        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     183        iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
     184        iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
     185        iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
     186        iomodel->FetchDataToInput(elements,VxEnum,0.);
     187        if(dakota_analysis)elements->InputDuplicate(VxEnum,QmuVxEnum);
     188        iomodel->FetchDataToInput(elements,VyEnum,0.);
     189        if(dakota_analysis)elements->InputDuplicate(VyEnum,QmuVyEnum);
     190        iomodel->FetchDataToInput(elements,LoadingforceXEnum);
     191        iomodel->FetchDataToInput(elements,LoadingforceYEnum);
     192        iomodel->FetchDataToInput(elements,DamageDEnum);
     193
     194        if(iomodel->meshtype==Mesh3DEnum){
     195                iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     196                iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     197                iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
     198                iomodel->FetchDataToInput(elements,FlowequationBorderFSEnum);
     199                iomodel->FetchDataToInput(elements,LoadingforceZEnum);
     200                iomodel->FetchDataToInput(elements,VzEnum,0.);
     201                if(dakota_analysis)elements->InputDuplicate(VzEnum,QmuVzEnum);
     202        }
     203        if(isFS){
     204                iomodel->FetchDataToInput(elements,MeshVertexonbedEnum);
     205                iomodel->FetchDataToInput(elements,PressureEnum,0.);
     206                if(dakota_analysis)elements->InputDuplicate(PressureEnum,QmuPressureEnum);
     207        }
     208
     209#ifdef _HAVE_ANDROID_
     210        elements->InputDuplicate(FrictionCoefficientEnum,AndroidFrictionCoefficientEnum);
     211#endif
     212
     213        /*Free data: */
     214        iomodel->DeleteData(1,FlowequationElementEquationEnum);
     215        xDelete<int>(finiteelement_list);
     216}/*}}}*/
     217void StressbalanceAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     218
     219        /*Intermediary*/
     220        bool isSSA,isL1L2,isHO,isFS,iscoupling;
     221        int  finiteelement=-1,approximation=-1;
     222
     223        /*Fetch parameters: */
     224        iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
     225        iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
     226        iomodel->Constant(&isHO,FlowequationIsHOEnum);
     227        iomodel->Constant(&isFS,FlowequationIsFSEnum);
     228
     229        /*Now, check that we have non SIA elements */
     230        if(!isSSA & !isL1L2 & !isHO & !isFS) return;
     231
     232        /*Do we have coupling*/
     233        if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
     234         iscoupling = true;
     235        else
     236         iscoupling = false;
     237
     238        /*If no coupling, call Regular CreateNodes, else, use P1 elements only*/
     239        if(!iscoupling){
     240
     241                /*Get finite element type*/
     242                if(isSSA){
     243                        approximation=SSAApproximationEnum;
     244                        iomodel->Constant(&finiteelement,FlowequationFeSSAEnum);
     245                }
     246                else if(isL1L2){
     247                        approximation = L1L2ApproximationEnum;
     248                        finiteelement = P1Enum;
     249                }
     250                else if(isHO){
     251                        approximation = HOApproximationEnum;
     252                        iomodel->Constant(&finiteelement,FlowequationFeHOEnum);
     253                }
     254                else if(isFS){
     255                        approximation = FSApproximationEnum;
     256                        iomodel->Constant(&finiteelement,FlowequationFeFSEnum);
     257                }
     258                iomodel->FetchData(3,FlowequationBorderSSAEnum,FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
     259                if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderFSEnum);
     260                ::CreateNodes(pnodes,iomodel,StressbalanceAnalysisEnum,finiteelement,approximation);
     261                iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
     262                                        FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
     263        }
     264        else{
     265                /*Coupling: we are going to create P1 Elements only*/
     266
     267                /*Recover nodes*/
     268                Nodes* nodes = *pnodes;
     269                Node*  node  = NULL;
     270                int    lid=0;
     271                if(!nodes) nodes = new Nodes();
     272
     273                iomodel->FetchData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
     274                                        FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
     275                if(isFS){
     276                        /*P1+ velocity*/
     277                        for(int i=0;i<iomodel->numberofvertices;i++){
     278                                if(iomodel->my_vertices[i]){
     279                                        approximation=reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i]);
     280                                        if(approximation==FSApproximationEnum)  approximation=FSvelocityEnum;
     281                                        nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,approximation));
     282                                }
     283                        }
     284                        for(int i=0;i<iomodel->numberofelements;i++){
     285                                if(iomodel->my_elements[i]){
     286                                        node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,StressbalanceAnalysisEnum,FSvelocityEnum);
     287                                        node->Deactivate();
     288                                        nodes->AddObject(node);
     289                                }
     290                        }
     291                        /*P1 pressure*/
     292                        for(int i=0;i<iomodel->numberofvertices;i++){
     293                                if(iomodel->my_vertices[i]){
     294                                        approximation=reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i]);
     295                                        node = new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,i,iomodel,StressbalanceAnalysisEnum,FSpressureEnum);
     296                                        if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum){
     297                                                node->Deactivate();
     298                                        }
     299                                        nodes->AddObject(node);
     300                                }
     301                        }
     302                }
     303                else{
     304                        for(int i=0;i<iomodel->numberofvertices;i++){
     305                                if(iomodel->my_vertices[i]){
     306                                        nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i])));
     307                                }
     308                        }
     309                }
     310                iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
     311                                        FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
     312
     313                /*Assign output pointer: */
     314                *pnodes=nodes;
     315        }
     316}/*}}}*/
     317void StressbalanceAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     318
     319        /*Intermediary*/
     320        int        i,j;
     321        int        count,finiteelement;
     322        IssmDouble g;
     323        IssmDouble rho_ice;
     324        IssmDouble FSreconditioning;
     325        bool       isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
     326        bool       spcpresent = false;
     327        int        Mx,Nx;
     328        int        My,Ny;
     329        int        Mz,Nz;
     330        IssmDouble *spcvx          = NULL;
     331        IssmDouble *spcvy          = NULL;
     332        IssmDouble *spcvz          = NULL;
     333        IssmDouble *nodeonSSA = NULL;
     334        IssmDouble *nodeonHO   = NULL;
     335        IssmDouble *nodeonFS   = NULL;
     336        IssmDouble *nodeonbed      = NULL;
     337        IssmDouble *groundedice_ls = NULL;
     338        IssmDouble *vertices_type  = NULL;
     339        IssmDouble *surface        = NULL;
     340        IssmDouble *z              = NULL;
     341        IssmDouble *timesx=NULL;
     342        IssmDouble *timesy=NULL;
     343        IssmDouble *timesz=NULL;
     344   IssmDouble* values=NULL;
     345
     346        /*Output*/
     347        Constraints *constraints      = NULL;
     348
     349        /*Fetch parameters: */
     350        iomodel->Constant(&g,ConstantsGEnum);
     351        iomodel->Constant(&rho_ice,MaterialsRhoIceEnum);
     352        iomodel->Constant(&FSreconditioning,StressbalanceFSreconditioningEnum);
     353        iomodel->Constant(&isSIA,FlowequationIsSIAEnum);
     354        iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
     355        iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
     356        iomodel->Constant(&isHO,FlowequationIsHOEnum);
     357        iomodel->Constant(&isFS,FlowequationIsFSEnum);
     358
     359        /*Recover pointer: */
     360        constraints=*pconstraints;
     361
     362        /*Now, is the flag macayaealHO on? otherwise, do nothing: */
     363        if(!isSSA && !isHO && !isFS && !isL1L2){
     364                *pconstraints=constraints;
     365                return;
     366        }
     367
     368        /*Do we have coupling*/
     369        if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
     370         iscoupling = true;
     371        else
     372         iscoupling = false;
     373
     374        /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/
     375        if(!iscoupling){
     376
     377                /*Get finite element type*/
     378                if(isSSA)       iomodel->Constant(&finiteelement,FlowequationFeSSAEnum);
     379                else if(isL1L2) finiteelement = P1Enum;
     380                else if(isHO)   iomodel->Constant(&finiteelement,FlowequationFeHOEnum);
     381                else if(isFS){  iomodel->Constant(&finiteelement,FlowequationFeFSEnum);
     382                        /*Deduce velocity interpolation from finite element*/
     383                        switch(finiteelement){
     384                                case P1P1Enum          : finiteelement = P1Enum;       break;
     385                                case P1P1GLSEnum       : finiteelement = P1Enum;       break;
     386                                case MINIcondensedEnum : finiteelement = P1bubbleEnum; break;
     387                                case MINIEnum          : finiteelement = P1bubbleEnum; break;
     388                                case TaylorHoodEnum    : finiteelement = P2Enum;       break;
     389                                default: _error_("finite element "<<finiteelement<<" not supported");
     390                        }
     391                }
     392                else{
     393                        _error_("model not supported yet");
     394                }
     395
     396                if(isFS){
     397
     398                        /*Constraint at the bedrock interface (v.n = vz = 0) (Coordinates will be updated according to the bed slope)*/
     399                        iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
     400                        iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum);
     401                        iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
     402                        iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum);
     403                        if(iomodel->meshtype==Mesh3DEnum){
     404                                iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum);
     405                        }
     406                        else if (iomodel->meshtype==Mesh2DverticalEnum){
     407                                iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvyEnum);
     408                        }
     409                        else{
     410                                _error_("not supported yet");
     411                        }
     412                        for(i=0;i<iomodel->numberofvertices;i++){
     413                                if(iomodel->my_vertices[i]){
     414                                        if(nodeonbed[i]>0. && groundedice_ls[i]>0. && nodeonFS[i]>0.){
     415                                                if(vertices_type[i] == FSApproximationEnum){
     416                                                        for(j=0;j<Nz;j++) spcvz[i*Nz+j] = 0.;
     417                                                }
     418                                                else{
     419                                                        _error_("not supported");
     420                                                }
     421                                        }
     422                                }
     423                        }
     424                        if(iomodel->meshtype==Mesh3DEnum){
     425                                IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1);
     426                                IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,2);
     427                                IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,3);
     428                                iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum);
     429                        }
     430                        else if (iomodel->meshtype==Mesh2DverticalEnum){
     431                                IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1);
     432                                IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,2);
     433                                iomodel->DeleteData(spcvz,StressbalanceSpcvyEnum);
     434                        }
     435                        else{
     436                                _error_("not supported yet");
     437                        }
     438                        iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
     439                        iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum);
     440                        iomodel->DeleteData(nodeonbed,MeshVertexonbedEnum);
     441                        iomodel->DeleteData(groundedice_ls,MaskGroundediceLevelsetEnum);
     442
     443                        /*Pressure spc*/
     444                        count = constraints->Size();
     445                        iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
     446                        iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum);
     447                        iomodel->FetchData(&z,NULL,NULL,MeshZEnum);
     448                        switch(finiteelement){
     449                                case P1bubbleEnum:
     450                                        for(i=0;i<iomodel->numberofvertices;i++){
     451                                                if(iomodel->my_vertices[i]){
     452                                                        if(reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
     453                                                                constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,1,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
     454                                                                count++;
     455                                                        }
     456                                                }
     457                                        }
     458                                        break;
     459                                case P2Enum:
     460                                        for(i=0;i<iomodel->numberofvertices;i++){
     461                                                if(iomodel->my_vertices[i]){
     462                                                        if(reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
     463                                                                constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+i+1,1,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
     464                                                                count++;
     465                                                        }
     466                                                }
     467                                        }
     468                                        break;
     469                                default:
     470                                        _error_("not implemented yet");
     471                        }
     472                        iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
     473                        iomodel->DeleteData(surface,SurfaceEnum);
     474                        iomodel->DeleteData(z,MeshZEnum);
     475                }
     476                else{
     477                        IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1);
     478                        IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,2);
     479                }
     480
     481                *pconstraints=constraints;
     482                return;
     483        }
     484
     485        /*Constraints: fetch data: */
     486        iomodel->FetchData(&spcvx,&Mx,&Nx,StressbalanceSpcvxEnum);
     487        iomodel->FetchData(&spcvy,&My,&Ny,StressbalanceSpcvyEnum);
     488        iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum);
     489        iomodel->FetchData(&nodeonSSA,NULL,NULL,FlowequationBorderSSAEnum);
     490        if(iomodel->meshtype==Mesh3DEnum)iomodel->FetchData(&nodeonHO,NULL,NULL,FlowequationBorderHOEnum);
     491        if(iomodel->meshtype==Mesh3DEnum)iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum);
     492        if(iomodel->meshtype==Mesh3DEnum)iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
     493        if(iomodel->meshtype==Mesh3DEnum)iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum);
     494        iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
     495        iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum);
     496        iomodel->FetchData(&z,NULL,NULL,MeshZEnum);
     497
     498        /*Initialize counter: */
     499        count=0;
     500
     501        /*figure out times: */
     502        timesx=xNew<IssmDouble>(Nx);
     503        for(j=0;j<Nx;j++){
     504                timesx[j]=spcvx[(Mx-1)*Nx+j];
     505        }
     506        /*figure out times: */
     507        timesy=xNew<IssmDouble>(Ny);
     508        for(j=0;j<Ny;j++){
     509                timesy[j]=spcvy[(My-1)*Ny+j];
     510        }
     511        /*figure out times: */
     512        timesz=xNew<IssmDouble>(Nz);
     513        for(j=0;j<Nz;j++){
     514                timesz[j]=spcvz[(Mz-1)*Nz+j];
     515        }
     516
     517        /*Create spcs from x,y,z, as well as the spc values on those spcs: */
     518        for(i=0;i<iomodel->numberofvertices;i++){
     519                if(iomodel->my_vertices[i]){
     520
     521                        /*Start with adding spcs of coupling: zero at the border SSA/HO for the appropriate dofs*/
     522                        if(reCast<int,IssmDouble>(vertices_type[i]==SSAHOApproximationEnum)){
     523                                /*If grionSSA, spc HO dofs: 3 & 4*/
     524                                        if (reCast<int,IssmDouble>(nodeonHO[i])){
     525                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     526                                                count++;
     527                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     528                                                count++;
     529                                                if (!xIsNan<IssmDouble>(spcvx[i])){
     530                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     531                                                        count++;
     532                                                }
     533                                                if (!xIsNan<IssmDouble>(spcvy[i])){
     534                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     535                                                        count++;
     536                                                }
     537
     538                                        }
     539                                        else if (reCast<int,IssmDouble>(nodeonSSA[i])){
     540                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     541                                                count++;
     542                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     543                                                count++;
     544                                                if (!xIsNan<IssmDouble>(spcvx[i])){
     545                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     546                                                        count++;
     547                                                }
     548                                                if (!xIsNan<IssmDouble>(spcvy[i])){
     549                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     550                                                        count++;
     551                                                }
     552
     553                                        }
     554                                        else _error_("if vertices_type is SSAHO, you shoud have nodeonHO or nodeonSSA");
     555                        }
     556                        /*Also add spcs of coupling: zero at the border HO/FS for the appropriate dofs*/
     557                        else if (reCast<int,IssmDouble>(vertices_type[i])==HOFSApproximationEnum){
     558                                /*If grion,HO spc FS dofs: 3 4 & 5*/
     559                                        if (reCast<int,IssmDouble>(nodeonHO[i])){
     560                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     561                                                count++;
     562                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     563                                                count++;
     564                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     565                                                count++;
     566                                                if (!xIsNan<IssmDouble>(spcvx[i])){
     567                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     568                                                        count++;
     569                                                }
     570                                                if (!xIsNan<IssmDouble>(spcvy[i])){
     571                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     572                                                        count++;
     573                                                }
     574
     575                                        }
     576                                        else if (reCast<int,IssmDouble>(nodeonFS[i])){ //spc HO nodes: 1 & 2
     577                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     578                                                count++;
     579                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     580                                                count++;
     581                                                if (!xIsNan<IssmDouble>(spcvx[i])){
     582                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     583                                                        count++;
     584                                                }
     585                                                if (!xIsNan<IssmDouble>(spcvy[i])){
     586                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     587                                                        count++;
     588                                                }
     589                                                if (!xIsNan<IssmDouble>(spcvz[i])){
     590                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     591                                                        count++;
     592                                                }
     593                                        }
     594                                        else _error_("if vertices_type is HOFS, you shoud have nodeonHO or nodeonFS");
     595                        }
     596                        /*Also add spcs of coupling: zero at the border HO/FS for the appropriate dofs*/
     597                        else if (reCast<int,IssmDouble>(vertices_type[i])==SSAFSApproximationEnum){
     598                                /*If grion,HO spc FS dofs: 3 4 & 5*/
     599                                        if (reCast<int,IssmDouble>(nodeonSSA[i])){
     600                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     601                                                count++;
     602                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     603                                                count++;
     604                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     605                                                count++;
     606                                                if (!xIsNan<IssmDouble>(spcvx[i])){
     607                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     608                                                        count++;
     609                                                }
     610                                                if (!xIsNan<IssmDouble>(spcvy[i])){
     611                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     612                                                        count++;
     613                                                }
     614
     615                                        }
     616                                        else if (reCast<int,IssmDouble>(nodeonFS[i])){ //spc SSA nodes: 1 & 2
     617                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     618                                                count++;
     619                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     620                                                count++;
     621                                                if (!xIsNan<IssmDouble>(spcvx[i])){
     622                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     623                                                        count++;
     624                                                }
     625                                                if (!xIsNan<IssmDouble>(spcvy[i])){
     626                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     627                                                        count++;
     628                                                }
     629                                                if (!xIsNan<IssmDouble>(spcvz[i])){
     630                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     631                                                        count++;
     632                                                }
     633                                        }
     634                                        else _error_("if vertices_type is SSAFS, you shoud have nodeonSSA or nodeonFS");
     635                        }
     636                        /*Now add the regular spcs*/
     637                        else{
     638                                if (Mx==iomodel->numberofvertices && !xIsNan<IssmDouble>(spcvx[i])){
     639                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     640                                        count++;
     641
     642                                }
     643                                else if (Mx==iomodel->numberofvertices+1) {
     644                                        /*figure out times and values: */
     645                                        values=xNew<IssmDouble>(Nx);
     646                                        spcpresent=false;
     647                                        for(j=0;j<Nx;j++){
     648                                                values[j]=spcvx[i*Nx+j];
     649                                                if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
     650                                        }
     651
     652                                        if(spcpresent){
     653                                                constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,Nx,timesx,values,StressbalanceAnalysisEnum));
     654                                                count++;
     655                                        }
     656                                        xDelete<IssmDouble>(values);
     657                                }
     658                                else if (vertices_type[i]==SIAApproximationEnum){
     659                                        constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,StressbalanceAnalysisEnum));
     660                                        count++;
     661                                }
     662
     663                                if (My==iomodel->numberofvertices && !xIsNan<IssmDouble>(spcvy[i])){
     664                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vy.
     665                                        count++;
     666                                }
     667                                else if (My==iomodel->numberofvertices+1){
     668                                        /*figure out times and values: */
     669                                        values=xNew<IssmDouble>(Ny);
     670                                        spcpresent=false;
     671                                        for(j=0;j<Ny;j++){
     672                                                values[j]=spcvy[i*Ny+j];
     673                                                if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
     674                                        }
     675                                        if(spcpresent){
     676                                                constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,Ny,timesy,values,StressbalanceAnalysisEnum));
     677                                                count++;
     678                                        }
     679                                        xDelete<IssmDouble>(values);
     680                                }
     681                                else if (vertices_type[i]==SIAApproximationEnum){
     682                                        constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,StressbalanceAnalysisEnum));
     683                                        count++;
     684                                }
     685
     686                                if (reCast<int,IssmDouble>(vertices_type[i])==FSApproximationEnum ||  (reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum)){
     687                                        if (Mz==iomodel->numberofvertices && !xIsNan<IssmDouble>(spcvz[i])){
     688                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
     689                                                count++;
     690                                        }
     691                                        else if (Mz==iomodel->numberofvertices+1){
     692                                                /*figure out times and values: */
     693                                                values=xNew<IssmDouble>(Nz);
     694                                                spcpresent=false;
     695                                                for(j=0;j<Nz;j++){
     696                                                        values[j]=spcvz[i*Nz+j];
     697                                                        if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
     698                                                }
     699                                                if(spcpresent){
     700                                                        constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,Nz,timesz,values,StressbalanceAnalysisEnum));
     701                                                        count++;
     702                                                }
     703                                                xDelete<IssmDouble>(values);
     704                                        }
     705
     706                                }
     707                                if (reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
     708                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,1,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
     709                                        count++;
     710                                }
     711                        }
     712
     713                        /*Constraint at the bedrock interface (v.n = vz = 0) (Coordinates will be updated according to the bed slope)*/
     714                        if (iomodel->meshtype==Mesh3DEnum) if(nodeonbed[i]>0. && groundedice_ls[i]>=0. && nodeonFS[i]>0.){
     715                                 switch(reCast<int,IssmDouble>(vertices_type[i])){
     716                                        case SSAFSApproximationEnum:
     717                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,0.,StressbalanceAnalysisEnum));
     718                                                count++;
     719                                                break;
     720                                        case HOFSApproximationEnum:
     721                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,0.,StressbalanceAnalysisEnum));
     722                                                count++;
     723                                                break;
     724                                        case FSApproximationEnum:
     725                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0.,StressbalanceAnalysisEnum));
     726                                                count++;
     727                                                break;
     728                                        default: _error_("Vertex approximation " << EnumToStringx(reCast<int,IssmDouble>(vertices_type[i])) << " not supported");
     729                                }
     730                        }
     731                }
     732        }
     733
     734        /*Free data: */
     735        iomodel->DeleteData(spcvx,StressbalanceSpcvxEnum);
     736        iomodel->DeleteData(spcvy,StressbalanceSpcvyEnum);
     737        iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum);
     738        iomodel->DeleteData(nodeonSSA,FlowequationBorderSSAEnum);
     739        if(iomodel->meshtype==Mesh3DEnum)iomodel->DeleteData(nodeonHO,FlowequationBorderHOEnum);
     740        if(iomodel->meshtype==Mesh3DEnum)iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum);
     741        if(iomodel->meshtype==Mesh3DEnum)iomodel->DeleteData(nodeonbed,MeshVertexonbedEnum);
     742        if(iomodel->meshtype==Mesh3DEnum)iomodel->DeleteData(groundedice_ls,MaskGroundediceLevelsetEnum);
     743        iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
     744        iomodel->DeleteData(surface,SurfaceEnum);
     745        iomodel->DeleteData(z,MeshZEnum);
     746
     747        /*Free resources:*/
     748        xDelete<IssmDouble>(timesx);
     749        xDelete<IssmDouble>(timesy);
     750        xDelete<IssmDouble>(timesz);
     751        xDelete<IssmDouble>(values);
     752
     753        /*Assign output pointer: */
     754        *pconstraints=constraints;
     755}/*}}}*/
     756void StressbalanceAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     757
     758        /*Intermediary*/
     759        const int   RIFTINFOSIZE = 12;
     760        int         i;
     761        int         count;
     762        int         penpair_ids[2];
     763        bool        isSSA,isL1L2,isHO,isFS;
     764        int         numpenalties,numrifts,numriftsegments;
     765        IssmDouble *riftinfo       = NULL;
     766        IssmDouble *penalties      = NULL;
     767        int         assert_int;
     768
     769        /*Fetch parameters: */
     770        iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
     771        iomodel->Constant(&isFS,FlowequationIsFSEnum);
     772        iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
     773        iomodel->Constant(&isHO,FlowequationIsHOEnum);
     774        iomodel->Constant(&numrifts,RiftsNumriftsEnum);
     775
     776        /*Recover pointer: */
     777        Loads* loads=*ploads;
     778
     779        /*Now, is the flag macayaealHO on? otherwise, do nothing: */
     780        if(!isSSA && !isHO && !isFS && !isL1L2) return;
     781
     782        /*Initialize counter: */
     783        count=0;
     784
     785        /*Create Penpair for penalties: */
     786        iomodel->FetchData(&penalties,&numpenalties,NULL,StressbalanceVertexPairingEnum);
     787
     788        for(i=0;i<numpenalties;i++){
     789
     790                if(iomodel->my_vertices[reCast<int,IssmDouble>(penalties[2*i+0]-1)]){
     791
     792                        /*In debugging mode, check that the second node is in the same cpu*/
     793                        assert_int=iomodel->my_vertices[reCast<int,IssmDouble>(penalties[2*i+1]-1)]; _assert_(assert_int);
     794
     795                        /*Get node ids*/
     796                        penpair_ids[0]=iomodel->nodecounter+reCast<int,IssmDouble>(penalties[2*i+0]);
     797                        penpair_ids[1]=iomodel->nodecounter+reCast<int,IssmDouble>(penalties[2*i+1]);
     798
     799                        /*Create Load*/
     800                        loads->AddObject(new Penpair(iomodel->loadcounter+count+1,&penpair_ids[0],StressbalanceAnalysisEnum));
     801                        count++;
     802                }
     803        }
     804
     805        /*free ressources: */
     806        iomodel->DeleteData(penalties,StressbalanceVertexPairingEnum);
     807
     808        /*Create Riffront loads for rifts: */
     809#ifdef _HAVE_RIFTS_
     810        if(numrifts){
     811                iomodel->FetchData(&riftinfo,&numriftsegments,NULL,RiftsRiftstructEnum);
     812                iomodel->FetchData(5,RiftsRiftstructEnum,ThicknessEnum,BedEnum,SurfaceEnum,MaskGroundediceLevelsetEnum);
     813                for(i=0;i<numriftsegments;i++){
     814                        if(iomodel->my_elements[reCast<int,IssmDouble>(*(riftinfo+RIFTINFOSIZE*i+2))-1]){
     815                                loads->AddObject(new Riftfront(iomodel->loadcounter+count+1,i,iomodel,StressbalanceAnalysisEnum));
     816                                count++;
     817                        }
     818                }
     819                iomodel->DeleteData(5,RiftsRiftstructEnum,ThicknessEnum,BedEnum,SurfaceEnum,MaskGroundediceLevelsetEnum);
     820                xDelete<IssmDouble>(riftinfo);
     821        }
     822#endif
     823
     824        /*Assign output pointer: */
     825        *ploads=loads;
     826}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int StressbalanceSIAAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  StressbalanceSIAAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 2;
    1010}/*}}}*/
     11void StressbalanceSIAAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12
     13        /*No specific parameters*/
     14
     15}/*}}}*/
     16void StressbalanceSIAAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     17
     18        /*Fetch data needed: */
     19        bool   isSIA;
     20        iomodel->Constant(&isSIA,FlowequationIsSIAEnum);
     21
     22        /*Now, is the flag SIA on? otherwise, do nothing: */
     23        if (!isSIA)return;
     24
     25        iomodel->FetchData(1,FlowequationElementEquationEnum);
     26
     27        /*Update elements: */
     28        int counter=0;
     29        for(int i=0;i<iomodel->numberofelements;i++){
     30                if(iomodel->my_elements[i]){
     31                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     32                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     33                        counter++;
     34                }
     35        }
     36
     37        iomodel->FetchDataToInput(elements,ThicknessEnum);
     38
     39        /*Free data: */
     40        iomodel->DeleteData(1,FlowequationElementEquationEnum);
     41}/*}}}*/
     42void StressbalanceSIAAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     43
     44        /*Intermediaries*/
     45        bool  isSIA;
     46        Node* node = NULL;
     47
     48        /*Fetch parameters: */
     49        iomodel->Constant(&isSIA,FlowequationIsSIAEnum);
     50
     51        /*Now, is the flag isSIA on? otherwise, do nothing: */
     52        if(!isSIA) return;
     53
     54        /*First create nodes*/
     55        Nodes* nodes=*pnodes;
     56        int    lid=0;
     57        if(!nodes) nodes = new Nodes();
     58
     59        iomodel->FetchData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
     60                                FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
     61
     62        for(int i=0;i<iomodel->numberofvertices;i++){
     63                if(iomodel->my_vertices[i]){
     64
     65                        /*Create new node if is in this processor's partition*/
     66                        node = new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceSIAAnalysisEnum,reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i]));
     67
     68                        /*Deactivate node if not SIA*/
     69                        if(reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i])!=SIAApproximationEnum){
     70                                node->Deactivate();
     71                        }
     72
     73                        /*Add to Nodes dataset*/
     74                        nodes->AddObject(node);
     75                }
     76        }
     77
     78        iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
     79                                FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
     80
     81        /*Assign output pointer: */
     82        *pnodes=nodes;
     83
     84}/*}}}*/
     85void StressbalanceSIAAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     86
     87        /*Intermediary*/
     88        int        count;
     89        IssmDouble yts;
     90        bool       isSIA;
     91
     92        /*Output*/
     93        Constraints* constraints = NULL;
     94
     95        /*Recover pointer: */
     96        constraints=*pconstraints;
     97
     98        /*Fetch parameters: */
     99        iomodel->Constant(&yts,ConstantsYtsEnum);
     100        iomodel->Constant(&isSIA,FlowequationIsSIAEnum);
     101
     102        /*Now, is the flag isSIA on? otherwise, do nothing: */
     103        if (!isSIA) return;
     104
     105        /*Fetch data: */
     106        iomodel->FetchData(3,StressbalanceSpcvxEnum,StressbalanceSpcvyEnum,FlowequationVertexEquationEnum);
     107
     108        /*Initialize conunter*/
     109        count=0;
     110
     111        /*vx and vy are spc'd if we are not on nodeonSIA: */
     112        for(int i=0;i<iomodel->numberofvertices;i++){
     113                /*keep only this partition's nodes:*/
     114                if((iomodel->my_vertices[i])){
     115                        if (reCast<int,IssmDouble>(iomodel->Data(FlowequationVertexEquationEnum)[i])!=SIAApproximationEnum){
     116
     117                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceSIAAnalysisEnum));
     118                                count++;
     119
     120                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceSIAAnalysisEnum));
     121                                count++;
     122                        }
     123                        else{
     124                                if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvxEnum)[i])){
     125                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,iomodel->Data(StressbalanceSpcvxEnum)[i]/yts,StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     126                                        count++;
     127                                }
     128
     129                                if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvyEnum)[i])){
     130                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,iomodel->Data(StressbalanceSpcvyEnum)[i]/yts,StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
     131                                        count++;
     132                                }
     133                        }
     134                }
     135        }
     136
     137        /*Free data: */
     138        iomodel->DeleteData(3,StressbalanceSpcvxEnum,StressbalanceSpcvyEnum,FlowequationVertexEquationEnum);
     139
     140        /*Assign output pointer: */
     141        *pconstraints=constraints;
     142}/*}}}*/
     143void StressbalanceSIAAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     144
     145        /*No loads*/
     146
     147}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int StressbalanceVerticalAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  StressbalanceVerticalAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void StressbalanceVerticalAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12
     13        /*No specific parameters*/
     14
     15}/*}}}*/
     16void StressbalanceVerticalAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     17
     18        /*return if not 3d mesh*/
     19        if(iomodel->meshtype!=Mesh3DEnum) return;
     20
     21        /*Update elements: */
     22        int counter=0;
     23        for(int i=0;i<iomodel->numberofelements;i++){
     24                if(iomodel->my_elements[i]){
     25                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     26                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     27                        counter++;
     28                }
     29        }
     30
     31        iomodel->FetchDataToInput(elements,ThicknessEnum);
     32        iomodel->FetchDataToInput(elements,SurfaceEnum);
     33        iomodel->FetchDataToInput(elements,BedEnum);
     34        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     35        iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     36        iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     37        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
     38        iomodel->FetchDataToInput(elements,VxEnum,0.);
     39        iomodel->FetchDataToInput(elements,VyEnum,0.);
     40}/*}}}*/
     41void StressbalanceVerticalAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     42
     43        /*return if not 3d mesh*/
     44        if(iomodel->meshtype!=Mesh3DEnum) return;
     45
     46        iomodel->FetchData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationVertexEquationEnum);
     47        ::CreateNodes(pnodes,iomodel,StressbalanceVerticalAnalysisEnum,P1Enum);
     48        iomodel->DeleteData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationVertexEquationEnum);
     49}/*}}}*/
     50void StressbalanceVerticalAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     51
     52        /*Intermediary*/
     53        int count;
     54        IssmDouble yts;
     55
     56        /*Fetch parameters: */
     57        iomodel->Constant(&yts,ConstantsYtsEnum);
     58
     59        /*Recover pointer: */
     60        Constraints* constraints=*pconstraints;
     61
     62        /*return if not 3d mesh*/
     63        if(iomodel->meshtype!=Mesh3DEnum) return;
     64
     65        /*Fetch data: */
     66        iomodel->FetchData(2,StressbalanceSpcvzEnum,FlowequationBorderFSEnum);
     67
     68        /*Initialize counter*/
     69        count=0;
     70
     71        /*Create spcs from x,y,z, as well as the spc values on those spcs: */
     72        for(int i=0;i<iomodel->numberofvertices;i++){
     73
     74                /*keep only this partition's nodes:*/
     75                if(iomodel->my_vertices[i]){
     76
     77                        if (reCast<int,IssmDouble>(iomodel->Data(FlowequationBorderFSEnum)[i])){
     78                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceVerticalAnalysisEnum)); //spc to zero as vertical velocity is done in Horiz for FS
     79                                count++;
     80                        }
     81                        else if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvzEnum)[i])){
     82                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,
     83                                                                iomodel->Data(StressbalanceSpcvzEnum)[i]/yts,StressbalanceVerticalAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     84                                count++;
     85
     86                        }
     87                }
     88        }
     89
     90        /*Free data: */
     91        iomodel->DeleteData(2,StressbalanceSpcvzEnum,FlowequationBorderFSEnum);
     92
     93        /*Assign output pointer: */
     94        *pconstraints=constraints;
     95}/*}}}*/
     96void StressbalanceVerticalAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     97
     98        /*No loads*/
     99
     100}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

  • issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp

    r16534 r16539  
    66
    77/*Model processing*/
    8 int ThermalAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
     8int  ThermalAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    99        return 1;
    1010}/*}}}*/
     11void ThermalAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12}/*}}}*/
     13void ThermalAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     14
     15        /*Now, is the model 3d? otherwise, do nothing: */
     16        if(iomodel->meshtype==Mesh2DhorizontalEnum)return;
     17
     18        /*Update elements: */
     19        int counter=0;
     20        for(int i=0;i<iomodel->numberofelements;i++){
     21                if(iomodel->my_elements[i]){
     22                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     23                        element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
     24                        counter++;
     25                }
     26        }
     27
     28        bool dakota_analysis;
     29        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
     30
     31        iomodel->FetchDataToInput(elements,ThicknessEnum);
     32        iomodel->FetchDataToInput(elements,SurfaceEnum);
     33        iomodel->FetchDataToInput(elements,BedEnum);
     34        iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
     35        iomodel->FetchDataToInput(elements,FrictionPEnum);
     36        iomodel->FetchDataToInput(elements,FrictionQEnum);
     37        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     38        iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
     39        iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     40        iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
     41        iomodel->FetchDataToInput(elements,FlowequationElementEquationEnum);
     42        iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
     43        iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
     44        iomodel->FetchDataToInput(elements,PressureEnum);
     45        iomodel->FetchDataToInput(elements,TemperatureEnum);
     46        iomodel->FetchDataToInput(elements,BasalforcingsGeothermalfluxEnum);
     47        iomodel->FetchDataToInput(elements,VxEnum);
     48        iomodel->FetchDataToInput(elements,VyEnum);
     49        iomodel->FetchDataToInput(elements,VzEnum);
     50        InputUpdateFromConstantx(elements,0.,VxMeshEnum);
     51        InputUpdateFromConstantx(elements,0.,VyMeshEnum);
     52        InputUpdateFromConstantx(elements,0.,VzMeshEnum);
     53        if(dakota_analysis){
     54                elements->InputDuplicate(TemperatureEnum,QmuTemperatureEnum);
     55                elements->InputDuplicate(BasalforcingsMeltingRateEnum,QmuMeltingEnum);
     56                elements->InputDuplicate(VxMeshEnum,QmuVxMeshEnum);
     57                elements->InputDuplicate(VxMeshEnum,QmuVyMeshEnum);
     58                elements->InputDuplicate(VxMeshEnum,QmuVzMeshEnum);
     59        }
     60}/*}}}*/
     61void ThermalAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/
     62
     63        if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     64        ::CreateNodes(pnodes,iomodel,ThermalAnalysisEnum,P1Enum);
     65        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     66}/*}}}*/
     67void ThermalAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/
     68
     69        /*Recover pointer: */
     70        Constraints* constraints=*pconstraints;
     71
     72        /*Only 3d mesh supported*/
     73        if(iomodel->meshtype==Mesh3DEnum){
     74                IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum,P1Enum);
     75        }
     76
     77        /*Assign output pointer: */
     78        *pconstraints=constraints;
     79}/*}}}*/
     80void ThermalAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/
     81
     82        /*Recover pointer: */
     83        Loads* loads=*ploads;
     84
     85        if(iomodel->meshtype==Mesh2DhorizontalEnum) _error_("2d meshes not supported yet");
     86
     87        /*create penalties for nodes: no node can have a temperature over the melting point*/
     88        iomodel->FetchData(1,ThermalSpctemperatureEnum);
     89        CreateSingleNodeToElementConnectivity(iomodel);
     90
     91        for(int i=0;i<iomodel->numberofvertices;i++){
     92
     93                /*keep only this partition's nodes:*/
     94                if(iomodel->my_vertices[i]){
     95                        if (xIsNan<IssmDouble>(iomodel->Data(ThermalSpctemperatureEnum)[i])){ //No penalty applied on spc nodes!
     96                                loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,ThermalAnalysisEnum));
     97                        }
     98                }
     99        }
     100        iomodel->DeleteData(1,ThermalSpctemperatureEnum);
     101
     102        /*Assign output pointer: */
     103        *ploads=loads;
     104}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/ThermalAnalysis.h

    r16534 r16539  
    1212
    1313        public:
    14                 int DofsPerNode(int** doflist,int meshtype,int approximation);
     14                int  DofsPerNode(int** doflist,int meshtype,int approximation);
     15                void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum);
     16                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     17                void CreateNodes(Nodes** pnodes,IoModel* iomodel);
     18                void CreateConstraints(Constraints** pconstraints,IoModel* iomodel);
     19                void CreateLoads(Loads** ploads, IoModel* iomodel);
    1520};
    1621#endif

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

    r16237 r16539  
    3636
    3737        /*data: */
     38        const int RIFTINFOSIZE = 12;
    3839        int    riftfront_node_ids[2];
    3940        int    riftfront_elem_ids[2];

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp

    r16533 r16539  
    1414#include "./ModelProcessorx.h"
    1515
    16 void CreateDataSets(Elements** pelements,Nodes** pnodes, Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads,Parameters** pparameters,IoModel* iomodel,FILE* toolkitfile,char* rootpath,const int solution_type,const int analysis_type,const int nummodels,int analysis_counter){
     16void CreateDataSets(Elements** pelements,Nodes** pnodes, Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads,Parameters** pparameters,IoModel* iomodel,FILE* toolkitfile,char* rootpath,const int solution_enum,const int analysis_enum,const int nummodels,int analysis_counter){
    1717
    1818        Elements   *elements   = NULL;
    … …  
    2525        ElementsAndVerticesPartitioning(&iomodel->my_elements,&iomodel->my_vertices,iomodel);
    2626
    27         /*Create elements, vertices and materials, independent of analysis_type: */
     27        /*Create elements, vertices and materials, independent of analysis_enum: */
    2828        CreateElementsVerticesAndMaterials(pelements, pvertices, pmaterials, iomodel,nummodels);
     29
     30        /*Create Parameters*/
     31        CreateParameters(pparameters,iomodel,rootpath,toolkitfile,solution_enum,analysis_enum,analysis_counter);
    2932
    3033        /*Recover elements and materials, for future update: */
    … …  
    3841
    3942        /*Now, branch onto analysis dependent model generation: */
    40         switch(analysis_type){
    41 
    42                 #ifdef _HAVE_STRESSBALANCE_
    43                 case StressbalanceAnalysisEnum:
    44                         CreateNodesStressbalance(pnodes,iomodel);
    45                         CreateConstraintsStressbalance(pconstraints,iomodel);
    46                         CreateLoadsStressbalance(ploads,iomodel);
    47                         UpdateElementsStressbalance(elements,iomodel,analysis_counter,analysis_type);
    48                         break;
    49 
    50                 case StressbalanceVerticalAnalysisEnum:
    51                         CreateNodesStressbalanceVertical(pnodes, iomodel);
    52                         CreateConstraintsStressbalanceVertical(pconstraints,iomodel);
    53                         CreateLoadsStressbalanceVertical(ploads,iomodel);
    54                         UpdateElementsStressbalanceVertical(elements,iomodel,analysis_counter,analysis_type);
    55                         break;
    56 
    57                 case StressbalanceSIAAnalysisEnum:
    58                         CreateNodesStressbalanceSIA(pnodes, iomodel);
    59                         CreateConstraintsStressbalanceSIA(pconstraints,iomodel);
    60                         CreateLoadsStressbalanceSIA(ploads,iomodel);
    61                         UpdateElementsStressbalanceSIA(elements,iomodel,analysis_counter,analysis_type);
    62                         break;
    63                 #endif
    64 
    65                 #ifdef _HAVE_HYDROLOGY_
    66                 case HydrologyShreveAnalysisEnum:
    67                         CreateNodesHydrologyShreve(pnodes, iomodel);
    68                         CreateConstraintsHydrologyShreve(pconstraints,iomodel);
    69                         CreateLoadsHydrologyShreve(ploads,iomodel);
    70                         UpdateElementsHydrologyShreve(elements,iomodel,analysis_counter,analysis_type);
    71                         break;
    72                 case HydrologyDCInefficientAnalysisEnum:
    73                         CreateNodesHydrologyDCInefficient(pnodes, iomodel);
    74                         CreateConstraintsHydrologyDCInefficient(pconstraints,iomodel);
    75                         CreateLoadsHydrologyDCInefficient(ploads,iomodel);
    76                         UpdateElementsHydrologyDCInefficient(elements,iomodel,analysis_counter,analysis_type);
    77                         break;
    78                 case HydrologyDCEfficientAnalysisEnum:
    79                         CreateNodesHydrologyDCEfficient(pnodes, iomodel);
    80                         CreateConstraintsHydrologyDCEfficient(pconstraints,iomodel);
    81                         CreateLoadsHydrologyDCEfficient(ploads,iomodel);
    82                         UpdateElementsHydrologyDCEfficient(elements,iomodel,analysis_counter,analysis_type);
    83                         break;
    84                 #endif
    85 
    86                 #ifdef _HAVE_THERMAL_
    87                 case ThermalAnalysisEnum:
    88                         CreateNodesThermal(pnodes, iomodel);
    89                         CreateConstraintsThermal(pconstraints,iomodel);
    90                         CreateLoadsThermal(ploads,iomodel);
    91                         UpdateElementsThermal(elements,iomodel,analysis_counter,analysis_type);
    92                         break;
    93 
    94                 case EnthalpyAnalysisEnum:
    95                         CreateNodesEnthalpy(pnodes, iomodel);
    96                         CreateConstraintsEnthalpy(pconstraints,iomodel);
    97                         CreateLoadsEnthalpy(ploads,iomodel);
    98                         UpdateElementsEnthalpy(elements,iomodel,analysis_counter,analysis_type);
    99                         break;
    100 
    101                 case MeltingAnalysisEnum:
    102                         CreateNodesMelting(pnodes, iomodel);
    103                         CreateConstraintsMelting(pconstraints,iomodel);
    104                         CreateLoadsMelting(ploads,iomodel);
    105                         UpdateElementsMelting(elements,iomodel,analysis_counter,analysis_type);
    106                         break;
    107                 #endif
    108 
    109                 #ifdef _HAVE_BALANCED_
    110                 case BalancethicknessAnalysisEnum:
    111                         CreateNodesBalancethickness(pnodes, iomodel);
    112                         CreateConstraintsBalancethickness(pconstraints,iomodel);
    113                         CreateLoadsBalancethickness(ploads,iomodel);
    114                         UpdateElementsBalancethickness(elements,iomodel,analysis_counter,analysis_type);
    115                         break;
    116                 case BalancevelocityAnalysisEnum:
    117                         CreateNodesBalancevelocity(pnodes, iomodel);
    118                         CreateConstraintsBalancevelocity(pconstraints,iomodel);
    119                         CreateLoadsBalancevelocity(ploads,iomodel);
    120                         UpdateElementsBalancevelocity(elements,iomodel,analysis_counter,analysis_type);
    121                         break;
    122                 case SmoothedSurfaceSlopeXAnalysisEnum:
    123                         iomodel->FetchData(1,MeshVertexonbedEnum);
    124                         CreateNodes(pnodes,iomodel,SmoothedSurfaceSlopeXAnalysisEnum,P1Enum);
    125                         iomodel->DeleteData(1,MeshVertexonbedEnum);
    126                         UpdateElementsL2ProjectionBase(elements,iomodel,analysis_counter,analysis_type);
    127                         break;
    128                 case SmoothedSurfaceSlopeYAnalysisEnum:
    129                         iomodel->FetchData(1,MeshVertexonbedEnum);
    130                         CreateNodes(pnodes,iomodel,SmoothedSurfaceSlopeYAnalysisEnum,P1Enum);
    131                         iomodel->DeleteData(1,MeshVertexonbedEnum);
    132                         UpdateElementsL2ProjectionBase(elements,iomodel,analysis_counter,analysis_type);
    133                         break;
    134                 #endif
    135 
    136                 #ifdef _HAVE_GIA_
    137                 case GiaAnalysisEnum:
    138                         CreateNodesGia(pnodes, iomodel);
    139                         CreateConstraintsGia(pconstraints,iomodel);
    140                         CreateLoadsGia(ploads,iomodel);
    141                         UpdateElementsGia(elements,iomodel,analysis_counter,analysis_type);
    142                         break;
    143                 #endif
    144 
    145                 #ifdef _HAVE_DAMAGE_
    146                 case DamageEvolutionAnalysisEnum:
    147                         CreateNodesDamage(pnodes, iomodel);
    148                         CreateConstraintsDamage(pconstraints,iomodel);
    149                         CreateLoadsDamage(ploads,iomodel);
    150                         UpdateElementsDamage(elements,iomodel,analysis_counter,analysis_type);
    151                         break;
    152                 #endif
    153 
    154 
    155                 #ifdef _HAVE_SLOPE_
    156                 case L2ProjectionBaseAnalysisEnum:
    157                         CreateNodesL2ProjectionBase(pnodes, iomodel);
    158                         CreateConstraintsL2ProjectionBase(pconstraints,iomodel);
    159                         CreateLoadsL2ProjectionBase(ploads,iomodel);
    160                         UpdateElementsL2ProjectionBase(elements,iomodel,analysis_counter,analysis_type);
    161                         break;
    162 
    163                 case L2ProjectionTopAnalysisEnum:
    164                         CreateNodesL2ProjectionTop(pnodes, iomodel);
    165                         CreateConstraintsL2ProjectionTop(pconstraints,iomodel);
    166                         CreateLoadsL2ProjectionTop(ploads,iomodel);
    167                         UpdateElementsL2ProjectionTop(elements,iomodel,analysis_counter,analysis_type);
    168                         break;
    169                 #endif
    170 
    171                 #ifdef _HAVE_MASSTRANSPORT_
    172                 case MasstransportAnalysisEnum:
    173                         CreateNodesMasstransport(pnodes, iomodel);
    174                         CreateConstraintsMasstransport(pconstraints,iomodel);
    175                         CreateLoadsMasstransport(ploads,iomodel);
    176                         UpdateElementsMasstransport(elements,iomodel,analysis_counter,analysis_type);
    177                         break;
    178                 case FreeSurfaceTopAnalysisEnum:
    179                         CreateNodesFreeSurfaceTop(pnodes, iomodel);
    180                         CreateConstraintsFreeSurfaceTop(pconstraints,iomodel);
    181                         CreateLoadsFreeSurfaceTop(ploads,iomodel);
    182                         UpdateElementsFreeSurfaceTop(elements,iomodel,analysis_counter,analysis_type);
    183                         break;
    184                 case FreeSurfaceBaseAnalysisEnum:
    185                         CreateNodesFreeSurfaceBase(pnodes, iomodel);
    186                         CreateConstraintsFreeSurfaceBase(pconstraints,iomodel);
    187                         CreateLoadsFreeSurfaceBase(ploads,iomodel);
    188                         UpdateElementsFreeSurfaceBase(elements,iomodel,analysis_counter,analysis_type);
    189                         break;
    190                 case ExtrudeFromBaseAnalysisEnum:
    191                         CreateNodesExtrudeFromBase(pnodes, iomodel);
    192                         CreateConstraintsExtrudeFromBase(pconstraints,iomodel);
    193                         CreateLoadsExtrudeFromBase(ploads,iomodel);
    194                         UpdateElementsExtrudeFromBase(elements,iomodel,analysis_counter,analysis_type);
    195                         break;
    196                 case ExtrudeFromTopAnalysisEnum:
    197                         CreateNodesExtrudeFromTop(pnodes, iomodel);
    198                         CreateConstraintsExtrudeFromTop(pconstraints,iomodel);
    199                         CreateLoadsExtrudeFromTop(ploads,iomodel);
    200                         UpdateElementsExtrudeFromTop(elements,iomodel,analysis_counter,analysis_type);
    201                         break;
    202                 #endif
    203 
    204                 default:
    205                         _error_("analysis_type: " << EnumToStringx(analysis_type) << " not supported yet!");
    206         }
     43        Analysis* analysis = EnumToAnalysis(analysis_enum);
     44        analysis->UpdateParameters(pparameters,iomodel,solution_enum,analysis_enum);
     45        analysis->CreateNodes(pnodes,iomodel);
     46        analysis->CreateConstraints(pconstraints,iomodel);
     47        analysis->CreateLoads(ploads,iomodel);
     48        analysis->UpdateElements(elements,iomodel,analysis_counter,analysis_enum);
     49        delete analysis;
    20750
    20851        /*Update Elements and Materials For Inversions*/
    … …  
    21659        #endif
    21760
    218         /*Generate objects that are not dependent on any analysis_type: */
    219         CreateParameters(pparameters,iomodel,rootpath,toolkitfile,solution_type,analysis_type,analysis_counter);
    22061
    22162        /*Update Elements in case we are running a transient solution: */
    22263        #ifdef _HAVE_TRANSIENT_
    22364        parameters=*pparameters;
    224         if(analysis_counter==(nummodels-1)&& solution_type==TransientSolutionEnum){
    225                 UpdateElementsTransient(elements,parameters,iomodel,analysis_counter,analysis_type);
     65        if(analysis_counter==(nummodels-1)&& solution_enum==TransientSolutionEnum){
     66                UpdateElementsTransient(elements,parameters,iomodel,analysis_counter,analysis_enum);
    22667        }
    22768        #endif

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r16470 r16539  
    2525        bool        ispdd,isdelta18o;
    2626
    27         /*parameters for mass flux: {{{*/
     27        /*parameters for mass flux:*/
    2828        int          mass_flux_num_profiles     = 0;
    2929        bool         qmu_mass_flux_present      = false;
    … …  
    3737        IssmDouble  *matrix                     = NULL;
    3838        int          count;
    39         /*}}}*/
    4039
    4140        if(*pparameters)return; //do not create parameters twice!
    … …  
    9493
    9594        /*For stress balance only*/
     95        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum));
     96        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum));
     97        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum));
     98        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum));
     99        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum));
     100        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum));
    96101        if(iomodel->meshtype==Mesh3DEnum)
    97102         parameters->AddObject(iomodel->CopyConstantObject(MeshNumberoflayersEnum));
    … …  
    221226        CreateOutputDefinitions(&parameters,iomodel);
    222227
    223         /*Solution specific parameters*/
    224         #ifdef _HAVE_HYDROLOGY_
    225         CreateParametersHydrologyShreve(&parameters,iomodel,solution_type,analysis_type);
    226         CreateParametersHydrologyDCInefficient(&parameters,iomodel,solution_type,analysis_type);
    227         CreateParametersHydrologyDCEfficient(&parameters,iomodel,solution_type,analysis_type);
    228         #endif
    229        
    230         #ifdef _HAVE_DAMAGE_
    231         CreateParametersDamage(&parameters,iomodel,solution_type,analysis_type);
    232         #endif
    233 
    234228        /*Before returning, create parameters in case we are running Qmu or control types runs: */
    235229        #ifdef _HAVE_CONTROL_
    … …  
    241235        #endif
    242236
    243         #ifdef _HAVE_STRESSBALANCE_
    244         CreateParametersStressbalance(&parameters,iomodel,solution_type,analysis_type);
    245         #endif
    246 
    247 
    248237        /*Now, deal with toolkits options, which need to be put into the parameters dataset: */
    249238        ParseToolkitsOptionsx(parameters,toolkitsoptionsfid);

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp

    r16334 r16539  
    1919        int i;
    2020
    21         int my_rank;
    22         int num_procs;
     21        const int RIFTINFOSIZE = 12;
    2322        int numberofelements2d;
    2423        int numberofvertices2d;
    … …  
    4140
    4241        /*Get my_rank:*/
    43         my_rank   = IssmComm::GetRank();
    44         num_procs = IssmComm::GetSize();
     42        int my_rank   = IssmComm::GetRank();
     43        int num_procs = IssmComm::GetSize();
    4544
    4645        /*Fetch parameters: */
    47 
    4846        iomodel->Constant(&numrifts,RiftsNumriftsEnum);
    4947
    5048        /*First, check that partitioning has not yet been carryed out. Just check whether my_elements pointers is not already assigned a value: */
    51         if (*pmy_elements)return;
     49        if(*pmy_elements)return;
    5250
    5351        /*Number of vertices per elements, needed to correctly retrieve data: */

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h

    r16534 r16539  
    66#define _MODEL_PROCESSORX_H_
    77
    8 #define RIFTINFOSIZE 12
    9 
    108#include "../../classes/classes.h"
     9#include "../../analyses/analyses.h"
    1110
    1211void ModelProcessorx(Elements** pelements, Nodes** pnodes, Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads, Parameters** pparameters, FILE* iomodel_handle,FILE* toolkitfile, char* rootpath,const int solution_type,const int nummodels,const int* analysis_type_listh);
    … …  
    1918void CreateParametersControl(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type);
    2019void CreateParametersDakota(Parameters** pparameters,IoModel* iomodel,char* rootpath,int solution_type,int analysis_type);
    21 void CreateParametersStressbalance(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type);
    22 void CreateParametersHydrologyShreve(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type);
    23 void CreateParametersHydrologyDCInefficient(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type);
    24 void CreateParametersHydrologyDCEfficient(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type);
    2520void CreateOutputDefinitions(Parameters** pparameters,IoModel* iomodel);
    2621void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel);
    2722void UpdateElementsAndMaterialsDakota(Elements* elements,Materials* materials, IoModel* iomodel);
     23void UpdateElementsTransient(Elements* elements,Parameters* parameters,IoModel* iomodel,int analysis_counter,int analysis_type);
    2824void CreateNodes(Nodes** pnodes, IoModel* iomodel,int analysis,int finite_element,int approximation=NoneApproximationEnum);
    29 
    30 /*Creation of fem datasets: specialised drivers: */
    31 
    32 /*stressbalance horizontal*/
    33 void CreateNodesStressbalance(Nodes** pnodes,IoModel* iomodel);
    34 void CreateConstraintsStressbalance(Constraints** pconstraints,IoModel* iomodel);
    35 void CreateLoadsStressbalance(Loads** ploads, IoModel* iomodel);
    36 void UpdateElementsStressbalance(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    37 
    38 /*stressbalance vertical*/
    39 void CreateNodesStressbalanceVertical(Nodes** pnodes,IoModel* iomodel);
    40 void CreateConstraintsStressbalanceVertical(Constraints** pconstraints,IoModel* iomodel);
    41 void CreateLoadsStressbalanceVertical(Loads** ploads, IoModel* iomodel);
    42 void UpdateElementsStressbalanceVertical(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    43 
    44 /*stressbalance SIA*/
    45 void CreateNodesStressbalanceSIA(Nodes** pnodes,IoModel* iomodel);
    46 void CreateConstraintsStressbalanceSIA(Constraints** pconstraints,IoModel* iomodel);
    47 void CreateLoadsStressbalanceSIA(Loads** ploads, IoModel* iomodel);
    48 void UpdateElementsStressbalanceSIA(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    49 
    50 #ifdef _HAVE_GIA_
    51 /*gia*/
    52 void CreateNodesGia(Nodes** pnodes,IoModel* iomodel);
    53 void CreateConstraintsGia(Constraints** pconstraints,IoModel* iomodel);
    54 void CreateLoadsGia(Loads** ploads, IoModel* iomodel);
    55 void UpdateElementsGia(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    56 #endif
    57 
    58 #ifdef _HAVE_DAMAGE_
    59 /*damage*/
    60 void CreateNodesDamage(Nodes** pnodes,IoModel* iomodel);
    61 void CreateConstraintsDamage(Constraints** pconstraints,IoModel* iomodel);
    62 void CreateLoadsDamage(Loads** ploads, IoModel* iomodel);
    63 void UpdateElementsDamage(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    64 void CreateParametersDamage(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type);
    65 #endif
    66 
    67 /*Extrude from base*/
    68 void CreateNodesExtrudeFromBase(Nodes** pnodes,IoModel* iomodel);
    69 void CreateConstraintsExtrudeFromBase(Constraints** pconstraints,IoModel* iomodel);
    70 void CreateLoadsExtrudeFromBase(Loads** ploads, IoModel* iomodel);
    71 void UpdateElementsExtrudeFromBase(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    72 
    73 /*Extrude from top*/
    74 void CreateNodesExtrudeFromTop(Nodes** pnodes,IoModel* iomodel);
    75 void CreateConstraintsExtrudeFromTop(Constraints** pconstraints,IoModel* iomodel);
    76 void CreateLoadsExtrudeFromTop(Loads** ploads, IoModel* iomodel);
    77 void UpdateElementsExtrudeFromTop(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    78 
    79 /*L2 projection base*/
    80 void CreateNodesL2ProjectionBase(Nodes** pnodes,IoModel* iomodel);
    81 void CreateConstraintsL2ProjectionBase(Constraints** pconstraints,IoModel* iomodel);
    82 void CreateLoadsL2ProjectionBase(Loads** ploads, IoModel* iomodel);
    83 void UpdateElementsL2ProjectionBase(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    84 
    85 /*L2 projection top*/
    86 void CreateNodesL2ProjectionTop(Nodes** pnodes,IoModel* iomodel);
    87 void CreateConstraintsL2ProjectionTop(Constraints** pconstraints,IoModel* iomodel);
    88 void CreateLoadsL2ProjectionTop(Loads** ploads, IoModel* iomodel);
    89 void UpdateElementsL2ProjectionTop(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    90 
    91 /*thermal:*/
    92 void CreateNodesThermal(Nodes** pnodes,IoModel* iomodel);
    93 void CreateConstraintsThermal(Constraints** pconstraints,IoModel* iomodel);
    94 void CreateLoadsThermal(Loads** ploads, IoModel* iomodel);
    95 void UpdateElementsThermal(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    96 
    97 /*enthalpy:*/
    98 void CreateNodesEnthalpy(Nodes** pnodes,IoModel* iomodel);
    99 void CreateConstraintsEnthalpy(Constraints** pconstraints,IoModel* iomodel);
    100 void CreateLoadsEnthalpy(Loads** ploads, IoModel* iomodel);
    101 void UpdateElementsEnthalpy(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    102 
    103 /*hydrology Shreve:*/
    104 void CreateNodesHydrologyShreve(Nodes** pnodes,IoModel* iomodel);
    105 void CreateConstraintsHydrologyShreve(Constraints** pconstraints,IoModel* iomodel);
    106 void CreateLoadsHydrologyShreve(Loads** ploads, IoModel* iomodel);
    107 void UpdateElementsHydrologyShreve(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    108 
    109 /*hydrology DC:*/
    110 void CreateNodesHydrologyDCInefficient(Nodes** pnodes,IoModel* iomodel);
    111 void CreateConstraintsHydrologyDCInefficient(Constraints** pconstraints,IoModel* iomodel);
    112 void CreateLoadsHydrologyDCInefficient(Loads** ploads, IoModel* iomodel);
    113 void UpdateElementsHydrologyDCInefficient(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    114 void CreateNodesHydrologyDCEfficient(Nodes** pnodes,IoModel* iomodel);
    115 void CreateConstraintsHydrologyDCEfficient(Constraints** pconstraints,IoModel* iomodel);
    116 void CreateLoadsHydrologyDCEfficient(Loads** ploads, IoModel* iomodel);
    117 void UpdateElementsHydrologyDCEfficient(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    118 
    119 /*melting:*/
    120 void CreateNodesMelting(Nodes** pnodes,IoModel* iomodel);
    121 void CreateConstraintsMelting(Constraints** pconstraints,IoModel* iomodel);
    122 void CreateLoadsMelting(Loads** ploads, IoModel* iomodel);
    123 void UpdateElementsMelting(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    124 
    125 /*masstransport:*/
    126 void CreateNodesMasstransport(Nodes** pnodes,IoModel* iomodel);
    127 void CreateConstraintsMasstransport(Constraints** pconstraints,IoModel* iomodel);
    128 void CreateLoadsMasstransport(Loads** ploads, IoModel* iomodel);
    129 void UpdateElementsMasstransport(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    130 void CreateNodesFreeSurfaceTop(Nodes** pnodes,IoModel* iomodel);
    131 void CreateConstraintsFreeSurfaceTop(Constraints** pconstraints,IoModel* iomodel);
    132 void CreateLoadsFreeSurfaceTop(Loads** ploads, IoModel* iomodel);
    133 void UpdateElementsFreeSurfaceTop(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    134 void CreateNodesFreeSurfaceBase(Nodes** pnodes,IoModel* iomodel);
    135 void CreateConstraintsFreeSurfaceBase(Constraints** pconstraints,IoModel* iomodel);
    136 void CreateLoadsFreeSurfaceBase(Loads** ploads, IoModel* iomodel);
    137 void UpdateElementsFreeSurfaceBase(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    138 
    139 /*balancedthickness:*/
    140 void CreateNodesBalancethickness(Nodes** pnodes,IoModel* iomodel);
    141 void CreateConstraintsBalancethickness(Constraints** pconstraints,IoModel* iomodel);
    142 void CreateLoadsBalancethickness(Loads** ploads, IoModel* iomodel);
    143 void UpdateElementsBalancethickness(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    144 
    145 /*balancedvelocity:*/
    146 void CreateNodesBalancevelocity(Nodes** pnodes,IoModel* iomodel);
    147 void CreateConstraintsBalancevelocity(Constraints** pconstraints,IoModel* iomodel);
    148 void CreateLoadsBalancevelocity(Loads** ploads, IoModel* iomodel);
    149 void UpdateElementsBalancevelocity(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    150 
    151 /*transient: */
    152 void UpdateElementsTransient(Elements* elements,Parameters* parameters,IoModel* iomodel,int analysis_counter,int analysis_type);
    15325
    15426/*partitioning: */
Note: See TracChangeset for help on using the changeset viewer.