Changeset 24131


Ignore:
Timestamp:
08/30/19 11:58:34 (6 years ago)
Author:
Mathieu Morlighem
Message:

CHG: not needed anymore

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

Legend:

Unmodified
Added
Removed

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

    r24091 r24131  
    225225        }
    226226
    227 }
    228 /*}}}*/
    229 void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
    230         /*Compute B  matrix. B=[phi1 phi2 -phi3 -phi4]
    231          *
    232          * and phi1=phi3 phi2=phi4
    233          *
    234          * We assume B has been allocated already, of size: 1x4
    235          */
    236 
    237         /*Fetch number of nodes for this finite element*/
    238         int numnodes = this->NumberofNodes(finiteelement);
    239 
    240         /*Get nodal functions*/
    241         IssmDouble* basis=xNew<IssmDouble>(numnodes);
    242         GetNodalFunctions(basis,gauss,finiteelement);
    243 
    244         /*Build B for this segment*/
    245         switch(finiteelement){
    246                 case P0DGEnum:
    247                         B[0] = +basis[0];
    248                         B[1] = -basis[0];
    249                         break;
    250                 case P1DGEnum:
    251                         B[0] = +basis[index1];
    252                         B[1] = +basis[index2];
    253                         B[2] = -basis[index1];
    254                         B[3] = -basis[index2];
    255                         break;
    256                 default:
    257                         _error_("not supported yet");
    258         }
    259 
    260         /*Clean-up*/
    261         xDelete<IssmDouble>(basis);
    262 }
    263 /*}}}*/
    264 void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
    265         /*Compute Bprime  matrix. Bprime=[phi1 phi2 phi3 phi4]
    266          *
    267          * and phi1=phi3 phi2=phi4
    268          *
    269          * We assume Bprime has been allocated already, of size: 1x4
    270          */
    271 
    272         /*Fetch number of nodes for this finite element*/
    273         int numnodes = this->NumberofNodes(finiteelement);
    274 
    275         /*Get nodal functions*/
    276         IssmDouble* basis=xNew<IssmDouble>(numnodes);
    277         GetNodalFunctions(basis,gauss,finiteelement);
    278 
    279         /*Build B'*/
    280         /*Build B for this segment*/
    281         switch(finiteelement){
    282                 case P0DGEnum:
    283                         Bprime[0] = basis[0];
    284                         Bprime[1] = basis[0];
    285                         break;
    286                 case P1DGEnum:
    287                         Bprime[0] = basis[index1];
    288                         Bprime[1] = basis[index2];
    289                         Bprime[2] = basis[index1];
    290                         Bprime[3] = basis[index2];
    291                         break;
    292                 default:
    293                         _error_("not supported yet");
    294         }
    295 
    296         /*Clean-up*/
    297         xDelete<IssmDouble>(basis);
    298227}
    299228/*}}}*/

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

    r23962 r24131  
    2525                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
    2626                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss,int finiteelement);
    27                 void GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement);
    28                 void GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement);
    2927                void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss);
    3028                void GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss, int index1,int index2,int finiteelement);

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

    r24130 r24131  
    574574
    575575        /* Intermediaries*/
    576         IssmDouble DL1,DL2,Jdet,vx,vy,UdotN;
     576        IssmDouble A1,A2,Jdet,vx,vy,UdotN;
    577577        IssmDouble xyz_list[NUMVERTICES][3];
    578578        IssmDouble normal[2];
    579579
    580580        /*Fetch number of nodes for this flux*/
    581         int numnodes = this->GetNumberOfNodes();
     581        int numnodes       = this->GetNumberOfNodes();
     582        int numnodes_plus  = this->GetNumberOfNodesOneSide();
     583        int numnodes_minus = numnodes_plus; /*For now we are not doing p-adaptive DG*/
     584        _assert_(numnodes==numnodes_plus+numnodes_minus);
    582585
    583586        /*Initialize variables*/
    584         ElementMatrix *Ke     = new ElementMatrix(nodes,numnodes,this->parameters);
    585         IssmDouble    *B      = xNew<IssmDouble>(numnodes);
    586         IssmDouble    *Bprime = xNew<IssmDouble>(numnodes);
     587        ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters);
     588        IssmDouble    *basis_plus  = xNew<IssmDouble>(numnodes_plus);
     589        IssmDouble    *basis_minus = xNew<IssmDouble>(numnodes_minus);
    587590
    588591        /*Retrieve all inputs and parameters*/
    … …  
    600603                gauss->GaussPoint(ig);
    601604
    602                 tria->GetSegmentBFlux(&B[0],gauss,index1,index2,tria->FiniteElement());
    603                 tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2,tria->FiniteElement());
     605                tria->GetSegmentNodalFunctions(&basis_plus[0] ,gauss,index1,index2,tria->FiniteElement());
     606                tria->GetSegmentNodalFunctions(&basis_minus[0],gauss,index1,index2,tria->FiniteElement());
    604607
    605608                vxaverage_input->GetInputValue(&vx,gauss);
    … …  
    607610                UdotN=vx*normal[0]+vy*normal[1];
    608611                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
    609                 DL1=gauss->weight*Jdet*UdotN/2;
    610                 DL2=gauss->weight*Jdet*fabs(UdotN)/2;
    611 
    612                 for(int i=0;i<numnodes;i++){
    613                         for(int j=0;j<numnodes;j++){
    614                                 Ke->values[i*numnodes+j]+=DL1*B[i]*Bprime[j];
    615                                 Ke->values[i*numnodes+j]+=DL2*B[i]*B[j];
    616                         }
    617                 }
    618         }
    619 
    620         /*Clean up and return*/
    621         xDelete<IssmDouble>(B);
    622         xDelete<IssmDouble>(Bprime);
    623         delete gauss;
    624         return Ke;
    625 }
    626 /*}}}*/
    627 ElementMatrix* Numericalflux::CreateKMatrixMasstransport(void){/*{{{*/
    628 
    629         switch(this->flux_type){
    630                 case InternalEnum:
    631                         return CreateKMatrixMasstransportInternal();
    632                 case BoundaryEnum:
    633                         return CreateKMatrixMasstransportBoundary();
    634                 default:
    635                         _error_("type not supported yet");
    636         }
    637 }
    638 /*}}}*/
    639 ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/
    640 
    641         /*Initialize Element matrix and return if necessary*/
    642         Tria*  tria=(Tria*)element;
    643         if(!tria->IsIceInElement()) return NULL;
    644 
    645         /* Intermediaries*/
    646         IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy;
    647         IssmDouble xyz_list[NUMVERTICES][3];
    648         IssmDouble normal[2];
    649 
    650         /*Retrieve all inputs and parameters*/
    651         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    652         IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
    653         Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
    654         Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
    655         GetNormal(&normal[0],xyz_list);
    656 
    657         /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    658         int index1=tria->GetVertexIndex(vertices[0]);
    659         int index2=tria->GetVertexIndex(vertices[1]);
    660 
    661         GaussTria* gauss=new GaussTria();
    662         gauss->GaussEdgeCenter(index1,index2);
    663         vxaverage_input->GetInputValue(&mean_vx,gauss);
    664         vyaverage_input->GetInputValue(&mean_vy,gauss);
    665         delete gauss;
    666 
    667         IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    668         if(UdotN<=0){
    669                 return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
    670         }
    671 
    672         /*Initialize Element vector and other vectors*/
    673    int            numnodes = this->GetNumberOfNodes();
    674    ElementMatrix *Ke       = new ElementMatrix(nodes,numnodes,this->parameters);
    675    IssmDouble    *basis    = xNew<IssmDouble>(numnodes);
    676 
    677         /* Start  looping on the number of gaussian points: */
    678         gauss=new GaussTria(index1,index2,2);
    679         for(int ig=gauss->begin();ig<gauss->end();ig++){
    680 
    681                 gauss->GaussPoint(ig);
    682 
    683                 tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
    684 
    685                 vxaverage_input->GetInputValue(&vx,gauss);
    686                 vyaverage_input->GetInputValue(&vy,gauss);
    687                 UdotN=vx*normal[0]+vy*normal[1];
    688                 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
    689                 DL=gauss->weight*Jdet*dt*UdotN;
    690 
    691                 for(int i=0;i<numnodes;i++){
    692                         for(int j=0;j<numnodes;j++){
    693                                 Ke->values[i*numnodes+j]+=DL*basis[i]*basis[j];
    694                         }
    695                 }
    696         }
    697 
    698         /*Clean up and return*/
    699         xDelete<IssmDouble>(basis);
    700         delete gauss;
    701         return Ke;
    702 }
    703 /*}}}*/
    704 ElementMatrix* Numericalflux::CreateKMatrixMasstransportInternal(void){/*{{{*/
    705 
    706         /*Initialize Element matrix and return if necessary*/
    707         Tria*  tria=(Tria*)element;
    708         if(!tria->IsIceInElement()) return NULL;
    709 
    710         /* Intermediaries*/
    711         IssmDouble A1,A2,Jdet,vx,vy,UdotN;
    712         IssmDouble xyz_list[NUMVERTICES][3];
    713         IssmDouble normal[2];
    714 
    715         /*Fetch number of nodes for this flux*/
    716         int numnodes       = this->GetNumberOfNodes();
    717         int numnodes_plus  = this->GetNumberOfNodesOneSide();
    718         int numnodes_minus = numnodes_plus; /*For now we are not doing p-adaptive DG*/
    719         _assert_(numnodes==numnodes_plus+numnodes_minus);
    720 
    721         /*Initialize variables*/
    722         ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters);
    723         IssmDouble    *basis_plus  = xNew<IssmDouble>(numnodes_plus);
    724         IssmDouble    *basis_minus = xNew<IssmDouble>(numnodes_minus);
    725 
    726         /*Retrieve all inputs and parameters*/
    727         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    728         IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
    729         Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
    730         Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
    731         GetNormal(&normal[0],xyz_list);
    732 
    733         /* Start  looping on the number of gaussian points: */
    734         int index1=tria->GetVertexIndex(vertices[0]);
    735         int index2=tria->GetVertexIndex(vertices[1]);
    736         GaussTria* gauss=new GaussTria(index1,index2,2);
    737         for(int ig=gauss->begin();ig<gauss->end();ig++){
    738 
    739                 gauss->GaussPoint(ig);
    740 
    741                 tria->GetSegmentNodalFunctions(&basis_plus[0] ,gauss,index1,index2,tria->FiniteElement());
    742                 tria->GetSegmentNodalFunctions(&basis_minus[0],gauss,index1,index2,tria->FiniteElement());
    743 
    744                 vxaverage_input->GetInputValue(&vx,gauss);
    745                 vyaverage_input->GetInputValue(&vy,gauss);
    746                 UdotN=vx*normal[0]+vy*normal[1];
    747                 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
    748                 A1=gauss->weight*Jdet*dt*UdotN/2;
    749                 A2=gauss->weight*Jdet*dt*fabs(UdotN)/2;
     612                A1=gauss->weight*Jdet*UdotN/2;
     613                A2=gauss->weight*Jdet*fabs(UdotN)/2;
    750614
    751615                /*Term 1 (numerical flux): {Hv}.[[phi]] = 0.5(H+v+ + H-v-)(phi+n+ + phi-n-)
    … …  
    798662}
    799663/*}}}*/
     664ElementMatrix* Numericalflux::CreateKMatrixMasstransport(void){/*{{{*/
     665
     666        switch(this->flux_type){
     667                case InternalEnum:
     668                        return CreateKMatrixMasstransportInternal();
     669                case BoundaryEnum:
     670                        return CreateKMatrixMasstransportBoundary();
     671                default:
     672                        _error_("type not supported yet");
     673        }
     674}
     675/*}}}*/
     676ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/
     677
     678        /*Initialize Element matrix and return if necessary*/
     679        Tria*  tria=(Tria*)element;
     680        if(!tria->IsIceInElement()) return NULL;
     681
     682        /* Intermediaries*/
     683        IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy;
     684        IssmDouble xyz_list[NUMVERTICES][3];
     685        IssmDouble normal[2];
     686
     687        /*Retrieve all inputs and parameters*/
     688        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     689        IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
     690        Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
     691        Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
     692        GetNormal(&normal[0],xyz_list);
     693
     694        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
     695        int index1=tria->GetVertexIndex(vertices[0]);
     696        int index2=tria->GetVertexIndex(vertices[1]);
     697
     698        GaussTria* gauss=new GaussTria();
     699        gauss->GaussEdgeCenter(index1,index2);
     700        vxaverage_input->GetInputValue(&mean_vx,gauss);
     701        vyaverage_input->GetInputValue(&mean_vy,gauss);
     702        delete gauss;
     703
     704        IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
     705        if(UdotN<=0){
     706                return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
     707        }
     708
     709        /*Initialize Element vector and other vectors*/
     710   int            numnodes = this->GetNumberOfNodes();
     711   ElementMatrix *Ke       = new ElementMatrix(nodes,numnodes,this->parameters);
     712   IssmDouble    *basis    = xNew<IssmDouble>(numnodes);
     713
     714        /* Start  looping on the number of gaussian points: */
     715        gauss=new GaussTria(index1,index2,2);
     716        for(int ig=gauss->begin();ig<gauss->end();ig++){
     717
     718                gauss->GaussPoint(ig);
     719
     720                tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
     721
     722                vxaverage_input->GetInputValue(&vx,gauss);
     723                vyaverage_input->GetInputValue(&vy,gauss);
     724                UdotN=vx*normal[0]+vy*normal[1];
     725                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
     726                DL=gauss->weight*Jdet*dt*UdotN;
     727
     728                for(int i=0;i<numnodes;i++){
     729                        for(int j=0;j<numnodes;j++){
     730                                Ke->values[i*numnodes+j]+=DL*basis[i]*basis[j];
     731                        }
     732                }
     733        }
     734
     735        /*Clean up and return*/
     736        xDelete<IssmDouble>(basis);
     737        delete gauss;
     738        return Ke;
     739}
     740/*}}}*/
     741ElementMatrix* Numericalflux::CreateKMatrixMasstransportInternal(void){/*{{{*/
     742
     743        /*Initialize Element matrix and return if necessary*/
     744        Tria*  tria=(Tria*)element;
     745        if(!tria->IsIceInElement()) return NULL;
     746
     747        /* Intermediaries*/
     748        IssmDouble A1,A2,Jdet,vx,vy,UdotN;
     749        IssmDouble xyz_list[NUMVERTICES][3];
     750        IssmDouble normal[2];
     751
     752        /*Fetch number of nodes for this flux*/
     753        int numnodes       = this->GetNumberOfNodes();
     754        int numnodes_plus  = this->GetNumberOfNodesOneSide();
     755        int numnodes_minus = numnodes_plus; /*For now we are not doing p-adaptive DG*/
     756        _assert_(numnodes==numnodes_plus+numnodes_minus);
     757
     758        /*Initialize variables*/
     759        ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters);
     760        IssmDouble    *basis_plus  = xNew<IssmDouble>(numnodes_plus);
     761        IssmDouble    *basis_minus = xNew<IssmDouble>(numnodes_minus);
     762
     763        /*Retrieve all inputs and parameters*/
     764        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     765        IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
     766        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
     767        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
     768        GetNormal(&normal[0],xyz_list);
     769
     770        /* Start  looping on the number of gaussian points: */
     771        int index1=tria->GetVertexIndex(vertices[0]);
     772        int index2=tria->GetVertexIndex(vertices[1]);
     773        GaussTria* gauss=new GaussTria(index1,index2,2);
     774        for(int ig=gauss->begin();ig<gauss->end();ig++){
     775
     776                gauss->GaussPoint(ig);
     777
     778                tria->GetSegmentNodalFunctions(&basis_plus[0] ,gauss,index1,index2,tria->FiniteElement());
     779                tria->GetSegmentNodalFunctions(&basis_minus[0],gauss,index1,index2,tria->FiniteElement());
     780
     781                vxaverage_input->GetInputValue(&vx,gauss);
     782                vyaverage_input->GetInputValue(&vy,gauss);
     783                UdotN=vx*normal[0]+vy*normal[1];
     784                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
     785                A1=gauss->weight*Jdet*dt*UdotN/2;
     786                A2=gauss->weight*Jdet*dt*fabs(UdotN)/2;
     787
     788                /*Term 1 (numerical flux): {Hv}.[[phi]] = 0.5(H+v+ + H-v-)(phi+n+ + phi-n-)
     789                 *                                      = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
     790                 *                                      = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
     791                 *
     792                 *Term 2 (stabilization)  |v.n|/2 [[H]].[[phi]] = |v.n|/2 (H+n+ + H-n-)(phi+n+ + phi-n-)
     793                 *                                      = |v.n|/2 (H+phi+ -H-phi+ -H+phi- +H-phi-)
     794                 *     | A++ | A+- |
     795                 * K = |-----------|
     796                 *     | A-+ | A-- |
     797                 *
     798                 *These 4 terms for each expressions are added independently*/
     799
     800                /*First term A++*/
     801                for(int i=0;i<numnodes_plus;i++){
     802                        for(int j=0;j<numnodes_plus;j++){
     803                                Ke->values[i*numnodes+j] += A1*(basis_plus[j]*basis_plus[i]);
     804                                Ke->values[i*numnodes+j] += A2*(basis_plus[j]*basis_plus[i]);
     805                        }
     806                }
     807                /*Second term A+-*/
     808                for(int i=0;i<numnodes_plus;i++){
     809                        for(int j=0;j<numnodes_minus;j++){
     810                                Ke->values[i*numnodes+numnodes_plus+j] +=  A1*(basis_minus[j]*basis_plus[i]);
     811                                Ke->values[i*numnodes+numnodes_plus+j] += -A2*(basis_minus[j]*basis_plus[i]);
     812                        }
     813                }
     814                /*Third term A-+*/
     815                for(int i=0;i<numnodes_minus;i++){
     816                        for(int j=0;j<numnodes_plus;j++){
     817                                Ke->values[(numnodes_plus+i)*numnodes+j] += -A1*(basis_plus[j]*basis_minus[i]);
     818                                Ke->values[(numnodes_plus+i)*numnodes+j] += -A2*(basis_plus[j]*basis_minus[i]);
     819                        }
     820                }
     821                /*Fourth term A-+*/
     822                for(int i=0;i<numnodes_minus;i++){
     823                        for(int j=0;j<numnodes_minus;j++){
     824                                Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] += -A1*(basis_minus[j]*basis_minus[i]);
     825                                Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] +=  A2*(basis_minus[j]*basis_minus[i]);
     826                        }
     827                }
     828        }
     829
     830        /*Clean up and return*/
     831   xDelete<IssmDouble>(basis_plus);
     832   xDelete<IssmDouble>(basis_minus);
     833        delete gauss;
     834        return Ke;
     835}
     836/*}}}*/
    800837ElementVector* Numericalflux::CreatePVectorAdjointBalancethickness(void){/*{{{*/
    801838
Note: See TracChangeset for help on using the changeset viewer.