source: issm/oecreview/Archive/23390-24306/ISSM-23800-23801.diff@ 24307

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

@@ -74 +74 @@
                 int            GetElementType(void);
                 void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
                 IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
+                IssmDouble              GetIcefrontArea();
                 void           GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
                 void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
                 void           GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented yet");};
                 void           GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues){_error_("not implemented yet");};
+                void                            GetLevelsetIntersectionBase(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level);
                 Node*          GetNode(int node_number);
                 int            GetNodeIndex(Node* node);
                 int            GetNumberOfNodes(void);
@@ -149 +151 @@
                 void           ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
                 void           ResetFSBasalBoundaryCondition(void);
                 void           ResetHooks();
+                void                            RignotMeltParameterization();
                 void           SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int N,int M);
                 void           SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
                 void           SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);

../trunk-jpl/src/c/classes/FemModel.cpp

@@ -1497 +1497 @@
 }/*}}}*/
 void FemModel::IcefrontAreax(){/*{{{*/
 
-        int numvertices      = this->GetElementsWidth();
         int numbasins;
-        IssmDouble* BasinId   = xNew<IssmDouble>(numvertices);
         this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
         IssmDouble* basin_icefront_area           = xNewZeroInit<IssmDouble>(numbasins);
 
@@ -1509 +1507 @@
 
                 for(int i=0;i<this->elements->Size();i++){
                         Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+                        if(!element->IsOnBase()) continue;
+                        int  numvertices = element->GetNumberOfVertices();
+                        IssmDouble* BasinId   = xNew<IssmDouble>(numvertices);
                         element->GetInputListOnVertices(BasinId,FrontalForcingsBasinIdEnum);
-                        for(int j=0;j<numvertices;j++){
+                        for(int j=0;j<3;j++){
                                 if(BasinId[j]==basin){
                                         local_icefront_area+=element->GetIcefrontArea();
                                         break;
                                 }
                         }
+                        xDelete<IssmDouble>(BasinId);
                 }
                 ISSM_MPI_Reduce(&local_icefront_area,&total_icefront_area,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
                 ISSM_MPI_Bcast(&total_icefront_area,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        
+                
                 basin_icefront_area[basin-1]=total_icefront_area;
         }
         
         this->parameters->AddObject(new DoubleVecParam(FrontalForcingsBasinIcefrontAreaEnum,basin_icefront_area,numbasins));
+        xDelete<IssmDouble>(basin_icefront_area);
         
-        xDelete<IssmDouble>(basin_icefront_area);
-        xDelete<IssmDouble>(BasinId);
 }/*}}}*/
 void FemModel::IceMassx(IssmDouble* pM, bool scaled){/*{{{*/
 

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

@@ -1069 +1069 @@
         return phi;
 }
 /*}}}*/
+IssmDouble Penta::GetIcefrontArea(){/*{{{*/
+        
+        IssmDouble  bed[NUMVERTICES]; //basinId[NUMVERTICES];
+        IssmDouble      Haverage,frontarea;
+        IssmDouble  x1,y1,x2,y2,distance;
+        IssmDouble lsf[NUMVERTICES], Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
+        int* indices=NULL;
+        IssmDouble* H=NULL;;
+        int nrfrontbed,numiceverts;
+
+        /*Retrieve all inputs and parameters*/
+        GetInputListOnVertices(&bed[0],BedEnum);
+        GetInputListOnVertices(&surfaces[0],SurfaceEnum);
+        GetInputListOnVertices(&bases[0],BaseEnum);
+        GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+
+        if(!this->IsOnBase()) return 0;
+        if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
+
+        nrfrontbed=0;
+        for(int i=0;i<NUMVERTICES2D;i++){
+                /*Find if bed<0*/
+                if(bed[i]<0.) nrfrontbed++;
+        }
+
+        if(nrfrontbed==3){
+                /*2. Find coordinates of where levelset crosses 0*/
+                int         numiceverts;
+                IssmDouble  s[2],x[2],y[2];
+                this->GetLevelsetIntersectionBase(&indices, &numiceverts,&s[0],MaskIceLevelsetEnum,0.);
+                _assert_(numiceverts); 
+
+                /*3 Write coordinates*/
+                IssmDouble  xyz_list[NUMVERTICES][3];
+                ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
+                int counter = 0;
+                if((numiceverts>0) && (numiceverts<NUMVERTICES2D)){
+                        for(int i=0;i<numiceverts;i++){
+                                for(int n=numiceverts;n<NUMVERTICES2D;n++){ // iterate over no-ice vertices
+                                        x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
+                                        y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
+                                        counter++;
+                                }
+                        }
+                }
+                else if(numiceverts==NUMVERTICES2D){ //NUMVERTICES ice vertices: calving front lies on element edge
+
+                        for(int i=0;i<NUMVERTICES2D;i++){
+                                if(lsf[indices[i]]==0.){
+                                        x[counter]=xyz_list[indices[i]][0];
+                                        y[counter]=xyz_list[indices[i]][1];
+                                        counter++;
+                                }
+                                if(counter==2) break;
+                        }
+                        if(counter==1){
+                                /*We actually have only 1 vertex on levelset, write a single point as a segment*/
+                                x[counter]=x[0];
+                                y[counter]=y[0];
+                                counter++;
+                        }
+                }
+                else{
+                        _error_("not sure what's going on here...");
+                }
+                x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
+                distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));
+                
+                int numthk=numiceverts+2;
+                H=xNew<IssmDouble>(numthk);
+                for(int iv=0;iv<NUMVERTICES2D;iv++) Haux[iv]=-bed[indices[iv]]; //sort bed in ice/noice
+
+                switch(numiceverts){
+                        case 1: // average over triangle
+                                H[0]=Haux[0];
+                                H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
+                                H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
+                                Haverage=(H[1]+H[2])/2;
+                                break;
+                        case 2: // average over quadrangle
+                                H[0]=Haux[0];
+                                H[1]=Haux[1];
+                                H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
+                                H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
+                                Haverage=(H[2]+H[3])/2;
+                                break;
+                        default:
+                                _error_("Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
+                                break;
+                }
+                frontarea=distance*Haverage;
+        }
+        else return 0;  
+        
+        xDelete<int>(indices);
+        xDelete<IssmDouble>(H);
+        
+        _assert_(frontarea>0);
+        return frontarea;
+}
+/*}}}*/
 void       Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
 
         /* Intermediaries */
@@ -1132 +1233 @@
         return lower_penta;
 }
 /*}}}*/
+void       Penta::GetLevelsetIntersectionBase(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level){/*{{{*/
+
+        /* GetLevelsetIntersection computes: 
+         * 1. indices of element, sorted in [iceverts, noiceverts] in counterclockwise fashion,
+         * 2. fraction of intersected triangle edges intersected by levelset, lying below level*/
+
+        /*Intermediaries*/
+        int i, numiceverts, numnoiceverts;
+        int ind0, ind1, lastindex;
+        int indices_ice[NUMVERTICES2D],indices_noice[NUMVERTICES2D];
+        IssmDouble lsf[NUMVERTICES];
+        int* indices = xNew<int>(NUMVERTICES2D);
+
+        /*Retrieve all inputs and parameters*/
+        GetInputListOnVertices(&lsf[0],levelset_enum);
+
+        /* Determine distribution of ice over element.
+         * Exploit: ice/no-ice parts are connected, so find starting vertex of segment*/
+        lastindex=0;
+        for(i=0;i<NUMVERTICES2D;i++){ // go backwards along vertices, and check for sign change
+                ind0=(NUMVERTICES2D-i)%NUMVERTICES2D;
+                ind1=(ind0-1+NUMVERTICES2D)%NUMVERTICES2D;
+                if((lsf[ind0]-level)*(lsf[ind1]-level)<=0.){ // levelset has been crossed, find last index belonging to segment
+                        if(lsf[ind1]==level) //if levelset intersects 2nd vertex, choose this vertex as last
+                                lastindex=ind1;
+                        else
+                                lastindex=ind0;
+                        break;
+                }
+        }
+
+        numiceverts=0;
+        numnoiceverts=0;
+        for(i=0;i<NUMVERTICES2D;i++){
+                ind0=(lastindex+i)%NUMVERTICES2D;
+                if(lsf[i]<=level){
+                        indices_ice[numiceverts]=i;
+                        numiceverts++;
+                }
+                else{
+                        indices_noice[numnoiceverts]=i;
+                        numnoiceverts++;
+                }
+        }
+        //merge indices 
+        for(i=0;i<numiceverts;i++){indices[i]=indices_ice[i];}
+        for(i=0;i<numnoiceverts;i++){indices[numiceverts+i]=indices_noice[i];}
+
+        switch (numiceverts){
+                case 0: // no vertex has ice: element is ice free, no intersection
+                        for(i=0;i<2;i++)
+                                fraction[i]=0.;
+                        break;
+                case 1: // one vertex has ice:
+                        for(i=0;i<2;i++){
+                                fraction[i]=(level-lsf[indices[0]])/(lsf[indices[numiceverts+i]]-lsf[indices[0]]);
+                        }
+                        break;
+                case 2: // two vertices have ice: fraction is computed from first ice vertex to last in CCW fashion
+                        for(i=0;i<2;i++){
+                                fraction[i]=(level-lsf[indices[i]])/(lsf[indices[numiceverts]]-lsf[indices[i]]);
+                        }
+                        break;
+                case NUMVERTICES2D: // all vertices have ice: return triangle area
+                        for(i=0;i<2;i++)
+                                fraction[i]=1.;
+                        break;
+                default:
+                        _error_("Wrong number of ice vertices in Penta::GetLevelsetIntersectionBase!");
+                        break;
+        }
+
+        *pindices=indices;
+        *pnumiceverts=numiceverts;
+}
+/*}}}*/
 Node*      Penta::GetNode(int node_number){/*{{{*/
         _assert_(node_number>=0); 
         _assert_(node_number<this->NumberofNodes(this->element_type)); 
@@ -2311 +2488 @@
 
 }
 /*}}}*/
+void       Penta::RignotMeltParameterization(){/*{{{*/
+
+   IssmDouble A, B, alpha, beta;
+        IssmDouble bed,qsg,qsg_basin,TF,yts;
+        int numbasins;
+        IssmDouble basinid[NUMVERTICES];
+        IssmDouble* basin_icefront_area=NULL;
+
+        /* Coefficients */
+        A    = 3e-4;        
+        B    = 0.15;        
+        alpha = 0.39;
+        beta = 1.18;
+        
+        /*Get inputs*/
+        Input* bed_input = this->GetInput(BedEnum);                     _assert_(bed_input);
+        Input* qsg_input = this->GetInput(FrontalForcingsSubglacialDischargeEnum);               _assert_(qsg_input);
+        Input* TF_input  = this->GetInput(FrontalForcingsThermalForcingEnum);          _assert_(TF_input);
+        GetInputListOnVertices(&basinid[0],FrontalForcingsBasinIdEnum); 
+        
+        this->FindParam(&yts, ConstantsYtsEnum);
+        this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
+        this->parameters->FindParam(&basin_icefront_area,&numbasins,FrontalForcingsBasinIcefrontAreaEnum);
+
+        IssmDouble meltrates[NUMVERTICES2D];  //frontal melt-rate
+        
+        /* Start looping on the number of vertices: */
+        GaussPenta* gauss=new GaussPenta();
+        for(int iv=0;iv<NUMVERTICES2D;iv++){
+                gauss->GaussVertex(iv);
+
+                /* Get variables */
+                bed_input->GetInputValue(&bed,gauss);
+                qsg_input->GetInputValue(&qsg,gauss);
+                TF_input->GetInputValue(&TF,gauss);
+
+                if(basin_icefront_area[reCast<int>(basinid[iv])-1]==0.) meltrates[iv]=0.;
+                else{
+                        /* change the unit of qsg (m^3/d -> m/d) with ice front area */
+                        qsg_basin=qsg/basin_icefront_area[reCast<int>(basinid[iv])-1];
+
+                        /* calculate melt rates */
+                        meltrates[iv]=(A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta)/86400; //[m/s]
+                }       
+
+                if(xIsNan<IssmDouble>(meltrates[iv])) _error_("NaN found in vector");
+                if(xIsInf<IssmDouble>(meltrates[iv])) _error_("Inf found in vector");
+        }
+
+        /*Add input*/
+        this->inputs->AddInput(new PentaInput(CalvingMeltingrateEnum,&meltrates[0],P1Enum));
+        
+        this->InputExtrude(CalvingMeltingrateEnum,-1);
+   
+        /*Cleanup and return*/
+        xDelete<IssmDouble>(basin_icefront_area);
+        delete gauss;
+}
+/*}}}*/
 void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int N, int M){/*{{{*/
 
         IssmDouble  values[NUMVERTICES];

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

@@ -3419 +3419 @@
         IssmDouble* basin_icefront_area=NULL;
 
         /* Coefficients */
-        A    = 3e-4;        //
-        B    = 0.15;        //
+        A    = 3e-4;        
+        B    = 0.15;      
         alpha = 0.39;
         beta = 1.18;
         
@@ -3457 +3457 @@
 
                 if(xIsNan<IssmDouble>(meltrates[iv])) _error_("NaN found in vector");
                 if(xIsInf<IssmDouble>(meltrates[iv])) _error_("Inf found in vector");
-                
         }
 
         /*Add input*/