Changeset 27909


Ignore:
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:
2 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r27852 r27909  
    14511451IssmDouble Penta::GetIcefrontArea(){/*{{{*/
    14521452
    1453         IssmDouble  bed[NUMVERTICES]; //basinId[NUMVERTICES];
    1454         IssmDouble      Haverage,frontarea;
    1455         IssmDouble  x1,y1,x2,y2,distance;
    1456         IssmDouble lsf[NUMVERTICES], Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
    1457         int* indices=NULL;
    1458         IssmDouble* H=NULL;;
    1459         int nrfrontbed,numiceverts;
    1460 
     1453        /*We need to be on base and cross the levelset*/
     1454        if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
    14611455        if(!this->IsOnBase()) return 0;
    1462         if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
    1463 
    1464         /*Retrieve all inputs and parameters*/
    1465         Element::GetInputListOnVertices(&bed[0],BedEnum);
    1466         Element::GetInputListOnVertices(&surfaces[0],SurfaceEnum);
    1467         Element::GetInputListOnVertices(&bases[0],BaseEnum);
    1468         Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
    1469 
    1470         nrfrontbed=0;
    1471         for(int i=0;i<NUMVERTICES2D;i++){
    1472                 /*Find if bed<0*/
    1473                 if(bed[i]<0.) nrfrontbed++;
    1474         }
    1475 
    1476         if(nrfrontbed==3){
    1477                 /*2. Find coordinates of where levelset crosses 0*/
    1478                 int         numiceverts;
    1479                 IssmDouble  s[2],x[2],y[2];
    1480                 this->GetLevelsetIntersectionBase(&indices, &numiceverts,&s[0],MaskIceLevelsetEnum,0.);
    1481                 _assert_(numiceverts);
    1482 
    1483                 /*3 Write coordinates*/
    1484                 IssmDouble  xyz_list[NUMVERTICES][3];
    1485                 ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
    1486                 int counter = 0;
    1487                 if((numiceverts>0) && (numiceverts<NUMVERTICES2D)){
    1488                         for(int i=0;i<numiceverts;i++){
    1489                                 for(int n=numiceverts;n<NUMVERTICES2D;n++){ // iterate over no-ice vertices
    1490                                         x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
    1491                                         y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
    1492                                         counter++;
    1493                                 }
    1494                         }
    1495                 }
    1496                 else if(numiceverts==NUMVERTICES2D){ //NUMVERTICES ice vertices: calving front lies on element edge
    1497 
    1498                         for(int i=0;i<NUMVERTICES2D;i++){
    1499                                 if(lsf[indices[i]]==0.){
    1500                                         x[counter]=xyz_list[indices[i]][0];
    1501                                         y[counter]=xyz_list[indices[i]][1];
    1502                                         counter++;
    1503                                 }
    1504                                 if(counter==2) break;
    1505                         }
    1506                         if(counter==1){
    1507                                 /*We actually have only 1 vertex on levelset, write a single point as a segment*/
    1508                                 x[counter]=x[0];
    1509                                 y[counter]=y[0];
    1510                                 counter++;
    1511                         }
    1512                 }
    1513                 else{
    1514                         _error_("not sure what's going on here...");
    1515                 }
    1516                 x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
    1517                 distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));
    1518 
    1519                 int numthk=numiceverts+2;
    1520                 H=xNew<IssmDouble>(numthk);
    1521                 for(int iv=0;iv<NUMVERTICES2D;iv++) Haux[iv]=-bed[indices[iv]]; //sort bed in ice/noice
    1522 
    1523                 switch(numiceverts){
    1524                         case 1: // average over triangle
    1525                                 H[0]=Haux[0];
    1526                                 H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
    1527                                 H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
    1528                                 Haverage=(H[1]+H[2])/2;
    1529                                 break;
    1530                         case 2: // average over quadrangle
    1531                                 H[0]=Haux[0];
    1532                                 H[1]=Haux[1];
    1533                                 H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
    1534                                 H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
    1535                                 Haverage=(H[2]+H[3])/2;
    1536                                 break;
    1537                         default:
    1538                                 _error_("Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
    1539                                 break;
    1540                 }
    1541                 frontarea=distance*Haverage;
    1542         }
    1543         else return 0;
    1544 
    1545         xDelete<int>(indices);
    1546         xDelete<IssmDouble>(H);
    1547 
    1548         _assert_(frontarea>0);
     1456
     1457        /*Spawn Tria element from the base of the Penta: */
     1458        Tria* tria=(Tria*)SpawnTria(0,1,2);
     1459        IssmDouble frontarea = tria->GetIcefrontArea();
     1460        delete tria->material; delete tria;
     1461
    15491462        return frontarea;
    1550 }
    1551 /*}}}*/
     1463}/*}}}*/
    15521464void       Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
    15531465
  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r27907 r27909  
    27792779
    27802780        /*Return if no ice front present*/
    2781         if(!this->IsIcefront()) return 0.;
     2781        if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
     2782        //if(!this->IsIcefront()) return 0.;
    27822783
    27832784        /*Retrieve all inputs and parameters*/
Note: See TracChangeset for help on using the changeset viewer.