Changeset 16805


Ignore:
Timestamp:
11/16/13 14:27:12 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: working on FS

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp

    r16782 r16805  
    5959}/*}}}*/
    6060ElementVector* L2ProjectionBaseAnalysis::CreatePVector(Element* element){/*{{{*/
    61 _error_("not implemented yet");
     61
     62        /*Intermediaries*/
     63        int      meshtype;
     64        Element* basalelement;
     65
     66        /*Get basal element*/
     67        element->FindParam(&meshtype,MeshTypeEnum);
     68        switch(meshtype){
     69                case Mesh2DhorizontalEnum:
     70                        basalelement = element;
     71                        break;
     72                case Mesh3DEnum:
     73                        if(!element->IsOnBed()) return NULL;
     74                        basalelement = element->SpawnBasalElement();
     75                        break;
     76                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     77        }
     78
     79        /*Intermediaries */
     80        int         input_enum;
     81        IssmDouble  Jdet,value,slopes[2];
     82        Input      *input     = NULL;
     83        Input      *input2    = NULL;
     84        IssmDouble *xyz_list  = NULL;
     85
     86        /*Fetch number of nodes and dof for this finite element*/
     87        int numnodes = basalelement->GetNumberOfNodes();
     88
     89        /*Initialize Element vector*/
     90        ElementVector* pe    = basalelement->NewElementVector();
     91        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     92
     93        /*Retrieve all inputs and parameters*/
     94        basalelement->GetVerticesCoordinates(&xyz_list);
     95        basalelement->FindParam(&input_enum,InputToL2ProjectEnum);
     96        switch(input_enum){
     97                case SurfaceSlopeXEnum: input2 = basalelement->GetInput(SurfaceEnum); _assert_(input2); break;
     98                case SurfaceSlopeYEnum: input2 = basalelement->GetInput(SurfaceEnum); _assert_(input2); break;
     99                case BedSlopeXEnum:     input2 = basalelement->GetInput(BedEnum);     _assert_(input2); break;
     100                case BedSlopeYEnum:     input2 = basalelement->GetInput(BedEnum);     _assert_(input2); break;
     101                default: input = element->GetInput(input_enum);
     102        }
     103
     104        /* Start  looping on the number of gaussian points: */
     105        Gauss* gauss=basalelement->NewGauss(2);
     106        for(int ig=gauss->begin();ig<gauss->end();ig++){
     107                gauss->GaussPoint(ig);
     108
     109                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     110                basalelement->NodalFunctions(basis,gauss);
     111
     112                if(input2) input2->GetInputDerivativeValue(&slopes[0],xyz_list,gauss);
     113                switch(input_enum){
     114                        case SurfaceSlopeXEnum: case BedSlopeXEnum: value = slopes[0]; break;
     115                        case SurfaceSlopeYEnum: case BedSlopeYEnum: value = slopes[1]; break;
     116                        default: input->GetInputValue(&value,gauss);
     117                }
     118
     119                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*value*basis[i];
     120        }
     121
     122        /*Clean up and return*/
     123        xDelete<IssmDouble>(xyz_list);
     124        xDelete<IssmDouble>(basis);
     125        delete gauss;
     126        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     127        return pe;
    62128}/*}}}*/
    63129void L2ProjectionBaseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r16804 r16805  
    814814                case HOApproximationEnum:
    815815                        return CreatePVectorHO(element);
     816                case FSApproximationEnum:
     817                        return CreatePVectorFS(element);
     818                case NoneApproximationEnum:
     819                        return NULL;
    816820                default:
    817821                        _error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
    818822        }
     823}/*}}}*/
     824ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
     825
     826        /*compute all load vectors for this element*/
     827        ElementVector* pe1=CreatePVectorFSViscous(element);
     828        ElementVector* pe2=CreatePVectorFSShelf(element);
     829        ElementVector* pe3=CreatePVectorFSFront(element);
     830        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
     831
     832        /*clean-up and return*/
     833        delete pe1;
     834        delete pe2;
     835        delete pe3;
     836        return pe;
     837}/*}}}*/
     838ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/
     839
     840        int         i,meshtype,dim;
     841        IssmDouble  Jdet,forcex,forcey,forcez;
     842        IssmDouble *xyz_list = NULL;
     843
     844        /*Get problem dimension*/
     845        element->FindParam(&meshtype,MeshTypeEnum);
     846        switch(meshtype){
     847                case Mesh2DverticalEnum: dim = 2; break;
     848                case Mesh3DEnum:         dim = 3; break;
     849                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     850        }
     851
     852        /*Fetch number of nodes and dof for this finite element*/
     853        int vnumnodes = element->GetNumberOfNodesVelocity();
     854        int pnumnodes = element->GetNumberOfNodesPressure();
     855
     856        /*Prepare coordinate system list*/
     857        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     858        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     859        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     860        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     861
     862        /*Initialize vectors*/
     863        ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
     864        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     865
     866        /*Retrieve all inputs and parameters*/
     867        element->GetVerticesCoordinates(&xyz_list);
     868        Input*      loadingforcex_input=element->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
     869        Input*      loadingforcey_input=element->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
     870        Input*      loadingforcez_input=element->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
     871        IssmDouble  rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
     872        IssmDouble  gravity =element->GetMaterialParameter(ConstantsGEnum);
     873
     874        /* Start  looping on the number of gaussian points: */
     875        Gauss* gauss=element->NewGauss(5);
     876        for(int ig=gauss->begin();ig<gauss->end();ig++){
     877                gauss->GaussPoint(ig);
     878
     879                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     880                element->NodalFunctionsVelocity(vbasis,gauss);
     881
     882                loadingforcex_input->GetInputValue(&forcex,gauss);
     883                loadingforcey_input->GetInputValue(&forcey,gauss);
     884                if(dim==3) loadingforcez_input->GetInputValue(&forcez,gauss);
     885
     886                for(i=0;i<vnumnodes;i++){
     887                        pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
     888                        pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
     889                        if(dim==3){
     890                                pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];
     891                                pe->values[i*dim+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
     892                        }
     893                        else{
     894                                pe->values[i*dim+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
     895                        }
     896                }
     897        }
     898
     899        /*Transform coordinate system*/
     900        element->TransformLoadVectorCoord(pe,cs_list);
     901
     902        /*Clean up and return*/
     903        delete gauss;
     904        xDelete<int>(cs_list);
     905        xDelete<IssmDouble>(vbasis);
     906        xDelete<IssmDouble>(xyz_list);
     907        return pe;
     908}/*}}}*/
     909ElementVector* StressbalanceAnalysis::CreatePVectorFSShelf(Element* element){/*{{{*/
     910
     911        return NULL;
     912}/*}}}*/
     913ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/
     914
     915        return NULL;
    819916}/*}}}*/
    820917ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/
    … …  
    9161013
    9171014        /*Initialize Element vector and vectors*/
    918         ElementVector* pe=element->NewElementVector(SSAApproximationEnum);
     1015        ElementVector* pe    = element->NewElementVector(SSAApproximationEnum);
    9191016        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    9201017
    … …  
    11461243        /*Prepare coordinate system list*/
    11471244        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    1148         if(dim==2){
    1149                 for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
    1150         }
    1151         else{
    1152                 for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
    1153         }
     1245        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     1246        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
    11541247        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    11551248

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h

    r16804 r16805  
    2323                ElementMatrix* CreateKMatrix(Element* element);
    2424                ElementVector* CreatePVector(Element* element);
     25                ElementVector* CreatePVectorFS(Element* element);
     26                ElementVector* CreatePVectorFSViscous(Element* element);
     27                ElementVector* CreatePVectorFSShelf(Element* element);
     28                ElementVector* CreatePVectorFSFront(Element* element);
    2529                ElementVector* CreatePVectorHO(Element* element);
    2630                ElementVector* CreatePVectorHODrivingStress(Element* element);

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

    r16803 r16805  
    5454                virtual IssmDouble GetMaterialParameter(int enum_in)=0;
    5555                virtual void   NodalFunctions(IssmDouble* basis,Gauss* gauss)=0;
     56                virtual void   NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss)=0;
     57                virtual void   NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss)=0;
    5658                virtual void   TransformLoadVectorCoord(ElementVector* pe,int transformenum)=0;
    5759                virtual void   TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list)=0;

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

    r16803 r16805  
    25432543        _assert_(gauss->Enum()==GaussPentaEnum);
    25442544        this->GetNodalFunctions(basis,(GaussPenta*)gauss);
     2545
     2546}
     2547/*}}}*/
     2548/*FUNCTION Penta::NodalFunctionsVelocity{{{*/
     2549void Penta::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){
     2550
     2551        _assert_(gauss->Enum()==GaussPentaEnum);
     2552        this->GetNodalFunctionsVelocity(basis,(GaussPenta*)gauss);
     2553
     2554}
     2555/*}}}*/
     2556/*FUNCTION Penta::NodalFunctionsPressure{{{*/
     2557void Penta::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){
     2558
     2559        _assert_(gauss->Enum()==GaussPentaEnum);
     2560        this->GetNodalFunctionsPressure(basis,(GaussPenta*)gauss);
    25452561
    25462562}

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

    r16803 r16805  
    247247                ElementVector* NewElementVector(int approximation_enum);
    248248                void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
     249                void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
     250                void           NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss);
    249251                IssmDouble     MinEdgeLength(IssmDouble xyz_list[6][3]);
    250252                void             SetClone(int* minranks);

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

    r16803 r16805  
    111111                void        JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
    112112                void        NodalFunctions(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
     113                void        NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
     114                void        NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
    113115                bool        NoIceInElement(){_error_("not implemented yet");};
    114116                int         NumberofNodesVelocity(void){_error_("not implemented yet");};

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

    r16801 r16805  
    22072207        _assert_(gauss->Enum()==GaussTriaEnum);
    22082208        this->GetNodalFunctions(basis,(GaussTria*)gauss);
     2209
     2210}
     2211/*}}}*/
     2212/*FUNCTION Tria::NodalFunctionsVelocity{{{*/
     2213void Tria::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){
     2214
     2215        _assert_(gauss->Enum()==GaussTriaEnum);
     2216        this->GetNodalFunctionsVelocity(basis,(GaussTria*)gauss);
     2217
     2218}
     2219/*}}}*/
     2220/*FUNCTION Tria::NodalFunctionsPressure{{{*/
     2221void Tria::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){
     2222
     2223        _assert_(gauss->Enum()==GaussTriaEnum);
     2224        this->GetNodalFunctionsPressure(basis,(GaussTria*)gauss);
    22092225
    22102226}

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

    r16803 r16805  
    274274                ElementVector* NewElementVector(int approximation_enum);
    275275                void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
     276                void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
     277                void           NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss);
    276278                void             SetClone(int* minranks);
    277279                Seg*             SpawnSeg(int index1,int index2);

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

    r16789 r16805  
    4141                for (i=0;i<femmodel->elements->Size();i++){
    4242                        element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
    43                         element->CreateKMatrix(Kff_temp,NULL);
     43                        ElementMatrix* Ke = element->CreateKMatrix();
     44                        ElementVector* pe = element->CreatePVector();
     45                        element->ReduceMatrices(Ke,pe);
     46                        if(Ke) Ke->AddToGlobal(Kff_temp,NULL);
     47                        delete Ke;
     48                        delete pe;
    4449                }
    4550
    … …  
    6873                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
    6974                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
    70                 //ElementVector* pe = analysis->CreatePVector(element);
     75                ElementVector* pe = analysis->CreatePVector(element);
    7176                //ElementMatrix* Ke = analysis->CreateKMatrix(element);
    72                 ElementVector* pe = element->CreatePVector();
     77                //ElementVector* pe = element->CreatePVector();
    7378                ElementMatrix* Ke = element->CreateKMatrix();
    7479                element->ReduceMatrices(Ke,pe);
Note: See TracChangeset for help on using the changeset viewer.