Changeset 9217


Ignore:
Timestamp:
08/09/11 11:30:31 (14 years ago)
Author:
Mathieu Morlighem
Message:

some cleanup up in elements: removed almost all Ke_gg_gaussian

Location:
issm/trunk/src/c/objects/Elements
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r9206 r9217  
    783783        ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
    784784        ElementMatrix* Ke2=new ElementMatrix(this->nodes     ,NUMVERTICES,this->parameters,PattynApproximationEnum);
    785         ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
     785        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    786786        delete Ke1; delete Ke2;
    787787
    … …  
    14511451        double     B[5][numdof];
    14521452        double     Bprime[5][numdof];
    1453         double     Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    14541453        Tria*      tria=NULL;
    14551454        GaussPenta *gauss=NULL;
    … …  
    14891488                                        &D[0][0],5,5,0,
    14901489                                        &Bprime[0][0],5,numdof,0,
    1491                                         &Ke_gg_gaussian[0][0],0);
    1492 
    1493                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j];
     1490                                        &Ke->values[0],1);
    14941491        }
    14951492
    … …  
    15171514        double    DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
    15181515        double    DL_scalar;
    1519         double    Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
    15201516        Friction  *friction = NULL;
    15211517        GaussPenta *gauss=NULL;
    … …  
    15651561                                        &DL[0][0],2,2,0,
    15661562                                        &L[0][0],2,numdof,0,
    1567                                         &Ke_gg_gaussian[0][0],0);
    1568 
    1569                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j];
     1563                                        &Ke->values[0],1);
    15701564        }
    15711565
    … …  
    17801774        double      Bprime[NDOF1][NUMVERTICES];
    17811775        double      DL_scalar;
    1782         double      Ke_gg[numdof][numdof]={0.0};
    17831776        GaussPenta  *gauss=NULL;
    17841777
    … …  
    18041797                                        &DL_scalar,1,1,0,
    18051798                                        &Bprime[0][0],1,NUMVERTICES,0,
    1806                                         &Ke_gg[0][0],0);
    1807 
    1808                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg[i][j];
     1799                                        &Ke->values[0],1);
    18091800        }
    18101801
    … …  
    18301821        double    basis[NUMVERTICES];
    18311822        GaussPenta *gauss=NULL;
    1832         double Ke_g[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    18331823
    18341824        /*Initialize Element matrix*/
    … …  
    18491839                GetNodalFunctionsP1(&basis[0], gauss);
    18501840
    1851                 /**********************Do not forget the sign**********************************/
    1852                 DL_scalar=- gauss->weight*Jdet2d*surface_normal[2];
    1853                 /******************************************************************************/
     1841                DL_scalar= - gauss->weight*Jdet2d*surface_normal[2];
    18541842
    18551843                TripleMultiply( basis,1,numdof,1,
    18561844                                        &DL_scalar,1,1,0,
    18571845                                        basis,1,numdof,0,
    1858                                         &Ke_g[0][0],0);
    1859 
    1860                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
     1846                                        &Ke->values[0],1);
    18611847        }
    18621848
    … …  
    19081894        double     D[3][3];
    19091895        double     K[2][2]={0.0};
    1910         double     Ke_gaussian_conduct[numdof][numdof];
    1911         double     Ke_gaussian_advec[numdof][numdof];
    1912         double     Ke_gaussian_artdiff[numdof][numdof];
    1913         double     Ke_gaussian_transient[numdof][numdof];
    19141896        Tria*      tria=NULL;
    19151897        GaussPenta *gauss=NULL;
    … …  
    19651947                                        &D[0][0],3,3,0,
    19661948                                        &B_conduct[0][0],3,numdof,0,
    1967                                         &Ke_gaussian_conduct[0][0],0);
     1949                                        &Ke->values[0],1);
    19681950
    19691951                /*Advection: */
    … …  
    19891971                                        &D[0][0],3,3,0,
    19901972                                        &Bprime_advec[0][0],3,numdof,0,
    1991                                         &Ke_gaussian_advec[0][0],0);
     1973                                        &Ke->values[0],1);
    19921974
    19931975                /*Transient: */
    … …  
    20011983                                                &D_scalar_trans,1,1,0,
    20021984                                                &L[0],1,numdof,0,
    2003                                                 &Ke_gaussian_transient[0][0],0);
    2004                 }
    2005                 else{
    2006                         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gaussian_transient[i][j]=0;
     1985                                                &Ke->values[0],1);
    20071986                }
    20081987
    … …  
    20212000                                                &K[0][0],2,2,0,
    20222001                                                &B_artdiff[0][0],2,numdof,0,
    2023                                                 &Ke_gaussian_artdiff[0][0],0);
     2002                                                &Ke->values[0],1);
    20242003                }
    20252004                else if(artdiff==2){
    … …  
    20312010                        for(i=0;i<numdof;i++){
    20322011                                for(j=0;j<numdof;j++){
    2033                                         Ke_gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*
     2012                                        Ke->values[i*numdof+j]=tau_parameter*D_scalar_advec*
    20342013                                          ((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i])*((u-um)*dbasis[0][j]+(v-vm)*dbasis[1][j]+(w-wm)*dbasis[2][j]);
    20352014                                }
    … …  
    20382017                                for(i=0;i<numdof;i++){
    20392018                                        for(j=0;j<numdof;j++){
    2040                                                 Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]);
     2019                                                Ke->values[i*numdof+j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]);
    20412020                                        }
    20422021                                }
    20432022                        }
    20442023                }
    2045                 else{
    2046                         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gaussian_artdiff[i][j]=0;
    2047                 }
    2048 
    2049                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j];
    20502024        }
    20512025
    … …  
    20702044        double    basis[NUMVERTICES];
    20712045        double    D_scalar;
    2072         double    Ke_gaussian[numdof][numdof]={0.0};
    20732046        GaussPenta *gauss=NULL;
    20742047
    … …  
    21022075                                        &D_scalar,1,1,0,
    21032076                                        &basis[0],1,numdof,0,
    2104                                         &Ke_gaussian[0][0],0);
    2105 
    2106                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian[i][j];
     2077                                        &Ke->values[0],1);
    21072078        }
    21082079       
    … …  
    21982169        double     D[3][3];
    21992170        double     K[2][2]={0.0};
    2200         double     Ke_gaussian_conduct[numdof][numdof];
    2201         double     Ke_gaussian_advec[numdof][numdof];
    2202         double     Ke_gaussian_artdiff[numdof][numdof];
    2203         double     Ke_gaussian_transient[numdof][numdof];
    22042171        Tria*      tria=NULL;
    22052172        GaussPenta *gauss=NULL;
    … …  
    22482215                                        &D[0][0],3,3,0,
    22492216                                        &B_conduct[0][0],3,numdof,0,
    2250                                         &Ke_gaussian_conduct[0][0],0);
     2217                                        &Ke->values[0],1);
    22512218
    22522219                /*Advection: */
    … …  
    22722239                                        &D[0][0],3,3,0,
    22732240                                        &Bprime_advec[0][0],3,numdof,0,
    2274                                         &Ke_gaussian_advec[0][0],0);
     2241                                        &Ke->values[0],1);
    22752242
    22762243                /*Transient: */
    … …  
    22842251                                                &D_scalar_trans,1,1,0,
    22852252                                                &L[0],1,numdof,0,
    2286                                                 &Ke_gaussian_transient[0][0],0);
    2287                 }
    2288                 else{
    2289                         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gaussian_transient[i][j]=0;
     2253                                                &Ke->values[0],1);
    22902254                }
    22912255
    … …  
    23042268                                                &K[0][0],2,2,0,
    23052269                                                &B_artdiff[0][0],2,numdof,0,
    2306                                                 &Ke_gaussian_artdiff[0][0],0);
     2270                                                &Ke->values[0],1);
    23072271                }
    23082272                else if(artdiff==2){
    … …  
    23142278                        for(i=0;i<numdof;i++){
    23152279                                for(j=0;j<numdof;j++){
    2316                                         Ke_gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*
     2280                                        Ke->values[i*numdof+j]=tau_parameter*D_scalar_advec*
    23172281                                          ((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i])*((u-um)*dbasis[0][j]+(v-vm)*dbasis[1][j]+(w-wm)*dbasis[2][j]);
    23182282                                }
    … …  
    23212285                                for(i=0;i<numdof;i++){
    23222286                                        for(j=0;j<numdof;j++){
    2323                                                 Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]);
     2287                                                Ke->values[i*numdof+j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]);
    23242288                                        }
    23252289                                }
    23262290                        }
    23272291                }
    2328                 else{
    2329                         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gaussian_artdiff[i][j]=0;
    2330                 }
    2331 
    2332                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j];
    23332292        }
    23342293
    … …  
    23542313        double    basis[NUMVERTICES];
    23552314        double    D_scalar;
    2356         double    Ke_gaussian[numdof][numdof]={0.0};
    23572315        GaussPenta *gauss=NULL;
    23582316
    … …  
    23862344                                        &D_scalar,1,1,0,
    23872345                                        &basis[0],1,numdof,0,
    2388                                         &Ke_gaussian[0][0],0);
    2389 
    2390                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian[i][j];
     2346                                        &Ke->values[0],1);
    23912347        }
    23922348       

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r9206 r9217  
    484484        double     DLprime[2][2]                    = {0.0};
    485485        double     DL_scalar;
    486         double     Ke_gg_gaussian[numdof][numdof]   = {0.0};
    487         double     Ke_gg_thickness1[numdof][numdof] = {0.0};
    488         double     Ke_gg_thickness2[numdof][numdof] = {0.0};
    489486        GaussTria *gauss                            = NULL;
    490487
    … …  
    548545                                        &DL[0][0],2,2,0,
    549546                                        &B[0][0],2,numdof,0,
    550                                         &Ke_gg_thickness1[0][0],0);
     547                                        &Ke->values[0],1);
    551548
    552549                TripleMultiply( &B[0][0],2,numdof,1,
    553550                                        &DLprime[0][0],2,2,0,
    554551                                        &Bprime[0][0],2,numdof,0,
    555                                         &Ke_gg_thickness2[0][0],0);
    556 
    557                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_thickness1[i][j];
    558                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_thickness2[i][j];
     552                                        &Ke->values[0],1);
    559553
    560554                if(artdiff){
    … …  
    565559                                                &KDL[0][0],2,2,0,
    566560                                                &Bprime[0][0],2,numdof,0,
    567                                                 &Ke_gg_gaussian[0][0],0);
    568 
    569                         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j];
     561                                                &Ke->values[0],1);
    570562                }
    571563        }
    … …  
    590582        double     DL[2][2]={0.0};
    591583        double     DL_scalar;
    592         double     Ke_gg[numdof][numdof]={0.0};
    593584        GaussTria  *gauss=NULL;
    594585
    … …  
    623614                                        &DL[0][0],2,2,0,
    624615                                        &Bprime[0][0],2,numdof,0,
    625                                         &Ke_gg[0][0],0);
    626 
    627                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg[i][j];
     616                                        &Ke->values[0],1);
    628617        }
    629618
    … …  
    655644        double  DLprime[2][2]={0.0};
    656645        double  DL_scalar;
    657         double  Ke_gg_gaussian[numdof][numdof]   = {0.0};
    658         double  Ke_gg_velocities[numdof][numdof] = {0.0};
    659646        GaussTria *gauss=NULL;
    660647
    … …  
    734721                                        &DLprime[0][0],2,2,0,
    735722                                        &Bprime[0][0],2,numdof,0,
    736                                         &Ke_gg_velocities[0][0],0);
    737 
    738                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_velocities[i][j];
     723                                        &Ke->values[0],1);
    739724
    740725                if(artdiff){
    … …  
    745730                                                &KDL[0][0],2,2,0,
    746731                                                &Bprime[0][0],2,numdof,0,
    747                                                 &Ke_gg_gaussian[0][0],0);
    748 
    749                         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j];
     732                                                &Ke->values[0],1);
    750733                }
    751734        }
    … …  
    786769        double     D[3][3]   = {0.0};
    787770        double     D_scalar;
    788         double     Ke_g[numdof][numdof];
    789771        GaussTria *gauss = NULL;
    790772
    … …  
    824806                                        &D[0][0],3,3,0,
    825807                                        &Bprime[0][0],3,numdof,0,
    826                                         &Ke_g[0][0],0);
    827 
    828                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
     808                                        &Ke->values[0],1);
    829809        }
    830810
    … …  
    849829        double     L[2][numdof];
    850830        double     DL[2][2]  = {{ 0,0 },{0,0}};
    851         double     Ke_g[numdof][numdof];
    852831        double     DL_scalar;
    853832        double     slope[2]  = {0.0,0.0};
    … …  
    894873                                        &DL[0][0],2,2,0,
    895874                                        &L[0][0],2,numdof,0,
    896                                         &Ke_g[0][0],0);
    897 
    898                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
     875                                        &Ke->values[0],1);
    899876        }
    900877
    … …  
    985962        double     dLy[NUMVERTICES];
    986963        double     K[2];
    987         double     Ke_gg_gaussian1[numdof][numdof]  ={0.0};
    988         double     Ke_gg_gaussian2[numdof][numdof]  ={0.0};
    989         double     Ke_gg_gaussian3[numdof][numdof]  ={0.0};
    990964        double     hydro_p,hydro_q,hydro_CR,hydro_n,rho_water,g,gamma;
    991965        double     DL_scalar1;
    … …  
    1019993                GetNodalFunctions(&L[0],gauss);
    1020994                GetNodalFunctionsDerivatives(&dL[0][0],&xyz_list[0][0],gauss);
    1021                 //printf("NUMVERTICES=%d, Jdettria=%f\n",NUMVERTICES,Jdettria); //BK
    1022995                for(i=0;i<NUMVERTICES;i++)dLx[i]=dL[0][i];
    1023996                for(i=0;i<NUMVERTICES;i++)dLy[i]=dL[1][i];
    1024                 /*    for(i=0;i<NUMVERTICES;i++)
    1025                                 { dLx[i]=dL[0][i];
    1026                                 printf("dLx[%d]=%f, L[%d]=%f\n",i,dLx[i]),i,L[i];}
    1027                                 for(i=0;i<NUMVERTICES;i++)
    1028                                 { dLy[i]=dL[1][i];
    1029                                 printf("dLy[%d]=%f\n",i,dLy[i]);} */
    1030997
    1031998                /*Get K parameter: */
    1032999                GetHydrologyK(&K[0],&xyz_list[0][0],gauss);
    1033                 //printf("K[0]=%g,K[1]=%g\n",K[0],K[1]);
    10341000
    10351001                if(dt){
    … …  
    10381004                                        &DL_scalar1,1,1,0,
    10391005                                        &L[0],1,numdof,0,
    1040                                         &Ke_gg_gaussian1[0][0],0);
     1006                                        &Ke->values[0],1);
    10411007                }
    10421008
    10431009                if(dt){
    1044                         DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]*dt; //don't forget the -
    1045                         DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]*dt; //don't forget the -
     1010                        DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]*dt;
     1011                        DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]*dt;
    10461012
    10471013                }
    10481014                else{
    1049                         DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]; //don't forget the -
    1050                         DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]; //don't forget the -
     1015                        DL_scalar2=-gauss->weight*Jdettria*gamma*K[0];
     1016                        DL_scalar3=-gauss->weight*Jdettria*gamma*K[1];
    10511017                }
    10521018               
    1053                 TripleMultiply( &dLx[0],1,numdof,1, &DL_scalar2,1,1,0, &L[0],1,numdof,0, &Ke_gg_gaussian2[0][0],0);
    1054                 TripleMultiply( &dLy[0],1,numdof,1, &DL_scalar3,1,1,0, &L[0],1,numdof,0, &Ke_gg_gaussian3[0][0],0);
    1055 
    1056                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian1[i][j];
    1057                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian2[i][j];
    1058                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian3[i][j];
    1059 
     1019                TripleMultiply(&dLx[0],1,numdof,1, &DL_scalar2,1,1,0, &L[0],1,numdof,0,&Ke->values[0],1);
     1020                TripleMultiply(&dLy[0],1,numdof,1, &DL_scalar3,1,1,0, &L[0],1,numdof,0,&Ke->values[0],1);
    10601021        }
    10611022
    … …  
    12101171        double     DLprime[2][2]={0.0};
    12111172        double     DL_scalar;
    1212         double     Ke_gg1[numdof][numdof]={0.0};
    1213         double     Ke_gg2[numdof][numdof]={0.0};
    12141173        GaussTria  *gauss=NULL;
    12151174
    … …  
    12491208                                        &DL_scalar,1,1,0,
    12501209                                        &L[0],1,numdof,0,
    1251                                         &Ke_gg1[0][0],0);
     1210                                        &Ke->values[0],1);
    12521211
    12531212                /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
    … …  
    12631222                                        &DLprime[0][0],2,2,0,
    12641223                                        &Bprime[0][0],2,numdof,0,
    1265                                         &Ke_gg2[0][0],0);
    1266 
    1267                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg1[i][j];
    1268                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg2[i][j];
     1224                                        &Ke->values[0],1);
    12691225        }
    12701226
    … …  
    12851241        double     xyz_list[NUMVERTICES][3];
    12861242        double     L[1][3];
    1287         double     Ke_g[numdof][numdof];
    12881243        GaussTria *gauss = NULL;
    12891244
    … …  
    13071262                                        &DL_scalar,1,1,0,
    13081263                                        &L[0][0],1,3,0,
    1309                                         &Ke_g[0][0],0);
    1310 
    1311                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
     1264                                        &Ke->values[0],1);
    13121265        }
    13131266
    … …  
    28682821        K[1]=pow(w,hydro_p)*pow((rho_ice*g*dsdy+(rho_water/rho_ice-1)*rho_ice*g*dbdy) - rho_ice * g * kn/w * (dsdy - dbdy ) * surface_slope,hydro_q);
    28692822       
    2870         //printf("K[0]=%g,K[1]=%g,w=%g\n",K[0],K[1],w); //bk
    2871         //if (w<0) {printf("negative w!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
    2872         //printf("K[0]=%g,K[1]=%g,w=%g\n",K[0],K[1],w);}
    28732823}
    28742824/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.