Changeset 5936


Ignore:
Timestamp:
09/22/10 08:28:23 (15 years ago)
Author:
Mathieu Morlighem
Message:

more element Matrix and Vectors

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

Legend:

Unmodified
Added
Removed

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

    r5934 r5936  
    688688
    689689        /*retrieve parameters: */
     690        ElementMatrix* Ke=NULL;
    690691        int analysis_type;
    691         ElementMatrix* Ke=NULL;
    692692        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    693693

  • issm/trunk/src/c/objects/Loads/Pengrid.cpp

    r5772 r5936  
    33 */
    44
    5 
     5/*Headers*/
     6/*{{{1*/
    67#ifdef HAVE_CONFIG_H
    78        #include "config.h"
    … …  
    1718#include "../../shared/shared.h"
    1819#include "../../Container/Container.h"
    19        
     20/*}}}*/
     21       
     22/*Element macros*/
     23#define NUMVERTICES   1
     24#define NDOF1 1
     25#define NDOF2 2
     26#define NDOF3 3
     27#define NDOF4 4
     28
    2029/*Pengrid constructors and destructor*/
    2130/*FUNCTION Pengrid::Pengrid(){{{1*/
    … …  
    299308void  Pengrid::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,double kmax){
    300309
     310        /*Retrieve parameters: */
     311        ElementMatrix* Ke=NULL;
    301312        int analysis_type;
    302 
    303         /*Retrieve parameters: */
    304313        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    305314
    306         if (analysis_type==DiagnosticHorizAnalysisEnum || analysis_type==AdjointHorizAnalysisEnum){
    307                 PenaltyCreateKMatrixDiagnosticStokes( Kgg,kmax);
    308         }
    309         else if (analysis_type==ThermalAnalysisEnum){
    310                 PenaltyCreateKMatrixThermal( Kgg,kmax);
    311         }
    312         else if (analysis_type==MeltingAnalysisEnum){
    313                 PenaltyCreateKMatrixMelting( Kgg,kmax);
    314         }
    315         else{
    316                 ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
    317         }
    318 
     315        switch(analysis_type){
     316                case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum:
     317                        Ke=PenaltyCreateKMatrixDiagnosticStokes(kmax);
     318                        break;
     319                case ThermalAnalysisEnum:
     320                        Ke=PenaltyCreateKMatrixThermal(kmax);
     321                        break;
     322                case MeltingAnalysisEnum:
     323                        Ke=PenaltyCreateKMatrixMelting(kmax);
     324                        break;
     325                default:
     326                        ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
     327        }
     328
     329        /*Add to global matrix*/
     330        if(Ke){
     331                Ke->AddToGlobal(Kgg,Kff,Kfs);
     332                delete Ke;
     333        }
    319334}
    320335/*}}}1*/
    … …  
    322337void  Pengrid::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){
    323338
     339        /*Retrieve parameters: */
     340        ElementVector* pe=NULL;
    324341        int analysis_type;
    325 
    326         /*Retrieve parameters: */
    327342        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    328343
    329         if (analysis_type==ThermalAnalysisEnum){
    330                 PenaltyCreatePVectorThermal( pg,kmax);
    331         }
    332         else if (analysis_type==MeltingAnalysisEnum){
    333                 PenaltyCreatePVectorMelting( pg,kmax);
    334         }
    335         else if (analysis_type==DiagnosticHorizAnalysisEnum || analysis_type==AdjointHorizAnalysisEnum){
    336 
    337                 /*No loads applied, do nothing: */
    338                 return;
    339 
    340         }
    341         else{
    342                 ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
    343         }
    344 
     344        switch(analysis_type){
     345                case ThermalAnalysisEnum:
     346                        pe=PenaltyCreatePVectorThermal(kmax);
     347                        break;
     348                case MeltingAnalysisEnum:
     349                        pe=PenaltyCreatePVectorMelting(kmax);
     350                        break;
     351                case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum:
     352                        break;
     353                default:
     354                        ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
     355        }
     356
     357        /*Add to global Vector*/
     358        if(pe){
     359                pe->AddToGlobal(pg,pf);
     360                delete pe;
     361        }
    345362}
    346363/*}}}1*/
    … …  
    541558/*}}}1*/
    542559/*FUNCTION Pengrid::PenaltyCreateKMatrixDiagnosticStokes {{{1*/
    543 void  Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax){
    544        
    545         const int numgrids=1;
    546         const int NDOF4=4;
    547         const int numdof=numgrids*NDOF4;
    548         int*         doflist=NULL;
    549 
    550         int dofs1[1]={0};
    551         int dofs2[1]={1};
    552         double slope[2];
    553         int found=0;
    554         double Ke[4][4]={0.0};
    555         double penalty_offset;
    556         int approximation;
    557 
    558         /*pointers: */
    559         Penta* penta=NULL;
    560 
    561         /*recover pointers: */
    562         penta=(Penta*)element;
    563 
    564         /*If not Stokes, return*/
     560ElementMatrix* Pengrid::PenaltyCreateKMatrixDiagnosticStokes(double kmax){
     561       
     562        const int numdof = NUMVERTICES *NDOF4;
     563        double    slope[2];
     564        double    penalty_offset;
     565        int       approximation;
     566
     567        Penta* penta=(Penta*)element;
     568
     569        /*Initialize Element vector and return if necessary*/
    565570        penta->inputs->GetParameterValue(&approximation,ApproximationEnum);
    566         if(approximation!=StokesApproximationEnum) return;
    567 
    568         /*Get dof list: */
    569         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    570 
    571         //recover slope: */
     571        if(approximation!=StokesApproximationEnum) return NULL;
     572        ElementMatrix* Ke=NewElementMatrix(&node,1,this->parameters,StokesApproximationEnum);
     573
     574        /*Retrieve all inputs and parameters*/
     575        parameters->FindParam(&penalty_offset,PenaltyOffsetEnum);
    572576        penta->GetParameterValue(&slope[0],node,BedSlopeXEnum);
    573577        penta->GetParameterValue(&slope[1],node,BedSlopeYEnum);
    574        
     578
     579        /*Create elementary matrix: add penalty to constrain wb (wb=ub*db/dx+vb*db/dy)*/
     580        Ke->values[2*NDOF4+0]=-slope[0]*kmax*pow((double)10.0,penalty_offset);
     581        Ke->values[2*NDOF4+1]=-slope[1]*kmax*pow((double)10.0,penalty_offset);
     582        Ke->values[2*NDOF4+2]= kmax*pow((double)10,penalty_offset);
     583
     584        /*Clean up and return*/
     585        return Ke;
     586}
     587/*}}}1*/
     588/*FUNCTION Pengrid::PenaltyCreateKMatrixMelting {{{1*/
     589ElementMatrix* Pengrid::PenaltyCreateKMatrixMelting(double kmax){
     590
     591        const int numdof=NUMVERTICES*NDOF1;
     592        double pressure,temperature,t_pmp;
     593        double penalty_offset;
     594
     595        Penta* penta=(Penta*)element;
     596
     597        /*check that pengrid is not a clone (penalty to be added only once)*/
     598        if (node->IsClone()) return NULL;
     599        ElementMatrix* Ke=NewElementMatrix(&node,1,this->parameters);
     600
     601        /*Retrieve all inputs and parameters*/
     602        penta->GetParameterValue(&pressure,node,PressureEnum);
     603        penta->GetParameterValue(&temperature,node,TemperatureEnum);
     604        parameters->FindParam(&penalty_offset,PenaltyOffsetEnum);
     605       
     606        /*Compute pressure melting point*/
     607        t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
     608
     609        /*Add penalty load*/
     610        if (temperature<t_pmp){ //If T<Tpmp, there must be no melting. Therefore, melting should be  constrained to 0 when T<Tpmp, instead of using spcs, use penalties
     611                Ke->values[0]=kmax*pow((double)10,penalty_offset);
     612        }
     613
     614        /*Clean up and return*/
     615        return Ke;
     616}
     617/*}}}1*/
     618/*FUNCTION Pengrid::PenaltyCreateKMatrixThermal {{{1*/
     619ElementMatrix* Pengrid::PenaltyCreateKMatrixThermal(double kmax){
     620
     621        const int numdof=NUMVERTICES*NDOF1;
     622        double    penalty_offset;
     623
     624        /*Initialize Element matrix and return if necessary*/
     625        if(!this->active) return NULL;
     626        ElementMatrix* Ke=NewElementMatrix(&node,NUMVERTICES,this->parameters);
     627
    575628        /*recover parameters: */
    576629        parameters->FindParam(&penalty_offset,PenaltyOffsetEnum);
    577630
    578         //Create elementary matrix: add penalty to constrain wb (wb=ub*db/dx+vb*db/dy)
    579         Ke[2][0]=-slope[0]*kmax*pow((double)10.0,penalty_offset);
    580         Ke[2][1]=-slope[1]*kmax*pow((double)10.0,penalty_offset);
    581         Ke[2][2]=kmax*pow((double)10,penalty_offset);
    582        
    583         /*Add Ke to global matrix Kgg: */
    584         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
    585        
    586         /*Free ressources:*/
    587         xfree((void**)&doflist);
    588 
    589 }
    590 /*}}}1*/
    591 /*FUNCTION Pengrid::PenaltyCreateKMatrixMelting {{{1*/
    592 void  Pengrid::PenaltyCreateKMatrixMelting(Mat Kgg,double kmax){
    593 
    594         int found=0;
    595         const int numgrids=1;
    596         const int NDOF1=1;
    597         const int numdof=numgrids*NDOF1;
    598         double Ke[numdof][numdof]={0.0};
    599         int       dofs1[1]={0};
    600         int*      doflist=NULL;
    601         double    meltingpoint;
    602 
    603         double pressure;
    604         double temperature;
    605         double beta,t_pmp;
    606         double penalty_offset;
    607 
    608         /*recover pointers: */
    609         Penta* penta=(Penta*)element;
    610 
    611         /*check that pengrid is not a clone (penalty to be added only once)*/
    612         if (node->IsClone()) return;
    613 
    614         //First recover pressure and temperature values, using the element: */
    615         penta->GetParameterValue(&pressure,node,PressureEnum);
    616         penta->GetParameterValue(&temperature,node,TemperatureEnum);
    617 
    618         /*recover parameters: */
    619         parameters->FindParam(&penalty_offset,PenaltyOffsetEnum);
    620 
    621         /*Get dof list: */
    622         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    623        
    624         //Compute pressure melting point
    625         meltingpoint=matpar->GetMeltingPoint();
    626         beta=matpar->GetBeta();
    627         t_pmp=meltingpoint-beta*pressure;
    628 
    629         //Add penalty load
    630         if (temperature<t_pmp){ //If T<Tpmp, there must be no melting. Therefore, melting should be  constrained to 0 when T<Tpmp, instead of using spcs, use penalties
    631                 Ke[0][0]=kmax*pow((double)10,penalty_offset);
    632         }
    633        
    634         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
    635 
    636         /*Clean up*/
    637         xfree((void**)&doflist);
    638 }
    639 /*}}}1*/
    640 /*FUNCTION Pengrid::PenaltyCreateKMatrixThermal {{{1*/
    641 void  Pengrid::PenaltyCreateKMatrixThermal(Mat Kgg,double kmax){
    642 
    643         int found=0;
    644        
    645         const int numgrids=1;
    646         const int NDOF1=1;
    647         const int numdof=numgrids*NDOF1;
    648         double Ke[numdof][numdof];
    649         int*         doflist=NULL;
    650         double    penalty_offset;
    651 
    652         if(!this->active)return;
    653 
    654         /*recover parameters: */
    655         parameters->FindParam(&penalty_offset,PenaltyOffsetEnum);
    656 
    657         /*Get dof list: */
    658         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    659 
    660         Ke[0][0]=kmax*pow((double)10,penalty_offset);
    661        
    662         /*Add Ke to global matrix Kgg: */
    663         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
    664 
    665         /*Free ressources:*/
    666         xfree((void**)&doflist);
     631        Ke->values[0]=kmax*pow((double)10,penalty_offset);
     632
     633        /*Clean up and return*/
     634        return Ke;
    667635}
    668636/*}}}1*/
    669637/*FUNCTION Pengrid::PenaltyCreatePVectorMelting {{{1*/
    670 void  Pengrid::PenaltyCreatePVectorMelting(Vec pg, double kmax){
    671        
    672         const int numgrids=1;
    673         const int NDOF1=1;
    674         const int numdof=numgrids*NDOF1;
    675         int*   doflist=NULL;
    676         double P_terms[numdof]={0.0};
    677         int    found=0;
    678         int    dofs1[1]={0};
     638ElementVector* Pengrid::PenaltyCreatePVectorMelting(double kmax){
     639       
     640        const int numdof=NUMVERTICES*NDOF1;
    679641        double pressure;
    680642        double temperature;
    681643        double melting_offset;
    682         double meltingpoint;
    683         double beta, heatcapacity;
    684         double latentheat;
    685644        double t_pmp;
    686645        double dt,penalty_offset;
    … …  
    690649
    691650        /*check that pengrid is not a clone (penalty to be added only once)*/
    692         if (node->IsClone()) return;
    693 
    694         /*Get dof list: */
    695         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    696 
    697         //First recover pressure and temperature values, using the element: */
     651        if (node->IsClone()) return NULL;
     652        ElementVector* pe=NewElementVector(&node,NUMVERTICES,this->parameters);
     653
     654        /*Retrieve all inputs and parameters*/
    698655        penta->GetParameterValue(&pressure,node,PressureEnum);
    699656        penta->GetParameterValue(&temperature,node,TemperatureEnum);
    … …  
    702659        parameters->FindParam(&penalty_offset,PenaltyOffsetEnum);
    703660
    704         meltingpoint=matpar->GetMeltingPoint();
    705         beta=matpar->GetBeta();
    706         heatcapacity=matpar->GetHeatCapacity();
    707         latentheat=matpar->GetLatentHeat();
    708 
    709         //Compute pressure melting point
    710         t_pmp=meltingpoint-beta*pressure;
    711 
    712         //Add penalty load
    713         //This time, the penalty must have the same value as the one used for the thermal computation
    714         //so that the corresponding melting can be computed correctly
    715         //In the thermal computation, we used kmax=melting_offset, and the same penalty_offset
     661        /*Compute pressure melting point*/
     662        t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
     663
     664        /*Add penalty load
     665          This time, the penalty must have the same value as the one used for the thermal computation
     666          so that the corresponding melting can be computed correctly
     667          In the thermal computation, we used kmax=melting_offset, and the same penalty_offset*/
    716668        if (temperature<t_pmp){ //%no melting
    717                 P_terms[0]=0;
     669                pe->values[0]=0;
    718670        }
    719671        else{
    720                 if (dt){
    721                         P_terms[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp)/dt;
    722                 }
    723                 else{
    724                         P_terms[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp);
    725                 }
    726         }
    727 
    728         /*Add P_terms to global vector pg: */
    729         VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
    730 
    731         /*Clean up*/
    732         xfree((void**)&doflist);
    733 
     672                if (dt) pe->values[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp)/dt;
     673                else    pe->values[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp);
     674        }
     675
     676        /*Clean up and return*/
     677        return pe;
    734678}
    735679/*}}}1*/
    736680/*FUNCTION Pengrid::PenaltyCreatePVectorThermal {{{1*/
    737 void  Pengrid::PenaltyCreatePVectorThermal(Vec pg,  double kmax){
    738 
    739         const int numgrids=1;
    740         const int NDOF1=1;
    741         const int numdof=numgrids*NDOF1;
    742         int*      doflist=NULL;
    743         double  P_terms[numdof]={0.0};
    744         int    found=0;
     681ElementVector* Pengrid::PenaltyCreatePVectorThermal(double kmax){
     682
     683        const int numdof=NUMVERTICES*NDOF1;
    745684        double pressure;
    746         int dofs1[1]={0};
    747         double meltingpoint;
    748         double beta;
    749685        double t_pmp;
    750686        double penalty_offset;
    751687
    752         /*recover pointers: */
    753688        Penta* penta=(Penta*)element;
    754689
    755         if(!this->active)return;
    756 
    757         /*Get dof list: */
    758         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    759 
    760         //First recover pressure  and penalty_offset
     690        /*Initialize Element matrix and return if necessary*/
     691        if(!this->active) return NULL;
     692        ElementVector* pe=NewElementVector(&node,1,this->parameters);
     693
     694        /*Retrieve all inputs and parameters*/
    761695        penta->GetParameterValue(&pressure,node,PressureEnum);
    762696        parameters->FindParam(&penalty_offset,PenaltyOffsetEnum);
    763697
    764         //Compute pressure melting point
    765         meltingpoint=matpar->GetMeltingPoint();
    766         beta=matpar->GetBeta();
    767         t_pmp=meltingpoint-beta*pressure;
    768 
    769         //Add penalty load
    770         P_terms[0]=kmax*pow((double)10,penalty_offset)*t_pmp;
    771 
    772         /*Add P_terms to global vector pg: */
    773         VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
    774 
    775         /*Clean up*/
    776         xfree((void**)&doflist);
    777 
     698        /*Compute pressure melting point*/
     699        t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
     700
     701        pe->values[0]=kmax*pow((double)10,penalty_offset)*t_pmp;
     702
     703        /*Clean up and return*/
     704        return pe;
    778705}
    779706/*}}}1*/

  • issm/trunk/src/c/objects/Loads/Pengrid.h

    r5772 r5936  
    7878                /*}}}*/
    7979                /*Pengrid management {{{1*/
    80                 void  PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax);
    8180                void  GetDofList(int** pdoflist,int approximation_enum,int setenum);
    82                 void  PenaltyCreateKMatrixThermal(Mat Kgg,double kmax);
    83                 void  PenaltyCreateKMatrixMelting(Mat Kgg,double kmax);
    84                 void  PenaltyCreatePVectorThermal(Vec pg, double kmax);
    85                 void  PenaltyCreatePVectorMelting(Vec pg, double kmax);
     81                ElementMatrix* PenaltyCreateKMatrixDiagnosticStokes(double kmax);
     82                ElementMatrix* PenaltyCreateKMatrixThermal(double kmax);
     83                ElementMatrix* PenaltyCreateKMatrixMelting(double kmax);
     84                ElementVector* PenaltyCreatePVectorThermal(double kmax);
     85                ElementVector* PenaltyCreatePVectorMelting(double kmax);
    8686                void  PenaltyConstrain(int* punstable);
    8787                void  PenaltyConstrainThermal(int* punstable);
Note: See TracChangeset for help on using the changeset viewer.