Changeset 27909
- Timestamp:
- 09/20/23 05:27:41 (18 months ago)
- 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 1451 1451 IssmDouble Penta::GetIcefrontArea(){/*{{{*/ 1452 1452 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; 1461 1455 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 1549 1462 return frontarea; 1550 } 1551 /*}}}*/ 1463 }/*}}}*/ 1552 1464 void Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/ 1553 1465 -
issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
r27907 r27909 2779 2779 2780 2780 /*Return if no ice front present*/ 2781 if(!this->IsIcefront()) return 0.; 2781 if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0; 2782 //if(!this->IsIcefront()) return 0.; 2782 2783 2783 2784 /*Retrieve all inputs and parameters*/
Note:
See TracChangeset
for help on using the changeset viewer.