Changeset 5625
- Timestamp:
- 08/30/10 17:01:51 (15 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 3 added
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Makefile.am
r5595 r5625 77 77 ./objects/Bamg/Mesh.cpp\ 78 78 ./objects/Bamg/Mesh.h\ 79 ./objects/Gauss/GaussTria.h\ 80 ./objects/Gauss/GaussTria.cpp\ 79 81 ./objects/Update.h\ 80 82 ./objects/Element.h\ … … 630 632 ./objects/Bamg/Mesh.h\ 631 633 ./objects/Bamg/Mesh.cpp\ 634 ./objects/Gauss/GaussTria.h\ 635 ./objects/Gauss/GaussTria.cpp\ 632 636 ./objects/Update.h\ 633 637 ./objects/Element.h\ -
issm/trunk/src/c/objects/Elements/Tria.cpp
r5624 r5625 4580 4580 void Tria::CreateKMatrixSlope(Mat Kgg){ 4581 4581 4582 /* local declarations */ 4583 int i,j; 4584 4585 /* node data: */ 4586 const int numgrids=3; 4582 /*constants: */ 4583 const int numnodes=3; 4587 4584 const int NDOF1=1; 4588 const int numdof=NDOF1*numgrids; 4589 double xyz_list[numgrids][3]; 4590 int* doflist=NULL; 4591 4592 /* gaussian points: */ 4593 int num_gauss,ig; 4594 double* first_gauss_area_coord = NULL; 4595 double* second_gauss_area_coord = NULL; 4596 double* third_gauss_area_coord = NULL; 4597 double* gauss_weights = NULL; 4598 double gauss_weight; 4599 double gauss_l1l2l3[3]; 4600 4601 /* matrices: */ 4602 double L[1][3]; 4603 double DL_scalar; 4604 4605 /* local element matrices: */ 4606 double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix 4607 double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. 4608 4609 double Jdet; 4610 4611 /* Get node coordinates and dof list: */ 4612 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 4585 const int numdof=NDOF1*numnodes; 4586 4587 /* Intermediaries */ 4588 int i,j,ig; 4589 double DL_scalar,Jdet; 4590 double xyz_list[numnodes][3]; 4591 double L[1][3]; 4592 double Ke[numdof][numdof] = {0.0}; 4593 double Ke_g[numdof][numdof]; 4594 GaussTria *gauss = NULL; 4595 int *doflist = NULL; 4596 4597 GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes); 4613 4598 GetDofList(&doflist); 4614 4599 4615 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 4616 GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); 4617 4618 /* Start looping on the number of gaussian points: */ 4619 for (ig=0; ig<num_gauss; ig++){ 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); 4625 4600 /* Start looping on the number of gaussian points: */ 4601 gauss=new GaussTria(2); 4602 for (ig=gauss->begin();ig<gauss->end();ig++){ 4603 4604 gauss->GaussPoint(ig); 4626 4605 4627 /*Get L matrix: */ 4628 GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1); 4629 4630 /* Get Jacobian determinant: */ 4631 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); 4632 4633 DL_scalar=gauss_weight*Jdet; 4634 4635 /* Do the triple producte tL*D*L: */ 4636 TripleMultiply( &L[0][0],1,3,1, 4606 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 4607 DL_scalar=gauss->weight*Jdet; 4608 4609 GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF1); 4610 4611 TripleMultiply(&L[0][0],1,3,1, 4637 4612 &DL_scalar,1,1,0, 4638 4613 &L[0][0],1,3,0, 4639 &Ke_gg_gaussian[0][0],0); 4640 4641 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ 4642 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j]; 4643 } //for (ig=0; ig<num_gauss; ig++ 4644 4645 /*Add Ke_gg to global matrix Kgg: */ 4646 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); 4614 &Ke_g[0][0],0); 4615 4616 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke[i][j]+=Ke_g[i][j]; 4617 } 4618 4619 /*Add Ke to global matrix Kgg: */ 4620 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES); 4647 4621 4648 xfree((void**)&first_gauss_area_coord); 4649 xfree((void**)&second_gauss_area_coord); 4650 xfree((void**)&third_gauss_area_coord); 4651 xfree((void**)&gauss_weights); 4622 /*Clean up*/ 4623 delete gauss; 4652 4624 xfree((void**)&doflist); 4653 4625 } -
issm/trunk/src/c/objects/Elements/TriaRef.cpp
r4921 r5625 180 180 } 181 181 /*}}}*/ 182 /*FUNCTION TriaRef::GetL {{{1*/182 /*FUNCTION TriaRef::GetL(double* L, double* xyz_list, double* gauss,int numdof) {{{1*/ 183 183 void TriaRef::GetL(double* L, double* xyz_list, double* gauss,int numdof){ 184 184 /*Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. … … 221 221 } 222 222 /*}}}*/ 223 /*FUNCTION TriaRef::GetJacobian {{{1*/ 223 /*FUNCTION TriaRef::GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof) {{{1*/ 224 void TriaRef::GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof){ 225 /*Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. 226 * For grid i, Li can be expressed in the actual coordinate system 227 * by: 228 * numdof=1: 229 * Li=h; 230 * numdof=2: 231 * Li=[ h 0 ] 232 * [ 0 h ] 233 * where h is the interpolation function for grid i. 234 * 235 * We assume L has been allocated already, of size: numgrids (numdof=1), or numdofx(numdof*numgrids) (numdof=2) 236 */ 237 238 int i; 239 const int NDOF2=2; 240 const int numgrids=3; 241 double l1l2l3[3]; 242 243 /*Get l1l2l3 in actual coordinate system: */ 244 GetNodalFunctions(l1l2l3,gauss); 245 246 /*Build L: */ 247 if(numdof==1){ 248 for (i=0;i<numgrids;i++){ 249 L[i]=l1l2l3[i]; 250 } 251 } 252 else{ 253 for (i=0;i<numgrids;i++){ 254 *(L+numdof*numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i]; 255 *(L+numdof*numgrids*0+numdof*i+1)=0; 256 *(L+numdof*numgrids*1+numdof*i)=0; 257 *(L+numdof*numgrids*1+numdof*i+1)=l1l2l3[i]; 258 } 259 } 260 } 261 /*}}}*/ 262 /*FUNCTION TriaRef::GetJacobian(double* J, double* xyz_list,double* gauss) {{{1*/ 224 263 void TriaRef::GetJacobian(double* J, double* xyz_list,double* gauss){ 225 264 /*The Jacobian is constant over the element, discard the gaussian points. … … 244 283 } 245 284 /*}}}*/ 246 /*FUNCTION TriaRef::GetJacobianDeterminant2d {{{1*/ 285 /*FUNCTION TriaRef::GetJacobian(double* J, double* xyz_list,GaussTria* gauss) {{{1*/ 286 void TriaRef::GetJacobian(double* J, double* xyz_list,GaussTria* gauss){ 287 /*The Jacobian is constant over the element, discard the gaussian points. 288 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ 289 290 const int NDOF2=2; 291 const int numgrids=3; 292 double x1,y1,x2,y2,x3,y3; 293 294 x1=*(xyz_list+numgrids*0+0); 295 y1=*(xyz_list+numgrids*0+1); 296 x2=*(xyz_list+numgrids*1+0); 297 y2=*(xyz_list+numgrids*1+1); 298 x3=*(xyz_list+numgrids*2+0); 299 y3=*(xyz_list+numgrids*2+1); 300 301 302 *(J+NDOF2*0+0)=0.5*(x2-x1); 303 *(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2); 304 *(J+NDOF2*0+1)=0.5*(y2-y1); 305 *(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2); 306 } 307 /*}}}*/ 308 /*FUNCTION TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss) {{{1*/ 247 309 void TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss){ 248 310 311 /*The Jacobian determinant is constant over the element, discard the gaussian points. 312 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ 313 double J[2][2]; 314 315 /*Get Jacobian*/ 316 GetJacobian(&J[0][0],xyz_list,gauss); 317 318 /*Get Determinant*/ 319 Matrix2x2Determinant(Jdet,&J[0][0]); 320 if(*Jdet<0) ISSMERROR("negative jacobian determinant!"); 321 322 } 323 /*}}}*/ 324 /*FUNCTION TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss) {{{1*/ 325 void TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss){ 249 326 /*The Jacobian determinant is constant over the element, discard the gaussian points. 250 327 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ … … 296 373 } 297 374 /*}}}*/ 298 /*FUNCTION TriaRef::GetNodalFunctions {{{1*/375 /*FUNCTION TriaRef::GetNodalFunctions(double* l1l2l3, double* gauss) {{{1*/ 299 376 void TriaRef::GetNodalFunctions(double* l1l2l3, double* gauss){ 300 377 /*This routine returns the values of the nodal functions at the gaussian point.*/ … … 303 380 l1l2l3[1]=gauss[1]; 304 381 l1l2l3[2]=gauss[2]; 382 383 } 384 /*}}}*/ 385 /*FUNCTION TriaRef::GetNodalFunctions(double* l1l2l3,GaussTria* gauss) {{{1*/ 386 void TriaRef::GetNodalFunctions(double* l1l2l3,GaussTria* gauss){ 387 /*This routine returns the values of the nodal functions at the gaussian point.*/ 388 389 l1l2l3[0]=gauss->coord1; 390 l1l2l3[1]=gauss->coord2; 391 l1l2l3[2]=gauss->coord3; 305 392 306 393 } -
issm/trunk/src/c/objects/Elements/TriaRef.h
r4921 r5625 7 7 #ifndef _TRIAREF_H_ 8 8 #define _TRIAREF_H_ 9 10 class GaussTria; 9 11 10 12 class TriaRef{ … … 27 29 void GetBprimePrognostic(double* Bprime_prog, double* xyz_list, double* gauss); 28 30 void GetBPrognostic(double* B_prog, double* xyz_list, double* gauss); 29 void GetL(double* L, double* xyz_list, double* gauss,int numdof); 31 void GetL(double* L, double* xyz_list,double* gauss,int numdof); 32 void GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof); 30 33 void GetJacobian(double* J, double* xyz_list,double* gauss); 34 void GetJacobian(double* J, double* xyz_list,GaussTria* gauss); 31 35 void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss); 36 void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss); 32 37 void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss); 33 38 void GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss); 34 39 void GetNodalFunctions(double* l1l2l3, double* gauss); 40 void GetNodalFunctions(double* l1l2l3,GaussTria* gauss); 35 41 void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, double* gauss); 36 42 void GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss); -
issm/trunk/src/c/objects/objects.h
r5578 r5625 22 22 /*Constraints: */ 23 23 #include "./Constraints/Spc.h" 24 25 /*Gauss*/ 26 #include "./Gauss/GaussTria.h" 24 27 25 28 /*Loads: */
Note:
See TracChangeset
for help on using the changeset viewer.