Changeset 15473


Ignore:
Timestamp:
07/09/13 16:23:42 (12 years ago)
Author:
seroussi
Message:

NEW: implementing Franca et al stabilization

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r15467 r15473  
    23082308
    23092309                /*Compute the length of this edge and compare it to the minimal length*/
    2310                 length=pow(pow(xyz_list[node0][0]-xyz_list[node1][0],2.0)+pow(xyz_list[node0][1]-xyz_list[node1][1],2.0)+pow(xyz_list[node0][2]-xyz_list[node1][2],2.0),0.5);
     2310                length=sqrt(pow(xyz_list[node0][0]-xyz_list[node1][0],2)+pow(xyz_list[node0][1]-xyz_list[node1][1],2)+pow(xyz_list[node0][2]-xyz_list[node1][2],2));
    23112311                if(length<minlength || minlength<0) minlength=length;
    23122312        }
    … …  
    68866886        /*Condensation*/
    68876887        ReduceMatrixStokes(Ke->values, &Ke_temp[0][0]);
     6888        //Ke->Echo();
     6889        //_error_("stop");
    68886890
    68896891        /*Transform Coordinate System*/
    … …  
    69156917        GaussPenta *gauss=NULL;
    69166918
     6919        /*Stabilization*/
     6920        bool       stabilization = false;
     6921        IssmDouble dbasis[3][6];
     6922        IssmDouble dmu[3];
     6923        IssmDouble mu;
     6924        IssmDouble dnodalbasis[6][6][3];
     6925        IssmDouble SU[6][4][4];
     6926        IssmDouble SW[6][4][4];
     6927        int p,q,ii;
     6928        int c=3; //index of pressure
     6929
    69176930        /*If on water or not Stokes, skip stiffness: */
    69186931        inputs->GetInputValue(&approximation,ApproximationEnum);
    … …  
    69326945
    69336946        /* Start  looping on the number of gaussian points: */
    6934         gauss=new GaussPenta(5,5);
     6947        gauss=new GaussPenta(3,2);
    69356948        for(int ig=gauss->begin();ig<gauss->end();ig++){
    69366949
    … …  
    69516964                                        &D[0][0],8,8,0,
    69526965                                        &B_prime[0][0],8,24,0,
    6953                                         &Ke_temp[0][0],1);
     6966                                        &Ke_temp[0][0],0);
    69546967
    69556968                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_temp[i][j];
    69566969
    69576970                /*Add stabilization*/
    6958                 D_scalar_stab=-gauss->weight*Jdet*1./3.*pow(diameter,2.)/(4.*1000*rigidity)*stokesreconditioning;
    6959                 GetBConduct(&B_stab[0][0],&xyz_list[0][0],gauss);
    6960 
    6961                 D_stab[0][0]=D_scalar_stab; D_stab[0][1]=0;             D_stab[0][2]=0;
    6962                 D_stab[1][0]=0;             D_stab[1][1]=D_scalar_stab; D_stab[1][2]=0;
    6963                 D_stab[2][0]=0;             D_stab[2][1]=0;             D_stab[2][2]=D_scalar_stab;
    6964 
    6965                 TripleMultiply(&B_stab[0][0],3,6,1,
    6966                                         &D_stab[0][0],3,3,0,
    6967                                         &B_stab[0][0],3,6,0,
    6968                                         &Ke_temp_stab[0][0],1);
    6969 
    6970                 for(i=0;i<NUMVERTICES;i++) for(j=0;j<NUMVERTICES;j++) Ke->values[numdof*(i*4+3)+j*4+3]+=Ke_temp_stab[i][j];
    6971 
    6972         }
    6973 
     6971                if(stabilization){
     6972                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
     6973                        dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
     6974                        mu = viscosity;
     6975                        for(i=0;i<6;i++){
     6976                                for(p=0;p<6;p++){
     6977                                        for(j=0;j<3;j++){
     6978                                                dnodalbasis[i][p][j] = dbasis[j][p];
     6979                                        }
     6980                                }
     6981                        }
     6982                        for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
     6983                                SU[p][i][j]=0.;
     6984                                SW[p][i][j]=0.;
     6985                        }
     6986                        for(p=0;p<6;p++){
     6987                                for(i=0;i<3;i++){
     6988                                        SU[p][i][c] += dbasis[i][p];
     6989                                        SW[p][c][i] += dbasis[i][p];
     6990                                        for(j=0;j<3;j++){
     6991                                                SU[p][i][i] += -dmu[j]*dbasis[j][p];
     6992                                                SU[p][i][j] += -dmu[i]*dbasis[i][p];
     6993                                                for(ii=0;ii<6;ii++){
     6994                                                        SU[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
     6995                                                        SU[p][i][j] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
     6996                                                }
     6997                                                SW[p][i][i] += -dmu[j]*dbasis[j][p];
     6998                                                SW[p][j][i] += -dmu[j]*dbasis[i][p];
     6999                                                for(ii=0;ii<6;ii++){
     7000                                                        SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
     7001                                                        SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
     7002                                                }
     7003                                        }
     7004                                }
     7005                        }
     7006                        IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity);
     7007                        for(p=0;p<6;p++){
     7008                                for(q=0;q<6;q++){
     7009                                        for(i=0;i<4;i++){
     7010                                                for(j=0;j<4;j++){
     7011                                                        Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][0]*SU[q][0][j];
     7012                                                        Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][1]*SU[q][1][j];
     7013                                                        Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][2]*SU[q][2][j];
     7014                                                }
     7015                                        }
     7016                                }
     7017                        }
     7018
     7019                }
     7020
     7021
     7022                //viscosity = 1000*rigidity;
     7023                //D_scalar_stab=-gauss->weight*Jdet*1./3.*pow(diameter,2)/(8.*2.*viscosity)*stokesreconditioning;
     7024                //if(id==24){
     7025                //printf("tau %g\n",1./3.*pow(diameter,2)/(8.*2.*viscosity));
     7026                //}
     7027                //GetBConduct(&B_stab[0][0],&xyz_list[0][0],gauss);
     7028
     7029                //D_stab[0][0]=D_scalar_stab; D_stab[0][1]=0;             D_stab[0][2]=0;
     7030                //D_stab[1][0]=0;             D_stab[1][1]=D_scalar_stab; D_stab[1][2]=0;
     7031                //D_stab[2][0]=0;             D_stab[2][1]=0;             D_stab[2][2]=D_scalar_stab;
     7032
     7033                //TripleMultiply(&B_stab[0][0],3,6,1,
     7034                //                      &D_stab[0][0],3,3,0,
     7035                //                      &B_stab[0][0],3,6,0,
     7036                //                      &Ke_temp_stab[0][0],0);
     7037
     7038                //for(i=0;i<NUMVERTICES;i++) for(j=0;j<NUMVERTICES;j++) Ke->values[numdof*(i*4+3)+j*4+3]+=Ke_temp_stab[i][j];
     7039
     7040        }
     7041
     7042//      if(id==24){
     7043//              _error_("");
     7044//      }
     7045
     7046        //Ke->Echo();
     7047        //_error_("stop");
    69747048        /*Transform Coordinate System*/
    69757049        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
     7050        if(id==24){
     7051                //printarray(Ke->values,24,24);
     7052        }
    69767053
    69777054        /*Clean up and return*/
    … …  
    77907867                for(i=0;i<NUMVERTICES+1;i++){
    77917868                        Pe_gaussian[i*NDOF4+2]+=-rho_ice*gravity*Jdet*gauss->weight*l1l7[i];
    7792                         Pe_gaussian[i*NDOF4+0]+=forcex*Jdet*gauss->weight*l1l7[i];
    7793                         Pe_gaussian[i*NDOF4+1]+=forcey*Jdet*gauss->weight*l1l7[i];
    7794                         Pe_gaussian[i*NDOF4+2]+=forcez*Jdet*gauss->weight*l1l7[i];
     7869                        Pe_gaussian[i*NDOF4+0]+=rho_ice*forcex*Jdet*gauss->weight*l1l7[i];
     7870                        Pe_gaussian[i*NDOF4+1]+=rho_ice*forcey*Jdet*gauss->weight*l1l7[i];
     7871                        Pe_gaussian[i*NDOF4+2]+=rho_ice*forcez*Jdet*gauss->weight*l1l7[i];
    77957872                }
    77967873
    … …  
    78287905        int        i,j;
    78297906        int        approximation;
    7830         IssmDouble Jdet,gravity,rho_ice,B,D_scalar_stab;
     7907        IssmDouble Jdet,gravity,rho_ice,B,D_scalar_stab,viscosity;
    78317908        IssmDouble forcex,forcey,forcez,diameter,stokesreconditioning;
    78327909        IssmDouble xyz_list[NUMVERTICES][3];
    … …  
    78347911        IssmDouble l1l6[6]; //for the six nodes and the bubble
    78357912        IssmDouble dh1dh6[3][NUMVERTICES];
    7836         IssmDouble Pe_gaussian[numdof]={0.0}; //for the six nodes and the bubble
    78377913        GaussPenta *gauss=NULL;
     7914
     7915        /*Stabilization*/
     7916        bool       stabilization = false;
     7917        IssmDouble dbasis[3][6];
     7918        IssmDouble dmu[3];
     7919        IssmDouble mu;
     7920        IssmDouble dnodalbasis[6][6][3];
     7921        IssmDouble SW[6][4][4];
     7922        int p,q,ii;
     7923        int c=3; //index of pressure
    78387924
    78397925        /*Initialize Element vector and return if necessary*/
    … …  
    78517937        Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
    78527938        Input* loadingforcez_input=inputs->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
     7939        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     7940        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     7941        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    78537942
    78547943        /*Find minimal length*/
    … …  
    78567945
    78577946        /* Start  looping on the number of gaussian points: */
    7858         gauss=new GaussPenta(5,5);
     7947        gauss=new GaussPenta(2,3);
    78597948        for(int ig=gauss->begin();ig<gauss->end();ig++){
    78607949
    … …  
    78637952                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    78647953                GetNodalFunctionsP1(&l1l6[0], gauss);
     7954                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     7955                material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    78657956
    78667957                loadingforcex_input->GetInputValue(&forcex, gauss);
    … …  
    78697960
    78707961                for(i=0;i<NUMVERTICES;i++){
    7871                         Pe_gaussian[i*NDOF4+2]+=-rho_ice*gravity*Jdet*gauss->weight*l1l6[i];
    7872                         Pe_gaussian[i*NDOF4+0]+=forcex*Jdet*gauss->weight*l1l6[i];
    7873                         Pe_gaussian[i*NDOF4+1]+=forcey*Jdet*gauss->weight*l1l6[i];
    7874                         Pe_gaussian[i*NDOF4+2]+=forcez*Jdet*gauss->weight*l1l6[i];
     7962                        pe->values[i*NDOF4+2]+=-rho_ice*gravity*Jdet*gauss->weight*l1l6[i];
     7963                        pe->values[i*NDOF4+0]+=rho_ice*forcex*Jdet*gauss->weight*l1l6[i];
     7964                        pe->values[i*NDOF4+1]+=rho_ice*forcey*Jdet*gauss->weight*l1l6[i];
     7965                        pe->values[i*NDOF4+2]+=rho_ice*forcez*Jdet*gauss->weight*l1l6[i];
    78757966                }
    78767967
    78777968                /*Add stabilization*/
    7878                 D_scalar_stab=-gauss->weight*Jdet*1./3.*pow(diameter,2.)/(4.*1000*B)*stokesreconditioning;
    7879                 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0],gauss);
    7880 
    7881                 for(i=0;i<NUMVERTICES;i++){
    7882                         Pe_gaussian[i*NDOF4+3]+=-rho_ice*gravity*D_scalar_stab*dh1dh6[2][i];
    7883                         Pe_gaussian[i*NDOF4+3]+=forcex*D_scalar_stab*dh1dh6[0][i];
    7884                         Pe_gaussian[i*NDOF4+3]+=forcey*D_scalar_stab*dh1dh6[1][i];
    7885                         Pe_gaussian[i*NDOF4+3]+=forcez*D_scalar_stab*dh1dh6[2][i];
    7886                 }
    7887         }
    7888 
    7889         for(i=0;i<numdof;i++) pe->values[i]+=Pe_gaussian[i];
     7969                if(stabilization){
     7970                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
     7971                        dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
     7972                        mu = viscosity;
     7973                        for(i=0;i<6;i++){
     7974                                for(p=0;p<6;p++){
     7975                                        for(j=0;j<3;j++){
     7976                                                dnodalbasis[i][p][j] = dbasis[j][p];
     7977                                        }
     7978                                }
     7979                        }
     7980                        for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
     7981                                SW[p][i][j]=0.;
     7982                        }
     7983                        for(p=0;p<6;p++){
     7984                                for(i=0;i<3;i++){
     7985                                        SW[p][c][i] += dbasis[i][p];
     7986                                        for(j=0;j<3;j++){
     7987                                                SW[p][i][i] += -dmu[j]*dbasis[j][p];
     7988                                                SW[p][j][i] += -dmu[j]*dbasis[i][p];
     7989                                                for(ii=0;ii<6;ii++){
     7990                                                        SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
     7991                                                        SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
     7992                                                }
     7993                                        }
     7994                                }
     7995                        }
     7996                        IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity);
     7997                        for(p=0;p<6;p++){
     7998                                for(j=0;j<4;j++){
     7999                                        pe->values[p*4+j] += gauss->weight*Jdet*tau*forcex*SW[p][j][0];
     8000                                        pe->values[p*4+j] += gauss->weight*Jdet*tau*forcey*SW[p][j][1];
     8001                                        pe->values[p*4+j] += gauss->weight*Jdet*tau*(forcez-rho_ice*gravity)*SW[p][j][2];
     8002                                }
     8003                        }
     8004
     8005                }
     8006                ///*Add stabilization*/
     8007                //D_scalar_stab=-gauss->weight*Jdet*1./3.*pow(diameter,2)/(8.*2.*viscosity)*stokesreconditioning;
     8008                //GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0],gauss);
     8009
     8010                //for(i=0;i<NUMVERTICES;i++){
     8011                //      Pe_gaussian[i*NDOF4+3]+=-rho_ice*gravity*D_scalar_stab*dh1dh6[2][i];
     8012                //      Pe_gaussian[i*NDOF4+3]+=forcex*D_scalar_stab*dh1dh6[0][i];
     8013                //      Pe_gaussian[i*NDOF4+3]+=forcey*D_scalar_stab*dh1dh6[1][i];
     8014                //      Pe_gaussian[i*NDOF4+3]+=forcez*D_scalar_stab*dh1dh6[2][i];
     8015                //}
     8016        }
    78908017
    78918018        /*Transform coordinate system*/
    78928019        TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
     8020
     8021        if(id==24){
     8022                printarray(pe->values,24,1);
     8023                _error_("STOP");
     8024        }
    78938025
    78948026        /*Clean up and return*/
Note: See TracChangeset for help on using the changeset viewer.