Changeset 16693


Ignore:
Timestamp:
11/08/13 16:32:04 (11 years ago)
Author:
seroussi
Message:

NEW: working on InputUpdateFromSolutionHoriz

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

Legend:

Unmodified
Added
Removed

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

    r16684 r16693  
    930930}/*}}}*/
    931931void StressbalanceAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    932         _error_("not implemented yet");
     932
     933        int approximation;
     934        element->GetInputValue(&approximation,ApproximationEnum);
     935        switch(approximation){
     936                case FSApproximationEnum: case NoneApproximationEnum:
     937                        //InputUpdateFromSolutionFS(solution,element);
     938                        return;
     939                case SSAApproximationEnum: case HOApproximationEnum: case SIAApproximationEnum:
     940                        InputUpdateFromSolutionHoriz(solution,element);
     941                        return;
     942                case L1L2ApproximationEnum:
     943                        //InputUpdateFromSolutionHoriz(solution,element);
     944                        return;
     945                case SSAHOApproximationEnum: case HOFSApproximationEnum: case SSAFSApproximationEnum:
     946                        /*the elements around will create the solution*/
     947                        return;
     948                default:
     949                        _error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
     950        }
    933951}/*}}}*/
     952void StressbalanceAnalysis::InputUpdateFromSolutionHoriz(IssmDouble* solution,Element* element){/*{{{*/
     953
     954        int        i;
     955        IssmDouble rho_ice,g;
     956        int*       doflist=NULL;
     957
     958        /*Fetch number of nodes and dof for this finite element*/
     959        int numnodes = element->GetNumberOfNodes();
     960        int numdof   = numnodes*2;
     961
     962        /*Fetch dof list and allocate solution vectors*/
     963        element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     964        IssmDouble* values    = xNew<IssmDouble>(numdof);
     965        IssmDouble* vx        = xNew<IssmDouble>(numnodes);
     966        IssmDouble* vy        = xNew<IssmDouble>(numnodes);
     967        IssmDouble* vz        = xNew<IssmDouble>(numnodes);
     968        IssmDouble* vel       = xNew<IssmDouble>(numnodes);
     969        IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
     970        IssmDouble* thickness = xNew<IssmDouble>(numnodes);
     971
     972        /*Use the dof list to index into the solution vector: */
     973        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
     974
     975        /*Transform solution in Cartesian Space*/
     976        element->TransformSolutionCoord(&values[0],XYEnum);
     977
     978        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
     979        for(i=0;i<numnodes;i++){
     980                vx[i]=values[i*NDOF2+0];
     981                vy[i]=values[i*NDOF2+1];
     982
     983                /*Check solution*/
     984                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
     985                if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
     986        }
     987
     988        /*Get Vz and compute vel*/
     989        element->GetInputListOnNodes(&vz[0],VzEnum,0.);
     990        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
     991
     992        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
     993         *so the pressure is just the pressure at the bedrock: */
     994        Matpar* matpar=element->GetMatparPointer();
     995        rho_ice=matpar->GetRhoIce();
     996        g=matpar->GetG();
     997        element->GetInputListOnNodes(&thickness[0],ThicknessEnum);
     998        for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*thickness[i];
     999
     1000        /*Now, we have to move the previous Vx and Vy inputs  to old
     1001         * status, otherwise, we'll wipe them off: */
     1002        element->InputChangeName(VxEnum,VxPicardEnum);
     1003        element->InputChangeName(VyEnum,VyPicardEnum);
     1004        element->InputChangeName(PressureEnum,PressurePicardEnum);
     1005
     1006        /*Add vx and vy as inputs to the tria element: */
     1007        element->AddInput(VxEnum,vx,P1Enum);
     1008        element->AddInput(VyEnum,vy,P1Enum);
     1009        element->AddInput(VelEnum,vel,P1Enum);
     1010        element->AddInput(PressureEnum,pressure,P1Enum);
     1011
     1012        /*Free ressources:*/
     1013        xDelete<IssmDouble>(thickness);
     1014        xDelete<IssmDouble>(pressure);
     1015        xDelete<IssmDouble>(vel);
     1016        xDelete<IssmDouble>(vz);
     1017        xDelete<IssmDouble>(vy);
     1018        xDelete<IssmDouble>(vx);
     1019        xDelete<IssmDouble>(values);
     1020        xDelete<int>(doflist);
     1021
     1022}/*}}}*/

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

    r16684 r16693  
    1919                void CreateLoads(Loads* loads, IoModel* iomodel);
    2020                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    21                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    2221                void GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
    2322                void GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element);
     23                void InputUpdateFromSolutionHoriz(IssmDouble* solution,Element* element);
     24                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    2425};
    2526#endif

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

    r16685 r16693  
    2121class Vertices;
    2222class Materials;
     23class Matpar;
    2324class Input;
    2425class Gauss;
    … …  
    3839                virtual void   SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0;
    3940                virtual void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;
     41                virtual void   AddInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
    4042                virtual void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>*  Kfs)=0;
    4143                virtual void   CreateDVector(Vector<IssmDouble>* df)=0;
    … …  
    4547                virtual void   FindParam(IssmDouble* pvalue,int paramenum)=0;
    4648                virtual int    FiniteElement(void)=0;
     49                virtual Matpar* GetMatparPointer(void)=0;
     50                virtual void   TransformSolutionCoord(IssmDouble* values,int transformenum)=0;
    4751                virtual void    GetDofList(int** pdoflist,int approximation_enum,int setenum)=0;
    4852                virtual void    GetDofListVelocity(int** pdoflist,int setenum)=0;
    … …  
    5862                virtual bool   IsNodeOnShelfFromFlags(IssmDouble* flags)=0;
    5963                virtual bool   IsOnBed()=0;
     64                virtual void   GetInputListOnNodes(IssmDouble* pvalue,int enumtype)=0;
     65                virtual void   GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue)=0;
    6066                virtual void   GetInputListOnVertices(IssmDouble* pvalue,int enumtype)=0;
    6167                virtual void   GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue)=0;
    … …  
    6975                virtual IssmDouble SurfaceArea(void)=0;
    7076                virtual void   InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum)=0;
     77                virtual void   InputChangeName(int enum_type,int enum_type_old)=0;
    7178                virtual void   ComputeBasalStress(Vector<IssmDouble>* sigma_b)=0;
    7279                virtual void   ComputeStrainRate(Vector<IssmDouble>* eps)=0;

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

    r16690 r16693  
    18151815}
    18161816/*}}}*/
     1817/*FUNCTION Penta::InputChangeName{{{*/
     1818void  Penta::InputChangeName(int new_enum,int original_enum){
     1819
     1820        /*Call inputs method*/
     1821        this->inputs->ChangeEnum(original_enum,new_enum);
     1822}
     1823/*}}}*/
    18171824/*FUNCTION Penta::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{*/
    18181825void Penta::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){
    … …  
    32273234        return dt;
    32283235}/*}}}*/
     3236/*FUNCTION Penta::TransformSolutionCoord{{{*/
     3237void Penta::TransformSolutionCoord(IssmDouble* values,int transformenum){
     3238
     3239        ::TransformSolutionCoord(values,this->nodes,this->NumberofNodes(),transformenum);
     3240
     3241}
     3242/*}}}*/
    32293243/*FUNCTION Penta::Update {{{*/
    32303244void Penta::Update(int index,IoModel* iomodel,int analysis_counter,int analysis_type,int finiteelement_type){
    … …  
    62866300
    62876301        /*Transform solution in Cartesian Space*/
    6288         TransformSolutionCoord(&vvalues[0],nodes,vnumnodes,XYZEnum);
     6302        ::TransformSolutionCoord(&vvalues[0],nodes,vnumnodes,XYZEnum);
    62896303
    62906304        /*fill in all arrays: */
    … …  
    97999813
    98009814        /*Transform solution in Cartesian Space*/
    9801         TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,XYEnum); /*2D: only the first 3 nodes are taken*/
     9815        ::TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,XYEnum); /*2D: only the first 3 nodes are taken*/
    98029816
    98039817        /*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */
    … …  
    98989912
    98999913        /*Transform solution in Cartesian Space*/
    9900         TransformSolutionCoord(&SSA_values[0],penta->nodes,NUMVERTICES,XYEnum);
    9901         TransformSolutionCoord(&HO_values[0],   this->nodes,NUMVERTICES,XYEnum);
     9914        ::TransformSolutionCoord(&SSA_values[0],penta->nodes,NUMVERTICES,XYEnum);
     9915        ::TransformSolutionCoord(&HO_values[0],   this->nodes,NUMVERTICES,XYEnum);
    99029916
    99039917        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    … …  
    999110005
    999210006        /*Transform solution in Cartesian Space*/
    9993         TransformSolutionCoord(&SSA_values[0],this->nodes,NUMVERTICES,XYEnum);
    9994         TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
     10007        ::TransformSolutionCoord(&SSA_values[0],this->nodes,NUMVERTICES,XYEnum);
     10008        ::TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
    999510009
    999610010        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    … …  
    1007310087
    1007410088        /*Transform solution in Cartesian Space*/
    10075         TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,XYEnum); /*2D: only the first 3 nodes are taken*/
     10089        ::TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,XYEnum); /*2D: only the first 3 nodes are taken*/
    1007610090
    1007710091        /*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */
    … …  
    1015510169
    1015610170        /*Transform solution in Cartesian Space*/
    10157         TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum);
     10171        ::TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum);
    1015810172        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    1015910173
    … …  
    1024510259
    1024610260        /*Transform solution in Cartesian Space*/
    10247         TransformSolutionCoord(&HO_values[0],this->nodes,NUMVERTICES,XYEnum);
    10248         TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
     10261        ::TransformSolutionCoord(&HO_values[0],this->nodes,NUMVERTICES,XYEnum);
     10262        ::TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
    1024910263
    1025010264        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    … …  
    1050610520
    1050710521        /*Transform solution in Cartesian Space*/
    10508         TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
     10522        ::TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
    1050910523
    1051010524        /*Ok, we have vx and vy in values, fill in all arrays: */

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

    r16685 r16693  
    9191                void   GetNodesLidList(int* lidlist);
    9292                int    GetNumberOfNodes(void);
     93                Matpar* GetMatparPointer(void){_error_("not implemented yet");};
    9394                void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
    9495                IssmDouble GetZcoord(GaussPenta* gauss);
    … …  
    179180                /*}}}*/
    180181                /*Penta specific routines:{{{*/
     182                void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
    181183                void             BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]);
    182184                ElementMatrix* CreateBasalMassMatrix(void);
    … …  
    218220                Penta*         GetLowerElement(void);
    219221                Penta*         GetBasalElement(void);
     222                void           InputChangeName(int input_enum, int enum_type_old);
    220223                void             InputExtrude(int enum_type,int object_type);
    221224                void           InputUpdateFromSolutionMasstransport(IssmDouble* solutiong);
    … …  
    235238                Tria*            SpawnTria(int location);
    236239                void             SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);
     240                void           TransformSolutionCoord(IssmDouble* values,int transformenum);
    237241
    238242                #ifdef _HAVE_STRESSBALANCE_

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

    r16685 r16693  
    6666                /*}}}*/
    6767                /*Element virtual functions definitions: {{{*/
     68                void        AddInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
    6869                void        ComputeBasalStress(Vector<IssmDouble>* sigma_b){_error_("not implemented yet");};
    6970                void        ComputeStrainRate(Vector<IssmDouble>* eps){_error_("not implemented yet");};
    … …  
    8485                void        GetDofListVelocity(int** pdoflist,int setenum){_error_("not implemented yet");};
    8586                void        GetDofListPressure(int** pdoflist,int setenum){_error_("not implemented yet");};
     87                void        GetInputListOnNodes(IssmDouble* pvalue,int enumtype){_error_("not implemented yet");};
     88                void        GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){_error_("not implemented yet");};
     89                Matpar*     GetMatparPointer(void){_error_("not implemented yet");};
    8690                int         GetNodeIndex(Node* node){_error_("not implemented yet");};
    8791                void        GetNodesSidList(int* sidlist){_error_("not implemented yet");};
    … …  
    8993                int         GetNumberOfNodes(void){_error_("not implemented yet");};
    9094                int         Sid(){_error_("not implemented yet");};
     95                void        InputChangeName(int input_enum, int enum_type_old){_error_("not implemented yet");};
    9196                bool        IsOnBed(){_error_("not implemented yet");};
    9297                bool        IsFloating(){_error_("not implemented yet");};
    … …  
    117122                void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){_error_("not implemented yet");};
    118123                void    ComputeEPLThickness(void){_error_("not implemented yet");};
    119                 void    UpdateConstraintsL2ProjectionEPL(){_error_("not implemented");};
    120124                #endif
    121125                void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type){_error_("not implemented yet");};
    … …  
    135139                void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement){_error_("not implemented yet");};
    136140                IssmDouble  TimeAdapt(){_error_("not implemented yet");};
     141                void        TransformSolutionCoord(IssmDouble* values,int transformenum){_error_("not implemented yet");};
    137142                void UpdateConstraintsExtrudeFromBase(){_error_("not implemented");};
    138143                void UpdateConstraintsExtrudeFromTop(){_error_("not implemented");};
     144                void UpdateConstraintsL2ProjectionEPL(){_error_("not implemented");};
    139145
    140146#ifdef _HAVE_RESPONSES_

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

    r16691 r16693  
    164164        *pd_nz=d_nz;
    165165        *po_nz=o_nz;
     166}
     167/*}}}*/
     168/*FUNCTION Tria::AddInput{{{*/
     169void  Tria::AddInput(int input_enum,IssmDouble* values, int interpolation_enum){
     170
     171        /*Call inputs method*/
     172        this->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
    166173}
    167174/*}}}*/
    … …  
    14901497}
    14911498/*}}}*/
     1499/*FUNCTION Tria::GetMatparPointer(void) {{{*/
     1500Matpar* Tria::GetMatparPointer(void){
     1501        return this->matpar;
     1502}
     1503/*}}}*/
    14921504/*FUNCTION Tria::GetVertexPidList {{{*/
    14931505void  Tria::GetVertexPidList(int* doflist){
    … …  
    15981610                inputs->DuplicateInput(original_enum,new_enum);
    15991611        }
     1612}
     1613/*}}}*/
     1614/*FUNCTION Tria::InputChangeName{{{*/
     1615void  Tria::InputChangeName(int new_enum,int original_enum){
     1616
     1617        /*Call inputs method*/
     1618        this->inputs->ChangeEnum(original_enum,new_enum);
    16001619}
    16011620/*}}}*/
    … …  
    27652784}
    27662785/*}}}*/
     2786/*FUNCTION Tria::TransformSolutionCoord{{{*/
     2787void Tria::TransformSolutionCoord(IssmDouble* values,int transformenum){
     2788
     2789        ::TransformSolutionCoord(values,this->nodes,this->NumberofNodes(),transformenum);
     2790
     2791}
     2792/*}}}*/
    27672793/*FUNCTION Tria::Update{{{*/
    27682794void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finiteelement_type){
    … …  
    43394365
    43404366        /*Transform solution in Cartesian Space*/
    4341         TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
     4367        ::TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
    43424368
    43434369        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    … …  
    44214447
    44224448        /*Transform solution in Cartesian Space*/
    4423         TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
     4449        ::TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
    44244450
    44254451        /*Ok, we have vx and vy in values, fill in all arrays: */
    … …  
    44884514
    44894515        /*Transform solution in Cartesian Space*/
    4490         TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
     4516        ::TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
    44914517
    44924518        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    … …  
    62586284
    62596285        /*Transform solution in Cartesian Space*/
    6260         TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
     6286        ::TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
    62616287
    62626288        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */

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

    r16685 r16693  
    191191                /*}}}*/
    192192                /*Tria specific routines:{{{*/
     193                void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
    193194                ElementMatrix* CreateKMatrix(void);
    194195                ElementMatrix* CreateKMatrixBalancethickness(void);
    … …  
    237238                IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
    238239                void           GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[2][3]);
     240      Matpar*        GetMatparPointer(void);
    239241                void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum);
    240242                Input*         GetInput(int inputenum);
    … …  
    250252                void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    251253                void           GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input);
     254                void             InputChangeName(int enum_type,int enum_type_old);
    252255                void             InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
    253256                void             InputUpdateFromSolutionMasstransport(IssmDouble* solution);
    … …  
    257260                Seg*             SpawnSeg(int index1,int index2);
    258261                void             SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);
     262                void           TransformSolutionCoord(IssmDouble* values,int transformenum);
    259263                #ifdef _HAVE_STRESSBALANCE_
    260264                ElementMatrix* CreateKMatrixStressbalanceSSA(void);
Note: See TracChangeset for help on using the changeset viewer.