![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 4769
- Timestamp:
- 07/22/10 15:46:39 (15 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 6 edited
-
Elements/Tria.cpp (modified) (3 diffs)
-
Elements/Tria.h (modified) (2 diffs)
-
Elements/TriaRef.cpp (modified) (1 diff)
-
Elements/TriaRef.h (modified) (1 diff)
-
Inputs/TriaVertexInput.cpp (modified) (2 diffs)
-
Inputs/TriaVertexInput.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Tria.cpp
r4765 r4769 4863 4863 } 4864 4864 /*}}}*/ 4865 /*FUNCTION Tria::GetB {{{1*/4866 4867 void Tria::GetB(double* B, double* xyz_list, double* gauss_l1l2l3){4868 4869 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.4870 * For grid i, Bi can be expressed in the actual coordinate system4871 * by:4872 * Bi=[ dh/dx 0 ]4873 * [ 0 dh/dy ]4874 * [ 1/2*dh/dy 1/2*dh/dx ]4875 * where h is the interpolation function for grid i.4876 *4877 * We assume B has been allocated already, of size: 3x(NDOF2*numgrids)4878 */4879 4880 int i;4881 const int NDOF2=2;4882 const int numgrids=3;4883 4884 double dh1dh3[NDOF2][numgrids];4885 4886 4887 /*Get dh1dh2dh3 in actual coordinate system: */4888 GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss_l1l2l3);4889 4890 /*Build B: */4891 for (i=0;i<numgrids;i++){4892 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];4893 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0;4894 *(B+NDOF2*numgrids*1+NDOF2*i)=0;4895 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh3[1][i];4896 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i];4897 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i];4898 }4899 }4900 /*}}}*/4901 4865 /*FUNCTION Tria::GetB_prog {{{1*/ 4902 4866 … … 5042 5006 } 5043 5007 /*}}}*/ 5044 /*FUNCTION Tria::GetJacobian {{{1*/5045 void Tria::GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3){5046 5047 /*The Jacobian is constant over the element, discard the gaussian points.5048 * J is assumed to have been allocated of size NDOF2xNDOF2.*/5049 5050 const int NDOF2=2;5051 const int numgrids=3;5052 double x1,y1,x2,y2,x3,y3;5053 5054 x1=*(xyz_list+numgrids*0+0);5055 y1=*(xyz_list+numgrids*0+1);5056 x2=*(xyz_list+numgrids*1+0);5057 y2=*(xyz_list+numgrids*1+1);5058 x3=*(xyz_list+numgrids*2+0);5059 y3=*(xyz_list+numgrids*2+1);5060 5061 5062 *(J+NDOF2*0+0)=0.5*(x2-x1);5063 *(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);5064 *(J+NDOF2*0+1)=0.5*(y2-y1);5065 *(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);5066 }5067 /*}}}*/5068 /*FUNCTION Tria::GetJacobianDeterminant2d {{{1*/5069 void Tria::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss_l1l2l3){5070 5071 /*The Jacobian determinant is constant over the element, discard the gaussian points.5072 * J is assumed to have been allocated of size NDOF2xNDOF2.*/5073 5074 double x1,x2,x3,y1,y2,y3;5075 5076 x1=*(xyz_list+3*0+0);5077 y1=*(xyz_list+3*0+1);5078 x2=*(xyz_list+3*1+0);5079 y2=*(xyz_list+3*1+1);5080 x3=*(xyz_list+3*2+0);5081 y3=*(xyz_list+3*2+1);5082 5083 *Jdet=SQRT3/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1));5084 5085 if(Jdet<0){5086 ISSMERROR("negative jacobian determinant!");5087 }5088 5089 }5090 /*}}}*/5091 /*FUNCTION Tria::GetJacobianDeterminant3d {{{1*/5092 void Tria::GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss_l1l2l3){5093 5094 /*The Jacobian determinant is constant over the element, discard the gaussian points.5095 * J is assumed to have been allocated of size NDOF2xNDOF2.*/5096 5097 double x1,x2,x3,y1,y2,y3,z1,z2,z3;5098 5099 x1=*(xyz_list+3*0+0);5100 y1=*(xyz_list+3*0+1);5101 z1=*(xyz_list+3*0+2);5102 x2=*(xyz_list+3*1+0);5103 y2=*(xyz_list+3*1+1);5104 z2=*(xyz_list+3*1+2);5105 x3=*(xyz_list+3*2+0);5106 y3=*(xyz_list+3*2+1);5107 z3=*(xyz_list+3*2+2);5108 5109 5110 *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);5111 5112 if(Jdet<0){5113 ISSMERROR("negative jacobian determinant!");5114 }5115 5116 }5117 /*}}}*/5118 /*FUNCTION Tria::GetJacobianInvert {{{1*/5119 void Tria::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss_l1l2l3){5120 5121 double Jdet;5122 const int NDOF2=2;5123 const int numgrids=3;5124 5125 /*Call Jacobian routine to get the jacobian:*/5126 GetJacobian(Jinv, xyz_list, gauss_l1l2l3);5127 5128 /*Invert Jacobian matrix: */5129 MatrixInverse(Jinv,NDOF2,NDOF2,NULL,0,&Jdet);5130 5131 }5132 /*}}}*/5133 5008 /*FUNCTION Tria::GetL {{{1*/ 5134 5009 … … 5178 5053 } 5179 5054 } 5180 }5181 /*}}}*/5182 /*FUNCTION Tria::GetNodalFunctions {{{1*/5183 void Tria::GetNodalFunctions(double* l1l2l3, double* gauss_l1l2l3){5184 5185 /*This routine returns the values of the nodal functions at the gaussian point.*/5186 5187 /*First nodal function: */5188 l1l2l3[0]=gauss_l1l2l3[0];5189 5190 /*Second nodal function: */5191 l1l2l3[1]=gauss_l1l2l3[1];5192 5193 /*Third nodal function: */5194 l1l2l3[2]=gauss_l1l2l3[2];5195 5196 }5197 /*}}}*/5198 /*FUNCTION Tria::GetNodalFunctionsDerivatives {{{1*/5199 void Tria::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss_l1l2l3){5200 5201 /*This routine returns the values of the nodal functions derivatives (with respect to the5202 * actual coordinate system: */5203 5204 int i;5205 const int NDOF2=2;5206 const int numgrids=3;5207 5208 double dh1dh3_ref[NDOF2][numgrids];5209 double Jinv[NDOF2][NDOF2];5210 5211 5212 /*Get derivative values with respect to parametric coordinate system: */5213 GetNodalFunctionsDerivativesReference(&dh1dh3_ref[0][0], gauss_l1l2l3);5214 5215 /*Get Jacobian invert: */5216 GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_l1l2l3);5217 5218 /*Build dh1dh3:5219 *5220 * [dhi/dx]= Jinv*[dhi/dr]5221 * [dhi/dy] [dhi/ds]5222 */5223 5224 for (i=0;i<numgrids;i++){5225 *(dh1dh3+numgrids*0+i)=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];5226 *(dh1dh3+numgrids*1+i)=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];5227 }5228 5229 }5230 /*}}}*/5231 /*FUNCTION Tria::GetNodalFunctionsDerivativesReference {{{1*/5232 void Tria::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3){5233 5234 /*This routine returns the values of the nodal functions derivatives (with respect to the5235 * natural coordinate system) at the gaussian point. */5236 5237 const int NDOF2=2;5238 const int numgrids=3;5239 5240 /*First nodal function: */5241 *(dl1dl3+numgrids*0+0)=-0.5;5242 *(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3);5243 5244 /*Second nodal function: */5245 *(dl1dl3+numgrids*0+1)=0.5;5246 *(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3);5247 5248 /*Third nodal function: */5249 *(dl1dl3+numgrids*0+2)=0;5250 *(dl1dl3+numgrids*1+2)=1.0/SQRT3;5251 5252 5055 } 5253 5056 /*}}}*/ -
issm/trunk/src/c/objects/Elements/Tria.h
r4765 r4769 136 136 double GetArea(void); 137 137 double GetAreaCoordinate(double x, double y, int which_one); 138 void GetB(double* B, double* xyz_list, double* gauss_l1l2l3);139 138 void GetBPrime(double* Bprime, double* xyz_list, double* gauss_l1l2l3); 140 139 void GetBPrime_prog(double* Bprime_prog, double* xyz_list, double* gauss_l1l2l3); … … 143 142 void GetDofList(int* doflist,int* pnumberofdofs); 144 143 void GetDofList1(int* doflist); 145 void GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3);146 void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss_l1l2l3);147 void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss_l1l2l3);148 void GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss_l1l2l3);149 144 void GetL(double* L, double* xyz_list, double* gauss_l1l2l3,int numdof); 150 void GetNodalFunctions(double* l1l2l3, double* gauss_l1l2l3);151 void GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss_l1l2l3);152 void GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3);153 145 void GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3); 154 146 void GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3); -
issm/trunk/src/c/objects/Elements/TriaRef.cpp
r4739 r4769 50 50 } 51 51 /*}}}*/ 52 53 /*Reference Element numerics*/ 54 /*FUNCTION TriaRef::GetB {{{1*/ 55 void TriaRef::GetB(double* B, double* xyz_list, double* gauss){ 56 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 57 * For grid i, Bi can be expressed in the actual coordinate system 58 * by: 59 * Bi=[ dh/dx 0 ] 60 * [ 0 dh/dy ] 61 * [ 1/2*dh/dy 1/2*dh/dx ] 62 * where h is the interpolation function for grid i. 63 * 64 * We assume B has been allocated already, of size: 3x(NDOF2*numgrids) 65 */ 66 67 int i; 68 const int NDOF2=2; 69 const int numgrids=3; 70 71 double dh1dh3[NDOF2][numgrids]; 72 73 74 /*Get dh1dh2dh3 in actual coordinate system: */ 75 GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss); 76 77 /*Build B: */ 78 for (i=0;i<numgrids;i++){ 79 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i]; 80 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0; 81 *(B+NDOF2*numgrids*1+NDOF2*i)=0; 82 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh3[1][i]; 83 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i]; 84 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i]; 85 } 86 } 87 /*}}}*/ 88 /*FUNCTION TriaRef::GetJacobian {{{1*/ 89 void TriaRef::GetJacobian(double* J, double* xyz_list,double* gauss){ 90 /*The Jacobian is constant over the element, discard the gaussian points. 91 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ 92 93 const int NDOF2=2; 94 const int numgrids=3; 95 double x1,y1,x2,y2,x3,y3; 96 97 x1=*(xyz_list+numgrids*0+0); 98 y1=*(xyz_list+numgrids*0+1); 99 x2=*(xyz_list+numgrids*1+0); 100 y2=*(xyz_list+numgrids*1+1); 101 x3=*(xyz_list+numgrids*2+0); 102 y3=*(xyz_list+numgrids*2+1); 103 104 105 *(J+NDOF2*0+0)=0.5*(x2-x1); 106 *(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2); 107 *(J+NDOF2*0+1)=0.5*(y2-y1); 108 *(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2); 109 } 110 /*}}}*/ 111 /*FUNCTION TriaRef::GetJacobianDeterminant2d {{{1*/ 112 void TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss){ 113 114 /*The Jacobian determinant is constant over the element, discard the gaussian points. 115 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ 116 117 double x1,x2,x3,y1,y2,y3; 118 119 x1=*(xyz_list+3*0+0); 120 y1=*(xyz_list+3*0+1); 121 x2=*(xyz_list+3*1+0); 122 y2=*(xyz_list+3*1+1); 123 x3=*(xyz_list+3*2+0); 124 y3=*(xyz_list+3*2+1); 125 126 *Jdet=SQRT3/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)); 127 128 if(Jdet<0){ 129 ISSMERROR("negative jacobian determinant!"); 130 } 131 132 } 133 /*}}}*/ 134 /*FUNCTION TriaRef::GetJacobianDeterminant3d {{{1*/ 135 void TriaRef::GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss){ 136 /*The Jacobian determinant is constant over the element, discard the gaussian points. 137 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ 138 139 double x1,x2,x3,y1,y2,y3,z1,z2,z3; 140 141 x1=*(xyz_list+3*0+0); 142 y1=*(xyz_list+3*0+1); 143 z1=*(xyz_list+3*0+2); 144 x2=*(xyz_list+3*1+0); 145 y2=*(xyz_list+3*1+1); 146 z2=*(xyz_list+3*1+2); 147 x3=*(xyz_list+3*2+0); 148 y3=*(xyz_list+3*2+1); 149 z3=*(xyz_list+3*2+2); 150 151 152 *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); 153 154 if(Jdet<0){ 155 ISSMERROR("negative jacobian determinant!"); 156 } 157 158 } 159 /*}}}*/ 160 /*FUNCTION TriaRef::GetJacobianInvert {{{1*/ 161 void TriaRef::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss){ 162 163 double Jdet; 164 const int NDOF2=2; 165 const int numgrids=3; 166 167 /*Call Jacobian routine to get the jacobian:*/ 168 GetJacobian(Jinv, xyz_list, gauss); 169 170 /*Invert Jacobian matrix: */ 171 MatrixInverse(Jinv,NDOF2,NDOF2,NULL,0,&Jdet); 172 173 } 174 /*}}}*/ 175 /*FUNCTION TriaRef::GetNodalFunctions {{{1*/ 176 void TriaRef::GetNodalFunctions(double* l1l2l3, double* gauss){ 177 178 /*This routine returns the values of the nodal functions at the gaussian point.*/ 179 180 /*First nodal function: */ 181 l1l2l3[0]=gauss[0]; 182 183 /*Second nodal function: */ 184 l1l2l3[1]=gauss[1]; 185 186 /*Third nodal function: */ 187 l1l2l3[2]=gauss[2]; 188 189 } 190 /*}}}*/ 191 /*FUNCTION TriaRef::GetNodalFunctionsDerivatives {{{1*/ 192 void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss){ 193 194 /*This routine returns the values of the nodal functions derivatives (with respect to the 195 * actual coordinate system: */ 196 int i; 197 const int NDOF2=2; 198 const int numgrids=3; 199 double dh1dh3_ref[NDOF2][numgrids]; 200 double Jinv[NDOF2][NDOF2]; 201 202 /*Get derivative values with respect to parametric coordinate system: */ 203 GetNodalFunctionsDerivativesReference(&dh1dh3_ref[0][0], gauss); 204 205 /*Get Jacobian invert: */ 206 GetJacobianInvert(&Jinv[0][0], xyz_list, gauss); 207 208 /*Build dh1dh3: 209 * 210 * [dhi/dx]= Jinv*[dhi/dr] 211 * [dhi/dy] [dhi/ds] 212 */ 213 for (i=0;i<numgrids;i++){ 214 dh1dh3[numgrids*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i]; 215 dh1dh3[numgrids*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i]; 216 } 217 218 } 219 /*}}}*/ 220 /*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference {{{1*/ 221 void TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss){ 222 223 /*This routine returns the values of the nodal functions derivatives (with respect to the 224 * natural coordinate system) at the gaussian point. */ 225 226 const int NDOF2=2; 227 const int numgrids=3; 228 229 /*First nodal function: */ 230 *(dl1dl3+numgrids*0+0)=-0.5; 231 *(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3); 232 233 /*Second nodal function: */ 234 *(dl1dl3+numgrids*0+1)=0.5; 235 *(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3); 236 237 /*Third nodal function: */ 238 *(dl1dl3+numgrids*0+2)=0; 239 *(dl1dl3+numgrids*1+2)=1.0/SQRT3; 240 241 } 242 /*}}}*/ -
issm/trunk/src/c/objects/Elements/TriaRef.h
r4739 r4769 19 19 ~TriaRef(); 20 20 21 /*Management*/ 21 22 void SetElementType(int type,int type_counter); 23 24 /*Numerics*/ 25 void GetB(double* B, double* xyz_list, double* gauss); 26 void GetJacobian(double* J, double* xyz_list,double* gauss); 27 void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss); 28 void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss); 29 void GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss); 30 void GetNodalFunctions(double* l1l2l3, double* gauss); 31 void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, double* gauss); 32 void GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss); 22 33 23 34 }; -
issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp
r4546 r4769 24 24 /*}}}*/ 25 25 /*FUNCTION TriaVertexInput::TriaVertexInput(int in_enum_type,double* values){{{1*/ 26 TriaVertexInput::TriaVertexInput(int in_enum_type,double* in_values){ 27 26 TriaVertexInput::TriaVertexInput(int in_enum_type,double* in_values) 27 :TriaRef(1) 28 { 29 30 /*Set TriaRef*/ 31 this->SetElementType(P1Enum,0); 32 this->element_type=P1Enum; 33 34 /*Set Enum*/ 28 35 enum_type=in_enum_type; 36 37 /*Set values*/ 29 38 values[0]=in_values[0]; 30 39 values[1]=in_values[1]; … … 311 320 312 321 /*Intermediary*/ 313 /*FUNCTION TriaVertexInput::GetNodalFunctions {{{1*/314 void TriaVertexInput::GetNodalFunctions(double* l1l2l3, double* gauss){315 316 /*This routine returns the values of the nodal functions at the gaussian point.*/317 318 /*First nodal function: */319 l1l2l3[0]=gauss[0];320 321 /*Second nodal function: */322 l1l2l3[1]=gauss[1];323 324 /*Third nodal function: */325 l1l2l3[2]=gauss[2];326 327 }328 /*}}}*/329 /*FUNCTION TriaVertexInput::GetB {{{1*/330 void TriaVertexInput::GetB(double* B, double* xyz_list, double* gauss){331 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.332 * For grid i, Bi can be expressed in the actual coordinate system333 * by:334 * Bi=[ dh/dx 0 ]335 * [ 0 dh/dy ]336 * [ 1/2*dh/dy 1/2*dh/dx ]337 * where h is the interpolation function for grid i.338 *339 * We assume B has been allocated already, of size: 3x(NDOF2*numgrids)340 */341 342 int i;343 const int NDOF2=2;344 const int numgrids=3;345 346 double dh1dh3[NDOF2][numgrids];347 348 349 /*Get dh1dh2dh3 in actual coordinate system: */350 GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);351 352 /*Build B: */353 for (i=0;i<numgrids;i++){354 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];355 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0;356 *(B+NDOF2*numgrids*1+NDOF2*i)=0;357 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh3[1][i];358 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i];359 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i];360 }361 }362 /*}}}*/363 /*FUNCTION TriaVertexInput::GetNodalFunctionsDerivatives {{{1*/364 void TriaVertexInput::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss){365 366 /*This routine returns the values of the nodal functions derivatives (with respect to the367 * actual coordinate system: */368 int i;369 const int NDOF2=2;370 const int numgrids=3;371 double dh1dh3_ref[NDOF2][numgrids];372 double Jinv[NDOF2][NDOF2];373 374 /*Get derivative values with respect to parametric coordinate system: */375 GetNodalFunctionsDerivativesReference(&dh1dh3_ref[0][0], gauss);376 377 /*Get Jacobian invert: */378 GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);379 380 /*Build dh1dh3:381 *382 * [dhi/dx]= Jinv*[dhi/dr]383 * [dhi/dy] [dhi/ds]384 */385 for (i=0;i<numgrids;i++){386 dh1dh3[numgrids*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];387 dh1dh3[numgrids*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];388 }389 390 }391 /*}}}*/392 /*FUNCTION TriaVertexInput::GetNodalFunctionsDerivativesReference {{{1*/393 void TriaVertexInput::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3){394 395 /*This routine returns the values of the nodal functions derivatives (with respect to the396 * natural coordinate system) at the gaussian point. */397 398 const int NDOF2=2;399 const int numgrids=3;400 401 /*First nodal function: */402 *(dl1dl3+numgrids*0+0)=-0.5;403 *(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3);404 405 /*Second nodal function: */406 *(dl1dl3+numgrids*0+1)=0.5;407 *(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3);408 409 /*Third nodal function: */410 *(dl1dl3+numgrids*0+2)=0;411 *(dl1dl3+numgrids*1+2)=1.0/SQRT3;412 413 }414 /*}}}*/415 /*FUNCTION TriaVertexInput::GetJacobian {{{1*/416 void TriaVertexInput::GetJacobian(double* J, double* xyz_list,double* gauss){417 418 /*The Jacobian is constant over the element, discard the gaussian points.419 * J is assumed to have been allocated of size NDOF2xNDOF2.*/420 421 const int NDOF2=2;422 const int numgrids=3;423 double x1,y1,x2,y2,x3,y3;424 425 x1=*(xyz_list+numgrids*0+0);426 y1=*(xyz_list+numgrids*0+1);427 x2=*(xyz_list+numgrids*1+0);428 y2=*(xyz_list+numgrids*1+1);429 x3=*(xyz_list+numgrids*2+0);430 y3=*(xyz_list+numgrids*2+1);431 432 433 *(J+NDOF2*0+0)=0.5*(x2-x1);434 *(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);435 *(J+NDOF2*0+1)=0.5*(y2-y1);436 *(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);437 }438 /*}}}*/439 /*FUNCTION TriaVertexInput::GetJacobianInvert {{{1*/440 void TriaVertexInput::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss){441 442 double Jdet;443 const int NDOF2=2;444 const int numgrids=3;445 446 /*Call Jacobian routine to get the jacobian:*/447 GetJacobian(Jinv, xyz_list, gauss);448 449 /*Invert Jacobian matrix: */450 MatrixInverse(Jinv,NDOF2,NDOF2,NULL,0,&Jdet);451 452 }453 /*}}}*/454 322 /*FUNCTION TriaVertexInput::SquareMin{{{1*/ 455 323 void TriaVertexInput::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){ -
issm/trunk/src/c/objects/Inputs/TriaVertexInput.h
r4704 r4769 65 65 void ChangeEnum(int newenumtype); 66 66 67 void GetNodalFunctions(double* l1l2l3, double* gauss);68 void GetB(double* B, double* xyz_list, double* gauss);69 void GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss);70 void GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3);71 void GetJacobian(double* J, double* xyz_list,double* gauss);72 void GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss);73 67 void SquareMin(double* psquaremin, bool process_units,Parameters* parameters); 74 68 void Scale(double scale_factor);
Note:
See TracChangeset
for help on using the changeset viewer.
