Changeset 16789


Ignore:
Timestamp:
11/15/13 15:41:42 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: AddToGlobal and condensation now done by SystemMatrices

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r16787 r16789  
    1111/*{{{*/
    1212#include "../../datastructures/datastructures.h"
     13#include "../../toolkits/toolkits.h"
    1314#include "../Update.h"
    14 
    1515class DataSet;
    1616class Parameters;
    … …  
    2626template <class doublematrix> class Matrix;
    2727template <class doubletype> class Vector;
    28 
    29 #include "../../toolkits/toolkits.h"
     28class ElementMatrix;
     29class ElementVector;
    3030/*}}}*/
    3131
    … …  
    4343                virtual void   AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
    4444                virtual void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>*  Kfs)=0;
     45                virtual ElementMatrix* CreateKMatrix(void)=0;
    4546                virtual void   CreateDVector(Vector<IssmDouble>* df)=0;
    4647                virtual void   CreatePVector(Vector<IssmDouble>* pf)=0;
     48                virtual ElementVector* CreatePVector(void)=0;
    4749                virtual void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0;
    4850                virtual void   DeleteMaterials(void)=0;
    … …  
    110112           virtual Element*   SpawnBasalElement(void)=0;
    111113                virtual IssmDouble TMeltingPoint(IssmDouble pressure)=0;
     114                virtual void   ReduceMatrices(ElementMatrix* Ke,ElementVector* pe)=0;
    112115                virtual void   ResetCoordinateSystem()=0;
    113116                virtual int    VelocityInterpolation()=0;

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r16788 r16789  
    26522652        /*clean-up*/
    26532653        delete gauss;
     2654}
     2655/*}}}*/
     2656/*FUNCTION Penta::ReduceMatrices) {{{*/
     2657void Penta::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){
     2658
     2659        int analysis_type;
     2660        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     2661
     2662        if(pe){
     2663                if(analysis_type==StressbalanceAnalysisEnum){
     2664                        if(this->element_type==MINIcondensedEnum){
     2665                                int approximation;
     2666                                inputs->GetInputValue(&approximation,ApproximationEnum);
     2667                                if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
     2668                                        //Do nothing, condensation already done in PVectorCoupling
     2669                                }
     2670                                else{
     2671                                        int indices[3]={18,19,20};
     2672                                        pe->StaticCondensation(Ke,3,&indices[0]);
     2673                                }
     2674                        }
     2675                        else if(this->element_type==P1bubblecondensedEnum){
     2676                                int size   = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     2677                                int offset = 0;
     2678                                for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     2679                                int* indices=xNew<int>(size);
     2680                                for(int i=0;i<size;i++) indices[i] = offset+i;
     2681                                pe->StaticCondensation(Ke,size,indices);
     2682                                xDelete<int>(indices);
     2683                                delete Ke;
     2684                        }
     2685                }
     2686        }
     2687
     2688        if(Ke){
     2689                if(analysis_type==StressbalanceAnalysisEnum){
     2690                        int approximation;
     2691                        inputs->GetInputValue(&approximation,ApproximationEnum);
     2692                        if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
     2693                                //Do nothing condensatino already done for Stokes part
     2694                        }
     2695                        else{
     2696                                if(this->element_type==MINIcondensedEnum){
     2697                                        int indices[3]={18,19,20};
     2698                                        Ke->StaticCondensation(3,&indices[0]);
     2699                                }
     2700                                else if(this->element_type==P1bubblecondensedEnum){
     2701                                        int size   = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     2702                                        int offset = 0;
     2703                                        for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     2704                                        int* indices=xNew<int>(size);
     2705                                        for(int i=0;i<size;i++) indices[i] = offset+i;
     2706                                        Ke->StaticCondensation(size,indices);
     2707                                        xDelete<int>(indices);
     2708                                }
     2709                        }
     2710                }
     2711        }
     2712
    26542713}
    26552714/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r16787 r16789  
    115115                void   ResultToVector(Vector<IssmDouble>* vector,int output_enum);
    116116                void   PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
     117                void   ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
    117118                void   ResetCoordinateSystem(void);
    118119                void   SmbGradients();

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r16787 r16789  
    7878                void        SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){_error_("not implemented yet");};
    7979                void        CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){_error_("not implemented yet");};
     80                ElementMatrix* CreateKMatrix(void){_error_("not implemented yet");};
    8081                void        CreateDVector(Vector<IssmDouble>* df){_error_("not implemented yet");};
    8182                void        CreatePVector(Vector<IssmDouble>* pf){_error_("not implemented yet");};
     83                ElementVector* CreatePVector(void){_error_("not implemented yet");};
    8284                ElementVector* CreatePVectorL2Projection(void);
    8385                void        CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("not implemented yet");};
    … …  
    147149                void        ResultToVector(Vector<IssmDouble>* vector,int output_enum){_error_("not implemented");};
    148150                void        ResetCoordinateSystem(void){_error_("not implemented yet");};
     151                void        ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){_error_("not implemented yet");};
    149152                void          SmbGradients(){_error_("not implemented yet");};
    150153                IssmDouble  TMeltingPoint(IssmDouble pressure){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r16788 r16789  
    23492349        /*clean-up*/
    23502350        delete gauss;
     2351}
     2352/*}}}*/
     2353/*FUNCTION Tria::ReduceMatrices {{{*/
     2354void Tria::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){
     2355
     2356        /*Static condensation if requested*/
     2357        if(pe){
     2358                if(this->element_type==MINIcondensedEnum){
     2359                        int indices[2]={6,7};
     2360                        pe->StaticCondensation(Ke,2,&indices[0]);
     2361                }
     2362                else if(this->element_type==P1bubblecondensedEnum){
     2363                        int size   = nodes[3]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     2364                        int offset = 0;
     2365                        for(int i=0;i<3;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     2366                        int* indices=xNew<int>(size);
     2367                        for(int i=0;i<size;i++) indices[i] = offset+i;
     2368                        pe->StaticCondensation(Ke,size,indices);
     2369                        xDelete<int>(indices);
     2370                }
     2371        }
     2372
     2373        if(Ke){
     2374                if(this->element_type==MINIcondensedEnum){
     2375                        int indices[2]={6,7};
     2376                        Ke->StaticCondensation(2,&indices[0]);
     2377                }
     2378                else if(this->element_type==P1bubblecondensedEnum){
     2379                        int size   = nodes[3]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     2380                        int offset = 0;
     2381                        for(int i=0;i<3;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     2382                        int* indices=xNew<int>(size);
     2383                        for(int i=0;i<size;i++) indices[i] = offset+i;
     2384                        Ke->StaticCondensation(size,indices);
     2385                        xDelete<int>(indices);
     2386                }
     2387        }
     2388
     2389
    23512390}
    23522391/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r16787 r16789  
    121121                void        ResultInterpolation(int* pinterpolation,int output_enum);
    122122                void        ResultToVector(Vector<IssmDouble>* vector,int output_enum);
     123                void        ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
    123124                void        ResetCoordinateSystem(void);
    124125                void          SmbGradients();

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

    r16126 r16789  
    1212        /*intermediary: */
    1313        int      i,M,N;
    14         int      configuration_type;
     14        int      configuration_type,analysisenum;
    1515        Element *element = NULL;
    1616        Load    *load    = NULL;
    … …  
    2626        if(VerboseModule()) _printf0_("   Generating matrices");
    2727
    28         /*retrive parameters: */
     28        /*retrieve parameters: */
    2929        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
     30        femmodel->parameters->FindParam(&analysisenum,AnalysisTypeEnum);
     31        Analysis* analysis = EnumToAnalysis(analysisenum);
    3032
    3133        /*First, we might need to do a dry run to get kmax if penalties are employed*/
    … …  
    6264        }
    6365
    64         /*Fill stiffness matrix from elements and loads */
     66        /*Fill stiffness matrix and load vector from elements*/
    6567        for (i=0;i<femmodel->elements->Size();i++){
    6668                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
    67                 element->CreateKMatrix(Kff,Kfs);
     69                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
     70                //ElementVector* pe = analysis->CreatePVector(element);
     71                //ElementMatrix* Ke = analysis->CreateKMatrix(element);
     72                ElementVector* pe = element->CreatePVector();
     73                ElementMatrix* Ke = element->CreateKMatrix();
     74                element->ReduceMatrices(Ke,pe);
     75                if(Ke) Ke->AddToGlobal(Kff,Kfs);
     76                if(pe) pe->AddToGlobal(pf);
     77                delete Ke;
     78                delete pe;
    6879        }
    6980
    70         for (i=0;i<femmodel->loads->Size();i++){
     81        /*Fill stiffness matrix and load vector from loads*/
     82        for(i=0;i<femmodel->loads->Size();i++){
    7183                load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
    72                 if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs);
    73         }
    74 
    75         /*Fill right hand side vector, from elements and loads */
    76         for (i=0;i<femmodel->elements->Size();i++){
    77                 element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
    78                 element->CreatePVector(pf);
    79         }
    80         for (i=0;i<femmodel->loads->Size();i++){
    81                 load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
    82                 if(load->InAnalysis(configuration_type)) load->CreatePVector(pf);
     84                if(load->InAnalysis(configuration_type)){
     85                        load->CreateKMatrix(Kff,Kfs);
     86                        load->CreatePVector(pf);
     87                }
    8388        }
    8489
    … …  
    8792                for (i=0;i<femmodel->loads->Size();i++){
    8893                        load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
    89                         if(load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
    90                 }
    91                 for (i=0;i<femmodel->loads->Size();i++){
    92                         load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
    93                         if(load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pf,kmax);
     94                        if(load->InAnalysis(configuration_type)){
     95                                load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
     96                                load->PenaltyCreatePVector(pf,kmax);
     97                        }
    9498                }
    9599        }
    … …  
    109113        //Kfs->AllocationInfo();
    110114
    111         /*Assign output pointers: */
     115        /*cleanu up and assign output pointers: */
     116        delete analysis;
    112117        if(pKff) *pKff=Kff;
    113118        else      delete Kff;

  • issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h

    r16126 r16789  
    66
    77#include "../../classes/classes.h"
     8#include "../../analyses/analyses.h"
    89
    910/* local prototypes: */
Note: See TracChangeset for help on using the changeset viewer.