Changeset 5911


Ignore:
Timestamp:
09/20/10 15:34:55 (15 years ago)
Author:
Mathieu Morlighem
Message:

Extended ElementMatrix to loads

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

Legend:

Unmodified
Added
Removed

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

    r5909 r5911  
    318318
    319319        /*Checks in debugging mode*/
     320        /*{{{2*/
    320321        ISSMASSERT(nodes);
    321322        ISSMASSERT(element);
    322323        ISSMASSERT(matpar);
     324        /*}}}*/
    323325
    324326        /*Retrieve parameters: */

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

    r5772 r5911  
    332332        int type;
    333333        int analysis_type;
     334        ElementMatrix* Ke=NULL;
    334335
    335336        /*recover some parameters*/
    … …  
    337338        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    338339
    339         if (analysis_type==PrognosticAnalysisEnum || analysis_type==BalancedthicknessAnalysisEnum || analysis_type==AdjointBalancedthicknessAnalysisEnum){
    340                 if      (type==InternalEnum) CreateKMatrixInternal(Kgg);
    341                 else if (type==BoundaryEnum) CreateKMatrixBoundary(Kgg);
    342                 else ISSMERROR("type not supported yet");
    343         }
    344         else{
    345                 ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type));
     340        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
     341        switch(analysis_type){
     342                case PrognosticAnalysisEnum: case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum:
     343                        switch(type){
     344                                case InternalEnum:
     345                                        Ke=CreateKMatrixInternal();
     346                                        break;
     347                                case BoundaryEnum:
     348                                        Ke=CreateKMatrixBoundary();
     349                                        break;
     350                                default:
     351                                        ISSMERROR("type not supported yet");
     352                        }
     353                        break;
     354                default:
     355                        ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
     356        }
     357
     358        /*Add to global matrix*/
     359        if(Ke){
     360                Ke->AddToGlobal(Kgg,Kff,Kfs);
     361                delete Ke;
    346362        }
    347363
    … …  
    353369        int type;
    354370        int analysis_type;
     371        ElementVector* pe=NULL;
    355372
    356373        /*recover some parameters*/
    … …  
    358375        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    359376
    360         if (analysis_type==PrognosticAnalysisEnum || analysis_type==BalancedthicknessAnalysisEnum){
    361                 if      (type==InternalEnum) CreatePVectorInternal(pg);
    362                 else if (type==BoundaryEnum) CreatePVectorBoundary(pg);
    363                 else ISSMERROR("type not supported yet");
    364         }
    365         else if (analysis_type==AdjointBalancedthicknessAnalysisEnum){
    366                 return;
    367         }
    368         else{
    369                 ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type));
    370         }
    371 
     377        switch(analysis_type){
     378                case PrognosticAnalysisEnum: case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum:
     379                        switch(type){
     380                                case InternalEnum:
     381                                        pe=CreatePVectorInternal();
     382                                        break;
     383                                case BoundaryEnum:
     384                                        pe=CreatePVectorBoundary();
     385                                        break;
     386                                default:
     387                                        ISSMERROR("type not supported yet");
     388                        }
     389                        break;
     390                default:
     391                        ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
     392        }
     393
     394        /*Add to global matrix*/
     395        if(pe){
     396                pe->AddToGlobal(pg,pf);
     397                delete pe;
     398        }
    372399
    373400}
    … …  
    398425/*Numericalflux management*/
    399426/*FUNCTION Numericalflux::CreateKMatrixInternal {{{1*/
    400 void  Numericalflux::CreateKMatrixInternal(Mat Kgg){
     427ElementMatrix* Numericalflux::CreateKMatrixInternal(void){
    401428
    402429        /* constants*/
    … …  
    412439        double     Ke_g1[numdof][numdof];
    413440        double     Ke_g2[numdof][numdof];
    414         double     Ke[numdof][numdof] = {0.0};
    415         int       *doflist = NULL;
    416441        GaussTria *gauss;
    417442
     443        /*Initialize Element matrix and return if necessary*/
     444        Tria*  tria=(Tria*)element;
     445        if(tria->IsOnWater()) return NULL;
     446        ElementMatrix* Ke=NewElementMatrix(nodes,NUMVERTICES_INTERNAL,this->parameters);
     447
     448        /*Retrieve all inputs and parameters*/
    418449        GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES_INTERNAL);
    419         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    420 
    421         /*Retrieve all inputs and parameters we will be needing: */
    422         Tria*  tria=(Tria*)element;
    423450        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
    424451        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
    … …  
    461488                                        &Ke_g2[0][0],0);
    462489
    463                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke[i][j]+=Ke_g1[i][j];
    464                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke[i][j]+=Ke_g2[i][j];
    465         }
    466 
    467         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
     490                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g1[i][j];
     491                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g2[i][j];
     492        }
    468493       
    469         /*Free ressources:*/
     494        /*Clean up and return*/
    470495        delete gauss;
    471         xfree((void**)&doflist);
    472 
     496        return Ke;
    473497}
    474498/*}}}*/
    475499/*FUNCTION Numericalflux::CreateKMatrixBoundary {{{1*/
    476 void  Numericalflux::CreateKMatrixBoundary(Mat Kgg){
     500ElementMatrix* Numericalflux::CreateKMatrixBoundary(void){
    477501
    478502        /* constants*/
    … …  
    486510        double     L[numdof];
    487511        double     Ke_g[numdof][numdof];
    488         double     Ke[numdof][numdof] = {0.0};
    489         int       *doflist = NULL;
    490512        GaussTria *gauss;
    491513
     514        /*Initialize Element matrix and return if necessary*/
     515        Tria*  tria=(Tria*)element;
     516        if(tria->IsOnWater()) return NULL;
     517        ElementMatrix* Ke=NewElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters);
     518
     519        /*Retrieve all inputs and parameters*/
    492520        GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES_BOUNDARY);
    493 
    494         /*Retrieve all inputs and parameters we will be needing: */
    495         Tria*  tria=(Tria*)element;
    496521        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
    497522        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
    … …  
    520545        if (UdotN<=0){
    521546                /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
    522                 return;
    523         }
    524 
    525         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     547                return NULL;
     548        }
    526549
    527550        /* Start  looping on the number of gaussian points: */
    … …  
    544567                                        &Ke_g[0][0],0);
    545568
    546                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke[i][j]+=Ke_g[i][j];
     569                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
    547570        }
    548571
    549         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
    550 
    551         /*Free ressources:*/
     572        /*Clean up and return*/
    552573        delete gauss;
    553         xfree((void**)&doflist);
     574        return Ke;
    554575}
    555576/*}}}*/
    556577/*FUNCTION Numericalflux::CreatePVectorInternal{{{1*/
    557 void  Numericalflux::CreatePVectorInternal(Vec pg){
     578ElementVector* Numericalflux::CreatePVectorInternal(void){
    558579
    559580        /*Nothing added to PVector*/
    560         return;
     581        return NULL;
    561582
    562583}
    563584/*}}}*/
    564585/*FUNCTION Numericalflux::CreatePVectorBoundary{{{1*/
    565 void  Numericalflux::CreatePVectorBoundary(Vec pg){
     586ElementVector* Numericalflux::CreatePVectorBoundary(void){
    566587
    567588        /* constants*/
    … …  
    574595        double     normal[2];
    575596        double     L[numdof];
    576         double     pe[numdof] = {0.0};
    577         int       *doflist = NULL;
    578597        GaussTria *gauss;
    579598
     599        /*Initialize Load Vectorand return if necessary*/
     600        Tria*  tria=(Tria*)element;
     601        if(tria->IsOnWater()) return NULL;
     602        ElementVector* pe=NewElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);
     603
     604        /*Retrieve all inputs and parameters*/
    580605        GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES_BOUNDARY);
    581 
    582         /*Retrieve all inputs and parameters we will be needing: */
    583         Tria*  tria=(Tria*)element;
    584606        Input* vxaverage_input=tria->inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input);
    585607        Input* vyaverage_input=tria->inputs->GetInput(VyEnum); ISSMASSERT(vyaverage_input);
    … …  
    614636        if (UdotN>0){
    615637                /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
    616                 return;
    617         }
    618 
    619         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     638                return NULL;
     639        }
    620640
    621641        /* Start  looping on the number of gaussian points: */
    … …  
    634654                DL= - gauss->weight*Jdet*dt*UdotN*thickness;
    635655
    636                 for(i=0;i<numdof;i++) pe[i] += DL*L[i];
    637         }
    638 
    639         /*Add pe_g to global matrix Kgg: */
    640         VecSetValues(pg,numdof,doflist,(const double*)pe,ADD_VALUES);
    641 
    642         /*Free ressources:*/
     656                for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
     657        }
     658
     659        /*Clean up and return*/
    643660        delete gauss;
    644         xfree((void**)&doflist);
    645 }
    646 /*}}}*/
    647 /*FUNCTION Numericalflux::GetDofList {{{1*/
    648 void  Numericalflux::GetDofList(int** pdoflist,int approximation,int setenum){
    649 
    650         int i,j;
    651         int numberofdofs=0;
    652         int count=0;
    653         int type;
    654         int numberofnodes;
    655         int* doflist=NULL;
    656        
    657         /*Some checks for debugging*/
    658         ISSMASSERT(nodes);
    659                
    660         /*How many nodes? :*/
    661         inputs->GetParameterValue(&type,TypeEnum);
    662         switch(type){
    663                 case InternalEnum:
    664                         numberofnodes=NUMVERTICES_INTERNAL; break;
    665                 case BoundaryEnum:
    666                         numberofnodes=NUMVERTICES_BOUNDARY; break;
    667                 default:
    668                         ISSMERROR("type %s not supported yet",type);
    669         }
    670        
    671         /*Figure out size of doflist: */
    672         for(i=0;i<numberofnodes;i++){
    673                 numberofdofs+=nodes[i]->GetNumberOfDofs(approximation,setenum);
    674         }
    675 
    676         /*Allocate: */
    677         doflist=(int*)xmalloc(numberofdofs*sizeof(int));
    678 
    679         /*Populate: */
    680         count=0;
    681         for(i=0;i<numberofnodes;i++){
    682                 nodes[i]->GetDofList(doflist+count,approximation,setenum);
    683                 count+=nodes[i]->GetNumberOfDofs(approximation,setenum);
    684         }
    685 
    686         /*Assign output pointers:*/
    687         *pdoflist=doflist;
     661        return pe;
    688662}
    689663/*}}}*/

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

    r5772 r5911  
    1313class Inputs;
    1414class IoModel;
    15 
     15class ElementMatrix;
     16class ElementVector;
    1617/*}}}*/
    1718
    … …  
    7273                /*}}}*/
    7374                /*Numericalflux management:{{{1*/
    74                 void  GetDofList(int** pdoflist,int approximation_enum,int setenum);
    7575                void  GetNormal(double* normal,double xyz_list[4][3]);
    76                 void  CreateKMatrixInternal(Mat Kgg);
    77                 void  CreateKMatrixBoundary(Mat Kgg);
    78                 void  CreatePVectorInternal(Vec pg);
    79                 void  CreatePVectorBoundary(Vec pg);
     76                ElementMatrix* CreateKMatrixInternal(void);
     77                ElementMatrix* CreateKMatrixBoundary(void);
     78                ElementVector* CreatePVectorInternal(void);
     79                ElementVector* CreatePVectorBoundary(void);
    8080                /*}}}*/
    8181

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

    r5772 r5911  
    378378void  Riftfront::CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs){
    379379        /*do nothing: */
     380        return;
    380381}
    381382/*}}}1*/
Note: See TracChangeset for help on using the changeset viewer.