Changeset 15538


Ignore:
Timestamp:
07/22/13 14:03:58 (12 years ago)
Author:
seroussi
Message:

NEW: added front in element for Pattyn and Stokes

Location:
issm/trunk-jpl/src/c
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r15535 r15538  
    784784}
    785785/*}}}*/
     786/*FUNCTION Penta::GetAreaCoordinates{{{*/
     787void Penta::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble xyz_zero[4][3],IssmDouble xyz_list[6][3],int numpoints){
     788        /*Computeportion of the element that is grounded*/
     789
     790        int         i,j,k;
     791        IssmDouble  area_init,area_portion;
     792        IssmDouble  xyz_bis[3][3];
     793
     794        area_init=fabs(xyz_list[1][0]*xyz_list[2][1] - xyz_list[1][1]*xyz_list[2][0] + xyz_list[0][0]*xyz_list[1][1] - xyz_list[0][1]*xyz_list[1][0] + xyz_list[2][0]*xyz_list[0][1] - xyz_list[2][1]*xyz_list[0][0])/2.;
     795
     796
     797        /*Initialize xyz_list with original xyz_list of triangle coordinates*/
     798        for(j=0;j<3;j++){
     799                for(k=0;k<3;k++){
     800                        xyz_bis[j][k]=xyz_list[j][k];
     801                }
     802        }
     803        for(i=0;i<numpoints;i++){
     804                for(j=0;j<3;j++){
     805                        for(k=0;k<3;k++){
     806                                /*Change appropriate line*/
     807                                xyz_bis[j][k]=xyz_zero[i][k];
     808                        }
     809
     810                        /*Compute area fraction*/
     811                        area_portion=fabs(xyz_bis[1][0]*xyz_bis[2][1] - xyz_bis[1][1]*xyz_bis[2][0] + xyz_bis[0][0]*xyz_bis[1][1] - xyz_bis[0][1]*xyz_bis[1][0] + xyz_bis[2][0]*xyz_bis[0][1] - xyz_bis[2][1]*xyz_bis[0][0])/2.;
     812                        *(area_coordinates+3*i+j)=area_portion/area_init;
     813
     814                        /*Reinitialize xyz_list*/
     815                        for(k=0;k<3;k++){
     816                                /*Reinitialize xyz_list with original coordinates*/
     817                                xyz_bis[j][k]=xyz_list[j][k];
     818                        }
     819                }
     820        }
     821}
     822/*}}}*/
    786823/*FUNCTION Penta::GetBasalElement{{{*/
    787824Penta* Penta::GetBasalElement(void){
    … …  
    11161153         *    = 4 * mu * eps_eff ^2*/
    11171154        *phi=4*pow(epsilon_eff,2.0)*viscosity;
     1155}
     1156/*}}}*/
     1157/*FUNCTION Penta::GetQuadNormal {{{*/
     1158void Penta:: GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){
     1159
     1160        /*Build unit outward pointing vector*/
     1161        IssmDouble AB[3];
     1162        IssmDouble AC[3];
     1163        IssmDouble norm;
     1164
     1165        AB[0]=xyz_list[1][0] - xyz_list[0][0];
     1166        AB[1]=xyz_list[1][1] - xyz_list[0][1];
     1167        AB[2]=xyz_list[1][2] - xyz_list[0][2];
     1168        AC[0]=xyz_list[2][0] - xyz_list[0][0];
     1169        AC[1]=xyz_list[2][1] - xyz_list[0][1];
     1170        AC[2]=xyz_list[2][2] - xyz_list[0][2];
     1171
     1172        cross(normal,AB,AC);
     1173        norm=sqrt(pow(normal[0],2.0)+pow(normal[1],2.0)+pow(normal[2],2.0));
     1174
     1175        for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
    11181176}
    11191177/*}}}*/
    … …  
    13111369
    13121370        return z;
     1371}
     1372/*}}}*/
     1373/*FUNCTION Penta::GetZeroLevelsetCoordinates{{{*/
     1374void Penta::GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[6][3],int levelsetenum){
     1375        /*Computeportion of the element that is grounded*/
     1376
     1377        int         normal_orientation;
     1378        IssmDouble  s1,s2;
     1379        IssmDouble  levelset[3];
     1380
     1381        /*Recover parameters and values*/
     1382        GetInputListOnVertices(&levelset[0],levelsetenum);
     1383
     1384        if(levelset[0]*levelset[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     1385                /*Portion of the segments*/
     1386                s1=levelset[2]/(levelset[2]-levelset[1]);
     1387                s2=levelset[2]/(levelset[2]-levelset[0]);
     1388
     1389                if(levelset[2]>0) normal_orientation=1;
     1390                /*New point 1*/
     1391                xyz_zero[3*normal_orientation+0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
     1392                xyz_zero[3*normal_orientation+1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]);
     1393                xyz_zero[3*normal_orientation+2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]);
     1394
     1395                /*New point 0*/
     1396                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]);
     1397                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]);
     1398                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]);
     1399
     1400                /*New point 3*/
     1401                xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[5][0]+s1*(xyz_list[4][0]-xyz_list[5][0]);
     1402                xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[5][1]+s1*(xyz_list[4][1]-xyz_list[5][1]);
     1403                xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[5][2]+s1*(xyz_list[4][2]-xyz_list[5][2]);
     1404
     1405                /*New point 4*/
     1406                xyz_zero[3*(2+normal_orientation)+0]=xyz_list[5][0]+s2*(xyz_list[3][0]-xyz_list[5][0]);
     1407                xyz_zero[3*(2+normal_orientation)+1]=xyz_list[5][1]+s2*(xyz_list[3][1]-xyz_list[5][1]);
     1408                xyz_zero[3*(2+normal_orientation)+2]=xyz_list[5][2]+s2*(xyz_list[3][2]-xyz_list[5][2]);
     1409        }
     1410        else if(levelset[1]*levelset[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
     1411                /*Portion of the segments*/
     1412                s1=levelset[0]/(levelset[0]-levelset[2]);
     1413                s2=levelset[0]/(levelset[0]-levelset[1]);
     1414
     1415                if(levelset[0]>0) normal_orientation=1;
     1416                /*New point 1*/
     1417                xyz_zero[3*normal_orientation+0]=xyz_list[0][0]+s1*(xyz_list[2][0]-xyz_list[0][0]);
     1418                xyz_zero[3*normal_orientation+1]=xyz_list[0][1]+s1*(xyz_list[2][1]-xyz_list[0][1]);
     1419                xyz_zero[3*normal_orientation+2]=xyz_list[0][2]+s1*(xyz_list[2][2]-xyz_list[0][2]);
     1420
     1421                /*New point 2*/
     1422                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0][0]+s2*(xyz_list[1][0]-xyz_list[0][0]);
     1423                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0][1]+s2*(xyz_list[1][1]-xyz_list[0][1]);
     1424                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0][2]+s2*(xyz_list[1][2]-xyz_list[0][2]);
     1425
     1426                /*New point 3*/
     1427                xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[3][0]+s1*(xyz_list[5][0]-xyz_list[3][0]);
     1428                xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[3][1]+s1*(xyz_list[5][1]-xyz_list[3][1]);
     1429                xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[3][2]+s1*(xyz_list[5][2]-xyz_list[3][2]);
     1430
     1431                /*New point 4*/
     1432                xyz_zero[3*(2+normal_orientation)+0]=xyz_list[3][0]+s2*(xyz_list[4][0]-xyz_list[3][0]);
     1433                xyz_zero[3*(2+normal_orientation)+1]=xyz_list[3][1]+s2*(xyz_list[4][1]-xyz_list[3][1]);
     1434                xyz_zero[3*(2+normal_orientation)+2]=xyz_list[3][2]+s2*(xyz_list[4][2]-xyz_list[3][2]);
     1435        }
     1436        else if(levelset[0]*levelset[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
     1437                /*Portion of the segments*/
     1438                s1=levelset[1]/(levelset[1]-levelset[0]);
     1439                s2=levelset[1]/(levelset[1]-levelset[2]);
     1440
     1441                if(levelset[1]>0) normal_orientation=1;
     1442                /*New point 0*/
     1443                xyz_zero[3*normal_orientation+0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
     1444                xyz_zero[3*normal_orientation+1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]);
     1445                xyz_zero[3*normal_orientation+2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]);
     1446
     1447                /*New point 2*/
     1448                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]);
     1449                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]);
     1450                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]);
     1451
     1452                /*New point 3*/
     1453                xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[4][0]+s1*(xyz_list[3][0]-xyz_list[4][0]);
     1454                xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[4][1]+s1*(xyz_list[3][1]-xyz_list[4][1]);
     1455                xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[4][2]+s1*(xyz_list[3][2]-xyz_list[4][2]);
     1456
     1457                /*New point 4*/
     1458                xyz_zero[3*(2+normal_orientation)+0]=xyz_list[4][0]+s2*(xyz_list[5][0]-xyz_list[4][0]);
     1459                xyz_zero[3*(2+normal_orientation)+1]=xyz_list[4][1]+s2*(xyz_list[5][1]-xyz_list[4][1]);
     1460                xyz_zero[3*(2+normal_orientation)+2]=xyz_list[4][2]+s2*(xyz_list[5][2]-xyz_list[4][2]);
     1461        }
     1462        else if(levelset[0]==0 && levelset[1]==0){ //front is on point 0 and 1
     1463                xyz_zero[3*0+0]=xyz_list[0][0];
     1464                xyz_zero[3*0+1]=xyz_list[0][1];
     1465                xyz_zero[3*0+2]=xyz_list[0][2];
     1466
     1467                /*New point 2*/
     1468                xyz_zero[3*1+0]=xyz_list[1][0];
     1469                xyz_zero[3*1+1]=xyz_list[1][1];
     1470                xyz_zero[3*1+2]=xyz_list[1][2];
     1471
     1472                /*New point 3*/
     1473                xyz_zero[3*2+0]=xyz_list[4][0];
     1474                xyz_zero[3*2+1]=xyz_list[4][1];
     1475                xyz_zero[3*2+2]=xyz_list[4][2];
     1476
     1477                /*New point 4*/
     1478                xyz_zero[3*3+0]=xyz_list[3][0];
     1479                xyz_zero[3*3+1]=xyz_list[3][1];
     1480                xyz_zero[3*3+2]=xyz_list[3][2];
     1481        }
     1482        else if(levelset[0]==0 && levelset[2]==0){ //front is on point 0 and 1
     1483                xyz_zero[3*0+0]=xyz_list[2][0];
     1484                xyz_zero[3*0+1]=xyz_list[2][1];
     1485                xyz_zero[3*0+2]=xyz_list[2][2];
     1486
     1487                /*New point 2*/
     1488                xyz_zero[3*1+0]=xyz_list[0][0];
     1489                xyz_zero[3*1+1]=xyz_list[0][1];
     1490                xyz_zero[3*1+2]=xyz_list[0][2];
     1491
     1492                /*New point 3*/
     1493                xyz_zero[3*2+0]=xyz_list[3][0];
     1494                xyz_zero[3*2+1]=xyz_list[3][1];
     1495                xyz_zero[3*2+2]=xyz_list[3][2];
     1496
     1497                /*New point 4*/
     1498                xyz_zero[3*3+0]=xyz_list[5][0];
     1499                xyz_zero[3*3+1]=xyz_list[5][1];
     1500                xyz_zero[3*3+2]=xyz_list[5][2];
     1501        }
     1502        else if(levelset[1]==0 && levelset[2]==0){ //front is on point 0 and 1
     1503                xyz_zero[3*0+0]=xyz_list[1][0];
     1504                xyz_zero[3*0+1]=xyz_list[1][1];
     1505                xyz_zero[3*0+2]=xyz_list[1][2];
     1506
     1507                /*New point 2*/
     1508                xyz_zero[3*1+0]=xyz_list[2][0];
     1509                xyz_zero[3*1+1]=xyz_list[2][1];
     1510                xyz_zero[3*1+2]=xyz_list[2][2];
     1511
     1512                /*New point 3*/
     1513                xyz_zero[3*2+0]=xyz_list[5][0];
     1514                xyz_zero[3*2+1]=xyz_list[5][1];
     1515                xyz_zero[3*2+2]=xyz_list[5][2];
     1516
     1517                /*New point 4*/
     1518                xyz_zero[3*3+0]=xyz_list[4][0];
     1519                xyz_zero[3*3+1]=xyz_list[4][1];
     1520                xyz_zero[3*3+2]=xyz_list[4][2];
     1521        }
     1522        else _error_("Case not covered");
    13131523}
    13141524/*}}}*/
    … …  
    76967906ElementVector* Penta::CreatePVectorDiagnosticPattyn(void){
    76977907
     7908        /*compute all load vectors for this element*/
     7909        ElementVector* pe1=CreatePVectorDiagnosticPattynDrivingStress();
     7910        ElementVector* pe2=CreatePVectorDiagnosticPattynFront();
     7911        ElementVector* pe =new ElementVector(pe1,pe2);
     7912
     7913        /*clean-up and return*/
     7914        delete pe1;
     7915        delete pe2;
     7916        return pe;
     7917}
     7918/*}}}*/
     7919/*FUNCTION Penta::CreatePVectorDiagnosticPattynDrivingStress{{{*/
     7920ElementVector* Penta::CreatePVectorDiagnosticPattynDrivingStress(void){
     7921
    76987922        /*Constants*/
    76997923        const int    numdof=NDOF2*NUMVERTICES;
    … …  
    77417965}
    77427966/*}}}*/
     7967/*FUNCTION Penta::CreatePVectorDiagnosticPattynFront{{{*/
     7968ElementVector* Penta::CreatePVectorDiagnosticPattynFront(void){
     7969
     7970        /*Intermediaries */
     7971        IssmDouble  ls[6];
     7972        IssmDouble  xyz_list[NUMVERTICES][3];
     7973        bool        isfront;
     7974
     7975        /*Retrieve all inputs and parameters*/
     7976        GetInputListOnVertices(&ls[0],IcelevelsetEnum);
     7977
     7978        /*If the level set is awlays <=0, there is no ice front here*/
     7979        isfront = false;
     7980        if(ls[0]>0. || ls[1]>0. || ls[2]>0.){
     7981                if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]==0.)){
     7982                        isfront = true;
     7983                }
     7984        }
     7985
     7986        /*If no front, return NULL*/
     7987        if(!isfront) return NULL;
     7988
     7989        /*Fetch number of nodes and dof for this finite element*/
     7990        int         numnodes = this->NumberofNodes();
     7991        int         numdof   = numnodes*NDOF2;
     7992        IssmDouble  rho_ice,rho_water,gravity;
     7993        IssmDouble  Jdet,surface,z_g,water_pressure,ice_pressure,air_pressure;
     7994        IssmDouble  surface_under_water,base_under_water,pressure;
     7995        GaussPenta*  gauss;
     7996        IssmDouble* basis = xNew<IssmDouble>(numnodes);
     7997        IssmDouble  xyz_list_front[4][3];
     7998        IssmDouble  area_coordinates[4][3];
     7999        IssmDouble  normal[3];
     8000
     8001        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     8002        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
     8003        Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     8004        rho_water=matpar->GetRhoWater();
     8005        rho_ice  =matpar->GetRhoIce();
     8006        gravity  =matpar->GetG();
     8007        GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,IcelevelsetEnum);
     8008        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4);
     8009        GetQuadNormal(&normal[0],xyz_list_front);
     8010
     8011        /*Start looping on Gaussian points*/
     8012        IssmDouble zmax=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
     8013        IssmDouble zmin=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
     8014        if(zmax>0 && zmin<0) gauss=new GaussPenta(area_coordinates,3,10); //refined in vertical because of the sea level discontinuity
     8015        else                 gauss=new GaussPenta(area_coordinates,3,3);
     8016
     8017        for(int ig=gauss->begin();ig<gauss->end();ig++){
     8018
     8019                gauss->GaussPoint(ig);
     8020                surface_input->GetInputValue(&surface,gauss);
     8021                z_g=GetZcoord(gauss);
     8022                GetNodalFunctions(basis,gauss);
     8023                GetQuadJacobianDeterminant(&Jdet,xyz_list_front,gauss);
     8024
     8025                water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
     8026                ice_pressure=rho_ice*gravity*(surface-z_g);
     8027                air_pressure=0;
     8028                pressure = ice_pressure + water_pressure + air_pressure;
     8029
     8030                for (int i=0;i<numnodes;i++){
     8031                        pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
     8032                        pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
     8033                }
     8034        }
     8035
     8036        /*Transform coordinate system*/
     8037        TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum);
     8038
     8039        /*Clean up and return*/
     8040        xDelete<IssmDouble>(basis);
     8041        delete gauss;
     8042        return pe;
     8043}
     8044/*}}}*/
    77438045/*FUNCTION Penta::CreatePVectorDiagnosticStokes {{{*/
    77448046ElementVector* Penta::CreatePVectorDiagnosticStokes(void){
    … …  
    77478049        ElementVector* pe1;
    77488050        ElementVector* pe2;
     8051        ElementVector* pe3;
    77498052        ElementVector* pe;
    77508053        parameters->FindParam(&fe_stokes,FlowequationFeStokesEnum);
    … …  
    77558058                        pe1=CreatePVectorDiagnosticStokesViscous();
    77568059                        pe2=CreatePVectorDiagnosticStokesShelf();
    7757                         pe =new ElementVector(pe1,pe2);
     8060                        pe3=CreatePVectorDiagnosticStokesFront();
     8061                        pe =new ElementVector(pe1,pe2,pe3);
    77588062                        break;
    77598063                case 1:
    … …  
    77618065                        pe1=CreatePVectorDiagnosticStokesGLSViscous();
    77628066                        pe2=CreatePVectorDiagnosticStokesShelf();
    7763                         pe =new ElementVector(pe1,pe2);
     8067                        pe3=CreatePVectorDiagnosticStokesFront();
     8068                        pe =new ElementVector(pe1,pe2,pe3);
    77648069                        break;
    77658070                default:
    … …  
    77718076        delete pe1;
    77728077        delete pe2;
     8078        delete pe3;
    77738079        return pe;
    77748080}
    … …  
    78608166
    78618167        /*Clean up and return*/
     8168        delete gauss;
     8169        return pe;
     8170}
     8171/*}}}*/
     8172/*FUNCTION Penta::CreatePVectorDiagnosticStokesFront{{{*/
     8173ElementVector* Penta::CreatePVectorDiagnosticStokesFront(void){
     8174
     8175        /*Intermediaries */
     8176        IssmDouble  ls[6];
     8177        IssmDouble  xyz_list[NUMVERTICES][3];
     8178        bool        isfront;
     8179
     8180        /*Retrieve all inputs and parameters*/
     8181        GetInputListOnVertices(&ls[0],IcelevelsetEnum);
     8182
     8183        /*If the level set is awlays <=0, there is no ice front here*/
     8184        isfront = false;
     8185        if(ls[0]>0. || ls[1]>0. || ls[2]>0.){
     8186                if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]==0.)){
     8187                        isfront = true;
     8188                }
     8189        }
     8190
     8191        /*If no front, return NULL*/
     8192        if(!isfront) return NULL;
     8193
     8194        /*Fetch number of nodes and dof for this finite element*/
     8195        int         numnodes = this->NumberofNodes();
     8196        int         numdof   = numnodes*NDOF4;
     8197        IssmDouble  rho_ice,rho_water,gravity;
     8198        IssmDouble  Jdet,z_g,water_pressure,air_pressure;
     8199        IssmDouble  surface_under_water,base_under_water,pressure;
     8200        GaussPenta* gauss;
     8201        IssmDouble* basis = xNew<IssmDouble>(numnodes);
     8202        IssmDouble  xyz_list_front[4][3];
     8203        IssmDouble  area_coordinates[4][3];
     8204        IssmDouble  normal[3];
     8205
     8206        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     8207        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
     8208        rho_water=matpar->GetRhoWater();
     8209        rho_ice  =matpar->GetRhoIce();
     8210        gravity  =matpar->GetG();
     8211        GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,IcelevelsetEnum);
     8212        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4);
     8213        GetQuadNormal(&normal[0],xyz_list_front);
     8214
     8215        /*Start looping on Gaussian points*/
     8216        IssmDouble zmax=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
     8217        IssmDouble zmin=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
     8218        if(zmax>0 && zmin<0) gauss=new GaussPenta(area_coordinates,3,10); //refined in vertical because of the sea level discontinuity
     8219        else                 gauss=new GaussPenta(area_coordinates,3,3);
     8220
     8221        for(int ig=gauss->begin();ig<gauss->end();ig++){
     8222
     8223                gauss->GaussPoint(ig);
     8224                z_g=GetZcoord(gauss);
     8225                GetNodalFunctions(basis,gauss);
     8226                GetQuadJacobianDeterminant(&Jdet,xyz_list_front,gauss);
     8227
     8228                water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
     8229                air_pressure=0;
     8230                pressure = water_pressure + air_pressure;
     8231
     8232                for (int i=0;i<numnodes;i++){
     8233                        pe->values[4*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
     8234                        pe->values[4*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
     8235                        pe->values[4*i+2]+= pressure*Jdet*gauss->weight*normal[2]*basis[i];
     8236                }
     8237        }
     8238
     8239        /*Transform coordinate system*/
     8240        TransformLoadVectorCoord(pe,nodes,numnodes,XYZPEnum);
     8241
     8242
     8243        /*Clean up and return*/
     8244        xDelete<IssmDouble>(basis);
    78628245        delete gauss;
    78638246        return pe;

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

    r15497 r15538  
    182182                ElementVector* CreatePVectorPrognostic(void);
    183183                ElementVector* CreatePVectorSlope(void);
     184                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[6][3],int numpoints);
    184185                void             GetDofList(int** pdoflist,int approximation_enum,int setenum);
    185186                void             GetVertexPidList(int* doflist);
    … …  
    193194                void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    194195                void             GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
     196                void           GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
    195197                void             GetSolutionFromInputsEnthalpy(Vector<IssmDouble>* solutiong);
    196198                IssmDouble     GetStabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa);
    … …  
    198200                void    GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
    199201                Penta*  GetUpperElement(void);
     202                void    GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[6][3],int levelsetenum);
    200203                Penta*  GetLowerElement(void);
    201204                Penta*  GetBasalElement(void);
    … …  
    284287                ElementVector* CreatePVectorDiagnosticL1L2(void);
    285288                ElementVector* CreatePVectorDiagnosticPattyn(void);
     289                ElementVector* CreatePVectorDiagnosticPattynDrivingStress(void);
     290                ElementVector* CreatePVectorDiagnosticPattynFront(void);
    286291                ElementVector* CreatePVectorDiagnosticPattynStokes(void);
    287292                ElementVector* CreatePVectorDiagnosticStokes(void);
     293                ElementVector* CreatePVectorDiagnosticStokesFront(void);
    288294                ElementVector* CreatePVectorDiagnosticStokesViscous(void);
    289295                ElementVector* CreatePVectorDiagnosticStokesGLSViscous(void);

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

    r15535 r15538  
    31653165        }
    31663166
    3167         return NULL;
    31683167        /*If no front, return NULL*/
    31693168        if(!isfront) return NULL;
    … …  
    32003199                thickness_input->GetInputValue(&thickness,gauss);
    32013200                bed_input->GetInputValue(&bed,gauss);
     3201                GetSegmentJacobianDeterminant(&Jdet,&xyz_list_front[0][0],gauss);
     3202                GetNodalFunctions(basis,gauss);
    32023203
    32033204                surface_under_water=min(0.,thickness+bed); // 0 if the top of the glacier is above water level
    … …  
    32063207                ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
    32073208                air_pressure=0;
    3208 
    32093209                pressure = ice_pressure + water_pressure + air_pressure;
    3210 
    3211                 GetSegmentJacobianDeterminant(&Jdet,&xyz_list_front[0][0],gauss);
    3212                 GetNodalFunctions(basis,gauss);
    32133210
    32143211                for (int i=0;i<numnodes;i++){

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp

    r15104 r15538  
    237237}
    238238/*}}}*/
     239/*FUNCTION GaussPenta::GaussPenta(IssmDouble area_coordinates[4][3],int order_horiz,int order_vert){{{*/
     240GaussPenta::GaussPenta(IssmDouble area_coordinates[4][3],int order_horiz,int order_vert){
     241
     242        /*Intermediaties*/
     243        IssmPDouble *seg_horiz_coords  = NULL;
     244        IssmPDouble *seg_horiz_weights = NULL;
     245        IssmPDouble *seg_vert_coords   = NULL;
     246        IssmPDouble *seg_vert_weights  = NULL;
     247
     248        /*get the gauss points using the product of two line rules*/
     249        GaussLegendreLinear(&seg_horiz_coords,&seg_horiz_weights,order_horiz);
     250        GaussLegendreLinear(&seg_vert_coords, &seg_vert_weights, order_vert);
     251
     252        /*Allocate GaussPenta fields*/
     253        numgauss=order_horiz*order_vert;
     254        coords1=xNew<IssmDouble>(numgauss);
     255        coords2=xNew<IssmDouble>(numgauss);
     256        coords3=xNew<IssmDouble>(numgauss);
     257        coords4=xNew<IssmDouble>(numgauss);
     258        weights=xNew<IssmDouble>(numgauss);
     259
     260        /*Quads: get the gauss points using the product of two line rules  */
     261        for(int i=0;i<order_horiz;i++){
     262                for(int j=0;j<order_vert;j++){
     263                        coords1[i*order_vert+j]=0.5*(area_coordinates[0][0]+area_coordinates[1][0]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][0]-area_coordinates[0][0]);
     264                        coords2[i*order_vert+j]=0.5*(area_coordinates[0][1]+area_coordinates[1][1]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][1]-area_coordinates[0][1]);
     265                        coords3[i*order_vert+j]=0.5*(area_coordinates[0][2]+area_coordinates[1][2]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][2]-area_coordinates[0][2]);
     266                        coords4[i*order_vert+j]=seg_vert_coords[j];
     267                        weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
     268                }
     269        }
     270
     271        /*clean-up*/
     272        xDelete<IssmPDouble>(seg_horiz_coords);
     273        xDelete<IssmPDouble>(seg_horiz_weights);
     274        xDelete<IssmPDouble>(seg_vert_coords);
     275        xDelete<IssmPDouble>(seg_vert_weights);
     276}
     277/*}}}*/
    239278/*FUNCTION GaussPenta::~GaussPenta(){{{*/
    240279GaussPenta::~GaussPenta(){

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h

    r15071 r15538  
    3535                GaussPenta(int index1, int index2, int index3, int order);
    3636                GaussPenta(int index1, int index2, int index3, int index4,int order_horiz,int order_vert);
     37                GaussPenta(IssmDouble area_coordinates[4][3],int order_horiz,int order_vert);
    3738                ~GaussPenta();
    3839

  • issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp

    r15517 r15538  
    134134
    135135        /*clean up*/
    136         xDelete<double>(seg_coords);
    137         xDelete<double>(seg_weights);
     136        xDelete<IssmPDouble>(seg_coords);
     137        xDelete<IssmPDouble>(seg_weights);
    138138}
    139139/*}}}*/

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp

    r15525 r15538  
    6666                /*Do not create ice front if Hutter or Stokes elements*/
    6767                if (reCast<int,IssmDouble>(*(elements_type+element))==HutterApproximationEnum) continue;
     68                if (reCast<int,IssmDouble>(*(elements_type+element))==MacAyealApproximationEnum) continue;
     69                if (reCast<int,IssmDouble>(*(elements_type+element))==L1L2ApproximationEnum) continue;
     70                if (reCast<int,IssmDouble>(*(elements_type+element))==PattynApproximationEnum) continue;
     71                if (reCast<int,IssmDouble>(*(elements_type+element))==StokesApproximationEnum) continue;
    6872
    6973                /*Create and  add load: */
    70                 if (reCast<int,IssmDouble>(*(elements_type+element))==(MacAyealApproximationEnum) && iomodel->dim==2){
    71                         loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal2dIceFrontEnum,DiagnosticHorizAnalysisEnum));
    72                         count++;
    73                 }
    74                 else if (reCast<int,IssmDouble>(*(elements_type+element))==(MacAyealApproximationEnum) && iomodel->dim==3){
    75                         loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum));
    76                         count++;
    77                 }
    78                 else if (reCast<int,IssmDouble>(*(elements_type+element))==(L1L2ApproximationEnum)){
    79                         loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum));
    80                         count++;
    81                 }
    82                 else if (reCast<int,IssmDouble>(*(elements_type+element))==(PattynApproximationEnum)){
    83                         loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,PattynIceFrontEnum,DiagnosticHorizAnalysisEnum));
    84                         count++;
    85                 }
    86                 else if (reCast<int,IssmDouble>(*(elements_type+element))==(L1L2ApproximationEnum)){
    87                         loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,L1L2IceFrontEnum,DiagnosticHorizAnalysisEnum));
    88                         count++;
    89                 }
    90                 else if (reCast<int,IssmDouble>(*(elements_type+element))==(StokesApproximationEnum)){
    91                         loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,StokesIceFrontEnum,DiagnosticHorizAnalysisEnum));
    92                         count++;
    93                 }
    94                 else if (reCast<int,IssmDouble>(*(elements_type+element))==(MacAyealPattynApproximationEnum)){
     74                if (reCast<int,IssmDouble>(*(elements_type+element))==(MacAyealPattynApproximationEnum)){
    9575                        loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum));
    9676                        count++;
Note: See TracChangeset for help on using the changeset viewer.