Changeset 5821


Ignore:
Timestamp:
09/15/10 10:18:38 (15 years ago)
Author:
Mathieu Morlighem
Message:

some cleaning

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

Legend:

Unmodified
Added
Removed

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

    r5803 r5821  
    12751275
    12761276        int    i;
    1277         bool   found = false;
    12781277        Input *input = NULL;
    12791278
    … …  
    23952394void  Tria::CreateKMatrixBalancedthickness_CG(Mat Kgg){
    23962395
    2397         /* local declarations */
    2398         int             i,j;
    2399 
    2400         /* node data: */
     2396        /*Constants*/
    24012397        const int    numdof=NDOF1*NUMVERTICES;
    2402         double       xyz_list[NUMVERTICES][3];
    2403         int*         doflist=NULL;
    2404 
    2405         /* gaussian points: */
    2406         int     ig;
    2407         GaussTria *gauss=NULL;
    2408 
    2409         /* matrices: */
    2410         double L[NUMVERTICES];
    2411         double B[2][NUMVERTICES];
    2412         double Bprime[2][NUMVERTICES];
    2413         double DL[2][2]={0.0};
    2414         double DLprime[2][2]={0.0};
    2415         double DL_scalar;
    2416         double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix
    2417         double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    2418         double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    2419         double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    2420 
    2421         double Jdettria;
    2422 
    2423         /*input parameters for structural analysis (diagnostic): */
    2424         double  dvx[2];
    2425         double  dvy[2];
    2426         double  vx,vy;
    2427         double  dvxdx,dvydy;
    2428         double  v_gauss[2]={0.0};
    2429 
    2430         double  K[2][2]={0.0};
    2431         double  KDL[2][2]={0.0};
    2432         int     found=0;
    2433         int     dim;
    2434 
    2435         /*inputs: */
    2436         Input* vxaverage_input=NULL;
    2437         Input* vyaverage_input=NULL;
    2438         bool artdiff;
     2398
     2399        /*Intermediaries */
     2400        bool       artdiff;
     2401        int        i,j,ig,dim;
     2402        double     Jdettria ,vx,vy,dvxdx,dvydy;
     2403        double     dvx[2],dvy[2];
     2404        double     xyz_list[NUMVERTICES][3];
     2405        double     L[NUMVERTICES];
     2406        double     B[2][NUMVERTICES];
     2407        double     Bprime[2][NUMVERTICES];
     2408        double     K[2][2]                          = {0.0};
     2409        double     KDL[2][2]                        = {0.0};
     2410        double     DL[2][2]                         = {0.0};
     2411        double     DLprime[2][2]                    = {0.0};
     2412        double     DL_scalar;
     2413        double     Ke_gg[numdof][numdof]            = {0.0};
     2414        double     Ke_gg_gaussian[numdof][numdof]   = {0.0};
     2415        double     Ke_gg_thickness1[numdof][numdof] = {0.0};
     2416        double     Ke_gg_thickness2[numdof][numdof] = {0.0};
     2417        int       *doflist                          = NULL;
     2418        GaussTria *gauss                            = NULL;
    24392419
    24402420        /* Get node coordinates and dof list: */
    … …  
    24422422        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    24432423
    2444         /*retrieve some parameters: */
     2424        /*Retrieve all Inputs and parameters: */
    24452425        this->parameters->FindParam(&artdiff,ArtDiffEnum);
    24462426        this->parameters->FindParam(&dim,DimEnum);
    2447 
    2448         /*Retrieve all inputs we will be needed*/
     2427        Input* vxaverage_input=NULL;
     2428        Input* vyaverage_input=NULL;
    24492429        if(dim==2){
    24502430                vxaverage_input=inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input);
    … …  
    24582438        //Create Artificial diffusivity once for all if requested
    24592439        if(artdiff){
    2460                 //Get the Jacobian determinant
    24612440                gauss=new GaussTria();
    24622441                gauss->GaussCenter();
    … …  
    24642443                delete gauss;
    24652444
    2466                 //Build K matrix (artificial diffusivity matrix)
    2467                 vxaverage_input->GetParameterAverage(&v_gauss[0]);
    2468                 vyaverage_input->GetParameterAverage(&v_gauss[1]);
    2469 
    2470                 K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]);
    2471                 K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]);
     2445                vxaverage_input->GetParameterAverage(&vx);
     2446                vyaverage_input->GetParameterAverage(&vy);
     2447                K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(vx);
     2448                K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(vy);
    24722449        }
    24732450
    … …  
    24782455                gauss->GaussPoint(ig);
    24792456
    2480                 /* Get Jacobian determinant: */
    24812457                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    2482 
    2483                 /*Get B  and B prime matrix: */
    24842458                GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
    24852459                GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
    24862460
    2487                 //Get vx, vy and their derivatives at gauss point
    24882461                vxaverage_input->GetParameterValue(&vx,gauss);
    24892462                vyaverage_input->GetParameterValue(&vy,gauss);
    2490 
    24912463                vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
    24922464                vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    … …  
    24942466                dvxdx=dvx[0];
    24952467                dvydy=dvy[1];
    2496 
    24972468                DL_scalar=gauss->weight*Jdettria;
    24982469
    2499                 //Create DL and DLprime matrix
    25002470                DL[0][0]=DL_scalar*dvxdx;
    25012471                DL[1][1]=DL_scalar*dvydy;
    … …  
    25032473                DLprime[0][0]=DL_scalar*vx;
    25042474                DLprime[1][1]=DL_scalar*vy;
    2505 
    2506                 //Do the triple product tL*D*L.
    2507                 //Ke_gg_thickness=B'*DLprime*Bprime;
    25082475
    25092476                TripleMultiply( &B[0][0],2,numdof,1,
    … …  
    25172484                                        &Ke_gg_thickness2[0][0],0);
    25182485
    2519                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    25202486                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness1[i][j];
    25212487                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
    25222488
    25232489                if(artdiff){
    2524 
    2525                         /* Compute artificial diffusivity */
    25262490                        KDL[0][0]=DL_scalar*K[0][0];
    25272491                        KDL[1][1]=DL_scalar*K[1][1];
    … …  
    25322496                                                &Ke_gg_gaussian[0][0],0);
    25332497
    2534                         /* Add artificial diffusivity matrix */
    25352498                        for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    2536 
    25372499                }
    25382500        }
    25392501
    2540         /*Add Ke_gg to global matrix Kgg: */
    25412502        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    25422503
    … …  
    25732534        /*input parameters for structural analysis (diagnostic): */
    25742535        double  vx,vy;
    2575         int     found;
    25762536        int     dim;
    25772537
    … …  
    26662626        double  K[2][2]={0.0};
    26672627        double  KDL[2][2]={0.0};
    2668         int     found=0;
    26692628        int     dim;
    26702629
    … …  
    27962755void  Tria::CreateKMatrixDiagnosticMacAyeal(Mat Kgg,Mat Kff, Mat Kfs){
    27972756       
    2798         const int      numdof         = 2 *NUMVERTICES;
     2757        /*Intermediaries*/
    27992758        double        *Ke_gg_viscous  = NULL;
    28002759        double        *Ke_gg_friction = NULL;
    … …  
    28052764        Ke_gg_friction = CreateKMatrixDiagnosticMacAyealFriction();
    28062765
    2807         /*Initialize element matrix: */
     2766        /*Add Ke_gg values to Ke element stifness matrix: */
    28082767        Ke=this->NewElementMatrix(MacAyealApproximationEnum);
    2809 
    2810         /*Add Ke_gg values to Ke element stifness matrix: */
    28112768        if(Ke_gg_viscous) Ke->AddValues(Ke_gg_viscous);
    28122769        if(Ke_gg_friction)Ke->AddValues(Ke_gg_friction);
    … …  
    28192776        xfree((void**)&Ke_gg_viscous);
    28202777        xfree((void**)&Ke_gg_friction);
    2821 
    2822 }
    2823 
     2778}
    28242779/*}}}*/
    28252780/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealViscous{{{1*/
    … …  
    35473502        double  K[2][2]={0.0};
    35483503        double  KDL[2][2]={0.0};
    3549         int     found;
    35503504        int     dim;
    35513505
    … …  
    37083662        /*input parameters for structural analysis (diagnostic): */
    37093663        double  vx,vy;
    3710         int     found;
    37113664        int     dim;
    37123665
    … …  
    38383791
    38393792        int i,j;
    3840         int found=0;
    38413793       
    38423794        /* node data: */
    … …  
    50985050void Tria::CreatePVectorThermalShelf( Vec pg){
    50995051
    5100         int i,found;
     5052        int i;
    51015053       
    51025054        const int  numdof=NUMVERTICES*NDOF1;
    … …  
    51865138void Tria::CreatePVectorThermalSheet( Vec pg){
    51875139
    5188         int i,found;
     5140        int i;
    51895141       
    51905142        const int  numdof=NUMVERTICES*NDOF1;
    … …  
    60165968        int* sinternaldoflist=NULL;
    60175969        int* sexternaldoflist=NULL;
    6018         bool symmetric=true;
     5970        bool square=true;
    60195971
    60205972        /*retrieve some parameters: */
    … …  
    60375989        }
    60385990
    6039         /*Use symmetric constructor for ElementMatrix: */
    6040         if(!kff) Ke=new ElementMatrix(gsize,symmetric,gexternaldoflist);
    6041         else     Ke=new ElementMatrix(gsize,symmetric,finternaldoflist,fexternaldoflist,fsize,sinternaldoflist,sexternaldoflist,ssize);
     5991        /*Use square constructor for ElementMatrix: */
     5992        if(!kff) Ke=new ElementMatrix(gsize,square,gexternaldoflist);
     5993        else     Ke=new ElementMatrix(gsize,square,finternaldoflist,fexternaldoflist,fsize,sinternaldoflist,sexternaldoflist,ssize);
    60425994
    60435995        /*Free ressources and return:*/
    … …  
    60816033        }
    60826034
    6083         /*Use symmetric constructor for ElementVector: */
     6035        /*constructor for ElementVector: */
    60846036        if(!kff)pe=new ElementVector(gsize,gexternaldoflist);
    60856037        else    pe=new ElementVector(gsize,finternaldoflist,fexternaldoflist,fsize);

  • issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp

    r5788 r5821  
    2525        this->nrows=0;
    2626        this->ncols=0;
    27         this->symmetric=false;
     27        this->values=NULL;
     28        this->square=false;
     29        this->kff=false;
     30
    2831        this->row_fsize=0;
    29         this->values=NULL;
    3032        this->row_finternaldoflist=NULL;
    3133        this->row_fexternaldoflist=NULL;
    … …  
    3335        this->row_sinternaldoflist=NULL;
    3436        this->row_sexternaldoflist=NULL;
     37
    3538        this->col_fsize=0;
    3639        this->col_finternaldoflist=NULL;
    … …  
    4245}
    4346/*}}}*/
    44 /*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool symmetric,int* in_gexternaldoflist){{{1*/
    45 ElementMatrix::ElementMatrix(int gsize,bool symmetric,int* in_gexternaldoflist){
    46 
    47         if(symmetric=false)ISSMERROR(" calling symmetric constructor with false symmetric flag!");
     47/*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool square,int* in_gexternaldoflist){{{1*/
     48ElementMatrix::ElementMatrix(int gsize,bool square,int* in_gexternaldoflist){
     49
     50        if(square=false)ISSMERROR(" calling square constructor with false square flag!");
    4851
    4952        this->nrows=gsize;
    5053        this->ncols=gsize;
    51         this->symmetric=true;
     54        this->square=true;
    5255        this->kff=false;
    5356
    … …  
    5962                memcpy(this->gexternaldoflist,in_gexternaldoflist,this->nrows*sizeof(int));
    6063        }
    61         else this->gexternaldoflist=NULL;
     64        else{
     65                this->gexternaldoflist=NULL;
     66        }
    6267
    6368        /*don't do rows and cols, because kff is false:*/
    … …  
    7883}
    7984/*}}}*/
    80 /*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool symmetric,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize){{{1*/
    81 ElementMatrix::ElementMatrix(int gsize,bool symmetric,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize){
    82 
    83         if(symmetric=false)ISSMERROR(" calling symmetric constructor with false symmetric flag!");
     85/*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool square,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize){{{1*/
     86ElementMatrix::ElementMatrix(int gsize,bool square,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize){
     87
     88        if(square=false)ISSMERROR(" calling square constructor with false square flag!");
    8489
    8590        this->nrows=gsize;
    8691        this->ncols=gsize;
    87         this->symmetric=true;
     92        this->square=true;
    8893        this->kff=true;
    8994
    … …  
    149154void ElementMatrix::AddValues(double* Ke_gg){
    150155
    151         int i,j;
    152 
    153156        if(Ke_gg){
    154                 for (i=0;i<this->nrows;i++){
    155                         for(j=0;j<this->ncols;j++){
     157                for (int i=0;i<this->nrows;i++){
     158                        for(int j=0;j<this->ncols;j++){
    156159                                *(this->values+this->ncols*i+j)+=*(Ke_gg+this->ncols*i+j);
    157160                        }
    … …  
    166169        double* internalvalues=NULL;
    167170
    168         if(this->symmetric){
     171        if(this->square){
    169172                /*only use row dofs to add values into global matrices: */
    170173               
    … …  
    207210        }
    208211        else{
    209                 ISSMERROR(" unsymmetric AddToGlobal routine not support yet!");
     212                ISSMERROR(" non square matrix AddToGlobal routine not support yet!");
    210213        }
    211214
    … …  
    219222        printf("   nrows: %i\n",nrows);
    220223        printf("   ncols: %i\n",ncols);
    221         printf("   symmetric: %s\n",symmetric?"true":"false");
     224        printf("   square: %s\n",square?"true":"false");
    222225        printf("   kff: %s\n",kff?"true":"false");
    223226
    … …  
    225228        for(i=0;i<nrows;i++){
    226229                printf("      %i: ",i);
    227                 for(j=0;j<ncols;j++){
    228                         printf("%g ",*(values+ncols*i+j));
    229                 }
     230                for(j=0;j<ncols;j++) printf("%g ",*(values+ncols*i+j));
    230231                printf("\n");
    231232        }
    … …  
    246247        if(row_sexternaldoflist)for(i=0;i<row_ssize;i++)printf("%i ",row_sexternaldoflist[i]); printf("\n");
    247248
    248         if(!symmetric){
    249 
     249        if(!square){
    250250                printf("   col_fsize: %i\n",col_fsize);
    251251                printf("   col_finternaldoflist (%p): ",col_finternaldoflist);

  • issm/trunk/src/c/objects/Numerics/ElementMatrix.h

    r5772 r5821  
    2222                int      ncols;
    2323                double*  values;
    24                 bool     symmetric;
     24                bool     square;
    2525                bool     kff;
    2626
    … …  
    5050                /*ElementMatrix constructors, destructors {{{1*/
    5151                ElementMatrix();
    52                 ElementMatrix(int gsize,bool symmetric,int* gexternaldoflist);
    53                 ElementMatrix(int gsize,bool symmetric,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize);
     52                ElementMatrix(int gsize,bool square,int* gexternaldoflist);
     53                ElementMatrix(int gsize,bool square,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize);
    5454                ~ElementMatrix();
    5555                /*}}}*/
Note: See TracChangeset for help on using the changeset viewer.