Changeset 5840


Ignore:
Timestamp:
09/16/10 10:14:40 (15 years ago)
Author:
Mathieu Morlighem
Message:

Some more ElementMatrix

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

Legend:

Unmodified
Added
Removed

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

    r5839 r5840  
    695695void  Tria::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){
    696696
    697         /*retrive parameters: */
     697        /*retreive parameters: */
     698        ElementMatrix* Ke=NULL;
    698699        int analysis_type;
    699700        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    … …  
    704705
    705706        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    706         if (analysis_type==DiagnosticHorizAnalysisEnum || analysis_type==AdjointHorizAnalysisEnum){
    707                 ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal();
     707        switch(analysis_type){
     708                case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum:
     709                        Ke=CreateKMatrixDiagnosticMacAyeal();
     710                        break;
     711                case DiagnosticHutterAnalysisEnum:
     712                        Ke=CreateKMatrixDiagnosticHutter();
     713                        break;
     714                case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
     715                        Ke=CreateKMatrixSlope();
     716                        break;
     717                case PrognosticAnalysisEnum:
     718                        Ke=CreateKMatrixPrognostic();
     719                        break;
     720                case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum:
     721                        Ke=CreateKMatrixBalancedthickness();
     722                        break;
     723                case BalancedvelocitiesAnalysisEnum:
     724                        Ke=CreateKMatrixBalancedvelocities();
     725                        break;
     726                default:
     727                        ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
     728        }
     729
     730        /*Add to global matrix*/
     731        if(Ke){
    708732                Ke->AddToGlobal(Kgg,Kff,Kfs);
    709733                delete Ke;
    710734        }
    711         else if (analysis_type==DiagnosticHutterAnalysisEnum){
    712                 CreateKMatrixDiagnosticHutter( Kgg);
    713         }
    714         else if (analysis_type==BedSlopeXAnalysisEnum || analysis_type==SurfaceSlopeXAnalysisEnum || analysis_type==BedSlopeYAnalysisEnum || analysis_type==SurfaceSlopeYAnalysisEnum){
    715                 CreateKMatrixSlope( Kgg);
    716         }
    717         else if (analysis_type==PrognosticAnalysisEnum){
    718                 ElementMatrix* Ke=CreateKMatrixPrognostic();
    719                 Ke->AddToGlobal(Kgg,Kff,Kfs);
    720                 delete Ke;
    721         }
    722         else if (analysis_type==BalancedthicknessAnalysisEnum || analysis_type==AdjointBalancedthicknessAnalysisEnum){
    723                 ElementMatrix* Ke=CreateKMatrixBalancedthickness();
    724                 Ke->AddToGlobal(Kgg,Kff,Kfs);
    725                 delete Ke;
    726         }
    727         else if (analysis_type==BalancedvelocitiesAnalysisEnum){
    728                 CreateKMatrixBalancedvelocities( Kgg);
    729         }
    730         else{
    731                 ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
    732         }
    733 
    734735}
    735736/*}}}*/
    … …  
    24942495                                        &Ke_gg_thickness2[0][0],0);
    24952496
    2496                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_thickness1[i][j];
    2497                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_thickness2[i][j];
     2497                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_thickness1[i][j];
     2498                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_thickness2[i][j];
    24982499
    24992500                if(artdiff){
    … …  
    25062507                                                &Ke_gg_gaussian[0][0],0);
    25072508
    2508                         for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j];
     2509                        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j];
    25092510                }
    25102511        }
    … …  
    25782579/*}}}*/
    25792580/*FUNCTION Tria::CreateKMatrixBalancedvelocities {{{1*/
    2580 void  Tria::CreateKMatrixBalancedvelocities(Mat Kgg){
     2581ElementMatrix* Tria::CreateKMatrixBalancedvelocities(void){
    25812582
    25822583        /* local declarations */
     2584        const int    numdof=NDOF1*NUMVERTICES;
    25832585        int             i,j;
    2584 
    2585         /* node data: */
    2586         const int    numdof=NDOF1*NUMVERTICES;
    25872586        double       xyz_list[NUMVERTICES][3];
    2588         int*         doflist=NULL;
    2589 
    2590         /* gaussian points: */
    25912587        int     ig;
    25922588        GaussTria *gauss=NULL;
    2593 
    2594         /* matrices: */
    25952589        double L[NUMVERTICES];
    25962590        double B[2][NUMVERTICES];
    … …  
    26042598        double Ke_gg_velocities2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    26052599        double Jdettria;
    2606 
    2607         /*input parameters for structural analysis (diagnostic): */
    26082600        double  surface_normal[3];
    26092601        double  surface_list[3];
    … …  
    26172609        double  KDL[2][2]={0.0};
    26182610        int     dim;
    2619 
    2620         /*inputs: */
     2611        bool artdiff;
     2612
     2613        /*Initialize Element matrix and return if necessary*/
     2614        if(IsOnWater()) return NULL;
     2615        ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum);
     2616
     2617        /*Retrieve all inputs and parameters*/
     2618        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     2619        this->parameters->FindParam(&artdiff,ArtDiffEnum);
     2620        this->parameters->FindParam(&dim,DimEnum);
     2621        Input* surface_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input);
    26212622        Input* vxaverage_input=NULL;
    26222623        Input* vyaverage_input=NULL;
    2623         bool artdiff;
    2624 
    2625         /*retrieve some parameters: */
    2626         this->parameters->FindParam(&artdiff,ArtDiffEnum);
    2627         this->parameters->FindParam(&dim,DimEnum);
    2628 
    2629         /* Get node coordinates and dof list: */
    2630         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2631         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    2632 
    2633         /*Retrieve all inputs we will be needed*/
    2634         Input* surface_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input);
    26352624        if(dim==2){
    26362625                vxaverage_input=inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input);
    … …  
    26642653        //Create Artificial diffusivity once for all if requested
    26652654        if(artdiff){
    2666                 //Get the Jacobian determinant
    26672655                gauss=new GaussTria();
    26682656                gauss->GaussCenter();
    … …  
    26702658                delete gauss;
    26712659
    2672                 //Build K matrix (artificial diffusivity matrix)
    26732660                vxaverage_input->GetParameterAverage(&v_gauss[0]);
    26742661                vyaverage_input->GetParameterAverage(&v_gauss[1]);
    … …  
    26842671                gauss->GaussPoint(ig);
    26852672
    2686                 /* Get Jacobian determinant: */
    2687                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    2688 
    2689                 /*Get B  and B prime matrix: */
    26902673                GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
    26912674                GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
    2692 
    2693                 //Get vx, vy and their derivatives at gauss point
     2675                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
     2676
    26942677                vxaverage_input->GetParameterValue(&vx,gauss);
    26952678                vyaverage_input->GetParameterValue(&vy,gauss);
    2696 
    26972679                vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
    26982680                vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    2699 
    27002681                dvxdx=dvx[0];
    27012682                dvydy=dvy[1];
    … …  
    27062687                DLprime[1][1]=DL_scalar*ny;
    27072688
    2708                 //Do the triple product tL*D*L.
    2709                 //Ke_gg_velocities=B'*DLprime*Bprime;
    27102689                TripleMultiply( &B[0][0],2,numdof,1,
    27112690                                        &DLprime[0][0],2,2,0,
    … …  
    27132692                                        &Ke_gg_velocities2[0][0],0);
    27142693
    2715                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    2716                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_velocities2[i][j];
     2694                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_velocities2[i][j];
    27172695
    27182696                if(artdiff){
    2719 
    2720                         /* Compute artificial diffusivity */
    27212697                        KDL[0][0]=DL_scalar*K[0][0];
    27222698                        KDL[1][1]=DL_scalar*K[1][1];
    … …  
    27272703                                                &Ke_gg_gaussian[0][0],0);
    27282704
    2729                         /* Add artificial diffusivity matrix */
    2730                         for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     2705                        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j];
    27312706
    27322707                }
    27332708
    27342709        }
    2735 
    2736         /*Add Ke_gg to global matrix Kgg: */
    2737         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    27382710
    27392711        /*Clean up and return*/
    27402712        delete gauss;
    2741         xfree((void**)&doflist);
     2713        return Ke;
    27422714}
    27432715/*}}}*/
    … …  
    31393111/*}}}*/
    31403112/*FUNCTION Tria::CreateKMatrixDiagnosticHutter{{{1*/
    3141 void  Tria::CreateKMatrixDiagnosticHutter(Mat Kgg){
    3142 
    3143         /*Collapsed formulation: */
    3144         int    i;
    3145         int    connectivity;
    3146         const int numdofs=NUMVERTICES*NDOF2;
    3147         int*         doflist=NULL;
    3148         double Ke_gg[numdofs][numdofs]={0.0};
    3149 
    3150         /*If on water, skip: */
    3151         if(IsOnWater())return;
    3152 
    3153         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    3154 
    3155         /*Spawn 3 sing elements: */
    3156         for(i=0;i<3;i++){
     3113ElementMatrix* Tria::CreateKMatrixDiagnosticHutter(void){
     3114
     3115        /*Intermediaries*/
     3116        const int numdof=NUMVERTICES*NDOF2;
     3117        int    i,connectivity;
     3118
     3119        /*Initialize Element matrix and return if necessary*/
     3120        if(IsOnWater()) return NULL;
     3121        ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum);
     3122
     3123        /*Create Element matrix*/
     3124        for(i=0;i<NUMVERTICES;i++){
    31573125                connectivity=nodes[i]->GetConnectivity();
    3158                 Ke_gg[2*i][2*i]=1/(double)connectivity;
    3159                 Ke_gg[2*i+1][2*i+1]=1/(double)connectivity;
    3160         }
    3161 
    3162         MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
    3163 
    3164         /*Free ressources:*/
    3165         xfree((void**)&doflist);
     3126                Ke->values[(2*i)*numdof  +(2*i)  ]=1/(double)connectivity;
     3127                Ke->values[(2*i+1)*numdof+(2*i+1)]=1/(double)connectivity;
     3128        }
     3129
     3130        /*Clean up and return*/
     3131        return Ke;
    31663132}
    31673133/*}}}*/
    … …  
    35703536/*}}}*/
    35713537/*FUNCTION Tria::CreateKMatrixSlope {{{1*/
    3572 void  Tria::CreateKMatrixSlope(Mat Kgg){
     3538ElementMatrix* Tria::CreateKMatrixSlope(void){
    35733539
    35743540        /*constants: */
    3575         const int    numnodes=3;
    3576         const int    numdof=NDOF1*numnodes;
     3541        const int    numdof=NDOF1*NUMVERTICES;
    35773542
    35783543        /* Intermediaries */
    35793544        int        i,j,ig;
    35803545        double     DL_scalar,Jdet;
    3581         double     xyz_list[numnodes][3];
     3546        double     xyz_list[NUMVERTICES][3];
    35823547        double     L[1][3];
    3583         double     Ke[numdof][numdof] = {0.0};
    35843548        double     Ke_g[numdof][numdof];
    35853549        GaussTria *gauss = NULL;
    3586         int       *doflist = NULL;
    3587 
    3588         GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes);
    3589         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     3550
     3551        /*Initialize Element matrix and return if necessary*/
     3552        if(IsOnWater()) return NULL;
     3553        ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum);
     3554
     3555        GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES);
    35903556
    35913557        /* Start looping on the number of gaussian points: */
    … …  
    36053571                                        &Ke_g[0][0],0);
    36063572
    3607                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke[i][j]+=Ke_g[i][j];
    3608         }
    3609 
    3610         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
    3611                
    3612         /*Clean up*/
     3573                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
     3574        }
     3575
     3576        /*Clean up and return*/
    36133577        delete gauss;
    3614         xfree((void**)&doflist);
     3578        return Ke;
    36153579}
    36163580/*}}}*/

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

    r5839 r5840  
    124124                ElementMatrix* CreateKMatrixBalancedthickness_DG(void);
    125125                ElementMatrix* CreateKMatrixBalancedthickness_CG(void);
    126                 void      CreateKMatrixBalancedvelocities(Mat Kgg);
     126                ElementMatrix* CreateKMatrixBalancedvelocities(void);
    127127                ElementMatrix* CreateKMatrixDiagnosticMacAyeal(void);
    128128                ElementMatrix* CreateKMatrixDiagnosticMacAyealViscous(void);
    … …  
    131131                void      CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg);
    132132                void      CreateKMatrixDiagnosticPattynFriction(Mat Kgg);
    133                 void      CreateKMatrixDiagnosticHutter(Mat Kgg);
     133                ElementMatrix* CreateKMatrixDiagnosticHutter(void);
    134134                void      CreateKMatrixDiagnosticSurfaceVert(Mat Kgg);
    135135                void      CreateKMatrixMelting(Mat Kgg);
    … …  
    137137                ElementMatrix* CreateKMatrixPrognostic_CG(void);
    138138                ElementMatrix* CreateKMatrixPrognostic_DG(void);
    139                 void      CreateKMatrixSlope(Mat Kgg);
     139                ElementMatrix* CreateKMatrixSlope(void);
    140140                void      CreateKMatrixThermal(Mat Kgg);
    141141                void      CreatePVectorBalancedthickness_DG(Vec pg);
Note: See TracChangeset for help on using the changeset viewer.