Changeset 24089


Ignore:
Timestamp:
07/15/19 04:40:13 (6 years ago)
Author:
Mathieu Morlighem
Message:

CHG: working on P0DG

Location:
issm/trunk-jpl/src/c/classes
Files:
18 edited

Legend:

Unmodified
Added
Removed

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

    r24049 r24089  
    13181318
    13191319}/*}}}*/
     1320Node*      Element::GetNode(int nodeindex){/*{{{*/
     1321   _assert_(nodeindex>=0);
     1322   _assert_(nodeindex<this->GetNumberOfNodes(this->element_type));
     1323   return this->nodes[nodeindex];
     1324}/*}}}*/
     1325int        Element::GetNodeIndex(Node* node){/*{{{*/
     1326
     1327        _assert_(this->nodes);
     1328        int numnodes = this->GetNumberOfNodes(this->element_type);
     1329
     1330        for(int i=0;i<numnodes;i++){
     1331                if(node==nodes[i]) return i;
     1332        }
     1333        _error_("Node provided not found among element nodes");
     1334
     1335}
     1336/*}}}*/
    13201337void       Element::GetNodesLidList(int* lidlist){/*{{{*/
    13211338

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

    r24049 r24089  
    9898                void               GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type);
    9999                void               GetInputsInterpolations(Vector<IssmDouble>* interps);
     100                Node*              GetNode(int nodeindex);
     101                int                GetNodeIndex(Node* node);
    100102                void               GetNodesLidList(int* lidlist);
    101103                void               GetNodesSidList(int* sidlist);
    … …  
    233235                virtual IssmDouble GetIcefrontArea(){_error_("not implemented");};
    234236                virtual void       GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum)=0;
    235                 virtual void       GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
     237                virtual void       GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){_error_("not implemented yet");};
     238                virtual void       GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
    236239                virtual void       GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level)=0;
    237240                virtual void       GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues)=0;
    238                 virtual Node*      GetNode(int node_number)=0;
    239                 virtual int        GetNodeIndex(Node* node)=0;
     241                virtual int        GetVertexIndex(Vertex* vertex){_error_("not implemented");};;
    240242                virtual int        GetNumberOfNodes(void)=0;
    241243                virtual int        GetNumberOfNodes(int enum_type)=0;

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

    r24021 r24089  
    12211221        if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
    12221222
     1223        int index = this->GetNodeIndex(node);
     1224
    12231225        GaussPenta* gauss=new GaussPenta();
    1224         gauss->GaussVertex(this->GetNodeIndex(node));
     1226        gauss->GaussNode(this->element_type,index);
     1227
     1228        input->GetInputValue(pvalue,gauss);
     1229        delete gauss;
     1230}
     1231/*}}}*/
     1232void       Penta::GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){/*{{{*/
     1233
     1234        Input* input=inputs->GetInput(enumtype);
     1235        if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
     1236
     1237        int index = this->GetVertexIndex(vertex);
     1238
     1239        GaussPenta* gauss=new GaussPenta();
     1240        gauss->GaussVertex(index);
    12251241
    12261242        input->GetInputValue(pvalue,gauss);
    … …  
    13131329}
    13141330/*}}}*/
    1315 Node*      Penta::GetNode(int node_number){/*{{{*/
    1316         _assert_(node_number>=0);
    1317         _assert_(node_number<this->NumberofNodes(this->element_type));
    1318         return this->nodes[node_number];
    1319 }
    1320 /*}}}*/
    1321 int        Penta::GetNodeIndex(Node* node){/*{{{*/
    1322 
    1323         _assert_(nodes);
    1324         int numnodes = this->NumberofNodes(this->element_type);
    1325 
    1326         for(int i=0;i<numnodes;i++){
    1327                 if(node==nodes[i]) return i;
    1328         }
    1329         _error_("Node provided not found among element nodes");
    1330 
     1331int        Penta::GetVertexIndex(Vertex* vertex){/*{{{*/
     1332
     1333        _assert_(vertices);
     1334        for(int i=0;i<NUMVERTICES;i++){
     1335                if(vertex==vertices[i])
     1336                 return i;
     1337        }
     1338        _error_("Vertex provided not found among element vertices");
    13311339}
    13321340/*}}}*/

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

    r24020 r24089  
    7777                IssmDouble              GetIcefrontArea();
    7878                void           GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
     79                void           GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype);
    7980                void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    8081                void           GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented yet");};
    8182                void           GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues){_error_("not implemented yet");};
    8283                void                            GetLevelsetIntersectionBase(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level);
    83                 Node*          GetNode(int node_number);
    84                 int            GetNodeIndex(Node* node);
    8584                int            GetNumberOfNodes(void);
    8685                int            GetNumberOfNodes(int enum_type);
    … …  
    9089                void           GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data);
    9190                void           GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,int offset);
     91                int            GetVertexIndex(Vertex* vertex);
    9292                void           GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
    9393                void           GetVerticesCoordinatesTop(IssmDouble** pxyz_list);

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

    r24020 r24089  
    6363                IssmDouble  GetGroundedPortion(IssmDouble* xyz_list);
    6464                void               GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
    65                 void        GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
     65                void        GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){_error_("not implemented yet");};
    6666                void               GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented");};
    6767                void        GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues){_error_("not implemented yet");};
    68                 Node*       GetNode(int node_number){_error_("Not implemented");};
    69                 int         GetNodeIndex(Node* node){_error_("not implemented yet");};
    7068                int         GetNumberOfNodes(void);
    7169                int         GetNumberOfNodes(int enum_type){_error_("not implemented yet");};

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

    r23644 r24089  
    256256        input->GetInputValue(pvalue,gauss);
    257257        delete gauss;
    258 }
    259 /*}}}*/
    260 int      Tetra::GetNodeIndex(Node* node){/*{{{*/
    261 
    262         _assert_(nodes);
    263         int numnodes = this->NumberofNodes(this->element_type);
    264 
    265         for(int i=0;i<numnodes;i++){
    266                 if(node==nodes[i]) return i;
    267         }
    268         _error_("Node provided not found among element nodes");
    269 
    270258}
    271259/*}}}*/

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

    r24020 r24089  
    7272                void               GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented yet");};
    7373                void        GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues){_error_("not implemented yet");};
    74                 Node*       GetNode(int node_number){_error_("Not implemented");};
    75                 int         GetNodeIndex(Node* node);
    7674                int         GetNumberOfNodes(void);
    7775                int         GetNumberOfNodes(int enum_type){_error_("not implemented yet");};

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

    r24088 r24089  
    18321832        if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
    18331833
     1834        int index = this->GetNodeIndex(node);
     1835
    18341836        GaussTria* gauss=new GaussTria();
    1835         gauss->GaussVertex(this->GetNodeIndex(node));
     1837        gauss->GaussNode(this->element_type,index);
     1838
     1839        input->GetInputValue(pvalue,gauss);
     1840        delete gauss;
     1841}
     1842/*}}}*/
     1843void       Tria::GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){/*{{{*/
     1844
     1845        Input* input=inputs->GetInput(enumtype);
     1846        if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
     1847
     1848        int index = this->GetVertexIndex(vertex);
     1849
     1850        GaussTria* gauss=new GaussTria();
     1851        gauss->GaussVertex(index);
    18361852
    18371853        input->GetInputValue(pvalue,gauss);
    … …  
    20072023}
    20082024/*}}}*/
    2009 Node*      Tria::GetNode(int node_number){/*{{{*/
    2010         _assert_(node_number>=0);
    2011         _assert_(node_number<this->NumberofNodes(this->element_type));
    2012         return this->nodes[node_number];
    2013 
    2014 }/*}}}*/
    2015 int        Tria::GetNodeIndex(Node* node){/*{{{*/
    2016 
    2017         _assert_(nodes);
    2018         for(int i=0;i<NUMVERTICES;i++){
    2019                 if(node==nodes[i])
    2020                  return i;
    2021         }
    2022         _error_("Node provided not found among element nodes");
    2023 }
    2024 /*}}}*/
    20252025int        Tria::GetVertexIndex(Vertex* vertex){/*{{{*/
    20262026
    … …  
    20302030                 return i;
    20312031        }
    2032         _error_("Vertex provided not found among element nodes");
     2032        _error_("Vertex provided not found among element vertices");
    20332033}
    20342034/*}}}*/

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

    r24088 r24089  
    8484                void          GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
    8585                void          GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level);
    86                 int         GetNodeIndex(Node* node);
    8786                int         GetVertexIndex(Vertex* vertex);
    8887                int         GetNumberOfNodes(void);
    … …  
    180179                int            GetElementType(void);
    181180                void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    182                 void            GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level);
     181                void           GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype);
     182                void                  GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level);
    183183                void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    184                 Node*          GetNode(int node_number);
    185184                void             InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
    186185                void             InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){_error_("not implemented yet");};

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

    r23962 r24089  
    157157                case NoneEnum:
    158158                        return;
    159                 case P0Enum:
     159                case P0Enum: case P0DGEnum:
    160160                        basis[0]=1.;
    161161                        return;
    … …  
    369369
    370370        switch(finiteelement){
    371                 case P0Enum:
     371                case P0Enum: case P0DGEnum:
    372372                        /*Nodal function 1*/
    373373                        dbasis[NUMNODESP0*0+0] = 0.;
    … …  
    512512                case NoneEnum:                return 0;
    513513                case P0Enum:                  return NUMNODESP0;
     514                case P0DGEnum:                return NUMNODESP0;
    514515                case P1Enum:                  return NUMNODESP1;
    515516                case P1DGEnum:                return NUMNODESP1;

  • TabularUnified issm/trunk-jpl/src/c/classes/Loads/Channel.cpp

    r24080 r24089  
    355355        if(!tria->IsIceInElement()) return NULL;
    356356        _assert_(tria->FiniteElement()==P1Enum);
    357         int index1=tria->GetNodeIndex(nodes[0]);
    358         int index2=tria->GetNodeIndex(nodes[1]);
     357        int index1=tria->GetVertexIndex(vertices[0]);
     358        int index2=tria->GetVertexIndex(vertices[1]);
    359359
    360360        /*Intermediaries */
    … …  
    488488        if(!tria->IsIceInElement()) return NULL;
    489489        _assert_(tria->FiniteElement()==P1Enum);
    490         int index1=tria->GetNodeIndex(nodes[0]);
    491         int index2=tria->GetNodeIndex(nodes[1]);
     490        int index1=tria->GetVertexIndex(vertices[0]);
     491        int index2=tria->GetVertexIndex(vertices[1]);
    492492
    493493        /*Intermediaries */
    … …  
    605605        }
    606606        _assert_(tria->FiniteElement()==P1Enum);
    607         int index1=tria->GetNodeIndex(nodes[0]);
    608         int index2=tria->GetNodeIndex(nodes[1]);
     607        int index1=tria->GetVertexIndex(vertices[0]);
     608        int index2=tria->GetVertexIndex(vertices[1]);
    609609
    610610        /*Intermediaries */

  • TabularUnified issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp

    r24069 r24089  
    2323        this->parameters=NULL;
    2424        this->hnode=NULL;
     25        this->hvertex=NULL;
    2526        this->node=NULL;
    2627        this->helement=NULL;
    … …  
    4748
    4849        this->hnode=new Hook(&pengrid_node_id,1);
     50        this->hvertex=new Hook(&pengrid_node_id,1);
    4951        this->helement=new Hook(&pengrid_element_id,1);
    5052
    … …  
    5254        this->parameters=NULL;
    5355        this->node=NULL;
     56        this->vertex=NULL;
    5457        this->element=NULL;
    5558}
    … …  
    5760Moulin::~Moulin(){/*{{{*/
    5861        delete hnode;
     62        delete hvertex;
    5963        delete helement;
    6064        return;
    … …  
    7781        /*now deal with hooks and objects: */
    7882        pengrid->hnode=(Hook*)this->hnode->copy();
     83        pengrid->hvertex=(Hook*)this->hvertex->copy();
    7984        pengrid->helement=(Hook*)this->helement->copy();
    8085
    8186        /*corresponding fields*/
    8287        pengrid->node  =(Node*)pengrid->hnode->delivers();
     88        pengrid->vertex=(Vertex*)pengrid->hvertex->delivers();
    8389        pengrid->element=(Element*)pengrid->helement->delivers();
    8490
    … …  
    9197        _printf_("   id: " << id << "\n");
    9298        hnode->DeepEcho();
     99        hvertex->DeepEcho();
    93100        helement->DeepEcho();
    94101        _printf_("   parameters\n");
    … …  
    112119        if(marshall_direction==MARSHALLING_BACKWARD){
    113120                this->hnode      = new Hook();
     121                this->hvertex      = new Hook();
    114122                this->helement   = new Hook();
    115123        }
    116124
    117125        this->hnode->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
     126        this->hvertex->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
    118127        this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
    119128
    120129        /*corresponding fields*/
    121130        node   =(Node*)this->hnode->delivers();
     131        vertex =(Vertex*)this->hvertex->delivers();
    122132        element=(Element*)this->helement->delivers();
    123133}
    … …  
    135145         * datasets, using internal ids and offsets hidden in hooks: */
    136146        hnode->configure(nodesin);
     147        hvertex->configure(verticesin);
    137148        helement->configure(elementsin);
    138149
    139150        /*Get corresponding fields*/
    140151        node=(Node*)hnode->delivers();
     152        vertex=(Vertex*)hvertex->delivers();
    141153        element=(Element*)helement->delivers();
    142154
    … …  
    252264        /*Get Element type*/
    253265        this->hnode->reset();
     266        this->hvertex->reset();
    254267        this->helement->reset();
    255268

  • TabularUnified issm/trunk-jpl/src/c/classes/Loads/Moulin.h

    r24053 r24089  
    2626
    2727                /*Hooks*/
    28                 Hook* hnode;  //hook to 1 node
     28                Hook* hnode;     //hook to 1 node
     29                Hook* hvertex;   //hook to 1 vertex
    2930                Hook* helement;  //hook to 1 element
    3031
    3132                /*Corresponding fields*/
    3233                Node    *node;
     34                Vertex  *vertex;
    3335                Element *element;
    3436

  • TabularUnified issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp

    r23970 r24089  
    364364
    365365        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    366         int index1=tria->GetNodeIndex(nodes[0]);
    367         int index2=tria->GetNodeIndex(nodes[1]);
     366        int index1=tria->GetVertexIndex(vertices[0]);
     367        int index2=tria->GetVertexIndex(vertices[1]);
    368368
    369369        /* Start  looping on the number of gaussian points: */
    … …  
    408408
    409409        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    410         int index1=tria->GetNodeIndex(nodes[0]);
    411         int index2=tria->GetNodeIndex(nodes[1]);
     410        int index1=tria->GetVertexIndex(vertices[0]);
     411        int index2=tria->GetVertexIndex(vertices[1]);
    412412
    413413        /* Start  looping on the number of gaussian points: */

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

    r24088 r24089  
    3232        /* Intermediary */
    3333        int pos1,pos2,pos3,pos4;
    34         int num_nodes;
     34        int numnodes;
    3535
    3636        /*numericalflux constructor data: */
    … …  
    5050        if(e2==-1){
    5151                /* Boundary edge, only one element */
    52                 num_nodes=2;
    5352                numericalflux_type=BoundaryEnum;
    5453                numericalflux_elem_ids[0]=e1;
    … …  
    5655        else{
    5756                /* internal edge: connected to 2 elements */
    58       num_nodes=4;
    5957                numericalflux_type=InternalEnum;
    6058                numericalflux_elem_ids[0]=e1;
    … …  
    6462   /*FIXME: hardcode element degree for now*/
    6563   this->flux_degree= P1DGEnum;
     64   //printf("-------------- file: Numericalflux.cpp line: %i\n",__LINE__);
     65   //this->flux_degree= P0DGEnum;
    6666
    6767        /*1: Get vertices ids*/
    … …  
    103103        }
    104104
     105   switch(this->flux_degree){
     106      case P0DGEnum:
     107         if(numericalflux_type==InternalEnum) numnodes = 2;
     108         else numnodes = 1;
     109                        for(int i=0;i<numnodes;i++) numericalflux_node_ids[i] = numericalflux_elem_ids[i];
     110         numericalflux_node_ids[1] = numericalflux_elem_ids[1];
     111         break;
     112      case P1DGEnum:
     113         if(numericalflux_type==InternalEnum) numnodes = 4;
     114         else numnodes = 2;
     115         for(int i=0;i<numnodes;i++) numericalflux_node_ids[i] = numericalflux_node_ids[i]; //FIXME: to be improved...
     116         break;
     117      default:
     118         _error_("not supported yet");
     119
     120   }
     121
    105122        /*Assign object fields: */
    106123        this->id          = numericalflux_id;
    … …  
    109126
    110127        /*Hooks: */
    111         this->hnodes    = new Hook(numericalflux_node_ids,num_nodes);
     128        this->hnodes    = new Hook(numericalflux_node_ids,numnodes);
    112129        this->hvertices = new Hook(&numericalflux_vertex_ids[0],2);
    113130        this->helement  = new Hook(numericalflux_elem_ids,1); // take only the first element for now
    … …  
    363380void  Numericalflux::ResetHooks(){/*{{{*/
    364381
    365         this->nodes=NULL;
    366         this->vertices=NULL;
    367         this->element=NULL;
    368         this->parameters=NULL;
     382        this->nodes      = NULL;
     383        this->vertices   = NULL;
     384        this->element    = NULL;
     385        this->parameters = NULL;
    369386
    370387        /*Get Element type*/
    … …  
    376393/*}}}*/
    377394void  Numericalflux::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
     395   /*Nothing to do :)*/
    378396
    379397}

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

    r23959 r24089  
    2323        this->parameters=NULL;
    2424        this->hnodes=NULL;
     25        this->hvertices=NULL;
    2526        this->helements=NULL;
    2627        this->nodes=NULL;
     28        this->vertices=NULL;
    2729        this->elements=NULL;
    2830}
    … …  
    6365        /*Hooks: */
    6466        this->hnodes=new Hook(riftfront_node_ids,2);
     67        this->hvertices=new Hook(riftfront_node_ids,2);
    6568        this->helements=new Hook(riftfront_elem_ids,2);
    6669
    … …  
    8891        this->parameters=NULL;
    8992        this->nodes= NULL;
     93        this->vertices= NULL;
    9094        this->elements= NULL;
    9195
    … …  
    9599        this->parameters=NULL;
    96100        delete hnodes;
     101        delete hvertices;
    97102        delete helements;
    98103}
    … …  
    119124        /*now deal with hooks and objects: */
    120125        riftfront->hnodes=(Hook*)this->hnodes->copy();
     126        riftfront->hvertices=(Hook*)this->hvertices->copy();
    121127        riftfront->helements=(Hook*)this->helements->copy();
    122128
    123129        /*corresponding fields*/
    124130        riftfront->nodes   =(Node**)riftfront->hnodes->deliverp();
     131        riftfront->vertices=(Vertex**)riftfront->hvertices->deliverp();
    125132        riftfront->elements=(Element**)riftfront->helements->deliverp();
    126133
    … …  
    147154        _printf_("   id: " << id << "\n");
    148155        hnodes->DeepEcho();
     156        hvertices->DeepEcho();
    149157        helements->DeepEcho();
    150158        _printf_("   parameters\n");
    … …  
    157165        _printf_("   id: " << id << "\n");
    158166        _printf_("   hnodes: " << hnodes << "\n");
     167        _printf_("   hvertices: " << hvertices << "\n");
    159168        _printf_("   helements: " << helements << "\n");
    160169        _printf_("   parameters: " << parameters << "\n");
    … …  
    193202        if(marshall_direction==MARSHALLING_BACKWARD){
    194203                this->hnodes      = new Hook();
     204                this->hvertices      = new Hook();
    195205                this->helements   = new Hook();
    196206        }
    197207
    198208        this->hnodes->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
     209        this->hvertices->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
    199210        this->helements->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
    200211
    201212        /*corresponding fields*/
    202213        nodes     =(Node**)this->hnodes->deliverp();
     214        vertices  =(Vertex**)this->hvertices->deliverp();
    203215        elements  =(Element**)this->helements->deliverp();
    204216
    … …  
    245257         * datasets, using internal ids and offsets hidden in hooks: */
    246258        hnodes->configure(nodesin);
     259        hvertices->configure(verticesin);
    247260        helements->configure(elementsin);
    248261
    249262        /*Initialize hooked fields*/
    250263        this->nodes   =(Node**)hnodes->deliverp();
     264        this->vertices=(Vertex**)hvertices->deliverp();
    251265        this->elements=(Element**)helements->deliverp();
    252266
    … …  
    344358
    345359        this->nodes=NULL;
     360        this->vertices=NULL;
    346361        this->elements=NULL;
    347362        this->parameters=NULL;
    … …  
    349364        /*Get Element type*/
    350365        this->hnodes->reset();
     366        this->hvertices->reset();
    351367        this->helements->reset();
    352368
    … …  
    420436        IssmDouble  penalty_offset;
    421437
    422         /*Objects: */
    423         Tria       *tria1               = NULL;
    424         Tria       *tria2               = NULL;
    425 
    426438        /*enum of element? */
    427439        if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
    428         tria1=(Tria*)elements[0];
    429         tria2=(Tria*)elements[1];
     440        Tria* tria1=(Tria*)elements[0];
     441        Tria* tria2=(Tria*)elements[1];
    430442
    431443        /*Initialize Element Matrix*/
    … …  
    435447        /*Get some parameters: */
    436448        this->parameters->FindParam(&penalty_offset,StressbalancePenaltyFactorEnum);
    437         tria1->GetInputValue(&h[0],nodes[0],ThicknessEnum);
    438         tria2->GetInputValue(&h[1],nodes[1],ThicknessEnum);
     449        tria1->GetInputValue(&h[0],vertices[0],ThicknessEnum);
     450        tria2->GetInputValue(&h[1],vertices[1],ThicknessEnum);
    439451        if (h[0]!=h[1])_error_("different thicknesses not supported for rift fronts");
    440452        thickness=h[0];
    … …  
    506518        IssmDouble pressure_water;
    507519
    508         /*Objects: */
    509         Tria *tria1 = NULL;
    510         Tria *tria2 = NULL;
    511 
    512520        /*enum of element? */
    513521        if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
    514         tria1=(Tria*)elements[0];
    515         tria2=(Tria*)elements[1];
     522        Tria* tria1=(Tria*)elements[0];
     523        Tria* tria2=(Tria*)elements[1];
    516524
    517525        /*Initialize Element Matrix*/
    … …  
    523531        rho_water=tria1->FindParam(MaterialsRhoSeawaterEnum);
    524532        gravity=tria1->FindParam(ConstantsGEnum);
    525         tria1->GetInputValue(&h[0],nodes[0],ThicknessEnum);
    526         tria2->GetInputValue(&h[1],nodes[1],ThicknessEnum);
     533        tria1->GetInputValue(&h[0],vertices[0],ThicknessEnum);
     534        tria2->GetInputValue(&h[1],vertices[1],ThicknessEnum);
    527535        if (h[0]!=h[1])_error_("different thicknesses not supported for rift fronts");
    528536        thickness=h[0];
    529         tria1->GetInputValue(&b[0],nodes[0],BaseEnum);
    530         tria2->GetInputValue(&b[1],nodes[1],BaseEnum);
     537        tria1->GetInputValue(&b[0],vertices[0],BaseEnum);
     538        tria2->GetInputValue(&b[1],vertices[1],BaseEnum);
    531539        if (b[0]!=b[1])_error_("different beds not supported for rift fronts");
    532540        bed=b[0];
    … …  
    571579
    572580        /*Ok, add contribution to first node, along the normal i==0: */
    573         for (j=0;j<2;j++){
     581        for(int j=0;j<2;j++){
    574582                pe->values[j]+=pressure*normal[j]*length;
    575583        }
    576584
    577585        /*Add contribution to second node, along the opposite normal: i==1 */
    578         for (j=0;j<2;j++){
     586        for(int j=0;j<2;j++){
    579587                pe->values[2+j]+= -pressure*normal[j]*length;
    580588        }       
    … …  
    624632
    625633        /*First recover velocity: */
    626         tria1->GetInputValue(&vx1,nodes[0],VxEnum);
    627         tria2->GetInputValue(&vx2,nodes[1],VxEnum);
    628         tria1->GetInputValue(&vy1,nodes[0],VyEnum);
    629         tria2->GetInputValue(&vy2,nodes[1],VyEnum);
     634        tria1->GetInputValue(&vx1,vertices[0],VxEnum);
     635        tria2->GetInputValue(&vx2,vertices[1],VxEnum);
     636        tria1->GetInputValue(&vy1,vertices[0],VyEnum);
     637        tria2->GetInputValue(&vy2,vertices[1],VyEnum);
    630638
    631639        /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */

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

    r23959 r24089  
    2828                /*hooks: */
    2929                Hook* hnodes;
     30                Hook* hvertices;
    3031                Hook* helements;
    3132
    3233                /*Corresponding fields*/
    3334                Node    **nodes;
     35                Vertex  **vertices;
    3436                Element **elements;
    3537

  • TabularUnified issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp

    r21892 r24089  
    497497        /*update static arrays*/
    498498        switch(finiteelement){
    499                 case P0Enum:
     499                case P0Enum: case P0DGEnum:
    500500                        switch(iv){
    501501                                case 0: coord1=1./3.; coord2=1./3.; coord3=1./3.; break;
Note: See TracChangeset for help on using the changeset viewer.