Changeset 13925


Ignore:
Timestamp:
11/09/12 15:25:17 (12 years ago)
Author:
Mathieu Morlighem
Message:

NEW: Finished stiffness allocation matrix (still turned off), also, as of PETSC 3.3, no new nonzero can be added after matrix assembly. We now do a dry run to get Kff_max if penalties are present

Location:
issm/trunk-jpl/src/c
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Container/Loads.cpp

    r13916 r13925  
    5151        }
    5252
     53}
     54/*}}}*/
     55/*FUNCTION Loads::IsPenalty{{{*/
     56bool Loads::IsPenalty(int analysis_type){
     57
     58        int ispenalty=0;
     59        int allispenalty=0;
     60
     61        /*Now go through all loads, and get how many nodes they own, unless they are clone nodes: */
     62        for(int i=0;i<this->Size();i++){
     63
     64                Load* load=dynamic_cast<Load*>(this->GetObjectByOffset(i));
     65                if (load->InAnalysis(analysis_type)){
     66                        if(load->IsPenalty()) ispenalty++;
     67                }
     68        }
     69
     70        /*Grab sum of all cpus: */
     71#ifdef _HAVE_MPI_
     72        MPI_Allreduce((void*)&ispenalty,(void*)&allispenalty,1,MPI_INT,MPI_SUM,IssmComm::GetComm());
     73        ispenalty=allispenalty;
     74#endif
     75
     76        if(ispenalty)
     77         return true;
     78        else
     79         return false;
    5380}
    5481/*}}}*/

  • issm/trunk-jpl/src/c/Container/Loads.h

    r13916 r13925  
    2626                /*numerics*/
    2727                void  Configure(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters);
     28                bool  IsPenalty(int analysis);
    2829                int   MaxNumNodes(int analysis);
    2930                int   NumberOfLoads(void);

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r13918 r13925  
    353353        Vector<IssmDouble> *df   = NULL;
    354354
    355         //bool oldalloc=false;
    356355        bool oldalloc=true;
    357356
    … …  
    369368
    370369        if(oldalloc){
    371                 Kff=new Matrix<IssmDouble>(fsize,fsize,connectivity,numberofdofspernode);
    372                 Kfs=new Matrix<IssmDouble>(fsize,ssize,connectivity,numberofdofspernode);
    373                 df =new Vector<IssmDouble>(fsize);
    374                 pf =new Vector<IssmDouble>(fsize);
     370                if(pKff) Kff=new Matrix<IssmDouble>(fsize,fsize,connectivity,numberofdofspernode);
     371                if(pKfs) Kfs=new Matrix<IssmDouble>(fsize,ssize,connectivity,numberofdofspernode);
     372                if(pdf)  df =new Vector<IssmDouble>(fsize);
     373                if(ppf)  pf =new Vector<IssmDouble>(fsize);
    375374        }
    376375        else{
    377                 /*Kff*/
    378                 m=flocalsize; n=flocalsize; /*local  sizes*/
    379                 M=fsize;      N=fsize;      /*global sizes*/
    380                 this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,FsetEnum);
    381                 Kff=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
    382                 xDelete<int>(d_nnz);
    383                 xDelete<int>(o_nnz);
    384 
    385                 /*Kfs*/
    386                 m=flocalsize; n=slocalsize; /*local  sizes*/
    387                 M=fsize;      N=ssize;      /*global sizes*/
    388                 this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,SsetEnum);
    389                 Kfs=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
    390                 xDelete<int>(d_nnz);
    391                 xDelete<int>(o_nnz);
    392 
    393                 /*Vectors*/
    394                 df =new Vector<IssmDouble>(flocalsize,fsize);
    395                 pf =new Vector<IssmDouble>(flocalsize,fsize);
     376                if(pKff){
     377                        m=flocalsize; n=flocalsize; /*local  sizes*/
     378                        M=fsize;      N=fsize;      /*global sizes*/
     379                        this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,FsetEnum);
     380                        Kff=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
     381                        xDelete<int>(d_nnz);
     382                        xDelete<int>(o_nnz);
     383                }
     384                if(pKfs){
     385                        m=flocalsize; n=slocalsize; /*local  sizes*/
     386                        M=fsize;      N=ssize;      /*global sizes*/
     387                        this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,SsetEnum);
     388                        Kfs=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
     389                        xDelete<int>(d_nnz);
     390                        xDelete<int>(o_nnz);
     391                }
     392                if(pdf) df =new Vector<IssmDouble>(flocalsize,fsize);
     393                if(ppf) pf =new Vector<IssmDouble>(flocalsize,fsize);
    396394        }
    397395
    398396        /*Allocate output pointers*/
    399         *pKff = Kff;
    400         *pKfs = Kfs;
    401         *pdf  = df;
    402         *ppf  = pf;
     397        if(pKff) *pKff = Kff;
     398        if(pKfs) *pKfs = Kfs;
     399        if(pdf)  *pdf  = df;
     400        if(ppf)  *ppf  = pf;
    403401}
    404402/*}}}*/
    … …  
    770768        Vector<IssmDouble> *pf   = NULL;
    771769        Vector<IssmDouble> *df   = NULL;
    772         IssmDouble          kmax;
     770        IssmDouble          kmax = 0;
    773771
    774772        /*Display message*/
    … …  
    778776        this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    779777
    780         /*Compute penalty free mstiffness matrix and load vector*/
     778        /*First, we might need to do a dry run to get kmax if penalties are employed*/
     779        if(loads->IsPenalty(configuration_type)){
     780
     781                /*Allocate Kff_temp*/
     782                Matrix<IssmDouble> *Kff_temp = NULL;
     783                this->AllocateSystemMatrices(&Kff_temp,NULL,NULL,NULL);
     784
     785                /*Get complete stiffness matrix without penalties*/
     786                for (i=0;i<this->elements->Size();i++){
     787                        element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
     788                        element->CreateKMatrix(Kff_temp,NULL,NULL);
     789                }
     790                for (i=0;i<this->loads->Size();i++){
     791                        load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
     792                        if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff_temp,NULL);
     793                }
     794                Kff_temp->Assemble();
     795
     796                /*Now, figure out maximum value of stiffness matrix, so that we can penalize it correctly: */
     797                kmax=Kff_temp->Norm(NORM_INF);
     798                delete Kff_temp;
     799        }
     800
     801        /*Allocate stiffness matrices and load vector*/
    781802        this->AllocateSystemMatrices(&Kff,&Kfs,&df,&pf);
    782803
    … …  
    788809        for (i=0;i<this->loads->Size();i++){
    789810                load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
    790                 if (load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs);
     811                if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs);
    791812        }
    792813
    … …  
    798819        for (i=0;i<this->loads->Size();i++){
    799820                load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
    800                 if (load->InAnalysis(configuration_type)) load->CreatePVector(pf);
    801         }
    802 
    803         /*Assemble matrices (required to get kmax)*/
    804         Kff->Assemble();
    805         Kfs->Assemble();
    806         df->Assemble();
    807 
    808         /*Now, figure out maximum value of stiffness matrix, so that we can penalize it correctly: */
    809         kmax=Kff->Norm(NORM_INF);
    810 
    811         /*Now that we have kmax, deal with penalties (only in loads)*/
    812         for (i=0;i<this->loads->Size();i++){
    813                 load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
    814                 if (load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
    815         }
    816         for (i=0;i<this->loads->Size();i++){
    817                 load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
    818                 if (load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pf,kmax);
     821                if(load->InAnalysis(configuration_type)) load->CreatePVector(pf);
     822        }
     823
     824        /*Now deal with penalties (only in loads)*/
     825        if(loads->IsPenalty(configuration_type)){
     826                for (i=0;i<this->loads->Size();i++){
     827                        load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
     828                        if(load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
     829                }
     830                for (i=0;i<this->loads->Size();i++){
     831                        load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
     832                        if(load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pf,kmax);
     833                }
    819834        }
    820835
    … …  
    823838        Kfs->Assemble();
    824839        pf->Assemble();
     840        df->Assemble();
    825841        //Kff->AllocationInfo();
    826842        //Kfs->AllocationInfo();

  • issm/trunk-jpl/src/c/classes/objects/Loads/Icefront.cpp

    r13915 r13925  
    316316        }
    317317
     318}
     319/*}}}*/
     320/*FUNCTION Icefront::IsPenalty{{{*/
     321bool Icefront::IsPenalty(void){
     322        return false;
    318323}
    319324/*}}}*/

  • issm/trunk-jpl/src/c/classes/objects/Loads/Icefront.h

    r13915 r13925  
    7373                int   GetNumberOfNodes(void);
    7474                void  GetNodesSidList(int* sidlist);
     75                bool  IsPenalty(void);
    7576                void  PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* kfs, IssmDouble kmax);
    7677                void  PenaltyCreatePVector(Vector<IssmDouble>*  pf, IssmDouble kmax);

  • issm/trunk-jpl/src/c/classes/objects/Loads/Load.h

    r13915 r13925  
    2525                virtual       ~Load(){};
    2626                virtual void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
     27                virtual bool  IsPenalty(void)=0;
    2728                virtual int   GetNumberOfNodes(void)=0;
    2829                virtual void  GetNodesSidList(int* sidlist)=0;

  • issm/trunk-jpl/src/c/classes/objects/Loads/Numericalflux.cpp

    r13915 r13925  
    340340        }
    341341
     342}
     343/*}}}*/
     344/*FUNCTION Numericalflux::IsPenalty{{{*/
     345bool Numericalflux::IsPenalty(void){
     346        return false;
    342347}
    343348/*}}}*/

  • issm/trunk-jpl/src/c/classes/objects/Loads/Numericalflux.h

    r13915 r13925  
    6767                int   GetNumberOfNodes(void);
    6868                void  CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("Not implemented yet");};
     69                bool  IsPenalty(void);
    6970                void  PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){_error_("Not implemented yet");};
    7071                void  PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* kfs, IssmDouble kmax);

  • issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.cpp

    r13915 r13925  
    293293        if (in_analysis_type==this->analysis_type)return true;
    294294        else return false;
     295}
     296/*}}}*/
     297/*FUNCTION Pengrid::IsPenalty{{{*/
     298bool Pengrid::IsPenalty(void){
     299        return true;
    295300}
    296301/*}}}*/

  • issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.h

    r13915 r13925  
    7474                void  GetNodesSidList(int* sidlist);
    7575                int   GetNumberOfNodes(void);
     76                bool  IsPenalty(void);
    7677                void  PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){_error_("Not implemented yet");};
    7778                void  PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* kfs, IssmDouble kmax);

  • issm/trunk-jpl/src/c/classes/objects/Loads/Penpair.cpp

    r13915 r13925  
    161161
    162162        return NUMVERTICES;
     163}
     164/*}}}*/
     165/*FUNCTION Penpair::IsPenalty{{{*/
     166bool Penpair::IsPenalty(void){
     167        return true;
    163168}
    164169/*}}}*/

  • issm/trunk-jpl/src/c/classes/objects/Loads/Penpair.h

    r13915 r13925  
    5959                void  GetNodesSidList(int* sidlist);
    6060                int   GetNumberOfNodes(void);
     61                bool  IsPenalty(void);
    6162                void  PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff,Matrix<IssmDouble>* Kfs,IssmDouble kmax);
    6263                void  PenaltyCreatePVector(Vector<IssmDouble>* pf, IssmDouble kmax);

  • issm/trunk-jpl/src/c/classes/objects/Loads/Riftfront.cpp

    r13915 r13925  
    293293        this->parameters=parametersin;
    294294
     295}
     296/*}}}*/
     297/*FUNCTION Riftfront::IsPenalty{{{*/
     298bool Riftfront::IsPenalty(void){
     299        return true;
    295300}
    296301/*}}}*/

  • issm/trunk-jpl/src/c/classes/objects/Loads/Riftfront.h

    r13915 r13925  
    8080                void  GetNodesSidList(int* sidlist);
    8181                int   GetNumberOfNodes(void);
     82                bool  IsPenalty(void);
    8283                void  PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){_error_("Not implemented yet");};
    8384                void  PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* kfs, IssmDouble kmax);

  • issm/trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp

    r13622 r13925  
    1212
    1313        /*variables: */
    14         int i;
    15         int configuration_type;
    16         int fsize;
    17         IssmDouble* ug_serial=NULL;
     14        int         configuration_type;
     15        int         fsize;
     16        IssmDouble *ug_serial = NULL;
     17        bool        oldalloc  = true;
    1818
    1919        /*first figure out fsize: */
    2020        parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    2121        fsize=nodes->NumberOfDofs(configuration_type,FsetEnum);
     22
     23        int    analysis_type;
     24        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     25        int flocalsize = nodes->NumberOfDofsLocal(analysis_type,FsetEnum);
    2226
    2327        if(fsize==0){
    … …  
    2630        else{
    2731                /*allocate: */
    28                 uf=new Vector<IssmDouble>(fsize);
     32                if(oldalloc)
     33                 uf=new Vector<IssmDouble>(fsize);
     34                else
     35                 uf=new Vector<IssmDouble>(flocalsize,fsize);
    2936
    3037                if(nodes->NumberOfNodes(configuration_type)){
    … …  
    3441
    3542                        /*Go through all nodes, and ask them to retrieve values from ug, and plug them into uf: */
    36                         for(i=0;i<nodes->Size();i++){
     43                        for(int i=0;i<nodes->Size();i++){
    3744
    3845                                Node* node=(Node*)nodes->GetObjectByOffset(i);

  • issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp

    r13893 r13925  
    4040}
    4141/*}}}*/
    42 /*FUNCTION PetscMat::PetscMat(int M,int N, IssmDouble sparsity){{{*/
     42/*FUNCTION PetscMat::PetscMat(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/
    4343PetscMat::PetscMat(int m,int n,int M,int N,int* d_nnz,int* o_nnz){
    4444
Note: See TracChangeset for help on using the changeset viewer.