Changeset 5647
- Timestamp:
- 09/01/10 10:45:00 (15 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Penta.cpp
r5640 r5647 2713 2713 2714 2714 /* 3d gaussian points: */ 2715 int num_gauss,ig; 2716 double* first_gauss_area_coord = NULL; 2717 double* second_gauss_area_coord = NULL; 2718 double* third_gauss_area_coord = NULL; 2719 double* fourth_gauss_vert_coord = NULL; 2720 double* area_gauss_weights = NULL; 2721 double* vert_gauss_weights = NULL; 2722 int ig1,ig2; 2723 double gauss_weight1,gauss_weight2; 2724 double gauss_coord[4]; 2725 int order_area_gauss; 2726 int num_vert_gauss; 2727 int num_area_gauss; 2728 double gauss_weight; 2729 2730 /* 2d gaussian point: */ 2731 int num_gauss2d; 2732 double* first_gauss_area_coord2d = NULL; 2733 double* second_gauss_area_coord2d = NULL; 2734 double* third_gauss_area_coord2d = NULL; 2735 double* gauss_weights2d=NULL; 2736 double gauss_l1l2l3[3]; 2715 int ig; 2716 GaussPenta *gauss=NULL; 2737 2717 2738 2718 /* material data: */ … … 2802 2782 GetDofList(&doflist,PattynApproximationEnum); 2803 2783 2804 /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore2805 get tria gaussian points as well as segment gaussian points. For tria gaussian2806 points, order of integration is 2, because we need to integrate the product tB*D*B'2807 which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian2808 points, same deal, which yields 3 gaussian points.*/2809 2810 order_area_gauss=5;2811 num_vert_gauss=5;2812 2813 gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);2814 2815 2784 /*Retrieve all inputs we will be needing: */ 2816 2785 vx_input=inputs->GetInput(VxEnum); … … 2820 2789 2821 2790 /* Start looping on the number of gaussian points: */ 2822 for (ig1=0; ig1<num_area_gauss; ig1++){ 2823 for (ig2=0; ig2<num_vert_gauss; ig2++){ 2824 2825 /*Pick up the gaussian point: */ 2826 gauss_weight1=*(area_gauss_weights+ig1); 2827 gauss_weight2=*(vert_gauss_weights+ig2); 2828 gauss_weight=gauss_weight1*gauss_weight2; 2829 2830 2831 gauss_coord[0]=*(first_gauss_area_coord+ig1); 2832 gauss_coord[1]=*(second_gauss_area_coord+ig1); 2833 gauss_coord[2]=*(third_gauss_area_coord+ig1); 2834 gauss_coord[3]=*(fourth_gauss_vert_coord+ig2); 2835 2836 2837 /*Get strain rate from velocity: */ 2838 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input); 2839 this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input); 2840 2841 /*Get viscosity: */ 2842 matice->GetViscosity3d(&viscosity, &epsilon[0]); 2843 matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); 2844 2845 /*Get B and Bprime matrices: */ 2846 GetBPattyn(&B[0][0], &xyz_list[0][0], gauss_coord); 2847 GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss_coord); 2848 2849 /* Get Jacobian determinant: */ 2850 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); 2851 2852 /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant 2853 onto this scalar matrix, so that we win some computational time: */ 2854 2855 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); 2856 D_scalar=2*newviscosity*gauss_weight*Jdet; 2857 for (i=0;i<5;i++){ 2858 D[i][i]=D_scalar; 2859 } 2860 2861 /* Do the triple product tB*D*Bprime: */ 2862 TripleMultiply( &B[0][0],5,numdof,1, 2863 &D[0][0],5,5,0, 2864 &Bprime[0][0],5,numdof,0, 2865 &Ke_gg_gaussian[0][0],0); 2866 2867 /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ 2868 for( i=0; i<numdof; i++){ 2869 for(j=0;j<numdof;j++){ 2870 Ke_gg[i][j]+=Ke_gg_gaussian[i][j]; 2871 } 2872 } 2873 } //for (ig2=0; ig2<num_vert_gauss; ig2++) 2874 } //for (ig1=0; ig1<num_area_gauss; ig1++) 2791 gauss=new GaussPenta(5,5); 2792 for (ig=gauss->begin();ig<gauss->end();ig++){ 2793 2794 gauss->GaussPoint(ig); 2795 2796 /*Get strain rate from velocity: */ 2797 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 2798 this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); 2799 2800 /*Get viscosity: */ 2801 matice->GetViscosity3d(&viscosity, &epsilon[0]); 2802 matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); 2803 2804 /*Get B and Bprime matrices: */ 2805 GetBPattyn(&B[0][0], &xyz_list[0][0], gauss); 2806 GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss); 2807 2808 /* Get Jacobian determinant: */ 2809 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); 2810 2811 /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant 2812 onto this scalar matrix, so that we win some computational time: */ 2813 2814 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); 2815 D_scalar=2*newviscosity*gauss->weight*Jdet; 2816 for (i=0;i<5;i++) D[i][i]=D_scalar; 2817 2818 /* Do the triple product tB*D*Bprime: */ 2819 TripleMultiply( &B[0][0],5,numdof,1, 2820 &D[0][0],5,5,0, 2821 &Bprime[0][0],5,numdof,0, 2822 &Ke_gg_gaussian[0][0],0); 2823 2824 /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ 2825 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j]; 2826 } 2875 2827 2876 2828 … … 2880 2832 //Deal with 2d friction at the bedrock interface 2881 2833 if((onbed && !shelf)){ 2882 2883 2834 /*Build a tria element using the 3 grids of the base of the penta. Then use 2884 2835 * the tria functionality to build a friction stiffness matrix on these 3 … … 2890 2841 } 2891 2842 2892 xfree((void**)&first_gauss_area_coord); 2893 xfree((void**)&second_gauss_area_coord); 2894 xfree((void**)&third_gauss_area_coord); 2895 xfree((void**)&fourth_gauss_vert_coord); 2896 xfree((void**)&area_gauss_weights); 2897 xfree((void**)&vert_gauss_weights); 2898 xfree((void**)&first_gauss_area_coord2d); 2899 xfree((void**)&second_gauss_area_coord2d); 2900 xfree((void**)&third_gauss_area_coord2d); 2901 xfree((void**)&gauss_weights2d); 2843 /*Clean up and return*/ 2844 delete gauss; 2902 2845 xfree((void**)&doflist); 2903 2846 } … … 5208 5151 /*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/ 5209 5152 void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input){ 5153 /*Compute the 3d Blatter/PattynStrain Rate (5 components): 5154 * 5155 * epsilon=[exx eyy exy exz eyz] 5156 * 5157 * with exz=1/2 du/dz 5158 * eyz=1/2 dv/dz 5159 * 5160 * the contribution of vz is neglected 5161 */ 5162 5163 int i; 5164 5165 double epsilonvx[5]; 5166 double epsilonvy[5]; 5167 5168 /*Check that both inputs have been found*/ 5169 if (!vx_input || !vy_input){ 5170 ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input); 5171 } 5172 5173 /*Get strain rate assuming that epsilon has been allocated*/ 5174 vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss); 5175 vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss); 5176 5177 /*Sum all contributions*/ 5178 for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; 5179 5180 } 5181 /*}}}*/ 5182 /*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/ 5183 void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input){ 5210 5184 /*Compute the 3d Blatter/PattynStrain Rate (5 components): 5211 5185 * -
issm/trunk/src/c/objects/Elements/Penta.h
r5635 r5647 168 168 void GetSolutionFromInputsThermal(Vec solutiong); 169 169 void GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input); 170 void GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input); 170 171 void GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input); 171 172 Penta* GetUpperElement(void); -
issm/trunk/src/c/objects/Elements/PentaRef.cpp
r5203 r5647 981 981 } 982 982 /*}}}*/ 983 984 /*FUNCTION PentaRef::GetBMacAyealPattyn {{{1*/ 985 void PentaRef::GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss){ 986 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 987 * For grid i, Bi can be expressed in the actual coordinate system 988 * by: 989 * Bi=[ dh/dx 0 ] 990 * [ 0 dh/dy ] 991 * [ 1/2*dh/dy 1/2*dh/dx ] 992 * where h is the interpolation function for grid i. 993 * 994 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids) 995 */ 996 997 int i; 998 const int numgrids=6; 999 const int NDOF3=3; 1000 const int NDOF2=2; 1001 1002 double dh1dh6[NDOF3][numgrids]; 1003 1004 /*Get dh1dh6 in actual coordinate system: */ 1005 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss); 1006 1007 /*Build B: */ 1008 for (i=0;i<numgrids;i++){ 1009 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i]; 1010 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0; 1011 1012 *(B+NDOF2*numgrids*1+NDOF2*i)=0.0; 1013 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i]; 1014 1015 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i]; 1016 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i]; 1017 1018 } 1019 1020 } 1021 /*}}}*/ 1022 /*FUNCTION PentaRef::GetBPattyn {{{1*/ 1023 void PentaRef::GetBPattyn(double* B, double* xyz_list, GaussPenta* gauss){ 1024 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 1025 * For grid i, Bi can be expressed in the actual coordinate system 1026 * by: 1027 * Bi=[ dh/dx 0 ] 1028 * [ 0 dh/dy ] 1029 * [ 1/2*dh/dy 1/2*dh/dx ] 1030 * [ 1/2*dh/dz 0 ] 1031 * [ 0 1/2*dh/dz ] 1032 * where h is the interpolation function for grid i. 1033 * 1034 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids) 1035 */ 1036 1037 int i; 1038 const int numgrids=6; 1039 const int NDOF3=3; 1040 const int NDOF2=2; 1041 1042 double dh1dh6[NDOF3][numgrids]; 1043 1044 /*Get dh1dh6 in actual coordinate system: */ 1045 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss); 1046 1047 /*Build B: */ 1048 for (i=0;i<numgrids;i++){ 1049 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i]; 1050 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0; 1051 1052 *(B+NDOF2*numgrids*1+NDOF2*i)=0.0; 1053 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i]; 1054 1055 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i]; 1056 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i]; 1057 1058 *(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh6[2][i]; 1059 *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0; 1060 1061 *(B+NDOF2*numgrids*4+NDOF2*i)=0.0; 1062 *(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh6[2][i]; 1063 } 1064 1065 } 1066 /*}}}*/ 1067 /*FUNCTION PentaRef::GetBprimePattyn {{{1*/ 1068 void PentaRef::GetBprimePattyn(double* B, double* xyz_list, GaussPenta* gauss_coord){ 1069 /*Compute B prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 1070 * For grid i, Bi can be expressed in the actual coordinate system 1071 * by: 1072 * Bi=[ 2*dh/dx dh/dy ] 1073 * [ dh/dx 2*dh/dy ] 1074 * [ dh/dy dh/dx ] 1075 * [ dh/dz 0 ] 1076 * [ 0 dh/dz ] 1077 * where h is the interpolation function for grid i. 1078 * 1079 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids) 1080 */ 1081 1082 int i; 1083 const int NDOF3=3; 1084 const int NDOF2=2; 1085 const int numgrids=6; 1086 1087 double dh1dh6[NDOF3][numgrids]; 1088 1089 /*Get dh1dh6 in actual coordinate system: */ 1090 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss_coord); 1091 1092 /*Build BPrime: */ 1093 for (i=0;i<numgrids;i++){ 1094 *(B+NDOF2*numgrids*0+NDOF2*i)=2.0*dh1dh6[0][i]; 1095 *(B+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh6[1][i]; 1096 1097 *(B+NDOF2*numgrids*1+NDOF2*i)=dh1dh6[0][i]; 1098 *(B+NDOF2*numgrids*1+NDOF2*i+1)=2.0*dh1dh6[1][i]; 1099 1100 *(B+NDOF2*numgrids*2+NDOF2*i)=dh1dh6[1][i]; 1101 *(B+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh6[0][i]; 1102 1103 *(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh6[2][i]; 1104 *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0; 1105 1106 *(B+NDOF2*numgrids*4+NDOF2*i)=0.0; 1107 *(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh6[2][i]; 1108 } 1109 } 1110 /*}}}*/ 1111 /*FUNCTION PentaRef::GetBStokes {{{1*/ 1112 void PentaRef::GetBStokes(double* B, double* xyz_list, GaussPenta* gauss){ 1113 1114 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID. 1115 * For grid i, Bi can be expressed in the actual coordinate system 1116 * by: Bi=[ dh/dx 0 0 0 ] 1117 * [ 0 dh/dy 0 0 ] 1118 * [ 0 0 dh/dy 0 ] 1119 * [ 1/2*dh/dy 1/2*dh/dx 0 0 ] 1120 * [ 1/2*dh/dz 0 1/2*dh/dx 0 ] 1121 * [ 0 1/2*dh/dz 1/2*dh/dy 0 ] 1122 * [ 0 0 0 h ] 1123 * [ dh/dx dh/dy dh/dz 0 ] 1124 * where h is the interpolation function for grid i. 1125 * Same thing for Bb except the last column that does not exist. 1126 */ 1127 1128 int i; 1129 const int calculationdof=3; 1130 const int numgrids=6; 1131 int DOFPERGRID=4; 1132 1133 double dh1dh7[calculationdof][numgrids+1]; 1134 double l1l6[numgrids]; 1135 1136 1137 /*Get dh1dh7 in actual coordinate system: */ 1138 GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss); 1139 GetNodalFunctionsP1(l1l6, gauss); 1140 1141 /*Build B: */ 1142 for (i=0;i<numgrids+1;i++){ 1143 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B[0][DOFPERGRID*i]=dh1dh6[0][i]; 1144 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0; 1145 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0; 1146 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0; 1147 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i]; 1148 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0; 1149 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0; 1150 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0; 1151 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i]; 1152 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7[1][i]; 1153 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7[0][i]; 1154 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0; 1155 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7[2][i]; 1156 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0; 1157 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7[0][i]; 1158 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0; 1159 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7[2][i]; 1160 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7[1][i]; 1161 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=0; 1162 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=0; 1163 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=0; 1164 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7[0][i]; 1165 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7[1][i]; 1166 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7[2][i]; 1167 } 1168 1169 for (i=0;i<numgrids;i++){ //last column not for the bubble function 1170 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0; 1171 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0; 1172 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0; 1173 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0; 1174 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0; 1175 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0; 1176 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=l1l6[i]; 1177 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=0; 1178 } 1179 1180 } 1181 /*}}}*/ 1182 /*FUNCTION PentaRef::GetBprimeStokes {{{1*/ 1183 void PentaRef::GetBprimeStokes(double* B_prime, double* xyz_list, GaussPenta* gauss){ 1184 /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. 1185 * For grid i, Bi' can be expressed in the actual coordinate system 1186 * by: 1187 * Bi'=[ dh/dx 0 0 0] 1188 * [ 0 dh/dy 0 0] 1189 * [ 0 0 dh/dz 0] 1190 * [ dh/dy dh/dx 0 0] 1191 * [ dh/dz 0 dh/dx 0] 1192 * [ 0 dh/dz dh/dy 0] 1193 * [ dh/dx dh/dy dh/dz 0] 1194 * [ 0 0 0 h] 1195 * where h is the interpolation function for grid i. 1196 * 1197 * Same thing for the bubble fonction except that there is no fourth column 1198 */ 1199 1200 int i; 1201 const int calculationdof=3; 1202 const int numgrids=6; 1203 int DOFPERGRID=4; 1204 1205 double dh1dh7[calculationdof][numgrids+1]; 1206 double l1l6[numgrids]; 1207 1208 /*Get dh1dh7 in actual coordinate system: */ 1209 GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss); 1210 1211 GetNodalFunctionsP1(l1l6, gauss); 1212 1213 /*B_primeuild B_prime: */ 1214 for (i=0;i<numgrids+1;i++){ 1215 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B_prime[0][DOFPERGRID*i]=dh1dh6[0][i]; 1216 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0; 1217 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0; 1218 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0; 1219 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i]; 1220 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0; 1221 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0; 1222 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0; 1223 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i]; 1224 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=dh1dh7[1][i]; 1225 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=dh1dh7[0][i]; 1226 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0; 1227 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=dh1dh7[2][i]; 1228 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0; 1229 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=dh1dh7[0][i]; 1230 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0; 1231 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=dh1dh7[2][i]; 1232 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=dh1dh7[1][i]; 1233 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=dh1dh7[0][i]; 1234 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=dh1dh7[1][i]; 1235 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=dh1dh7[2][i]; 1236 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=0; 1237 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=0; 1238 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=0; 1239 } 1240 1241 for (i=0;i<numgrids;i++){ //last column not for the bubble function 1242 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0; 1243 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0; 1244 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0; 1245 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0; 1246 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0; 1247 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0; 1248 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=0; 1249 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=l1l6[i]; 1250 } 1251 1252 } 1253 /*}}}*/ 1254 /*FUNCTION PentaRef::GetBArtdiff {{{1*/ 1255 void PentaRef::GetBArtdiff(double* B_artdiff, double* xyz_list, GaussPenta* gauss){ 1256 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 1257 * For grid i, Bi' can be expressed in the actual coordinate system 1258 * by: 1259 * Bi_artdiff=[ dh/dx ] 1260 * [ dh/dy ] 1261 * where h is the interpolation function for grid i. 1262 * 1263 * We assume B has been allocated already, of size: 2x(DOFPERGRID*numgrids) 1264 */ 1265 1266 int i; 1267 const int calculationdof=3; 1268 const int numgrids=6; 1269 int DOFPERGRID=1; 1270 1271 /*Same thing in the actual coordinate system: */ 1272 double dh1dh6[calculationdof][numgrids]; 1273 1274 /*Get dh1dh6 in actual coordinates system : */ 1275 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss); 1276 1277 /*Build B': */ 1278 for (i=0;i<numgrids;i++){ 1279 *(B_artdiff+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i]; 1280 *(B_artdiff+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i]; 1281 } 1282 } 1283 /*}}}*/ 1284 /*FUNCTION PentaRef::GetBAdvec{{{1*/ 1285 void PentaRef::GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss){ 1286 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 1287 * For grid i, Bi' can be expressed in the actual coordinate system 1288 * by: 1289 * Bi_advec =[ h ] 1290 * [ h ] 1291 * [ h ] 1292 * where h is the interpolation function for grid i. 1293 * 1294 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids) 1295 */ 1296 1297 int i; 1298 int calculationdof=3; 1299 int numgrids=6; 1300 int DOFPERGRID=1; 1301 1302 /*Same thing in the actual coordinate system: */ 1303 double l1l6[6]; 1304 1305 /*Get dh1dh2dh3 in actual coordinates system : */ 1306 GetNodalFunctionsP1(l1l6, gauss); 1307 1308 /*Build B': */ 1309 for (i=0;i<numgrids;i++){ 1310 *(B_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=l1l6[i]; 1311 *(B_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=l1l6[i]; 1312 *(B_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=l1l6[i]; 1313 } 1314 } 1315 /*}}}*/ 1316 /*FUNCTION PentaRef::GetBConduct{{{1*/ 1317 void PentaRef::GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss){ 1318 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 1319 * For grid i, Bi' can be expressed in the actual coordinate system 1320 * by: 1321 * Bi_conduct=[ dh/dx ] 1322 * [ dh/dy ] 1323 * [ dh/dz ] 1324 * where h is the interpolation function for grid i. 1325 * 1326 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids) 1327 */ 1328 1329 int i; 1330 const int calculationdof=3; 1331 const int numgrids=6; 1332 int DOFPERGRID=1; 1333 1334 /*Same thing in the actual coordinate system: */ 1335 double dh1dh6[calculationdof][numgrids]; 1336 1337 /*Get dh1dh2dh3 in actual coordinates system : */ 1338 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss); 1339 1340 /*Build B': */ 1341 for (i=0;i<numgrids;i++){ 1342 *(B_conduct+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i]; 1343 *(B_conduct+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i]; 1344 *(B_conduct+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i]; 1345 } 1346 } 1347 /*}}}*/ 1348 /*FUNCTION PentaRef::GetBVert{{{1*/ 1349 void PentaRef::GetBVert(double* B, double* xyz_list, GaussPenta* gauss){ 1350 /* Compute B matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz]; 1351 where hi is the interpolation function for grid i.*/ 1352 1353 int i; 1354 const int NDOF3=3; 1355 const int numgrids=6; 1356 double dh1dh6[NDOF3][numgrids]; 1357 1358 /*Get dh1dh6 in actual coordinate system: */ 1359 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss); 1360 1361 /*Build B: */ 1362 for (i=0;i<numgrids;i++){ 1363 B[i]=dh1dh6[2][i]; 1364 } 1365 1366 } 1367 /*}}}*/ 1368 /*FUNCTION PentaRef::GetBprimeAdvec{{{1*/ 1369 void PentaRef::GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss){ 1370 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 1371 * For grid i, Bi' can be expressed in the actual coordinate system 1372 * by: 1373 * Biprime_advec=[ dh/dx ] 1374 * [ dh/dy ] 1375 * [ dh/dz ] 1376 * where h is the interpolation function for grid i. 1377 * 1378 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids) 1379 */ 1380 1381 int i; 1382 const int calculationdof=3; 1383 const int numgrids=6; 1384 int DOFPERGRID=1; 1385 1386 /*Same thing in the actual coordinate system: */ 1387 double dh1dh6[calculationdof][numgrids]; 1388 1389 /*Get dh1dh2dh3 in actual coordinates system : */ 1390 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss); 1391 1392 /*Build B': */ 1393 for (i=0;i<numgrids;i++){ 1394 *(Bprime_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i]; 1395 *(Bprime_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i]; 1396 *(Bprime_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i]; 1397 } 1398 } 1399 /*}}}*/ 1400 /*FUNCTION PentaRef::GetBprimeVert{{{1*/ 1401 void PentaRef::GetBprimeVert(double* B, double* xyz_list, GaussPenta* gauss){ 1402 /* Compute Bprime matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for grid i*/ 1403 1404 GetNodalFunctionsP1(B, gauss); 1405 1406 } 1407 /*}}}*/ 1408 /*FUNCTION PentaRef::GetLprimeStokes {{{1*/ 1409 void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_tria, GaussPenta* gauss){ 1410 1411 /* 1412 * Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 1413 * For grid i, Lpi can be expressed in the actual coordinate system 1414 * by: 1415 * Lpi=[ h 0 0 0] 1416 * [ 0 h 0 0] 1417 * [ h 0 0 0] 1418 * [ 0 h 0 0] 1419 * [ 0 0 h 0] 1420 * [ 0 0 h 0] 1421 * [ 0 0 dh/dz 0] 1422 * [ 0 0 dh/dz 0] 1423 * [ 0 0 dh/dz 0] 1424 * [dh/dz 0 dh/dx 0] 1425 * [ 0 dh/dz dh/dy 0] 1426 * [ 0 0 0 h] 1427 * [ 0 0 0 h] 1428 * [ 0 0 0 h] 1429 * where h is the interpolation function for grid i. 1430 */ 1431 int i; 1432 const int numgrids2d=3; 1433 int num_dof=4; 1434 1435 double l1l2l3[numgrids2d]; 1436 double dh1dh6[3][6]; 1437 1438 /*Get l1l2l3 in actual coordinate system: */ 1439 l1l2l3[0]=gauss_tria[0]; 1440 l1l2l3[1]=gauss_tria[1]; 1441 l1l2l3[2]=gauss_tria[2]; 1442 1443 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss); 1444 1445 /*Build LprimeStokes: */ 1446 for (i=0;i<3;i++){ 1447 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i]; 1448 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+1)=0; 1449 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+2)=0; 1450 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+3)=0; 1451 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i)=0; 1452 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i]; 1453 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+2)=0; 1454 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+3)=0; 1455 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i)=l1l2l3[i]; 1456 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+1)=0; 1457 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+2)=0; 1458 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+3)=0; 1459 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i)=0; 1460 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+1)=l1l2l3[i]; 1461 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+2)=0; 1462 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+3)=0; 1463 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i)=0; 1464 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+1)=0; 1465 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+2)=l1l2l3[i]; 1466 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+3)=0; 1467 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i)=0; 1468 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+1)=0; 1469 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+2)=l1l2l3[i]; 1470 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+3)=0; 1471 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i)=0; 1472 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+1)=0; 1473 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+2)=dh1dh6[2][i]; 1474 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+3)=0; 1475 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i)=0; 1476 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+1)=0; 1477 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+2)=dh1dh6[2][i]; 1478 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+3)=0; 1479 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i)=0; 1480 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+1)=0; 1481 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+2)=dh1dh6[2][i]; 1482 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+3)=0; 1483 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i)=dh1dh6[2][i]; 1484 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+1)=0; 1485 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+2)=dh1dh6[0][i]; 1486 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+3)=0; 1487 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i)=0; 1488 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+1)=dh1dh6[2][i]; 1489 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+2)=dh1dh6[1][i]; 1490 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+3)=0; 1491 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i)=0; 1492 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+1)=0; 1493 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+2)=0; 1494 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+3)=l1l2l3[i]; 1495 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i)=0; 1496 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+1)=0; 1497 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+2)=0; 1498 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+3)=l1l2l3[i]; 1499 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i)=0; 1500 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+1)=0; 1501 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+2)=0; 1502 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+3)=l1l2l3[i]; 1503 1504 } 1505 } 1506 /*}}}*/ 1507 /*FUNCTION PentaRef::GetJacobian {{{1*/ 1508 void PentaRef::GetJacobian(double* J, double* xyz_list,GaussPenta* gauss){ 1509 1510 const int NDOF3=3; 1511 int i,j; 1512 1513 /*The Jacobian is constant over the element, discard the gaussian points. 1514 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ 1515 1516 double A1,A2,A3; //area coordinates 1517 double xi,eta,zi; //parametric coordinates 1518 1519 double x1,x2,x3,x4,x5,x6; 1520 double y1,y2,y3,y4,y5,y6; 1521 double z1,z2,z3,z4,z5,z6; 1522 1523 /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */ 1524 A1=gauss->coord1; 1525 A2=gauss->coord2; 1526 A3=gauss->coord3; 1527 1528 xi=A2-A1; 1529 eta=SQRT3*A3; 1530 zi=gauss->coord4; 1531 1532 x1=*(xyz_list+3*0+0); 1533 x2=*(xyz_list+3*1+0); 1534 x3=*(xyz_list+3*2+0); 1535 x4=*(xyz_list+3*3+0); 1536 x5=*(xyz_list+3*4+0); 1537 x6=*(xyz_list+3*5+0); 1538 1539 y1=*(xyz_list+3*0+1); 1540 y2=*(xyz_list+3*1+1); 1541 y3=*(xyz_list+3*2+1); 1542 y4=*(xyz_list+3*3+1); 1543 y5=*(xyz_list+3*4+1); 1544 y6=*(xyz_list+3*5+1); 1545 1546 z1=*(xyz_list+3*0+2); 1547 z2=*(xyz_list+3*1+2); 1548 z3=*(xyz_list+3*2+2); 1549 z4=*(xyz_list+3*3+2); 1550 z5=*(xyz_list+3*4+2); 1551 z6=*(xyz_list+3*5+2); 1552 1553 *(J+NDOF3*0+0)=0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5); 1554 *(J+NDOF3*1+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6); 1555 *(J+NDOF3*2+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4); 1556 1557 *(J+NDOF3*0+1)=0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5); 1558 *(J+NDOF3*1+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6); 1559 *(J+NDOF3*2+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2); 1560 1561 *(J+NDOF3*0+2)=0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5); 1562 *(J+NDOF3*1+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6); 1563 *(J+NDOF3*2+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4); 1564 1565 } 1566 /*}}}*/ 1567 /*FUNCTION PentaRef::GetJacobianDeterminant {{{1*/ 1568 void PentaRef::GetJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss){ 1569 /*On a penta, Jacobian varies according to coordinates. We need to get the Jacobian, and take 1570 * the determinant of it: */ 1571 double J[3][3]; 1572 1573 /*Get Jacobian*/ 1574 GetJacobian(&J[0][0],xyz_list,gauss); 1575 1576 /*Get Determinant*/ 1577 Matrix3x3Determinant(Jdet,&J[0][0]); 1578 if(*Jdet<0) ISSMERROR("negative jacobian determinant!"); 1579 1580 } 1581 /*}}}*/ 1582 /*FUNCTION PentaRef::GetJacobianInvert {{{1*/ 1583 void PentaRef::GetJacobianInvert(double* Jinv, double* xyz_list,GaussPenta* gauss){ 1584 1585 /*Jacobian*/ 1586 double J[3][3]; 1587 1588 /*Call Jacobian routine to get the jacobian:*/ 1589 GetJacobian(&J[0][0], xyz_list, gauss); 1590 1591 /*Invert Jacobian matrix: */ 1592 Matrix3x3Invert(Jinv,&J[0][0]); 1593 } 1594 /*}}}*/ 1595 /*FUNCTION PentaRef::GetNodalFunctionsMINI{{{1*/ 1596 void PentaRef::GetNodalFunctionsMINI(double* l1l7, GaussPenta* gauss){ 1597 /*This routine returns the values of the nodal functions at the gaussian point.*/ 1598 1599 l1l7[0]=gauss->coord1*(1.0-gauss->coord4)/2.0; 1600 l1l7[1]=gauss->coord2*(1.0-gauss->coord4)/2.0; 1601 l1l7[2]=gauss->coord3*(1.0-gauss->coord4)/2.0; 1602 l1l7[3]=gauss->coord1*(1.0+gauss->coord4)/2.0; 1603 l1l7[4]=gauss->coord2*(1.0+gauss->coord4)/2.0; 1604 l1l7[5]=gauss->coord3*(1.0+gauss->coord4)/2.0; 1605 l1l7[6]=27*gauss->coord1*gauss->coord2*gauss->coord3*(1.0+gauss->coord4)*(1.0-gauss->coord4); 1606 1607 } 1608 /*}}}*/ 1609 /*FUNCTION PentaRef::GetNodalFunctionsMINIDerivatives{{{1*/ 1610 void PentaRef::GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, GaussPenta* gauss){ 1611 1612 /*This routine returns the values of the nodal functions derivatives (with respect to the 1613 * actual coordinate system): */ 1614 1615 int i; 1616 const int numgrids = 7; 1617 double dh1dh7_ref[3][numgrids]; 1618 double Jinv[3][3]; 1619 1620 /*Get derivative values with respect to parametric coordinate system: */ 1621 GetNodalFunctionsMINIDerivativesReference(&dh1dh7_ref[0][0], gauss); 1622 1623 /*Get Jacobian invert: */ 1624 GetJacobianInvert(&Jinv[0][0], xyz_list, gauss); 1625 1626 /*Build dh1dh6: 1627 * 1628 * [dhi/dx]= Jinv'*[dhi/dr] 1629 * [dhi/dy] [dhi/ds] 1630 * [dhi/dz] [dhi/dzeta] 1631 */ 1632 1633 for (i=0;i<numgrids;i++){ 1634 *(dh1dh7+numgrids*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i]; 1635 *(dh1dh7+numgrids*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i]; 1636 *(dh1dh7+numgrids*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i]; 1637 } 1638 1639 } 1640 /*}}}*/ 1641 /*FUNCTION PentaRef::GetNodalFunctionsMINIDerivativesReference{{{1*/ 1642 void PentaRef::GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,GaussPenta* gauss){ 1643 1644 /*This routine returns the values of the nodal functions derivatives (with respect to the 1645 * natural coordinate system) at the gaussian point. */ 1646 1647 int numgrids=7; //six plus bubble grids 1648 1649 double r=gauss->coord2-gauss->coord1; 1650 double s=-3.0/SQRT3*(gauss->coord1+gauss->coord2-2.0/3.0); 1651 double zeta=gauss->coord4; 1652 1653 /*First nodal function: */ 1654 *(dl1dl7+numgrids*0+0)=-0.5*(1.0-zeta)/2.0; 1655 *(dl1dl7+numgrids*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0; 1656 *(dl1dl7+numgrids*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD); 1657 1658 /*Second nodal function: */ 1659 *(dl1dl7+numgrids*0+1)=0.5*(1.0-zeta)/2.0; 1660 *(dl1dl7+numgrids*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0; 1661 *(dl1dl7+numgrids*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD); 1662 1663 /*Third nodal function: */ 1664 *(dl1dl7+numgrids*0+2)=0; 1665 *(dl1dl7+numgrids*1+2)=SQRT3/3.0*(1.0-zeta)/2.0; 1666 *(dl1dl7+numgrids*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD); 1667 1668 /*Fourth nodal function: */ 1669 *(dl1dl7+numgrids*0+3)=-0.5*(1.0+zeta)/2.0; 1670 *(dl1dl7+numgrids*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0; 1671 *(dl1dl7+numgrids*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD); 1672 1673 /*Fith nodal function: */ 1674 *(dl1dl7+numgrids*0+4)=0.5*(1.0+zeta)/2.0; 1675 *(dl1dl7+numgrids*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0; 1676 *(dl1dl7+numgrids*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD); 1677 1678 /*Sixth nodal function: */ 1679 *(dl1dl7+numgrids*0+5)=0; 1680 *(dl1dl7+numgrids*1+5)=SQRT3/3.0*(1.0+zeta)/2.0; 1681 *(dl1dl7+numgrids*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD); 1682 1683 /*Seventh nodal function: */ 1684 *(dl1dl7+numgrids*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0); 1685 *(dl1dl7+numgrids*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0)); 1686 *(dl1dl7+numgrids*2+6)=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta); 1687 1688 } 1689 /*}}}*/ 1690 /*FUNCTION PentaRef::GetNodalFunctionsP1 {{{1*/ 1691 void PentaRef::GetNodalFunctionsP1(double* l1l6, GaussPenta* gauss){ 1692 /*This routine returns the values of the nodal functions at the gaussian point.*/ 1693 1694 l1l6[0]=gauss->coord1*(1-gauss->coord4)/2.0; 1695 l1l6[1]=gauss->coord2*(1-gauss->coord4)/2.0; 1696 l1l6[2]=gauss->coord3*(1-gauss->coord4)/2.0; 1697 l1l6[3]=gauss->coord1*(1+gauss->coord4)/2.0; 1698 l1l6[4]=gauss->coord2*(1+gauss->coord4)/2.0; 1699 l1l6[5]=gauss->coord3*(1+gauss->coord4)/2.0; 1700 1701 } 1702 /*}}}*/ 1703 /*FUNCTION PentaRef::GetNodalFunctionsP1Derivatives {{{1*/ 1704 void PentaRef::GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, GaussPenta* gauss){ 1705 1706 /*This routine returns the values of the nodal functions derivatives (with respect to the 1707 * actual coordinate system): */ 1708 int i; 1709 const int NDOF3 = 3; 1710 const int numgrids = 6; 1711 double dh1dh6_ref[NDOF3][numgrids]; 1712 double Jinv[NDOF3][NDOF3]; 1713 1714 /*Get derivative values with respect to parametric coordinate system: */ 1715 GetNodalFunctionsP1DerivativesReference(&dh1dh6_ref[0][0], gauss); 1716 1717 /*Get Jacobian invert: */ 1718 GetJacobianInvert(&Jinv[0][0], xyz_list, gauss); 1719 1720 /*Build dh1dh3: 1721 * 1722 * [dhi/dx]= Jinv*[dhi/dr] 1723 * [dhi/dy] [dhi/ds] 1724 * [dhi/dz] [dhi/dn] 1725 */ 1726 1727 for (i=0;i<numgrids;i++){ 1728 *(dh1dh6+numgrids*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i]; 1729 *(dh1dh6+numgrids*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i]; 1730 *(dh1dh6+numgrids*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i]; 1731 } 1732 1733 } 1734 /*}}}*/ 1735 /*FUNCTION PentaRef::GetNodalFunctionsP1DerivativesReference {{{1*/ 1736 void PentaRef::GetNodalFunctionsP1DerivativesReference(double* dl1dl6,GaussPenta* gauss){ 1737 1738 /*This routine returns the values of the nodal functions derivatives (with respect to the 1739 * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */ 1740 1741 const int numgrids=6; 1742 double A1,A2,A3,z; 1743 1744 A1=gauss->coord1; ISSMASSERT(A1>=0 && A1<=1);//first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3); 1745 A2=gauss->coord2; ISSMASSERT(A2>=0 && A2<=1);//second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3); 1746 A3=gauss->coord3; ISSMASSERT(A3>=0 && A3<=1);//third area coordinate value In term of xi and eta: A3=y/SQRT3; 1747 z =gauss->coord4; ISSMASSERT(z>=-1 && z<=1);//fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2 1748 1749 /*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/ 1750 *(dl1dl6+numgrids*0+0)=-0.5*(1.0-z)/2.0; 1751 *(dl1dl6+numgrids*1+0)=-0.5/SQRT3*(1.0-z)/2.0; 1752 *(dl1dl6+numgrids*2+0)=-0.5*A1; 1753 1754 /*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/ 1755 *(dl1dl6+numgrids*0+1)=0.5*(1.0-z)/2.0; 1756 *(dl1dl6+numgrids*1+1)=-0.5/SQRT3*(1.0-z)/2.0; 1757 *(dl1dl6+numgrids*2+1)=-0.5*A2; 1758 1759 /*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/ 1760 *(dl1dl6+numgrids*0+2)=0.0; 1761 *(dl1dl6+numgrids*1+2)=1.0/SQRT3*(1.0-z)/2.0; 1762 *(dl1dl6+numgrids*2+2)=-0.5*A3; 1763 1764 /*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/ 1765 *(dl1dl6+numgrids*0+3)=-0.5*(1.0+z)/2.0; 1766 *(dl1dl6+numgrids*1+3)=-0.5/SQRT3*(1.0+z)/2.0; 1767 *(dl1dl6+numgrids*2+3)=0.5*A1; 1768 1769 /*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/ 1770 *(dl1dl6+numgrids*0+4)=0.5*(1.0+z)/2.0; 1771 *(dl1dl6+numgrids*1+4)=-0.5/SQRT3*(1.0+z)/2.0; 1772 *(dl1dl6+numgrids*2+4)=0.5*A2; 1773 1774 /*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/ 1775 *(dl1dl6+numgrids*0+5)=0.0; 1776 *(dl1dl6+numgrids*1+5)=1.0/SQRT3*(1.0+z)/2.0; 1777 *(dl1dl6+numgrids*2+5)=0.5*A3; 1778 } 1779 /*}}}*/ 1780 /*FUNCTION PentaRef::GetParameterValue{{{1*/ 1781 void PentaRef::GetParameterValue(double* pvalue,double* plist,GaussPenta* gauss){ 1782 /*P1 interpolation on Gauss point*/ 1783 1784 /*intermediary*/ 1785 double l1l6[6]; 1786 1787 /*nodal functions: */ 1788 GetNodalFunctionsP1(&l1l6[0],gauss); 1789 1790 /*Assign output pointers:*/ 1791 *pvalue=l1l6[0]*plist[0]+l1l6[1]*plist[1]+l1l6[2]*plist[2]+l1l6[3]*plist[3]+l1l6[4]*plist[4]+l1l6[5]*plist[5]; 1792 1793 } 1794 /*}}}*/ 1795 /*FUNCTION PentaRef::GetParameterDerivativeValue{{{1*/ 1796 void PentaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, GaussPenta* gauss){ 1797 /*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_coord: 1798 * dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx; 1799 * dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy; 1800 * dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz; 1801 * 1802 * p is a vector of size 3x1 already allocated. 1803 */ 1804 double dh1dh6[3][6]; 1805 1806 /*Get nodal funnctions derivatives in actual coordinate system: */ 1807 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss); 1808 1809 /*Assign output*/ 1810 p[0]=plist[0]*dh1dh6[0][0]+plist[1]*dh1dh6[0][1]+plist[2]*dh1dh6[0][2]+plist[3]*dh1dh6[0][3]+plist[4]*dh1dh6[0][4]+plist[5]*dh1dh6[0][5]; 1811 p[1]=plist[0]*dh1dh6[1][0]+plist[1]*dh1dh6[1][1]+plist[2]*dh1dh6[1][2]+plist[3]*dh1dh6[1][3]+plist[4]*dh1dh6[1][4]+plist[5]*dh1dh6[1][5]; 1812 p[2]=plist[0]*dh1dh6[2][0]+plist[1]*dh1dh6[2][1]+plist[2]*dh1dh6[2][2]+plist[3]*dh1dh6[2][3]+plist[4]*dh1dh6[2][4]+plist[5]*dh1dh6[2][5]; 1813 1814 } 1815 /*}}}*/ -
issm/trunk/src/c/objects/Elements/PentaRef.h
r5629 r5647 50 50 void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");}; 51 51 52 void GetNodalFunctionsP1(double* l1l6, GaussPenta* gauss); 53 void GetNodalFunctionsMINI(double* l1l7, GaussPenta* gauss); 54 void GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, GaussPenta* gauss); 55 void GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, GaussPenta* gauss); 56 void GetNodalFunctionsP1DerivativesReference(double* dl1dl6,GaussPenta* gauss); 57 void GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,GaussPenta* gauss); 58 void GetJacobian(double* J, double* xyz_list,GaussPenta* gauss); 59 void GetJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss); 60 void GetJacobianInvert(double* Jinv, double* xyz_list,GaussPenta* gauss); 61 void GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss); 62 void GetBPattyn(double* B, double* xyz_list, GaussPenta* gauss); 63 void GetBprimePattyn(double* B, double* xyz_list, GaussPenta* gauss); 64 void GetBStokes(double* B, double* xyz_list, GaussPenta* gauss); 65 void GetBprimeStokes(double* B_prime, double* xyz_list, GaussPenta* gauss); 66 void GetBprimeVert(double* B, double* xyz_list, GaussPenta* gauss); 67 void GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss); 68 void GetBArtdiff(double* B_artdiff, double* xyz_list, GaussPenta* gauss); 69 void GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss); 70 void GetBVert(double* B, double* xyz_list, GaussPenta* gauss); 71 void GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss); 72 //void GetLStokes(double* LStokes, double* gauss_tria); 73 void GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_tria, GaussPenta* gauss); 74 void GetParameterValue(double* pvalue,double* plist,GaussPenta* gauss); 75 //void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");}; 76 void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussPenta* gauss); 77 //void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");}; 78 52 79 }; 53 80 #endif -
issm/trunk/src/c/objects/Inputs/BoolInput.cpp
r5629 r5647 172 172 void BoolInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");} 173 173 /*}}}*/ 174 /*FUNCTION BoolInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/ 175 void BoolInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");} 176 /*}}}*/ 174 177 /*FUNCTION BoolInput::GetParameterValues{{{1*/ 175 178 void BoolInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){ISSMERROR(" not supported yet!");} … … 180 183 /*FUNCTION BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/ 181 184 void BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR(" not supported yet!");} 185 /*}}}*/ 186 /*FUNCTION BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){{{1*/ 187 void BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){ISSMERROR(" not supported yet!");} 182 188 /*}}}*/ 183 189 /*FUNCTION BoolInput::ChangeEnum{{{1*/ -
issm/trunk/src/c/objects/Inputs/BoolInput.h
r5631 r5647 49 49 void GetParameterValue(double* pvalue,double* gauss); 50 50 void GetParameterValue(double* pvalue,GaussTria* gauss); 51 void GetParameterValue(double* pvalue,GaussPenta* gauss); 51 52 void GetParameterValues(double* values,double* gauss_pointers, int numgauss); 52 53 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss); 53 54 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); 55 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss); 54 56 void GetParameterAverage(double* pvalue){ISSMERROR("not implemented yet");}; 55 57 void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");}; … … 62 64 void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; 63 65 void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; 64 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};65 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};66 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};67 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};68 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};66 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 67 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 68 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 69 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 70 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 69 71 void ChangeEnum(int newenumtype); 70 72 void SquareMin(double* psquaremin, bool process_units,Parameters* parameters); -
issm/trunk/src/c/objects/Inputs/DoubleInput.cpp
r5629 r5647 185 185 void DoubleInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");} 186 186 /*}}}*/ 187 /*FUNCTION DoubleInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/ 188 void DoubleInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");} 189 /*}}}*/ 187 190 /*FUNCTION DoubleInput::GetParameterValues{{{1*/ 188 191 void DoubleInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){ISSMERROR(" not supported yet!");} … … 193 196 /*FUNCTION DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/ 194 197 void DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR(" not supported yet!");} 198 /*}}}*/ 199 /*FUNCTION DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){{{1*/ 200 void DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){ISSMERROR(" not supported yet!");} 195 201 /*}}}*/ 196 202 /*FUNCTION DoubleInput::ChangeEnum{{{1*/ -
issm/trunk/src/c/objects/Inputs/DoubleInput.h
r5631 r5647 48 48 void GetParameterValue(double* pvalue,double* gauss); 49 49 void GetParameterValue(double* pvalue,GaussTria* gauss); 50 void GetParameterValue(double* pvalue,GaussPenta* gauss); 50 51 void GetParameterValues(double* values,double* gauss_pointers, int numgauss); 51 52 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss); 52 53 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); 54 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss); 53 55 void GetParameterAverage(double* pvalue); 54 56 void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");}; … … 61 63 void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; 62 64 void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; 63 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};64 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};65 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};66 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};67 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};65 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 66 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 67 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 68 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 69 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 68 70 void ChangeEnum(int newenumtype); 69 71 void SquareMin(double* psquaremin, bool process_units,Parameters* parameters); -
issm/trunk/src/c/objects/Inputs/Input.h
r5631 r5647 27 27 virtual void GetParameterValue(double* pvalue,double* gauss)=0; 28 28 virtual void GetParameterValue(double* pvalue,GaussTria* gauss)=0; 29 virtual void GetParameterValue(double* pvalue,GaussPenta* gauss)=0; 29 30 virtual void GetParameterValues(double* values,double* gauss_pointers, int numgauss)=0; 30 31 virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss)=0; 31 32 virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss)=0; 33 virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss)=0; 32 34 virtual void GetParameterAverage(double* pvalue)=0; 33 35 virtual void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss)=0; … … 40 42 virtual void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss)=0; 41 43 virtual void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss)=0; 42 virtual void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, Gauss Tria* gauss)=0;43 virtual void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, Gauss Tria* gauss)=0;44 virtual void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, Gauss Tria* gauss)=0;45 virtual void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, Gauss Tria* gauss)=0;46 virtual void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, Gauss Tria* gauss)=0;44 virtual void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss)=0; 45 virtual void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss)=0; 46 virtual void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss)=0; 47 virtual void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss)=0; 48 virtual void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss)=0; 47 49 virtual void ChangeEnum(int newenumtype)=0; 48 50 -
issm/trunk/src/c/objects/Inputs/IntInput.cpp
r5629 r5647 173 173 void IntInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");} 174 174 /*}}}*/ 175 /*FUNCTION IntInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/ 176 void IntInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");} 177 /*}}}*/ 175 178 /*FUNCTION IntInput::GetParameterValues{{{1*/ 176 179 void IntInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){ISSMERROR(" not supported yet!");} … … 181 184 /*FUNCTION IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/ 182 185 void IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR(" not supported yet!");} 186 /*}}}*/ 187 /*FUNCTION IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){{{1*/ 188 void IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){ISSMERROR(" not supported yet!");} 183 189 /*}}}*/ 184 190 /*FUNCTION IntInput::ChangeEnum{{{1*/ -
issm/trunk/src/c/objects/Inputs/IntInput.h
r5631 r5647 49 49 void GetParameterValue(double* pvalue,double* gauss); 50 50 void GetParameterValue(double* pvalue,GaussTria* gauss); 51 void GetParameterValue(double* pvalue,GaussPenta* gauss); 51 52 void GetParameterValues(double* values,double* gauss_pointers, int numgauss); 52 53 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss); 53 54 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); 55 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss); 54 56 void GetParameterAverage(double* pvalue){ISSMERROR("not implemented yet");}; 55 57 void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");}; … … 62 64 void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; 63 65 void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; 64 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};65 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};66 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};67 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};68 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};66 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 67 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 68 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 69 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 70 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 69 71 void ChangeEnum(int newenumtype); 70 72 void SquareMin(double* psquaremin, bool process_units,Parameters* parameters); -
issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp
r5629 r5647 184 184 } 185 185 /*}}}*/ 186 /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,Gauss Tria* gauss){{{1*/187 void PentaVertexInput::GetParameterValue(double* pvalue,Gauss Tria* gauss){186 /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/ 187 void PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ 188 188 189 189 /*Call PentaRef function*/ … … 216 216 } 217 217 /*}}}*/ 218 /*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, Gauss Tria* gauss){{{1*/219 void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, Gauss Tria* gauss){218 /*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussPenta* gauss){{{1*/ 219 void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussPenta* gauss){ 220 220 221 221 /*Call PentaRef function*/ … … 386 386 /*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn{{{1*/ 387 387 void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss){ 388 389 int i; 390 const int numgrids=6; 391 const int NDOF2=2; 392 double B[5][NDOF2*numgrids]; 393 double velocity[numgrids][NDOF2]; 394 395 /*Get B matrix: */ 396 GetBPattyn(&B[0][0], xyz_list, gauss); 397 398 /*Here, we are computing the strain rate of (0,vy)*/ 399 for(i=0;i<numgrids;i++){ 400 velocity[i][0]=0.0; 401 velocity[i][1]=this->values[i]; 402 } 403 404 /*Multiply B by velocity, to get strain rate: */ 405 MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0, 406 &velocity[0][0],NDOF2*numgrids,1,0, 407 epsilonvy,0); 408 409 } 410 /*}}}*/ 411 /*FUNCTION PentaVertexInput::GetVxStrainRate3d{{{1*/ 412 void PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ 413 int i,j; 414 415 const int numgrids=6; 416 const int DOFVELOCITY=3; 417 double B[8][27]; 418 double B_reduced[6][DOFVELOCITY*numgrids]; 419 double velocity[numgrids][DOFVELOCITY]; 420 421 /*Get B matrix: */ 422 GetBStokes(&B[0][0], xyz_list, gauss); 423 /*Create a reduced matrix of B to get rid of pressure */ 424 for (i=0;i<6;i++){ 425 for (j=0;j<3;j++){ 426 B_reduced[i][j]=B[i][j]; 427 } 428 for (j=4;j<7;j++){ 429 B_reduced[i][j-1]=B[i][j]; 430 } 431 for (j=8;j<11;j++){ 432 B_reduced[i][j-2]=B[i][j]; 433 } 434 for (j=12;j<15;j++){ 435 B_reduced[i][j-3]=B[i][j]; 436 } 437 for (j=16;j<19;j++){ 438 B_reduced[i][j-4]=B[i][j]; 439 } 440 for (j=20;j<23;j++){ 441 B_reduced[i][j-5]=B[i][j]; 442 } 443 } 444 445 /*Here, we are computing the strain rate of (vx,0,0)*/ 446 for(i=0;i<numgrids;i++){ 447 velocity[i][0]=this->values[i]; 448 velocity[i][1]=0.0; 449 velocity[i][2]=0.0; 450 } 451 /*Multiply B by velocity, to get strain rate: */ 452 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvx,0); 453 454 } 455 /*}}}*/ 456 /*FUNCTION PentaVertexInput::GetVyStrainRate3d{{{1*/ 457 void PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ 458 int i,j; 459 460 const int numgrids=6; 461 const int DOFVELOCITY=3; 462 double B[8][27]; 463 double B_reduced[6][DOFVELOCITY*numgrids]; 464 double velocity[numgrids][DOFVELOCITY]; 465 466 /*Get B matrix: */ 467 GetBStokes(&B[0][0], xyz_list, gauss); 468 /*Create a reduced matrix of B to get rid of pressure */ 469 for (i=0;i<6;i++){ 470 for (j=0;j<3;j++){ 471 B_reduced[i][j]=B[i][j]; 472 } 473 for (j=4;j<7;j++){ 474 B_reduced[i][j-1]=B[i][j]; 475 } 476 for (j=8;j<11;j++){ 477 B_reduced[i][j-2]=B[i][j]; 478 } 479 for (j=12;j<15;j++){ 480 B_reduced[i][j-3]=B[i][j]; 481 } 482 for (j=16;j<19;j++){ 483 B_reduced[i][j-4]=B[i][j]; 484 } 485 for (j=20;j<23;j++){ 486 B_reduced[i][j-5]=B[i][j]; 487 } 488 } 489 490 /*Here, we are computing the strain rate of (0,vy,0)*/ 491 for(i=0;i<numgrids;i++){ 492 velocity[i][0]=0.0; 493 velocity[i][1]=this->values[i]; 494 velocity[i][2]=0.0; 495 } 496 /*Multiply B by velocity, to get strain rate: */ 497 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvy,0); 498 499 } 500 /*}}}*/ 501 /*FUNCTION PentaVertexInput::GetVzStrainRate3d{{{1*/ 502 void PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){ 503 int i,j; 504 505 const int numgrids=6; 506 const int DOFVELOCITY=3; 507 double B[8][27]; 508 double B_reduced[6][DOFVELOCITY*numgrids]; 509 double velocity[numgrids][DOFVELOCITY]; 510 511 /*Get B matrix: */ 512 GetBStokes(&B[0][0], xyz_list, gauss); 513 /*Create a reduced matrix of B to get rid of pressure */ 514 for (i=0;i<6;i++){ 515 for (j=0;j<3;j++){ 516 B_reduced[i][j]=B[i][j]; 517 } 518 for (j=4;j<7;j++){ 519 B_reduced[i][j-1]=B[i][j]; 520 } 521 for (j=8;j<11;j++){ 522 B_reduced[i][j-2]=B[i][j]; 523 } 524 for (j=12;j<15;j++){ 525 B_reduced[i][j-3]=B[i][j]; 526 } 527 for (j=16;j<19;j++){ 528 B_reduced[i][j-4]=B[i][j]; 529 } 530 for (j=20;j<23;j++){ 531 B_reduced[i][j-5]=B[i][j]; 532 } 533 } 534 535 /*Here, we are computing the strain rate of (0,0,vz)*/ 536 for(i=0;i<numgrids;i++){ 537 velocity[i][0]=0.0; 538 velocity[i][1]=0.0; 539 velocity[i][2]=this->values[i]; 540 } 541 542 /*Multiply B by velocity, to get strain rate: */ 543 MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvz,0); 544 545 } 546 /*}}}*/ 547 /*FUNCTION PentaVertexInput::GetVxStrainRate3dPattyn{{{1*/ 548 void PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ 549 550 int i; 551 const int numgrids=6; 552 const int NDOF2=2; 553 double B[5][NDOF2*numgrids]; 554 double velocity[numgrids][NDOF2]; 555 556 /*Get B matrix: */ 557 GetBPattyn(&B[0][0], xyz_list, gauss); 558 559 /*Here, we are computing the strain rate of (vx,0)*/ 560 for(i=0;i<numgrids;i++){ 561 velocity[i][0]=this->values[i]; 562 velocity[i][1]=0.0; 563 } 564 565 /*Multiply B by velocity, to get strain rate: */ 566 MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0, 567 &velocity[0][0],NDOF2*numgrids,1,0, 568 epsilonvx,0); 569 570 } 571 /*}}}*/ 572 /*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn{{{1*/ 573 void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ 388 574 389 575 int i; -
issm/trunk/src/c/objects/Inputs/PentaVertexInput.h
r5631 r5647 48 48 void GetParameterValue(double* pvalue){ISSMERROR("not implemented yet");}; 49 49 void GetParameterValue(double* pvalue,double* gauss); 50 void GetParameterValue(double* pvalue,GaussTria* gauss); 50 void GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR("not implemented yet");}; 51 void GetParameterValue(double* pvalue,GaussPenta* gauss); 51 52 void GetParameterValues(double* values,double* gauss_pointers, int numgauss); 52 53 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss); 53 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); 54 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; 55 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss); 54 56 void GetParameterAverage(double* pvalue); 55 57 void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");}; … … 62 64 void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; 63 65 void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; 64 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};65 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};66 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};67 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};68 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};66 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss); 67 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss); 68 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss); 69 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss); 70 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss); 69 71 void ChangeEnum(int newenumtype); 70 72 -
issm/trunk/src/c/objects/Inputs/TriaVertexInput.h
r5631 r5647 49 49 void GetParameterValue(double* pvalue,double* gauss); 50 50 void GetParameterValue(double* pvalue,GaussTria* gauss); 51 void GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR("not implemented yet");}; 51 52 void GetParameterValues(double* values,double* gauss_pointers, int numgauss); 52 53 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss); 53 54 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); 55 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 54 56 void GetParameterAverage(double* pvalue); 55 57 void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss); … … 62 64 void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss); 63 65 void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss); 64 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};65 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};66 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};67 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};68 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, Gauss Tria* gauss){ISSMERROR("not implemented yet");};66 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 67 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 68 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 69 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 70 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; 69 71 void ChangeEnum(int newenumtype); 70 72
Note:
See TracChangeset
for help on using the changeset viewer.