Changeset 5638
- Timestamp:
- 08/31/10 14:26:29 (15 years ago)
- Location:
- issm/trunk/src/c/objects/Elements
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Tria.cpp
r5637 r5638 1007 1007 1008 1008 /* gaussian points: */ 1009 int num_gauss,ig; 1010 double* first_gauss_area_coord = NULL; 1011 double* second_gauss_area_coord = NULL; 1012 double* third_gauss_area_coord = NULL; 1013 double* gauss_weights = NULL; 1014 double gauss_weight; 1015 double gauss_l1l2l3[3]; 1009 int ig; 1010 GaussTria *gauss=NULL; 1016 1011 1017 1012 /* parameters: */ … … 1035 1030 double l1l2l3[3]; 1036 1031 double alpha2complement_list[numvertices]; //TO BE DELETED 1037 double gauss [numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED1032 double gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED 1038 1033 1039 1034 /* strain rate: */ … … 1075 1070 /*COMPUT alpha2_list (TO BE DELETED)*/ 1076 1071 for(i=0;i<numvertices;i++){ 1077 friction->GetAlphaComplement(&alpha2complement_list[i],&gauss[i][0],VxEnum,VyEnum); 1078 } 1079 1080 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 1081 GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4); 1072 friction->GetAlphaComplement(&alpha2complement_list[i],&gaussgrids[i][0],VxEnum,VyEnum); 1073 } 1082 1074 1083 1075 /*Retrieve all inputs we will be needing: */ … … 1089 1081 1090 1082 /* Start looping on the number of gaussian points: */ 1091 for (ig=0; ig<num_gauss; ig++){ 1092 /*Pick up the gaussian point: */ 1093 gauss_weight=*(gauss_weights+ig); 1094 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 1095 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 1096 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 1083 gauss=new GaussTria(4); 1084 for (ig=gauss->begin();ig<gauss->end();ig++){ 1085 1086 gauss->GaussPoint(ig); 1097 1087 1098 1088 /*Build alpha_complement_list: */ 1099 //if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss _l1l2l3,VxEnum,VyEnum); // TO BE UNCOMMENTED1089 //if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); // TO BE UNCOMMENTED 1100 1090 //else alpha_complement=0; 1101 TriaRef::GetParameterValue(&alpha_complement,&alpha2complement_list[0],gauss _l1l2l3); // TO BE DELETED1091 TriaRef::GetParameterValue(&alpha_complement,&alpha2complement_list[0],gauss); // TO BE DELETED 1102 1092 1103 1093 /*Recover alpha_complement and k: */ 1104 dragcoefficient_input->GetParameterValue(&drag, gauss _l1l2l3);1094 dragcoefficient_input->GetParameterValue(&drag, gauss); 1105 1095 1106 1096 /*recover lambda and mu: */ 1107 adjointx_input->GetParameterValue(&lambda, gauss _l1l2l3);1108 adjointy_input->GetParameterValue(&mu, gauss _l1l2l3);1097 adjointx_input->GetParameterValue(&lambda, gauss); 1098 adjointy_input->GetParameterValue(&mu, gauss); 1109 1099 1110 1100 /*recover vx and vy: */ 1111 inputs->GetParameterValue(&vx, gauss_l1l2l3,VxEnum);1112 inputs->GetParameterValue(&vy, gauss_l1l2l3,VyEnum);1101 vx_input->GetParameterValue(&vx,gauss); 1102 vy_input->GetParameterValue(&vy,gauss); 1113 1103 1114 1104 /* Get Jacobian determinant: */ 1115 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss _l1l2l3);1105 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 1116 1106 1117 1107 /* Get nodal functions value at gaussian point:*/ 1118 GetNodalFunctions(l1l2l3, gauss _l1l2l3);1108 GetNodalFunctions(l1l2l3, gauss); 1119 1109 1120 1110 /*Get nodal functions derivatives*/ 1121 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss _l1l2l3);1111 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss); 1122 1112 1123 1113 /*Get k derivative: dk/dx */ 1124 dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0], &gauss_l1l2l3[0]);1114 dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); 1125 1115 1126 1116 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ … … 1128 1118 1129 1119 //standard term dJ/dki 1130 grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss _weight*l1l2l3[i];1120 grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*l1l2l3[i]; 1131 1121 1132 1122 //noise dampening d/dki(1/2*(dk/dx)^2) 1133 grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss _weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);1123 grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); 1134 1124 } 1135 1125 … … 1141 1131 VecSetValues(gradient,numvertices,doflist1,(const double*)grade_g,ADD_VALUES); 1142 1132 1143 xfree((void**)&first_gauss_area_coord); 1144 xfree((void**)&second_gauss_area_coord); 1145 xfree((void**)&third_gauss_area_coord); 1146 xfree((void**)&gauss_weights); 1133 /*Clean up and return*/ 1134 delete gauss; 1147 1135 delete friction; 1148 1136 … … 3672 3660 3673 3661 /* gaussian points: */ 3674 int num_gauss,ig; 3675 double* first_gauss_area_coord = NULL; 3676 double* second_gauss_area_coord = NULL; 3677 double* third_gauss_area_coord = NULL; 3678 double* gauss_weights = NULL; 3679 double gauss_weight; 3680 double gauss_l1l2l3[3]; 3681 3662 int ig; 3663 GaussTria *gauss=NULL; 3682 3664 3683 3665 /* surface normal: */ … … 3704 3686 GetDofList(&doflist); 3705 3687 3706 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */3707 GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);3708 3709 3688 /*Build normal vector to the surface:*/ 3710 3711 3689 x4=xyz_list[0][0]; 3712 3690 y4=xyz_list[0][1]; … … 3737 3715 3738 3716 /* Start looping on the number of gaussian points: */ 3739 for (ig=0; ig<num_gauss; ig++){ 3740 /*Pick up the gaussian point: */ 3741 gauss_weight=*(gauss_weights+ig); 3742 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 3743 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 3744 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 3717 gauss=new GaussTria(2); 3718 for (ig=gauss->begin();ig<gauss->end();ig++){ 3719 3720 gauss->GaussPoint(ig); 3745 3721 3746 3722 /* Get Jacobian determinant: */ 3747 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss _l1l2l3);3723 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss); 3748 3724 3749 3725 //Get L matrix if viscous basal drag present: 3750 GetL(&L[0], &xyz_list[0][0], gauss _l1l2l3,NDOF1);3726 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); 3751 3727 3752 3728 /**********************Do not forget the sign**********************************/ 3753 DL_scalar=- gauss _weight*Jdet*nz;3729 DL_scalar=- gauss->weight*Jdet*nz; 3754 3730 /******************************************************************************/ 3755 3731 … … 3769 3745 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); 3770 3746 3771 xfree((void**)&first_gauss_area_coord); 3772 xfree((void**)&second_gauss_area_coord); 3773 xfree((void**)&third_gauss_area_coord); 3774 xfree((void**)&gauss_weights); 3747 /*Clean up and return*/ 3748 delete gauss; 3775 3749 xfree((void**)&doflist); 3776 3750 } … … 3794 3768 3795 3769 /* gaussian points: */ 3796 int num_area_gauss,ig; 3797 double* gauss_weights = NULL; 3798 double* first_gauss_area_coord = NULL; 3799 double* second_gauss_area_coord = NULL; 3800 double* third_gauss_area_coord = NULL; 3801 double gauss_weight; 3802 double gauss_coord[3]; 3770 int ig; 3771 GaussTria *gauss=NULL; 3803 3772 3804 3773 /*matrices: */ … … 3818 3787 GetDofList(&doflist); 3819 3788 3820 /* Get gaussian points and weights: */3821 GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);3822 3789 3823 3790 /* Start looping on the number of gauss (nodes on the bedrock) */ 3824 for (ig=0; ig<num_area_gauss; ig++){ 3825 gauss_weight=*(gauss_weights+ig); 3826 gauss_coord[0]=*(first_gauss_area_coord+ig); 3827 gauss_coord[1]=*(second_gauss_area_coord+ig); 3828 gauss_coord[2]=*(third_gauss_area_coord+ig); 3829 3830 //Get the Jacobian determinant 3831 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss_coord); 3791 gauss=new GaussTria(2); 3792 for (ig=gauss->begin();ig<gauss->end();ig++){ 3793 3794 gauss->GaussPoint(ig); 3795 3796 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss); 3832 3797 3833 3798 /*Get L matrix : */ 3834 GetL(&L[0], &xyz_list[0][0], gauss _coord,NDOF1);3799 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); 3835 3800 3836 3801 /*Calculate DL on gauss point */ 3837 D_scalar=latentheat/heatcapacity*gauss _weight*Jdet;3802 D_scalar=latentheat/heatcapacity*gauss->weight*Jdet; 3838 3803 3839 3804 /* Do the triple product tL*D*L: */ … … 3851 3816 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES); 3852 3817 3853 xfree((void**)&first_gauss_area_coord); 3854 xfree((void**)&second_gauss_area_coord); 3855 xfree((void**)&third_gauss_area_coord); 3856 xfree((void**)&gauss_weights); 3818 /*Clean up and return*/ 3819 delete gauss; 3857 3820 xfree((void**)&doflist); 3858 3821 } … … 3873 3836 3874 3837 /* gaussian points: */ 3875 int num_gauss,ig; 3876 double* first_gauss_area_coord = NULL; 3877 double* second_gauss_area_coord = NULL; 3878 double* third_gauss_area_coord = NULL; 3879 double* gauss_weights = NULL; 3880 double gauss_weight; 3881 double gauss_l1l2l3[3]; 3838 int ig; 3839 GaussTria *gauss=NULL; 3882 3840 3883 3841 /* matrices: */ … … 3933 3891 if(artdiff){ 3934 3892 //Get the Jacobian determinant 3935 gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; 3936 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); 3893 gauss=new GaussTria(); 3894 gauss->GaussCenter(); 3895 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); 3896 delete gauss; 3937 3897 3938 3898 //Build K matrix (artificial diffusivity matrix) … … 3944 3904 } 3945 3905 3946 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */3947 GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);3948 3949 3906 /* Start looping on the number of gaussian points: */ 3950 for (ig=0; ig<num_gauss; ig++){ 3951 3952 /*Pick up the gaussian point: */ 3953 gauss_weight=*(gauss_weights+ig); 3954 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 3955 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 3956 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 3907 gauss=new GaussTria(2); 3908 for (ig=gauss->begin();ig<gauss->end();ig++){ 3909 3910 gauss->GaussPoint(ig); 3957 3911 3958 3912 /* Get Jacobian determinant: */ 3959 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss _l1l2l3);3913 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); 3960 3914 3961 3915 /*Get L matrix: */ 3962 GetL(&L[0], &xyz_list[0][0], gauss _l1l2l3,numberofdofspernode);3963 3964 DL_scalar=gauss _weight*Jdettria;3916 GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); 3917 3918 DL_scalar=gauss->weight*Jdettria; 3965 3919 3966 3920 /* Do the triple product tL*D*L: */ … … 3971 3925 3972 3926 /*Get B and B prime matrix: */ 3973 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss _l1l2l3);3974 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss _l1l2l3);3927 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss); 3928 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); 3975 3929 3976 3930 //Get vx, vy and their derivatives at gauss point 3977 vxaverage_input->GetParameterValue(&vx, &gauss_l1l2l3[0]);3978 vyaverage_input->GetParameterValue(&vy, &gauss_l1l2l3[0]);3979 3980 vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0], &gauss_l1l2l3[0]);3981 vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0], &gauss_l1l2l3[0]);3931 vxaverage_input->GetParameterValue(&vx,gauss); 3932 vyaverage_input->GetParameterValue(&vy,gauss); 3933 3934 vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); 3935 vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); 3982 3936 3983 3937 dvxdx=dvx[0]; 3984 3938 dvydy=dvy[1]; 3985 3939 3986 DL_scalar=dt*gauss _weight*Jdettria;3940 DL_scalar=dt*gauss->weight*Jdettria; 3987 3941 3988 3942 //Create DL and DLprime matrix … … 4032 3986 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); 4033 3987 4034 xfree((void**)&first_gauss_area_coord); 4035 xfree((void**)&second_gauss_area_coord); 4036 xfree((void**)&third_gauss_area_coord); 4037 xfree((void**)&gauss_weights); 3988 /*Clean up and return*/ 3989 delete gauss; 4038 3990 xfree((void**)&doflist); 4039 3991 } … … 4054 4006 4055 4007 /* gaussian points: */ 4056 int num_gauss,ig; 4057 double* first_gauss_area_coord = NULL; 4058 double* second_gauss_area_coord = NULL; 4059 double* third_gauss_area_coord = NULL; 4060 double* gauss_weights = NULL; 4061 double gauss_weight; 4062 double gauss_l1l2l3[3]; 4008 int ig; 4009 GaussTria *gauss=NULL; 4063 4010 4064 4011 /* matrices: */ … … 4092 4039 GetDofList(&doflist); 4093 4040 4094 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */4095 GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);4096 4097 4041 /*Retrieve all inputs we will be needed*/ 4098 4042 if(dim==2){ … … 4106 4050 4107 4051 /* Start looping on the number of gaussian points: */ 4108 for (ig=0; ig<num_gauss; ig++){ 4109 4110 /*Pick up the gaussian point: */ 4111 gauss_weight=*(gauss_weights+ig); 4112 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 4113 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 4114 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 4052 gauss=new GaussTria(2); 4053 for (ig=gauss->begin();ig<gauss->end();ig++){ 4054 4055 gauss->GaussPoint(ig); 4115 4056 4116 4057 /* Get Jacobian determinant: */ 4117 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss _l1l2l3);4058 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); 4118 4059 4119 4060 /*Get L matrix: */ 4120 GetL(&L[0], &xyz_list[0][0], gauss _l1l2l3,numberofdofspernode);4121 4122 DL_scalar=gauss _weight*Jdettria;4061 GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); 4062 4063 DL_scalar=gauss->weight*Jdettria; 4123 4064 4124 4065 /* Do the triple product tL*D*L: */ … … 4130 4071 /*Get B and B prime matrix: */ 4131 4072 /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/ 4132 GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss _l1l2l3);4133 GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss _l1l2l3);4073 GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); 4074 GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss); 4134 4075 4135 4076 //Get vx, vy and their derivatives at gauss point 4136 vxaverage_input->GetParameterValue(&vx, &gauss_l1l2l3[0]);4137 vyaverage_input->GetParameterValue(&vy, &gauss_l1l2l3[0]);4138 4139 DL_scalar=-dt*gauss _weight*Jdettria;4077 vxaverage_input->GetParameterValue(&vx,gauss); 4078 vyaverage_input->GetParameterValue(&vy,gauss); 4079 4080 DL_scalar=-dt*gauss->weight*Jdettria; 4140 4081 4141 4082 DLprime[0][0]=DL_scalar*vx; … … 4157 4098 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); 4158 4099 4159 xfree((void**)&first_gauss_area_coord); 4160 xfree((void**)&second_gauss_area_coord); 4161 xfree((void**)&third_gauss_area_coord); 4162 xfree((void**)&gauss_weights); 4100 /*Clean up and return*/ 4101 delete gauss; 4163 4102 xfree((void**)&doflist); 4164 4103 } … … 4231 4170 double heatcapacity; 4232 4171 4233 int num_gauss,ig; 4234 double* first_gauss_area_coord = NULL; 4235 double* second_gauss_area_coord = NULL; 4236 double* third_gauss_area_coord = NULL; 4237 double* gauss_weights = NULL; 4238 double gauss_weight; 4239 double gauss_coord[3]; 4172 int ig; 4173 GaussTria *gauss=NULL; 4240 4174 4241 4175 /*matrices: */ … … 4264 4198 heatcapacity=matpar->GetHeatCapacity(); 4265 4199 4266 GaussLegendreTria (&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);4267 4268 4200 /* Start looping on the number of gauss (nodes on the bedrock) */ 4269 for (ig=0; ig<num_gauss; ig++){ 4270 gauss_weight=*(gauss_weights+ig); 4271 gauss_coord[0]=*(first_gauss_area_coord+ig); 4272 gauss_coord[1]=*(second_gauss_area_coord+ig); 4273 gauss_coord[2]=*(third_gauss_area_coord+ig); 4201 gauss=new GaussTria(2); 4202 for (ig=gauss->begin();ig<gauss->end();ig++){ 4203 4204 gauss->GaussPoint(ig); 4274 4205 4275 4206 //Get the Jacobian determinant 4276 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss _coord);4207 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss); 4277 4208 4278 4209 /*Get nodal functions values: */ 4279 GetNodalFunctions(&l1l2l3[0], gauss _coord);4210 GetNodalFunctions(&l1l2l3[0], gauss); 4280 4211 4281 4212 /*Calculate DL on gauss point */ 4282 D_scalar=gauss _weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice);4213 D_scalar=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice); 4283 4214 if(dt){ 4284 4215 D_scalar=dt*D_scalar; … … 4299 4230 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES); 4300 4231 4301 xfree((void**)&first_gauss_area_coord); 4302 xfree((void**)&second_gauss_area_coord); 4303 xfree((void**)&third_gauss_area_coord); 4304 xfree((void**)&gauss_weights); 4232 4233 /*Clean up and return*/ 4234 delete gauss; 4305 4235 xfree((void**)&doflist); 4306 4236 } … … 4321 4251 4322 4252 /* gaussian points: */ 4323 int num_gauss,ig; 4324 double* first_gauss_area_coord = NULL; 4325 double* second_gauss_area_coord = NULL; 4326 double* third_gauss_area_coord = NULL; 4327 double* gauss_weights = NULL; 4328 double gauss_weight; 4329 double gauss_l1l2l3[3]; 4253 int ig; 4254 GaussTria* gauss=NULL; 4330 4255 4331 4256 /* matrix */ … … 4347 4272 GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices); 4348 4273 GetDofList(&doflist); 4349 4350 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */4351 GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);4352 4274 4353 4275 /*retrieve inputs :*/ … … 4357 4279 4358 4280 /* Start looping on the number of gaussian points: */ 4359 for (ig=0; ig<num_gauss; ig++){ 4360 /*Pick up the gaussian point: */ 4361 gauss_weight=*(gauss_weights+ig); 4362 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 4363 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 4364 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 4281 gauss=new GaussTria(2); 4282 for(ig=gauss->begin();ig<gauss->end();ig++){ 4283 4284 gauss->GaussPoint(ig); 4365 4285 4366 4286 /* Get Jacobian determinant: */ 4367 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss _l1l2l3);4287 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); 4368 4288 4369 4289 /*Get L matrix: */ 4370 GetL(&L[0], &xyz_list[0][0], gauss _l1l2l3,numberofdofspernode);4290 GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); 4371 4291 4372 4292 /* Get accumulation, melting at gauss point */ 4373 accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]);4374 melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]);4375 dhdt_input->GetParameterValue(&dhdt_g, &gauss_l1l2l3[0]);4293 accumulation_input->GetParameterValue(&accumulation_g,gauss); 4294 melting_input->GetParameterValue(&melting_g,gauss); 4295 dhdt_input->GetParameterValue(&dhdt_g,gauss); 4376 4296 4377 4297 /* Add value into pe_g: */ 4378 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss _weight*(accumulation_g-melting_g-dhdt_g)*L[i];4298 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g-dhdt_g)*L[i]; 4379 4299 4380 4300 } // for (ig=0; ig<num_gauss; ig++) … … 4383 4303 VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); 4384 4304 4385 xfree((void**)&first_gauss_area_coord); 4386 xfree((void**)&second_gauss_area_coord); 4387 xfree((void**)&third_gauss_area_coord); 4388 xfree((void**)&gauss_weights); 4305 /*Clean up and return*/ 4306 delete gauss; 4389 4307 xfree((void**)&doflist); 4390 4308 } … … 4405 4323 4406 4324 /* gaussian points: */ 4407 int num_gauss,ig; 4408 double* first_gauss_area_coord = NULL; 4409 double* second_gauss_area_coord = NULL; 4410 double* third_gauss_area_coord = NULL; 4411 double* gauss_weights = NULL; 4412 double gauss_weight; 4413 double gauss_l1l2l3[3]; 4325 int ig; 4326 GaussTria* gauss=NULL; 4414 4327 4415 4328 /* matrix */ … … 4432 4345 GetDofList(&doflist); 4433 4346 4434 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */4435 GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);4436 4437 4347 /*Retrieve all inputs we will be needing: */ 4438 4348 accumulation_input=inputs->GetInput(AccumulationRateEnum); … … 4441 4351 4442 4352 /* Start looping on the number of gaussian points: */ 4443 for (ig=0; ig<num_gauss; ig++){ 4444 /*Pick up the gaussian point: */ 4445 gauss_weight=*(gauss_weights+ig); 4446 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 4447 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 4448 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 4353 gauss=new GaussTria(2); 4354 for(ig=gauss->begin();ig<gauss->end();ig++){ 4355 4356 gauss->GaussPoint(ig); 4449 4357 4450 4358 /* Get Jacobian determinant: */ 4451 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss _l1l2l3);4359 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); 4452 4360 4453 4361 /*Get L matrix: */ 4454 GetL(&L[0], &xyz_list[0][0], gauss _l1l2l3,numberofdofspernode);4362 GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); 4455 4363 4456 4364 /* Get accumulation, melting and thickness at gauss point */ 4457 accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]);4458 melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]);4459 dhdt_input->GetParameterValue(&dhdt_g, &gauss_l1l2l3[0]);4365 accumulation_input->GetParameterValue(&accumulation_g,gauss); 4366 melting_input->GetParameterValue(&melting_g,gauss); 4367 dhdt_input->GetParameterValue(&dhdt_g,gauss); 4460 4368 4461 4369 /* Add value into pe_g: */ 4462 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss _weight*(accumulation_g-melting_g-dhdt_g)*L[i];4370 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g-dhdt_g)*L[i]; 4463 4371 4464 4372 } // for (ig=0; ig<num_gauss; ig++) … … 4467 4375 VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); 4468 4376 4469 xfree((void**)&first_gauss_area_coord); 4470 xfree((void**)&second_gauss_area_coord); 4471 xfree((void**)&third_gauss_area_coord); 4472 xfree((void**)&gauss_weights); 4377 /*Clean up and return*/ 4378 delete gauss; 4473 4379 xfree((void**)&doflist); 4474 4380 } … … 4489 4395 4490 4396 /* gaussian points: */ 4491 int num_gauss,ig; 4492 double* first_gauss_area_coord = NULL; 4493 double* second_gauss_area_coord = NULL; 4494 double* third_gauss_area_coord = NULL; 4495 double* gauss_weights = NULL; 4496 double gauss_weight; 4497 double gauss_l1l2l3[3]; 4397 int ig; 4398 GaussTria* gauss=NULL; 4498 4399 4499 4400 /* matrix */ … … 4514 4415 GetDofList(&doflist); 4515 4416 4516 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */4517 GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);4518 4519 4417 /*Retrieve all inputs we will be needing: */ 4520 4418 accumulation_input=inputs->GetInput(AccumulationRateEnum); … … 4522 4420 4523 4421 /* Start looping on the number of gaussian points: */ 4524 for (ig=0; ig<num_gauss; ig++){ 4525 /*Pick up the gaussian point: */ 4526 gauss_weight=*(gauss_weights+ig); 4527 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 4528 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 4529 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 4422 gauss=new GaussTria(2); 4423 for(ig=gauss->begin();ig<gauss->end();ig++){ 4424 4425 gauss->GaussPoint(ig); 4530 4426 4531 4427 /* Get Jacobian determinant: */ 4532 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss _l1l2l3);4428 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); 4533 4429 4534 4430 /*Get L matrix: */ 4535 GetL(&L[0], &xyz_list[0][0], gauss _l1l2l3,numberofdofspernode);4431 GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); 4536 4432 4537 4433 /* Get accumulation, melting at gauss point */ 4538 accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]);4539 melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]);4434 accumulation_input->GetParameterValue(&accumulation_g,gauss); 4435 melting_input->GetParameterValue(&melting_g,gauss); 4540 4436 4541 4437 /* Add value into pe_g: */ 4542 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss _weight*(accumulation_g-melting_g)*L[i];4438 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g)*L[i]; 4543 4439 4544 4440 } // for (ig=0; ig<num_gauss; ig++) … … 4547 4443 VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); 4548 4444 4549 /*Free ressources:*/ 4550 xfree((void**)&first_gauss_area_coord); 4551 xfree((void**)&second_gauss_area_coord); 4552 xfree((void**)&third_gauss_area_coord); 4553 xfree((void**)&gauss_weights); 4445 /*Clean up and return*/ 4446 delete gauss; 4554 4447 xfree((void**)&doflist); 4555 4448 } … … 4568 4461 4569 4462 /* gaussian points: */ 4570 int num_gauss,ig; 4571 double* first_gauss_area_coord = NULL; 4572 double* second_gauss_area_coord = NULL; 4573 double* third_gauss_area_coord = NULL; 4574 double* gauss_weights = NULL; 4575 double gauss_weight; 4576 double gauss_l1l2l3[3]; 4463 int ig; 4464 GaussTria* gauss=NULL; 4577 4465 4578 4466 /* Jacobian: */ … … 4605 4493 GetDofList(&doflist); 4606 4494 4607 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */4608 GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);4609 4610 4495 /*Retrieve all inputs we will be needing: */ 4611 4496 melting_input=inputs->GetInput(MeltingRateEnum); … … 4616 4501 /*For icesheets: */ 4617 4502 /* Start looping on the number of gaussian points: */ 4618 for (ig=0; ig<num_gauss; ig++){ 4619 4620 /*Pick up the gaussian point: */ 4621 gauss_weight=*(gauss_weights+ig); 4622 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 4623 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 4624 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 4503 gauss=new GaussTria(2); 4504 for(ig=gauss->begin();ig<gauss->end();ig++){ 4505 4506 gauss->GaussPoint(ig); 4625 4507 4626 4508 /*Get melting at gaussian point: */ 4627 melting_input->GetParameterValue(&meltingvalue, &gauss_l1l2l3[0]);4509 melting_input->GetParameterValue(&meltingvalue, gauss); 4628 4510 4629 4511 /*Get velocity at gaussian point: */ 4630 vx_input->GetParameterValue(&vx, &gauss_l1l2l3[0]);4631 vy_input->GetParameterValue(&vy, &gauss_l1l2l3[0]);4512 vx_input->GetParameterValue(&vx, gauss); 4513 vy_input->GetParameterValue(&vy, gauss); 4632 4514 4633 4515 /*Get bed slope: */ 4634 bed_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0], &gauss_l1l2l3[0]);4516 bed_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); 4635 4517 dbdx=slope[0]; 4636 4518 dbdy=slope[1]; 4637 4519 4638 4520 /* Get Jacobian determinant: */ 4639 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss _l1l2l3);4521 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss); 4640 4522 4641 4523 //Get L matrix if viscous basal drag present: 4642 GetL(&L[0], &xyz_list[0][0], gauss _l1l2l3,NDOF1);4524 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); 4643 4525 4644 4526 /*Build gaussian vector: */ 4645 4527 for(i=0;i<numvertices;i++){ 4646 pe_g_gaussian[i]=-Jdet*gauss _weight*(vx*dbdx+vy*dbdy-meltingvalue)*L[i];4528 pe_g_gaussian[i]=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-meltingvalue)*L[i]; 4647 4529 } 4648 4530 … … 4655 4537 VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); 4656 4538 4657 /*Free ressources:*/ 4658 xfree((void**)&first_gauss_area_coord); 4659 xfree((void**)&second_gauss_area_coord); 4660 xfree((void**)&third_gauss_area_coord); 4661 xfree((void**)&gauss_weights); 4539 /*Clean up and return*/ 4540 delete gauss; 4662 4541 xfree((void**)&doflist); 4663 4542 } -
issm/trunk/src/c/objects/Elements/TriaRef.cpp
r5636 r5638 487 487 } 488 488 /*}}}*/ 489 /*FUNCTION TriaRef::GetJacobianDeterminant3d {{{1*/ 490 void TriaRef::GetJacobianDeterminant3d(double* Jdet, double* xyz_list,GaussTria* gauss){ 491 /*The Jacobian determinant is constant over the element, discard the gaussian points. 492 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ 493 494 double x1,x2,x3,y1,y2,y3,z1,z2,z3; 495 496 x1=*(xyz_list+3*0+0); 497 y1=*(xyz_list+3*0+1); 498 z1=*(xyz_list+3*0+2); 499 x2=*(xyz_list+3*1+0); 500 y2=*(xyz_list+3*1+1); 501 z2=*(xyz_list+3*1+2); 502 x3=*(xyz_list+3*2+0); 503 y3=*(xyz_list+3*2+1); 504 z3=*(xyz_list+3*2+2); 505 506 *Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5); 507 if(*Jdet<0) ISSMERROR("negative jacobian determinant!"); 508 509 } 510 /*}}}*/ 489 511 /*FUNCTION TriaRef::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss) {{{1*/ 490 512 void TriaRef::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss){ -
issm/trunk/src/c/objects/Elements/TriaRef.h
r5636 r5638 40 40 void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss); 41 41 void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss); 42 void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,GaussTria* gauss); 42 43 void GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss); 43 44 void GetJacobianInvert(double* Jinv, double* xyz_list,GaussTria* gauss);
Note:
See TracChangeset
for help on using the changeset viewer.