Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 25933) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 25934) @@ -299,22 +299,22 @@ /* Start looping on the number of vertices: */ - GaussPenta* gauss=new GaussPenta(); + GaussPenta gauss; for(int iv=0;iv<3;iv++){ - gauss->GaussVertex(iv); + gauss.GaussVertex(iv); /*Get velocity components and thickness*/ - B_input->GetInputValue(&B,gauss); - n_input->GetInputValue(&n,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - gr_input->GetInputValue(&groundedice,gauss); - bs_input->GetInputValue(&bed,gauss); - smax_fl_input->GetInputValue(&sigma_max_floating,gauss); - smax_gr_input->GetInputValue(&sigma_max_grounded,gauss); + B_input->GetInputValue(&B,&gauss); + n_input->GetInputValue(&n,&gauss); + vx_input->GetInputValue(&vx,&gauss); + vy_input->GetInputValue(&vy,&gauss); + gr_input->GetInputValue(&groundedice,&gauss); + bs_input->GetInputValue(&bed,&gauss); + smax_fl_input->GetInputValue(&sigma_max_floating,&gauss); + smax_gr_input->GetInputValue(&sigma_max_grounded,&gauss); vel=sqrt(vx*vx+vy*vy)+1.e-14; - sl_input->GetInputValue(&sealevel,gauss); + sl_input->GetInputValue(&sealevel,&gauss); /*Compute strain rate and viscosity: */ - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); /*Get Eigen values*/ @@ -359,7 +359,4 @@ this->InputExtrude(CalvingCalvingrateEnum,-1); this->InputExtrude(SigmaVMEnum,-1); - - /*Clean up and return*/ - delete gauss; } /*}}}*/ @@ -382,15 +379,15 @@ /* Start looping on the number of vertices: */ - GaussPenta* gauss=new GaussPenta(); - for (int iv=0;ivGaussVertex(iv); + GaussPenta gauss; + for(int iv=0;ivGetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); + vx_input->GetInputValue(&vx,&gauss); + vy_input->GetInputValue(&vy,&gauss); vel=vx*vx+vy*vy; - strainparallel_input->GetInputValue(&strainparallel,gauss); - strainperpendicular_input->GetInputValue(&strainperpendicular,gauss); - levermanncoeff_input->GetInputValue(&propcoeff,gauss); + strainparallel_input->GetInputValue(&strainparallel,&gauss); + strainperpendicular_input->GetInputValue(&strainperpendicular,&gauss); + levermanncoeff_input->GetInputValue(&propcoeff,&gauss); /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */ @@ -407,7 +404,4 @@ this->AddBasalInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum); this->AddBasalInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); - - /*Clean up and return*/ - delete gauss; }/*}}}*/ void Penta::CalvingFluxLevelset(){/*{{{*/ @@ -761,5 +755,4 @@ IssmDouble tau_yz[NUMVERTICES]; IssmDouble tau_eff[NUMVERTICES]; - GaussPenta* gauss=NULL; /* Get node coordinates and dof list: */ @@ -772,11 +765,11 @@ /* Start looping on the number of vertices: */ - gauss=new GaussPenta(); + GaussPenta gauss; for (int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /*Compute strain rate viscosity and pressure: */ - this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); + this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); /*Compute Stress*/ @@ -802,7 +795,4 @@ this->AddInput(DeviatoricStresszzEnum,&tau_zz[0],P1DGEnum); this->AddInput(DeviatoricStresseffectiveEnum,&tau_eff[0],P1DGEnum); - - /*Clean up and return*/ - delete gauss; } /*}}}*/ @@ -818,5 +808,4 @@ IssmDouble sigma_xz[NUMVERTICES]; IssmDouble sigma_yz[NUMVERTICES]; - GaussPenta* gauss=NULL; /* Get node coordinates and dof list: */ @@ -830,12 +819,12 @@ /* Start looping on the number of vertices: */ - gauss=new GaussPenta(); - for (int iv=0;ivGaussVertex(iv); + GaussPenta gauss; + for(int iv=0;ivStrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - pressure_input->GetInputValue(&pressure,gauss); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); + this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); + pressure_input->GetInputValue(&pressure,&gauss); /*Compute Stress*/ @@ -855,7 +844,4 @@ this->AddInput(StressTensoryzEnum,&sigma_yz[0],P1DGEnum); this->AddInput(StressTensorzzEnum,&sigma_zz[0],P1DGEnum); - - /*Clean up and return*/ - delete gauss; } /*}}}*/ @@ -952,15 +938,14 @@ control_gradient->Serve(NUMVERTICES,&lidlist[0]); - GaussPenta* gauss=new GaussPenta(); + GaussPenta gauss; for (int iv=0;ivGaussVertex(iv); - - control_value->GetInputValue(&value,gauss); - control_gradient->GetInputValue(&gradient,gauss); + gauss.GaussVertex(iv); + + control_value->GetInputValue(&value,&gauss); + control_gradient->GetInputValue(&gradient,&gauss); values[iv] = reCast(value)/reCast(connectivity[iv]); gradients[iv] = reCast(gradient)/reCast(connectivity[iv]); } - delete gauss; vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL); @@ -1139,11 +1124,11 @@ /*1. Recover stresses at the base*/ - GaussPenta* gauss=new GaussPenta(); + GaussPenta gauss; for (int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /*Compute strain rate viscosity and pressure: */ - this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); + this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); /*FIXME: this is for Hongju only*/ //pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]); @@ -1180,7 +1165,4 @@ } } - - /*clean up*/ - delete gauss; } /*}}}*/ @@ -1631,9 +1613,7 @@ int index = this->GetNodeIndex(node); - GaussPenta* gauss=new GaussPenta(); - gauss->GaussNode(this->element_type,index); - - input->GetInputValue(pvalue,gauss); - delete gauss; + GaussPenta gauss; + gauss.GaussNode(this->element_type,index); + input->GetInputValue(pvalue,&gauss); } /*}}}*/ @@ -1645,9 +1625,7 @@ int index = this->GetVertexIndex(vertex); - GaussPenta* gauss=new GaussPenta(); - gauss->GaussVertex(index); - - input->GetInputValue(pvalue,gauss); - delete gauss; + GaussPenta gauss; + gauss.GaussVertex(index); + input->GetInputValue(pvalue,&gauss); } /*}}}*/ @@ -2851,12 +2829,12 @@ /* Start looping on the number of vertices: */ - GaussPenta* gauss=new GaussPenta(); + GaussPenta gauss; for(int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /* Advection */ - vx_input->GetInputValue(&v[0],gauss); - vy_input->GetInputValue(&v[1],gauss); - gr_input->GetInputValue(&groundedice,gauss); + vx_input->GetInputValue(&v[0],&gauss); + vy_input->GetInputValue(&v[1],&gauss); + gr_input->GetInputValue(&groundedice,&gauss); /*Get calving speed*/ @@ -2864,8 +2842,8 @@ case DefaultCalvingEnum: case CalvingVonmisesEnum: - lsf_slopex_input->GetInputValue(&dlsf[0],gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); - calvingrate_input->GetInputValue(&calvingrate,gauss); - meltingrate_input->GetInputValue(&meltingrate,gauss); + lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); + calvingrate_input->GetInputValue(&calvingrate,&gauss); + meltingrate_input->GetInputValue(&meltingrate,&gauss); if(groundedice<0) meltingrate = 0.; @@ -2886,7 +2864,7 @@ case CalvingLevermannEnum: - calvingratex_input->GetInputValue(&c[0],gauss); - if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss); - meltingrate_input->GetInputValue(&meltingrate,gauss); + calvingratex_input->GetInputValue(&c[0],&gauss); + if(dim==2) calvingratey_input->GetInputValue(&c[1],&gauss); + meltingrate_input->GetInputValue(&meltingrate,&gauss); norm_calving=0.; for(i=0;iGetInputValue(&dlsf[0],gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); - meltingrate_input->GetInputValue(&meltingrate,gauss); + lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); + meltingrate_input->GetInputValue(&meltingrate,&gauss); norm_dlsf=0.; @@ -2917,7 +2895,7 @@ case CalvingHabEnum: - lsf_slopex_input->GetInputValue(&dlsf[0],gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); - meltingrate_input->GetInputValue(&meltingrate,gauss); + lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); + meltingrate_input->GetInputValue(&meltingrate,&gauss); norm_dlsf=0.; @@ -2938,7 +2916,7 @@ case CalvingCrevasseDepthEnum: - lsf_slopex_input->GetInputValue(&dlsf[0],gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); - meltingrate_input->GetInputValue(&meltingrate,gauss); + lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); + meltingrate_input->GetInputValue(&meltingrate,&gauss); if(groundedice<0) meltingrate = 0.; @@ -2962,9 +2940,9 @@ case CalvingDev2Enum: { - lsf_slopex_input->GetInputValue(&dlsf[0],gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); - calvingrate_input->GetInputValue(&calvingrate,gauss); - meltingrate_input->GetInputValue(&meltingrate,gauss); - gr_input->GetInputValue(&groundedice,gauss); + lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); + calvingrate_input->GetInputValue(&calvingrate,&gauss); + meltingrate_input->GetInputValue(&meltingrate,&gauss); + gr_input->GetInputValue(&groundedice,&gauss); //idea: no retreat on ice above critical calving height "calvinghaf" . Limit using regularized Heaviside function. @@ -3021,6 +2999,4 @@ this->InputExtrude(MovingFrontalVxEnum,-1); this->InputExtrude(MovingFrontalVyEnum,-1); - /*Clean up and return*/ - delete gauss; } /*}}}*/ @@ -3142,5 +3118,5 @@ int found=0; IssmDouble value; - GaussPenta* gauss=NULL; + GaussPenta gauss; /*First, serarch the input: */ @@ -3154,6 +3130,6 @@ if(data){ /*ok, we are good. retrieve value of input at vertex :*/ - gauss=new GaussPenta(); gauss->GaussVertex(i); - data->GetInputValue(&value,gauss); + gauss.GaussVertex(i); + data->GetInputValue(&value,&gauss); found=1; break; @@ -3162,5 +3138,4 @@ } - delete gauss; if(found)*pvalue=value; return found; @@ -3369,11 +3344,11 @@ /*Loop over basal nodes and update their CS*/ - GaussPenta* gauss = new GaussPenta(); + GaussPenta gauss; for(int i=0;iGaussNode(this->VelocityInterpolation(),indices[i]); - - slopex_input->GetInputValue(&slopex,gauss); - slopey_input->GetInputValue(&slopey,gauss); - groundedicelevelset_input->GetInputValue(&groundedice,gauss); + gauss.GaussNode(this->VelocityInterpolation(),indices[i]); + + slopex_input->GetInputValue(&slopex,&gauss); + slopey_input->GetInputValue(&slopey,&gauss); + groundedicelevelset_input->GetInputValue(&groundedice,&gauss); /*New X axis New Z axis*/ @@ -3406,5 +3381,4 @@ /*cleanup*/ xDelete(indices); - delete gauss; } /*}}}*/ @@ -3455,12 +3429,12 @@ /* Start looping on the number of vertices: */ - GaussPenta* gauss=new GaussPenta(); + GaussPenta gauss; for(int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /* Get variables */ - bed_input->GetInputValue(&bed,gauss); - qsg_input->GetInputValue(&qsg,gauss); - TF_input->GetInputValue(&TF,gauss); + bed_input->GetInputValue(&bed,&gauss); + qsg_input->GetInputValue(&qsg,&gauss); + TF_input->GetInputValue(&TF,&gauss); if(basinid[iv]==0 || basin_icefront_area[reCast(basinid[iv])-1]==0.) meltrates[iv]=0.; @@ -3484,5 +3458,4 @@ /*Cleanup and return*/ xDelete(basin_icefront_area); - delete gauss; } /*}}}*/ @@ -3692,5 +3665,4 @@ IssmDouble epsilon[6]; - GaussPenta* gauss=NULL; IssmDouble vx,vy,vel; IssmDouble strainxx; @@ -3709,15 +3681,15 @@ /* Start looping on the number of vertices: */ - gauss=new GaussPenta(); + GaussPenta gauss; for (int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /* Get the value we need*/ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); + vx_input->GetInputValue(&vx,&gauss); + vy_input->GetInputValue(&vy,&gauss); vel=vx*vx+vy*vy; /*Compute strain rate and viscosity: */ - this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); strainxx=epsilon[0]; strainyy=epsilon[1]; @@ -3730,7 +3702,4 @@ /*Add input*/ this->AddInput(StrainRateparallelEnum,&strainparallel[0],P1DGEnum); - - /*Clean up and return*/ - delete gauss; } /*}}}*/ @@ -3738,5 +3707,4 @@ IssmDouble epsilon[6]; - GaussPenta* gauss=NULL; IssmDouble vx,vy,vel; IssmDouble strainxx; @@ -3755,15 +3723,15 @@ /* Start looping on the number of vertices: */ - gauss=new GaussPenta(); + GaussPenta gauss; for (int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /* Get the value we need*/ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); + vx_input->GetInputValue(&vx,&gauss); + vy_input->GetInputValue(&vy,&gauss); vel=vx*vx+vy*vy; /*Compute strain rate and viscosity: */ - this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); strainxx=epsilon[0]; strainyy=epsilon[1]; @@ -3776,7 +3744,4 @@ /*Add input*/ this->AddInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1DGEnum); - - /*Clean up and return*/ - delete gauss; } /*}}}*/ @@ -4779,17 +4744,14 @@ Input* onbase = this->GetInput(MeshVertexonbaseEnum); _assert_(onbase); - GaussPenta* gauss=new GaussPenta(); + GaussPenta gauss; for(int iv=0;ivNumberofNodes(this->element_type);iv++){ - gauss->GaussNode(this->element_type,iv); - onbase->GetInputValue(&isonbase,gauss); + gauss.GaussNode(this->element_type,iv); + onbase->GetInputValue(&isonbase,&gauss); if(isonbase==1.){ - input->GetInputValue(&value,gauss); + input->GetInputValue(&value,&gauss); this->nodes[iv]->ApplyConstraint(0,value); } } - delete gauss; - -} -/*}}}*/ +}/*}}}*/ void Penta::UpdateConstraintsExtrudeFromTop(void){/*{{{*/ @@ -4803,11 +4765,10 @@ Input* input = this->GetInput(extrusioninput); _assert_(extrusioninput); - GaussPenta* gauss=new GaussPenta(); + GaussPenta gauss; for(int i=0;i<3;i++){ - gauss->GaussNode(P1Enum,indices[i]); - input->GetInputValue(&value,gauss); + gauss.GaussNode(P1Enum,indices[i]); + input->GetInputValue(&value,&gauss); this->nodes[indices[i]]->ApplyConstraint(0,value); } - delete gauss; } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25933) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25934) @@ -339,22 +339,22 @@ /* Start looping on the number of vertices: */ - GaussTria* gauss=new GaussTria(); + GaussTria gauss; for(int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /*Get velocity components and thickness*/ - B_input->GetInputValue(&B,gauss); - n_input->GetInputValue(&n,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - gr_input->GetInputValue(&groundedice,gauss); - bs_input->GetInputValue(&bed,gauss); - smax_fl_input->GetInputValue(&sigma_max_floating,gauss); - smax_gr_input->GetInputValue(&sigma_max_grounded,gauss); + B_input->GetInputValue(&B,&gauss); + n_input->GetInputValue(&n,&gauss); + vx_input->GetInputValue(&vx,&gauss); + vy_input->GetInputValue(&vy,&gauss); + gr_input->GetInputValue(&groundedice,&gauss); + bs_input->GetInputValue(&bed,&gauss); + smax_fl_input->GetInputValue(&sigma_max_floating,&gauss); + smax_gr_input->GetInputValue(&sigma_max_grounded,&gauss); vel=sqrt(vx*vx+vy*vy)+1.e-14; - sl_input->GetInputValue(&sealevel,gauss); + sl_input->GetInputValue(&sealevel,&gauss); /*Compute strain rate and viscosity: */ - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); /*Get Eigen values*/ @@ -399,7 +399,4 @@ this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); this->AddInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum); - - /*Clean up and return*/ - delete gauss; } /*}}}*/ @@ -437,21 +434,21 @@ /*Loop over all elements of this partition*/ - GaussTria* gauss=new GaussTria(); + GaussTria gauss; for (int iv=0;ivGaussVertex(iv); - - H_input->GetInputValue(&thickness,gauss); - bed_input->GetInputValue(&bed,gauss); - surface_input->GetInputValue(&surface,gauss); - strainrateparallel_input->GetInputValue(&strainparallel,gauss); - strainrateeffective_input->GetInputValue(&straineffective,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - waterheight_input->GetInputValue(&water_height,gauss); - s_xx_input->GetInputValue(&s_xx,gauss); - s_xy_input->GetInputValue(&s_xy,gauss); - s_yy_input->GetInputValue(&s_yy,gauss); - B_input->GetInputValue(&B,gauss); - n_input->GetInputValue(&n,gauss); + gauss.GaussVertex(iv); + + H_input->GetInputValue(&thickness,&gauss); + bed_input->GetInputValue(&bed,&gauss); + surface_input->GetInputValue(&surface,&gauss); + strainrateparallel_input->GetInputValue(&strainparallel,&gauss); + strainrateeffective_input->GetInputValue(&straineffective,&gauss); + vx_input->GetInputValue(&vx,&gauss); + vy_input->GetInputValue(&vy,&gauss); + waterheight_input->GetInputValue(&water_height,&gauss); + s_xx_input->GetInputValue(&s_xx,&gauss); + s_xy_input->GetInputValue(&s_xy,&gauss); + s_yy_input->GetInputValue(&s_yy,&gauss); + B_input->GetInputValue(&B,&gauss); + n_input->GetInputValue(&n,&gauss); vel=sqrt(vx*vx+vy*vy)+1.e-14; @@ -494,6 +491,4 @@ this->AddInput(BasalCrevasseEnum,&basal_crevasse[0],P1DGEnum); this->AddInput(CrevasseDepthEnum,&crevasse_depth[0],P1DGEnum); - - delete gauss; } /*}}}*/ @@ -517,16 +512,16 @@ /* Start looping on the number of vertices: */ - GaussTria* gauss=new GaussTria(); + GaussTria gauss; for (int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /* Get the value we need*/ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); + vx_input->GetInputValue(&vx,&gauss); + vy_input->GetInputValue(&vy,&gauss); vel=vx*vx+vy*vy; - strainparallel_input->GetInputValue(&strainparallel,gauss); - strainperpendicular_input->GetInputValue(&strainperpendicular,gauss); - levermanncoeff_input->GetInputValue(&propcoeff,gauss); - bs_input->GetInputValue(&bed,gauss); + strainparallel_input->GetInputValue(&strainparallel,&gauss); + strainperpendicular_input->GetInputValue(&strainperpendicular,&gauss); + levermanncoeff_input->GetInputValue(&propcoeff,&gauss); + bs_input->GetInputValue(&bed,&gauss); /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */ @@ -545,7 +540,4 @@ this->AddInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum); this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); - - /*Clean up and return*/ - delete gauss; }/*}}}*/ void Tria::CalvingFluxLevelset(){/*{{{*/ @@ -852,5 +844,4 @@ IssmDouble tau_1[NUMVERTICES]; IssmDouble tau_2[NUMVERTICES]; - GaussTria* gauss=NULL; int domaintype,dim=2; @@ -868,19 +859,19 @@ /* Start looping on the number of vertices: */ - gauss=new GaussTria(); + GaussTria gauss; for (int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /*Compute strain rate and viscosity: */ - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); switch(approximation){ case SSAApproximationEnum: - this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input); + this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0],&gauss,vx_input,vy_input); break; case HOApproximationEnum: - this->material->ViscosityHO(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input); + this->material->ViscosityHO(&viscosity,dim,&xyz_list[0][0],&gauss,vx_input,vy_input); break; case FSApproximationEnum: - this->material->ViscosityFS(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input,NULL); + this->material->ViscosityFS(&viscosity,dim,&xyz_list[0][0],&gauss,vx_input,vy_input,NULL); break; default: @@ -908,7 +899,4 @@ this->AddInput(DeviatoricStress1Enum,&tau_1[0],P1DGEnum); this->AddInput(DeviatoricStress2Enum,&tau_2[0],P1DGEnum); - - /*Clean up and return*/ - delete gauss; } /*}}}*/ @@ -921,5 +909,4 @@ IssmDouble strain_xy[NUMVERTICES]; IssmDouble vorticity_xy[NUMVERTICES]; - GaussTria* gauss=NULL; /* Get node coordinates and dof list: */ @@ -931,10 +918,10 @@ /* Start looping on the number of vertices: */ - gauss=new GaussTria(); + GaussTria gauss; for (int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /*Compute strain rate and vorticity rate: */ - this->StrainRateESA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateESA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); /*Compute Stress*/ @@ -950,7 +937,4 @@ this->AddInput(EsaStrainratexyEnum,&strain_xy[0],P1DGEnum); this->AddInput(EsaRotationrateEnum,&vorticity_xy[0],P1DGEnum); - - /*Clean up and return*/ - delete gauss; } /*}}}*/ @@ -1026,5 +1010,4 @@ IssmDouble sigma_xz[NUMVERTICES]={0,0,0}; IssmDouble sigma_yz[NUMVERTICES]={0,0,0}; - GaussTria* gauss=NULL; int domaintype,dim=2; @@ -1040,12 +1023,12 @@ /* Start looping on the number of vertices: */ - gauss=new GaussTria(); + GaussTria gauss; for (int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /*Compute strain rate viscosity and pressure: */ - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input); - pressure_input->GetInputValue(&pressure,gauss); + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); + this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0],&gauss,vx_input,vy_input); + pressure_input->GetInputValue(&pressure,&gauss); /*Compute Stress*/ @@ -1062,7 +1045,4 @@ this->AddInput(StressTensoryzEnum,&sigma_yz[0],P1DGEnum); this->AddInput(StressTensorzzEnum,&sigma_zz[0],P1DGEnum); - - /*Clean up and return*/ - delete gauss; } /*}}}*/ @@ -1176,15 +1156,14 @@ control_gradient->Serve(NUMVERTICES,&lidlist[0]); - GaussTria* gauss=new GaussTria(); + GaussTria gauss; for (int iv=0;ivGaussVertex(iv); - - control_value->GetInputValue(&value,gauss); - control_gradient->GetInputValue(&gradient,gauss); + gauss.GaussVertex(iv); + + control_value->GetInputValue(&value,&gauss); + control_gradient->GetInputValue(&gradient,&gauss); values[iv] = reCast(value)/reCast(connectivity[iv]); gradients[iv] = reCast(gradient)/reCast(connectivity[iv]); } - delete gauss; vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL); @@ -1414,14 +1393,14 @@ /*1. Recover stresses at the base*/ - GaussTria* gauss=new GaussTria(); + GaussTria gauss; for (int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /*Compute strain rate viscosity and pressure: */ - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->material->ViscosityFS(&viscosity,2,&xyz_list[0][0],gauss,vx_input,vy_input,NULL); + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); + this->material->ViscosityFS(&viscosity,2,&xyz_list[0][0],&gauss,vx_input,vy_input,NULL); /*FIXME: this is for Hongju only*/ - // pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]); - // if (pressure[iv]/pressureice[iv]>1) pressure[iv]=pressureice[iv]; + // pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]); + // if (pressure[iv]/pressureice[iv]>1) pressure[iv]=pressureice[iv]; /*Compute Stress*/ @@ -1450,7 +1429,4 @@ } } - - /*clean up*/ - delete gauss; } /*}}}*/ @@ -2155,9 +2131,7 @@ int index = this->GetNodeIndex(node); - GaussTria* gauss=new GaussTria(); - gauss->GaussNode(this->element_type,index); - - input->GetInputValue(pvalue,gauss); - delete gauss; + GaussTria gauss; + gauss.GaussNode(this->element_type,index); + input->GetInputValue(pvalue,&gauss); } /*}}}*/ @@ -2169,9 +2143,7 @@ int index = this->GetVertexIndex(vertex); - GaussTria* gauss=new GaussTria(); - gauss->GaussVertex(index); - - input->GetInputValue(pvalue,gauss); - delete gauss; + GaussTria gauss; + gauss.GaussVertex(index); + input->GetInputValue(pvalue,&gauss); } /*}}}*/ @@ -3316,8 +3288,8 @@ /*get area coordinates of 0 and 1 locations: */ - GaussTria* gauss_1=new GaussTria(); - gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]); - GaussTria* gauss_2=new GaussTria(); - gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]); + GaussTria gauss_1; + gauss_1.GaussFromCoords(x1,y1,&xyz_list[0][0]); + GaussTria gauss_2; + gauss_2.GaussFromCoords(x2,y2,&xyz_list[0][0]); /*Get segment length and normal (needs to be properly oriented)*/ @@ -3340,10 +3312,10 @@ } - thickness_input->GetInputValue(&h1, gauss_1); - thickness_input->GetInputValue(&h2, gauss_2); - vx_input->GetInputValue(&vx1,gauss_1); - vx_input->GetInputValue(&vx2,gauss_2); - vy_input->GetInputValue(&vy1,gauss_1); - vy_input->GetInputValue(&vy2,gauss_2); + thickness_input->GetInputValue(&h1, &gauss_1); + thickness_input->GetInputValue(&h2, &gauss_2); + vx_input->GetInputValue(&vx1,&gauss_1); + vx_input->GetInputValue(&vx2,&gauss_2); + vy_input->GetInputValue(&vy1,&gauss_1); + vy_input->GetInputValue(&vy2,&gauss_2); mass_flux= rho_ice*length*( @@ -3353,6 +3325,4 @@ /*clean up and return:*/ - delete gauss_1; - delete gauss_2; return mass_flux; } @@ -3540,12 +3510,12 @@ /* Start looping on the number of vertices: */ - GaussTria* gauss=new GaussTria(); + GaussTria gauss; for(int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /* Advection */ - vx_input->GetInputValue(&v[0],gauss); - vy_input->GetInputValue(&v[1],gauss); - gr_input->GetInputValue(&groundedice,gauss); + vx_input->GetInputValue(&v[0],&gauss); + vy_input->GetInputValue(&v[1],&gauss); + gr_input->GetInputValue(&groundedice,&gauss); /*Get calving speed*/ @@ -3553,8 +3523,8 @@ case DefaultCalvingEnum: case CalvingVonmisesEnum: - lsf_slopex_input->GetInputValue(&dlsf[0],gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); - calvingrate_input->GetInputValue(&calvingrate,gauss); - meltingrate_input->GetInputValue(&meltingrate,gauss); + lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); + calvingrate_input->GetInputValue(&calvingrate,&gauss); + meltingrate_input->GetInputValue(&meltingrate,&gauss); if(groundedice<0) meltingrate = 0.; @@ -3574,7 +3544,7 @@ case CalvingLevermannEnum: - calvingratex_input->GetInputValue(&c[0],gauss); - if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss); - meltingrate_input->GetInputValue(&meltingrate,gauss); + calvingratex_input->GetInputValue(&c[0],&gauss); + if(dim==2) calvingratey_input->GetInputValue(&c[1],&gauss); + meltingrate_input->GetInputValue(&meltingrate,&gauss); norm_calving=0.; for(i=0;iGetInputValue(&dlsf[0],gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); - meltingrate_input->GetInputValue(&meltingrate,gauss); + lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); + meltingrate_input->GetInputValue(&meltingrate,&gauss); norm_dlsf=0.; @@ -3605,7 +3575,7 @@ case CalvingHabEnum: - lsf_slopex_input->GetInputValue(&dlsf[0],gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); - meltingrate_input->GetInputValue(&meltingrate,gauss); + lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); + meltingrate_input->GetInputValue(&meltingrate,&gauss); norm_dlsf=0.; @@ -3626,7 +3596,7 @@ case CalvingCrevasseDepthEnum: - lsf_slopex_input->GetInputValue(&dlsf[0],gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); - meltingrate_input->GetInputValue(&meltingrate,gauss); + lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); + meltingrate_input->GetInputValue(&meltingrate,&gauss); if(groundedice<0) meltingrate = 0.; @@ -3650,9 +3620,9 @@ case CalvingDev2Enum: { - lsf_slopex_input->GetInputValue(&dlsf[0],gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); - calvingrate_input->GetInputValue(&calvingrate,gauss); - meltingrate_input->GetInputValue(&meltingrate,gauss); - gr_input->GetInputValue(&groundedice,gauss); + lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); + calvingrate_input->GetInputValue(&calvingrate,&gauss); + meltingrate_input->GetInputValue(&meltingrate,&gauss); + gr_input->GetInputValue(&groundedice,&gauss); //idea: no retreat on ice above critical calving height "calvinghaf" . Limit using regularized Heaviside function. @@ -3705,7 +3675,4 @@ this->AddInput(MovingFrontalVxEnum,&movingfrontvx[0],P1DGEnum); this->AddInput(MovingFrontalVyEnum,&movingfrontvy[0],P1DGEnum); - - /*Clean up and return*/ - delete gauss; } /*}}}*/ @@ -3828,5 +3795,5 @@ int found = 0; IssmDouble value; - GaussTria *gauss = NULL; + GaussTria gauss; /*First, serarch the input: */ @@ -3840,6 +3807,6 @@ if(data){ /*ok, we are good. retrieve value of input at vertex :*/ - gauss=new GaussTria(); gauss->GaussVertex(i); - data->GetInputValue(&value,gauss); + gauss.GaussVertex(i); + data->GetInputValue(&value,&gauss); found=1; break; @@ -3847,7 +3814,4 @@ } } - - /*clean-up*/ - delete gauss; if(found)*pvalue=value; @@ -4006,11 +3970,11 @@ /*Loop over basal nodes and update their CS*/ - GaussTria* gauss = new GaussTria(); + GaussTria gauss; for(int i=0;iNumberofNodesVelocity();i++){ if(vertexonbase[i]==1){ - gauss->GaussNode(this->VelocityInterpolation(),i); - slope_input->GetInputValue(&slope,gauss); - groundedicelevelset_input->GetInputValue(&groundedice,gauss); + gauss.GaussNode(this->VelocityInterpolation(),i); + slope_input->GetInputValue(&slope,&gauss); + groundedicelevelset_input->GetInputValue(&groundedice,&gauss); IssmDouble theta = atan(slope); @@ -4033,5 +3997,4 @@ /*cleanup*/ xDelete(vertexonbase); - delete gauss; } /*}}}*/ @@ -4079,12 +4042,12 @@ /* Start looping on the number of vertices: */ - GaussTria* gauss=new GaussTria(); + GaussTria gauss; for(int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /* Get variables */ - bed_input->GetInputValue(&bed,gauss); - qsg_input->GetInputValue(&qsg,gauss); - TF_input->GetInputValue(&TF,gauss); + bed_input->GetInputValue(&bed,&gauss); + qsg_input->GetInputValue(&qsg,&gauss); + TF_input->GetInputValue(&TF,&gauss); if(basin_icefront_area[reCast(basinid[iv])-1]==0.) meltrates[iv]=0.; @@ -4106,5 +4069,4 @@ /*Cleanup and return*/ xDelete(basin_icefront_area); - delete gauss; } /*}}}*/ @@ -4294,5 +4256,4 @@ IssmDouble epsilon[3]; - GaussTria* gauss=NULL; IssmDouble vx,vy,vel; IssmDouble strainxx; @@ -4310,15 +4271,15 @@ /* Start looping on the number of vertices: */ - gauss=new GaussTria(); + GaussTria gauss; for (int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /* Get the value we need*/ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); + vx_input->GetInputValue(&vx,&gauss); + vy_input->GetInputValue(&vy,&gauss); vel=vx*vx+vy*vy; /*Compute strain rate viscosity and pressure: */ - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); strainxx=epsilon[0]; strainyy=epsilon[1]; @@ -4331,12 +4292,8 @@ /*Add input*/ this->AddInput(StrainRateparallelEnum,&strainparallel[0],P1DGEnum); - - /*Clean up and return*/ - delete gauss; } /*}}}*/ void Tria::StrainRateperpendicular(){/*{{{*/ - GaussTria* gauss=NULL; IssmDouble epsilon[3]; IssmDouble vx,vy,vel; @@ -4355,15 +4312,15 @@ /* Start looping on the number of vertices: */ - gauss=new GaussTria(); + GaussTria gauss; for (int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /* Get the value we need*/ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); + vx_input->GetInputValue(&vx,&gauss); + vy_input->GetInputValue(&vy,&gauss); vel=vx*vx+vy*vy; /*Compute strain rate viscosity and pressure: */ - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); strainxx=epsilon[0]; strainyy=epsilon[1]; @@ -4376,7 +4333,4 @@ /*Add input*/ this->AddInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1DGEnum); - - /*Clean up and return*/ - delete gauss; } /*}}}*/ @@ -4996,14 +4950,13 @@ Input* onbase = this->GetInput(MeshVertexonbaseEnum); _assert_(onbase); - GaussTria* gauss=new GaussTria(); + GaussTria gauss; for(int iv=0;ivNumberofNodes(this->element_type);iv++){ - gauss->GaussNode(this->element_type,iv); - onbase->GetInputValue(&isonbase,gauss); + gauss.GaussNode(this->element_type,iv); + onbase->GetInputValue(&isonbase,&gauss); if(isonbase==1.){ - input->GetInputValue(&value,gauss); + input->GetInputValue(&value,&gauss); this->nodes[iv]->ApplyConstraint(0,value); } } - delete gauss; } @@ -5020,14 +4973,13 @@ Input* onsurf = this->GetInput(MeshVertexonsurfaceEnum); _assert_(onsurf); - GaussTria* gauss=new GaussTria(); + GaussTria gauss; for(int iv=0;ivNumberofNodes(this->element_type);iv++){ - gauss->GaussNode(this->element_type,iv); - onsurf->GetInputValue(&isonsurface,gauss); + gauss.GaussNode(this->element_type,iv); + onsurf->GetInputValue(&isonsurface,&gauss); if(isonsurface==1.){ - input->GetInputValue(&value,gauss); + input->GetInputValue(&value,&gauss); this->nodes[iv]->ApplyConstraint(0,value); } } - delete gauss; } /*}}}*/