Changeset 442


Ignore:
Timestamp:
05/15/09 11:25:19 (16 years ago)
Author:
seroussi
Message:

fixed problems with stokes

Location:
issm/trunk/src/c/objects
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Icefront.cpp

    r416 r442  
    686686        QuadPressureLoadStokes(pe_g,matpar->GetRhoWater(),matpar->GetRhoIce(),matpar->GetG(),thickness_list_quad,bed_list_quad,normal1,normal2,normal3,normal4,&xyz_list_quad[0][0],fill);
    687687
    688 
    689688        #ifdef _DEBUGELEMENTS_
    690689        if(my_rank==RANK && eid==ELID){
    … …  
    702701                printf("fill %i\n",fill);
    703702                printf("pe_g->terms\n");
    704                 for(i=0;i<numgridsload*NDOF4;i++){
     703                for(i=0;i<numdofloads*NDOF4;i++){
    705704                        printf("%g ",*(pe_g->terms+i));
    706705                }
    … …  
    11801179        double J[4];
    11811180        double z_g[4];
    1182         double ice_pressure_tria[4];
    11831181        double air_pressure_tria;
    11841182        double water_level_above_g_tria;
    11851183        double water_pressure_tria;
    1186         double pressure_tria[4];
     1184        double pressure_tria;
    11871185
    11881186        /*To use tria functionalities: */
    … …  
    13311329                        //Loop on the grids of the quad
    13321330                        //Calculate the ice pressure
    1333                         for(j=0;j<4;j++){
    1334                                 ice_pressure_tria[j]=gravity*rho_ice*(surface_list[j]-z_g[i]);
    1335                         }
    13361331                        air_pressure_tria=0;
    13371332
    … …  
    13481343                        }
    13491344
    1350                         //Add pressure from triangle i
    1351                         //Loop on the grids of the quad
    1352                         for(j=0;j<4;j++){
    1353                                 pressure_tria[j] = ice_pressure_tria[j] + water_pressure_tria + air_pressure_tria;
    1354                         }
    1355 
    1356 
    1357                         pe_g[0]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*nx[i];
    1358                         pe_g[1]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*ny[i];
    1359                         pe_g[2]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*nz[i];
     1345                        //Add pressure
     1346                        pressure_tria = water_pressure_tria + air_pressure_tria;
     1347
     1348                        pe_g[0]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][0]*nx[i];
     1349                        pe_g[1]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][0]*ny[i];
     1350                        pe_g[2]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][0]*nz[i];
    13601351                        pe_g[3]+= 0;
    1361                         pe_g[4]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*nx[i];
    1362                         pe_g[5]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*ny[i];
    1363                         pe_g[6]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*nz[i];
     1352                        pe_g[4]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][1]*nx[i];
     1353                        pe_g[5]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][1]*ny[i];
     1354                        pe_g[6]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][1]*nz[i];
    13641355                        pe_g[7]+= 0;
    1365                         pe_g[8]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*nx[i];
    1366                         pe_g[9]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*ny[i];
    1367                         pe_g[10]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*nz[i];
     1356                        pe_g[8]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][2]*nx[i];
     1357                        pe_g[9]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][2]*ny[i];
     1358                        pe_g[10]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][2]*nz[i];
    13681359                        pe_g[11]+= 0;
    1369                         pe_g[12]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*nx[i];
    1370                         pe_g[13]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*ny[i];
    1371                         pe_g[14]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*nz[i];
     1360                        pe_g[12]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][3]*nx[i];
     1361                        pe_g[13]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][3]*ny[i];
     1362                        pe_g[14]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][3]*nz[i];
    13721363                        pe_g[15]+= 0;
    1373 
    1374                        
    13751364
    13761365                } //for(i=0;i<4;i++)

  • issm/trunk/src/c/objects/Matice.cpp

    r394 r442  
    333333        double eps0;
    334334
    335         eps0=10^-27;
     335        eps0=pow(10,-27);
    336336       
    337337        if (n==1){
    … …  
    362362                        else{
    363363                                e=(n-1)/2/n;
    364                        
    365364                                viscosity3d=2*B/(2*pow(A,e));
    366365                        }

  • issm/trunk/src/c/objects/Pengrid.cpp

    r394 r442  
    102102void  Pengrid::Demarshall(char** pmarshalled_dataset){
    103103
    104         int i;
    105104        char* marshalled_dataset=NULL;
    106105
    … …  
    207206        double slope[2];
    208207        int found=0;
    209         double Ke[4][4]={{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}};
     208        double Ke[4][4]={0.0};
    210209       
    211210        ParameterInputs* inputs=NULL;
    212        
    213211
    214212        /*recover pointers: */
    … …  
    229227        Ke[2][2]=kmax*pow(10,penalty_offset);
    230228       
    231        
    232229        /*Add Ke to global matrix Kgg: */
    233230        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
    … …  
    237234#define __FUNCT__ "Pengrid::PenaltyCreatePVector"
    238235void  Pengrid::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){
    239         throw ErrorException(__FUNCT__," not implemented yet!");
     236
     237        /*No penalty applied, do nothing: */
     238        return;
     239
    240240}
    241241

  • issm/trunk/src/c/objects/Penta.cpp

    r437 r442  
    679679                        /* Get Jacobian determinant: */
    680680                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
    681 
    682681                        DL_scalar=gauss_weight*Jdet;
    683682
    … …  
    746745
    747746        /*matrices: */
    748         double     Ke_temp[27][27]; //for the six nodes and the bubble
     747        double     Ke_temp[27][27]={0.0}; //for the six nodes and the bubble
    749748        double     Ke_reduced[numdof][numdof]; //for the six nodes only
    750749        double     Ke_gaussian[27][27];
    … …  
    756755        double     Jdet;
    757756        double     Jdet2d;
    758         double     D[8][8]={{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 }};         
     757        double     D[8][8]={0.0};
    759758        double     D_scalar;
    760759        double     tBD[numdof][8];
    761         double     DLStokes[14][14];
     760        double     DLStokes[14][14]={0.0};
    762761        double     tLDStokes[numdof2d][14];
    763762       
    … …  
    801800        gravity=matpar->GetG();
    802801
    803 
    804802        /* Get node coordinates and dof list: */
    805803        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    … …  
    808806        /*recover extra inputs from users, at current convergence iteration: */
    809807        inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
    810 
    811         /*Initialize Ke_temp to zero befor adding terms */
    812         for (i=0;i<27;i++){
    813                 for (j=0;j<27;j++){
    814                         Ke_temp[i][j]=0;
    815                 }
    816         }
    817808
    818809        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    … …  
    821812           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
    822813           points, same deal, which yields 3 gaussian points.*/
    823 
    824814
    825815        area_order=5;
    … …  
    841831                        gauss_coord[3]=*(vert_gauss_coord+igvert);
    842832
    843 
    844833                        /*Compute thickness: */
    845834                        GetParameterValue(&thickness, &h[0],gauss_coord);
    … …  
    868857                        }
    869858
    870 
    871859                        /*  Do the triple product tB*D*Bprime: */
    872860                        MatrixMultiply(&B[0][0],8,27,1,&D[0][0],8,8,0,&tBD[0][0],0);
    … …  
    883871
    884872        if((onbed==1) && (shelf==0)){
    885 
    886873
    887874                /*Build alpha2_list used by drag stiffness matrix*/
    … …  
    915902                }
    916903               
    917                
    918904                GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
    919905
    … …  
    974960                        for(i=0;i<numdof2d;i++){
    975961                                for(j=0;j<numdof2d;j++){
    976                                         K_terms[i][j]+=Ke_drag_gaussian[i][j];
     962                                        Ke_temp[i][j]+=Ke_drag_gaussian[i][j];
    977963                                }
    978964                        }
    … …  
    993979       
    994980        cleanup_and_return:
    995 
    996981
    997982        return;
    … …  
    24082393       
    24092394        #ifdef _DEBUG_
    2410         for (i=0;i<6;i++){
    2411                 printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i]);
     2395        for (i=0;i<7;i++){
     2396                printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i],dh1dh7_basic[2][i]);
    24122397        }
    24132398
    … …  
    27802765
    27812766        for (i=0;i<numgrids;i++){
    2782                 *(dh1dh7_basic+numgrids*0+i)=Jinv[0][0]*dh1dh7_param[0][i]+Jinv[1][0]*dh1dh7_param[1][i]+Jinv[2][0]*dh1dh7_param[2][i];
    2783                 *(dh1dh7_basic+numgrids*1+i)=Jinv[0][1]*dh1dh7_param[0][i]+Jinv[1][1]*dh1dh7_param[1][i]+Jinv[2][1]*dh1dh7_param[2][i];
    2784                 *(dh1dh7_basic+numgrids*2+i)=Jinv[0][2]*dh1dh7_param[0][i]+Jinv[1][2]*dh1dh7_param[1][i]+Jinv[2][2]*dh1dh7_param[2][i];
     2767                *(dh1dh7_basic+numgrids*0+i)=Jinv[0][0]*dh1dh7_param[0][i]+Jinv[0][1]*dh1dh7_param[1][i]+Jinv[0][2]*dh1dh7_param[2][i];
     2768                *(dh1dh7_basic+numgrids*1+i)=Jinv[1][0]*dh1dh7_param[0][i]+Jinv[1][1]*dh1dh7_param[1][i]+Jinv[1][2]*dh1dh7_param[2][i];
     2769                *(dh1dh7_basic+numgrids*2+i)=Jinv[2][0]*dh1dh7_param[0][i]+Jinv[2][1]*dh1dh7_param[1][i]+Jinv[2][2]*dh1dh7_param[2][i];
    27852770        }
    27862771
    … …  
    28902875
    28912876        /*matrices: */
    2892         double     Pe_temp[27]; //for the six nodes and the bubble
    2893         double     Pe_gaussian[27]; //for the six nodes and the bubble
    2894         double     Ke_temp[27][3]; //for the six nodes and the bubble
     2877        double     Pe_temp[27]={0.0}; //for the six nodes and the bubble
     2878        double     Pe_gaussian[27]={0.0}; //for the six nodes and the bubble
     2879        double     Ke_temp[27][3]={0.0}; //for the six nodes and the bubble
    28952880        double     Pe_reduced[numdof]; //for the six nodes only
    28962881        double     Ke_gaussian[27][3];
    … …  
    29022887        double     Jdet;
    29032888        double     Jdet2d;
    2904         double     D[8][8]={{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 }};         
     2889        double     D[8][8]={0.0};
    29052890        double     D_scalar;
    29062891        double     tBD[27][8];
    … …  
    29302915        GetDofList(&doflist[0],&numberofdofspernode);
    29312916
    2932         /*Initialize Ke_temp to zero befor adding terms */
    2933         for (i=0;i<27;i++){
    2934                 Pe_temp[i]=0;
    2935                 for (j=0;j<3;j++){
    2936                         Ke_temp[i][j]=0;
    2937                 }
    2938         }
    2939 
    2940        
    2941        
    29422917        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    29432918           get tria gaussian points as well as segment gaussian points. For tria gaussian
    … …  
    29452920           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
    29462921           points, same deal, which yields 3 gaussian points.*/
    2947 
    2948        
    29492922
    29502923        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    … …  
    30172990        }
    30182991
    3019 
    30202992        /*Deal with 2d friction at the bedrock interface: */
    30212993        if ( (onbed==1) && (shelf==1)){
    … …  
    30403012                        gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
    30413013                        gauss_coord_tria[2]=*(third_gauss_area_coord+igarea);
    3042                        
    30433014
    30443015                        /*Get the Jacobian determinant */
    … …  
    30763047        }
    30773048
    3078 
    30793049        /*Add P_terms to global vector pg: */
    30803050        VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
Note: See TracChangeset for help on using the changeset viewer.