Changeset 5932


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

more ElementVector, almost done

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

Legend:

Unmodified
Added
Removed

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

    r5931 r5932  
    789789                        break;
    790790                case ThermalAnalysisEnum:
    791                         CreatePVectorThermal( pg);
     791                        pe=CreatePVectorThermal();
     792                        if(pe) pe->AddToGlobal(pg,pf);
     793                        if(pe) delete pe;
    792794                        break;
    793795                case MeltingAnalysisEnum:
    794                         CreatePVectorMelting( pg);
     796                        pe=CreatePVectorMelting();
     797                        if(pe) pe->AddToGlobal(pg,pf);
     798                        if(pe) delete pe;
    795799                        break;
    796800                default:
    … …  
    37363740/*}}}*/
    37373741/*FUNCTION Penta::CreatePVectorMelting {{{1*/
    3738 void Penta::CreatePVectorMelting( Vec pg){
    3739         return;
     3742ElementVector* Penta::CreatePVectorMelting(void){
     3743        return NULL;
    37403744}
    37413745/*}}}*/
    … …  
    37773781/*}}}*/
    37783782/*FUNCTION Penta::CreatePVectorThermal {{{1*/
    3779 void Penta::CreatePVectorThermal( Vec pg){
    3780 
    3781         /*indexing: */
     3783ElementVector* Penta::CreatePVectorThermal(void){
     3784
     3785        /*compute all load vectorsfor this element*/
     3786        ElementVector* pe1=CreatePVectorThermalVolume();
     3787        ElementVector* pe2=CreatePVectorThermalSheet();
     3788        ElementVector* pe3=CreatePVectorThermalShelf();
     3789        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
     3790
     3791        /*clean-up and return*/
     3792        delete pe1;
     3793        delete pe2;
     3794        delete pe3;
     3795        return pe;
     3796}
     3797/*}}}*/
     3798/*FUNCTION Penta::CreatePVectorThermalVolume {{{1*/
     3799ElementVector* Penta::CreatePVectorThermalVolume(void){
     3800
     3801        const int  numdof=NUMVERTICES*NDOF1;
    37823802        int i,j;
     3803        int     ig;
    37833804        int found=0;
    3784 
    3785         const int  numdof=NUMVERTICES*NDOF1;
    3786         int*       doflist=NULL;
    3787 
    3788         /*Grid data: */
    37893805        double        xyz_list[NUMVERTICES][3];
    3790 
    3791         /* gaussian points: */
    3792         int     ig;
    37933806        GaussPenta *gauss=NULL;
    3794 
    37953807        double temperature_list[NUMVERTICES];
    37963808        double temperature;
    3797 
    3798         /*Material properties: */
    37993809        double gravity,rho_ice,rho_water;
    38003810        double mixed_layer_capacity,heatcapacity;
    38013811        double beta,meltingpoint,thermal_exchange_velocity;
    3802 
    3803         /* element parameters: */
    38043812        int    friction_type;
    3805 
    3806         /*matrices: */
    38073813        double P_terms[numdof]={0.0};
    38083814        double L[numdof];
    … …  
    38123818        double epsilon_sqr[3][3];
    38133819        double epsilon_matrix[3][3];
    3814 
    38153820        double Jdet;
    38163821        double viscosity;
    … …  
    38223827        double scalar_ocean;
    38233828        double scalar_transient;
    3824 
    3825         /*Collapsed formulation: */
    38263829        Tria*  tria=NULL;
    38273830        double dt;
    38283831
    3829         /*retrieve some parameters: */
     3832        /*Initialize Element vector and return if necessary*/
     3833        if(IsOnWater()) return NULL;
     3834        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters);
     3835
     3836        /*Retrieve all inputs and parameters*/
     3837        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    38303838        this->parameters->FindParam(&dt,DtEnum);
    3831 
    3832         /*If on water, skip: */
    3833         if(IsOnWater())return;
    3834 
    3835         /* Get node coordinates and dof list: */
    3836         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3837         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    3838 
    3839         /*recovre material parameters: */
    38403839        rho_water=matpar->GetRhoWater();
    38413840        rho_ice=matpar->GetRhoIce();
    … …  
    38443843        beta=matpar->GetBeta();
    38453844        meltingpoint=matpar->GetMeltingPoint();
    3846 
    3847         /*Retrieve all inputs we will be needing: */
    38483845        Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
    38493846        Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
    … …  
    38603857                gauss->GaussPoint(ig);
    38613858
    3862                 /*Compute strain rate and viscosity: */
    38633859                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    38643860                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    3865 
    3866                 /* Get Jacobian determinant: */
    38673861                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3868 
    3869                 /* Get nodal functions */
    38703862                GetNodalFunctionsP1(&L[0], gauss);
    3871 
    3872                 /*Build deformational heating: */
    38733863                GetPhi(&phi, &epsilon[0], viscosity);
    38743864
    3875                 /*Build pe_gaussian */
    38763865                scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight;
    38773866                if(dt) scalar_def=scalar_def*dt;
    38783867
    3879                 for(i=0;i<NUMVERTICES;i++) P_terms[i]+=scalar_def*L[i];
     3868                for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_def*L[i];
    38803869
    38813870                /* Build transient now */
    … …  
    38833872                        temperature_input->GetParameterValue(&temperature, gauss);
    38843873                        scalar_transient=temperature*Jdet*gauss->weight;
    3885                         for(i=0;i<NUMVERTICES;i++) P_terms[i]+=scalar_transient*L[i];
    3886                 }
    3887         }
    3888 
    3889         /*Add pe_g to global vector pg: */
    3890         VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
     3874                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_transient*L[i];
     3875                }
     3876        }
     3877
     3878        /*Clean up and return*/
     3879        delete gauss;
     3880        return pe;
     3881}
     3882/*}}}*/
     3883/*FUNCTION Penta::CreatePVectorThermalShelf {{{1*/
     3884ElementVector* Penta::CreatePVectorThermalShelf(void){
    38913885
    38923886        /* Ice/ocean heat exchange flux on ice shelf base */
    3893         if(IsOnBed() && IsOnShelf()){
    3894 
    3895                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    3896                 tria->CreatePVectorThermalShelf(pg);
    3897                 delete tria->matice; delete tria;
    3898         }
     3887        if (!IsOnBed() || !IsOnShelf() || IsOnWater()) return NULL;
     3888
     3889        /*Call Tria function*/
     3890        Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     3891        ElementVector* pe=tria->CreatePVectorThermalShelf();
     3892        delete tria->matice; delete tria;
     3893
     3894        /*Clean up and return*/
     3895        return pe;
     3896}
     3897/*}}}*/
     3898/*FUNCTION Penta::CreatePVectorThermalSheet {{{1*/
     3899ElementVector* Penta::CreatePVectorThermalSheet(void){
    38993900
    39003901        /* Geothermal flux on ice sheet base and basal friction */
    3901         if(IsOnBed() && !IsOnShelf()){
    3902 
    3903                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    3904                 tria->CreatePVectorThermalSheet(pg);
    3905                 delete tria->matice; delete tria;
    3906         }
    3907 
    3908         xfree((void**)&doflist);
    3909         delete gauss;
    3910 
     3902        if (!IsOnBed() || IsOnShelf() || IsOnWater()) return NULL;
     3903
     3904        /*Call Tria function*/
     3905        Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     3906        ElementVector* pe=tria->CreatePVectorThermalSheet();
     3907        delete tria->matice; delete tria;
     3908
     3909        /*Clean up and return*/
     3910        return pe;
    39113911}
    39123912/*}}}*/

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

    r5931 r5932  
    169169                ElementVector* CreatePVectorAdjointStokes(void);
    170170                void      CreatePVectorDiagnosticVert( Vec pg);
    171                 void      CreatePVectorMelting( Vec pg);
     171                ElementVector* CreatePVectorMelting(void);
    172172                ElementVector* CreatePVectorPrognostic(void);
    173173                ElementVector* CreatePVectorSlope(void);
    174                 void      CreatePVectorThermal( Vec pg);
     174                ElementVector* CreatePVectorThermal(void);
     175                ElementVector* CreatePVectorThermalVolume(void);
     176                ElementVector* CreatePVectorThermalShelf(void);
     177                ElementVector* CreatePVectorThermalSheet(void);
    175178                void      GetDofList(int** pdoflist,int approximation_enum,int setenum);
    176179                void      GetDofList1(int* doflist);

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

    r5926 r5932  
    43884388/*}}}*/
    43894389/*FUNCTION Tria::CreatePVectorThermalShelf {{{1*/
    4390 void Tria::CreatePVectorThermalShelf( Vec pg){
    4391 
    4392         int i;
     4390ElementVector* Tria::CreatePVectorThermalShelf(void){
    43934391       
    43944392        const int  numdof=NUMVERTICES*NDOF1;
    4395         int*         doflist=NULL;
     4393        int i;
    43964394        double       xyz_list[NUMVERTICES][3];
    4397 
    43984395        double mixed_layer_capacity;
    43994396        double thermal_exchange_velocity;
    … …  
    44034400        double beta;
    44044401        double meltingpoint;
    4405 
    4406         /*inputs: */
    44074402        double dt;
    44084403        double pressure;
    4409 
    4410         /* gaussian points: */
    44114404        int     ig;
    44124405        GaussTria* gauss=NULL;
    4413 
    4414         /*matrices: */
    44154406        double  Jdet;
    4416         double  P_terms[numdof]={0.0};
    44174407        double  l1l2l3[NUMVERTICES];
    44184408        double  t_pmp;
    44194409        double  scalar_ocean;
    44204410
    4421         /* Get node coordinates and dof list: */
     4411        /*Initialize Element vector and return if necessary*/
     4412        if(IsOnWater()) return NULL;
     4413        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters);
     4414
     4415        /*Retrieve all inputs and parameters*/
    44224416        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4423         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    4424 
    4425         //recover material parameters
    44264417        mixed_layer_capacity=matpar->GetMixedLayerCapacity();
    44274418        thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
    … …  
    44314422        beta=matpar->GetBeta();
    44324423        meltingpoint=matpar->GetMeltingPoint();
    4433        
    4434         /*retrieve some solution parameters: */
    44354424        this->parameters->FindParam(&dt,DtEnum);
    4436 
    4437         /* Ice/ocean heat exchange flux on ice shelf base */
    4438 
    4439         /*Retrieve all inputs we will be needing: */
    44404425        Input* pressure_input=inputs->GetInput(PressureEnum); ISSMASSERT(pressure_input);
    44414426
    … …  
    44464431                gauss->GaussPoint(ig);
    44474432
    4448                 //Get the Jacobian determinant
    44494433                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss);
    4450 
    4451                 /*Get nodal functions values: */
    44524434                GetNodalFunctions(&l1l2l3[0], gauss);
    44534435
    4454                 /*Get geothermal flux and basal friction */
    44554436                pressure_input->GetParameterValue(&pressure,gauss);
    44564437                t_pmp=meltingpoint-beta*pressure;
    44574438
    4458                 /*Calculate scalar parameter*/
    44594439                scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice);
    44604440                if(dt){
    … …  
    44624442                }
    44634443
    4464                 for(i=0;i<3;i++){
    4465                         P_terms[i]+=scalar_ocean*l1l2l3[i];
    4466                 }
    4467         }
    4468 
    4469         /*Add pe_g to global vector pg: */
    4470         VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
     4444                for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*l1l2l3[i];
     4445        }
    44714446
    44724447        /*Clean up and return*/
    44734448        delete gauss;
    4474         xfree((void**)&doflist);
     4449        return pe;
    44754450}
    44764451/*}}}*/
    44774452/*FUNCTION Tria::CreatePVectorThermalSheet {{{1*/
    4478 void Tria::CreatePVectorThermalSheet( Vec pg){
    4479 
     4453ElementVector* Tria::CreatePVectorThermalSheet(void){
     4454
     4455        const int  numdof=NUMVERTICES*NDOF1;
    44804456        int i;
    4481        
    4482         const int  numdof=NUMVERTICES*NDOF1;
    4483         int*       doflist=NULL;
     4457        int     ig;
    44844458        double     xyz_list[NUMVERTICES][3];
    4485 
    44864459        double rho_ice;
    44874460        double heatcapacity;
    4488 
    4489         /*inputs: */
    44904461        double dt;
    44914462        double pressure_list[3];
    … …  
    44964467        double alpha2,vx,vy;
    44974468        double geothermalflux_value;
    4498 
    4499         /* gaussian points: */
    4500         int     ig;
    45014469        GaussTria* gauss=NULL;
    4502 
    4503         /*matrices: */
    45044470        double  Jdet;
    45054471        double  P_terms[numdof]={0.0};
    45064472        double  l1l2l3[NUMVERTICES];
    45074473        double  scalar;
    4508 
    45094474        int analysis_type;
    45104475
    4511         /*retrive parameters: */
     4476        /*Initialize Element vector and return if necessary*/
     4477        if(IsOnWater()) return NULL;
     4478        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters);
     4479
     4480        /*Retrieve all inputs and parameters*/
     4481        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    45124482        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    4513 
    4514         /* Get node coordinates and dof list: */
    4515         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4516         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    4517 
    4518         //recover material parameters
    45194483        rho_ice=matpar->GetRhoIce();
    45204484        heatcapacity=matpar->GetHeatCapacity();
    4521 
    4522         /*retrieve some parameters: */
    45234485        this->parameters->FindParam(&dt,DtEnum);
     4486        Input* vx_input=inputs->GetInput(VxEnum);                         ISSMASSERT(vx_input);
     4487        Input* vy_input=inputs->GetInput(VyEnum);                         ISSMASSERT(vy_input);
     4488        Input* vz_input=inputs->GetInput(VzEnum);                         ISSMASSERT(vz_input);
     4489        Input* geothermalflux_input=inputs->GetInput(GeothermalFluxEnum); ISSMASSERT(geothermalflux_input);
    45244490
    45254491        /*Build frictoin element, needed later: */
    … …  
    45284494        friction=new Friction("3d",inputs,matpar,analysis_type);
    45294495
    4530         /*Retrieve all inputs we will be needing: */
    4531         Input* vx_input=inputs->GetInput(VxEnum);                         ISSMASSERT(vx_input);
    4532         Input* vy_input=inputs->GetInput(VyEnum);                         ISSMASSERT(vy_input);
    4533         Input* vz_input=inputs->GetInput(VzEnum);                         ISSMASSERT(vz_input);
    4534         Input* geothermalflux_input=inputs->GetInput(GeothermalFluxEnum); ISSMASSERT(geothermalflux_input);
    4535 
    4536         /* Ice/ocean heat exchange flux on ice shelf base */
    4537 
    45384496        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    45394497        gauss=new GaussTria(2);
    … …  
    45424500                gauss->GaussPoint(ig);
    45434501
    4544                 //Get the Jacobian determinant
    45454502                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss);
    4546 
    4547                 /*Get nodal functions values: */
    45484503                GetNodalFunctions(&l1l2l3[0], gauss);
    45494504
    4550                 /*Get geothermal flux and basal friction */
    45514505                geothermalflux_input->GetParameterValue(&geothermalflux_value,gauss);
    4552        
    4553                 /*Friction: */
    45544506                friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
    45554507               
    4556                 /*Calculate scalar parameter*/
    45574508                scalar=gauss->weight*Jdet*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
    45584509                if(dt){
    … …  
    45604511                }
    45614512
    4562                 for(i=0;i<3;i++){
    4563                         P_terms[i]+=scalar*l1l2l3[i];
    4564                 }
    4565 
    4566         }
    4567 
    4568         /*Add pe_g to global vector pg: */
    4569         VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
     4513                for(i=0;i<numdof;i++) pe->values[i]+=scalar*l1l2l3[i];
     4514        }
    45704515
    45714516        /*Clean up and return*/
    45724517        delete gauss;
    45734518        delete friction;
    4574         xfree((void**)&doflist);
     4519        return pe;
    45754520}
    45764521/*}}}*/

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

    r5926 r5932  
    152152                ElementVector* CreatePVectorPrognostic_DG(void);
    153153                ElementVector* CreatePVectorSlope(void);
    154                 void      CreatePVectorThermalSheet( Vec pg);
    155                 void      CreatePVectorThermalShelf( Vec pg);
     154                ElementVector* CreatePVectorThermalSheet(void);
     155                ElementVector* CreatePVectorThermalShelf(void);
    156156                double  GetArea(void);
    157157                int     GetElementType(void);
Note: See TracChangeset for help on using the changeset viewer.