![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 25934
- Timestamp:
- 01/07/21 19:54:55 (4 years ago)
- Location:
- issm/trunk-jpl/src/c/classes/Elements
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
r25839 r25934 299 299 300 300 /* Start looping on the number of vertices: */ 301 GaussPenta * gauss=new GaussPenta();301 GaussPenta gauss; 302 302 for(int iv=0;iv<3;iv++){ 303 gauss ->GaussVertex(iv);303 gauss.GaussVertex(iv); 304 304 305 305 /*Get velocity components and thickness*/ 306 B_input->GetInputValue(&B, gauss);307 n_input->GetInputValue(&n, gauss);308 vx_input->GetInputValue(&vx, gauss);309 vy_input->GetInputValue(&vy, gauss);310 gr_input->GetInputValue(&groundedice, gauss);311 bs_input->GetInputValue(&bed, gauss);312 smax_fl_input->GetInputValue(&sigma_max_floating, gauss);313 smax_gr_input->GetInputValue(&sigma_max_grounded, gauss);306 B_input->GetInputValue(&B,&gauss); 307 n_input->GetInputValue(&n,&gauss); 308 vx_input->GetInputValue(&vx,&gauss); 309 vy_input->GetInputValue(&vy,&gauss); 310 gr_input->GetInputValue(&groundedice,&gauss); 311 bs_input->GetInputValue(&bed,&gauss); 312 smax_fl_input->GetInputValue(&sigma_max_floating,&gauss); 313 smax_gr_input->GetInputValue(&sigma_max_grounded,&gauss); 314 314 vel=sqrt(vx*vx+vy*vy)+1.e-14; 315 sl_input->GetInputValue(&sealevel, gauss);315 sl_input->GetInputValue(&sealevel,&gauss); 316 316 317 317 /*Compute strain rate and viscosity: */ 318 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0], gauss,vx_input,vy_input);318 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); 319 319 320 320 /*Get Eigen values*/ … … 359 359 this->InputExtrude(CalvingCalvingrateEnum,-1); 360 360 this->InputExtrude(SigmaVMEnum,-1); 361 362 /*Clean up and return*/363 delete gauss;364 361 } 365 362 /*}}}*/ … … 382 379 383 380 /* Start looping on the number of vertices: */ 384 GaussPenta * gauss=new GaussPenta();385 for 386 gauss ->GaussVertex(iv);381 GaussPenta gauss; 382 for(int iv=0;iv<NUMVERTICES;iv++){ 383 gauss.GaussVertex(iv); 387 384 388 385 /* Get the value we need*/ 389 vx_input->GetInputValue(&vx, gauss);390 vy_input->GetInputValue(&vy, gauss);386 vx_input->GetInputValue(&vx,&gauss); 387 vy_input->GetInputValue(&vy,&gauss); 391 388 vel=vx*vx+vy*vy; 392 strainparallel_input->GetInputValue(&strainparallel, gauss);393 strainperpendicular_input->GetInputValue(&strainperpendicular, gauss);394 levermanncoeff_input->GetInputValue(&propcoeff, gauss);389 strainparallel_input->GetInputValue(&strainparallel,&gauss); 390 strainperpendicular_input->GetInputValue(&strainperpendicular,&gauss); 391 levermanncoeff_input->GetInputValue(&propcoeff,&gauss); 395 392 396 393 /*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 404 this->AddBasalInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum); 408 405 this->AddBasalInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); 409 410 /*Clean up and return*/411 delete gauss;412 406 }/*}}}*/ 413 407 void Penta::CalvingFluxLevelset(){/*{{{*/ … … 761 755 IssmDouble tau_yz[NUMVERTICES]; 762 756 IssmDouble tau_eff[NUMVERTICES]; 763 GaussPenta* gauss=NULL;764 757 765 758 /* Get node coordinates and dof list: */ … … 772 765 773 766 /* Start looping on the number of vertices: */ 774 gauss=new GaussPenta();767 GaussPenta gauss; 775 768 for (int iv=0;iv<NUMVERTICES;iv++){ 776 gauss ->GaussVertex(iv);769 gauss.GaussVertex(iv); 777 770 778 771 /*Compute strain rate viscosity and pressure: */ 779 this->StrainRateFS(&epsilon[0],&xyz_list[0][0], gauss,vx_input,vy_input,vz_input);780 this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0], gauss,vx_input,vy_input,vz_input);772 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); 773 this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); 781 774 782 775 /*Compute Stress*/ … … 802 795 this->AddInput(DeviatoricStresszzEnum,&tau_zz[0],P1DGEnum); 803 796 this->AddInput(DeviatoricStresseffectiveEnum,&tau_eff[0],P1DGEnum); 804 805 /*Clean up and return*/806 delete gauss;807 797 } 808 798 /*}}}*/ … … 818 808 IssmDouble sigma_xz[NUMVERTICES]; 819 809 IssmDouble sigma_yz[NUMVERTICES]; 820 GaussPenta* gauss=NULL;821 810 822 811 /* Get node coordinates and dof list: */ … … 830 819 831 820 /* Start looping on the number of vertices: */ 832 gauss=new GaussPenta();833 for 834 gauss ->GaussVertex(iv);821 GaussPenta gauss; 822 for(int iv=0;iv<NUMVERTICES;iv++){ 823 gauss.GaussVertex(iv); 835 824 836 825 /*Compute strain rate viscosity and pressure: */ 837 this->StrainRateFS(&epsilon[0],&xyz_list[0][0], gauss,vx_input,vy_input,vz_input);838 this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0], gauss,vx_input,vy_input,vz_input);839 pressure_input->GetInputValue(&pressure, gauss);826 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); 827 this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); 828 pressure_input->GetInputValue(&pressure,&gauss); 840 829 841 830 /*Compute Stress*/ … … 855 844 this->AddInput(StressTensoryzEnum,&sigma_yz[0],P1DGEnum); 856 845 this->AddInput(StressTensorzzEnum,&sigma_zz[0],P1DGEnum); 857 858 /*Clean up and return*/859 delete gauss;860 846 } 861 847 /*}}}*/ … … 952 938 control_gradient->Serve(NUMVERTICES,&lidlist[0]); 953 939 954 GaussPenta * gauss=new GaussPenta();940 GaussPenta gauss; 955 941 for (int iv=0;iv<NUMVERTICES;iv++){ 956 gauss ->GaussVertex(iv);957 958 control_value->GetInputValue(&value, gauss);959 control_gradient->GetInputValue(&gradient, gauss);942 gauss.GaussVertex(iv); 943 944 control_value->GetInputValue(&value,&gauss); 945 control_gradient->GetInputValue(&gradient,&gauss); 960 946 961 947 values[iv] = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]); 962 948 gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]); 963 949 } 964 delete gauss;965 950 966 951 vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL); … … 1139 1124 1140 1125 /*1. Recover stresses at the base*/ 1141 GaussPenta * gauss=new GaussPenta();1126 GaussPenta gauss; 1142 1127 for (int iv=0;iv<NUMVERTICES;iv++){ 1143 gauss ->GaussVertex(iv);1128 gauss.GaussVertex(iv); 1144 1129 1145 1130 /*Compute strain rate viscosity and pressure: */ 1146 this->StrainRateFS(&epsilon[0],&xyz_list[0][0], gauss,vx_input,vy_input,vz_input);1147 this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0], gauss,vx_input,vy_input,vz_input);1131 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); 1132 this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); 1148 1133 /*FIXME: this is for Hongju only*/ 1149 1134 //pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]); … … 1180 1165 } 1181 1166 } 1182 1183 /*clean up*/1184 delete gauss;1185 1167 } 1186 1168 /*}}}*/ … … 1631 1613 int index = this->GetNodeIndex(node); 1632 1614 1633 GaussPenta* gauss=new GaussPenta(); 1634 gauss->GaussNode(this->element_type,index); 1635 1636 input->GetInputValue(pvalue,gauss); 1637 delete gauss; 1615 GaussPenta gauss; 1616 gauss.GaussNode(this->element_type,index); 1617 input->GetInputValue(pvalue,&gauss); 1638 1618 } 1639 1619 /*}}}*/ … … 1645 1625 int index = this->GetVertexIndex(vertex); 1646 1626 1647 GaussPenta* gauss=new GaussPenta(); 1648 gauss->GaussVertex(index); 1649 1650 input->GetInputValue(pvalue,gauss); 1651 delete gauss; 1627 GaussPenta gauss; 1628 gauss.GaussVertex(index); 1629 input->GetInputValue(pvalue,&gauss); 1652 1630 } 1653 1631 /*}}}*/ … … 2851 2829 2852 2830 /* Start looping on the number of vertices: */ 2853 GaussPenta * gauss=new GaussPenta();2831 GaussPenta gauss; 2854 2832 for(int iv=0;iv<NUMVERTICES;iv++){ 2855 gauss ->GaussVertex(iv);2833 gauss.GaussVertex(iv); 2856 2834 2857 2835 /* Advection */ 2858 vx_input->GetInputValue(&v[0], gauss);2859 vy_input->GetInputValue(&v[1], gauss);2860 gr_input->GetInputValue(&groundedice, gauss);2836 vx_input->GetInputValue(&v[0],&gauss); 2837 vy_input->GetInputValue(&v[1],&gauss); 2838 gr_input->GetInputValue(&groundedice,&gauss); 2861 2839 2862 2840 /*Get calving speed*/ … … 2864 2842 case DefaultCalvingEnum: 2865 2843 case CalvingVonmisesEnum: 2866 lsf_slopex_input->GetInputValue(&dlsf[0], gauss);2867 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1], gauss);2868 calvingrate_input->GetInputValue(&calvingrate, gauss);2869 meltingrate_input->GetInputValue(&meltingrate, gauss);2844 lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); 2845 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); 2846 calvingrate_input->GetInputValue(&calvingrate,&gauss); 2847 meltingrate_input->GetInputValue(&meltingrate,&gauss); 2870 2848 2871 2849 if(groundedice<0) meltingrate = 0.; … … 2886 2864 2887 2865 case CalvingLevermannEnum: 2888 calvingratex_input->GetInputValue(&c[0], gauss);2889 if(dim==2) calvingratey_input->GetInputValue(&c[1], gauss);2890 meltingrate_input->GetInputValue(&meltingrate, gauss);2866 calvingratex_input->GetInputValue(&c[0],&gauss); 2867 if(dim==2) calvingratey_input->GetInputValue(&c[1],&gauss); 2868 meltingrate_input->GetInputValue(&meltingrate,&gauss); 2891 2869 norm_calving=0.; 2892 2870 for(i=0;i<dim;i++) norm_calving+=pow(c[i],2); … … 2896 2874 2897 2875 case CalvingMinthicknessEnum: 2898 lsf_slopex_input->GetInputValue(&dlsf[0], gauss);2899 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1], gauss);2900 meltingrate_input->GetInputValue(&meltingrate, gauss);2876 lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); 2877 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); 2878 meltingrate_input->GetInputValue(&meltingrate,&gauss); 2901 2879 2902 2880 norm_dlsf=0.; … … 2917 2895 2918 2896 case CalvingHabEnum: 2919 lsf_slopex_input->GetInputValue(&dlsf[0], gauss);2920 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1], gauss);2921 meltingrate_input->GetInputValue(&meltingrate, gauss);2897 lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); 2898 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); 2899 meltingrate_input->GetInputValue(&meltingrate,&gauss); 2922 2900 2923 2901 norm_dlsf=0.; … … 2938 2916 2939 2917 case CalvingCrevasseDepthEnum: 2940 lsf_slopex_input->GetInputValue(&dlsf[0], gauss);2941 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1], gauss);2942 meltingrate_input->GetInputValue(&meltingrate, gauss);2918 lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); 2919 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); 2920 meltingrate_input->GetInputValue(&meltingrate,&gauss); 2943 2921 2944 2922 if(groundedice<0) meltingrate = 0.; … … 2962 2940 case CalvingDev2Enum: 2963 2941 { 2964 lsf_slopex_input->GetInputValue(&dlsf[0], gauss);2965 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1], gauss);2966 calvingrate_input->GetInputValue(&calvingrate, gauss);2967 meltingrate_input->GetInputValue(&meltingrate, gauss);2968 gr_input->GetInputValue(&groundedice, gauss);2942 lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); 2943 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); 2944 calvingrate_input->GetInputValue(&calvingrate,&gauss); 2945 meltingrate_input->GetInputValue(&meltingrate,&gauss); 2946 gr_input->GetInputValue(&groundedice,&gauss); 2969 2947 2970 2948 //idea: no retreat on ice above critical calving height "calvinghaf" . Limit using regularized Heaviside function. … … 3021 2999 this->InputExtrude(MovingFrontalVxEnum,-1); 3022 3000 this->InputExtrude(MovingFrontalVyEnum,-1); 3023 /*Clean up and return*/3024 delete gauss;3025 3001 } 3026 3002 /*}}}*/ … … 3142 3118 int found=0; 3143 3119 IssmDouble value; 3144 GaussPenta * gauss=NULL;3120 GaussPenta gauss; 3145 3121 3146 3122 /*First, serarch the input: */ … … 3154 3130 if(data){ 3155 3131 /*ok, we are good. retrieve value of input at vertex :*/ 3156 gauss =new GaussPenta(); gauss->GaussVertex(i);3157 data->GetInputValue(&value, gauss);3132 gauss.GaussVertex(i); 3133 data->GetInputValue(&value,&gauss); 3158 3134 found=1; 3159 3135 break; … … 3162 3138 } 3163 3139 3164 delete gauss;3165 3140 if(found)*pvalue=value; 3166 3141 return found; … … 3369 3344 3370 3345 /*Loop over basal nodes and update their CS*/ 3371 GaussPenta * gauss = new GaussPenta();3346 GaussPenta gauss; 3372 3347 for(int i=0;i<numindices;i++){//FIXME 3373 gauss ->GaussNode(this->VelocityInterpolation(),indices[i]);3374 3375 slopex_input->GetInputValue(&slopex, gauss);3376 slopey_input->GetInputValue(&slopey, gauss);3377 groundedicelevelset_input->GetInputValue(&groundedice, gauss);3348 gauss.GaussNode(this->VelocityInterpolation(),indices[i]); 3349 3350 slopex_input->GetInputValue(&slopex,&gauss); 3351 slopey_input->GetInputValue(&slopey,&gauss); 3352 groundedicelevelset_input->GetInputValue(&groundedice,&gauss); 3378 3353 3379 3354 /*New X axis New Z axis*/ … … 3406 3381 /*cleanup*/ 3407 3382 xDelete<int>(indices); 3408 delete gauss;3409 3383 } 3410 3384 /*}}}*/ … … 3455 3429 3456 3430 /* Start looping on the number of vertices: */ 3457 GaussPenta * gauss=new GaussPenta();3431 GaussPenta gauss; 3458 3432 for(int iv=0;iv<NUMVERTICES2D;iv++){ 3459 gauss ->GaussVertex(iv);3433 gauss.GaussVertex(iv); 3460 3434 3461 3435 /* Get variables */ 3462 bed_input->GetInputValue(&bed, gauss);3463 qsg_input->GetInputValue(&qsg, gauss);3464 TF_input->GetInputValue(&TF, gauss);3436 bed_input->GetInputValue(&bed,&gauss); 3437 qsg_input->GetInputValue(&qsg,&gauss); 3438 TF_input->GetInputValue(&TF,&gauss); 3465 3439 3466 3440 if(basinid[iv]==0 || basin_icefront_area[reCast<int>(basinid[iv])-1]==0.) meltrates[iv]=0.; … … 3484 3458 /*Cleanup and return*/ 3485 3459 xDelete<IssmDouble>(basin_icefront_area); 3486 delete gauss;3487 3460 } 3488 3461 /*}}}*/ … … 3692 3665 3693 3666 IssmDouble epsilon[6]; 3694 GaussPenta* gauss=NULL;3695 3667 IssmDouble vx,vy,vel; 3696 3668 IssmDouble strainxx; … … 3709 3681 3710 3682 /* Start looping on the number of vertices: */ 3711 gauss=new GaussPenta();3683 GaussPenta gauss; 3712 3684 for (int iv=0;iv<NUMVERTICES;iv++){ 3713 gauss ->GaussVertex(iv);3685 gauss.GaussVertex(iv); 3714 3686 3715 3687 /* Get the value we need*/ 3716 vx_input->GetInputValue(&vx, gauss);3717 vy_input->GetInputValue(&vy, gauss);3688 vx_input->GetInputValue(&vx,&gauss); 3689 vy_input->GetInputValue(&vy,&gauss); 3718 3690 vel=vx*vx+vy*vy; 3719 3691 3720 3692 /*Compute strain rate and viscosity: */ 3721 this->StrainRateFS(&epsilon[0],&xyz_list[0][0], gauss,vx_input,vy_input,vz_input);3693 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); 3722 3694 strainxx=epsilon[0]; 3723 3695 strainyy=epsilon[1]; … … 3730 3702 /*Add input*/ 3731 3703 this->AddInput(StrainRateparallelEnum,&strainparallel[0],P1DGEnum); 3732 3733 /*Clean up and return*/3734 delete gauss;3735 3704 } 3736 3705 /*}}}*/ … … 3738 3707 3739 3708 IssmDouble epsilon[6]; 3740 GaussPenta* gauss=NULL;3741 3709 IssmDouble vx,vy,vel; 3742 3710 IssmDouble strainxx; … … 3755 3723 3756 3724 /* Start looping on the number of vertices: */ 3757 gauss=new GaussPenta();3725 GaussPenta gauss; 3758 3726 for (int iv=0;iv<NUMVERTICES;iv++){ 3759 gauss ->GaussVertex(iv);3727 gauss.GaussVertex(iv); 3760 3728 3761 3729 /* Get the value we need*/ 3762 vx_input->GetInputValue(&vx, gauss);3763 vy_input->GetInputValue(&vy, gauss);3730 vx_input->GetInputValue(&vx,&gauss); 3731 vy_input->GetInputValue(&vy,&gauss); 3764 3732 vel=vx*vx+vy*vy; 3765 3733 3766 3734 /*Compute strain rate and viscosity: */ 3767 this->StrainRateFS(&epsilon[0],&xyz_list[0][0], gauss,vx_input,vy_input,vz_input);3735 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input,vz_input); 3768 3736 strainxx=epsilon[0]; 3769 3737 strainyy=epsilon[1]; … … 3776 3744 /*Add input*/ 3777 3745 this->AddInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1DGEnum); 3778 3779 /*Clean up and return*/3780 delete gauss;3781 3746 } 3782 3747 /*}}}*/ … … 4779 4744 Input* onbase = this->GetInput(MeshVertexonbaseEnum); _assert_(onbase); 4780 4745 4781 GaussPenta * gauss=new GaussPenta();4746 GaussPenta gauss; 4782 4747 for(int iv=0;iv<this->NumberofNodes(this->element_type);iv++){ 4783 gauss ->GaussNode(this->element_type,iv);4784 onbase->GetInputValue(&isonbase, gauss);4748 gauss.GaussNode(this->element_type,iv); 4749 onbase->GetInputValue(&isonbase,&gauss); 4785 4750 if(isonbase==1.){ 4786 input->GetInputValue(&value, gauss);4751 input->GetInputValue(&value,&gauss); 4787 4752 this->nodes[iv]->ApplyConstraint(0,value); 4788 4753 } 4789 4754 } 4790 delete gauss; 4791 4792 } 4793 /*}}}*/ 4755 }/*}}}*/ 4794 4756 void Penta::UpdateConstraintsExtrudeFromTop(void){/*{{{*/ 4795 4757 … … 4803 4765 Input* input = this->GetInput(extrusioninput); _assert_(extrusioninput); 4804 4766 4805 GaussPenta * gauss=new GaussPenta();4767 GaussPenta gauss; 4806 4768 for(int i=0;i<3;i++){ 4807 gauss ->GaussNode(P1Enum,indices[i]);4808 input->GetInputValue(&value, gauss);4769 gauss.GaussNode(P1Enum,indices[i]); 4770 input->GetInputValue(&value,&gauss); 4809 4771 this->nodes[indices[i]]->ApplyConstraint(0,value); 4810 4772 } 4811 delete gauss;4812 4773 4813 4774 } -
issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
r25905 r25934 339 339 340 340 /* Start looping on the number of vertices: */ 341 GaussTria * gauss=new GaussTria();341 GaussTria gauss; 342 342 for(int iv=0;iv<NUMVERTICES;iv++){ 343 gauss ->GaussVertex(iv);343 gauss.GaussVertex(iv); 344 344 345 345 /*Get velocity components and thickness*/ 346 B_input->GetInputValue(&B, gauss);347 n_input->GetInputValue(&n, gauss);348 vx_input->GetInputValue(&vx, gauss);349 vy_input->GetInputValue(&vy, gauss);350 gr_input->GetInputValue(&groundedice, gauss);351 bs_input->GetInputValue(&bed, gauss);352 smax_fl_input->GetInputValue(&sigma_max_floating, gauss);353 smax_gr_input->GetInputValue(&sigma_max_grounded, gauss);346 B_input->GetInputValue(&B,&gauss); 347 n_input->GetInputValue(&n,&gauss); 348 vx_input->GetInputValue(&vx,&gauss); 349 vy_input->GetInputValue(&vy,&gauss); 350 gr_input->GetInputValue(&groundedice,&gauss); 351 bs_input->GetInputValue(&bed,&gauss); 352 smax_fl_input->GetInputValue(&sigma_max_floating,&gauss); 353 smax_gr_input->GetInputValue(&sigma_max_grounded,&gauss); 354 354 vel=sqrt(vx*vx+vy*vy)+1.e-14; 355 sl_input->GetInputValue(&sealevel, gauss);355 sl_input->GetInputValue(&sealevel,&gauss); 356 356 357 357 /*Compute strain rate and viscosity: */ 358 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0], gauss,vx_input,vy_input);358 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); 359 359 360 360 /*Get Eigen values*/ … … 399 399 this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); 400 400 this->AddInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum); 401 402 /*Clean up and return*/403 delete gauss;404 401 } 405 402 /*}}}*/ … … 437 434 438 435 /*Loop over all elements of this partition*/ 439 GaussTria * gauss=new GaussTria();436 GaussTria gauss; 440 437 for (int iv=0;iv<NUMVERTICES;iv++){ 441 gauss ->GaussVertex(iv);442 443 H_input->GetInputValue(&thickness, gauss);444 bed_input->GetInputValue(&bed, gauss);445 surface_input->GetInputValue(&surface, gauss);446 strainrateparallel_input->GetInputValue(&strainparallel, gauss);447 strainrateeffective_input->GetInputValue(&straineffective, gauss);448 vx_input->GetInputValue(&vx, gauss);449 vy_input->GetInputValue(&vy, gauss);450 waterheight_input->GetInputValue(&water_height, gauss);451 s_xx_input->GetInputValue(&s_xx, gauss);452 s_xy_input->GetInputValue(&s_xy, gauss);453 s_yy_input->GetInputValue(&s_yy, gauss);454 B_input->GetInputValue(&B, gauss);455 n_input->GetInputValue(&n, gauss);438 gauss.GaussVertex(iv); 439 440 H_input->GetInputValue(&thickness,&gauss); 441 bed_input->GetInputValue(&bed,&gauss); 442 surface_input->GetInputValue(&surface,&gauss); 443 strainrateparallel_input->GetInputValue(&strainparallel,&gauss); 444 strainrateeffective_input->GetInputValue(&straineffective,&gauss); 445 vx_input->GetInputValue(&vx,&gauss); 446 vy_input->GetInputValue(&vy,&gauss); 447 waterheight_input->GetInputValue(&water_height,&gauss); 448 s_xx_input->GetInputValue(&s_xx,&gauss); 449 s_xy_input->GetInputValue(&s_xy,&gauss); 450 s_yy_input->GetInputValue(&s_yy,&gauss); 451 B_input->GetInputValue(&B,&gauss); 452 n_input->GetInputValue(&n,&gauss); 456 453 457 454 vel=sqrt(vx*vx+vy*vy)+1.e-14; … … 494 491 this->AddInput(BasalCrevasseEnum,&basal_crevasse[0],P1DGEnum); 495 492 this->AddInput(CrevasseDepthEnum,&crevasse_depth[0],P1DGEnum); 496 497 delete gauss;498 493 } 499 494 /*}}}*/ … … 517 512 518 513 /* Start looping on the number of vertices: */ 519 GaussTria * gauss=new GaussTria();514 GaussTria gauss; 520 515 for (int iv=0;iv<NUMVERTICES;iv++){ 521 gauss ->GaussVertex(iv);516 gauss.GaussVertex(iv); 522 517 523 518 /* Get the value we need*/ 524 vx_input->GetInputValue(&vx, gauss);525 vy_input->GetInputValue(&vy, gauss);519 vx_input->GetInputValue(&vx,&gauss); 520 vy_input->GetInputValue(&vy,&gauss); 526 521 vel=vx*vx+vy*vy; 527 strainparallel_input->GetInputValue(&strainparallel, gauss);528 strainperpendicular_input->GetInputValue(&strainperpendicular, gauss);529 levermanncoeff_input->GetInputValue(&propcoeff, gauss);530 bs_input->GetInputValue(&bed, gauss);522 strainparallel_input->GetInputValue(&strainparallel,&gauss); 523 strainperpendicular_input->GetInputValue(&strainperpendicular,&gauss); 524 levermanncoeff_input->GetInputValue(&propcoeff,&gauss); 525 bs_input->GetInputValue(&bed,&gauss); 531 526 532 527 /*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 540 this->AddInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum); 546 541 this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); 547 548 /*Clean up and return*/549 delete gauss;550 542 }/*}}}*/ 551 543 void Tria::CalvingFluxLevelset(){/*{{{*/ … … 852 844 IssmDouble tau_1[NUMVERTICES]; 853 845 IssmDouble tau_2[NUMVERTICES]; 854 GaussTria* gauss=NULL;855 846 int domaintype,dim=2; 856 847 … … 868 859 869 860 /* Start looping on the number of vertices: */ 870 gauss=new GaussTria();861 GaussTria gauss; 871 862 for (int iv=0;iv<NUMVERTICES;iv++){ 872 gauss ->GaussVertex(iv);863 gauss.GaussVertex(iv); 873 864 874 865 /*Compute strain rate and viscosity: */ 875 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0], gauss,vx_input,vy_input);866 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); 876 867 switch(approximation){ 877 868 case SSAApproximationEnum: 878 this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0], gauss,vx_input,vy_input);869 this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0],&gauss,vx_input,vy_input); 879 870 break; 880 871 case HOApproximationEnum: 881 this->material->ViscosityHO(&viscosity,dim,&xyz_list[0][0], gauss,vx_input,vy_input);872 this->material->ViscosityHO(&viscosity,dim,&xyz_list[0][0],&gauss,vx_input,vy_input); 882 873 break; 883 874 case FSApproximationEnum: 884 this->material->ViscosityFS(&viscosity,dim,&xyz_list[0][0], gauss,vx_input,vy_input,NULL);875 this->material->ViscosityFS(&viscosity,dim,&xyz_list[0][0],&gauss,vx_input,vy_input,NULL); 885 876 break; 886 877 default: … … 908 899 this->AddInput(DeviatoricStress1Enum,&tau_1[0],P1DGEnum); 909 900 this->AddInput(DeviatoricStress2Enum,&tau_2[0],P1DGEnum); 910 911 /*Clean up and return*/912 delete gauss;913 901 } 914 902 /*}}}*/ … … 921 909 IssmDouble strain_xy[NUMVERTICES]; 922 910 IssmDouble vorticity_xy[NUMVERTICES]; 923 GaussTria* gauss=NULL;924 911 925 912 /* Get node coordinates and dof list: */ … … 931 918 932 919 /* Start looping on the number of vertices: */ 933 gauss=new GaussTria();920 GaussTria gauss; 934 921 for (int iv=0;iv<NUMVERTICES;iv++){ 935 gauss ->GaussVertex(iv);922 gauss.GaussVertex(iv); 936 923 937 924 /*Compute strain rate and vorticity rate: */ 938 this->StrainRateESA(&epsilon[0],&xyz_list[0][0], gauss,vx_input,vy_input);925 this->StrainRateESA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); 939 926 940 927 /*Compute Stress*/ … … 950 937 this->AddInput(EsaStrainratexyEnum,&strain_xy[0],P1DGEnum); 951 938 this->AddInput(EsaRotationrateEnum,&vorticity_xy[0],P1DGEnum); 952 953 /*Clean up and return*/954 delete gauss;955 939 } 956 940 /*}}}*/ … … 1026 1010 IssmDouble sigma_xz[NUMVERTICES]={0,0,0}; 1027 1011 IssmDouble sigma_yz[NUMVERTICES]={0,0,0}; 1028 GaussTria* gauss=NULL;1029 1012 int domaintype,dim=2; 1030 1013 … … 1040 1023 1041 1024 /* Start looping on the number of vertices: */ 1042 gauss=new GaussTria();1025 GaussTria gauss; 1043 1026 for (int iv=0;iv<NUMVERTICES;iv++){ 1044 gauss ->GaussVertex(iv);1027 gauss.GaussVertex(iv); 1045 1028 1046 1029 /*Compute strain rate viscosity and pressure: */ 1047 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0], gauss,vx_input,vy_input);1048 this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0], gauss,vx_input,vy_input);1049 pressure_input->GetInputValue(&pressure, gauss);1030 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); 1031 this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0],&gauss,vx_input,vy_input); 1032 pressure_input->GetInputValue(&pressure,&gauss); 1050 1033 1051 1034 /*Compute Stress*/ … … 1062 1045 this->AddInput(StressTensoryzEnum,&sigma_yz[0],P1DGEnum); 1063 1046 this->AddInput(StressTensorzzEnum,&sigma_zz[0],P1DGEnum); 1064 1065 /*Clean up and return*/1066 delete gauss;1067 1047 } 1068 1048 /*}}}*/ … … 1176 1156 control_gradient->Serve(NUMVERTICES,&lidlist[0]); 1177 1157 1178 GaussTria * gauss=new GaussTria();1158 GaussTria gauss; 1179 1159 for (int iv=0;iv<NUMVERTICES;iv++){ 1180 gauss ->GaussVertex(iv);1181 1182 control_value->GetInputValue(&value, gauss);1183 control_gradient->GetInputValue(&gradient, gauss);1160 gauss.GaussVertex(iv); 1161 1162 control_value->GetInputValue(&value,&gauss); 1163 control_gradient->GetInputValue(&gradient,&gauss); 1184 1164 1185 1165 values[iv] = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]); 1186 1166 gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]); 1187 1167 } 1188 delete gauss;1189 1168 1190 1169 vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL); … … 1414 1393 1415 1394 /*1. Recover stresses at the base*/ 1416 GaussTria * gauss=new GaussTria();1395 GaussTria gauss; 1417 1396 for (int iv=0;iv<NUMVERTICES;iv++){ 1418 gauss ->GaussVertex(iv);1397 gauss.GaussVertex(iv); 1419 1398 1420 1399 /*Compute strain rate viscosity and pressure: */ 1421 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0], gauss,vx_input,vy_input);1422 this->material->ViscosityFS(&viscosity,2,&xyz_list[0][0], gauss,vx_input,vy_input,NULL);1400 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); 1401 this->material->ViscosityFS(&viscosity,2,&xyz_list[0][0],&gauss,vx_input,vy_input,NULL); 1423 1402 /*FIXME: this is for Hongju only*/ 1424 // pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]);1425 // if (pressure[iv]/pressureice[iv]>1) pressure[iv]=pressureice[iv];1403 // pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]); 1404 // if (pressure[iv]/pressureice[iv]>1) pressure[iv]=pressureice[iv]; 1426 1405 1427 1406 /*Compute Stress*/ … … 1450 1429 } 1451 1430 } 1452 1453 /*clean up*/1454 delete gauss;1455 1431 } 1456 1432 /*}}}*/ … … 2155 2131 int index = this->GetNodeIndex(node); 2156 2132 2157 GaussTria* gauss=new GaussTria(); 2158 gauss->GaussNode(this->element_type,index); 2159 2160 input->GetInputValue(pvalue,gauss); 2161 delete gauss; 2133 GaussTria gauss; 2134 gauss.GaussNode(this->element_type,index); 2135 input->GetInputValue(pvalue,&gauss); 2162 2136 } 2163 2137 /*}}}*/ … … 2169 2143 int index = this->GetVertexIndex(vertex); 2170 2144 2171 GaussTria* gauss=new GaussTria(); 2172 gauss->GaussVertex(index); 2173 2174 input->GetInputValue(pvalue,gauss); 2175 delete gauss; 2145 GaussTria gauss; 2146 gauss.GaussVertex(index); 2147 input->GetInputValue(pvalue,&gauss); 2176 2148 } 2177 2149 /*}}}*/ … … 3316 3288 3317 3289 /*get area coordinates of 0 and 1 locations: */ 3318 GaussTria * gauss_1=new GaussTria();3319 gauss_1 ->GaussFromCoords(x1,y1,&xyz_list[0][0]);3320 GaussTria * gauss_2=new GaussTria();3321 gauss_2 ->GaussFromCoords(x2,y2,&xyz_list[0][0]);3290 GaussTria gauss_1; 3291 gauss_1.GaussFromCoords(x1,y1,&xyz_list[0][0]); 3292 GaussTria gauss_2; 3293 gauss_2.GaussFromCoords(x2,y2,&xyz_list[0][0]); 3322 3294 3323 3295 /*Get segment length and normal (needs to be properly oriented)*/ … … 3340 3312 } 3341 3313 3342 thickness_input->GetInputValue(&h1, gauss_1);3343 thickness_input->GetInputValue(&h2, gauss_2);3344 vx_input->GetInputValue(&vx1, gauss_1);3345 vx_input->GetInputValue(&vx2, gauss_2);3346 vy_input->GetInputValue(&vy1, gauss_1);3347 vy_input->GetInputValue(&vy2, gauss_2);3314 thickness_input->GetInputValue(&h1, &gauss_1); 3315 thickness_input->GetInputValue(&h2, &gauss_2); 3316 vx_input->GetInputValue(&vx1,&gauss_1); 3317 vx_input->GetInputValue(&vx2,&gauss_2); 3318 vy_input->GetInputValue(&vy1,&gauss_1); 3319 vy_input->GetInputValue(&vy2,&gauss_2); 3348 3320 3349 3321 mass_flux= rho_ice*length*( … … 3353 3325 3354 3326 /*clean up and return:*/ 3355 delete gauss_1;3356 delete gauss_2;3357 3327 return mass_flux; 3358 3328 } … … 3540 3510 3541 3511 /* Start looping on the number of vertices: */ 3542 GaussTria * gauss=new GaussTria();3512 GaussTria gauss; 3543 3513 for(int iv=0;iv<NUMVERTICES;iv++){ 3544 gauss ->GaussVertex(iv);3514 gauss.GaussVertex(iv); 3545 3515 3546 3516 /* Advection */ 3547 vx_input->GetInputValue(&v[0], gauss);3548 vy_input->GetInputValue(&v[1], gauss);3549 gr_input->GetInputValue(&groundedice, gauss);3517 vx_input->GetInputValue(&v[0],&gauss); 3518 vy_input->GetInputValue(&v[1],&gauss); 3519 gr_input->GetInputValue(&groundedice,&gauss); 3550 3520 3551 3521 /*Get calving speed*/ … … 3553 3523 case DefaultCalvingEnum: 3554 3524 case CalvingVonmisesEnum: 3555 lsf_slopex_input->GetInputValue(&dlsf[0], gauss);3556 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1], gauss);3557 calvingrate_input->GetInputValue(&calvingrate, gauss);3558 meltingrate_input->GetInputValue(&meltingrate, gauss);3525 lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); 3526 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); 3527 calvingrate_input->GetInputValue(&calvingrate,&gauss); 3528 meltingrate_input->GetInputValue(&meltingrate,&gauss); 3559 3529 if(groundedice<0) meltingrate = 0.; 3560 3530 … … 3574 3544 3575 3545 case CalvingLevermannEnum: 3576 calvingratex_input->GetInputValue(&c[0], gauss);3577 if(dim==2) calvingratey_input->GetInputValue(&c[1], gauss);3578 meltingrate_input->GetInputValue(&meltingrate, gauss);3546 calvingratex_input->GetInputValue(&c[0],&gauss); 3547 if(dim==2) calvingratey_input->GetInputValue(&c[1],&gauss); 3548 meltingrate_input->GetInputValue(&meltingrate,&gauss); 3579 3549 norm_calving=0.; 3580 3550 for(i=0;i<dim;i++) norm_calving+=pow(c[i],2); … … 3584 3554 3585 3555 case CalvingMinthicknessEnum: 3586 lsf_slopex_input->GetInputValue(&dlsf[0], gauss);3587 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1], gauss);3588 meltingrate_input->GetInputValue(&meltingrate, gauss);3556 lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); 3557 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); 3558 meltingrate_input->GetInputValue(&meltingrate,&gauss); 3589 3559 3590 3560 norm_dlsf=0.; … … 3605 3575 3606 3576 case CalvingHabEnum: 3607 lsf_slopex_input->GetInputValue(&dlsf[0], gauss);3608 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1], gauss);3609 meltingrate_input->GetInputValue(&meltingrate, gauss);3577 lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); 3578 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); 3579 meltingrate_input->GetInputValue(&meltingrate,&gauss); 3610 3580 3611 3581 norm_dlsf=0.; … … 3626 3596 3627 3597 case CalvingCrevasseDepthEnum: 3628 lsf_slopex_input->GetInputValue(&dlsf[0], gauss);3629 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1], gauss);3630 meltingrate_input->GetInputValue(&meltingrate, gauss);3598 lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); 3599 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); 3600 meltingrate_input->GetInputValue(&meltingrate,&gauss); 3631 3601 3632 3602 if(groundedice<0) meltingrate = 0.; … … 3650 3620 case CalvingDev2Enum: 3651 3621 { 3652 lsf_slopex_input->GetInputValue(&dlsf[0], gauss);3653 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1], gauss);3654 calvingrate_input->GetInputValue(&calvingrate, gauss);3655 meltingrate_input->GetInputValue(&meltingrate, gauss);3656 gr_input->GetInputValue(&groundedice, gauss);3622 lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); 3623 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); 3624 calvingrate_input->GetInputValue(&calvingrate,&gauss); 3625 meltingrate_input->GetInputValue(&meltingrate,&gauss); 3626 gr_input->GetInputValue(&groundedice,&gauss); 3657 3627 3658 3628 //idea: no retreat on ice above critical calving height "calvinghaf" . Limit using regularized Heaviside function. … … 3705 3675 this->AddInput(MovingFrontalVxEnum,&movingfrontvx[0],P1DGEnum); 3706 3676 this->AddInput(MovingFrontalVyEnum,&movingfrontvy[0],P1DGEnum); 3707 3708 /*Clean up and return*/3709 delete gauss;3710 3677 } 3711 3678 /*}}}*/ … … 3828 3795 int found = 0; 3829 3796 IssmDouble value; 3830 GaussTria *gauss = NULL;3797 GaussTria gauss; 3831 3798 3832 3799 /*First, serarch the input: */ … … 3840 3807 if(data){ 3841 3808 /*ok, we are good. retrieve value of input at vertex :*/ 3842 gauss =new GaussTria(); gauss->GaussVertex(i);3843 data->GetInputValue(&value, gauss);3809 gauss.GaussVertex(i); 3810 data->GetInputValue(&value,&gauss); 3844 3811 found=1; 3845 3812 break; … … 3847 3814 } 3848 3815 } 3849 3850 /*clean-up*/3851 delete gauss;3852 3816 3853 3817 if(found)*pvalue=value; … … 4006 3970 4007 3971 /*Loop over basal nodes and update their CS*/ 4008 GaussTria * gauss = new GaussTria();3972 GaussTria gauss; 4009 3973 for(int i=0;i<this->NumberofNodesVelocity();i++){ 4010 3974 4011 3975 if(vertexonbase[i]==1){ 4012 gauss ->GaussNode(this->VelocityInterpolation(),i);4013 slope_input->GetInputValue(&slope, gauss);4014 groundedicelevelset_input->GetInputValue(&groundedice, gauss);3976 gauss.GaussNode(this->VelocityInterpolation(),i); 3977 slope_input->GetInputValue(&slope,&gauss); 3978 groundedicelevelset_input->GetInputValue(&groundedice,&gauss); 4015 3979 IssmDouble theta = atan(slope); 4016 3980 … … 4033 3997 /*cleanup*/ 4034 3998 xDelete<IssmDouble>(vertexonbase); 4035 delete gauss;4036 3999 } 4037 4000 /*}}}*/ … … 4079 4042 4080 4043 /* Start looping on the number of vertices: */ 4081 GaussTria * gauss=new GaussTria();4044 GaussTria gauss; 4082 4045 for(int iv=0;iv<NUMVERTICES;iv++){ 4083 gauss ->GaussVertex(iv);4046 gauss.GaussVertex(iv); 4084 4047 4085 4048 /* Get variables */ 4086 bed_input->GetInputValue(&bed, gauss);4087 qsg_input->GetInputValue(&qsg, gauss);4088 TF_input->GetInputValue(&TF, gauss);4049 bed_input->GetInputValue(&bed,&gauss); 4050 qsg_input->GetInputValue(&qsg,&gauss); 4051 TF_input->GetInputValue(&TF,&gauss); 4089 4052 4090 4053 if(basin_icefront_area[reCast<int>(basinid[iv])-1]==0.) meltrates[iv]=0.; … … 4106 4069 /*Cleanup and return*/ 4107 4070 xDelete<IssmDouble>(basin_icefront_area); 4108 delete gauss;4109 4071 } 4110 4072 /*}}}*/ … … 4294 4256 4295 4257 IssmDouble epsilon[3]; 4296 GaussTria* gauss=NULL;4297 4258 IssmDouble vx,vy,vel; 4298 4259 IssmDouble strainxx; … … 4310 4271 4311 4272 /* Start looping on the number of vertices: */ 4312 gauss=new GaussTria();4273 GaussTria gauss; 4313 4274 for (int iv=0;iv<NUMVERTICES;iv++){ 4314 gauss ->GaussVertex(iv);4275 gauss.GaussVertex(iv); 4315 4276 4316 4277 /* Get the value we need*/ 4317 vx_input->GetInputValue(&vx, gauss);4318 vy_input->GetInputValue(&vy, gauss);4278 vx_input->GetInputValue(&vx,&gauss); 4279 vy_input->GetInputValue(&vy,&gauss); 4319 4280 vel=vx*vx+vy*vy; 4320 4281 4321 4282 /*Compute strain rate viscosity and pressure: */ 4322 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0], gauss,vx_input,vy_input);4283 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); 4323 4284 strainxx=epsilon[0]; 4324 4285 strainyy=epsilon[1]; … … 4331 4292 /*Add input*/ 4332 4293 this->AddInput(StrainRateparallelEnum,&strainparallel[0],P1DGEnum); 4333 4334 /*Clean up and return*/4335 delete gauss;4336 4294 } 4337 4295 /*}}}*/ 4338 4296 void Tria::StrainRateperpendicular(){/*{{{*/ 4339 4297 4340 GaussTria* gauss=NULL;4341 4298 IssmDouble epsilon[3]; 4342 4299 IssmDouble vx,vy,vel; … … 4355 4312 4356 4313 /* Start looping on the number of vertices: */ 4357 gauss=new GaussTria();4314 GaussTria gauss; 4358 4315 for (int iv=0;iv<NUMVERTICES;iv++){ 4359 gauss ->GaussVertex(iv);4316 gauss.GaussVertex(iv); 4360 4317 4361 4318 /* Get the value we need*/ 4362 vx_input->GetInputValue(&vx, gauss);4363 vy_input->GetInputValue(&vy, gauss);4319 vx_input->GetInputValue(&vx,&gauss); 4320 vy_input->GetInputValue(&vy,&gauss); 4364 4321 vel=vx*vx+vy*vy; 4365 4322 4366 4323 /*Compute strain rate viscosity and pressure: */ 4367 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0], gauss,vx_input,vy_input);4324 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],&gauss,vx_input,vy_input); 4368 4325 strainxx=epsilon[0]; 4369 4326 strainyy=epsilon[1]; … … 4376 4333 /*Add input*/ 4377 4334 this->AddInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1DGEnum); 4378 4379 /*Clean up and return*/4380 delete gauss;4381 4335 } 4382 4336 /*}}}*/ … … 4996 4950 Input* onbase = this->GetInput(MeshVertexonbaseEnum); _assert_(onbase); 4997 4951 4998 GaussTria * gauss=new GaussTria();4952 GaussTria gauss; 4999 4953 for(int iv=0;iv<this->NumberofNodes(this->element_type);iv++){ 5000 gauss ->GaussNode(this->element_type,iv);5001 onbase->GetInputValue(&isonbase, gauss);4954 gauss.GaussNode(this->element_type,iv); 4955 onbase->GetInputValue(&isonbase,&gauss); 5002 4956 if(isonbase==1.){ 5003 input->GetInputValue(&value, gauss);4957 input->GetInputValue(&value,&gauss); 5004 4958 this->nodes[iv]->ApplyConstraint(0,value); 5005 4959 } 5006 4960 } 5007 delete gauss;5008 4961 5009 4962 } … … 5020 4973 Input* onsurf = this->GetInput(MeshVertexonsurfaceEnum); _assert_(onsurf); 5021 4974 5022 GaussTria * gauss=new GaussTria();4975 GaussTria gauss; 5023 4976 for(int iv=0;iv<this->NumberofNodes(this->element_type);iv++){ 5024 gauss ->GaussNode(this->element_type,iv);5025 onsurf->GetInputValue(&isonsurface, gauss);4977 gauss.GaussNode(this->element_type,iv); 4978 onsurf->GetInputValue(&isonsurface,&gauss); 5026 4979 if(isonsurface==1.){ 5027 input->GetInputValue(&value, gauss);4980 input->GetInputValue(&value,&gauss); 5028 4981 this->nodes[iv]->ApplyConstraint(0,value); 5029 4982 } 5030 4983 } 5031 delete gauss;5032 4984 } 5033 4985 /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.
