![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 15436
- Timestamp:
- 07/05/13 09:22:25 (12 years ago)
- Location:
- issm/trunk-jpl/src/c/classes/Elements
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified issm/trunk-jpl/src/c/classes/Elements/Penta.cpp ¶
r15433 r15436 6816 6816 /*compute all stiffness matrices for this element*/ 6817 6817 Ke1=CreateKMatrixDiagnosticStokesGLSViscous(); 6818 Ke2=CreateKMatrixDiagnosticStokes GLSFriction();6818 Ke2=CreateKMatrixDiagnosticStokesFriction(); 6819 6819 Ke =new ElementMatrix(Ke1,Ke2); 6820 6820 break; … … 6894 6894 ElementMatrix* Penta::CreateKMatrixDiagnosticStokesGLSViscous(void){ 6895 6895 6896 int numdof = NDOF4*NUMVERTICES; 6897 6896 6898 /*Intermediaries */ 6897 int i, approximation;6899 int i,j,approximation; 6898 6900 IssmDouble Jdet,viscosity,stokesreconditioning; 6899 6901 IssmDouble xyz_list[NUMVERTICES][3]; … … 6903 6905 IssmDouble D_scalar; 6904 6906 IssmDouble D[8][8]={0.0}; 6905 IssmDouble Ke_temp[24][24]={ 1.0}; //for the six nodes6907 IssmDouble Ke_temp[24][24]={0.0}; //for the six nodes 6906 6908 GaussPenta *gauss=NULL; 6907 6909 … … 6939 6941 &B_prime[0][0],8,24,0, 6940 6942 &Ke_temp[0][0],1); 6943 6944 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_temp[i][j]; 6941 6945 } 6942 6946 … … 6951 6955 /*FUNCTION Penta::CreateKMatrixDiagnosticStokesFriction{{{*/ 6952 6956 ElementMatrix* Penta::CreateKMatrixDiagnosticStokesFriction(void){ 6953 6954 /*Constants*/6955 const int numdof=NUMVERTICES*NDOF4;6956 const int numdof2d=NUMVERTICES2D*NDOF4;6957 6958 /*Intermediaries */6959 int i,j;6960 int analysis_type,approximation;6961 IssmDouble alpha2,Jdet2d;6962 IssmDouble stokesreconditioning,viscosity;6963 IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/6964 IssmDouble xyz_list[NUMVERTICES][3];6965 IssmDouble xyz_list_tria[NUMVERTICES2D][3];6966 IssmDouble LStokes[2][numdof2d];6967 IssmDouble DLStokes[2][2]={0.0};6968 IssmDouble Ke_drag_gaussian[numdof2d][numdof2d];6969 Friction* friction=NULL;6970 GaussPenta *gauss=NULL;6971 6972 /*If on water or not Stokes, skip stiffness: */6973 inputs->GetInputValue(&approximation,ApproximationEnum);6974 if(IsFloating() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum)) return NULL;6975 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);6976 6977 /*Retrieve all inputs and parameters*/6978 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);6979 parameters->FindParam(&analysis_type,AnalysisTypeEnum);6980 parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);6981 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);6982 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);6983 Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);6984 for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];6985 6986 /*build friction object, used later on: */6987 friction=new Friction("3d",inputs,matpar,analysis_type);6988 6989 /* Start looping on the number of gaussian points: */6990 gauss=new GaussPenta(0,1,2,2);6991 for(int ig=gauss->begin();ig<gauss->end();ig++){6992 6993 gauss->GaussPoint(ig);6994 6995 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);6996 GetLStokes(&LStokes[0][0], gauss);6997 6998 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);6999 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);7000 7001 friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);7002 7003 DLStokes[0][0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx7004 DLStokes[1][1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy7005 7006 TripleMultiply( &LStokes[0][0],2,numdof2d,1,7007 &DLStokes[0][0],2,2,0,7008 &LStokes[0][0],2,numdof2d,0,7009 &Ke_drag_gaussian[0][0],0);7010 7011 for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof+j]+=Ke_drag_gaussian[i][j];7012 }7013 7014 /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/7015 //TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);7016 7017 /*Clean up and return*/7018 delete gauss;7019 delete friction;7020 return Ke;7021 }7022 /*}}}*/7023 /*FUNCTION Penta::CreateKMatrixDiagnosticStokesGLSFriction{{{*/7024 ElementMatrix* Penta::CreateKMatrixDiagnosticStokesGLSFriction(void){7025 6957 7026 6958 /*Constants*/ … … 7740 7672 ElementVector* Penta::CreatePVectorDiagnosticStokes(void){ 7741 7673 7742 /*compute all load vectors for this element*/ 7743 ElementVector* pe1=CreatePVectorDiagnosticStokesViscous(); 7744 ElementVector* pe2=CreatePVectorDiagnosticStokesShelf(); 7745 ElementVector* pe =new ElementVector(pe1,pe2); 7674 int fe_stokes; 7675 ElementVector* pe1; 7676 ElementVector* pe2; 7677 ElementVector* pe; 7678 parameters->FindParam(&fe_stokes,FlowequationFeStokesEnum); 7679 7680 switch(fe_stokes){ 7681 case 0: 7682 /*compute all stiffness matrices for this element*/ 7683 pe1=CreatePVectorDiagnosticStokesViscous(); 7684 pe2=CreatePVectorDiagnosticStokesShelf(); 7685 pe =new ElementVector(pe1,pe2); 7686 break; 7687 case 1: 7688 /*compute all stiffness matrices for this element*/ 7689 pe1=CreatePVectorDiagnosticStokesGLSViscous(); 7690 pe2=CreatePVectorDiagnosticStokesShelf(); 7691 pe =new ElementVector(pe1,pe2); 7692 break; 7693 default: 7694 _error_("Finite element" << fe_stokes << " not supported yet"); 7695 } 7696 7746 7697 7747 7698 /*clean-up and return*/ … … 7830 7781 } 7831 7782 7783 if(IsOnBed() || IsOnSurface()){ 7784 IssmDouble xyz_list_tria[NUMVERTICES2D][3]; 7785 IssmDouble pi=3.141592653589793; 7786 IssmDouble x,y,z; 7787 IssmDouble basis[6]; //for the six nodes of the penta 7788 GaussPenta *gauss=NULL; 7789 Input* x_input=inputs->GetInput(MeshXEnum); _assert_(x_input); 7790 Input* y_input=inputs->GetInput(MeshYEnum); _assert_(y_input); 7791 Input* z_input=inputs->GetInput(MeshZEnum); _assert_(z_input); 7792 7793 7794 if(IsOnBed()){ 7795 for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j]; 7796 gauss=new GaussPenta(0,1,2,2); 7797 } 7798 if(IsOnSurface()){ 7799 for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i+3][j]; 7800 gauss=new GaussPenta(3,4,5,2); 7801 } 7802 7803 7804 for(int ig=gauss->begin();ig<gauss->end();ig++){ 7805 7806 gauss->GaussPoint(ig); 7807 7808 GetTriaJacobianDeterminant(&Jdet, &xyz_list_tria[0][0], gauss); 7809 GetNodalFunctionsP1(basis, gauss); 7810 x_input->GetInputValue(&x, gauss); 7811 y_input->GetInputValue(&y, gauss); 7812 z_input->GetInputValue(&z, gauss); 7813 7814 7815 forcex=-((cos(2*pi*z)-1)*2*pi*sin(2*pi*y)*cos(2*pi*x)-2*(cos(2*pi*x)-1)*2*pi*sin(2*pi*y)*cos(2*pi*z))*1; 7816 forcey=-((cos(2*pi*z)-1)*2*pi*sin(2*pi*x)*cos(2*pi*y)+(cos(2*pi*y)-1)*2*pi*sin(2*pi*x)*cos(2*pi*z))*1; 7817 forcez=-(-2*pi*2*sin(2*pi*x)*sin(2*pi*y)*sin(2*pi*z)+sin(2*pi*x)*sin(2*pi*y)*sin(2*pi*z)); 7818 7819 if(IsOnBed()){ 7820 for(i=0;i<3;i++){ 7821 Pe_gaussian[i*NDOF4+0]+=forcex*Jdet*gauss->weight*basis[i]*-1.; 7822 Pe_gaussian[i*NDOF4+1]+=forcey*Jdet*gauss->weight*basis[i]*-1.; 7823 Pe_gaussian[i*NDOF4+2]+=forcez*Jdet*gauss->weight*basis[i]*-1.; 7824 } 7825 } 7826 if(IsOnSurface()){ 7827 for(i=3;i<6;i++){ 7828 Pe_gaussian[i*NDOF4+0]+=forcex*Jdet*gauss->weight*basis[i]*1.; 7829 Pe_gaussian[i*NDOF4+1]+=forcey*Jdet*gauss->weight*basis[i]*1.; 7830 Pe_gaussian[i*NDOF4+2]+=forcez*Jdet*gauss->weight*basis[i]*1.; 7831 } 7832 } 7833 } 7834 } 7835 7832 7836 /*Condensation*/ 7833 7837 ReduceVectorStokes(pe->values, &Ke_temp[0][0], &Pe_gaussian[0]); 7838 7839 /*Transform coordinate system*/ 7840 TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum); 7841 7842 /*Clean up and return*/ 7843 delete gauss; 7844 return pe; 7845 } 7846 /*}}}*/ 7847 /*FUNCTION Penta::CreatePVectorDiagnosticStokesGLSViscous {{{*/ 7848 ElementVector* Penta::CreatePVectorDiagnosticStokesGLSViscous(void){ 7849 7850 /*Constants*/ 7851 const int numdof=NDOF4*NUMVERTICES; 7852 7853 /*Intermediaries*/ 7854 int i,j; 7855 int approximation; 7856 IssmDouble Jdet,gravity,rho_ice; 7857 IssmDouble forcex,forcey,forcez; 7858 IssmDouble xyz_list[NUMVERTICES][3]; 7859 IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ 7860 IssmDouble l1l6[6]; //for the six nodes and the bubble 7861 IssmDouble Pe_gaussian[numdof]={0.0}; //for the six nodes and the bubble 7862 GaussPenta *gauss=NULL; 7863 7864 /*Initialize Element vector and return if necessary*/ 7865 inputs->GetInputValue(&approximation,ApproximationEnum); 7866 if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL; 7867 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); 7868 7869 /*Retrieve all inputs and parameters*/ 7870 rho_ice=matpar->GetRhoIce(); 7871 gravity=matpar->GetG(); 7872 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 7873 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); 7874 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); 7875 Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); 7876 Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum); _assert_(loadingforcex_input); 7877 Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input); 7878 Input* loadingforcez_input=inputs->GetInput(LoadingforceZEnum); _assert_(loadingforcez_input); 7879 7880 /* Start looping on the number of gaussian points: */ 7881 gauss=new GaussPenta(5,5); 7882 for(int ig=gauss->begin();ig<gauss->end();ig++){ 7883 7884 gauss->GaussPoint(ig); 7885 7886 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); 7887 GetNodalFunctionsP1(&l1l6[0], gauss); 7888 7889 loadingforcex_input->GetInputValue(&forcex, gauss); 7890 loadingforcey_input->GetInputValue(&forcey, gauss); 7891 loadingforcez_input->GetInputValue(&forcez, gauss); 7892 7893 for(i=0;i<NUMVERTICES;i++){ 7894 Pe_gaussian[i*NDOF4+2]+=-rho_ice*gravity*Jdet*gauss->weight*l1l6[i]; 7895 Pe_gaussian[i*NDOF4+0]+=forcex*Jdet*gauss->weight*l1l6[i]; 7896 Pe_gaussian[i*NDOF4+1]+=forcey*Jdet*gauss->weight*l1l6[i]; 7897 Pe_gaussian[i*NDOF4+2]+=forcez*Jdet*gauss->weight*l1l6[i]; 7898 } 7899 } 7900 7901 for(i=0;i<NUMVERTICES;i++) for(j=0;j<3;j++) pe->values[i*NDOF4+j]+=Pe_gaussian[i*NDOF4+j]; 7834 7902 7835 7903 /*Transform coordinate system*/ -
TabularUnified issm/trunk-jpl/src/c/classes/Elements/Penta.h ¶
r15433 r15436 287 287 ElementVector* CreatePVectorDiagnosticStokes(void); 288 288 ElementVector* CreatePVectorDiagnosticStokesViscous(void); 289 ElementVector* CreatePVectorDiagnosticStokesGLSViscous(void); 289 290 ElementVector* CreatePVectorDiagnosticStokesShelf(void); 290 291 ElementVector* CreatePVectorDiagnosticVert(void); -
TabularUnified issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp ¶
r15425 r15436 275 275 /*Build B: */ 276 276 for (i=0;i<NUMNODESMINI;i++){ 277 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh7[0][i]; //B[0][NDOF4*i]=dh1dh6[0][i];278 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0;279 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;280 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0;281 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh7[1][i];282 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;283 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0;284 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0;285 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=dh1dh7[2][i];286 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=(float).5*dh1dh7[1][i];287 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=(float).5*dh1dh7[0][i];288 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;289 *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i)=(float).5*dh1dh7[2][i];290 *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1)=0;291 *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2)=(float).5*dh1dh7[0][i];292 *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i)=0;293 *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1)=(float).5*dh1dh7[2][i];294 *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2)=(float).5*dh1dh7[1][i];295 *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i)=0;296 *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1)=0;297 *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2)=0;298 *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i)=dh1dh7[0][i];299 *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1)=dh1dh7[1][i];300 *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2)=dh1dh7[2][i];277 B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dh1dh7[0][i+0]; //B[0][NDOF4*i+0] = dh1dh6[0][i+0]; 278 B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = 0.; 279 B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.; 280 B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = 0.; 281 B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dh1dh7[1][i+0]; 282 B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.; 283 B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = 0.; 284 B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = 0.; 285 B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = dh1dh7[2][i+0]; 286 B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = (float).5*dh1dh7[1][i+0]; 287 B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = (float).5*dh1dh7[0][i+0]; 288 B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2] = 0.; 289 B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+0] = (float).5*dh1dh7[2][i+0]; 290 B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1] = 0.; 291 B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2] = (float).5*dh1dh7[0][i+0]; 292 B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+0] = 0.; 293 B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1] = (float).5*dh1dh7[2][i+0]; 294 B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2] = (float).5*dh1dh7[1][i+0]; 295 B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+0] = 0.; 296 B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1] = 0.; 297 B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2] = 0.; 298 B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+0] = dh1dh7[0][i+0]; 299 B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1] = dh1dh7[1][i+0]; 300 B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2] = dh1dh7[2][i+0]; 301 301 } 302 302 … … 338 338 /*Get dh1dh7 in actual coordinate system: */ 339 339 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss); 340 GetNodalFunctionsP1( l1l6, gauss);340 GetNodalFunctionsP1(&l1l6[0], gauss); 341 341 342 342 /*Build B: */ 343 343 for (i=0;i<NUMNODESP1;i++){ 344 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh6[0][i]; //B[0][NDOF4*i]=dh1dh6[0][i]; 345 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0.; 346 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0.; 347 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0.; 348 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh6[1][i]; 349 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0.; 350 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0.; 351 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0.; 352 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=dh1dh6[2][i]; 353 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=.5*dh1dh6[1][i]; 354 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=.5*dh1dh6[0][i]; 355 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0.; 356 *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i)=.5*dh1dh6[2][i]; 357 *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1)=0.; 358 *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2)=.5*dh1dh6[0][i]; 359 *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i)=0.; 360 *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1)=.5*dh1dh6[2][i]; 361 *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2)=.5*dh1dh6[1][i]; 362 *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i)=0.; 363 *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1)=0.; 364 *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2)=0.; 365 *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i)=dh1dh6[0][i]; 366 *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1)=dh1dh6[1][i]; 367 *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2)=dh1dh6[2][i]; 344 *(B+(NDOF4*NUMNODESP1)*0+NDOF4*i)=dh1dh6[0][i]; //B[0][NDOF4*i]=dh1dh6[0][i]; 345 *(B+(NDOF4*NUMNODESP1)*0+NDOF4*i+1)=0.; 346 *(B+(NDOF4*NUMNODESP1)*0+NDOF4*i+2)=0.; 347 *(B+(NDOF4*NUMNODESP1)*1+NDOF4*i)=0.; 348 *(B+(NDOF4*NUMNODESP1)*1+NDOF4*i+1)=dh1dh6[1][i]; 349 *(B+(NDOF4*NUMNODESP1)*1+NDOF4*i+2)=0.; 350 *(B+(NDOF4*NUMNODESP1)*2+NDOF4*i)=0.; 351 *(B+(NDOF4*NUMNODESP1)*2+NDOF4*i+1)=0.; 352 *(B+(NDOF4*NUMNODESP1)*2+NDOF4*i+2)=dh1dh6[2][i]; 353 *(B+(NDOF4*NUMNODESP1)*3+NDOF4*i)=.5*dh1dh6[1][i]; 354 *(B+(NDOF4*NUMNODESP1)*3+NDOF4*i+1)=.5*dh1dh6[0][i]; 355 *(B+(NDOF4*NUMNODESP1)*3+NDOF4*i+2)=0.; 356 *(B+(NDOF4*NUMNODESP1)*4+NDOF4*i)=.5*dh1dh6[2][i]; 357 *(B+(NDOF4*NUMNODESP1)*4+NDOF4*i+1)=0.; 358 *(B+(NDOF4*NUMNODESP1)*4+NDOF4*i+2)=.5*dh1dh6[0][i]; 359 *(B+(NDOF4*NUMNODESP1)*5+NDOF4*i)=0.; 360 *(B+(NDOF4*NUMNODESP1)*5+NDOF4*i+1)=.5*dh1dh6[2][i]; 361 *(B+(NDOF4*NUMNODESP1)*5+NDOF4*i+2)=.5*dh1dh6[1][i]; 362 *(B+(NDOF4*NUMNODESP1)*6+NDOF4*i)=0.; 363 *(B+(NDOF4*NUMNODESP1)*6+NDOF4*i+1)=0.; 364 *(B+(NDOF4*NUMNODESP1)*6+NDOF4*i+2)=0.; 365 *(B+(NDOF4*NUMNODESP1)*7+NDOF4*i)=dh1dh6[0][i]; 366 *(B+(NDOF4*NUMNODESP1)*7+NDOF4*i+1)=dh1dh6[1][i]; 367 *(B+(NDOF4*NUMNODESP1)*7+NDOF4*i+2)=dh1dh6[2][i]; 368 _assert_((NDOF4*NUMNODESP1)*7+NDOF4*i+2<8*24); 368 369 } 369 370 370 371 for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function 371 *(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0.; 372 *(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0.; 373 *(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0.; 374 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0.; 375 *(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0.; 376 *(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0.; 377 *(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=l1l6[i]; 378 *(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=0.; 372 B[(NDOF4*NUMNODESP1)*0+NDOF4*i+3]=0.; 373 B[(NDOF4*NUMNODESP1)*1+NDOF4*i+3]=0.; 374 B[(NDOF4*NUMNODESP1)*2+NDOF4*i+3]=0.; 375 B[(NDOF4*NUMNODESP1)*3+NDOF4*i+3]=0.; 376 B[(NDOF4*NUMNODESP1)*4+NDOF4*i+3]=0.; 377 B[(NDOF4*NUMNODESP1)*5+NDOF4*i+3]=0.; 378 B[(NDOF4*NUMNODESP1)*6+NDOF4*i+3]=l1l6[i]; 379 B[(NDOF4*NUMNODESP1)*7+NDOF4*i+3]=0.; 380 _assert_(((NDOF4*NUMNODESP1)*7+NDOF4*i+3)<8*24); 379 381 } 380 382 … … 476 478 /*B_primeuild B_prime: */ 477 479 for (i=0;i<NUMNODESP1;i++){ 478 *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh6[0][i]; //B_prime[0][NDOF4*i]=dh1dh6[0][i]; 479 *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0.; 480 *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0.; 481 *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0.; 482 *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh6[1][i]; 483 *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0.; 484 *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0.; 485 *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0.; 486 *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=dh1dh6[2][i]; 487 *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=dh1dh6[1][i]; 488 *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=dh1dh6[0][i]; 489 *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0.; 490 *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i)=dh1dh6[2][i]; 491 *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1)=0.; 492 *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2)=dh1dh6[0][i]; 493 *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i)=0.; 494 *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1)=dh1dh6[2][i]; 495 *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2)=dh1dh6[1][i]; 496 *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i)=dh1dh6[0][i]; 497 *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1)=dh1dh6[1][i]; 498 *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2)=dh1dh6[2][i]; 499 *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i)=0.; 500 *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1)=0.; 501 *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2)=0.; 480 *(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i)=dh1dh6[0][i]; //B_prime[0][NDOF4*i]=dh1dh6[0][i]; 481 *(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i+1)=0.; 482 *(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i+2)=0.; 483 *(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i)=0.; 484 *(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i+1)=dh1dh6[1][i]; 485 *(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i+2)=0.; 486 *(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i)=0.; 487 *(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i+1)=0.; 488 *(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i+2)=dh1dh6[2][i]; 489 *(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i)=dh1dh6[1][i]; 490 *(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i+1)=dh1dh6[0][i]; 491 *(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i+2)=0.; 492 *(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i)=dh1dh6[2][i]; 493 *(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i+1)=0.; 494 *(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i+2)=dh1dh6[0][i]; 495 *(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i)=0.; 496 *(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i+1)=dh1dh6[2][i]; 497 *(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i+2)=dh1dh6[1][i]; 498 *(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i)=dh1dh6[0][i]; 499 *(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i+1)=dh1dh6[1][i]; 500 *(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i+2)=dh1dh6[2][i]; 501 *(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i)=0.; 502 *(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i+1)=0.; 503 *(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i+2)=0.; 504 _assert_((NDOF4*NUMNODESP1)*7+NDOF4*i+2<8*24); 502 505 } 503 506 504 507 for (i=0;i<NUMNODESP1;i++){ //last column 505 *(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0.; 506 *(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0.; 507 *(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0.; 508 *(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0.; 509 *(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0.; 510 *(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0.; 511 *(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=0.; 512 *(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=l1l6[i]; 508 *(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i+3)=0.; 509 *(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i+3)=0.; 510 *(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i+3)=0.; 511 *(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i+3)=0.; 512 *(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i+3)=0.; 513 *(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i+3)=0.; 514 *(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i+3)=0.; 515 *(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i+3)=l1l6[i]; 516 _assert_((NDOF4*NUMNODESP1)*7+NDOF4*i+2<8*24); 513 517 } 514 518 … … 1087 1091 if(*Jdet<0) _error_("negative jacobian determinant!"); 1088 1092 } 1089 /*}}}* /1093 /*}}}* 1090 1094 /*FUNCTION PentaRef::GetSegmentJacobianDeterminant{{{*/ 1091 1095 void PentaRef::GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussPenta* gauss){
Note:
See TracChangeset
for help on using the changeset viewer.
