Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16837) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16838) @@ -1110,4 +1110,63 @@ ElementVector* StressbalanceAnalysis::CreatePVectorSSAFront(Element* element){/*{{{*/ + /*If no front, return NULL*/ + if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL; + + /*Intermediaries*/ + IssmDouble rho_ice,rho_water,gravity; + IssmDouble Jdet,thickness,bed,water_pressure,ice_pressure; + IssmDouble surface_under_water,base_under_water,pressure; + IssmDouble *xyz_list = NULL; + IssmDouble *xyz_list_front = NULL; + IssmDouble normal[2]; + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementVector* pe = element->NewElementVector(SSAApproximationEnum); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* bed_input =element->GetInput(BedEnum); _assert_(bed_input); + rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum); + rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + gravity = element->GetMaterialParameter(ConstantsGEnum); + element->GetVerticesCoordinates(&xyz_list); + element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); + element->NormalBase(&normal[0],xyz_list_front); + + /*Start looping on Gaussian points*/ + Gauss* gauss=element->NewGauss(xyz_list,xyz_list_front,3); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + thickness_input->GetInputValue(&thickness,gauss); + bed_input->GetInputValue(&bed,gauss); + element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss); + element->NodalFunctions(basis,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*(surface_under_water*surface_under_water - base_under_water*base_under_water); + ice_pressure = 1.0/2.0*gravity*rho_ice*thickness*thickness; + pressure = ice_pressure + water_pressure; + + 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]; + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,XYEnum); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(xyz_list_front); + xDelete(basis); + delete gauss; + return pe; return NULL; }/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16837) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16838) @@ -82,4 +82,5 @@ virtual void GetDofListPressure(int** pdoflist,int setenum)=0; virtual void JacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss)=0; + virtual void JacobianDeterminantSurface(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss)=0; virtual void JacobianDeterminantBase(IssmDouble* Jdet,IssmDouble* xyz_list_base,Gauss* gauss)=0; virtual void JacobianDeterminantTop(IssmDouble* Jdet,IssmDouble* xyz_list_base,Gauss* gauss)=0; @@ -130,4 +131,5 @@ virtual Gauss* NewGauss(void)=0; virtual Gauss* NewGauss(int order)=0; + virtual Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order)=0; virtual Gauss* NewGaussBase(int order)=0; virtual Gauss* NewGaussTop(int order)=0; @@ -152,4 +154,5 @@ virtual int PressureInterpolation()=0; virtual bool IsZeroLevelset(int levelset_enum)=0; + virtual void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum)=0; #ifdef _HAVE_RESPONSES_ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16837) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16838) @@ -884,5 +884,5 @@ /*}}}*/ /*FUNCTION Penta::GetAreaCoordinates{{{*/ -void Penta::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble xyz_zero[4][3],IssmDouble xyz_list[6][3],int numpoints){ +void Penta::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){ /*Computeportion of the element that is grounded*/ @@ -891,10 +891,10 @@ IssmDouble xyz_bis[3][3]; - 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.; + area_init=fabs(xyz_list[1*3+0]*xyz_list[2*3+1] - xyz_list[1*3+1]*xyz_list[2*3+0] + xyz_list[0*3+0]*xyz_list[1*3+1] - xyz_list[0*3+1]*xyz_list[1*3+0] + xyz_list[2*3+0]*xyz_list[0*3+1] - xyz_list[2*3+1]*xyz_list[0*3+0])/2.; /*Initialize xyz_list with original xyz_list of triangle coordinates*/ for(j=0;j<3;j++){ for(k=0;k<3;k++){ - xyz_bis[j][k]=xyz_list[j][k]; + xyz_bis[j][k]=xyz_list[j*3+k]; } } @@ -903,5 +903,5 @@ for(k=0;k<3;k++){ /*Change appropriate line*/ - xyz_bis[j][k]=xyz_zero[i][k]; + xyz_bis[j][k]=xyz_zero[i*3+k]; } @@ -913,5 +913,5 @@ for(k=0;k<3;k++){ /*Reinitialize xyz_list with original coordinates*/ - xyz_bis[j][k]=xyz_list[j][k]; + xyz_bis[j][k]=xyz_list[j*3+k]; } } @@ -1512,5 +1512,5 @@ /*}}}*/ /*FUNCTION Penta::GetQuadNormal {{{*/ -void Penta:: GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){ +void Penta:: GetQuadNormal(IssmDouble* normal,IssmDouble* xyz_list){ /*Build unit outward pointing vector*/ @@ -1519,10 +1519,10 @@ IssmDouble norm; - AB[0]=xyz_list[1][0] - xyz_list[0][0]; - AB[1]=xyz_list[1][1] - xyz_list[0][1]; - AB[2]=xyz_list[1][2] - xyz_list[0][2]; - AC[0]=xyz_list[2][0] - xyz_list[0][0]; - AC[1]=xyz_list[2][1] - xyz_list[0][1]; - AC[2]=xyz_list[2][2] - xyz_list[0][2]; + AB[0]=xyz_list[1*3+0] - xyz_list[0*3+0]; + AB[1]=xyz_list[1*3+1] - xyz_list[0*3+1]; + AB[2]=xyz_list[1*3+2] - xyz_list[0*3+2]; + AC[0]=xyz_list[2*3+0] - xyz_list[0*3+0]; + AC[1]=xyz_list[2*3+1] - xyz_list[0*3+1]; + AC[2]=xyz_list[2*3+2] - xyz_list[0*3+2]; cross(normal,AB,AC); @@ -1648,6 +1648,6 @@ } /*}}}*/ -/*FUNCTION Penta::GetZeroLevelsetCoordinates{{{*/ -void Penta::GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[6][3],int levelsetenum){ +/*FUNCTION Penta::ZeroLevelsetCoordinates{{{*/ +void Penta::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){ /*Computeportion of the element that is grounded*/ @@ -1655,4 +1655,5 @@ IssmDouble s1,s2; IssmDouble levelset[NUMVERTICES]; + IssmDouble* xyz_zero = xNew(4*3); /*Recover parameters and values*/ @@ -1666,22 +1667,22 @@ if(levelset[2]>0.) normal_orientation=1; /*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]); + xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]); + xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]); + xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+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]); + xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]); + xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]); + xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]); /*New point 3*/ - xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[5][0]+s1*(xyz_list[4][0]-xyz_list[5][0]); - xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[5][1]+s1*(xyz_list[4][1]-xyz_list[5][1]); - xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[5][2]+s1*(xyz_list[4][2]-xyz_list[5][2]); + xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[5*3+0]+s1*(xyz_list[4*3+0]-xyz_list[5*3+0]); + xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[5*3+1]+s1*(xyz_list[4*3+1]-xyz_list[5*3+1]); + xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[5*3+2]+s1*(xyz_list[4*3+2]-xyz_list[5*3+2]); /*New point 4*/ - xyz_zero[3*(2+normal_orientation)+0]=xyz_list[5][0]+s2*(xyz_list[3][0]-xyz_list[5][0]); - xyz_zero[3*(2+normal_orientation)+1]=xyz_list[5][1]+s2*(xyz_list[3][1]-xyz_list[5][1]); - xyz_zero[3*(2+normal_orientation)+2]=xyz_list[5][2]+s2*(xyz_list[3][2]-xyz_list[5][2]); + xyz_zero[3*(2+normal_orientation)+0]=xyz_list[5*3+0]+s2*(xyz_list[3*3+0]-xyz_list[5*3+0]); + xyz_zero[3*(2+normal_orientation)+1]=xyz_list[5*3+1]+s2*(xyz_list[3*3+1]-xyz_list[5*3+1]); + xyz_zero[3*(2+normal_orientation)+2]=xyz_list[5*3+2]+s2*(xyz_list[3*3+2]-xyz_list[5*3+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 @@ -1692,22 +1693,22 @@ if(levelset[0]>0.) normal_orientation=1; /*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]); + xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]); + xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]); + xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+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]); + xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]); + xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]); + xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]); /*New point 3*/ - xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[3][0]+s1*(xyz_list[5][0]-xyz_list[3][0]); - xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[3][1]+s1*(xyz_list[5][1]-xyz_list[3][1]); - xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[3][2]+s1*(xyz_list[5][2]-xyz_list[3][2]); + xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[3*3+0]+s1*(xyz_list[5*3+0]-xyz_list[3*3+0]); + xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[3*3+1]+s1*(xyz_list[5*3+1]-xyz_list[3*3+1]); + xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[3*3+2]+s1*(xyz_list[5*3+2]-xyz_list[3*3+2]); /*New point 4*/ - xyz_zero[3*(2+normal_orientation)+0]=xyz_list[3][0]+s2*(xyz_list[4][0]-xyz_list[3][0]); - xyz_zero[3*(2+normal_orientation)+1]=xyz_list[3][1]+s2*(xyz_list[4][1]-xyz_list[3][1]); - xyz_zero[3*(2+normal_orientation)+2]=xyz_list[3][2]+s2*(xyz_list[4][2]-xyz_list[3][2]); + xyz_zero[3*(2+normal_orientation)+0]=xyz_list[3*3+0]+s2*(xyz_list[4*3+0]-xyz_list[3*3+0]); + xyz_zero[3*(2+normal_orientation)+1]=xyz_list[3*3+1]+s2*(xyz_list[4*3+1]-xyz_list[3*3+1]); + xyz_zero[3*(2+normal_orientation)+2]=xyz_list[3*3+2]+s2*(xyz_list[4*3+2]-xyz_list[3*3+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 @@ -1718,84 +1719,87 @@ if(levelset[1]>0.) normal_orientation=1; /*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]); + xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]); + xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]); + xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+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]); + xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]); + xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]); + xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]); /*New point 3*/ - xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[4][0]+s1*(xyz_list[3][0]-xyz_list[4][0]); - xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[4][1]+s1*(xyz_list[3][1]-xyz_list[4][1]); - xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[4][2]+s1*(xyz_list[3][2]-xyz_list[4][2]); + xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[4*3+0]+s1*(xyz_list[3*3+0]-xyz_list[4*3+0]); + xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[4*3+1]+s1*(xyz_list[3*3+1]-xyz_list[4*3+1]); + xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[4*3+2]+s1*(xyz_list[3*3+2]-xyz_list[4*3+2]); /*New point 4*/ - xyz_zero[3*(2+normal_orientation)+0]=xyz_list[4][0]+s2*(xyz_list[5][0]-xyz_list[4][0]); - xyz_zero[3*(2+normal_orientation)+1]=xyz_list[4][1]+s2*(xyz_list[5][1]-xyz_list[4][1]); - xyz_zero[3*(2+normal_orientation)+2]=xyz_list[4][2]+s2*(xyz_list[5][2]-xyz_list[4][2]); + xyz_zero[3*(2+normal_orientation)+0]=xyz_list[4*3+0]+s2*(xyz_list[5*3+0]-xyz_list[4*3+0]); + xyz_zero[3*(2+normal_orientation)+1]=xyz_list[4*3+1]+s2*(xyz_list[5*3+1]-xyz_list[4*3+1]); + xyz_zero[3*(2+normal_orientation)+2]=xyz_list[4*3+2]+s2*(xyz_list[5*3+2]-xyz_list[4*3+2]); } else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1 - xyz_zero[3*0+0]=xyz_list[0][0]; - xyz_zero[3*0+1]=xyz_list[0][1]; - xyz_zero[3*0+2]=xyz_list[0][2]; + xyz_zero[3*0+0]=xyz_list[0*3+0]; + xyz_zero[3*0+1]=xyz_list[0*3+1]; + xyz_zero[3*0+2]=xyz_list[0*3+2]; /*New point 2*/ - xyz_zero[3*1+0]=xyz_list[1][0]; - xyz_zero[3*1+1]=xyz_list[1][1]; - xyz_zero[3*1+2]=xyz_list[1][2]; + xyz_zero[3*1+0]=xyz_list[1*3+0]; + xyz_zero[3*1+1]=xyz_list[1*3+1]; + xyz_zero[3*1+2]=xyz_list[1*3+2]; /*New point 3*/ - xyz_zero[3*2+0]=xyz_list[4][0]; - xyz_zero[3*2+1]=xyz_list[4][1]; - xyz_zero[3*2+2]=xyz_list[4][2]; + xyz_zero[3*2+0]=xyz_list[4*3+0]; + xyz_zero[3*2+1]=xyz_list[4*3+1]; + xyz_zero[3*2+2]=xyz_list[4*3+2]; /*New point 4*/ - xyz_zero[3*3+0]=xyz_list[3][0]; - xyz_zero[3*3+1]=xyz_list[3][1]; - xyz_zero[3*3+2]=xyz_list[3][2]; + xyz_zero[3*3+0]=xyz_list[3*3+0]; + xyz_zero[3*3+1]=xyz_list[3*3+1]; + xyz_zero[3*3+2]=xyz_list[3*3+2]; } else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1 - xyz_zero[3*0+0]=xyz_list[2][0]; - xyz_zero[3*0+1]=xyz_list[2][1]; - xyz_zero[3*0+2]=xyz_list[2][2]; + xyz_zero[3*0+0]=xyz_list[2*3+0]; + xyz_zero[3*0+1]=xyz_list[2*3+1]; + xyz_zero[3*0+2]=xyz_list[2*3+2]; /*New point 2*/ - xyz_zero[3*1+0]=xyz_list[0][0]; - xyz_zero[3*1+1]=xyz_list[0][1]; - xyz_zero[3*1+2]=xyz_list[0][2]; + xyz_zero[3*1+0]=xyz_list[0*3+0]; + xyz_zero[3*1+1]=xyz_list[0*3+1]; + xyz_zero[3*1+2]=xyz_list[0*3+2]; /*New point 3*/ - xyz_zero[3*2+0]=xyz_list[3][0]; - xyz_zero[3*2+1]=xyz_list[3][1]; - xyz_zero[3*2+2]=xyz_list[3][2]; + xyz_zero[3*2+0]=xyz_list[3*3+0]; + xyz_zero[3*2+1]=xyz_list[3*3+1]; + xyz_zero[3*2+2]=xyz_list[3*3+2]; /*New point 4*/ - xyz_zero[3*3+0]=xyz_list[5][0]; - xyz_zero[3*3+1]=xyz_list[5][1]; - xyz_zero[3*3+2]=xyz_list[5][2]; + xyz_zero[3*3+0]=xyz_list[5*3+0]; + xyz_zero[3*3+1]=xyz_list[5*3+1]; + xyz_zero[3*3+2]=xyz_list[5*3+2]; } else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1 - xyz_zero[3*0+0]=xyz_list[1][0]; - xyz_zero[3*0+1]=xyz_list[1][1]; - xyz_zero[3*0+2]=xyz_list[1][2]; + xyz_zero[3*0+0]=xyz_list[1*3+0]; + xyz_zero[3*0+1]=xyz_list[1*3+1]; + xyz_zero[3*0+2]=xyz_list[1*3+2]; /*New point 2*/ - xyz_zero[3*1+0]=xyz_list[2][0]; - xyz_zero[3*1+1]=xyz_list[2][1]; - xyz_zero[3*1+2]=xyz_list[2][2]; + xyz_zero[3*1+0]=xyz_list[2*3+0]; + xyz_zero[3*1+1]=xyz_list[2*3+1]; + xyz_zero[3*1+2]=xyz_list[2*3+2]; /*New point 3*/ - xyz_zero[3*2+0]=xyz_list[5][0]; - xyz_zero[3*2+1]=xyz_list[5][1]; - xyz_zero[3*2+2]=xyz_list[5][2]; + xyz_zero[3*2+0]=xyz_list[5*3+0]; + xyz_zero[3*2+1]=xyz_list[5*3+1]; + xyz_zero[3*2+2]=xyz_list[5*3+2]; /*New point 4*/ - xyz_zero[3*3+0]=xyz_list[4][0]; - xyz_zero[3*3+1]=xyz_list[4][1]; - xyz_zero[3*3+2]=xyz_list[4][2]; + xyz_zero[3*3+0]=xyz_list[4*3+0]; + xyz_zero[3*3+1]=xyz_list[4*3+1]; + xyz_zero[3*3+2]=xyz_list[4*3+2]; } else _error_("Case not covered"); + + /*Assign output pointer*/ + *pxyz_zero= xyz_zero; } /*}}}*/ @@ -8764,5 +8768,5 @@ IssmDouble surface_under_water,base_under_water,pressure; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_front[4][3]; + IssmDouble* xyz_list_front = NULL; IssmDouble area_coordinates[4][3]; IssmDouble normal[3]; @@ -8782,6 +8786,6 @@ rho_ice =matpar->GetRhoIce(); gravity =matpar->GetG(); - GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum); - GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4); + ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum); + GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],4); GetQuadNormal(&normal[0],xyz_list_front); @@ -8854,5 +8858,5 @@ IssmDouble surface_under_water,base_under_water,pressure; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_front[4][3]; + IssmDouble* xyz_list_front = NULL; IssmDouble area_coordinates[4][3]; IssmDouble normal[3]; @@ -8869,14 +8873,14 @@ /*Initialize Element matrix and vectors*/ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); ElementVector* pe = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum); IssmDouble* vbasis = xNew(vnumnodes); /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); rho_water=matpar->GetRhoWater(); rho_ice =matpar->GetRhoIce(); gravity =matpar->GetG(); - GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum); - GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4); + ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum); + GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],4); GetQuadNormal(&normal[0],xyz_list_front); @@ -8910,4 +8914,5 @@ /*Clean up and return*/ xDelete(cs_list); + xDelete(xyz_list_front); xDelete(vbasis); delete gauss; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16837) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16838) @@ -114,4 +114,5 @@ int PressureInterpolation(); bool IsZeroLevelset(int levelset_enum); + void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); void ResultInterpolation(int* pinterpolation,int output_enum); @@ -211,5 +212,5 @@ IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure); IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure); - void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[6][3],int numpoints); + void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints); void GetVertexPidList(int* doflist); @@ -231,9 +232,8 @@ void GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype); void GetPhi(IssmDouble* phi, IssmDouble* epsilon, IssmDouble viscosity); - void GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]); + void GetQuadNormal(IssmDouble* normal,IssmDouble* xyz_list); void GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input); void GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input, Input* vz_input); Penta* GetUpperElement(void); - void GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[6][3],int levelsetenum); Penta* GetLowerElement(void); void InputChangeName(int input_enum, int enum_type_old); @@ -247,4 +247,5 @@ bool IsNodeOnShelfFromFlags(IssmDouble* flags); void JacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); + void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); @@ -252,4 +253,5 @@ Gauss* NewGauss(void); Gauss* NewGauss(int order); + Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");}; Gauss* NewGaussBase(int order); Gauss* NewGaussTop(int order); Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 16837) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 16838) @@ -1945,21 +1945,21 @@ /*}}}*/ /*FUNCTION PentaRef::GetQuadJacobianDeterminant{{{*/ -void PentaRef::GetQuadJacobianDeterminant(IssmDouble* Jdet,IssmDouble xyz_list[4][3],Gauss* gauss){ +void PentaRef::GetQuadJacobianDeterminant(IssmDouble* Jdet,IssmDouble* xyz_list,Gauss* gauss){ /*This routine returns the values of the nodal functions at the gaussian point.*/ IssmDouble x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4; - x1=xyz_list[0][0]; - y1=xyz_list[0][1]; - z1=xyz_list[0][2]; - x2=xyz_list[1][0]; - y2=xyz_list[1][1]; - z2=xyz_list[1][2]; - x3=xyz_list[2][0]; - y3=xyz_list[2][1]; - z3=xyz_list[2][2]; - x4=xyz_list[3][0]; - y4=xyz_list[3][1]; - z4=xyz_list[3][2]; + x1=xyz_list[0*3+0]; + y1=xyz_list[0*3+1]; + z1=xyz_list[0*3+2]; + x2=xyz_list[1*3+0]; + y2=xyz_list[1*3+1]; + z2=xyz_list[1*3+2]; + x3=xyz_list[2*3+0]; + y3=xyz_list[2*3+1]; + z3=xyz_list[2*3+2]; + x4=xyz_list[3*3+0]; + y4=xyz_list[3*3+1]; + z4=xyz_list[3*3+2]; /*Jdet = (Area of the trapezoid)/(Area trapezoid ref) with AreaRef = 4*/ Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h (revision 16837) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h (revision 16838) @@ -35,5 +35,5 @@ void GetNodalFunctionsP1DerivativesReference(IssmDouble* dl1dl6,Gauss* gauss); void GetNodalFunctionsMINIDerivativesReference(IssmDouble* dl1dl7,Gauss* gauss); - void GetQuadJacobianDeterminant(IssmDouble* Jdet, IssmDouble xyz_list[4][3],Gauss* gauss); + void GetQuadJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss); void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16837) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16838) @@ -114,4 +114,5 @@ bool IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");}; void JacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; + void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; @@ -142,4 +143,5 @@ Gauss* NewGauss(void){_error_("not implemented yet");}; Gauss* NewGauss(int order){_error_("not implemented yet");}; + Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");}; Gauss* NewGaussBase(int order){_error_("not implemented yet");}; Gauss* NewGaussTop(int order){_error_("not implemented yet");}; @@ -151,4 +153,5 @@ void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; bool IsZeroLevelset(int levelset_enum){_error_("not implemented");}; + void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");}; #ifdef _HAVE_THERMAL_ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16837) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16838) @@ -879,5 +879,5 @@ /*}}}*/ /*FUNCTION Tria::GetAreaCoordinates{{{*/ -void Tria::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints){ +void Tria::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){ /*Computeportion of the element that is grounded*/ @@ -891,5 +891,5 @@ for(j=0;j<3;j++){ for(k=0;k<3;k++){ - xyz_bis[j][k]=xyz_list[j][k]; + xyz_bis[j][k]=xyz_list[j*3+k]; } } @@ -898,5 +898,5 @@ for(k=0;k<3;k++){ /*Change appropriate line*/ - xyz_bis[j][k]=xyz_zero[i][k]; + xyz_bis[j][k]=xyz_zero[i*3+k]; } @@ -908,5 +908,5 @@ for(k=0;k<3;k++){ /*Reinitialize xyz_list with original coordinates*/ - xyz_bis[j][k]=xyz_list[j][k]; + xyz_bis[j][k]=xyz_list[j*3+k]; } } @@ -1184,7 +1184,6 @@ } /*}}}*/ -/*FUNCTION Tria::GetZeroLevelsetCoordinates{{{*/ -void Tria::GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum){ - /*Computeportion of the element that is grounded*/ +/*FUNCTION Tria::ZeroLevelsetCoordinates{{{*/ +void Tria::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){ int normal_orientation=0; @@ -1193,4 +1192,5 @@ /*Recover parameters and values*/ + IssmDouble* xyz_zero = xNew(2*3); GetInputListOnVertices(&levelset[0],levelsetenum); @@ -1202,12 +1202,12 @@ if(levelset[2]>0.) normal_orientation=1; /*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]); + xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]); + xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]); + xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+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]); + xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]); + xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]); + xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+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 @@ -1218,12 +1218,12 @@ if(levelset[0]>0.) normal_orientation=1; /*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]); + xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]); + xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]); + xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+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]); + xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]); + xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]); + xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+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 @@ -1234,44 +1234,47 @@ if(levelset[1]>0.) normal_orientation=1; /*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]); + xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]); + xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]); + xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+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]); + xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]); + xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]); + xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]); } else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1 - xyz_zero[3*0+0]=xyz_list[0][0]; - xyz_zero[3*0+1]=xyz_list[0][1]; - xyz_zero[3*0+2]=xyz_list[0][2]; + xyz_zero[3*0+0]=xyz_list[0*3+0]; + xyz_zero[3*0+1]=xyz_list[0*3+1]; + xyz_zero[3*0+2]=xyz_list[0*3+2]; /*New point 2*/ - xyz_zero[3*1+0]=xyz_list[1][0]; - xyz_zero[3*1+1]=xyz_list[1][1]; - xyz_zero[3*1+2]=xyz_list[1][2]; + xyz_zero[3*1+0]=xyz_list[1*3+0]; + xyz_zero[3*1+1]=xyz_list[1*3+1]; + xyz_zero[3*1+2]=xyz_list[1*3+2]; } else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1 - xyz_zero[3*0+0]=xyz_list[2][0]; - xyz_zero[3*0+1]=xyz_list[2][1]; - xyz_zero[3*0+2]=xyz_list[2][2]; + xyz_zero[3*0+0]=xyz_list[2*3+0]; + xyz_zero[3*0+1]=xyz_list[2*3+1]; + xyz_zero[3*0+2]=xyz_list[2*3+2]; /*New point 2*/ - xyz_zero[3*1+0]=xyz_list[0][0]; - xyz_zero[3*1+1]=xyz_list[0][1]; - xyz_zero[3*1+2]=xyz_list[0][2]; + xyz_zero[3*1+0]=xyz_list[0*3+0]; + xyz_zero[3*1+1]=xyz_list[0*3+1]; + xyz_zero[3*1+2]=xyz_list[0*3+2]; } else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1 - xyz_zero[3*0+0]=xyz_list[1][0]; - xyz_zero[3*0+1]=xyz_list[1][1]; - xyz_zero[3*0+2]=xyz_list[1][2]; + xyz_zero[3*0+0]=xyz_list[1*3+0]; + xyz_zero[3*0+1]=xyz_list[1*3+1]; + xyz_zero[3*0+2]=xyz_list[1*3+2]; /*New point 2*/ - xyz_zero[3*1+0]=xyz_list[2][0]; - xyz_zero[3*1+1]=xyz_list[2][1]; - xyz_zero[3*1+2]=xyz_list[2][2]; + xyz_zero[3*1+0]=xyz_list[2*3+0]; + xyz_zero[3*1+1]=xyz_list[2*3+1]; + xyz_zero[3*1+2]=xyz_list[2*3+2]; } else _error_("Case not covered"); + + /*Assign output pointer*/ + *pxyz_zero= xyz_zero; } /*}}}*/ @@ -2014,4 +2017,12 @@ } /*}}}*/ +/*FUNCTION Tria::JacobianDeterminantSurface{{{*/ +void Tria::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){ + + _assert_(gauss->Enum()==GaussTriaEnum); + this->GetSegmentJacobianDeterminant(pJdet,xyz_list,(GaussTria*)gauss); + +} +/*}}}*/ /*FUNCTION Tria::HasEdgeOnBed {{{*/ bool Tria::HasEdgeOnBed(){ @@ -2207,4 +2218,14 @@ Gauss* Tria::NewGauss(int order){ return new GaussTria(order); +} +/*}}}*/ +/*FUNCTION Tria::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){{{*/ +Gauss* Tria::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){ + + IssmDouble area_coordinates[2][3]; + + GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2); + + return new GaussTria(area_coordinates,order); } /*}}}*/ @@ -3930,5 +3951,5 @@ IssmDouble Jdet,water_pressure,air_pressure,pressure; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_front[2][3]; + IssmDouble* xyz_list_front = NULL; IssmDouble area_coordinates[2][3]; IssmDouble normal[2]; @@ -3953,11 +3974,11 @@ rho_ice = matpar->GetRhoIce(); gravity = matpar->GetG(); - GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum); - GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2); - NormalBase(&normal[0],&xyz_list_front[0][0]); + ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum); + GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],2); + NormalBase(&normal[0],xyz_list_front); /*Start looping on Gaussian points*/ - IssmDouble ymax=max(xyz_list_front[0][1],xyz_list_front[1][1]); - IssmDouble ymin=min(xyz_list_front[0][1],xyz_list_front[1][1]); + IssmDouble ymax=max(xyz_list_front[0*3+1],xyz_list_front[1*3+1]); + IssmDouble ymin=min(xyz_list_front[0*3+1],xyz_list_front[1*3+1]); if(ymax>0. && ymin<0.) gauss=new GaussTria(area_coordinates,30); //refined in vertical because of the sea level discontinuity else gauss=new GaussTria(area_coordinates,3); @@ -3968,5 +3989,5 @@ y_g=GetYcoord(gauss); GetNodalFunctionsVelocity(vbasis,gauss); - GetSegmentJacobianDeterminant(&Jdet,&xyz_list_front[0][0],gauss); + GetSegmentJacobianDeterminant(&Jdet,xyz_list_front,gauss); water_pressure = rho_water*gravity*min(0.,y_g); //0 if the gaussian point is above water level @@ -3980,5 +4001,4 @@ } - /*Transform coordinate system*/ ::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list); @@ -3986,4 +4006,5 @@ /*Clean up and return*/ xDelete(cs_list); + xDelete(xyz_list_front); xDelete(vbasis); delete gauss; @@ -4201,5 +4222,5 @@ IssmDouble surface_under_water,base_under_water,pressure; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_front[2][3]; + IssmDouble *xyz_list_front = NULL; IssmDouble area_coordinates[2][3]; IssmDouble normal[2]; @@ -4219,7 +4240,7 @@ rho_ice = matpar->GetRhoIce(); gravity = matpar->GetG(); - GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum); - GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2); - NormalBase(&normal[0],&xyz_list_front[0][0]); + ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum); + GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],2); + NormalBase(&normal[0],xyz_list_front); /*Start looping on Gaussian points*/ @@ -4230,5 +4251,5 @@ thickness_input->GetInputValue(&thickness,gauss); bed_input->GetInputValue(&bed,gauss); - GetSegmentJacobianDeterminant(&Jdet,&xyz_list_front[0][0],gauss); + GetSegmentJacobianDeterminant(&Jdet,xyz_list_front,gauss); GetNodalFunctions(basis,gauss); @@ -4251,4 +4272,5 @@ /*Clean up and return*/ xDelete(basis); + xDelete(xyz_list_front); delete gauss; return pe; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16837) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16838) @@ -135,4 +135,5 @@ void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement); IssmDouble TimeAdapt(); + void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); bool IsZeroLevelset(int levelset_enum); @@ -248,5 +249,5 @@ IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");}; IssmDouble GetArea(void); - void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints); + void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints); int GetElementType(void); @@ -259,5 +260,4 @@ void NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; IssmDouble GetMaterialParameter(int enum_in); - void GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum); Input* GetInput(int inputenum); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype); @@ -277,4 +277,5 @@ bool IsInput(int name); void JacobianDeterminant(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); + void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; @@ -282,4 +283,5 @@ Gauss* NewGauss(void); Gauss* NewGauss(int order); + Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order); Gauss* NewGaussBase(int order){_error_("not implemented yet");}; Gauss* NewGaussTop(int order){_error_("not implemented yet");};