Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15537) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15538) @@ -784,4 +784,41 @@ } /*}}}*/ +/*FUNCTION Penta::GetAreaCoordinates{{{*/ +void Penta::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble xyz_zero[4][3],IssmDouble xyz_list[6][3],int numpoints){ + /*Computeportion of the element that is grounded*/ + + int i,j,k; + IssmDouble area_init,area_portion; + 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.; + + + /*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]; + } + } + for(i=0;i0){ //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=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]); + + /*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]); + + /*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]); + + /*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]); + } + 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=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]); + + /*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]); + + /*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]); + + /*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]); + } + 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=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]); + + /*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]); + + /*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]); + + /*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]); + } + 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]; + + /*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]; + + /*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]; + + /*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]; + } + 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]; + + /*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]; + + /*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]; + + /*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]; + } + 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]; + + /*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]; + + /*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]; + + /*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]; + } + else _error_("Case not covered"); } /*}}}*/ @@ -7696,4 +7906,18 @@ ElementVector* Penta::CreatePVectorDiagnosticPattyn(void){ + /*compute all load vectors for this element*/ + ElementVector* pe1=CreatePVectorDiagnosticPattynDrivingStress(); + ElementVector* pe2=CreatePVectorDiagnosticPattynFront(); + ElementVector* pe =new ElementVector(pe1,pe2); + + /*clean-up and return*/ + delete pe1; + delete pe2; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticPattynDrivingStress{{{*/ +ElementVector* Penta::CreatePVectorDiagnosticPattynDrivingStress(void){ + /*Constants*/ const int numdof=NDOF2*NUMVERTICES; @@ -7741,4 +7965,82 @@ } /*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticPattynFront{{{*/ +ElementVector* Penta::CreatePVectorDiagnosticPattynFront(void){ + + /*Intermediaries */ + IssmDouble ls[6]; + IssmDouble xyz_list[NUMVERTICES][3]; + bool isfront; + + /*Retrieve all inputs and parameters*/ + GetInputListOnVertices(&ls[0],IcelevelsetEnum); + + /*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.)){ + 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,surface,z_g,water_pressure,ice_pressure,air_pressure; + IssmDouble surface_under_water,base_under_water,pressure; + GaussPenta* gauss; + IssmDouble* basis = xNew(numnodes); + IssmDouble xyz_list_front[4][3]; + IssmDouble area_coordinates[4][3]; + IssmDouble normal[3]; + + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum); + Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_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,4); + GetQuadNormal(&normal[0],xyz_list_front); + + /*Start looping on Gaussian points*/ + IssmDouble zmax=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2]; + IssmDouble zmin=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]0 && zmin<0) gauss=new GaussPenta(area_coordinates,3,10); //refined in vertical because of the sea level discontinuity + else gauss=new GaussPenta(area_coordinates,3,3); + + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + surface_input->GetInputValue(&surface,gauss); + z_g=GetZcoord(gauss); + GetNodalFunctions(basis,gauss); + GetQuadJacobianDeterminant(&Jdet,xyz_list_front,gauss); + + water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level + ice_pressure=rho_ice*gravity*(surface-z_g); + air_pressure=0; + pressure = ice_pressure + water_pressure + air_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*/ + TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum); + + /*Clean up and return*/ + xDelete(basis); + delete gauss; + return pe; +} +/*}}}*/ /*FUNCTION Penta::CreatePVectorDiagnosticStokes {{{*/ ElementVector* Penta::CreatePVectorDiagnosticStokes(void){ @@ -7747,4 +8049,5 @@ ElementVector* pe1; ElementVector* pe2; + ElementVector* pe3; ElementVector* pe; parameters->FindParam(&fe_stokes,FlowequationFeStokesEnum); @@ -7755,5 +8058,6 @@ pe1=CreatePVectorDiagnosticStokesViscous(); pe2=CreatePVectorDiagnosticStokesShelf(); - pe =new ElementVector(pe1,pe2); + pe3=CreatePVectorDiagnosticStokesFront(); + pe =new ElementVector(pe1,pe2,pe3); break; case 1: @@ -7761,5 +8065,6 @@ pe1=CreatePVectorDiagnosticStokesGLSViscous(); pe2=CreatePVectorDiagnosticStokesShelf(); - pe =new ElementVector(pe1,pe2); + pe3=CreatePVectorDiagnosticStokesFront(); + pe =new ElementVector(pe1,pe2,pe3); break; default: @@ -7771,4 +8076,5 @@ delete pe1; delete pe2; + delete pe3; return pe; } @@ -7860,4 +8166,81 @@ /*Clean up and return*/ + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticStokesFront{{{*/ +ElementVector* Penta::CreatePVectorDiagnosticStokesFront(void){ + + /*Intermediaries */ + IssmDouble ls[6]; + IssmDouble xyz_list[NUMVERTICES][3]; + bool isfront; + + /*Retrieve all inputs and parameters*/ + GetInputListOnVertices(&ls[0],IcelevelsetEnum); + + /*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.)){ + 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*NDOF4; + IssmDouble rho_ice,rho_water,gravity; + IssmDouble Jdet,z_g,water_pressure,air_pressure; + IssmDouble surface_under_water,base_under_water,pressure; + GaussPenta* gauss; + IssmDouble* basis = xNew(numnodes); + IssmDouble xyz_list_front[4][3]; + IssmDouble area_coordinates[4][3]; + IssmDouble normal[3]; + + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum); + 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,4); + GetQuadNormal(&normal[0],xyz_list_front); + + /*Start looping on Gaussian points*/ + IssmDouble zmax=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2]; + IssmDouble zmin=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]0 && zmin<0) gauss=new GaussPenta(area_coordinates,3,10); //refined in vertical because of the sea level discontinuity + else gauss=new GaussPenta(area_coordinates,3,3); + + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + z_g=GetZcoord(gauss); + GetNodalFunctions(basis,gauss); + GetQuadJacobianDeterminant(&Jdet,xyz_list_front,gauss); + + water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level + air_pressure=0; + pressure = water_pressure + air_pressure; + + for (int i=0;ivalues[4*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; + pe->values[4*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; + pe->values[4*i+2]+= pressure*Jdet*gauss->weight*normal[2]*basis[i]; + } + } + + /*Transform coordinate system*/ + TransformLoadVectorCoord(pe,nodes,numnodes,XYZPEnum); + + + /*Clean up and return*/ + xDelete(basis); delete gauss; return pe; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15537) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15538) @@ -182,4 +182,5 @@ ElementVector* CreatePVectorPrognostic(void); ElementVector* CreatePVectorSlope(void); + void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[6][3],int numpoints); void GetDofList(int** pdoflist,int approximation_enum,int setenum); void GetVertexPidList(int* doflist); @@ -193,4 +194,5 @@ void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype); void GetPhi(IssmDouble* phi, IssmDouble* epsilon, IssmDouble viscosity); + void GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]); void GetSolutionFromInputsEnthalpy(Vector* solutiong); IssmDouble GetStabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa); @@ -198,4 +200,5 @@ void GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* 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); Penta* GetBasalElement(void); @@ -284,6 +287,9 @@ ElementVector* CreatePVectorDiagnosticL1L2(void); ElementVector* CreatePVectorDiagnosticPattyn(void); + ElementVector* CreatePVectorDiagnosticPattynDrivingStress(void); + ElementVector* CreatePVectorDiagnosticPattynFront(void); ElementVector* CreatePVectorDiagnosticPattynStokes(void); ElementVector* CreatePVectorDiagnosticStokes(void); + ElementVector* CreatePVectorDiagnosticStokesFront(void); ElementVector* CreatePVectorDiagnosticStokesViscous(void); ElementVector* CreatePVectorDiagnosticStokesGLSViscous(void); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15537) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15538) @@ -3165,5 +3165,4 @@ } - return NULL; /*If no front, return NULL*/ if(!isfront) return NULL; @@ -3200,4 +3199,6 @@ thickness_input->GetInputValue(&thickness,gauss); bed_input->GetInputValue(&bed,gauss); + GetSegmentJacobianDeterminant(&Jdet,&xyz_list_front[0][0],gauss); + GetNodalFunctions(basis,gauss); surface_under_water=min(0.,thickness+bed); // 0 if the top of the glacier is above water level @@ -3206,9 +3207,5 @@ 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;i(numgauss); + coords2=xNew(numgauss); + coords3=xNew(numgauss); + coords4=xNew(numgauss); + weights=xNew(numgauss); + + /*Quads: get the gauss points using the product of two line rules */ + for(int i=0;i(seg_horiz_coords); + xDelete(seg_horiz_weights); + xDelete(seg_vert_coords); + xDelete(seg_vert_weights); +} +/*}}}*/ /*FUNCTION GaussPenta::~GaussPenta(){{{*/ GaussPenta::~GaussPenta(){ Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h (revision 15537) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h (revision 15538) @@ -35,4 +35,5 @@ GaussPenta(int index1, int index2, int index3, int order); GaussPenta(int index1, int index2, int index3, int index4,int order_horiz,int order_vert); + GaussPenta(IssmDouble area_coordinates[4][3],int order_horiz,int order_vert); ~GaussPenta(); Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 15537) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 15538) @@ -134,6 +134,6 @@ /*clean up*/ - xDelete(seg_coords); - xDelete(seg_weights); + xDelete(seg_coords); + xDelete(seg_weights); } /*}}}*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 15537) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 15538) @@ -66,31 +66,11 @@ /*Do not create ice front if Hutter or Stokes elements*/ if (reCast(*(elements_type+element))==HutterApproximationEnum) continue; + if (reCast(*(elements_type+element))==MacAyealApproximationEnum) continue; + if (reCast(*(elements_type+element))==L1L2ApproximationEnum) continue; + if (reCast(*(elements_type+element))==PattynApproximationEnum) continue; + if (reCast(*(elements_type+element))==StokesApproximationEnum) continue; /*Create and add load: */ - if (reCast(*(elements_type+element))==(MacAyealApproximationEnum) && iomodel->dim==2){ - loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal2dIceFrontEnum,DiagnosticHorizAnalysisEnum)); - count++; - } - else if (reCast(*(elements_type+element))==(MacAyealApproximationEnum) && iomodel->dim==3){ - loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum)); - count++; - } - else if (reCast(*(elements_type+element))==(L1L2ApproximationEnum)){ - loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum)); - count++; - } - else if (reCast(*(elements_type+element))==(PattynApproximationEnum)){ - loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,PattynIceFrontEnum,DiagnosticHorizAnalysisEnum)); - count++; - } - else if (reCast(*(elements_type+element))==(L1L2ApproximationEnum)){ - loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,L1L2IceFrontEnum,DiagnosticHorizAnalysisEnum)); - count++; - } - else if (reCast(*(elements_type+element))==(StokesApproximationEnum)){ - loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,StokesIceFrontEnum,DiagnosticHorizAnalysisEnum)); - count++; - } - else if (reCast(*(elements_type+element))==(MacAyealPattynApproximationEnum)){ + if (reCast(*(elements_type+element))==(MacAyealPattynApproximationEnum)){ loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum)); count++;