Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15516) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15517) @@ -6949,5 +6949,5 @@ /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(3,2); + gauss=new GaussPenta(5,5); for(int ig=gauss->begin();igend();ig++){ @@ -7927,5 +7927,5 @@ /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(3,2); + gauss=new GaussPenta(5,5); for(int ig=gauss->begin();igend();ig++){ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15516) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15517) @@ -692,4 +692,40 @@ } /*}}}*/ +/*FUNCTION Tria::GetAreaCoordinates{{{*/ +void Tria::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints){ + /*Computeportion of the element that is grounded*/ + + int i,j,k; + IssmDouble area_init,area_portion; + IssmDouble xyz_bis[3][3]; + + GetJacobianDeterminant(&area_init, &xyz_list[0][0],NULL); + + /*Initialize xyz_list with original xyz_list of triangle coordinates*/ + for(j=0;j<3;j++){ + for(k=0;k<3;j++){ + xyz_bis[j][k]=xyz_list[j][k]; + } + } + for(i=0;i0 && gl[1]>0 && gl[2]>0){ // All grounded + point=0; + f1=1.; + f2=1.; + } + else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating + point=0; + f1=0.; + f2=0.; + } + else{ + if(gl[0]*gl[1]*gl[2]<0) floating=false; + + if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 + point=2; + f1=gl[2]/(gl[2]-gl[0]); + f2=gl[2]/(gl[2]-gl[1]); + } + else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 + point=0; + f1=gl[0]/(gl[0]-gl[1]); + f2=gl[0]/(gl[0]-gl[2]); + } + else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 + point=1; + f1=gl[1]/(gl[1]-gl[2]); + f2=gl[1]/(gl[1]-gl[0]); + } + } + *point1=point; + *fraction1=f1; + *fraction2=f2; + *mainlyfloating=floating; +} +/*}}}*/ /*FUNCTION Tria::GetGroundedPortion{{{*/ IssmDouble Tria::GetGroundedPortion(IssmDouble* xyz_list){ @@ -826,56 +916,80 @@ } /*}}}*/ -/*FUNCTION Tria::GetGroundedPart{{{*/ -void Tria::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating){ +/*FUNCTION Tria::GetSegmentNormal {{{*/ +void Tria:: GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[2][3]){ + + /*Build unit outward pointing vector*/ + IssmDouble vector[2]; + IssmDouble norm; + + vector[0]=xyz_list[1][0] - xyz_list[0][0]; + vector[1]=xyz_list[1][1] - xyz_list[0][1]; + + norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0)); + + normal[0]= + vector[1]/norm; + normal[1]= - vector[0]/norm; +} +/*}}}*/ +/*FUNCTION Tria::GetZeroLevelsetCoordinates{{{*/ +void Tria::GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum){ /*Computeportion of the element that is grounded*/ - bool floating=true; - int point; - const IssmPDouble epsilon= 1.e-15; - IssmDouble gl[3]; - IssmDouble f1,f2; + int normal_orientation; + IssmDouble s1,s2; + IssmDouble levelset[3]; /*Recover parameters and values*/ - GetInputListOnVertices(&gl[0],GLlevelsetEnum); - - /*Be sure that values are not zero*/ - if(gl[0]==0) gl[0]=gl[0]+epsilon; - if(gl[1]==0) gl[1]=gl[1]+epsilon; - if(gl[2]==0) gl[2]=gl[2]+epsilon; - - /*Check that not all nodes are grounded or floating*/ - if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded - point=0; - f1=1.; - f2=1.; - } - else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating - point=0; - f1=0.; - f2=0.; - } - else{ - if(gl[0]*gl[1]*gl[2]<0) floating=false; - - if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 - point=2; - f1=gl[2]/(gl[2]-gl[0]); - f2=gl[2]/(gl[2]-gl[1]); - } - else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 - point=0; - f1=gl[0]/(gl[0]-gl[1]); - f2=gl[0]/(gl[0]-gl[2]); - } - else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 - point=1; - f1=gl[1]/(gl[1]-gl[2]); - f2=gl[1]/(gl[1]-gl[0]); - } - } - *point1=point; - *fraction1=f1; - *fraction2=f2; - *mainlyfloating=floating; + GetInputListOnVertices(&levelset[0],levelsetenum); + + if(levelset[0]*levelset[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 + /*Portion of the segments*/ + s1=levelset[2]/(levelset[2]-levelset[1]); + s2=levelset[2]/(levelset[2]-levelset[0]); + + if(levelset[2]>0) normal_orientation=0; + /*New point 1*/ + *(xyz_zero+3*normal_orientation+0)=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]); + *(xyz_zero+3*normal_orientation+1)=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]); + *(xyz_zero+3*normal_orientation+2)=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]); + + /*New point 0*/ + *(xyz_zero+3*(1-normal_orientation)+0)=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]); + *(xyz_zero+3*(1-normal_orientation)+1)=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]); + *(xyz_zero+3*(1-normal_orientation)+2)=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]); + } + 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 + /*Portion of the segments*/ + s1=levelset[0]/(levelset[0]-levelset[2]); + s2=levelset[0]/(levelset[0]-levelset[1]); + + if(levelset[0]>0) normal_orientation=0; + /*New point 1*/ + *(xyz_zero+3*normal_orientation+0)=xyz_list[0][0]+s1*(xyz_list[2][0]-xyz_list[0][0]); + *(xyz_zero+3*normal_orientation+1)=xyz_list[0][1]+s1*(xyz_list[2][1]-xyz_list[0][1]); + *(xyz_zero+3*normal_orientation+2)=xyz_list[0][2]+s1*(xyz_list[2][2]-xyz_list[0][2]); + + /*New point 2*/ + *(xyz_zero+3*(1-normal_orientation)+0)=xyz_list[0][0]+s2*(xyz_list[1][0]-xyz_list[0][0]); + *(xyz_zero+3*(1-normal_orientation)+1)=xyz_list[0][1]+s2*(xyz_list[1][1]-xyz_list[0][1]); + *(xyz_zero+3*(1-normal_orientation)+2)=xyz_list[0][2]+s2*(xyz_list[1][2]-xyz_list[0][2]); + } + 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 + /*Portion of the segments*/ + s1=levelset[1]/(levelset[1]-levelset[0]); + s2=levelset[1]/(levelset[1]-levelset[2]); + + if(levelset[1]>0) normal_orientation=0; + /*New point 0*/ + *(xyz_zero+3*normal_orientation+0)=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]); + *(xyz_zero+3*normal_orientation+1)=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]); + *(xyz_zero+3*normal_orientation+2)=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]); + + /*New point 2*/ + *(xyz_zero+3*(1-normal_orientation)+0)=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]); + *(xyz_zero+3*(1-normal_orientation)+1)=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]); + *(xyz_zero+3*(1-normal_orientation)+2)=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]); + } + } /*}}}*/ @@ -3006,25 +3120,70 @@ /*Intermediaries */ - int i,j; - IssmDouble ls[3]; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - int numdof = numnodes*NDOF2; - Icefront *icefront=NULL; + IssmDouble ls[3]; + IssmDouble xyz_list[NUMVERTICES][3]; + bool isfront; return NULL; + /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); GetInputListOnVertices(&ls[0],IcelevelsetEnum); - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum); - GaussTria* gauss = new GaussTria(2); - IssmDouble* basis = xNew(numnodes); - - /*Create Ice Front if necessary*/ + + /*If the level set is awlays <0, there is no ice front here*/ + isfront = false; if(ls[0]>0. || ls[1]>0. || ls[2]>0.){ if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]==0.)){ - //icefront=new Icefront("2d",inputs,matpar,MacAyealApproximationEnum,analysis_type); + isfront = true; + } + } + + /*If no front, return NULL*/ + if(!isfront) return NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF2; + IssmDouble rho_ice,rho_water,gravity; + IssmDouble Jdet,thickness,bed,water_pressure,ice_pressure,air_pressure; + IssmDouble surface_under_water,base_under_water,pressure; + GaussTria* gauss; + IssmDouble* basis = xNew(numnodes); + IssmDouble xyz_list_front[2][3]; + IssmDouble area_coordinates[2][3]; + IssmDouble normal[2]; + + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum); + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* bed_input =inputs->GetInput(BedEnum); _assert_(bed_input); + rho_water=matpar->GetRhoWater(); + rho_ice =matpar->GetRhoIce(); + gravity =matpar->GetG(); + GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,IcelevelsetEnum); + GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2); + GetSegmentNormal(&normal[0],xyz_list_front); + + /*Start looping on Gaussian points*/ + gauss=new GaussTria(area_coordinates,3); + + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + thickness_input->GetInputValue(&thickness,gauss); + bed_input->GetInputValue(&bed,gauss); + + surface_under_water=min(0.,thickness+bed); // 0 if the top of the glacier is above water level + base_under_water=min(0.,bed); // 0 if the bottom of the glacier is above water level + water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2)); + ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); + air_pressure=0; + + pressure = ice_pressure + water_pressure + air_pressure; + + GetSegmentJacobianDeterminant(&Jdet,&xyz_list_front[0][0],gauss); + GetNodalFunctions(basis,gauss); + + for (int i=0;ivalues[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; + pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; } } @@ -3037,4 +3196,5 @@ delete gauss; return pe; + } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15516) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15517) @@ -195,4 +195,5 @@ ElementVector* CreatePVectorSlope(void); IssmDouble GetArea(void); + void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints); int GetElementType(void); void GetDofList(int** pdoflist,int approximation_enum,int setenum); @@ -202,4 +203,6 @@ IssmDouble GetGroundedPortion(IssmDouble* xyz_list); void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating); + void GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[2][3]); + void GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue); Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 15516) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 15517) @@ -87,4 +87,41 @@ coords2[i]=0.5*(b1+b2) + 0.5*seg_coords[i]*(b2-b1); coords3[i]=0.5*(c1+c2) + 0.5*seg_coords[i]*(c2-c1); + weights[i]=seg_weights[i]; + } + + /*Initialize static fields as undefined*/ + weight=UNDEF; + coord1=UNDEF; + coord2=UNDEF; + coord3=UNDEF; + + /*clean up*/ + xDelete(seg_coords); + xDelete(seg_weights); +} +/*}}}*/ +/*FUNCTION GaussTria::GaussTria(IssmDouble area_coordinates,int order) {{{*/ +GaussTria::GaussTria(IssmDouble area_coordinates[2][3],int order){ + + /*Intermediaties*/ + IssmPDouble *seg_coords = NULL; + IssmPDouble *seg_weights = NULL; + int i,index3; + + /*Get Segment gauss points*/ + numgauss=order; + GaussLegendreLinear(&seg_coords,&seg_weights,numgauss); + + /*Allocate GaussTria fields*/ + coords1=xNew(numgauss); + coords2=xNew(numgauss); + coords3=xNew(numgauss); + weights=xNew(numgauss); + + /*Build Triangle Gauss point*/ + for(i=0;i