Changeset 7000
- Timestamp:
- 01/07/11 16:34:33 (14 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Penta.cpp
r6988 r7000 719 719 const int numdofm=NDOF2*NUMVERTICES2D; 720 720 const int numdofs=NDOF4*NUMVERTICES; 721 const int numdoftotal= numdofm+numdofs;721 const int numdoftotal=2*numdofm+numdofs; 722 722 723 723 /*Intermediaries */ 724 724 int i,j,ig; 725 725 double Jdet; 726 double viscosity ; //viscosity727 double epsilon[ 5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/726 double viscosity,stokesreconditioning; //viscosity 727 double epsilon[6]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ 728 728 double xyz_list[NUMVERTICES][3]; 729 double B[4][numdofs ];729 double B[4][numdofs+3]; 730 730 double Bprime[4][numdofm]; 731 731 double B2[3][numdofm]; 732 double Bprime2[3][numdofs ];732 double Bprime2[3][numdofs+3]; 733 733 double D[4][4]={0.0}; // material matrix, simple scalar matrix. 734 734 double D2[3][3]={0.0}; // material matrix, simple scalar matrix. … … 736 736 double Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix 737 737 double Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix 738 double Ke_gg_gaussian[numdofs ][numdofm]; //stiffness matrix evaluated at the gaussian point.739 double Ke_gg_gaussian2[numdofm][numdofs ]; //stiffness matrix evaluated at the gaussian point.738 double Ke_gg_gaussian[numdofs+3][numdofm]; //stiffness matrix evaluated at the gaussian point. 739 double Ke_gg_gaussian2[numdofm][numdofs+3]; //stiffness matrix evaluated at the gaussian point. 740 740 GaussPenta *gauss=NULL; 741 741 GaussTria *gauss_tria=NULL; … … 754 754 /* Get node coordinates and dof list: */ 755 755 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 756 parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); 756 757 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); 757 758 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); … … 777 778 D_scalar=2*viscosity*gauss->weight*Jdet; 778 779 for (i=0;i<3;i++) D[i][i]=D_scalar; 779 D[ 4][4]=(-1)*gauss->weight*Jdet;780 D[3][3]=-gauss->weight*Jdet*stokesreconditioning; 780 781 for (i=0;i<3;i++) D2[i][i]=D_scalar; 781 782 782 TripleMultiply( &B[0][0],4,numdofs ,1,783 TripleMultiply( &B[0][0],4,numdofs+3,1, 783 784 &D[0][0],4,4,0, 784 785 &Bprime[0][0],4,numdofm,0, … … 787 788 TripleMultiply( &B2[0][0],3,numdofm,1, 788 789 &D2[0][0],3,3,0, 789 &Bprime2[0][0],3,numdofs ,0,790 &Bprime2[0][0],3,numdofs+3,0, 790 791 &Ke_gg_gaussian2[0][0],0); 791 792 … … 794 795 } 795 796 for(i=0;i<numdofs;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j]; 796 for(i=0;i<numdofm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[j][i]; 797 _error_("coupling not finished yet"); 797 for(i=0;i<numdofm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[i][j]; 798 798 799 799 /*Clean-up and return*/ … … 990 990 double viscosity, oldviscosity, newviscosity, viscosity_overshoot; 991 991 double epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ 992 double epsilons[6]; //6 for stokes 992 993 double xyz_list[NUMVERTICES][3]; 993 994 double B[3][numdof2d]; … … 1040 1041 } 1041 1042 else if (approximation==MacAyealStokesApproximationEnum){ 1042 this->GetStrainRate3d(&epsilon [0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);1043 matice->GetViscosity3dStokes(&newviscosity,&epsilon [0]);1043 this->GetStrainRate3d(&epsilons[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 1044 matice->GetViscosity3dStokes(&newviscosity,&epsilons[0]); 1044 1045 } 1045 1046 else _error_("approximation %i (%s) not supported yet",approximation,EnumToString(approximation)); -
issm/trunk/src/c/objects/Elements/PentaRef.cpp
r6987 r7000 134 134 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=l1l6[i]; 135 135 } 136 137 136 } 138 137 /*}}}*/ … … 245 244 *(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=dh1dh7[0][i]; 246 245 *(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=0; 247 *(Bprime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=0;248 *(Bprime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=0;249 *(Bprime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;250 246 } 251 247 -
issm/trunk/src/c/objects/Node.cpp
r6412 r7000 104 104 if (iomodel->vertices_type[io_index]==MacAyealStokesApproximationEnum && iomodel->gridonmacayeal[io_index]){ 105 105 if(!iomodel->gridonbed[io_index]){ 106 for(k=1;k<= gsize;k++) this->FreezeDof(k);106 for(k=1;k<=2;k++) this->FreezeDof(k); 107 107 } 108 108 }
Note:
See TracChangeset
for help on using the changeset viewer.