Changeset 16382


Ignore:
Timestamp:
10/11/13 10:13:16 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: added bed slope calculation (needs debugging still)

Location:
issm/trunk-jpl/src/c
Files:
2 added
35 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r16378 r16382  
    9292                                        ./classes/Inputs/Input.h\
    9393                                        ./classes/Inputs/InputLocal.h\
     94                                        ./classes/Inputs/SegInput.h\
     95                                        ./classes/Inputs/SegInput.cpp\
    9496                                        ./classes/Inputs/TriaInput.h\
    9597                                        ./classes/Inputs/TriaInput.cpp\

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

    r15849 r16382  
    1414        /*parameters: */
    1515        bool save_results;
     16        int  meshtype;
    1617
    1718        /*Recover some parameters: */
    1819        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
     20        femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
    1921
    2022        if(VerboseSolution()) _printf0_("   computing slope\n");
    … …  
    2325        femmodel->SetCurrentConfiguration(BedSlopeAnalysisEnum,BedSlopeXAnalysisEnum);
    2426        solutionsequence_linear(femmodel);
    25         femmodel->SetCurrentConfiguration(BedSlopeAnalysisEnum,BedSlopeYAnalysisEnum);
    26         solutionsequence_linear(femmodel);
     27        if(meshtype!=Mesh2DverticalEnum){
     28                femmodel->SetCurrentConfiguration(BedSlopeAnalysisEnum,BedSlopeYAnalysisEnum);
     29                solutionsequence_linear(femmodel);
     30        }
    2731
    2832        if(save_results){
    2933                if(VerboseSolution()) _printf0_("   saving results\n");
    3034                InputToResultx(femmodel,BedSlopeXEnum);
    31                 InputToResultx(femmodel,BedSlopeYEnum);
     35                if(meshtype!=Mesh2DverticalEnum) InputToResultx(femmodel,BedSlopeYEnum);
    3236        }
    3337

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

    r16363 r16382  
    4949        /*Compute slopes: */
    5050        if(isSIA) surfaceslope_core(femmodel);
    51         if(isFS && meshtype==Mesh3DEnum){
     51        if(isFS){
    5252                bedslope_core(femmodel);
    5353                femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);

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

    r15737 r16382  
    128128}
    129129/*}}}*/
     130/*FUNCTION ElementHook::SpawnSegHook{{{*/
     131void ElementHook::SpawnSegHook(ElementHook* triahook,int index1,int index2){
     132
     133        triahook->numanalyses=this->numanalyses;
     134
     135        int indices[2];
     136        indices[0]=index1;
     137        indices[1]=index2;
     138
     139        /*Spawn nodes hook*/
     140        triahook->hnodes=new Hook*[this->numanalyses];
     141        for(int i=0;i<this->numanalyses;i++){
     142                /*Do not do anything if Hook is empty*/
     143                if (!this->hnodes[i] || this->hnodes[i]->GetNum()==0){
     144                        triahook->hnodes[i]=NULL;
     145                }
     146                else{
     147                        triahook->hnodes[i]=this->hnodes[i]->Spawn(indices,2);
     148                }
     149        }
     150
     151        /*do not spawn hmaterial. material will be taken care of by Tria*/
     152        triahook->hmaterial=NULL;
     153        triahook->hvertices=(Hook*)this->hvertices->Spawn(indices,2);
     154        triahook->hmatpar=(Hook*)this->hmatpar->copy();
     155}
     156/*}}}*/

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

    r15737 r16382  
    2626                void SetHookNodes(int* node_ids,int numnodes,int analysis_counter);
    2727                void SpawnTriaHook(ElementHook* triahook,int location); //3d only
     28                void SpawnSegHook(ElementHook* triahook,int ndex1,int index2); //2d only
    2829                void InitHookNeighbors(int* element_ids);               //3d only
    2930};

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

    r16343 r16382  
    139139}
    140140/*}}}*/
     141
     142/*FUNCTION Seg::CreateMassMatrix {{{*/
     143ElementMatrix* Seg::CreateMassMatrix(void){
     144
     145        /* Intermediaries */
     146        IssmDouble  D,Jdet;
     147        IssmDouble  xyz_list[NUMVERTICES][2];
     148
     149        /*Fetch number of nodes and dof for this finite element*/
     150        int numnodes = this->NumberofNodes();
     151
     152        /*Initialize Element matrix and vectors*/
     153        ElementMatrix* Ke    = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     154        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     155
     156        /*Retrieve all inputs and parameters*/
     157        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     158
     159        /* Start looping on the number of gaussian points: */
     160        GaussSeg* gauss=new GaussSeg(2);
     161        for(int ig=gauss->begin();ig<gauss->end();ig++){
     162
     163                gauss->GaussPoint(ig);
     164
     165                GetNodalFunctions(basis,gauss);
     166                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     167                D=gauss->weight*Jdet;
     168
     169                TripleMultiply(basis,1,numnodes,1,
     170                                        &D,1,1,0,
     171                                        basis,1,numnodes,0,
     172                                        &Ke->values[0],1);
     173        }
     174
     175        /*Clean up and return*/
     176        delete gauss;
     177        xDelete<IssmDouble>(basis);
     178        return Ke;
     179}
     180/*}}}*/
     181/*FUNCTION Seg::CreatePVectorSlope {{{*/
     182ElementVector* Seg::CreatePVectorSlope(void){
     183
     184        /*Intermediaries */
     185        int        i,analysis_type;
     186        IssmDouble Jdet;
     187        IssmDouble xyz_list[NUMVERTICES][2];
     188        IssmDouble slope;
     189
     190        /*Fetch number of nodes and dof for this finite element*/
     191        int numnodes = this->NumberofNodes();
     192
     193        /*Initialize Element vector*/
     194        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     195        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     196
     197        /*Retrieve all inputs and parameters*/
     198        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     199        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     200        Input* slope_input=NULL;
     201        if(analysis_type==SurfaceSlopeXAnalysisEnum){
     202                slope_input=inputs->GetInput(SurfaceEnum); _assert_(slope_input);
     203        }
     204        else if(analysis_type==BedSlopeXAnalysisEnum){
     205                slope_input=inputs->GetInput(BedEnum);     _assert_(slope_input);
     206        }
     207        else{
     208                _error_("Analysis "<<EnumToStringx(analysis_type)<<" not inplemented yet");
     209        }
     210
     211        /* Start  looping on the number of gaussian points: */
     212        GaussSeg* gauss=new GaussSeg(2);
     213        for(int ig=gauss->begin();ig<gauss->end();ig++){
     214
     215                gauss->GaussPoint(ig);
     216
     217                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     218                GetNodalFunctions(basis,gauss);
     219
     220                slope_input->GetInputDerivativeValue(&slope,&xyz_list[0][0],gauss);
     221
     222                for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*slope*basis[i];
     223        }
     224
     225        /*Clean up and return*/
     226        xDelete<IssmDouble>(basis);
     227        delete gauss;
     228        return pe;
     229}
     230/*}}}*/

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

    r16343 r16382  
    7676                void        CreateDVector(Vector<IssmDouble>* df){_error_("not implemented yet");};
    7777                void        CreatePVector(Vector<IssmDouble>* pf){_error_("not implemented yet");};
     78                ElementVector* CreatePVectorSlope(void);
    7879                void        CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("not implemented yet");};
    7980                void        Delta18oParameterization(void){_error_("not implemented yet");};
    … …  
    8788                bool        IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");};
    8889                bool        NoIceInElement(){_error_("not implemented yet");};
     90                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype){_error_("not implemented yet");};
     91                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){_error_("not implemented yet");};
     92                void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
     93                void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
     94                #ifdef _HAVE_THERMAL_
     95                void UpdateBasalConstraintsEnthalpy(void){_error_("not implemented yet");};
     96                void ComputeBasalMeltingrate(void){_error_("not implemented yet");};
     97                void DrainWaterfraction(void){_error_("not implemented yet");};
     98                #endif
     99                #ifdef _HAVE_HYDROLOGY_
     100                void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){_error_("not implemented yet");};
     101                void    GetHydrologyTransfer(Vector<IssmDouble>* transfer){_error_("not implemented yet");};
     102                void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){_error_("not implemented yet");};
     103                void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){_error_("not implemented yet");};
     104                #endif
    89105                void        GetSolutionFromInputs(Vector<IssmDouble>* solution){_error_("not implemented yet");};
    90106                void        GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum){_error_("not implemented yet");};
    … …  
    173189#endif
    174190                /*}}}*/
     191                /*Seg specific routines:*/
     192                ElementMatrix* CreateMassMatrix(void);
    175193};
    176194#endif  /* _SEG_H */

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

    r16378 r16382  
    4848
    4949/*Reference Element numerics*/
    50 /*FUNCTION SegRef::GetNodalFunctions{{{*/
     50/*FUNCTION SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss){{{*/
    5151void SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss){
    5252        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    … …  
    5454        _assert_(basis);
    5555
    56         switch(this->element_type){
     56        GetNodalFunctions(basis,gauss,this->element_type);
     57}
     58/*}}}*/
     59/*FUNCTION SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss,int finiteelement){{{*/
     60void SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss,int finiteelement){
     61        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     62
     63        _assert_(basis);
     64
     65        switch(element_type){
    5766                case P1Enum: case P1DGEnum:
    5867                        basis[0]=(1.-gauss->coord1)/2.;
    … …  
    6069                        return;
    6170                default:
     71                        _error_("Element type "<<EnumToStringx(element_type)<<" not supported yet");
     72        }
     73}
     74/*}}}*/
     75/*FUNCTION SegRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss){{{*/
     76void SegRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss){
     77
     78        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,this->element_type);
     79
     80}
     81/*}}}*/
     82/*FUNCTION SegRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement){{{*/
     83void SegRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement){
     84
     85        /*This routine returns the values of the nodal functions derivatives  (with respect to the
     86         * actual coordinate system): */
     87        IssmDouble    Jinv;
     88
     89        /*Fetch number of nodes for this finite element*/
     90        int numnodes = this->NumberofNodes(finiteelement);
     91
     92        /*Get nodal functions derivatives in reference triangle*/
     93        IssmDouble* dbasis_ref=xNew<IssmDouble>(numnodes);
     94        GetNodalFunctionsDerivativesReference(dbasis_ref,gauss,finiteelement);
     95
     96        /*Get Jacobian invert: */
     97        GetJacobianInvert(&Jinv, xyz_list, gauss);
     98
     99        /*Build dbasis:
     100         * [dhi/dx]= Jinv*[dhi/dr]
     101         */
     102        for(int i=0;i<numnodes;i++){
     103                dbasis[i] = Jinv*dbasis_ref[i];
     104        }
     105
     106        /*Clean up*/
     107        xDelete<IssmDouble>(dbasis_ref);
     108
     109}
     110/*}}}*/
     111/*FUNCTION SegRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss){{{*/
     112void SegRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss){
     113        /*This routine returns the values of the nodal functions derivatives  (with respect to the
     114         * natural coordinate system) at the gaussian point. */
     115
     116        GetNodalFunctionsDerivativesReference(dbasis,gauss,this->element_type);
     117
     118}
     119/*}}}*/
     120/*FUNCTION SegRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss,int finiteelement){{{*/
     121void SegRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss,int finiteelement){
     122        /*This routine returns the values of the nodal functions derivatives  (with respect to the
     123         * natural coordinate system) at the gaussian point. */
     124
     125        _assert_(dbasis && gauss);
     126
     127        switch(finiteelement){
     128                case P1Enum: case P1DGEnum:
     129                        /*Nodal function 1*/
     130                        dbasis[0] = -0.5;
     131                        /*Nodal function 2*/
     132                        dbasis[1] = 0.5;
     133                        return;
     134                default:
    62135                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
    63136        }
     137
     138}
     139/*}}}*/
     140/*FUNCTION SegRef::GetInputDerivativeValue{{{*/
     141void SegRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss){
     142
     143        /*From node values of parameter p (plist[0],plist[1]), return parameter derivative value at gaussian
     144         * point specified by gauss_basis:
     145         *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx
     146         *
     147         * p is a vector already allocated.
     148         */
     149
     150        /*Output*/
     151        IssmDouble dpx;
     152
     153        /*Fetch number of nodes for this finite element*/
     154        int numnodes = this->NumberofNodes();
     155
     156        /*Get nodal functions derivatives*/
     157        IssmDouble* dbasis=xNew<IssmDouble>(1*numnodes);
     158        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     159
     160        /*Calculate parameter for this Gauss point*/
     161        for(int i=0;i<numnodes;i++) dpx += dbasis[i]*plist[i];
     162
     163        /*Assign values*/
     164        xDelete<IssmDouble>(dbasis);
     165        *p=dpx;
     166
     167}
     168/*}}}*/
     169/*FUNCTION SegRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss){{{*/
     170void SegRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss){
     171
     172        GetInputValue(p,plist,gauss,this->element_type);
     173}
     174/*}}}*/
     175/*FUNCTION SegRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss,int finiteelement){{{*/
     176void SegRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss,int finiteelement){
     177
     178        /*Output*/
     179        IssmDouble value =0.;
     180
     181        /*Fetch number of nodes for this finite element*/
     182        int numnodes = this->NumberofNodes(finiteelement);
     183
     184        /*Get nodal functions*/
     185        IssmDouble* basis=xNew<IssmDouble>(numnodes);
     186        GetNodalFunctions(basis, gauss,finiteelement);
     187
     188        /*Calculate parameter for this Gauss point*/
     189        for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
     190
     191        /*Assign output pointer*/
     192        xDelete<IssmDouble>(basis);
     193        *p = value;
    64194}
    65195/*}}}*/
    … …  
    103233}
    104234/*}}}*/
     235/*FUNCTION SegRef::NumberofNodes(){{{*/
     236int SegRef::NumberofNodes(void){
     237
     238        return this->NumberofNodes(this->element_type);
     239}
     240/*}}}*/
     241/*FUNCTION SegRef::NumberofNodes(int finiteelement){{{*/
     242int SegRef::NumberofNodes(int finiteelement){
     243
     244        switch(finiteelement){
     245                case P1Enum:                return NUMNODESP1;
     246                case P1DGEnum:              return NUMNODESP1;
     247                default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     248        }
     249
     250        return -1;
     251}
     252/*}}}*/

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

    r16377 r16382  
    2626                void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussSeg* gauss);
    2727                void GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss);
     28                void GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss,int finiteelement);
     29                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss);
     30                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement);
     31                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss);
     32                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss,int finiteelement);
     33                void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss);
     34                void GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss);
     35                void GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss,int finiteelement);
     36
     37                int  NumberofNodes(void);
     38                int  NumberofNodes(int finiteelement);
    2839};
    2940#endif

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

    r16373 r16382  
    246246                #endif
    247247                case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
    248                         return CreateMassMatrix();
     248                        int meshtype;
     249                        parameters->FindParam(&meshtype,MeshTypeEnum);
     250                        if(meshtype==Mesh2DverticalEnum){
     251                                return CreateBasalMassMatrix();
     252                        }
     253                        else{
     254                                return CreateMassMatrix();
     255                        }
    249256                        break;
    250257                #ifdef _HAVE_MASSTRANSPORT_
    … …  
    330337}
    331338/*}}}*/
     339/*FUNCTION Tria::CreateBasalMassMatrix{{{*/
     340ElementMatrix* Tria::CreateBasalMassMatrix(void){
     341
     342        if(!HasEdgeOnBed()) return NULL;
     343
     344        int index1,index2;
     345        this->EdgeOnBedIndices(&index1,&index2);
     346
     347        Seg* seg=(Seg*)SpawnSeg(index1,index2);
     348        ElementMatrix* Ke=seg->CreateMassMatrix();
     349        delete seg->material; delete seg;
     350
     351        /*clean up and return*/
     352        return Ke;
     353}
     354/*}}}*/
    332355/*FUNCTION Tria::CreateDVector {{{*/
    333356void  Tria::CreateDVector(Vector<IssmDouble>* df){
    … …  
    413436                #endif
    414437                case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
    415                         return CreatePVectorSlope();
     438                        int meshtype;
     439                        parameters->FindParam(&meshtype,MeshTypeEnum);
     440                        if(meshtype==Mesh2DverticalEnum){
     441                                return CreateBasalPVectorSlope();
     442                        }
     443                        else{
     444                                return CreatePVectorSlope();
     445                        }
    416446                        break;
    417447                #ifdef _HAVE_MASSTRANSPORT_
    … …  
    459489        /*make compiler happy*/
    460490        return NULL;
     491}
     492/*}}}*/
     493/*FUNCTION Tria::CreateBasalPVectorSlope{{{*/
     494ElementVector* Tria::CreateBasalPVectorSlope(void){
     495
     496        if(!HasEdgeOnBed()) return NULL;
     497
     498        int index1,index2;
     499        this->EdgeOnBedIndices(&index1,&index2);
     500
     501        Seg* seg=(Seg*)SpawnSeg(index1,index2);
     502        ElementVector* pe=seg->CreatePVectorSlope();
     503        delete seg->material; delete seg;
     504
     505        /*clean up and return*/
     506        return pe;
    461507}
    462508/*}}}*/
    … …  
    24592505        else this->nodes=NULL;
    24602506
     2507}
     2508/*}}}*/
     2509/*FUNCTION Tria::SpawnSeg {{{*/
     2510Seg*  Tria::SpawnSeg(int index1,int index2){
     2511
     2512        int analysis_counter;
     2513
     2514        /*go into parameters and get the analysis_counter: */
     2515        this->parameters->FindParam(&analysis_counter,AnalysisCounterEnum);
     2516
     2517        /*Create Seg*/
     2518        Seg* seg=new Seg();
     2519        seg->id=this->id;
     2520        seg->inputs=(Inputs*)this->inputs->SpawnSegInputs(index1,index2);
     2521        seg->parameters=this->parameters;
     2522        seg->element_type=P1Enum; //Only P1 CG for now (TO BE CHANGED)
     2523        this->SpawnSegHook(dynamic_cast<ElementHook*>(seg),index1,index2);
     2524
     2525        /*Spawn material*/
     2526        seg->material=(Material*)this->material->copy();
     2527        delete seg->material->inputs;
     2528        seg->material->inputs=(Inputs*)this->material->inputs->SpawnSegInputs(index1,index2);
     2529
     2530        /*recover nodes, material and matpar: */
     2531        seg->nodes    = (Node**)seg->hnodes[analysis_counter]->deliverp();
     2532        seg->vertices = (Vertex**)seg->hvertices->deliverp();
     2533        seg->matpar   = (Matpar*)seg->hmatpar->delivers();
     2534
     2535        /*Return new Seg*/
     2536        return seg;
    24612537}
    24622538/*}}}*/

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

    r16364 r16382  
    1818class Material;
    1919class Matpar;
     20class Seg;
    2021class ElementMatrix;
    2122class ElementVector;
    … …  
    190191                ElementMatrix* CreateKMatrixFreeSurfaceBase(void);
    191192                ElementMatrix* CreateMassMatrix(void);
     193                ElementMatrix* CreateBasalMassMatrix(void);
    192194                ElementVector* CreatePVector(void);
    193195                ElementVector* CreatePVectorBalancethickness(void);
    … …  
    203205                ElementVector* CreatePVectorFreeSurfaceBase(void);
    204206                ElementVector* CreatePVectorSlope(void);
     207                ElementVector* CreateBasalPVectorSlope(void);
    205208                IssmDouble     GetArea(void);
    206209                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints);
    … …  
    228231                bool             IsInput(int name);
    229232                void             SetClone(int* minranks);
     233                Seg*             SpawnSeg(int index1,int index2);
    230234                void             SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);
    231235

  • issm/trunk-jpl/src/c/classes/Inputs/BoolInput.cpp

    r15737 r16382  
    7373/*FUNCTION BoolInput::SpawnTriaInput{{{*/
    7474Input* BoolInput::SpawnTriaInput(int location){
     75
     76                /*output*/
     77                BoolInput* outinput=new BoolInput();
     78
     79                /*only copy current value*/
     80                outinput->enum_type=this->enum_type;
     81                outinput->value=this->value;
     82
     83                /*Assign output*/
     84                return outinput;
     85
     86}
     87/*}}}*/
     88/*FUNCTION BoolInput::SpawnSegInput{{{*/
     89Input* BoolInput::SpawnSegInput(int index1,int index2){
    7590
    7691                /*output*/

  • issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h

    r16238 r16382  
    3535                int   InstanceEnum();
    3636                Input* SpawnTriaInput(int location);
     37                Input* SpawnSegInput(int index1,int index2);
    3738                Input* PointwiseDivide(Input* inputB){_error_("not implemented yet");};
    3839                Input* PointwiseMin(Input* inputB){_error_("not implemented yet");};
    … …  
    5253                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");};
    5354                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error_("not implemented yet");};
     55                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
    5456                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss);
    5557                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss);

  • issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp

    r16169 r16382  
    186186        return values->SpawnTriaInput(location);
    187187}/*}}}*/
     188/*FUNCTION ControlInput::SpawnSegInput{{{*/
     189Input* ControlInput::SpawnSegInput(int index1,int index2){
     190        return values->SpawnSegInput(index1,index2);
     191}/*}}}*/
    188192/*FUNCTION ControlInput::SpawnGradient{{{*/
    189193ElementResult* ControlInput::SpawnGradient(int step, IssmDouble time){

  • issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h

    r16238 r16382  
    3939                int    InstanceEnum();
    4040                Input* SpawnTriaInput(int location);
     41                Input* SpawnSegInput(int index1,int index2);
    4142                Input* PointwiseDivide(Input* inputB){_error_("not implemented yet");};
    4243                Input* PointwiseMin(Input* inputB){_error_("not implemented yet");};
    … …  
    5758                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");};
    5859                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error_("not implemented yet");};
     60                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
    5961                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss);
    6062                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss);

  • issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.cpp

    r16315 r16382  
    9696}
    9797/*}}}*/
     98/*FUNCTION DatasetInput::SpawnSegInput{{{*/
     99Input* DatasetInput::SpawnSegInput(int index1,int index2){
     100
     101        /*output*/
     102        DatasetInput* outinput=NULL;
     103
     104        /*Create new Datasetinput (copy of current input)*/
     105        outinput=new DatasetInput();
     106        outinput->enum_type=this->enum_type;
     107        outinput->inputs=dynamic_cast<Inputs*>(this->inputs->SpawnSegInputs(index1,index2));
     108        outinput->numids=this->numids;
     109        outinput->ids=xNew<int>(this->numids);
     110        xMemCpy(outinput->ids,this->ids,this->numids);
     111
     112        /*Assign output*/
     113        return outinput;
     114}
     115/*}}}*/
    98116
    99117/*DatasetInput management*/

  • issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h

    r16313 r16382  
    3737                int    InstanceEnum();
    3838                Input* SpawnTriaInput(int location);
     39                Input* SpawnSegInput(int index1,int index2);
    3940                Input* PointwiseDivide(Input* inputB){_error_("not implemented yet");};
    4041                Input* PointwiseMin(Input* inputB){_error_("not implemented yet");};
    … …  
    5455                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index);
    5556                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index);
     57                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
    5658                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");};
    5759                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.cpp

    r15737 r16382  
    8686}
    8787/*}}}*/
     88/*FUNCTION DoubleInput::SpawnSegInput{{{*/
     89Input* DoubleInput::SpawnSegInput(int index1,int index2){
     90
     91        /*output*/
     92        DoubleInput* outinput=new DoubleInput();
     93
     94        /*only copy current value*/
     95        outinput->enum_type=this->enum_type;
     96        outinput->value=this->value;
     97
     98        /*Assign output*/
     99        return outinput;
     100
     101}
     102/*}}}*/
    88103/*FUNCTION DoubleInput::SpawnResult{{{*/
    89104ElementResult* DoubleInput::SpawnResult(int step, IssmDouble time){

  • issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h

    r16238 r16382  
    3434                int   InstanceEnum();
    3535                Input* SpawnTriaInput(int location);
     36                Input* SpawnSegInput(int index1,int index2);
    3637                Input* PointwiseDivide(Input* inputB);
    3738                Input* PointwiseMin(Input* inputB);
    … …  
    5152                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");};
    5253                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error_("not implemented yet");};
     54                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
    5355                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss);
    5456                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss);

  • issm/trunk-jpl/src/c/classes/Inputs/Input.h

    r16238 r16382  
    1313class ElementResult;
    1414class GaussTria;
     15class GaussSeg;
    1516class Parameters;
    1617class GaussPenta;
    … …  
    3435                virtual void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index)=0;
    3536                virtual void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int index)=0;
     37                virtual void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss)=0;
    3638                virtual void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss)=0;
    3739                virtual void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss)=0;
    … …  
    6567
    6668                virtual Input* SpawnTriaInput(int location)=0;
     69                virtual Input* SpawnSegInput(int index1,int index2)=0;
    6770                virtual Input* PointwiseDivide(Input* inputB)=0;
    6871                virtual Input* PointwiseMax(Input* inputmax)=0;

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

    r16149 r16382  
    395395}
    396396/*}}}*/
     397/*FUNCTION Inputs::SpawnSegInputs{{{*/
     398Inputs* Inputs::SpawnSegInputs(int index1,int index2){
     399
     400        /*Intermediary*/
     401        vector<Object*>::iterator object;
     402        Input* inputin=NULL;
     403        Input* inputout=NULL;
     404
     405        /*Output*/
     406        Inputs* newinputs=new Inputs();
     407
     408        /*Go through inputs and call Spawn function*/
     409        for ( object=objects.begin() ; object < objects.end(); object++ ){
     410
     411                /*Create new input*/
     412                inputin=dynamic_cast<Input*>(*object);
     413                inputout=inputin->SpawnSegInput(index1,index2);
     414
     415                /*Add input to new inputs*/
     416                newinputs->AddObject(inputout);
     417        }
     418
     419        /*Assign output pointer*/
     420        return newinputs;
     421}
     422/*}}}*/
    397423/*FUNCTION Inputs::AXPY{{{*/
    398424void  Inputs::AXPY(int inputy_enum, IssmDouble scalar, int inputx_enum){

  • issm/trunk-jpl/src/c/classes/Inputs/Inputs.h

    r16149 r16382  
    2929                Input*      GetInput(int enum_name);
    3030                Inputs*     SpawnTriaInputs(int position);
     31                Inputs*     SpawnSegInputs(int index1,int index2);
     32                Inputs*     SpawnSegInputs(int position);
    3133                void        AXPY(int inputy_enum, IssmDouble scalar, int inputx_enum);
    3234                IssmDouble  InfinityNorm(int enumtype);

  • issm/trunk-jpl/src/c/classes/Inputs/IntInput.cpp

    r15737 r16382  
    7373/*FUNCTION IntInput::SpawnTriaInput{{{*/
    7474Input* IntInput::SpawnTriaInput(int location){
     75
     76        /*output*/
     77        IntInput* outinput=new IntInput();
     78
     79        /*only copy current value*/
     80        outinput->enum_type=this->enum_type;
     81        outinput->value=this->value;
     82
     83        /*Assign output*/
     84        return outinput;
     85}
     86/*}}}*/
     87/*FUNCTION IntInput::SpawnSegInput{{{*/
     88Input* IntInput::SpawnSegInput(int index1,int index2){
    7589
    7690        /*output*/

  • issm/trunk-jpl/src/c/classes/Inputs/IntInput.h

    r16238 r16382  
    3535                int   InstanceEnum();
    3636                Input* SpawnTriaInput(int location);
     37                Input* SpawnSegInput(int index1,int index2);
    3738                Input* PointwiseDivide(Input* inputB){_error_("not implemented yet");};
    3839                Input* PointwiseMin(Input* inputB){_error_("not implemented yet");};
    … …  
    5253                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");};
    5354                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error_("not implemented yet");};
     55                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
    5456                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss);
    5557                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss);

  • issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp

    r15737 r16382  
    124124}
    125125/*}}}*/
     126/*FUNCTION PentaInput::SpawnSegInput{{{*/
     127Input* PentaInput::SpawnSegInput(int index1,int index2){
     128
     129        _error_("not supported");
     130}
     131/*}}}*/
    126132/*FUNCTION PentaInput::SpawnResult{{{*/
    127133ElementResult* PentaInput::SpawnResult(int step, IssmDouble time){

  • issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h

    r16238 r16382  
    3535                int   InstanceEnum();
    3636                Input* SpawnTriaInput(int location);
     37                Input* SpawnSegInput(int index1,int index2);
    3738                Input* PointwiseDivide(Input* inputB);
    3839                Input* PointwiseMin(Input* inputB);
    … …  
    5152                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");};
    5253                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error_("not implemented yet");};
     54                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
    5355                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");};
    5456                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss);

  • issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp

    r16151 r16382  
    120120        xMemCpy(outinput->timesteps,this->timesteps,this->numtimesteps);
    121121        outinput->inputs=(Inputs*)this->inputs->SpawnTriaInputs(location);
     122        outinput->parameters=this->parameters;
     123
     124        /*Assign output*/
     125        return outinput;
     126
     127}
     128/*}}}*/
     129/*FUNCTION TransientInput::SpawnSegInput{{{*/
     130Input* TransientInput::SpawnSegInput(int index1,int index2){
     131
     132        /*output*/
     133        TransientInput* outinput=NULL;
     134
     135        /*Create new Transientinput (copy of current input)*/
     136        outinput=new TransientInput();
     137        outinput->enum_type=this->enum_type;
     138        outinput->numtimesteps=this->numtimesteps;
     139        outinput->timesteps=xNew<IssmDouble>(this->numtimesteps);
     140        xMemCpy(outinput->timesteps,this->timesteps,this->numtimesteps);
     141        outinput->inputs=(Inputs*)this->inputs->SpawnSegInputs(index1,index2);
    122142        outinput->parameters=this->parameters;
    123143

  • issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h

    r16238 r16382  
    3939                int    InstanceEnum();
    4040                Input* SpawnTriaInput(int location);
     41                Input* SpawnSegInput(int index1,int index2);
    4142                Input* PointwiseDivide(Input* forcingB){_error_("not implemented yet");};
    4243                Input* PointwiseMin(Input* forcingB){_error_("not implemented yet");};
    … …  
    5556                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");};
    5657                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error_("not implemented yet");};
     58                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
    5759                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss);
    5860                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp

    r16238 r16382  
    9797}
    9898/*}}}*/
     99/*FUNCTION SegInput::SpawnSegInput{{{*/
     100Input* TriaInput::SpawnSegInput(int index1,int index2){
     101
     102        /*output*/
     103        SegInput* outinput=NULL;
     104        IssmDouble newvalues[2]; //Assume P1 interpolation only for now
     105
     106        /*Create arrow of indices depending on location (0=base 1=surface)*/
     107
     108        newvalues[0]=this->values[index1];
     109        newvalues[1]=this->values[index2];
     110
     111        /*Create new Seg input*/
     112        outinput=new SegInput(this->enum_type,&newvalues[0],P1Enum);
     113
     114        /*Assign output*/
     115        return outinput;
     116
     117}
     118/*}}}*/
    99119/*FUNCTION TriaInput::SpawnResult{{{*/
    100120ElementResult* TriaInput::SpawnResult(int step, IssmDouble time){

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h

    r16238 r16382  
    3535                int    InstanceEnum();
    3636                Input* SpawnTriaInput(int location);
     37                Input* SpawnSegInput(int index1,int index2);
    3738                Input* PointwiseDivide(Input* inputB);
    3839                Input* PointwiseMin(Input* inputB);
    … …  
    5253                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");};
    5354                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int index){_error_("not implemented yet");};
     55                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
    5456                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss);
    5557                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");};

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

    r16343 r16382  
    9595                                analysis_type==HydrologyDCEfficientAnalysisEnum
    9696                                ){
    97                 if(iomodel->meshtype==Mesh3DEnum){
     97                if(iomodel->meshtype==Mesh3DEnum || iomodel->meshtype==Mesh2DverticalEnum){
    9898                        /*On a 3d mesh, we may have collapsed elements, hence dead nodes. Freeze them out: */
    9999                        _assert_(iomodel->Data(MeshVertexonbedEnum));

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

    r16343 r16382  
    5959#include "./Inputs/PentaInput.h"
    6060#include "./Inputs/TriaInput.h"
     61#include "./Inputs/SegInput.h"
    6162#include "./Inputs/ControlInput.h"
    6263#include "./Inputs/DatasetInput.h"

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/BedSlope/CreateNodesBedSlope.cpp

    r16291 r16382  
    1111void    CreateNodesBedSlope(Nodes** pnodes, IoModel* iomodel){
    1212
    13         if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     13        if(iomodel->meshtype==Mesh3DEnum){
     14                iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     15        }
     16        else if(iomodel->meshtype==Mesh2DverticalEnum){
     17                iomodel->FetchData(1,MeshVertexonbedEnum);
     18        }
    1419        CreateNodes(pnodes,iomodel,BedSlopeAnalysisEnum,P1Enum);
    1520        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/BedSlope/UpdateElementsBedSlope.cpp

    r16291 r16382  
    2828                iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
    2929        }
     30        if(iomodel->meshtype==Mesh2DverticalEnum){
     31                iomodel->FetchDataToInput(elements,MeshVertexonbedEnum);
     32        }
    3033}

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/SurfaceSlope/CreateNodesSurfaceSlope.cpp

    r16291 r16382  
    1111void    CreateNodesSurfaceSlope(Nodes** pnodes, IoModel* iomodel){
    1212
    13         if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     13        if(iomodel->meshtype==Mesh3DEnum){
     14                iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
     15        }
     16        else if(iomodel->meshtype==Mesh2DverticalEnum){
     17                iomodel->FetchData(1,MeshVertexonbedEnum);
     18        }
    1419        CreateNodes(pnodes,iomodel,SurfaceSlopeAnalysisEnum,P1Enum);
    1520        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
Note: See TracChangeset for help on using the changeset viewer.