Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 16538) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 16539) @@ -376,25 +376,5 @@ #}}} #Masstransport sources {{{ -masstransport_sources = ./modules/ModelProcessorx/Masstransport/UpdateElementsMasstransport.cpp\ - ./modules/ModelProcessorx/Masstransport/CreateNodesMasstransport.cpp\ - ./modules/ModelProcessorx/Masstransport/CreateConstraintsMasstransport.cpp\ - ./modules/ModelProcessorx/Masstransport/CreateLoadsMasstransport.cpp\ - ./modules/ModelProcessorx/FreeSurfaceTop/UpdateElementsFreeSurfaceTop.cpp\ - ./modules/ModelProcessorx/FreeSurfaceTop/CreateNodesFreeSurfaceTop.cpp\ - ./modules/ModelProcessorx/FreeSurfaceTop/CreateConstraintsFreeSurfaceTop.cpp\ - ./modules/ModelProcessorx/FreeSurfaceTop/CreateLoadsFreeSurfaceTop.cpp\ - ./modules/ModelProcessorx/FreeSurfaceBase/UpdateElementsFreeSurfaceBase.cpp\ - ./modules/ModelProcessorx/FreeSurfaceBase/CreateNodesFreeSurfaceBase.cpp\ - ./modules/ModelProcessorx/FreeSurfaceBase/CreateConstraintsFreeSurfaceBase.cpp\ - ./modules/ModelProcessorx/FreeSurfaceBase/CreateLoadsFreeSurfaceBase.cpp\ - ./modules/ModelProcessorx/ExtrudeFromBase/UpdateElementsExtrudeFromBase.cpp\ - ./modules/ModelProcessorx/ExtrudeFromBase/CreateNodesExtrudeFromBase.cpp \ - ./modules/ModelProcessorx/ExtrudeFromBase/CreateConstraintsExtrudeFromBase.cpp\ - ./modules/ModelProcessorx/ExtrudeFromBase/CreateLoadsExtrudeFromBase.cpp\ - ./modules/ModelProcessorx/ExtrudeFromTop/UpdateElementsExtrudeFromTop.cpp\ - ./modules/ModelProcessorx/ExtrudeFromTop/CreateNodesExtrudeFromTop.cpp \ - ./modules/ModelProcessorx/ExtrudeFromTop/CreateConstraintsExtrudeFromTop.cpp\ - ./modules/ModelProcessorx/ExtrudeFromTop/CreateLoadsExtrudeFromTop.cpp\ - ./analyses/ExtrudeFromBaseAnalysis.h\ +masstransport_sources = ./analyses/ExtrudeFromBaseAnalysis.h\ ./analyses/ExtrudeFromBaseAnalysis.cpp\ ./analyses/ExtrudeFromTopAnalysis.h\ @@ -412,17 +392,5 @@ #}}} #Thermal sources {{{ -thermal_sources = ./modules/ModelProcessorx/Thermal/UpdateElementsThermal.cpp\ - ./modules/ModelProcessorx/Thermal/CreateNodesThermal.cpp\ - ./modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp\ - ./modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp\ - ./modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp\ - ./modules/ModelProcessorx/Enthalpy/CreateNodesEnthalpy.cpp\ - ./modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp\ - ./modules/ModelProcessorx/Enthalpy/CreateLoadsEnthalpy.cpp\ - ./modules/ModelProcessorx/Melting/UpdateElementsMelting.cpp\ - ./modules/ModelProcessorx/Melting/CreateNodesMelting.cpp\ - ./modules/ModelProcessorx/Melting/CreateConstraintsMelting.cpp\ - ./modules/ModelProcessorx/Melting/CreateLoadsMelting.cpp\ - ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.h\ +thermal_sources = ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.h\ ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.cpp\ ./analyses/ThermalAnalysis.h\ @@ -493,20 +461,5 @@ #}}} #Hydrology sources {{{ -hydrology_sources = ./modules/ModelProcessorx/HydrologyShreve/UpdateElementsHydrologyShreve.cpp\ - ./modules/ModelProcessorx/HydrologyShreve/CreateNodesHydrologyShreve.cpp\ - ./modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp\ - ./modules/ModelProcessorx/HydrologyShreve/CreateLoadsHydrologyShreve.cpp \ - ./modules/ModelProcessorx/HydrologyShreve/CreateParametersHydrologyShreve.cpp \ - ./modules/ModelProcessorx/HydrologyDCInefficient/UpdateElementsHydrologyDCInefficient.cpp\ - ./modules/ModelProcessorx/HydrologyDCInefficient/CreateNodesHydrologyDCInefficient.cpp\ - ./modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp\ - ./modules/ModelProcessorx/HydrologyDCInefficient/CreateLoadsHydrologyDCInefficient.cpp \ - ./modules/ModelProcessorx/HydrologyDCInefficient/CreateParametersHydrologyDCInefficient.cpp \ - ./modules/ModelProcessorx/HydrologyDCEfficient/UpdateElementsHydrologyDCEfficient.cpp\ - ./modules/ModelProcessorx/HydrologyDCEfficient/CreateNodesHydrologyDCEfficient.cpp\ - ./modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp\ - ./modules/ModelProcessorx/HydrologyDCEfficient/CreateLoadsHydrologyDCEfficient.cpp \ - ./modules/ModelProcessorx/HydrologyDCEfficient/CreateParametersHydrologyDCEfficient.cpp \ - ./analyses/HydrologyDCEfficientAnalysis.h\ +hydrology_sources = ./analyses/HydrologyDCEfficientAnalysis.h\ ./analyses/HydrologyDCEfficientAnalysis.cpp\ ./analyses/HydrologyDCInefficientAnalysis.h\ @@ -518,18 +471,5 @@ #}}} #Stressbalance sources {{{ -stressbalance_sources = ./modules/ModelProcessorx/Stressbalance/UpdateElementsStressbalance.cpp\ - ./modules/ModelProcessorx/Stressbalance/CreateNodesStressbalance.cpp \ - ./modules/ModelProcessorx/Stressbalance/CreateConstraintsStressbalance.cpp \ - ./modules/ModelProcessorx/Stressbalance/CreateLoadsStressbalance.cpp\ - ./modules/ModelProcessorx/Stressbalance/CreateParametersStressbalance.cpp\ - ./modules/ModelProcessorx/StressbalanceVertical/UpdateElementsStressbalanceVertical.cpp\ - ./modules/ModelProcessorx/StressbalanceVertical/CreateNodesStressbalanceVertical.cpp \ - ./modules/ModelProcessorx/StressbalanceVertical/CreateConstraintsStressbalanceVertical.cpp \ - ./modules/ModelProcessorx/StressbalanceVertical/CreateLoadsStressbalanceVertical.cpp\ - ./modules/ModelProcessorx/StressbalanceSIA/UpdateElementsStressbalanceSIA.cpp\ - ./modules/ModelProcessorx/StressbalanceSIA/CreateNodesStressbalanceSIA.cpp \ - ./modules/ModelProcessorx/StressbalanceSIA/CreateConstraintsStressbalanceSIA.cpp \ - ./modules/ModelProcessorx/StressbalanceSIA/CreateLoadsStressbalanceSIA.cpp \ - ./analyses/StressbalanceAnalysis.h\ +stressbalance_sources = ./analyses/StressbalanceAnalysis.h\ ./analyses/StressbalanceAnalysis.cpp\ ./analyses/StressbalanceSIAAnalysis.h\ @@ -541,13 +481,5 @@ #}}} #Balanced sources {{{ -balanced_sources = ./modules/ModelProcessorx/Balancethickness/UpdateElementsBalancethickness.cpp\ - ./modules/ModelProcessorx/Balancethickness/CreateNodesBalancethickness.cpp\ - ./modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp\ - ./modules/ModelProcessorx/Balancethickness/CreateLoadsBalancethickness.cpp\ - ./modules/ModelProcessorx/Balancevelocity/UpdateElementsBalancevelocity.cpp\ - ./modules/ModelProcessorx/Balancevelocity/CreateNodesBalancevelocity.cpp\ - ./modules/ModelProcessorx/Balancevelocity/CreateConstraintsBalancevelocity.cpp\ - ./modules/ModelProcessorx/Balancevelocity/CreateLoadsBalancevelocity.cpp\ - ./analyses/BalancevelocityAnalysis.h\ +balanced_sources = ./analyses/BalancevelocityAnalysis.h\ ./analyses/BalancevelocityAnalysis.cpp\ ./analyses/SmoothedSurfaceSlopeXAnalysis.h\ @@ -564,13 +496,5 @@ #}}} #Slope sources {{{ -slope_sources = ./modules/ModelProcessorx/L2ProjectionBase/UpdateElementsL2ProjectionBase.cpp\ - ./modules/ModelProcessorx/L2ProjectionBase/CreateNodesL2ProjectionBase.cpp \ - ./modules/ModelProcessorx/L2ProjectionBase/CreateConstraintsL2ProjectionBase.cpp\ - ./modules/ModelProcessorx/L2ProjectionBase/CreateLoadsL2ProjectionBase.cpp\ - ./modules/ModelProcessorx/L2ProjectionTop/UpdateElementsL2ProjectionTop.cpp\ - ./modules/ModelProcessorx/L2ProjectionTop/CreateNodesL2ProjectionTop.cpp \ - ./modules/ModelProcessorx/L2ProjectionTop/CreateConstraintsL2ProjectionTop.cpp\ - ./modules/ModelProcessorx/L2ProjectionTop/CreateLoadsL2ProjectionTop.cpp\ - ./analyses/L2ProjectionBaseAnalysis.h\ +slope_sources = ./analyses/L2ProjectionBaseAnalysis.h\ ./analyses/L2ProjectionBaseAnalysis.cpp\ ./cores/surfaceslope_core.cpp\ @@ -578,9 +502,5 @@ #}}} #MeshDeformation sources {{{ -meshdeformation_sources = ./modules/ModelProcessorx/MeshDeformation/UpdateElementsMeshDeformation.cpp\ - ./modules/ModelProcessorx/MeshDeformation/CreateNodesMeshDeformation.cpp \ - ./modules/ModelProcessorx/MeshDeformation/CreateConstraintsMeshDeformation.cpp\ - ./modules/ModelProcessorx/MeshDeformation/CreateLoadsMeshDeformation.cpp\ - ./analyses/MeshdeformationAnalysis.h\ +meshdeformation_sources = ./analyses/MeshdeformationAnalysis.h\ ./analyses/MeshdeformationAnalysis.cpp\ ./cores/meshdeformation_core.cpp @@ -590,8 +510,4 @@ ./analyses/GiaAnalysis.h\ ./analyses/GiaAnalysis.cpp\ - ./modules/ModelProcessorx/Gia/UpdateElementsGia.cpp\ - ./modules/ModelProcessorx/Gia/CreateNodesGia.cpp \ - ./modules/ModelProcessorx/Gia/CreateConstraintsGia.cpp\ - ./modules/ModelProcessorx/Gia/CreateLoadsGia.cpp\ ./modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp\ ./modules/GiaDeflectionCorex/GiaDeflectionCorex.h\ @@ -610,9 +526,4 @@ ./analyses/DamageEvolutionAnalysis.h\ ./analyses/DamageEvolutionAnalysis.cpp\ - ./modules/ModelProcessorx/Damage/UpdateElementsDamage.cpp\ - ./modules/ModelProcessorx/Damage/CreateNodesDamage.cpp \ - ./modules/ModelProcessorx/Damage/CreateConstraintsDamage.cpp\ - ./modules/ModelProcessorx/Damage/CreateParametersDamage.cpp\ - ./modules/ModelProcessorx/Damage/CreateLoadsDamage.cpp\ ./solutionsequences/solutionsequence_damage_nonlinear.cpp Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (revision 16539) @@ -6,5 +6,20 @@ /*Model processor*/ -int AdjointBalancethicknessAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int AdjointBalancethicknessAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void AdjointBalancethicknessAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void AdjointBalancethicknessAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void AdjointBalancethicknessAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void AdjointBalancethicknessAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void AdjointBalancethicknessAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 16539) @@ -6,5 +6,20 @@ /*Model processing*/ -int AdjointHorizAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int AdjointHorizAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ _error_("not implemented"); }/*}}}*/ +void AdjointHorizAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void AdjointHorizAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void AdjointHorizAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void AdjointHorizAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void AdjointHorizAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/Analysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/Analysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/Analysis.h (revision 16539) @@ -6,10 +6,22 @@ #define _ANALYSIS_H_ +class Parameters; +class IoModel; +class Elements; +class Nodes; +class Constraints; +class Loads; + class Analysis{ public: - virtual ~Analysis(){}; - virtual int DofsPerNode(int** doflist,int meshtype,int approximation)=0; + virtual ~Analysis(){}; + virtual int DofsPerNode(int** doflist,int meshtype,int approximation)=0; + virtual void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum)=0; + virtual void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type)=0; + virtual void CreateNodes(Nodes** pnodes,IoModel* iomodel)=0; + virtual void CreateConstraints(Constraints** pconstraints,IoModel* iomodel)=0; + virtual void CreateLoads(Loads** ploads, IoModel* iomodel)=0; }; #endif Index: /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp (revision 16539) @@ -6,5 +6,119 @@ /*Model processing*/ -int BalancethicknessAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int BalancethicknessAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void BalancethicknessAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void BalancethicknessAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + int stabilization,finiteelement; + + /*Fetch data needed: */ + iomodel->Constant(&stabilization,BalancethicknessStabilizationEnum); + + /*Finite element type*/ + finiteelement = P1Enum; + if(stabilization==3){ + finiteelement = P1DGEnum; + } + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement); + counter++; + } + } + + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,VxEnum); + iomodel->FetchDataToInput(elements,VyEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum); + iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum); + + if(iomodel->meshtype==Mesh3DEnum){ + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + } +}/*}}}*/ +void BalancethicknessAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + int stabilization; + iomodel->Constant(&stabilization,BalancethicknessStabilizationEnum); + + /*Check in 3d*/ + if(stabilization==3 && iomodel->meshtype==Mesh3DEnum) _error_("DG 3d not implemented yet"); + + /*First fetch data: */ + if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); + if(stabilization!=3){ + ::CreateNodes(pnodes,iomodel,BalancethicknessAnalysisEnum,P1Enum); + } + else{ + ::CreateNodes(pnodes,iomodel,BalancethicknessAnalysisEnum,P1DGEnum); + } + iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); +}/*}}}*/ +void BalancethicknessAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + + /*Fetch parameters: */ + int stabilization; + iomodel->Constant(&stabilization,BalancethicknessStabilizationEnum); + + /*Recover pointer: */ + Constraints* constraints=*pconstraints; + + /*Do not add constraints in DG*/ + if(stabilization!=3){ + IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,BalancethicknessAnalysisEnum,P1Enum); + } + + /*Assign output pointer: */ + *pconstraints=constraints; +}/*}}}*/ +void BalancethicknessAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*Intermediary*/ + int element; + int stabilization; + + /*Fetch parameters: */ + iomodel->Constant(&stabilization,BalancethicknessStabilizationEnum); + + /*Recover pointer: */ + Loads* loads=*ploads; + + /*Loads only in DG*/ + if (stabilization==3){ + + /*Get faces and elements*/ + CreateFaces(iomodel); + iomodel->FetchData(1,ThicknessEnum); + + /*First load data:*/ + for(int i=0;inumberoffaces;i++){ + + /*Get left and right elements*/ + element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2] + + /*Now, if this element is not in the partition, pass: */ + if(!iomodel->my_elements[element]) continue; + + /* Add load */ + loads->AddObject(new Numericalflux(iomodel->loadcounter+i+1,i,i,iomodel,BalancethicknessAnalysisEnum)); + } + + /*Free data: */ + iomodel->DeleteData(1,ThicknessEnum); + } + + /*Assign output pointer: */ + *ploads=loads; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/BalancethicknessSoftAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/BalancethicknessSoftAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/BalancethicknessSoftAnalysis.cpp (revision 16539) @@ -6,5 +6,20 @@ /*Model processing*/ -int BalancethicknessSoftAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int BalancethicknessSoftAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ _error_("not implemented"); }/*}}}*/ +void BalancethicknessSoftAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void BalancethicknessSoftAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void BalancethicknessSoftAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void BalancethicknessSoftAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void BalancethicknessSoftAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/BalancethicknessSoftAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/BalancethicknessSoftAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/BalancethicknessSoftAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp (revision 16539) @@ -6,5 +6,59 @@ /*Model processing*/ -int BalancevelocityAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int BalancevelocityAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void BalancevelocityAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void BalancevelocityAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,VxEnum); + iomodel->FetchDataToInput(elements,VyEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum); + iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum); + + if(iomodel->meshtype==Mesh3DEnum){ + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + } +}/*}}}*/ +void BalancevelocityAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + /*Check in 3d*/ + if(iomodel->meshtype==Mesh3DEnum) _error_("DG 3d not implemented yet"); + + /*First fetch data: */ + iomodel->FetchData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationVertexEquationEnum); + ::CreateNodes(pnodes,iomodel,BalancevelocityAnalysisEnum,P1Enum); + iomodel->DeleteData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationVertexEquationEnum); +}/*}}}*/ +void BalancevelocityAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + + /*Recover pointer: */ + Constraints* constraints=*pconstraints; + + /*No constraints for now*/ + //IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,BalancevelocityAnalysisEnum,P1Enum); + + /*Assign output pointer: */ + *pconstraints=constraints; +}/*}}}*/ +void BalancevelocityAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*No loads*/ +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 16539) @@ -6,5 +6,76 @@ /*Model processing*/ -int DamageEvolutionAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int DamageEvolutionAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void DamageEvolutionAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void DamageEvolutionAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + /*What input do I need to run my damage evolution model?*/ + iomodel->FetchDataToInput(elements,VxEnum); + iomodel->FetchDataToInput(elements,VyEnum); + iomodel->FetchDataToInput(elements,VzEnum); + iomodel->FetchDataToInput(elements,DamageDEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,PressureEnum); + + bool dakota_analysis; + iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum); + if(dakota_analysis){ + elements->InputDuplicate(DamageDEnum, QmuDamageDEnum); + } +}/*}}}*/ +void DamageEvolutionAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + iomodel->FetchData(1,MeshVertexonbedEnum); + ::CreateNodes(pnodes,iomodel,DamageEvolutionAnalysisEnum,P1Enum); + iomodel->DeleteData(1,MeshVertexonbedEnum); +}/*}}}*/ +void DamageEvolutionAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + + /*Fetch parameters: */ + int stabilization; + iomodel->Constant(&stabilization,DamageStabilizationEnum); + + /*Recover pointer: */ + Constraints* constraints=*pconstraints; + + IoModelToConstraintsx(constraints,iomodel,DamageSpcdamageEnum,DamageEvolutionAnalysisEnum,P1Enum); + + /*Assign output pointer: */ + *pconstraints=constraints; +}/*}}}*/ +void DamageEvolutionAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*Recover pointer: */ + Loads* loads=*ploads; + + /*create penalties for nodes: no node can have a damage > 1*/ + iomodel->FetchData(1,DamageSpcdamageEnum); + CreateSingleNodeToElementConnectivity(iomodel); + + for(int i=0;inumberofvertices;i++){ + + /*keep only this partition's nodes:*/ + if(iomodel->my_vertices[i]){ + if (xIsNan(iomodel->Data(DamageSpcdamageEnum)[i])){ //No penalty applied on spc nodes! + loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,DamageEvolutionAnalysisEnum)); + } + } + } + iomodel->DeleteData(1,DamageSpcdamageEnum); + + /*Assign output pointer: */ + *ploads=loads; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 16539) @@ -6,5 +6,172 @@ /*Model processing*/ -int EnthalpyAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int EnthalpyAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void EnthalpyAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void EnthalpyAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + bool dakota_analysis; + bool isenthalpy; + + /*Now, is the model 3d? otherwise, do nothing: */ + if(iomodel->meshtype==Mesh2DhorizontalEnum)return; + + /*Is enthalpy requested?*/ + iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum); + if(!isenthalpy) return; + + /*Fetch data needed: */ + iomodel->FetchData(3,TemperatureEnum,WaterfractionEnum,PressureEnum); + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum); + + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,FrictionCoefficientEnum); + iomodel->FetchDataToInput(elements,FrictionPEnum); + iomodel->FetchDataToInput(elements,FrictionQEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum); + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + iomodel->FetchDataToInput(elements,FlowequationElementEquationEnum); + iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum); + iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum); + iomodel->FetchDataToInput(elements,PressureEnum); + iomodel->FetchDataToInput(elements,TemperatureEnum); + iomodel->FetchDataToInput(elements,WaterfractionEnum); + iomodel->FetchDataToInput(elements,EnthalpyEnum); + iomodel->FetchDataToInput(elements,BasalforcingsGeothermalfluxEnum); + iomodel->FetchDataToInput(elements,WatercolumnEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,VxEnum); + iomodel->FetchDataToInput(elements,VyEnum); + iomodel->FetchDataToInput(elements,VzEnum); + InputUpdateFromConstantx(elements,0.,VxMeshEnum); + InputUpdateFromConstantx(elements,0.,VyMeshEnum); + InputUpdateFromConstantx(elements,0.,VzMeshEnum); + if(dakota_analysis){ + elements->InputDuplicate(TemperatureEnum,QmuTemperatureEnum); + elements->InputDuplicate(BasalforcingsMeltingRateEnum,QmuMeltingEnum); + elements->InputDuplicate(VxMeshEnum,QmuVxMeshEnum); + elements->InputDuplicate(VxMeshEnum,QmuVyMeshEnum); + elements->InputDuplicate(VxMeshEnum,QmuVzMeshEnum); + } + + /*Free data: */ + iomodel->DeleteData(3,TemperatureEnum,WaterfractionEnum,PressureEnum); +}/*}}}*/ +void EnthalpyAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); + ::CreateNodes(pnodes,iomodel,EnthalpyAnalysisEnum,P1Enum); + iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); +}/*}}}*/ +void EnthalpyAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + + /*Intermediary*/ + int count; + int M,N; + bool spcpresent = false; + IssmDouble heatcapacity; + IssmDouble referencetemperature; + + /*Output*/ + IssmDouble *spcvector = NULL; + IssmDouble* times=NULL; + IssmDouble* values=NULL; + + /*Fetch parameters: */ + iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum); + iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum); + + /*Recover pointer: */ + Constraints* constraints=*pconstraints; + + /*return if 2d mesh*/ + if(iomodel->meshtype==Mesh2DhorizontalEnum) return; + + /*Fetch data: */ + iomodel->FetchData(&spcvector,&M,&N,ThermalSpctemperatureEnum); + + //FIX ME: SHOULD USE IOMODELCREATECONSTRAINTS + /*Transient or static?:*/ + if(M==iomodel->numberofvertices){ + /*static: just create Constraints objects*/ + count=0; + + for(int i=0;inumberofvertices;i++){ + /*keep only this partition's nodes:*/ + if((iomodel->my_vertices[i])){ + + if (!xIsNan(spcvector[i])){ + + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,heatcapacity*(spcvector[i]-referencetemperature),EnthalpyAnalysisEnum)); + count++; + + } + } + } + } + else if (M==(iomodel->numberofvertices+1)){ + /*transient: create transient SpcTransient objects. Same logic, except we need to retrieve + * various times and values to initialize an SpcTransient object: */ + count=0; + + /*figure out times: */ + times=xNew(N); + for(int j=0;jnumberofvertices;i++){ + + /*keep only this partition's nodes:*/ + if((iomodel->my_vertices[i])){ + + /*figure out times and values: */ + values=xNew(N); + spcpresent=false; + for(int j=0;j(values[j]))spcpresent=true; //NaN means no spc by default + } + + if(spcpresent){ + constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,N,times,values,EnthalpyAnalysisEnum)); + count++; + } + xDelete(values); + } + } + } + else{ + _error_("Size of field " << EnumToStringx(ThermalSpctemperatureEnum) << " not supported"); + } + + /*Free ressources:*/ + iomodel->DeleteData(spcvector,ThermalSpctemperatureEnum); + xDelete(times); + xDelete(values); + + /*Assign output pointer: */ + *pconstraints=constraints; +}/*}}}*/ +void EnthalpyAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*No loads */ +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.h (revision 16539) @@ -7,2 +7,8 @@ #endif + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp (revision 16539) @@ -6,5 +6,31 @@ /*Model processing*/ -int ExtrudeFromBaseAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int ExtrudeFromBaseAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void ExtrudeFromBaseAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void ExtrudeFromBaseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + if(iomodel->meshtype==Mesh2DverticalEnum){ + iomodel->FetchDataToInput(elements,MeshVertexonbedEnum); + } +}/*}}}*/ +void ExtrudeFromBaseAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + ::CreateNodes(pnodes,iomodel,ExtrudeFromBaseAnalysisEnum,P1Enum); + +}/*}}}*/ +void ExtrudeFromBaseAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void ExtrudeFromBaseAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp (revision 16539) @@ -6,5 +6,31 @@ /*Model processing*/ -int ExtrudeFromTopAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int ExtrudeFromTopAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void ExtrudeFromTopAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void ExtrudeFromTopAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + if(iomodel->meshtype==Mesh2DverticalEnum){ + iomodel->FetchDataToInput(elements,MeshVertexonbedEnum); + } +}/*}}}*/ +void ExtrudeFromTopAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + ::CreateNodes(pnodes,iomodel,ExtrudeFromBaseAnalysisEnum,P1Enum); + +}/*}}}*/ +void ExtrudeFromTopAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void ExtrudeFromTopAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 16539) @@ -6,5 +6,90 @@ /*Model processing*/ -int FreeSurfaceBaseAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int FreeSurfaceBaseAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void FreeSurfaceBaseAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void FreeSurfaceBaseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Now, is the model 3d? otherwise, do nothing: */ + if (iomodel->meshtype==Mesh2DhorizontalEnum)return; + + /*Finite element type*/ + int finiteelement = P1Enum; + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement); + counter++; + } + } + + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateCorrectionEnum,0.); + iomodel->FetchDataToInput(elements,VxEnum); + iomodel->FetchDataToInput(elements,VyEnum); + if(iomodel->meshtype==Mesh3DEnum){ + iomodel->FetchDataToInput(elements,VzEnum); + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + } +}/*}}}*/ +void FreeSurfaceBaseAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); + ::CreateNodes(pnodes,iomodel,FreeSurfaceBaseAnalysisEnum,P1Enum); + iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); +}/*}}}*/ +void FreeSurfaceBaseAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void FreeSurfaceBaseAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*Intermediaries*/ + int penpair_ids[2]; + int count=0; + int numvertex_pairing; + + /*Recover pointer: */ + Loads* loads=*ploads; + + /*Create Penpair for vertex_pairing: */ + IssmDouble *vertex_pairing=NULL; + IssmDouble *nodeonbed=NULL; + iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); + iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum); + for(int i=0;imy_vertices[reCast(vertex_pairing[2*i+0])-1]){ + + /*In debugging mode, check that the second node is in the same cpu*/ + _assert_(iomodel->my_vertices[reCast(vertex_pairing[2*i+1])-1]); + + /*Skip if one of the two is not on the bed*/ + if(!(reCast(nodeonbed[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonbed[reCast(vertex_pairing[2*i+1])-1]))) continue; + + /*Get node ids*/ + penpair_ids[0]=iomodel->nodecounter+reCast(vertex_pairing[2*i+0]); + penpair_ids[1]=iomodel->nodecounter+reCast(vertex_pairing[2*i+1]); + + /*Create Load*/ + loads->AddObject(new Penpair( + iomodel->loadcounter+count+1, + &penpair_ids[0], + FreeSurfaceBaseAnalysisEnum)); + count++; + } + } + + /*free ressources: */ + iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum); + iomodel->DeleteData(nodeonbed,MeshVertexonbedEnum); + + /*Assign output pointer: */ + *ploads=loads; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp (revision 16539) @@ -6,5 +6,88 @@ /*Model processing*/ -int FreeSurfaceTopAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int FreeSurfaceTopAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void FreeSurfaceTopAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void FreeSurfaceTopAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Now, is the model 3d? otherwise, do nothing: */ + if (iomodel->meshtype==Mesh2DhorizontalEnum)return; + + int finiteelement = P1Enum; + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement); + counter++; + } + } + + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum,0.); + iomodel->FetchDataToInput(elements,VxEnum); + iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum); + if(iomodel->meshtype==Mesh3DEnum){ + iomodel->FetchDataToInput(elements,VzEnum); + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + } +}/*}}}*/ +void FreeSurfaceTopAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); + ::CreateNodes(pnodes,iomodel,FreeSurfaceTopAnalysisEnum,P1Enum); + iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); +}/*}}}*/ +void FreeSurfaceTopAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void FreeSurfaceTopAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*Intermediaries*/ + int penpair_ids[2]; + int count=0; + int numvertex_pairing; + + /*Recover pointer: */ + Loads* loads=*ploads; + + /*Create Penpair for vertex_pairing: */ + IssmDouble *vertex_pairing=NULL; + IssmDouble *nodeonsurface=NULL; + iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); + iomodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum); + for(int i=0;imy_vertices[reCast(vertex_pairing[2*i+0])-1]){ + + /*In debugging mode, check that the second node is in the same cpu*/ + _assert_(iomodel->my_vertices[reCast(vertex_pairing[2*i+1])-1]); + + /*Skip if one of the two is not on the bed*/ + if(!(reCast(nodeonsurface[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonsurface[reCast(vertex_pairing[2*i+1])-1]))) continue; + + /*Get node ids*/ + penpair_ids[0]=iomodel->nodecounter+reCast(vertex_pairing[2*i+0]); + penpair_ids[1]=iomodel->nodecounter+reCast(vertex_pairing[2*i+1]); + + /*Create Load*/ + loads->AddObject(new Penpair( + iomodel->loadcounter+count+1, + &penpair_ids[0], + FreeSurfaceTopAnalysisEnum)); + count++; + } + } + + /*free ressources: */ + iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum); + iomodel->DeleteData(nodeonsurface,MeshVertexonsurfaceEnum); + + /*Assign output pointer: */ + *ploads=loads; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/GiaAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/GiaAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/GiaAnalysis.cpp (revision 16539) @@ -6,5 +6,32 @@ /*Model processing*/ -int GiaAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int GiaAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void GiaAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void GiaAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,GiaMantleViscosityEnum); + iomodel->FetchDataToInput(elements,GiaLithosphereThicknessEnum); +}/*}}}*/ +void GiaAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + ::CreateNodes(pnodes,iomodel,GiaAnalysisEnum,P1Enum); +}/*}}}*/ +void GiaAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + /*No constraints*/ +}/*}}}*/ +void GiaAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + /*No loads*/ +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/GiaAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/GiaAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/GiaAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 16539) @@ -6,5 +6,123 @@ /*Model processing*/ -int HydrologyDCEfficientAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int HydrologyDCEfficientAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void HydrologyDCEfficientAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + Parameters *parameters = NULL; + int hydrology_model; + bool isefficientlayer; + + /*Get parameters: */ + parameters=*pparameters; + + /*retrieve some parameters: */ + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + + /*Now, do we really want DC?*/ + if(hydrology_model!=HydrologydcEnum){ + *pparameters=parameters; + return; + } + + /*Do we want an efficient layer*/ + iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if(!isefficientlayer){ + *pparameters=parameters; + return; + } + + /*Nothing for now*/ + + /*Assign output pointer: */ + *pparameters=parameters; +}/*}}}*/ +void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + bool isefficientlayer; + int hydrology_model; + + /*Now, do we really want DC?*/ + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + if(hydrology_model!=HydrologydcEnum) return; + + /*Do we want an efficient layer*/ + iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if(!isefficientlayer) return; + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,EplHeadEnum); +}/*}}}*/ +void HydrologyDCEfficientAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + /*Now, do we really want DC?*/ + int hydrology_model; + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + if(hydrology_model!=HydrologydcEnum) return; + + /*Do we want an efficient layer*/ + bool isefficientlayer; + iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if(!isefficientlayer) return; + + iomodel->FetchData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationVertexEquationEnum); + ::CreateNodes(pnodes,iomodel,HydrologyDCEfficientAnalysisEnum,P1Enum); + iomodel->DeleteData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationVertexEquationEnum); +}/*}}}*/ +void HydrologyDCEfficientAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + + /*Recover pointer: */ + Constraints* constraints=*pconstraints; + + /*Do we really want DC?*/ + int hydrology_model; + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + if(hydrology_model!=HydrologydcEnum) return; + + /*Do we want an efficient layer*/ + bool isefficientlayer; + iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if(!isefficientlayer) return; + + IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum,P1Enum); + + /*Assign output pointer: */ + *pconstraints=constraints; +}/*}}}*/ +void HydrologyDCEfficientAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*Recover pointer: */ + Loads* loads=*ploads; + + /*Do we really want DC?*/ + int hydrology_model; + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + if(hydrology_model!=HydrologydcEnum) return; + + /*Do we want an efficient layer*/ + bool isefficientlayer; + iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if(!isefficientlayer) return; + + /*Nothing for now*/ + + /*Assign output pointer: */ + *ploads=loads; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 16539) @@ -6,5 +6,148 @@ /*Model processing*/ -int HydrologyDCInefficientAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int HydrologyDCInefficientAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void HydrologyDCInefficientAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + Parameters *parameters = NULL; + int hydrology_model; + int sedimentlimit_flag; + int transfer_flag; + bool isefficientlayer; + IssmDouble sedimentlimit; + IssmDouble penalty_factor; + IssmDouble leakagefactor; + IssmDouble rel_tol; + + /*Get parameters: */ + parameters=*pparameters; + + /*retrieve some parameters: */ + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + + /*Now, do we really want DC?*/ + if(hydrology_model!=HydrologydcEnum){ + *pparameters=parameters; + return; + } + + iomodel->FetchData(&isefficientlayer,HydrologydcIsefficientlayerEnum); + iomodel->FetchData(&sedimentlimit_flag,HydrologydcSedimentlimitFlagEnum); + iomodel->FetchData(&transfer_flag,HydrologydcTransferFlagEnum); + iomodel->FetchData(&penalty_factor,HydrologydcPenaltyFactorEnum); + iomodel->FetchData(&rel_tol,HydrologydcRelTolEnum); + + if(sedimentlimit_flag==1){ + iomodel->FetchData(&sedimentlimit,HydrologydcSedimentlimitEnum); + parameters->AddObject(new DoubleParam(HydrologydcSedimentlimitEnum,sedimentlimit)); + } + + if(transfer_flag==1){ + iomodel->FetchData(&leakagefactor,HydrologydcLeakageFactorEnum); + parameters->AddObject(new DoubleParam(HydrologydcLeakageFactorEnum,leakagefactor)); + } + + parameters->AddObject(new DoubleParam(HydrologydcPenaltyFactorEnum,penalty_factor)); + parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model)); + parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer)); + parameters->AddObject(new IntParam(HydrologydcSedimentlimitFlagEnum,sedimentlimit_flag)); + parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag)); + parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol)); + + /*Assign output pointer: */ + *pparameters=parameters; +}/*}}}*/ +void HydrologyDCInefficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + bool isefficientlayer; + int hydrology_model; + + /*Fetch data needed: */ + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + + /*Now, do we really want DC?*/ + if(hydrology_model!=HydrologydcEnum) return; + + /*Fetch data needed: */ + iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum); + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,SedimentHeadEnum); + if(isefficientlayer)iomodel->FetchDataToInput(elements,HydrologydcMaskEplactiveEnum); +}/*}}}*/ +void HydrologyDCInefficientAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + /*Fetch parameters: */ + int hydrology_model; + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + + /*Now, do we really want DC?*/ + if(hydrology_model!=HydrologydcEnum) return; + + if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); + ::CreateNodes(pnodes,iomodel,HydrologyDCInefficientAnalysisEnum,P1Enum); + iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); +}/*}}}*/ +void HydrologyDCInefficientAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + + /*Recover pointer: */ + Constraints* constraints=*pconstraints; + + /*retrieve some parameters: */ + int hydrology_model; + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + if(hydrology_model!=HydrologydcEnum) return; + + IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyDCInefficientAnalysisEnum,P1Enum); + + /*Assign output pointer: */ + *pconstraints=constraints; +}/*}}}*/ +void HydrologyDCInefficientAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*Recover pointer: */ + Loads* loads=*ploads; + + /*Fetch parameters: */ + int hydrology_model; + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + if(hydrology_model!=HydrologydcEnum) return; + + iomodel->FetchData(1,MeshVertexonbedEnum); + + //create penalties for nodes: no node can have a temperature over the melting point + CreateSingleNodeToElementConnectivity(iomodel); + for(int i=0;inumberofvertices;i++){ + if (iomodel->meshtype==Mesh3DEnum){ + /*keep only this partition's nodes:*/ + if(iomodel->my_vertices[i]){ + loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum)); + } + } + else if(reCast(iomodel->Data(MeshVertexonbedEnum)[i])){ + if(iomodel->my_vertices[i]){ + loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum)); + } + } + } + /*Assign output pointer: */ + *ploads=loads; + iomodel->DeleteData(1,MeshVertexonbedEnum); +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 16539) @@ -6,5 +6,91 @@ /*Model processing*/ -int HydrologyShreveAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int HydrologyShreveAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void HydrologyShreveAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + Parameters *parameters = NULL; + int hydrology_model; + + /*Get parameters: */ + parameters=*pparameters; + + /*retrieve some parameters: */ + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + + /*Now, do we really want Shreve?*/ + if(hydrology_model!=HydrologyshreveEnum){ + *pparameters=parameters; + return; + } + + parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model)); + parameters->AddObject(iomodel->CopyConstantObject(HydrologyshreveStabilizationEnum)); + + /*Assign output pointer: */ + *pparameters=parameters; +}/*}}}*/ +void HydrologyShreveAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Fetch data needed: */ + int hydrology_model; + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + + /*Now, do we really want Shreve?*/ + if(hydrology_model!=HydrologyshreveEnum) return; + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,WatercolumnEnum); + + elements->InputDuplicate(WatercolumnEnum,WaterColumnOldEnum); +}/*}}}*/ +void HydrologyShreveAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + /*Fetch parameters: */ + int hydrology_model; + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + + /*Now, do we really want Shreve?*/ + if(hydrology_model!=HydrologyshreveEnum) return; + + if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); + ::CreateNodes(pnodes,iomodel,HydrologyShreveAnalysisEnum,P1Enum); + iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); +}/*}}}*/ +void HydrologyShreveAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + + /*Recover pointer: */ + Constraints* constraints=*pconstraints; + + /*retrieve some parameters: */ + int hydrology_model; + iomodel->Constant(&hydrology_model,HydrologyModelEnum); + + if(hydrology_model!=HydrologyshreveEnum) return; + + IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum,P1Enum); + + /*Assign output pointer: */ + *pconstraints=constraints; +}/*}}}*/ +void HydrologyShreveAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + /*No loads*/ +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 16539) @@ -6,5 +6,49 @@ /*Model processing*/ -int L2ProjectionBaseAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int L2ProjectionBaseAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void L2ProjectionBaseAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void L2ProjectionBaseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + if(iomodel->meshtype==Mesh3DEnum){ + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + } + if(iomodel->meshtype==Mesh2DverticalEnum){ + iomodel->FetchDataToInput(elements,MeshVertexonbedEnum); + } +}/*}}}*/ +void L2ProjectionBaseAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + if(iomodel->meshtype==Mesh3DEnum){ + iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); + } + else if(iomodel->meshtype==Mesh2DverticalEnum){ + iomodel->FetchData(1,MeshVertexonbedEnum); + } + ::CreateNodes(pnodes,iomodel,L2ProjectionBaseAnalysisEnum,P1Enum); + iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); +}/*}}}*/ +void L2ProjectionBaseAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + + /*No constraints*/ +}/*}}}*/ +void L2ProjectionBaseAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*No loads*/ +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 16539) @@ -6,5 +6,208 @@ /*Model processing*/ -int MasstransportAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int MasstransportAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void MasstransportAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + /*Get parameters: */ + Parameters *parameters=*pparameters; + + parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum)); + + /*Assign output pointer: */ + *pparameters = parameters; +}/*}}}*/ +void MasstransportAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + int stabilization,finiteelement; + bool dakota_analysis; + bool issmbgradients; + bool ispdd; + bool isdelta18o; + bool isgroundingline; + + /*Fetch data needed: */ + iomodel->Constant(&stabilization,MasstransportStabilizationEnum); + iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum); + iomodel->Constant(&ispdd,SurfaceforcingsIspddEnum); + iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum); + iomodel->Constant(&issmbgradients,SurfaceforcingsIssmbgradientsEnum); + iomodel->Constant(&isgroundingline,TransientIsgroundinglineEnum); + + /*Finite element type*/ + finiteelement = P1Enum; + if(stabilization==3){ + finiteelement = P1DGEnum; + } + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement); + counter++; + } + } + + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateCorrectionEnum,0.); + iomodel->FetchDataToInput(elements,VxEnum); + iomodel->FetchDataToInput(elements,VyEnum); + + if(stabilization==3){ + iomodel->FetchDataToInput(elements,MasstransportSpcthicknessEnum); //for DG, we need the spc in the element + } + + if(dakota_analysis){ + elements->InputDuplicate(BedEnum,QmuBedEnum); + elements->InputDuplicate(ThicknessEnum,QmuThicknessEnum); + elements->InputDuplicate(SurfaceEnum,QmuSurfaceEnum); + elements->InputDuplicate(MaskIceLevelsetEnum,QmuMaskIceLevelsetEnum); + if(isgroundingline) elements->InputDuplicate(MaskGroundediceLevelsetEnum,QmuMaskGroundediceLevelsetEnum); + } + + if(iomodel->meshtype==Mesh3DEnum){ + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + } + + if(issmbgradients){ + iomodel->FetchDataToInput(elements,SurfaceforcingsHrefEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsSmbrefEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsBPosEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsBNegEnum); + } + if(ispdd){ + iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum); + if(isdelta18o){ + iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesLgmEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesPresentdayEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsPresentdayEnum); + } + else{ + iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsMonthlytemperaturesEnum); + } + } + if(~ispdd && ~issmbgradients){ + iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum,0.); + } +}/*}}}*/ +void MasstransportAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + /*Fetch parameters: */ + int stabilization; + iomodel->Constant(&stabilization,MasstransportStabilizationEnum); + + /*Check in 3d*/ + if(stabilization==3 && iomodel->meshtype==Mesh3DEnum) _error_("DG 3d not implemented yet"); + + /*Create Nodes either DG or CG depending on stabilization*/ + if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); + if(stabilization!=3){ + ::CreateNodes(pnodes,iomodel,MasstransportAnalysisEnum,P1Enum); + } + else{ + ::CreateNodes(pnodes,iomodel,MasstransportAnalysisEnum,P1DGEnum); + } + iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); +}/*}}}*/ +void MasstransportAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + + /*Fetch parameters: */ + int stabilization; + iomodel->Constant(&stabilization,MasstransportStabilizationEnum); + + /*Recover pointer: */ + Constraints* constraints=*pconstraints; + + /*Do not add constraints in DG, they are weakly imposed*/ + if(stabilization!=3){ + IoModelToConstraintsx(constraints,iomodel,MasstransportSpcthicknessEnum,MasstransportAnalysisEnum,P1Enum); + } + + /*Assign output pointer: */ + *pconstraints=constraints; +}/*}}}*/ +void MasstransportAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*Intermediaries*/ + int element; + int penpair_ids[2]; + int count=0; + int stabilization; + int numvertex_pairing; + + /*Fetch parameters: */ + iomodel->Constant(&stabilization,MasstransportStabilizationEnum); + + /*Recover pointer: */ + Loads* loads=*ploads; + + /*Loads only in DG*/ + if (stabilization==3){ + + /*Get faces and elements*/ + CreateFaces(iomodel); + iomodel->FetchData(1,ThicknessEnum); + + /*First load data:*/ + for(int i=0;inumberoffaces;i++){ + + /*Get left and right elements*/ + element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2] + + /*Now, if this element is not in the partition, pass: */ + if(!iomodel->my_elements[element]) continue; + + /* Add load */ + loads->AddObject(new Numericalflux(iomodel->loadcounter+i+1,i,i,iomodel,MasstransportAnalysisEnum)); + } + + /*Free data: */ + iomodel->DeleteData(1,ThicknessEnum); + } + + /*Create Penpair for vertex_pairing: */ + IssmDouble *vertex_pairing=NULL; + IssmDouble *nodeonbed=NULL; + iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); + iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum); + + for(int i=0;imy_vertices[reCast(vertex_pairing[2*i+0])-1]){ + + /*In debugging mode, check that the second node is in the same cpu*/ + _assert_(iomodel->my_vertices[reCast(vertex_pairing[2*i+1])-1]); + + /*Skip if one of the two is not on the bed*/ + if(!(reCast(nodeonbed[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonbed[reCast(vertex_pairing[2*i+1])-1]))) continue; + + /*Get node ids*/ + penpair_ids[0]=iomodel->nodecounter+reCast(vertex_pairing[2*i+0]); + penpair_ids[1]=iomodel->nodecounter+reCast(vertex_pairing[2*i+1]); + + /*Create Load*/ + loads->AddObject(new Penpair( + iomodel->loadcounter+count+1, + &penpair_ids[0], + MasstransportAnalysisEnum)); + count++; + } + } + + /*free ressources: */ + iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum); + iomodel->DeleteData(nodeonbed,MeshVertexonbedEnum); + + /*Assign output pointer: */ + *ploads=loads; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp (revision 16539) @@ -6,5 +6,71 @@ /*Model processing*/ -int MeltingAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int MeltingAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void MeltingAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void MeltingAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Now, is the model 3d? otherwise, do nothing: */ + if(iomodel->meshtype==Mesh2DhorizontalEnum)return; + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + /*Create inputs: */ + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,FrictionCoefficientEnum); + iomodel->FetchDataToInput(elements,FrictionPEnum); + iomodel->FetchDataToInput(elements,FrictionQEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + iomodel->FetchDataToInput(elements,FlowequationElementEquationEnum); + iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum); + iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,PressureEnum); +}/*}}}*/ +void MeltingAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); + ::CreateNodes(pnodes,iomodel,MeltingAnalysisEnum,P1Enum); + iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); +}/*}}}*/ +void MeltingAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + /*No Constraints*/ +}/*}}}*/ +void MeltingAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*if 2d: Error*/ + if(iomodel->meshtype==Mesh2DhorizontalEnum) _error_("2d meshes not supported yet"); + + /*Recover pointer: */ + Loads* loads=*ploads; + + //create penalties for nodes: no node can have a temperature over the melting point + iomodel->FetchData(1,MeshVertexonbedEnum); + CreateSingleNodeToElementConnectivity(iomodel); + + for(int i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + if (reCast(iomodel->Data(MeshVertexonbedEnum)[i])){ + loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,MeltingAnalysisEnum)); + } + } + } + iomodel->DeleteData(1,MeshVertexonbedEnum); + + /*Assign output pointer: */ + *ploads=loads; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/MeshdeformationAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MeshdeformationAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/MeshdeformationAnalysis.cpp (revision 16539) @@ -6,5 +6,20 @@ /*Model processing*/ -int MeshdeformationAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int MeshdeformationAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ _error_("not implemented"); }/*}}}*/ +void MeshdeformationAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void MeshdeformationAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void MeshdeformationAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void MeshdeformationAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void MeshdeformationAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/MeshdeformationAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/MeshdeformationAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/MeshdeformationAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.cpp (revision 16539) @@ -6,5 +6,39 @@ /*Model processing*/ -int SmoothedSurfaceSlopeXAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int SmoothedSurfaceSlopeXAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void SmoothedSurfaceSlopeXAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void SmoothedSurfaceSlopeXAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + if(iomodel->meshtype==Mesh3DEnum){ + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + } + if(iomodel->meshtype==Mesh2DverticalEnum){ + iomodel->FetchDataToInput(elements,MeshVertexonbedEnum); + } +}/*}}}*/ +void SmoothedSurfaceSlopeXAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + iomodel->FetchData(1,MeshVertexonbedEnum); + ::CreateNodes(pnodes,iomodel,SmoothedSurfaceSlopeXAnalysisEnum,P1Enum); + iomodel->DeleteData(1,MeshVertexonbedEnum); +}/*}}}*/ +void SmoothedSurfaceSlopeXAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void SmoothedSurfaceSlopeXAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.cpp (revision 16539) @@ -6,5 +6,39 @@ /*Model processing*/ -int SmoothedSurfaceSlopeYAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int SmoothedSurfaceSlopeYAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void SmoothedSurfaceSlopeYAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void SmoothedSurfaceSlopeYAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + if(iomodel->meshtype==Mesh3DEnum){ + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + } + if(iomodel->meshtype==Mesh2DverticalEnum){ + iomodel->FetchDataToInput(elements,MeshVertexonbedEnum); + } +}/*}}}*/ +void SmoothedSurfaceSlopeYAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + iomodel->FetchData(1,MeshVertexonbedEnum); + ::CreateNodes(pnodes,iomodel,SmoothedSurfaceSlopeYAnalysisEnum,P1Enum); + iomodel->DeleteData(1,MeshVertexonbedEnum); +}/*}}}*/ +void SmoothedSurfaceSlopeYAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ +}/*}}}*/ +void SmoothedSurfaceSlopeYAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16539) @@ -6,5 +6,5 @@ /*Model processing*/ -int StressbalanceAnalysis::DofsPerNode(int** pdoftype,int meshtype,int approximation){/*{{{*/ +int StressbalanceAnalysis::DofsPerNode(int** pdoftype,int meshtype,int approximation){/*{{{*/ /*output*/ @@ -66,2 +66,761 @@ return numdofs; }/*}}}*/ +void StressbalanceAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + /*Intermediaries*/ + int numoutputs; + char** requestedoutputs = NULL; + + /*Get parameters: */ + Parameters *parameters=*pparameters; + + parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSIAEnum)); + parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSSAEnum)); + parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsL1L2Enum)); + parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsHOEnum)); + parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum)); + parameters->AddObject(iomodel->CopyConstantObject(FlowequationFeFSEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceIsnewtonEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalancePenaltyFactorEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceFSreconditioningEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceShelfDampeningEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceViscosityOvershootEnum)); + + /*Requested outputs*/ + iomodel->FetchData(&requestedoutputs,&numoutputs,StressbalanceRequestedOutputsEnum); + parameters->AddObject(new IntParam(StressbalanceNumRequestedOutputsEnum,numoutputs)); + if(numoutputs)parameters->AddObject(new StringArrayParam(StressbalanceRequestedOutputsEnum,requestedoutputs,numoutputs)); + iomodel->DeleteData(&requestedoutputs,numoutputs,StressbalanceRequestedOutputsEnum); + + /*Assign output pointer: */ + *pparameters = parameters; +}/*}}}*/ +void StressbalanceAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Intermediaries*/ + int materials_type,finiteelement; + int approximation; + int* finiteelement_list=NULL; + bool isSSA,isL1L2,isHO,isFS,iscoupling; + bool control_analysis; + bool dakota_analysis; + + /*Fetch constants needed: */ + iomodel->Constant(&isSSA,FlowequationIsSSAEnum); + iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum); + iomodel->Constant(&isHO,FlowequationIsHOEnum); + iomodel->Constant(&isFS,FlowequationIsFSEnum); + iomodel->Constant(&control_analysis,InversionIscontrolEnum); + iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum); + iomodel->Constant(&materials_type,MaterialsEnum); + + /*return if no processing required*/ + if(!isSSA & !isL1L2 & !isHO & !isFS) return; + + /*Fetch data needed and allocate vectors: */ + iomodel->FetchData(1,FlowequationElementEquationEnum); + finiteelement_list=xNewZeroInit(iomodel->numberofelements); + + /*Do we have coupling*/ + if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) + iscoupling = true; + else + iscoupling = false; + + /*Get finite element type*/ + if(!iscoupling){ + if(isSSA) iomodel->Constant(&finiteelement,FlowequationFeSSAEnum); + else if(isL1L2) finiteelement = P1Enum; + else if(isHO) iomodel->Constant(&finiteelement,FlowequationFeHOEnum); + else if(isFS) iomodel->Constant(&finiteelement,FlowequationFeFSEnum); + for(int i=0;inumberofelements;i++){ + finiteelement_list[i]=finiteelement; + } + } + else{ + if(isFS){ + for(int i=0;inumberofelements;i++){ + approximation=reCast(iomodel->Data(FlowequationElementEquationEnum)[i]); + if(approximation==FSApproximationEnum || approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){ + finiteelement_list[i]=MINIcondensedEnum; + } + else{ + finiteelement_list[i]=P1Enum; + } + } + } + else{ + finiteelement = P1Enum; + for(int i=0;inumberofelements;i++){ + finiteelement_list[i]=finiteelement; + } + } + } + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement_list[i]); + counter++; + } + } + + /*Create inputs: */ + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,FrictionCoefficientEnum); + iomodel->FetchDataToInput(elements,FrictionPEnum); + iomodel->FetchDataToInput(elements,FrictionQEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum); + iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum); + iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum); + iomodel->FetchDataToInput(elements,VxEnum,0.); + if(dakota_analysis)elements->InputDuplicate(VxEnum,QmuVxEnum); + iomodel->FetchDataToInput(elements,VyEnum,0.); + if(dakota_analysis)elements->InputDuplicate(VyEnum,QmuVyEnum); + iomodel->FetchDataToInput(elements,LoadingforceXEnum); + iomodel->FetchDataToInput(elements,LoadingforceYEnum); + iomodel->FetchDataToInput(elements,DamageDEnum); + + if(iomodel->meshtype==Mesh3DEnum){ + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,FlowequationBorderFSEnum); + iomodel->FetchDataToInput(elements,LoadingforceZEnum); + iomodel->FetchDataToInput(elements,VzEnum,0.); + if(dakota_analysis)elements->InputDuplicate(VzEnum,QmuVzEnum); + } + if(isFS){ + iomodel->FetchDataToInput(elements,MeshVertexonbedEnum); + iomodel->FetchDataToInput(elements,PressureEnum,0.); + if(dakota_analysis)elements->InputDuplicate(PressureEnum,QmuPressureEnum); + } + +#ifdef _HAVE_ANDROID_ + elements->InputDuplicate(FrictionCoefficientEnum,AndroidFrictionCoefficientEnum); +#endif + + /*Free data: */ + iomodel->DeleteData(1,FlowequationElementEquationEnum); + xDelete(finiteelement_list); +}/*}}}*/ +void StressbalanceAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + /*Intermediary*/ + bool isSSA,isL1L2,isHO,isFS,iscoupling; + int finiteelement=-1,approximation=-1; + + /*Fetch parameters: */ + iomodel->Constant(&isSSA,FlowequationIsSSAEnum); + iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum); + iomodel->Constant(&isHO,FlowequationIsHOEnum); + iomodel->Constant(&isFS,FlowequationIsFSEnum); + + /*Now, check that we have non SIA elements */ + if(!isSSA & !isL1L2 & !isHO & !isFS) return; + + /*Do we have coupling*/ + if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) + iscoupling = true; + else + iscoupling = false; + + /*If no coupling, call Regular CreateNodes, else, use P1 elements only*/ + if(!iscoupling){ + + /*Get finite element type*/ + if(isSSA){ + approximation=SSAApproximationEnum; + iomodel->Constant(&finiteelement,FlowequationFeSSAEnum); + } + else if(isL1L2){ + approximation = L1L2ApproximationEnum; + finiteelement = P1Enum; + } + else if(isHO){ + approximation = HOApproximationEnum; + iomodel->Constant(&finiteelement,FlowequationFeHOEnum); + } + else if(isFS){ + approximation = FSApproximationEnum; + iomodel->Constant(&finiteelement,FlowequationFeFSEnum); + } + iomodel->FetchData(3,FlowequationBorderSSAEnum,FlowequationVertexEquationEnum,StressbalanceReferentialEnum); + if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderFSEnum); + ::CreateNodes(pnodes,iomodel,StressbalanceAnalysisEnum,finiteelement,approximation); + iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, + FlowequationVertexEquationEnum,StressbalanceReferentialEnum); + } + else{ + /*Coupling: we are going to create P1 Elements only*/ + + /*Recover nodes*/ + Nodes* nodes = *pnodes; + Node* node = NULL; + int lid=0; + if(!nodes) nodes = new Nodes(); + + iomodel->FetchData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, + FlowequationVertexEquationEnum,StressbalanceReferentialEnum); + if(isFS){ + /*P1+ velocity*/ + for(int i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + approximation=reCast(iomodel->Data(FlowequationVertexEquationEnum)[i]); + if(approximation==FSApproximationEnum) approximation=FSvelocityEnum; + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,approximation)); + } + } + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,StressbalanceAnalysisEnum,FSvelocityEnum); + node->Deactivate(); + nodes->AddObject(node); + } + } + /*P1 pressure*/ + for(int i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + approximation=reCast(iomodel->Data(FlowequationVertexEquationEnum)[i]); + node = new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,i,iomodel,StressbalanceAnalysisEnum,FSpressureEnum); + if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum){ + node->Deactivate(); + } + nodes->AddObject(node); + } + } + } + else{ + for(int i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,reCast(iomodel->Data(FlowequationVertexEquationEnum)[i]))); + } + } + } + iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, + FlowequationVertexEquationEnum,StressbalanceReferentialEnum); + + /*Assign output pointer: */ + *pnodes=nodes; + } +}/*}}}*/ +void StressbalanceAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + + /*Intermediary*/ + int i,j; + int count,finiteelement; + IssmDouble g; + IssmDouble rho_ice; + IssmDouble FSreconditioning; + bool isSIA,isSSA,isL1L2,isHO,isFS,iscoupling; + bool spcpresent = false; + int Mx,Nx; + int My,Ny; + int Mz,Nz; + IssmDouble *spcvx = NULL; + IssmDouble *spcvy = NULL; + IssmDouble *spcvz = NULL; + IssmDouble *nodeonSSA = NULL; + IssmDouble *nodeonHO = NULL; + IssmDouble *nodeonFS = NULL; + IssmDouble *nodeonbed = NULL; + IssmDouble *groundedice_ls = NULL; + IssmDouble *vertices_type = NULL; + IssmDouble *surface = NULL; + IssmDouble *z = NULL; + IssmDouble *timesx=NULL; + IssmDouble *timesy=NULL; + IssmDouble *timesz=NULL; + IssmDouble* values=NULL; + + /*Output*/ + Constraints *constraints = NULL; + + /*Fetch parameters: */ + iomodel->Constant(&g,ConstantsGEnum); + iomodel->Constant(&rho_ice,MaterialsRhoIceEnum); + iomodel->Constant(&FSreconditioning,StressbalanceFSreconditioningEnum); + iomodel->Constant(&isSIA,FlowequationIsSIAEnum); + iomodel->Constant(&isSSA,FlowequationIsSSAEnum); + iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum); + iomodel->Constant(&isHO,FlowequationIsHOEnum); + iomodel->Constant(&isFS,FlowequationIsFSEnum); + + /*Recover pointer: */ + constraints=*pconstraints; + + /*Now, is the flag macayaealHO on? otherwise, do nothing: */ + if(!isSSA && !isHO && !isFS && !isL1L2){ + *pconstraints=constraints; + return; + } + + /*Do we have coupling*/ + if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) + iscoupling = true; + else + iscoupling = false; + + /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/ + if(!iscoupling){ + + /*Get finite element type*/ + if(isSSA) iomodel->Constant(&finiteelement,FlowequationFeSSAEnum); + else if(isL1L2) finiteelement = P1Enum; + else if(isHO) iomodel->Constant(&finiteelement,FlowequationFeHOEnum); + else if(isFS){ iomodel->Constant(&finiteelement,FlowequationFeFSEnum); + /*Deduce velocity interpolation from finite element*/ + switch(finiteelement){ + case P1P1Enum : finiteelement = P1Enum; break; + case P1P1GLSEnum : finiteelement = P1Enum; break; + case MINIcondensedEnum : finiteelement = P1bubbleEnum; break; + case MINIEnum : finiteelement = P1bubbleEnum; break; + case TaylorHoodEnum : finiteelement = P2Enum; break; + default: _error_("finite element "<FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum); + iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum); + iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum); + iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum); + if(iomodel->meshtype==Mesh3DEnum){ + iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum); + } + else if (iomodel->meshtype==Mesh2DverticalEnum){ + iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvyEnum); + } + else{ + _error_("not supported yet"); + } + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + if(nodeonbed[i]>0. && groundedice_ls[i]>0. && nodeonFS[i]>0.){ + if(vertices_type[i] == FSApproximationEnum){ + for(j=0;jmeshtype==Mesh3DEnum){ + IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1); + IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,2); + IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,3); + iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum); + } + else if (iomodel->meshtype==Mesh2DverticalEnum){ + IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1); + IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,2); + iomodel->DeleteData(spcvz,StressbalanceSpcvyEnum); + } + else{ + _error_("not supported yet"); + } + iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum); + iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum); + iomodel->DeleteData(nodeonbed,MeshVertexonbedEnum); + iomodel->DeleteData(groundedice_ls,MaskGroundediceLevelsetEnum); + + /*Pressure spc*/ + count = constraints->Size(); + iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum); + iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum); + iomodel->FetchData(&z,NULL,NULL,MeshZEnum); + switch(finiteelement){ + case P1bubbleEnum: + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + if(reCast(vertices_type[i])==NoneApproximationEnum){ + constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,1,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); + count++; + } + } + } + break; + case P2Enum: + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + if(reCast(vertices_type[i])==NoneApproximationEnum){ + constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+i+1,1,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); + count++; + } + } + } + break; + default: + _error_("not implemented yet"); + } + iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum); + iomodel->DeleteData(surface,SurfaceEnum); + iomodel->DeleteData(z,MeshZEnum); + } + else{ + IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1); + IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,2); + } + + *pconstraints=constraints; + return; + } + + /*Constraints: fetch data: */ + iomodel->FetchData(&spcvx,&Mx,&Nx,StressbalanceSpcvxEnum); + iomodel->FetchData(&spcvy,&My,&Ny,StressbalanceSpcvyEnum); + iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum); + iomodel->FetchData(&nodeonSSA,NULL,NULL,FlowequationBorderSSAEnum); + if(iomodel->meshtype==Mesh3DEnum)iomodel->FetchData(&nodeonHO,NULL,NULL,FlowequationBorderHOEnum); + if(iomodel->meshtype==Mesh3DEnum)iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum); + if(iomodel->meshtype==Mesh3DEnum)iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum); + if(iomodel->meshtype==Mesh3DEnum)iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum); + iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum); + iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum); + iomodel->FetchData(&z,NULL,NULL,MeshZEnum); + + /*Initialize counter: */ + count=0; + + /*figure out times: */ + timesx=xNew(Nx); + for(j=0;j(Ny); + for(j=0;j(Nz); + for(j=0;jnumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + + /*Start with adding spcs of coupling: zero at the border SSA/HO for the appropriate dofs*/ + if(reCast(vertices_type[i]==SSAHOApproximationEnum)){ + /*If grionSSA, spc HO dofs: 3 & 4*/ + if (reCast(nodeonHO[i])){ + 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. + count++; + 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. + count++; + if (!xIsNan(spcvx[i])){ + 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. + count++; + } + if (!xIsNan(spcvy[i])){ + 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. + count++; + } + + } + else if (reCast(nodeonSSA[i])){ + 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. + count++; + 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. + count++; + if (!xIsNan(spcvx[i])){ + 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. + count++; + } + if (!xIsNan(spcvy[i])){ + 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. + count++; + } + + } + else _error_("if vertices_type is SSAHO, you shoud have nodeonHO or nodeonSSA"); + } + /*Also add spcs of coupling: zero at the border HO/FS for the appropriate dofs*/ + else if (reCast(vertices_type[i])==HOFSApproximationEnum){ + /*If grion,HO spc FS dofs: 3 4 & 5*/ + if (reCast(nodeonHO[i])){ + 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. + count++; + 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. + count++; + 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. + count++; + if (!xIsNan(spcvx[i])){ + 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. + count++; + } + if (!xIsNan(spcvy[i])){ + 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. + count++; + } + + } + else if (reCast(nodeonFS[i])){ //spc HO nodes: 1 & 2 + 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. + count++; + 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. + count++; + if (!xIsNan(spcvx[i])){ + 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. + count++; + } + if (!xIsNan(spcvy[i])){ + 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. + count++; + } + if (!xIsNan(spcvz[i])){ + 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. + count++; + } + } + else _error_("if vertices_type is HOFS, you shoud have nodeonHO or nodeonFS"); + } + /*Also add spcs of coupling: zero at the border HO/FS for the appropriate dofs*/ + else if (reCast(vertices_type[i])==SSAFSApproximationEnum){ + /*If grion,HO spc FS dofs: 3 4 & 5*/ + if (reCast(nodeonSSA[i])){ + 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. + count++; + 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. + count++; + 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. + count++; + if (!xIsNan(spcvx[i])){ + 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. + count++; + } + if (!xIsNan(spcvy[i])){ + 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. + count++; + } + + } + else if (reCast(nodeonFS[i])){ //spc SSA nodes: 1 & 2 + 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. + count++; + 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. + count++; + if (!xIsNan(spcvx[i])){ + 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. + count++; + } + if (!xIsNan(spcvy[i])){ + 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. + count++; + } + if (!xIsNan(spcvz[i])){ + 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. + count++; + } + } + else _error_("if vertices_type is SSAFS, you shoud have nodeonSSA or nodeonFS"); + } + /*Now add the regular spcs*/ + else{ + if (Mx==iomodel->numberofvertices && !xIsNan(spcvx[i])){ + 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. + count++; + + } + else if (Mx==iomodel->numberofvertices+1) { + /*figure out times and values: */ + values=xNew(Nx); + spcpresent=false; + for(j=0;j(values[j]))spcpresent=true; //NaN means no spc by default + } + + if(spcpresent){ + constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,Nx,timesx,values,StressbalanceAnalysisEnum)); + count++; + } + xDelete(values); + } + else if (vertices_type[i]==SIAApproximationEnum){ + constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,StressbalanceAnalysisEnum)); + count++; + } + + if (My==iomodel->numberofvertices && !xIsNan(spcvy[i])){ + 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. + count++; + } + else if (My==iomodel->numberofvertices+1){ + /*figure out times and values: */ + values=xNew(Ny); + spcpresent=false; + for(j=0;j(values[j]))spcpresent=true; //NaN means no spc by default + } + if(spcpresent){ + constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,Ny,timesy,values,StressbalanceAnalysisEnum)); + count++; + } + xDelete(values); + } + else if (vertices_type[i]==SIAApproximationEnum){ + constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,StressbalanceAnalysisEnum)); + count++; + } + + if (reCast(vertices_type[i])==FSApproximationEnum || (reCast(vertices_type[i])==NoneApproximationEnum)){ + if (Mz==iomodel->numberofvertices && !xIsNan(spcvz[i])){ + 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 + count++; + } + else if (Mz==iomodel->numberofvertices+1){ + /*figure out times and values: */ + values=xNew(Nz); + spcpresent=false; + for(j=0;j(values[j]))spcpresent=true; //NaN means no spc by default + } + if(spcpresent){ + constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,Nz,timesz,values,StressbalanceAnalysisEnum)); + count++; + } + xDelete(values); + } + + } + if (reCast(vertices_type[i])==NoneApproximationEnum){ + 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 + count++; + } + } + + /*Constraint at the bedrock interface (v.n = vz = 0) (Coordinates will be updated according to the bed slope)*/ + if (iomodel->meshtype==Mesh3DEnum) if(nodeonbed[i]>0. && groundedice_ls[i]>=0. && nodeonFS[i]>0.){ + switch(reCast(vertices_type[i])){ + case SSAFSApproximationEnum: + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,0.,StressbalanceAnalysisEnum)); + count++; + break; + case HOFSApproximationEnum: + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,0.,StressbalanceAnalysisEnum)); + count++; + break; + case FSApproximationEnum: + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0.,StressbalanceAnalysisEnum)); + count++; + break; + default: _error_("Vertex approximation " << EnumToStringx(reCast(vertices_type[i])) << " not supported"); + } + } + } + } + + /*Free data: */ + iomodel->DeleteData(spcvx,StressbalanceSpcvxEnum); + iomodel->DeleteData(spcvy,StressbalanceSpcvyEnum); + iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum); + iomodel->DeleteData(nodeonSSA,FlowequationBorderSSAEnum); + if(iomodel->meshtype==Mesh3DEnum)iomodel->DeleteData(nodeonHO,FlowequationBorderHOEnum); + if(iomodel->meshtype==Mesh3DEnum)iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum); + if(iomodel->meshtype==Mesh3DEnum)iomodel->DeleteData(nodeonbed,MeshVertexonbedEnum); + if(iomodel->meshtype==Mesh3DEnum)iomodel->DeleteData(groundedice_ls,MaskGroundediceLevelsetEnum); + iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum); + iomodel->DeleteData(surface,SurfaceEnum); + iomodel->DeleteData(z,MeshZEnum); + + /*Free resources:*/ + xDelete(timesx); + xDelete(timesy); + xDelete(timesz); + xDelete(values); + + /*Assign output pointer: */ + *pconstraints=constraints; +}/*}}}*/ +void StressbalanceAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*Intermediary*/ + const int RIFTINFOSIZE = 12; + int i; + int count; + int penpair_ids[2]; + bool isSSA,isL1L2,isHO,isFS; + int numpenalties,numrifts,numriftsegments; + IssmDouble *riftinfo = NULL; + IssmDouble *penalties = NULL; + int assert_int; + + /*Fetch parameters: */ + iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum); + iomodel->Constant(&isFS,FlowequationIsFSEnum); + iomodel->Constant(&isSSA,FlowequationIsSSAEnum); + iomodel->Constant(&isHO,FlowequationIsHOEnum); + iomodel->Constant(&numrifts,RiftsNumriftsEnum); + + /*Recover pointer: */ + Loads* loads=*ploads; + + /*Now, is the flag macayaealHO on? otherwise, do nothing: */ + if(!isSSA && !isHO && !isFS && !isL1L2) return; + + /*Initialize counter: */ + count=0; + + /*Create Penpair for penalties: */ + iomodel->FetchData(&penalties,&numpenalties,NULL,StressbalanceVertexPairingEnum); + + for(i=0;imy_vertices[reCast(penalties[2*i+0]-1)]){ + + /*In debugging mode, check that the second node is in the same cpu*/ + assert_int=iomodel->my_vertices[reCast(penalties[2*i+1]-1)]; _assert_(assert_int); + + /*Get node ids*/ + penpair_ids[0]=iomodel->nodecounter+reCast(penalties[2*i+0]); + penpair_ids[1]=iomodel->nodecounter+reCast(penalties[2*i+1]); + + /*Create Load*/ + loads->AddObject(new Penpair(iomodel->loadcounter+count+1,&penpair_ids[0],StressbalanceAnalysisEnum)); + count++; + } + } + + /*free ressources: */ + iomodel->DeleteData(penalties,StressbalanceVertexPairingEnum); + + /*Create Riffront loads for rifts: */ +#ifdef _HAVE_RIFTS_ + if(numrifts){ + iomodel->FetchData(&riftinfo,&numriftsegments,NULL,RiftsRiftstructEnum); + iomodel->FetchData(5,RiftsRiftstructEnum,ThicknessEnum,BedEnum,SurfaceEnum,MaskGroundediceLevelsetEnum); + for(i=0;imy_elements[reCast(*(riftinfo+RIFTINFOSIZE*i+2))-1]){ + loads->AddObject(new Riftfront(iomodel->loadcounter+count+1,i,iomodel,StressbalanceAnalysisEnum)); + count++; + } + } + iomodel->DeleteData(5,RiftsRiftstructEnum,ThicknessEnum,BedEnum,SurfaceEnum,MaskGroundediceLevelsetEnum); + xDelete(riftinfo); + } +#endif + + /*Assign output pointer: */ + *ploads=loads; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (revision 16539) @@ -6,5 +6,142 @@ /*Model processing*/ -int StressbalanceSIAAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int StressbalanceSIAAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 2; }/*}}}*/ +void StressbalanceSIAAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + /*No specific parameters*/ + +}/*}}}*/ +void StressbalanceSIAAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Fetch data needed: */ + bool isSIA; + iomodel->Constant(&isSIA,FlowequationIsSIAEnum); + + /*Now, is the flag SIA on? otherwise, do nothing: */ + if (!isSIA)return; + + iomodel->FetchData(1,FlowequationElementEquationEnum); + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + iomodel->FetchDataToInput(elements,ThicknessEnum); + + /*Free data: */ + iomodel->DeleteData(1,FlowequationElementEquationEnum); +}/*}}}*/ +void StressbalanceSIAAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + /*Intermediaries*/ + bool isSIA; + Node* node = NULL; + + /*Fetch parameters: */ + iomodel->Constant(&isSIA,FlowequationIsSIAEnum); + + /*Now, is the flag isSIA on? otherwise, do nothing: */ + if(!isSIA) return; + + /*First create nodes*/ + Nodes* nodes=*pnodes; + int lid=0; + if(!nodes) nodes = new Nodes(); + + iomodel->FetchData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, + FlowequationVertexEquationEnum,StressbalanceReferentialEnum); + + for(int i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + + /*Create new node if is in this processor's partition*/ + node = new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceSIAAnalysisEnum,reCast(iomodel->Data(FlowequationVertexEquationEnum)[i])); + + /*Deactivate node if not SIA*/ + if(reCast(iomodel->Data(FlowequationVertexEquationEnum)[i])!=SIAApproximationEnum){ + node->Deactivate(); + } + + /*Add to Nodes dataset*/ + nodes->AddObject(node); + } + } + + iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, + FlowequationVertexEquationEnum,StressbalanceReferentialEnum); + + /*Assign output pointer: */ + *pnodes=nodes; + +}/*}}}*/ +void StressbalanceSIAAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + + /*Intermediary*/ + int count; + IssmDouble yts; + bool isSIA; + + /*Output*/ + Constraints* constraints = NULL; + + /*Recover pointer: */ + constraints=*pconstraints; + + /*Fetch parameters: */ + iomodel->Constant(&yts,ConstantsYtsEnum); + iomodel->Constant(&isSIA,FlowequationIsSIAEnum); + + /*Now, is the flag isSIA on? otherwise, do nothing: */ + if (!isSIA) return; + + /*Fetch data: */ + iomodel->FetchData(3,StressbalanceSpcvxEnum,StressbalanceSpcvyEnum,FlowequationVertexEquationEnum); + + /*Initialize conunter*/ + count=0; + + /*vx and vy are spc'd if we are not on nodeonSIA: */ + for(int i=0;inumberofvertices;i++){ + /*keep only this partition's nodes:*/ + if((iomodel->my_vertices[i])){ + if (reCast(iomodel->Data(FlowequationVertexEquationEnum)[i])!=SIAApproximationEnum){ + + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceSIAAnalysisEnum)); + count++; + + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceSIAAnalysisEnum)); + count++; + } + else{ + if (!xIsNan(iomodel->Data(StressbalanceSpcvxEnum)[i])){ + 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. + count++; + } + + if (!xIsNan(iomodel->Data(StressbalanceSpcvyEnum)[i])){ + 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 + count++; + } + } + } + } + + /*Free data: */ + iomodel->DeleteData(3,StressbalanceSpcvxEnum,StressbalanceSpcvyEnum,FlowequationVertexEquationEnum); + + /*Assign output pointer: */ + *pconstraints=constraints; +}/*}}}*/ +void StressbalanceSIAAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*No loads*/ + +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 16539) @@ -6,5 +6,95 @@ /*Model processing*/ -int StressbalanceVerticalAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int StressbalanceVerticalAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void StressbalanceVerticalAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + /*No specific parameters*/ + +}/*}}}*/ +void StressbalanceVerticalAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*return if not 3d mesh*/ + if(iomodel->meshtype!=Mesh3DEnum) return; + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,VxEnum,0.); + iomodel->FetchDataToInput(elements,VyEnum,0.); +}/*}}}*/ +void StressbalanceVerticalAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + /*return if not 3d mesh*/ + if(iomodel->meshtype!=Mesh3DEnum) return; + + iomodel->FetchData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationVertexEquationEnum); + ::CreateNodes(pnodes,iomodel,StressbalanceVerticalAnalysisEnum,P1Enum); + iomodel->DeleteData(3,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationVertexEquationEnum); +}/*}}}*/ +void StressbalanceVerticalAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + + /*Intermediary*/ + int count; + IssmDouble yts; + + /*Fetch parameters: */ + iomodel->Constant(&yts,ConstantsYtsEnum); + + /*Recover pointer: */ + Constraints* constraints=*pconstraints; + + /*return if not 3d mesh*/ + if(iomodel->meshtype!=Mesh3DEnum) return; + + /*Fetch data: */ + iomodel->FetchData(2,StressbalanceSpcvzEnum,FlowequationBorderFSEnum); + + /*Initialize counter*/ + count=0; + + /*Create spcs from x,y,z, as well as the spc values on those spcs: */ + for(int i=0;inumberofvertices;i++){ + + /*keep only this partition's nodes:*/ + if(iomodel->my_vertices[i]){ + + if (reCast(iomodel->Data(FlowequationBorderFSEnum)[i])){ + 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 + count++; + } + else if (!xIsNan(iomodel->Data(StressbalanceSpcvzEnum)[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1, + iomodel->Data(StressbalanceSpcvzEnum)[i]/yts,StressbalanceVerticalAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + + } + } + } + + /*Free data: */ + iomodel->DeleteData(2,StressbalanceSpcvzEnum,FlowequationBorderFSEnum); + + /*Assign output pointer: */ + *pconstraints=constraints; +}/*}}}*/ +void StressbalanceVerticalAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*No loads*/ + +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 16539) @@ -6,5 +6,99 @@ /*Model processing*/ -int ThermalAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ +int ThermalAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/ return 1; }/*}}}*/ +void ThermalAnalysis::UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ +}/*}}}*/ +void ThermalAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Now, is the model 3d? otherwise, do nothing: */ + if(iomodel->meshtype==Mesh2DhorizontalEnum)return; + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + bool dakota_analysis; + iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum); + + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,FrictionCoefficientEnum); + iomodel->FetchDataToInput(elements,FrictionPEnum); + iomodel->FetchDataToInput(elements,FrictionQEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum); + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + iomodel->FetchDataToInput(elements,FlowequationElementEquationEnum); + iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum); + iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum); + iomodel->FetchDataToInput(elements,PressureEnum); + iomodel->FetchDataToInput(elements,TemperatureEnum); + iomodel->FetchDataToInput(elements,BasalforcingsGeothermalfluxEnum); + iomodel->FetchDataToInput(elements,VxEnum); + iomodel->FetchDataToInput(elements,VyEnum); + iomodel->FetchDataToInput(elements,VzEnum); + InputUpdateFromConstantx(elements,0.,VxMeshEnum); + InputUpdateFromConstantx(elements,0.,VyMeshEnum); + InputUpdateFromConstantx(elements,0.,VzMeshEnum); + if(dakota_analysis){ + elements->InputDuplicate(TemperatureEnum,QmuTemperatureEnum); + elements->InputDuplicate(BasalforcingsMeltingRateEnum,QmuMeltingEnum); + elements->InputDuplicate(VxMeshEnum,QmuVxMeshEnum); + elements->InputDuplicate(VxMeshEnum,QmuVyMeshEnum); + elements->InputDuplicate(VxMeshEnum,QmuVzMeshEnum); + } +}/*}}}*/ +void ThermalAnalysis::CreateNodes(Nodes** pnodes,IoModel* iomodel){/*{{{*/ + + if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); + ::CreateNodes(pnodes,iomodel,ThermalAnalysisEnum,P1Enum); + iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum); +}/*}}}*/ +void ThermalAnalysis::CreateConstraints(Constraints** pconstraints,IoModel* iomodel){/*{{{*/ + + /*Recover pointer: */ + Constraints* constraints=*pconstraints; + + /*Only 3d mesh supported*/ + if(iomodel->meshtype==Mesh3DEnum){ + IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum,P1Enum); + } + + /*Assign output pointer: */ + *pconstraints=constraints; +}/*}}}*/ +void ThermalAnalysis::CreateLoads(Loads** ploads, IoModel* iomodel){/*{{{*/ + + /*Recover pointer: */ + Loads* loads=*ploads; + + if(iomodel->meshtype==Mesh2DhorizontalEnum) _error_("2d meshes not supported yet"); + + /*create penalties for nodes: no node can have a temperature over the melting point*/ + iomodel->FetchData(1,ThermalSpctemperatureEnum); + CreateSingleNodeToElementConnectivity(iomodel); + + for(int i=0;inumberofvertices;i++){ + + /*keep only this partition's nodes:*/ + if(iomodel->my_vertices[i]){ + if (xIsNan(iomodel->Data(ThermalSpctemperatureEnum)[i])){ //No penalty applied on spc nodes! + loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,ThermalAnalysisEnum)); + } + } + } + iomodel->DeleteData(1,ThermalSpctemperatureEnum); + + /*Assign output pointer: */ + *ploads=loads; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.h (revision 16538) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.h (revision 16539) @@ -12,5 +12,10 @@ public: - int DofsPerNode(int** doflist,int meshtype,int approximation); + int DofsPerNode(int** doflist,int meshtype,int approximation); + void UpdateParameters(Parameters** pparameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes** pnodes,IoModel* iomodel); + void CreateConstraints(Constraints** pconstraints,IoModel* iomodel); + void CreateLoads(Loads** ploads, IoModel* iomodel); }; #endif Index: /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp (revision 16539) @@ -36,4 +36,5 @@ /*data: */ + const int RIFTINFOSIZE = 12; int riftfront_node_ids[2]; int riftfront_elem_ids[2]; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 16539) @@ -14,5 +14,5 @@ #include "./ModelProcessorx.h" -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){ +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_enum,const int analysis_enum,const int nummodels,int analysis_counter){ Elements *elements = NULL; @@ -25,6 +25,9 @@ ElementsAndVerticesPartitioning(&iomodel->my_elements,&iomodel->my_vertices,iomodel); - /*Create elements, vertices and materials, independent of analysis_type: */ + /*Create elements, vertices and materials, independent of analysis_enum: */ CreateElementsVerticesAndMaterials(pelements, pvertices, pmaterials, iomodel,nummodels); + + /*Create Parameters*/ + CreateParameters(pparameters,iomodel,rootpath,toolkitfile,solution_enum,analysis_enum,analysis_counter); /*Recover elements and materials, for future update: */ @@ -38,171 +41,11 @@ /*Now, branch onto analysis dependent model generation: */ - switch(analysis_type){ - - #ifdef _HAVE_STRESSBALANCE_ - case StressbalanceAnalysisEnum: - CreateNodesStressbalance(pnodes,iomodel); - CreateConstraintsStressbalance(pconstraints,iomodel); - CreateLoadsStressbalance(ploads,iomodel); - UpdateElementsStressbalance(elements,iomodel,analysis_counter,analysis_type); - break; - - case StressbalanceVerticalAnalysisEnum: - CreateNodesStressbalanceVertical(pnodes, iomodel); - CreateConstraintsStressbalanceVertical(pconstraints,iomodel); - CreateLoadsStressbalanceVertical(ploads,iomodel); - UpdateElementsStressbalanceVertical(elements,iomodel,analysis_counter,analysis_type); - break; - - case StressbalanceSIAAnalysisEnum: - CreateNodesStressbalanceSIA(pnodes, iomodel); - CreateConstraintsStressbalanceSIA(pconstraints,iomodel); - CreateLoadsStressbalanceSIA(ploads,iomodel); - UpdateElementsStressbalanceSIA(elements,iomodel,analysis_counter,analysis_type); - break; - #endif - - #ifdef _HAVE_HYDROLOGY_ - case HydrologyShreveAnalysisEnum: - CreateNodesHydrologyShreve(pnodes, iomodel); - CreateConstraintsHydrologyShreve(pconstraints,iomodel); - CreateLoadsHydrologyShreve(ploads,iomodel); - UpdateElementsHydrologyShreve(elements,iomodel,analysis_counter,analysis_type); - break; - case HydrologyDCInefficientAnalysisEnum: - CreateNodesHydrologyDCInefficient(pnodes, iomodel); - CreateConstraintsHydrologyDCInefficient(pconstraints,iomodel); - CreateLoadsHydrologyDCInefficient(ploads,iomodel); - UpdateElementsHydrologyDCInefficient(elements,iomodel,analysis_counter,analysis_type); - break; - case HydrologyDCEfficientAnalysisEnum: - CreateNodesHydrologyDCEfficient(pnodes, iomodel); - CreateConstraintsHydrologyDCEfficient(pconstraints,iomodel); - CreateLoadsHydrologyDCEfficient(ploads,iomodel); - UpdateElementsHydrologyDCEfficient(elements,iomodel,analysis_counter,analysis_type); - break; - #endif - - #ifdef _HAVE_THERMAL_ - case ThermalAnalysisEnum: - CreateNodesThermal(pnodes, iomodel); - CreateConstraintsThermal(pconstraints,iomodel); - CreateLoadsThermal(ploads,iomodel); - UpdateElementsThermal(elements,iomodel,analysis_counter,analysis_type); - break; - - case EnthalpyAnalysisEnum: - CreateNodesEnthalpy(pnodes, iomodel); - CreateConstraintsEnthalpy(pconstraints,iomodel); - CreateLoadsEnthalpy(ploads,iomodel); - UpdateElementsEnthalpy(elements,iomodel,analysis_counter,analysis_type); - break; - - case MeltingAnalysisEnum: - CreateNodesMelting(pnodes, iomodel); - CreateConstraintsMelting(pconstraints,iomodel); - CreateLoadsMelting(ploads,iomodel); - UpdateElementsMelting(elements,iomodel,analysis_counter,analysis_type); - break; - #endif - - #ifdef _HAVE_BALANCED_ - case BalancethicknessAnalysisEnum: - CreateNodesBalancethickness(pnodes, iomodel); - CreateConstraintsBalancethickness(pconstraints,iomodel); - CreateLoadsBalancethickness(ploads,iomodel); - UpdateElementsBalancethickness(elements,iomodel,analysis_counter,analysis_type); - break; - case BalancevelocityAnalysisEnum: - CreateNodesBalancevelocity(pnodes, iomodel); - CreateConstraintsBalancevelocity(pconstraints,iomodel); - CreateLoadsBalancevelocity(ploads,iomodel); - UpdateElementsBalancevelocity(elements,iomodel,analysis_counter,analysis_type); - break; - case SmoothedSurfaceSlopeXAnalysisEnum: - iomodel->FetchData(1,MeshVertexonbedEnum); - CreateNodes(pnodes,iomodel,SmoothedSurfaceSlopeXAnalysisEnum,P1Enum); - iomodel->DeleteData(1,MeshVertexonbedEnum); - UpdateElementsL2ProjectionBase(elements,iomodel,analysis_counter,analysis_type); - break; - case SmoothedSurfaceSlopeYAnalysisEnum: - iomodel->FetchData(1,MeshVertexonbedEnum); - CreateNodes(pnodes,iomodel,SmoothedSurfaceSlopeYAnalysisEnum,P1Enum); - iomodel->DeleteData(1,MeshVertexonbedEnum); - UpdateElementsL2ProjectionBase(elements,iomodel,analysis_counter,analysis_type); - break; - #endif - - #ifdef _HAVE_GIA_ - case GiaAnalysisEnum: - CreateNodesGia(pnodes, iomodel); - CreateConstraintsGia(pconstraints,iomodel); - CreateLoadsGia(ploads,iomodel); - UpdateElementsGia(elements,iomodel,analysis_counter,analysis_type); - break; - #endif - - #ifdef _HAVE_DAMAGE_ - case DamageEvolutionAnalysisEnum: - CreateNodesDamage(pnodes, iomodel); - CreateConstraintsDamage(pconstraints,iomodel); - CreateLoadsDamage(ploads,iomodel); - UpdateElementsDamage(elements,iomodel,analysis_counter,analysis_type); - break; - #endif - - - #ifdef _HAVE_SLOPE_ - case L2ProjectionBaseAnalysisEnum: - CreateNodesL2ProjectionBase(pnodes, iomodel); - CreateConstraintsL2ProjectionBase(pconstraints,iomodel); - CreateLoadsL2ProjectionBase(ploads,iomodel); - UpdateElementsL2ProjectionBase(elements,iomodel,analysis_counter,analysis_type); - break; - - case L2ProjectionTopAnalysisEnum: - CreateNodesL2ProjectionTop(pnodes, iomodel); - CreateConstraintsL2ProjectionTop(pconstraints,iomodel); - CreateLoadsL2ProjectionTop(ploads,iomodel); - UpdateElementsL2ProjectionTop(elements,iomodel,analysis_counter,analysis_type); - break; - #endif - - #ifdef _HAVE_MASSTRANSPORT_ - case MasstransportAnalysisEnum: - CreateNodesMasstransport(pnodes, iomodel); - CreateConstraintsMasstransport(pconstraints,iomodel); - CreateLoadsMasstransport(ploads,iomodel); - UpdateElementsMasstransport(elements,iomodel,analysis_counter,analysis_type); - break; - case FreeSurfaceTopAnalysisEnum: - CreateNodesFreeSurfaceTop(pnodes, iomodel); - CreateConstraintsFreeSurfaceTop(pconstraints,iomodel); - CreateLoadsFreeSurfaceTop(ploads,iomodel); - UpdateElementsFreeSurfaceTop(elements,iomodel,analysis_counter,analysis_type); - break; - case FreeSurfaceBaseAnalysisEnum: - CreateNodesFreeSurfaceBase(pnodes, iomodel); - CreateConstraintsFreeSurfaceBase(pconstraints,iomodel); - CreateLoadsFreeSurfaceBase(ploads,iomodel); - UpdateElementsFreeSurfaceBase(elements,iomodel,analysis_counter,analysis_type); - break; - case ExtrudeFromBaseAnalysisEnum: - CreateNodesExtrudeFromBase(pnodes, iomodel); - CreateConstraintsExtrudeFromBase(pconstraints,iomodel); - CreateLoadsExtrudeFromBase(ploads,iomodel); - UpdateElementsExtrudeFromBase(elements,iomodel,analysis_counter,analysis_type); - break; - case ExtrudeFromTopAnalysisEnum: - CreateNodesExtrudeFromTop(pnodes, iomodel); - CreateConstraintsExtrudeFromTop(pconstraints,iomodel); - CreateLoadsExtrudeFromTop(ploads,iomodel); - UpdateElementsExtrudeFromTop(elements,iomodel,analysis_counter,analysis_type); - break; - #endif - - default: - _error_("analysis_type: " << EnumToStringx(analysis_type) << " not supported yet!"); - } + Analysis* analysis = EnumToAnalysis(analysis_enum); + analysis->UpdateParameters(pparameters,iomodel,solution_enum,analysis_enum); + analysis->CreateNodes(pnodes,iomodel); + analysis->CreateConstraints(pconstraints,iomodel); + analysis->CreateLoads(ploads,iomodel); + analysis->UpdateElements(elements,iomodel,analysis_counter,analysis_enum); + delete analysis; /*Update Elements and Materials For Inversions*/ @@ -216,12 +59,10 @@ #endif - /*Generate objects that are not dependent on any analysis_type: */ - CreateParameters(pparameters,iomodel,rootpath,toolkitfile,solution_type,analysis_type,analysis_counter); /*Update Elements in case we are running a transient solution: */ #ifdef _HAVE_TRANSIENT_ parameters=*pparameters; - if(analysis_counter==(nummodels-1)&& solution_type==TransientSolutionEnum){ - UpdateElementsTransient(elements,parameters,iomodel,analysis_counter,analysis_type); + if(analysis_counter==(nummodels-1)&& solution_enum==TransientSolutionEnum){ + UpdateElementsTransient(elements,parameters,iomodel,analysis_counter,analysis_enum); } #endif Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 16539) @@ -25,5 +25,5 @@ bool ispdd,isdelta18o; - /*parameters for mass flux: {{{*/ + /*parameters for mass flux:*/ int mass_flux_num_profiles = 0; bool qmu_mass_flux_present = false; @@ -37,5 +37,4 @@ IssmDouble *matrix = NULL; int count; - /*}}}*/ if(*pparameters)return; //do not create parameters twice! @@ -94,4 +93,10 @@ /*For stress balance only*/ + parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum)); if(iomodel->meshtype==Mesh3DEnum) parameters->AddObject(iomodel->CopyConstantObject(MeshNumberoflayersEnum)); @@ -221,15 +226,4 @@ CreateOutputDefinitions(¶meters,iomodel); - /*Solution specific parameters*/ - #ifdef _HAVE_HYDROLOGY_ - CreateParametersHydrologyShreve(¶meters,iomodel,solution_type,analysis_type); - CreateParametersHydrologyDCInefficient(¶meters,iomodel,solution_type,analysis_type); - CreateParametersHydrologyDCEfficient(¶meters,iomodel,solution_type,analysis_type); - #endif - - #ifdef _HAVE_DAMAGE_ - CreateParametersDamage(¶meters,iomodel,solution_type,analysis_type); - #endif - /*Before returning, create parameters in case we are running Qmu or control types runs: */ #ifdef _HAVE_CONTROL_ @@ -241,9 +235,4 @@ #endif - #ifdef _HAVE_STRESSBALANCE_ - CreateParametersStressbalance(¶meters,iomodel,solution_type,analysis_type); - #endif - - /*Now, deal with toolkits options, which need to be put into the parameters dataset: */ ParseToolkitsOptionsx(parameters,toolkitsoptionsfid); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp (revision 16538) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp (revision 16539) @@ -19,6 +19,5 @@ int i; - int my_rank; - int num_procs; + const int RIFTINFOSIZE = 12; int numberofelements2d; int numberofvertices2d; @@ -41,13 +40,12 @@ /*Get my_rank:*/ - my_rank = IssmComm::GetRank(); - num_procs = IssmComm::GetSize(); + int my_rank = IssmComm::GetRank(); + int num_procs = IssmComm::GetSize(); /*Fetch parameters: */ - iomodel->Constant(&numrifts,RiftsNumriftsEnum); /*First, check that partitioning has not yet been carryed out. Just check whether my_elements pointers is not already assigned a value: */ - if (*pmy_elements)return; + if(*pmy_elements)return; /*Number of vertices per elements, needed to correctly retrieve data: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 16538) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 16539) @@ -6,7 +6,6 @@ #define _MODEL_PROCESSORX_H_ -#define RIFTINFOSIZE 12 - #include "../../classes/classes.h" +#include "../../analyses/analyses.h" void 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); @@ -19,136 +18,9 @@ void CreateParametersControl(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); void CreateParametersDakota(Parameters** pparameters,IoModel* iomodel,char* rootpath,int solution_type,int analysis_type); -void CreateParametersStressbalance(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); -void CreateParametersHydrologyShreve(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); -void CreateParametersHydrologyDCInefficient(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); -void CreateParametersHydrologyDCEfficient(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); void CreateOutputDefinitions(Parameters** pparameters,IoModel* iomodel); void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel); void UpdateElementsAndMaterialsDakota(Elements* elements,Materials* materials, IoModel* iomodel); +void UpdateElementsTransient(Elements* elements,Parameters* parameters,IoModel* iomodel,int analysis_counter,int analysis_type); void CreateNodes(Nodes** pnodes, IoModel* iomodel,int analysis,int finite_element,int approximation=NoneApproximationEnum); - -/*Creation of fem datasets: specialised drivers: */ - -/*stressbalance horizontal*/ -void CreateNodesStressbalance(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsStressbalance(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsStressbalance(Loads** ploads, IoModel* iomodel); -void UpdateElementsStressbalance(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -/*stressbalance vertical*/ -void CreateNodesStressbalanceVertical(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsStressbalanceVertical(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsStressbalanceVertical(Loads** ploads, IoModel* iomodel); -void UpdateElementsStressbalanceVertical(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -/*stressbalance SIA*/ -void CreateNodesStressbalanceSIA(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsStressbalanceSIA(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsStressbalanceSIA(Loads** ploads, IoModel* iomodel); -void UpdateElementsStressbalanceSIA(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -#ifdef _HAVE_GIA_ -/*gia*/ -void CreateNodesGia(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsGia(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsGia(Loads** ploads, IoModel* iomodel); -void UpdateElementsGia(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); -#endif - -#ifdef _HAVE_DAMAGE_ -/*damage*/ -void CreateNodesDamage(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsDamage(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsDamage(Loads** ploads, IoModel* iomodel); -void UpdateElementsDamage(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); -void CreateParametersDamage(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); -#endif - -/*Extrude from base*/ -void CreateNodesExtrudeFromBase(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsExtrudeFromBase(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsExtrudeFromBase(Loads** ploads, IoModel* iomodel); -void UpdateElementsExtrudeFromBase(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -/*Extrude from top*/ -void CreateNodesExtrudeFromTop(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsExtrudeFromTop(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsExtrudeFromTop(Loads** ploads, IoModel* iomodel); -void UpdateElementsExtrudeFromTop(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -/*L2 projection base*/ -void CreateNodesL2ProjectionBase(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsL2ProjectionBase(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsL2ProjectionBase(Loads** ploads, IoModel* iomodel); -void UpdateElementsL2ProjectionBase(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -/*L2 projection top*/ -void CreateNodesL2ProjectionTop(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsL2ProjectionTop(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsL2ProjectionTop(Loads** ploads, IoModel* iomodel); -void UpdateElementsL2ProjectionTop(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -/*thermal:*/ -void CreateNodesThermal(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsThermal(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsThermal(Loads** ploads, IoModel* iomodel); -void UpdateElementsThermal(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -/*enthalpy:*/ -void CreateNodesEnthalpy(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsEnthalpy(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsEnthalpy(Loads** ploads, IoModel* iomodel); -void UpdateElementsEnthalpy(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -/*hydrology Shreve:*/ -void CreateNodesHydrologyShreve(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsHydrologyShreve(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsHydrologyShreve(Loads** ploads, IoModel* iomodel); -void UpdateElementsHydrologyShreve(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -/*hydrology DC:*/ -void CreateNodesHydrologyDCInefficient(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsHydrologyDCInefficient(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsHydrologyDCInefficient(Loads** ploads, IoModel* iomodel); -void UpdateElementsHydrologyDCInefficient(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); -void CreateNodesHydrologyDCEfficient(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsHydrologyDCEfficient(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsHydrologyDCEfficient(Loads** ploads, IoModel* iomodel); -void UpdateElementsHydrologyDCEfficient(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -/*melting:*/ -void CreateNodesMelting(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsMelting(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsMelting(Loads** ploads, IoModel* iomodel); -void UpdateElementsMelting(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -/*masstransport:*/ -void CreateNodesMasstransport(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsMasstransport(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsMasstransport(Loads** ploads, IoModel* iomodel); -void UpdateElementsMasstransport(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); -void CreateNodesFreeSurfaceTop(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsFreeSurfaceTop(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsFreeSurfaceTop(Loads** ploads, IoModel* iomodel); -void UpdateElementsFreeSurfaceTop(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); -void CreateNodesFreeSurfaceBase(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsFreeSurfaceBase(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsFreeSurfaceBase(Loads** ploads, IoModel* iomodel); -void UpdateElementsFreeSurfaceBase(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -/*balancedthickness:*/ -void CreateNodesBalancethickness(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsBalancethickness(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsBalancethickness(Loads** ploads, IoModel* iomodel); -void UpdateElementsBalancethickness(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -/*balancedvelocity:*/ -void CreateNodesBalancevelocity(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsBalancevelocity(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsBalancevelocity(Loads** ploads, IoModel* iomodel); -void UpdateElementsBalancevelocity(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); - -/*transient: */ -void UpdateElementsTransient(Elements* elements,Parameters* parameters,IoModel* iomodel,int analysis_counter,int analysis_type); /*partitioning: */