Changeset 24088


Ignore:
Timestamp:
07/14/19 19:14:42 (6 years ago)
Author:
Mathieu Morlighem
Message:

CHG: improved numerical flux calculation and working on adding P0DG

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

Legend:

Unmodified
Added
Removed

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

    r24049 r24088  
    20212021        }
    20222022        _error_("Node provided not found among element nodes");
     2023}
     2024/*}}}*/
     2025int        Tria::GetVertexIndex(Vertex* vertex){/*{{{*/
     2026
     2027        _assert_(vertices);
     2028        for(int i=0;i<NUMVERTICES;i++){
     2029                if(vertex==vertices[i])
     2030                 return i;
     2031        }
     2032        _error_("Vertex provided not found among element nodes");
    20232033}
    20242034/*}}}*/
    … …  
    44294439        /*Recover nodes ids needed to initialize the node hook.*/
    44304440        switch(finiteelement_type){
     4441                case P0DGEnum:
     4442                        numnodes        = 1;
     4443                        tria_node_ids   = xNew<int>(numnodes);
     4444                        tria_node_ids[0]= index + 1;
     4445                        break;
    44314446                case P1Enum:
    44324447                        numnodes        = 3;

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

    r24049 r24088  
    8585                void          GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level);
    8686                int         GetNodeIndex(Node* node);
     87                int         GetVertexIndex(Vertex* vertex);
    8788                int         GetNumberOfNodes(void);
    8889                int         GetNumberOfNodes(int enum_type);

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

    r23959 r24088  
    1717/*Load macros*/
    1818#define NUMVERTICES 2
    19 #define NUMNODES_INTERNAL 4
    20 #define NUMNODES_BOUNDARY 2
    2119
    2220/*Numericalflux constructors and destructor*/
    … …  
    3331
    3432        /* Intermediary */
    35         int  j;
    36         int  pos1,pos2,pos3,pos4;
    37         int  num_nodes;
     33        int pos1,pos2,pos3,pos4;
     34        int num_nodes;
    3835
    3936        /*numericalflux constructor data: */
    40         int   numericalflux_elem_ids[2];
    41         int   numericalflux_vertex_ids[2];
    42         int   numericalflux_node_ids[4];
    43         int   numericalflux_type;
     37        int numericalflux_elem_ids[2];
     38        int numericalflux_vertex_ids[2];
     39        int numericalflux_node_ids[4];
     40        int numericalflux_type;
     41   int numericalflux_degree;
    4442
    4543        /*Get edge*/
    … …  
    5856        else{
    5957                /* internal edge: connected to 2 elements */
    60                 num_nodes=4;
     58      num_nodes=4;
    6159                numericalflux_type=InternalEnum;
    6260                numericalflux_elem_ids[0]=e1;
    … …  
    6462        }
    6563
     64   /*FIXME: hardcode element degree for now*/
     65   this->flux_degree= P1DGEnum;
     66
    6667        /*1: Get vertices ids*/
    6768        numericalflux_vertex_ids[0]=i1;
    … …  
    6970
    7071        /*2: Get node ids*/
    71         if (numericalflux_type==InternalEnum){
    72 
    73                 /*Now, we must get the nodes of the 4 nodes located on the edge*/
    74 
    75                 /*2: Get the column where these ids are located in the index*/
     72        if(numericalflux_type==InternalEnum){
     73                /*Get the column where these ids are located in the index*/
    7674                pos1=pos2=pos3=pos4=UNDEF;
    77                 for(j=0;j<3;j++){
     75                for(int j=0;j<3;j++){
    7876                        if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
    7977                        if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
    … …  
    8381                _assert_(pos1!=UNDEF && pos2!=UNDEF && pos3!=UNDEF && pos4!=UNDEF);
    8482
    85                 /*3: We have the id of the elements and the position of the vertices in the index
     83                /* We have the id of the elements and the position of the vertices in the index
    8684                 * we can compute their dofs!*/
    8785                numericalflux_node_ids[0]=3*(e1-1)+pos1;
    … …  
    9189        }
    9290        else{
    93 
    94                 /*2: Get the column where these ids are located in the index*/
     91                /*Get the column where these ids are located in the index*/
    9592                pos1=pos2=UNDEF;
    96                 for(j=0;j<3;j++){
     93                for(int j=0;j<3;j++){
    9794                        if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
    9895                        if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
    … …  
    10097                _assert_(pos1!=UNDEF && pos2!=UNDEF);
    10198
    102                 /*3: We have the id of the elements and the position of the vertices in the index
     99                /* We have the id of the elements and the position of the vertices in the index
    103100                 * we can compute their dofs!*/
    104101                numericalflux_node_ids[0]=3*(e1-1)+pos1;
    … …  
    106103        }
    107104
    108         /*Ok, we have everything to build the object: */
    109         this->id=numericalflux_id;
    110         this->flux_type = numericalflux_type;
     105        /*Assign object fields: */
     106        this->id          = numericalflux_id;
     107        this->flux_type   = numericalflux_type;
     108   this->flux_degree = numericalflux_degree;
    111109
    112110        /*Hooks: */
    113         this->hnodes    =new Hook(numericalflux_node_ids,num_nodes);
    114         this->hvertices =new Hook(&numericalflux_vertex_ids[0],2);
    115         this->helement  =new Hook(numericalflux_elem_ids,1); // take only the first element for now
    116 
    117         //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
    118         this->parameters=NULL;
    119         this->element=NULL;
    120         this->nodes=NULL;
     111        this->hnodes    = new Hook(numericalflux_node_ids,num_nodes);
     112        this->hvertices = new Hook(&numericalflux_vertex_ids[0],2);
     113        this->helement  = new Hook(numericalflux_elem_ids,1); // take only the first element for now
     114
     115        /*other fields*/
     116        this->parameters = NULL;
     117        this->element    = NULL;
     118        this->nodes      = NULL;
    121119}
    122120/*}}}*/
    … …  
    139137        numericalflux->id=this->id;
    140138        numericalflux->flux_type=this->flux_type;
     139        numericalflux->flux_degree=this->flux_degree;
    141140
    142141        /*point parameters: */
    … …  
    161160        _printf_("   id: " << id << "\n");
    162161        _printf_("   flux_type: " << this->flux_type<< "\n");
     162        _printf_("   flux_degree: " << this->flux_degree<< "\n");
    163163        hnodes->DeepEcho();
    164164        hvertices->DeepEcho();
    … …  
    175175        _printf_("   id: " << id << "\n");
    176176        _printf_("   flux_type: " << this->flux_type<< "\n");
     177        _printf_("   flux_degree: " << this->flux_degree<< "\n");
    177178        hnodes->Echo();
    178179        hvertices->Echo();
    … …  
    193194        MARSHALLING(id);
    194195        MARSHALLING(flux_type);
     196        MARSHALLING(flux_degree);
    195197
    196198        if(marshall_direction==MARSHALLING_BACKWARD){
    … …  
    300302        _assert_(nodes);
    301303
    302         switch(this->flux_type){
    303                 case InternalEnum:
    304                         for(int i=0;i<NUMNODES_INTERNAL;i++) lidlist[i]=nodes[i]->Lid();
    305                         return;
    306                 case BoundaryEnum:
    307                         for(int i=0;i<NUMNODES_BOUNDARY;i++) lidlist[i]=nodes[i]->Lid();
    308                         return;
    309                 default:
    310                         _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
    311         }
     304        int numnodes = this->GetNumberOfNodes();
     305        for(int i=0;i<numnodes;i++) lidlist[i]=nodes[i]->Lid();
    312306}
    313307/*}}}*/
    … …  
    317311        _assert_(nodes);
    318312
    319         switch(this->flux_type){
    320                 case InternalEnum:
    321                         for(int i=0;i<NUMNODES_INTERNAL;i++) sidlist[i]=nodes[i]->Sid();
    322                         return;
    323                 case BoundaryEnum:
    324                         for(int i=0;i<NUMNODES_BOUNDARY;i++) sidlist[i]=nodes[i]->Sid();
    325                         return;
    326                 default:
    327                         _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
    328         }
     313        int numnodes = this->GetNumberOfNodes();
     314        for(int i=0;i<numnodes;i++) sidlist[i]=nodes[i]->Sid();
    329315}
    330316/*}}}*/
    331317int   Numericalflux::GetNumberOfNodes(void){/*{{{*/
    332318
    333         switch(this->flux_type){
    334                 case InternalEnum:
    335                         return NUMNODES_INTERNAL;
    336                 case BoundaryEnum:
    337                         return NUMNODES_BOUNDARY;
    338                 default:
    339                         _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
     319        if(this->flux_degree==P0DGEnum){
     320                switch(this->flux_type){
     321                        case InternalEnum:
     322                                return 2;
     323                        case BoundaryEnum:
     324                                return 1;
     325                        default:
     326                                _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
     327                }
     328        }
     329        else if(this->flux_degree==P1DGEnum){
     330                switch(this->flux_type){
     331                        case InternalEnum:
     332                                return 4;
     333                        case BoundaryEnum:
     334                                return 2;
     335                        default:
     336                                _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
     337                }
     338        }
     339        else{
     340                _error_("Numericalflux " << EnumToStringx(this->flux_degree) << " not supported yet");
    340341        }
    341342
    … …  
    474475ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessBoundary(void){/*{{{*/
    475476
    476         /* constants*/
    477         const int numdof=NDOF1*NUMNODES_BOUNDARY;
     477        /*Initialize Element matrix and return if necessary*/
     478        Tria* tria=(Tria*)element;
     479        if(!tria->IsIceInElement()) return NULL;
    478480
    479481        /* Intermediaries*/
    480         int        i,j,ig,index1,index2;
    481         IssmDouble     DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN;
    482         IssmDouble     xyz_list[NUMVERTICES][3];
    483         IssmDouble     normal[2];
    484         IssmDouble     L[numdof];
    485         IssmDouble     Ke_g[numdof][numdof];
    486         GaussTria *gauss;
    487 
    488         /*Initialize Element matrix and return if necessary*/
    489         ElementMatrix* Ke = NULL;
    490         Tria*  tria=(Tria*)element;
    491         if(!tria->IsIceInElement()) return NULL;
     482        IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy;
     483        IssmDouble xyz_list[NUMVERTICES][3];
     484        IssmDouble normal[2];
    492485
    493486        /*Retrieve all inputs and parameters*/
    … …  
    498491
    499492        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    500         index1=tria->GetNodeIndex(nodes[0]);
    501         index2=tria->GetNodeIndex(nodes[1]);
    502 
    503         gauss=new GaussTria();
     493        int index1=tria->GetVertexIndex(vertices[0]);
     494        int index2=tria->GetVertexIndex(vertices[1]);
     495
     496        GaussTria* gauss=new GaussTria();
    504497        gauss->GaussEdgeCenter(index1,index2);
    505498        vxaverage_input->GetInputValue(&mean_vx,gauss);
    … …  
    507500        delete gauss;
    508501
    509         UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    510         if (UdotN<=0){
     502        IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
     503        if(UdotN<=0){
    511504                return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
    512505        }
    513         else{
    514                 Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters);
    515         }
     506
     507        /*Initialize Element vector and other vectors*/
     508   int            numnodes = this->GetNumberOfNodes();
     509   ElementMatrix *Ke       = new ElementMatrix(nodes,numnodes,this->parameters);
     510   IssmDouble    *basis    = xNew<IssmDouble>(numnodes);
    516511
    517512        /* Start  looping on the number of gaussian points: */
    518513        gauss=new GaussTria(index1,index2,2);
    519         for(ig=gauss->begin();ig<gauss->end();ig++){
     514        for(int ig=gauss->begin();ig<gauss->end();ig++){
    520515
    521516                gauss->GaussPoint(ig);
    522517
    523                 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
     518                tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
    524519
    525520                vxaverage_input->GetInputValue(&vx,gauss);
    … …  
    529524                DL=gauss->weight*Jdet*UdotN;
    530525
    531                 TripleMultiply(&L[0],1,numdof,1,
    532                                         &DL,1,1,0,
    533                                         &L[0],1,numdof,0,
    534                                         &Ke_g[0][0],0);
    535 
    536                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
     526                for(int i=0;i<numnodes;i++){
     527                        for(int j=0;j<numnodes;j++){
     528                                Ke->values[i*numnodes+j]+=DL*basis[i]*basis[j];
     529                        }
     530                }
    537531        }
    538532
    539533        /*Clean up and return*/
     534        xDelete<IssmDouble>(basis);
    540535        delete gauss;
    541536        return Ke;
    … …  
    543538/*}}}*/
    544539ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessInternal(void){/*{{{*/
    545 
    546         /* constants*/
    547         const int numdof=NDOF1*NUMNODES_INTERNAL;
    548 
    549         /* Intermediaries*/
    550         int        i,j,ig,index1,index2;
    551         IssmDouble     DL1,DL2,Jdet,vx,vy,UdotN;
    552         IssmDouble     xyz_list[NUMVERTICES][3];
    553         IssmDouble     normal[2];
    554         IssmDouble     B[numdof];
    555         IssmDouble     Bprime[numdof];
    556         IssmDouble     Ke_g1[numdof][numdof];
    557         IssmDouble     Ke_g2[numdof][numdof];
    558         GaussTria *gauss;
    559540
    560541        /*Initialize Element matrix and return if necessary*/
    561542        Tria*  tria=(Tria*)element;
    562543        if(!tria->IsIceInElement()) return NULL;
    563         ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES_INTERNAL,this->parameters);
     544
     545        /* Intermediaries*/
     546        IssmDouble DL1,DL2,Jdet,vx,vy,UdotN;
     547        IssmDouble xyz_list[NUMVERTICES][3];
     548        IssmDouble normal[2];
     549
     550        /*Fetch number of nodes for this flux*/
     551        int numnodes = this->GetNumberOfNodes();
     552
     553        /*Initialize variables*/
     554        ElementMatrix *Ke     = new ElementMatrix(nodes,numnodes,this->parameters);
     555        IssmDouble    *B      = xNew<IssmDouble>(numnodes);
     556        IssmDouble    *Bprime = xNew<IssmDouble>(numnodes);
    564557
    565558        /*Retrieve all inputs and parameters*/
    … …  
    570563
    571564        /* Start  looping on the number of gaussian points: */
    572         index1=tria->GetNodeIndex(nodes[0]);
    573         index2=tria->GetNodeIndex(nodes[1]);
    574         gauss=new GaussTria(index1,index2,2);
    575         for(ig=gauss->begin();ig<gauss->end();ig++){
     565        int index1=tria->GetVertexIndex(vertices[0]);
     566        int index2=tria->GetVertexIndex(vertices[1]);
     567        GaussTria* gauss=new GaussTria(index1,index2,2);
     568        for(int ig=gauss->begin();ig<gauss->end();ig++){
    576569
    577570                gauss->GaussPoint(ig);
    … …  
    587580                DL2=gauss->weight*Jdet*fabs(UdotN)/2;
    588581
    589                 TripleMultiply(&B[0],1,numdof,1,
    590                                         &DL1,1,1,0,
    591                                         &Bprime[0],1,numdof,0,
    592                                         &Ke_g1[0][0],0);
    593                 TripleMultiply(&B[0],1,numdof,1,
    594                                         &DL2,1,1,0,
    595                                         &B[0],1,numdof,0,
    596                                         &Ke_g2[0][0],0);
    597 
    598                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g1[i][j];
    599                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g2[i][j];
     582                for(int i=0;i<numnodes;i++){
     583                        for(int j=0;j<numnodes;j++){
     584                                Ke->values[i*numnodes+j]+=DL1*B[i]*Bprime[j];
     585                                Ke->values[i*numnodes+j]+=DL2*B[i]*B[j];
     586                        }
     587                }
    600588        }
    601589
    602590        /*Clean up and return*/
     591        xDelete<IssmDouble>(B);
     592        xDelete<IssmDouble>(Bprime);
    603593        delete gauss;
    604594        return Ke;
    … …  
    619609ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/
    620610
    621         /* constants*/
    622         const int numdof=NDOF1*NUMNODES_BOUNDARY;
    623 
    624         /* Intermediaries*/
    625         int        i,j,ig,index1,index2;
    626         IssmDouble     DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN;
    627         IssmDouble     xyz_list[NUMVERTICES][3];
    628         IssmDouble     normal[2];
    629         IssmDouble     L[numdof];
    630         IssmDouble     Ke_g[numdof][numdof];
    631         GaussTria *gauss;
    632 
    633611        /*Initialize Element matrix and return if necessary*/
    634         ElementMatrix* Ke = NULL;
    635612        Tria*  tria=(Tria*)element;
    636613        if(!tria->IsIceInElement()) return NULL;
    637614
     615        /* Intermediaries*/
     616        IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy;
     617        IssmDouble xyz_list[NUMVERTICES][3];
     618        IssmDouble normal[2];
     619
    638620        /*Retrieve all inputs and parameters*/
    639621        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    640         parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     622        IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
    641623        Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
    642624        Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
    … …  
    644626
    645627        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    646         index1=tria->GetNodeIndex(nodes[0]);
    647         index2=tria->GetNodeIndex(nodes[1]);
    648 
    649         gauss=new GaussTria();
     628        int index1=tria->GetVertexIndex(vertices[0]);
     629        int index2=tria->GetVertexIndex(vertices[1]);
     630
     631        GaussTria* gauss=new GaussTria();
    650632        gauss->GaussEdgeCenter(index1,index2);
    651633        vxaverage_input->GetInputValue(&mean_vx,gauss);
    … …  
    653635        delete gauss;
    654636
    655         UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    656         if (UdotN<=0){
     637        IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
     638        if(UdotN<=0){
    657639                return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
    658640        }
    659         else{
    660                 Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters);
    661         }
     641
     642        /*Initialize Element vector and other vectors*/
     643   int            numnodes = this->GetNumberOfNodes();
     644   ElementMatrix *Ke       = new ElementMatrix(nodes,numnodes,this->parameters);
     645   IssmDouble    *basis    = xNew<IssmDouble>(numnodes);
    662646
    663647        /* Start  looping on the number of gaussian points: */
    664648        gauss=new GaussTria(index1,index2,2);
    665         for(ig=gauss->begin();ig<gauss->end();ig++){
     649        for(int ig=gauss->begin();ig<gauss->end();ig++){
    666650
    667651                gauss->GaussPoint(ig);
    668652
    669                 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
     653                tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
    670654
    671655                vxaverage_input->GetInputValue(&vx,gauss);
    … …  
    675659                DL=gauss->weight*Jdet*dt*UdotN;
    676660
    677                 TripleMultiply(&L[0],1,numdof,1,
    678                                         &DL,1,1,0,
    679                                         &L[0],1,numdof,0,
    680                                         &Ke_g[0][0],0);
    681 
    682                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
     661                for(int i=0;i<numnodes;i++){
     662                        for(int j=0;j<numnodes;j++){
     663                                Ke->values[i*numnodes+j]+=DL*basis[i]*basis[j];
     664                        }
     665                }
    683666        }
    684667
    685668        /*Clean up and return*/
     669        xDelete<IssmDouble>(basis);
    686670        delete gauss;
    687671        return Ke;
    … …  
    689673/*}}}*/
    690674ElementMatrix* Numericalflux::CreateKMatrixMasstransportInternal(void){/*{{{*/
    691 
    692         /* constants*/
    693         const int numdof=NDOF1*NUMNODES_INTERNAL;
    694 
    695         /* Intermediaries*/
    696         int        i,j,ig,index1,index2;
    697         IssmDouble     DL1,DL2,Jdet,dt,vx,vy,UdotN;
    698         IssmDouble     xyz_list[NUMVERTICES][3];
    699         IssmDouble     normal[2];
    700         IssmDouble     B[numdof];
    701         IssmDouble     Bprime[numdof];
    702         IssmDouble     Ke_g1[numdof][numdof];
    703         IssmDouble     Ke_g2[numdof][numdof];
    704         GaussTria *gauss;
    705675
    706676        /*Initialize Element matrix and return if necessary*/
    707677        Tria*  tria=(Tria*)element;
    708678        if(!tria->IsIceInElement()) return NULL;
    709         ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES_INTERNAL,this->parameters);
     679
     680        /* Intermediaries*/
     681        IssmDouble DL1,DL2,Jdet,vx,vy,UdotN;
     682        IssmDouble xyz_list[NUMVERTICES][3];
     683        IssmDouble normal[2];
     684
     685        /*Fetch number of nodes for this flux*/
     686        int numnodes = this->GetNumberOfNodes();
     687
     688        /*Initialize variables*/
     689        ElementMatrix *Ke     = new ElementMatrix(nodes,numnodes,this->parameters);
     690        IssmDouble    *B      = xNew<IssmDouble>(numnodes);
     691        IssmDouble    *Bprime = xNew<IssmDouble>(numnodes);
    710692
    711693        /*Retrieve all inputs and parameters*/
    712694        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    713         parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     695        IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
    714696        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
    715697        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
    … …  
    717699
    718700        /* Start  looping on the number of gaussian points: */
    719         index1=tria->GetNodeIndex(nodes[0]);
    720         index2=tria->GetNodeIndex(nodes[1]);
    721         gauss=new GaussTria(index1,index2,2);
    722         for(ig=gauss->begin();ig<gauss->end();ig++){
     701        int index1=tria->GetVertexIndex(vertices[0]);
     702        int index2=tria->GetVertexIndex(vertices[1]);
     703        GaussTria* gauss=new GaussTria(index1,index2,2);
     704        for(int ig=gauss->begin();ig<gauss->end();ig++){
    723705
    724706                gauss->GaussPoint(ig);
    … …  
    734716                DL2=gauss->weight*Jdet*dt*fabs(UdotN)/2;
    735717
    736                 TripleMultiply(&B[0],1,numdof,1,
    737                                         &DL1,1,1,0,
    738                                         &Bprime[0],1,numdof,0,
    739                                         &Ke_g1[0][0],0);
    740                 TripleMultiply(&B[0],1,numdof,1,
    741                                         &DL2,1,1,0,
    742                                         &B[0],1,numdof,0,
    743                                         &Ke_g2[0][0],0);
    744 
    745                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g1[i][j];
    746                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g2[i][j];
     718                for(int i=0;i<numnodes;i++){
     719                        for(int j=0;j<numnodes;j++){
     720                                Ke->values[i*numnodes+j]+=DL1*B[i]*Bprime[j];
     721                                Ke->values[i*numnodes+j]+=DL2*B[i]*B[j];
     722                        }
     723                }
    747724        }
    748725
    749726        /*Clean up and return*/
     727   xDelete<IssmDouble>(B);
     728   xDelete<IssmDouble>(Bprime);
    750729        delete gauss;
    751730        return Ke;
    … …  
    772751ElementVector* Numericalflux::CreatePVectorBalancethicknessBoundary(void){/*{{{*/
    773752
    774         /* constants*/
    775         const int numdof=NDOF1*NUMNODES_BOUNDARY;
     753        /*Initialize Load Vector and return if necessary*/
     754        Tria*  tria=(Tria*)element;
     755        if(!tria->IsIceInElement()) return NULL;
    776756
    777757        /* Intermediaries*/
    778         int        i,ig,index1,index2;
    779         IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN,thickness;
     758        IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,thickness;
    780759        IssmDouble xyz_list[NUMVERTICES][3];
    781760        IssmDouble normal[2];
    782         IssmDouble L[numdof];
    783         GaussTria *gauss;
    784 
    785         /*Initialize Load Vector and return if necessary*/
    786         ElementVector* pe = NULL;
    787         Tria*  tria=(Tria*)element;
    788         if(!tria->IsIceInElement()) return NULL;
    789761
    790762        /*Retrieve all inputs and parameters*/
    791763        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    792         Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
    793         Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
    794         Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     764        Input* vxaverage_input = tria->inputs->GetInput(VxEnum);        _assert_(vxaverage_input);
     765        Input* vyaverage_input = tria->inputs->GetInput(VyEnum);        _assert_(vyaverage_input);
     766        Input* thickness_input = tria->inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    795767        GetNormal(&normal[0],xyz_list);
    796768
    797769        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    798         index1=tria->GetNodeIndex(nodes[0]);
    799         index2=tria->GetNodeIndex(nodes[1]);
    800 
    801         gauss=new GaussTria();
     770        int index1=tria->GetVertexIndex(vertices[0]);
     771        int index2=tria->GetVertexIndex(vertices[1]);
     772        GaussTria* gauss=new GaussTria();
    802773        gauss->GaussEdgeCenter(index1,index2);
    803774        vxaverage_input->GetInputValue(&mean_vx,gauss);
    804775        vyaverage_input->GetInputValue(&mean_vy,gauss);
    805776        delete gauss;
    806         UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    807         if (UdotN>0){
     777        IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
     778        if(UdotN>0){
    808779                return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
    809780        }
    810         else{
    811                 pe=new ElementVector(nodes,NUMNODES_BOUNDARY,this->parameters);
    812         }
     781
     782        /*Initialize Load Vector */
     783        int            numnodes = this->GetNumberOfNodes();
     784        ElementVector *pe       = new ElementVector(nodes,numnodes,this->parameters);
     785        IssmDouble    *basis    = xNew<IssmDouble>(numnodes);
    813786
    814787        /* Start  looping on the number of gaussian points: */
    815788        gauss=new GaussTria(index1,index2,2);
    816         for(ig=gauss->begin();ig<gauss->end();ig++){
     789        for(int ig=gauss->begin();ig<gauss->end();ig++){
    817790
    818791                gauss->GaussPoint(ig);
    819792
    820                 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
     793                tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
    821794
    822795                vxaverage_input->GetInputValue(&vx,gauss);
    … …  
    828801                DL= - gauss->weight*Jdet*UdotN*thickness;
    829802
    830                 for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
     803                for(int i=0;i<numnodes;i++) pe->values[i] += DL*basis[i];
    831804        }
    832805
    833806        /*Clean up and return*/
     807   xDelete<IssmDouble>(basis);
    834808        delete gauss;
    835809        return pe;
    … …  
    857831ElementVector* Numericalflux::CreatePVectorMasstransportBoundary(void){/*{{{*/
    858832
    859         /* constants*/
    860         const int numdof=NDOF1*NUMNODES_BOUNDARY;
     833        /*Initialize Load Vector and return if necessary*/
     834        Tria* tria=(Tria*)element;
     835        if(!tria->IsIceInElement()) return NULL;
    861836
    862837        /* Intermediaries*/
    863         int        i,ig,index1,index2;
    864         IssmDouble     DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN,thickness;
    865         IssmDouble     xyz_list[NUMVERTICES][3];
    866         IssmDouble     normal[2];
    867         IssmDouble     L[numdof];
    868         GaussTria *gauss;
    869 
    870         /*Initialize Load Vector and return if necessary*/
    871         ElementVector* pe = NULL;
    872         Tria*  tria=(Tria*)element;
    873         if(!tria->IsIceInElement()) return NULL;
     838        IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,thickness;
     839        IssmDouble xyz_list[NUMVERTICES][3];
     840        IssmDouble normal[2];
    874841
    875842        /*Retrieve all inputs and parameters*/
    876843        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    877         parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    878         Input* vxaverage_input   =tria->inputs->GetInput(VxEnum);                     _assert_(vxaverage_input);
    879         Input* vyaverage_input   =tria->inputs->GetInput(VyEnum);                     _assert_(vyaverage_input);
    880         Input* spcthickness_input=tria->inputs->GetInput(MasstransportSpcthicknessEnum); _assert_(spcthickness_input);
     844        IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
     845        Input* vxaverage_input    = tria->inputs->GetInput(VxEnum);                        _assert_(vxaverage_input);
     846        Input* vyaverage_input    = tria->inputs->GetInput(VyEnum);                        _assert_(vyaverage_input);
     847        Input* spcthickness_input = tria->inputs->GetInput(MasstransportSpcthicknessEnum); _assert_(spcthickness_input);
    881848        GetNormal(&normal[0],xyz_list);
    882849
    883850        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    884         index1=tria->GetNodeIndex(nodes[0]);
    885         index2=tria->GetNodeIndex(nodes[1]);
    886 
    887         gauss=new GaussTria();
     851        int index1=tria->GetVertexIndex(vertices[0]);
     852        int index2=tria->GetVertexIndex(vertices[1]);
     853        GaussTria* gauss=new GaussTria();
    888854        gauss->GaussEdgeCenter(index1,index2);
    889855        vxaverage_input->GetInputValue(&mean_vx,gauss);
    890856        vyaverage_input->GetInputValue(&mean_vy,gauss);
    891857        delete gauss;
    892 
    893         UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    894         if (UdotN>0){
     858        IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
     859        if(UdotN>0){
    895860                return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
    896861        }
    897         else{
    898                 pe=new ElementVector(nodes,NUMNODES_BOUNDARY,this->parameters);
    899         }
     862
     863        /*Initialize Load Vector */
     864        int            numnodes = this->GetNumberOfNodes();
     865        ElementVector *pe       = new ElementVector(nodes,numnodes,this->parameters);
     866        IssmDouble    *basis    = xNew<IssmDouble>(numnodes);
    900867
    901868        /* Start  looping on the number of gaussian points: */
    902869        gauss=new GaussTria(index1,index2,2);
    903         for(ig=gauss->begin();ig<gauss->end();ig++){
     870        for(int ig=gauss->begin();ig<gauss->end();ig++){
    904871
    905872                gauss->GaussPoint(ig);
    906873
    907                 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
     874                tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
    908875
    909876                vxaverage_input->GetInputValue(&vx,gauss);
    … …  
    916883                DL= - gauss->weight*Jdet*dt*UdotN*thickness;
    917884
    918                 for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
     885                for(int i=0;i<numnodes;i++) pe->values[i] += DL*basis[i];
    919886        }
    920887
    921888        /*Clean up and return*/
     889   xDelete<IssmDouble>(basis);
    922890        delete gauss;
    923891        return pe;
    … …  
    935903        /*Build unit outward pointing vector*/
    936904        IssmDouble vector[2];
    937         IssmDouble norm;
    938905
    939906        vector[0]=xyz_list[1][0] - xyz_list[0][0];
    940907        vector[1]=xyz_list[1][1] - xyz_list[0][1];
    941908
    942         norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
     909        IssmDouble norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
    943910
    944911        normal[0]= + vector[1]/norm;

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

    r23959 r24088  
    2121                int id;
    2222                int flux_type;
     23                int flux_degree;
    2324
    2425                /*Hooks*/

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

    r23640 r24088  
    5050        if(iomodel->meshelementtype==TriaEnum){
    5151                switch(finite_element){
     52                        case P0DGEnum:
     53                                numnodes = iomodel->numberofelements;
     54                                break;
    5255                        case P1Enum:
    5356                                numnodes = iomodel->numberofvertices;
    … …  
    173176                if(iomodel->meshelementtype==TriaEnum){
    174177                        switch(finite_element){
     178                                case P0DGEnum:
     179                                        element_numnodes=1;
     180                                        element_node_ids[0]=i;
     181                                        break;
    175182                                case P1Enum:
    176183                                        element_numnodes=3;

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r24082 r24088  
    11351135syn keyword cConstant OptionEnum
    11361136syn keyword cConstant P0ArrayEnum
     1137syn keyword cConstant P0DGEnum
    11371138syn keyword cConstant P1DGEnum
    11381139syn keyword cConstant P1P1Enum
    … …  
    12801281syn keyword cType Cfsurfacesquare
    12811282syn keyword cType Channel
    1282 syn keyword cType classes
    12831283syn keyword cType Constraint
    12841284syn keyword cType Constraints
    … …  
    12871287syn keyword cType ControlInput
    12881288syn keyword cType Covertree
     1289syn keyword cType DataSetParam
    12891290syn keyword cType DatasetInput
    1290 syn keyword cType DataSetParam
    12911291syn keyword cType Definition
    12921292syn keyword cType DependentObject
    … …  
    13011301syn keyword cType ElementHook
    13021302syn keyword cType ElementMatrix
     1303syn keyword cType ElementVector
    13031304syn keyword cType Elements
    1304 syn keyword cType ElementVector
    13051305syn keyword cType ExponentialVariogram
    13061306syn keyword cType ExternalResult
    … …  
    13091309syn keyword cType Friction
    13101310syn keyword cType Gauss
    1311 syn keyword cType GaussianVariogram
    1312 syn keyword cType gaussobjects
    13131311syn keyword cType GaussPenta
    13141312syn keyword cType GaussSeg
    13151313syn keyword cType GaussTetra
    13161314syn keyword cType GaussTria
     1315syn keyword cType GaussianVariogram
    13171316syn keyword cType GenericExternalResult
    13181317syn keyword cType GenericOption
    … …  
    13221321syn keyword cType Input
    13231322syn keyword cType Inputs
    1324 syn keyword cType IntArrayInput
    13251323syn keyword cType IntInput
    13261324syn keyword cType IntMatParam
    … …  
    13301328syn keyword cType IssmDirectApplicInterface
    13311329syn keyword cType IssmParallelDirectApplicInterface
    1332 syn keyword cType krigingobjects
    13331330syn keyword cType Load
    13341331syn keyword cType Loads
    … …  
    13411338syn keyword cType Matice
    13421339syn keyword cType Matlitho
    1343 syn keyword cType matrixobjects
    13441340syn keyword cType MatrixParam
    13451341syn keyword cType Misfit
    … …  
    13541350syn keyword cType Observations
    13551351syn keyword cType Option
     1352syn keyword cType OptionUtilities
    13561353syn keyword cType Options
    1357 syn keyword cType OptionUtilities
    13581354syn keyword cType Param
    13591355syn keyword cType Parameters
    … …  
    13681364syn keyword cType Regionaloutput
    13691365syn keyword cType Results
     1366syn keyword cType RiftStruct
    13701367syn keyword cType Riftfront
    1371 syn keyword cType RiftStruct
    13721368syn keyword cType Seg
    13731369syn keyword cType SegInput
     1370syn keyword cType SegRef
    13741371syn keyword cType Segment
    1375 syn keyword cType SegRef
    13761372syn keyword cType SpcDynamic
    13771373syn keyword cType SpcStatic
    … …  
    13931389syn keyword cType Vertex
    13941390syn keyword cType Vertices
     1391syn keyword cType classes
     1392syn keyword cType gaussobjects
     1393syn keyword cType krigingobjects
     1394syn keyword cType matrixobjects
    13951395syn keyword cType AdjointBalancethickness2Analysis
    13961396syn keyword cType AdjointBalancethicknessAnalysis
    … …  
    14111411syn keyword cType FreeSurfaceBaseAnalysis
    14121412syn keyword cType FreeSurfaceTopAnalysis
     1413syn keyword cType GLheightadvectionAnalysis
    14131414syn keyword cType GiaIvinsAnalysis
    1414 syn keyword cType GLheightadvectionAnalysis
    14151415syn keyword cType HydrologyDCEfficientAnalysis
    14161416syn keyword cType HydrologyDCInefficientAnalysis

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r24082 r24088  
    11331133        OptionEnum,
    11341134        P0ArrayEnum,
     1135        P0DGEnum,
    11351136        P1DGEnum,
    11361137        P1P1Enum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r24082 r24088  
    11371137                case OptionEnum : return "Option";
    11381138                case P0ArrayEnum : return "P0Array";
     1139                case P0DGEnum : return "P0DG";
    11391140                case P1DGEnum : return "P1DG";
    11401141                case P1P1Enum : return "P1P1";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r24082 r24088  
    11641164              else if (strcmp(name,"Option")==0) return OptionEnum;
    11651165              else if (strcmp(name,"P0Array")==0) return P0ArrayEnum;
     1166              else if (strcmp(name,"P0DG")==0) return P0DGEnum;
    11661167              else if (strcmp(name,"P1DG")==0) return P1DGEnum;
    11671168              else if (strcmp(name,"P1P1")==0) return P1P1Enum;
    … …  
    12431244              else if (strcmp(name,"StringParam")==0) return StringParamEnum;
    12441245              else if (strcmp(name,"SubelementFriction1")==0) return SubelementFriction1Enum;
    1245               else if (strcmp(name,"SubelementFriction2")==0) return SubelementFriction2Enum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum;
     1249              if (strcmp(name,"SubelementFriction2")==0) return SubelementFriction2Enum;
     1250              else if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum;
    12501251              else if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum;
    12511252              else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
Note: See TracChangeset for help on using the changeset viewer.