Changeset 16539 for issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

Timestamp:

10/24/13 10:12:44 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: moving model processor stuff to analyses

File:

• issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (modified) (1 diff)

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

@@ -6 +6 @@
 
 /*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;i<iomodel->numberofelements;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;i<iomodel->numberoffaces;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;i<numvertex_pairing;i++){
+
+                if(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+0])-1]){
+
+                        /*In debugging mode, check that the second node is in the same cpu*/
+                        _assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);
+
+                        /*Skip if one of the two is not on the bed*/
+                        if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
+
+                        /*Get node ids*/
+                        penpair_ids[0]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+0]);
+                        penpair_ids[1]=iomodel->nodecounter+reCast<int>(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;
+}/*}}}*/