Changeset 5933


Ignore:
Timestamp:
09/21/10 16:56:47 (15 years ago)
Author:
Mathieu Morlighem
Message:

Final ElementVector

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

Legend:

Unmodified
Added
Removed

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

    r5932 r5933  
    752752                case DiagnosticHorizAnalysisEnum:
    753753                        pe=CreatePVectorDiagnosticHoriz();
    754                         if(pe) pe->AddToGlobal(pg,pf);
    755                         if(pe) delete pe;
    756754                        break;
    757755                case AdjointHorizAnalysisEnum:
    758756                        pe=CreatePVectorAdjointHoriz();
    759                         if(pe) pe->AddToGlobal(pg,pf);
    760                         if(pe) delete pe;
    761757                        break;
    762758                case DiagnosticHutterAnalysisEnum:
    763759                        pe=CreatePVectorDiagnosticHutter();
    764                         if(pe) pe->AddToGlobal(pg,pf);
    765                         if(pe) delete pe;
    766760                        break;
    767761                case DiagnosticVertAnalysisEnum:
    768                         CreatePVectorDiagnosticVert( pg);
     762                        pe=CreatePVectorDiagnosticVert();
    769763                        break;
    770764                case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
    771765                        pe=CreatePVectorSlope();
    772                         if(pe) pe->AddToGlobal(pg,pf);
    773                         if(pe) delete pe;
    774766                        break;
    775767                case PrognosticAnalysisEnum:
    776768                        pe=CreatePVectorPrognostic();
    777                         if(pe) pe->AddToGlobal(pg,pf);
    778                         if(pe) delete pe;
    779769                        break;
    780770                case BalancedthicknessAnalysisEnum:
    781771                        pe=CreatePVectorBalancedthickness();
    782                         if(pe) pe->AddToGlobal(pg,pf);
    783                         if(pe) delete pe;
    784772                        break;
    785773                case BalancedvelocitiesAnalysisEnum:
    786774                        pe=CreatePVectorBalancedvelocities();
    787                         if(pe) pe->AddToGlobal(pg,pf);
    788                         if(pe) delete pe;
    789775                        break;
    790776                case ThermalAnalysisEnum:
    791777                        pe=CreatePVectorThermal();
    792                         if(pe) pe->AddToGlobal(pg,pf);
    793                         if(pe) delete pe;
    794778                        break;
    795779                case MeltingAnalysisEnum:
    796780                        pe=CreatePVectorMelting();
    797                         if(pe) pe->AddToGlobal(pg,pf);
    798                         if(pe) delete pe;
    799781                        break;
    800782                default:
    … …  
    802784        }
    803785
     786        /*Add to global Vector*/
     787        if(pe){
     788                pe->AddToGlobal(pg,pf);
     789                delete pe;
     790        }
    804791}
    805792/*}}}*/
    … …  
    30583045ElementVector* Penta::CreatePVectorCouplingPattynStokes(void){
    30593046
    3060         /*compute all load vectorsfor this element*/
     3047        /*compute all load vectors for this element*/
    30613048        ElementVector* pe1=CreatePVectorCouplingPattynStokesViscous();
    30623049        ElementVector* pe2=CreatePVectorCouplingPattynStokesFriction();
    … …  
    32533240ElementVector* Penta::CreatePVectorDiagnosticMacAyealPattyn(void){
    32543241
    3255         /*compute all load vectorsfor this element*/
     3242        /*compute all load vectors for this element*/
    32563243        ElementVector* pe1=CreatePVectorDiagnosticMacAyeal();
    32573244        ElementVector* pe2=CreatePVectorDiagnosticPattyn();
    … …  
    32673254ElementVector* Penta::CreatePVectorDiagnosticPattynStokes(void){
    32683255
    3269         /*compute all load vectorsfor this element*/
     3256        /*compute all load vectors for this element*/
    32703257        ElementVector* pe1=CreatePVectorDiagnosticPattyn();
    32713258        ElementVector* pe2=CreatePVectorDiagnosticStokes();
    … …  
    34353422ElementVector* Penta::CreatePVectorDiagnosticStokes(void){
    34363423
    3437         /*compute all load vectorsfor this element*/
     3424        /*compute all load vectors for this element*/
    34383425        ElementVector* pe1=CreatePVectorDiagnosticStokesViscous();
    34393426        ElementVector* pe2=CreatePVectorDiagnosticStokesShelf();
    … …  
    36423629/*}}}*/
    36433630/*FUNCTION Penta::CreatePVectorDiagnosticVert {{{1*/
    3644 void  Penta::CreatePVectorDiagnosticVert( Vec pg){
    3645 
    3646         int i;
    3647 
    3648         /* node data: */
     3631ElementVector* Penta::CreatePVectorDiagnosticVert(void){
     3632
     3633        /*compute all load vectors for this element*/
     3634        ElementVector* pe1=CreatePVectorDiagnosticVertVolume();
     3635        ElementVector* pe2=CreatePVectorDiagnosticVertBase();
     3636        ElementVector* pe =new ElementVector(pe1,pe2);
     3637
     3638        /*clean-up and return*/
     3639        delete pe1;
     3640        delete pe2;
     3641        return pe;
     3642}
     3643/*}}}*/
     3644/*FUNCTION Penta::CreatePVectorDiagnosticVertVolume {{{1*/
     3645ElementVector* Penta::CreatePVectorDiagnosticVertVolume(void){
     3646
    36493647        const int    numdof=NDOF1*NUMVERTICES;
    36503648        double       xyz_list[NUMVERTICES][3];
    36513649        int*         doflist=NULL;
    3652 
    3653         /* gaussian points: */
     3650        int i;
    36543651        int     ig;
    36553652        GaussPenta *gauss=NULL;
    3656 
    3657         /* Jacobian: */
    36583653        double Jdet;
    3659 
    3660         /*element vector at the gaussian points: */
    3661         double  pe_g[numdof]={0.0};
    3662         double  pe_g_gaussian[numdof];
    36633654        double l1l6[6];
    3664 
    3665         /*Spawning: */
    36663655        Tria* tria=NULL;
    3667 
    3668         /*input parameters for structural analysis (diagnostic): */
    36693656        double du[3];
    36703657        double dv[3];
    36713658        double dw[3];
    36723659        double dudx,dvdy,dwdz;
    3673 
    3674         /*Get some parameters*/
    36753660        int approximation;
    36763661
    3677         /*If on water, skip: */
    3678         if(IsOnWater())return;
     3662        /*Initialize Element vector and return if necessary*/
     3663        if(IsOnWater()) return NULL;
     3664        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters);
     3665
     3666        /*Retrieve all inputs and parameters*/
     3667        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    36793668        inputs->GetParameterValue(&approximation,ApproximationEnum);
    3680 
    3681         /*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the
    3682          *diagnostic base vertical stifness: */
    3683         if(IsOnBed()){
    3684                 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock
    3685                 tria->CreatePVectorDiagnosticBaseVert(pg);
    3686                 delete tria->matice; delete tria;
    3687         }
    3688 
    3689         /*Now, handle the standard penta element: */
    3690         /* Get node coordinates and dof list: */
    3691         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3692         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    3693 
    3694         /*Retrieve all inputs we will be needing: */
    36953669        Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
    36963670        Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
    … …  
    37063680                gauss->GaussPoint(ig);
    37073681
    3708                 /*Get velocity derivative, with respect to x and y: */
    37093682                vx_input->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss);
    37103683                vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss);
    … …  
    37173690                dvdy=dv[1];
    37183691
    3719                 /* Get Jacobian determinant: */
    37203692                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3721 
    3722                 /*Get nodal functions: */
    37233693                GetNodalFunctionsP1(l1l6, gauss);
    37243694
    3725                 /*Build pe_g_gaussian vector: */
    3726                 for (i=0;i<NUMVERTICES;i++){
    3727                         pe_g_gaussian[i]=(dudx+dvdy+dwdz)*Jdet*gauss->weight*l1l6[i];
    3728                 }
    3729 
    3730                 /*Add pe_g_gaussian vector to pe_g: */
    3731                 for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
    3732         }
    3733 
    3734         /*Add pe_g to global vector pg: */
    3735         VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    3736 
    3737         xfree((void**)&doflist);
     3695                for (i=0;i<numdof;i++) pe->values[i] += (dudx+dvdy+dwdz)*Jdet*gauss->weight*l1l6[i];
     3696        }
     3697
     3698        /*Clean up and return*/
    37383699        delete gauss;
     3700        return pe;
     3701}
     3702/*}}}*/
     3703/*FUNCTION Penta::CreatePVectorDiagnosticVertBase {{{1*/
     3704ElementVector* Penta::CreatePVectorDiagnosticVertBase(void){
     3705
     3706        if (!IsOnBed() || IsOnWater()) return NULL;
     3707
     3708        /*Call Tria function*/
     3709        Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     3710        ElementVector* pe=tria->CreatePVectorDiagnosticBaseVert();
     3711        delete tria->matice; delete tria;
     3712
     3713        /*Clean up and return*/
     3714        return pe;
    37393715}
    37403716/*}}}*/
    … …  
    37833759ElementVector* Penta::CreatePVectorThermal(void){
    37843760
    3785         /*compute all load vectorsfor this element*/
     3761        /*compute all load vectors for this element*/
    37863762        ElementVector* pe1=CreatePVectorThermalVolume();
    37873763        ElementVector* pe2=CreatePVectorThermalSheet();

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

    r5932 r5933  
    168168                ElementVector* CreatePVectorDiagnosticStokesShelf(void);
    169169                ElementVector* CreatePVectorAdjointStokes(void);
    170                 void      CreatePVectorDiagnosticVert( Vec pg);
     170                ElementVector* CreatePVectorDiagnosticVert(void);
     171                ElementVector* CreatePVectorDiagnosticVertVolume(void);
     172                ElementVector* CreatePVectorDiagnosticVertBase(void);
    171173                ElementVector* CreatePVectorMelting(void);
    172174                ElementVector* CreatePVectorPrognostic(void);

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

    r5932 r5933  
    35953595/*}}}*/
    35963596/*FUNCTION Tria::CreatePVectorDiagnosticBaseVert {{{1*/
    3597 void  Tria::CreatePVectorDiagnosticBaseVert(Vec pg){
     3597ElementVector* Tria::CreatePVectorDiagnosticBaseVert(void){
    35983598
    35993599        /*Constants*/
    … …  
    36033603        int        i,j,ig;
    36043604        int        approximation;
    3605         int*       doflist=NULL;
    36063605        double     xyz_list[NUMVERTICES][3];
    36073606        double     Jdet;
    … …  
    36093608        double     slope[2];
    36103609        double     L[NUMVERTICES];
    3611         double     pe_g[numdof]={0.0};
    3612         double     pe_g_gaussian[numdof];
    36133610        GaussTria* gauss=NULL;
    36143611
    3615         /* Get node coordinates and dof list: */
     3612        /*Initialize Element vector and return if necessary*/
     3613        if(IsOnWater()) return NULL;
     3614        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters);
     3615
     3616        /*Retrieve all inputs and parameters*/
    36163617        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3617         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    3618 
    3619         /*Retrieve all inputs we will be needing: */
    36203618        inputs->GetParameterValue(&approximation,ApproximationEnum);
    36213619        Input* bed_input=inputs->GetInput(BedEnum);             ISSMASSERT(bed_input);
    … …  
    36493647                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    36503648
    3651                 for(i=0;i<NUMVERTICES;i++){
    3652                         pe_g_gaussian[i]=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-meltingvalue)*L[i];
    3653                 }
    3654 
    3655                 for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
    3656         }
    3657 
    3658         VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
     3649                for(i=0;i<numdof;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-meltingvalue)*L[i];
     3650        }
    36593651
    36603652        /*Clean up and return*/
    36613653        delete gauss;
    3662         xfree((void**)&doflist);
     3654        return pe;
    36633655}
    36643656/*}}}*/

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

    r5932 r5933  
    142142                ElementVector* CreatePVectorBalancedthickness_CG(void);
    143143                ElementVector* CreatePVectorBalancedvelocities(void);
    144                 void      CreatePVectorDiagnosticBaseVert(Vec pg);
     144                ElementVector* CreatePVectorDiagnosticBaseVert(void);
    145145                ElementVector* CreatePVectorDiagnosticMacAyeal(void);
    146146                ElementVector* CreatePVectorAdjointHoriz(void);
Note: See TracChangeset for help on using the changeset viewer.