Changeset 27909

Timestamp:

09/20/23 05:27:41 (18 months ago)

Author:

Mathieu Morlighem

Message:

CHG: fixing NR

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• Penta.cpp (modified) (1 diff)
• Tria.cpp (modified) (1 diff)

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

@@ -1451 +1451 @@
 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;
-
+        /*We need to be on base and cross the levelset*/
+        if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
         if(!this->IsOnBase()) return 0;
-        if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
-
-        /*Retrieve all inputs and parameters*/
-        Element::GetInputListOnVertices(&bed[0],BedEnum);
-        Element::GetInputListOnVertices(&surfaces[0],SurfaceEnum);
-        Element::GetInputListOnVertices(&bases[0],BaseEnum);
-        Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
-
-        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);
+
+        /*Spawn Tria element from the base of the Penta: */
+        Tria* tria=(Tria*)SpawnTria(0,1,2);
+        IssmDouble frontarea = tria->GetIcefrontArea();
+        delete tria->material; delete tria;
+
         return frontarea;
-}
-/*}}}*/
+}/*}}}*/
 void       Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
 

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

@@ -2779 +2779 @@
 
         /*Return if no ice front present*/
-        if(!this->IsIcefront()) return 0.;
+        if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
+        //if(!this->IsIcefront()) return 0.;
 
         /*Retrieve all inputs and parameters*/