Changeset 6026


Ignore:
Timestamp:
09/24/10 11:37:17 (15 years ago)
Author:
Mathieu Morlighem
Message:

split CreateKMatrix and PVector depending on analysis: it is too dangerous to keep only one method

Location:
issm/trunk/src/c/objects/Loads
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified issm/trunk/src/c/objects/Loads/Numericalflux.cpp

    r5989 r6026  
    330330void  Numericalflux::CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs){
    331331
    332         int type;
     332        /*recover some parameters*/
     333        ElementMatrix* Ke=NULL;
    333334        int analysis_type;
    334         ElementMatrix* Ke=NULL;
    335 
    336         /*recover some parameters*/
    337         inputs->GetParameterValue(&type,TypeEnum);
    338335        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    339336
    340337        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    341338        switch(analysis_type){
    342                 case PrognosticAnalysisEnum: case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum:
    343                         switch(type){
    344                                 case InternalEnum:
    345                                         Ke=CreateKMatrixInternal();
    346                                         break;
    347                                 case BoundaryEnum:
    348                                         Ke=CreateKMatrixBoundary();
    349                                         break;
    350                                 default:
    351                                         ISSMERROR("type not supported yet");
    352                         }
     339                case PrognosticAnalysisEnum:
     340                        Ke=CreateKMatrixPrognostic();
     341                        break;
     342                case BalancedthicknessAnalysisEnum:
     343                        Ke=CreateKMatrixBalancedthickness();
     344                        break;
     345                case AdjointBalancedthicknessAnalysisEnum:
     346                        Ke=CreateKMatrixAdjointBalancedthickness();
    353347                        break;
    354348                default:
    … …  
    367361void  Numericalflux::CreatePVector(Vec pg,Vec pf){
    368362
    369         int type;
     363        /*recover some parameters*/
     364        ElementVector* pe=NULL;
    370365        int analysis_type;
    371         ElementVector* pe=NULL;
    372 
    373         /*recover some parameters*/
    374         inputs->GetParameterValue(&type,TypeEnum);
    375366        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    376367
    377368        switch(analysis_type){
    378                 case PrognosticAnalysisEnum: case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum:
    379                         switch(type){
    380                                 case InternalEnum:
    381                                         pe=CreatePVectorInternal();
    382                                         break;
    383                                 case BoundaryEnum:
    384                                         pe=CreatePVectorBoundary();
    385                                         break;
    386                                 default:
    387                                         ISSMERROR("type not supported yet");
    388                         }
     369                case PrognosticAnalysisEnum:
     370                        pe=CreatePVectorPrognostic();
     371                        break;
     372                case BalancedthicknessAnalysisEnum:
     373                        pe=CreatePVectorBalancedthickness();
     374                        break;
     375                case AdjointBalancedthicknessAnalysisEnum:
     376                        pe=CreatePVectorAdjointBalancedthickness();
    389377                        break;
    390378                default:
    … …  
    424412
    425413/*Numericalflux management*/
    426 /*FUNCTION Numericalflux::CreateKMatrixInternal {{{1*/
    427 ElementMatrix* Numericalflux::CreateKMatrixInternal(void){
     414/*FUNCTION Numericalflux::CreateKMatrixPrognostic{{{1*/
     415ElementMatrix* Numericalflux::CreateKMatrixPrognostic(void){
     416
     417        int type;
     418        inputs->GetParameterValue(&type,TypeEnum);
     419
     420        switch(type){
     421                case InternalEnum:
     422                        return CreateKMatrixPrognosticInternal();
     423                case BoundaryEnum:
     424                        return CreateKMatrixPrognosticBoundary();
     425                default:
     426                        ISSMERROR("type not supported yet");
     427        }
     428}
     429/*}}}*/
     430/*FUNCTION Numericalflux::CreateKMatrixPrognosticInternal {{{1*/
     431ElementMatrix* Numericalflux::CreateKMatrixPrognosticInternal(void){
    428432
    429433        /* constants*/
    … …  
    431435
    432436        /* Intermediaries*/
    433         int        i,j,ig,index1,index2,analysis_type;
     437        int        i,j,ig,index1,index2;
    434438        double     DL1,DL2,Jdet,dt,vx,vy,UdotN;
    435439        double     xyz_list[NUMVERTICES_INTERNAL][3];
    … …  
    448452        /*Retrieve all inputs and parameters*/
    449453        GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES_INTERNAL);
     454        parameters->FindParam(&dt,DtEnum);
    450455        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
    451456        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
    452         this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    453457        GetNormal(&normal[0],xyz_list);
    454         switch(analysis_type){
    455                 case PrognosticAnalysisEnum:
    456                         parameters->FindParam(&dt,DtEnum); break;
    457                 case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum:
    458                         dt=1; break;/*No transient term is involved*/
    459                 default:
    460                         ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type));
    461         }
    462458
    463459        /* Start  looping on the number of gaussian points: */
    … …  
    497493}
    498494/*}}}*/
    499 /*FUNCTION Numericalflux::CreateKMatrixBoundary {{{1*/
    500 ElementMatrix* Numericalflux::CreateKMatrixBoundary(void){
     495/*FUNCTION Numericalflux::CreateKMatrixPrognosticBoundary {{{1*/
     496ElementMatrix* Numericalflux::CreateKMatrixPrognosticBoundary(void){
    501497
    502498        /* constants*/
    … …  
    504500
    505501        /* Intermediaries*/
    506         int        i,j,ig,index1,index2,analysis_type;
     502        int        i,j,ig,index1,index2;
    507503        double     DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN;
    508504        double     xyz_list[NUMVERTICES_BOUNDARY][3];
    … …  
    519515        /*Retrieve all inputs and parameters*/
    520516        GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES_BOUNDARY);
     517        parameters->FindParam(&dt,DtEnum);
    521518        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
    522519        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
    523         this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    524520        GetNormal(&normal[0],xyz_list);
    525         switch(analysis_type){
    526                 case PrognosticAnalysisEnum:
    527                         parameters->FindParam(&dt,DtEnum); break;
    528                 case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum:
    529                         dt=1; break;/*No transient term is involved*/
    530                 default:
    531                         ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type));
    532         }
    533521
    534522        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    … …  
    575563}
    576564/*}}}*/
    577 /*FUNCTION Numericalflux::CreatePVectorInternal{{{1*/
    578 ElementVector* Numericalflux::CreatePVectorInternal(void){
     565/*FUNCTION Numericalflux::CreateKMatrixBalancedthickness{{{1*/
     566ElementMatrix* Numericalflux::CreateKMatrixBalancedthickness(void){
     567
     568        int type;
     569        inputs->GetParameterValue(&type,TypeEnum);
     570
     571        switch(type){
     572                case InternalEnum:
     573                        return CreateKMatrixBalancedthicknessInternal();
     574                case BoundaryEnum:
     575                        return CreateKMatrixBalancedthicknessBoundary();
     576                default:
     577                        ISSMERROR("type not supported yet");
     578        }
     579}
     580/*}}}*/
     581/*FUNCTION Numericalflux::CreateKMatrixBalancedthicknessInternal {{{1*/
     582ElementMatrix* Numericalflux::CreateKMatrixBalancedthicknessInternal(void){
     583
     584        /* constants*/
     585        const int numdof=NDOF1*NUMVERTICES_INTERNAL;
     586
     587        /* Intermediaries*/
     588        int        i,j,ig,index1,index2;
     589        double     DL1,DL2,Jdet,vx,vy,UdotN;
     590        double     xyz_list[NUMVERTICES_INTERNAL][3];
     591        double     normal[2];
     592        double     B[numdof];
     593        double     Bprime[numdof];
     594        double     Ke_g1[numdof][numdof];
     595        double     Ke_g2[numdof][numdof];
     596        GaussTria *gauss;
     597
     598        /*Initialize Element matrix and return if necessary*/
     599        Tria*  tria=(Tria*)element;
     600        if(tria->IsOnWater()) return NULL;
     601        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES_INTERNAL,this->parameters);
     602
     603        /*Retrieve all inputs and parameters*/
     604        GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES_INTERNAL);
     605        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
     606        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
     607        GetNormal(&normal[0],xyz_list);
     608
     609        /* Start  looping on the number of gaussian points: */
     610        index1=tria->GetNodeIndex(nodes[0]);
     611        index2=tria->GetNodeIndex(nodes[1]);
     612        gauss=new GaussTria(index1,index2,2);
     613        for(ig=gauss->begin();ig<gauss->end();ig++){
     614
     615                gauss->GaussPoint(ig);
     616
     617                tria->GetSegmentBFlux(&B[0],gauss,index1,index2);
     618                tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2);
     619
     620                vxaverage_input->GetParameterValue(&vx,gauss);
     621                vyaverage_input->GetParameterValue(&vy,gauss);
     622                UdotN=vx*normal[0]+vy*normal[1];
     623                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
     624                DL1=gauss->weight*Jdet*UdotN/2;
     625                DL2=gauss->weight*Jdet*fabs(UdotN)/2;
     626
     627                TripleMultiply(&B[0],1,numdof,1,
     628                                        &DL1,1,1,0,
     629                                        &Bprime[0],1,numdof,0,
     630                                        &Ke_g1[0][0],0);
     631                TripleMultiply(&B[0],1,numdof,1,
     632                                        &DL2,1,1,0,
     633                                        &B[0],1,numdof,0,
     634                                        &Ke_g2[0][0],0);
     635
     636                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g1[i][j];
     637                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g2[i][j];
     638        }
     639
     640        /*Clean up and return*/
     641        delete gauss;
     642        return Ke;
     643}
     644/*}}}*/
     645/*FUNCTION Numericalflux::CreateKMatrixBalancedthicknessBoundary {{{1*/
     646ElementMatrix* Numericalflux::CreateKMatrixBalancedthicknessBoundary(void){
     647
     648        /* constants*/
     649        const int numdof=NDOF1*NUMVERTICES_BOUNDARY;
     650
     651        /* Intermediaries*/
     652        int        i,j,ig,index1,index2;
     653        double     DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN;
     654        double     xyz_list[NUMVERTICES_BOUNDARY][3];
     655        double     normal[2];
     656        double     L[numdof];
     657        double     Ke_g[numdof][numdof];
     658        GaussTria *gauss;
     659
     660        /*Initialize Element matrix and return if necessary*/
     661        Tria*  tria=(Tria*)element;
     662        if(tria->IsOnWater()) return NULL;
     663        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters);
     664
     665        /*Retrieve all inputs and parameters*/
     666        GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES_BOUNDARY);
     667        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
     668        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
     669        GetNormal(&normal[0],xyz_list);
     670
     671        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
     672        index1=tria->GetNodeIndex(nodes[0]);
     673        index2=tria->GetNodeIndex(nodes[1]);
     674
     675        gauss=new GaussTria();
     676        gauss->GaussEdgeCenter(index1,index2);
     677        vxaverage_input->GetParameterValue(&mean_vx,gauss);
     678        vyaverage_input->GetParameterValue(&mean_vy,gauss);
     679        delete gauss;
     680
     681        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
     682        if (UdotN<=0){
     683                /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
     684                return NULL;
     685        }
     686
     687        /* Start  looping on the number of gaussian points: */
     688        gauss=new GaussTria(index1,index2,2);
     689        for(ig=gauss->begin();ig<gauss->end();ig++){
     690
     691                gauss->GaussPoint(ig);
     692
     693                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
     694
     695                vxaverage_input->GetParameterValue(&vx,gauss);
     696                vyaverage_input->GetParameterValue(&vy,gauss);
     697                UdotN=vx*normal[0]+vy*normal[1];
     698                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
     699                DL=gauss->weight*Jdet*UdotN;
     700
     701                TripleMultiply(&L[0],1,numdof,1,
     702                                        &DL,1,1,0,
     703                                        &L[0],1,numdof,0,
     704                                        &Ke_g[0][0],0);
     705
     706                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
     707        }
     708
     709        /*Clean up and return*/
     710        delete gauss;
     711        return Ke;
     712}
     713/*}}}*/
     714/*FUNCTION Numericalflux::CreateKMatrixAdjointBalancedthickness{{{1*/
     715ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancedthickness(void){
     716
     717        int type;
     718        inputs->GetParameterValue(&type,TypeEnum);
     719
     720        switch(type){
     721                case InternalEnum:
     722                        return CreateKMatrixAdjointBalancedthicknessInternal();
     723                case BoundaryEnum:
     724                        return CreateKMatrixAdjointBalancedthicknessBoundary();
     725                default:
     726                        ISSMERROR("type not supported yet");
     727        }
     728}
     729/*}}}*/
     730/*FUNCTION Numericalflux::CreateKMatrixAdjointBalancedthicknessInternal {{{1*/
     731ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancedthicknessInternal(void){
     732
     733        return CreateKMatrixBalancedthicknessInternal();
     734}
     735/*}}}*/
     736/*FUNCTION Numericalflux::CreateKMatrixAdjointBalancedthicknessBoundary {{{1*/
     737ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancedthicknessBoundary(void){
     738
     739        return CreateKMatrixBalancedthicknessBoundary();
     740}
     741/*}}}*/
     742/*FUNCTION Numericalflux::CreatePVectorPrognostic{{{1*/
     743ElementVector* Numericalflux::CreatePVectorPrognostic(void){
     744
     745        int type;
     746        inputs->GetParameterValue(&type,TypeEnum);
     747
     748        switch(type){
     749                case InternalEnum:
     750                        return CreatePVectorPrognosticInternal();
     751                case BoundaryEnum:
     752                        return CreatePVectorPrognosticBoundary();
     753                default:
     754                        ISSMERROR("type not supported yet");
     755        }
     756}
     757/*}}}*/
     758/*FUNCTION Numericalflux::CreatePVectorPrognosticInternal{{{1*/
     759ElementVector* Numericalflux::CreatePVectorPrognosticInternal(void){
    579760
    580761        /*Nothing added to PVector*/
    … …  
    583764}
    584765/*}}}*/
    585 /*FUNCTION Numericalflux::CreatePVectorBoundary{{{1*/
    586 ElementVector* Numericalflux::CreatePVectorBoundary(void){
     766/*FUNCTION Numericalflux::CreatePVectorPrognosticBoundary{{{1*/
     767ElementVector* Numericalflux::CreatePVectorPrognosticBoundary(void){
    587768
    588769        /* constants*/
    … …  
    590771
    591772        /* Intermediaries*/
    592         int        i,j,ig,index1,index2,analysis_type;
     773        int        i,j,ig,index1,index2;
    593774        double     DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN,thickness;
    594775        double     xyz_list[NUMVERTICES_BOUNDARY][3];
    … …  
    597778        GaussTria *gauss;
    598779
    599         /*Initialize Load Vectorand return if necessary*/
     780        /*Initialize Load Vector and return if necessary*/
    600781        Tria*  tria=(Tria*)element;
    601782        if(tria->IsOnWater()) return NULL;
    … …  
    604785        /*Retrieve all inputs and parameters*/
    605786        GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES_BOUNDARY);
     787        parameters->FindParam(&dt,DtEnum);
    606788        Input* vxaverage_input=tria->inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input);
    607789        Input* vyaverage_input=tria->inputs->GetInput(VyEnum); ISSMASSERT(vyaverage_input);
    608         Input* thickness_input=tria->inputs->GetInput(ThicknessObsEnum);
    609 
    610         /*Here, as it is a forcing, we have H=Hobs by default (for control methods)*/
    611         if (!thickness_input){
    612                 thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
    613         }
    614 
    615         this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     790        Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
    616791        GetNormal(&normal[0],xyz_list);
    617         switch(analysis_type){
    618                 case PrognosticAnalysisEnum:
    619                         parameters->FindParam(&dt,DtEnum); break;
    620                 case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum:
    621                         dt=1; break;/*No transient term is involved*/
    622                 default:
    623                         ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type));
    624         }
    625792
    626793        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    … …  
    662829}
    663830/*}}}*/
     831/*FUNCTION Numericalflux::CreatePVectorBalancedthickness{{{1*/
     832ElementVector* Numericalflux::CreatePVectorBalancedthickness(void){
     833
     834        int type;
     835        inputs->GetParameterValue(&type,TypeEnum);
     836
     837        switch(type){
     838                case InternalEnum:
     839                        return CreatePVectorBalancedthicknessInternal();
     840                case BoundaryEnum:
     841                        return CreatePVectorBalancedthicknessBoundary();
     842                default:
     843                        ISSMERROR("type not supported yet");
     844        }
     845}
     846/*}}}*/
     847/*FUNCTION Numericalflux::CreatePVectorBalancedthicknessInternal{{{1*/
     848ElementVector* Numericalflux::CreatePVectorBalancedthicknessInternal(void){
     849
     850        /*Nothing added to PVector*/
     851        return NULL;
     852
     853}
     854/*}}}*/
     855/*FUNCTION Numericalflux::CreatePVectorBalancedthicknessBoundary{{{1*/
     856ElementVector* Numericalflux::CreatePVectorBalancedthicknessBoundary(void){
     857
     858        /* constants*/
     859        const int numdof=NDOF1*NUMVERTICES_BOUNDARY;
     860
     861        /* Intermediaries*/
     862        int        i,j,ig,index1,index2;
     863        double     DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN,thickness;
     864        double     xyz_list[NUMVERTICES_BOUNDARY][3];
     865        double     normal[2];
     866        double     L[numdof];
     867        GaussTria *gauss;
     868
     869        /*Initialize Load Vector and return if necessary*/
     870        Tria*  tria=(Tria*)element;
     871        if(tria->IsOnWater()) return NULL;
     872        ElementVector* pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);
     873
     874        /*Retrieve all inputs and parameters*/
     875        GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES_BOUNDARY);
     876        Input* vxaverage_input=tria->inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input);
     877        Input* vyaverage_input=tria->inputs->GetInput(VyEnum); ISSMASSERT(vyaverage_input);
     878        Input* thickness_input=tria->inputs->GetInput(ThicknessObsEnum);
     879
     880        /*Here, as it is a forcing, we have H=Hobs by default (for control methods)*/
     881        if (!thickness_input){
     882                thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
     883        }
     884        GetNormal(&normal[0],xyz_list);
     885
     886        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
     887        index1=tria->GetNodeIndex(nodes[0]);
     888        index2=tria->GetNodeIndex(nodes[1]);
     889
     890        gauss=new GaussTria();
     891        gauss->GaussEdgeCenter(index1,index2);
     892        vxaverage_input->GetParameterValue(&mean_vx,gauss);
     893        vyaverage_input->GetParameterValue(&mean_vy,gauss);
     894        delete gauss;
     895        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
     896        if (UdotN>0){
     897                /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
     898                return NULL;
     899        }
     900
     901        /* Start  looping on the number of gaussian points: */
     902        gauss=new GaussTria(index1,index2,2);
     903        for(ig=gauss->begin();ig<gauss->end();ig++){
     904
     905                gauss->GaussPoint(ig);
     906
     907                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
     908
     909                vxaverage_input->GetParameterValue(&vx,gauss);
     910                vyaverage_input->GetParameterValue(&vy,gauss);
     911                thickness_input->GetParameterValue(&thickness,gauss);
     912                UdotN=vx*normal[0]+vy*normal[1];
     913                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
     914                DL= - gauss->weight*Jdet*UdotN*thickness;
     915
     916                for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
     917        }
     918
     919        /*Clean up and return*/
     920        delete gauss;
     921        return pe;
     922}
     923/*}}}*/
     924/*FUNCTION Numericalflux::CreatePVectorAdjointBalancedthickness{{{1*/
     925ElementVector* Numericalflux::CreatePVectorAdjointBalancedthickness(void){
     926
     927        int type;
     928        inputs->GetParameterValue(&type,TypeEnum);
     929
     930        switch(type){
     931                case InternalEnum:
     932                        return CreatePVectorAdjointBalancedthicknessInternal();
     933                case BoundaryEnum:
     934                        return CreatePVectorAdjointBalancedthicknessBoundary();
     935                default:
     936                        ISSMERROR("type not supported yet");
     937        }
     938}
     939/*}}}*/
     940/*FUNCTION Numericalflux::CreatePVectorAdjointBalancedthicknessInternal{{{1*/
     941ElementVector* Numericalflux::CreatePVectorAdjointBalancedthicknessInternal(void){
     942
     943        /*Nothing added to PVector*/
     944        return NULL;
     945
     946}
     947/*}}}*/
     948/*FUNCTION Numericalflux::CreatePVectorAdjointBalancedthicknessBoundary{{{1*/
     949ElementVector* Numericalflux::CreatePVectorAdjointBalancedthicknessBoundary(void){
     950
     951        /* constants*/
     952        const int numdof=NDOF1*NUMVERTICES_BOUNDARY;
     953
     954        /* Intermediaries*/
     955        int        i,j,ig,index1,index2;
     956        double     DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN,thickness;
     957        double     xyz_list[NUMVERTICES_BOUNDARY][3];
     958        double     normal[2];
     959        double     L[numdof];
     960        GaussTria *gauss;
     961
     962        /*Initialize Load Vector and return if necessary*/
     963        Tria*  tria=(Tria*)element;
     964        if(tria->IsOnWater()) return NULL;
     965        ElementVector* pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);
     966
     967        /*Retrieve all inputs and parameters*/
     968        GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES_BOUNDARY);
     969        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);           ISSMASSERT(vxaverage_input);
     970        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);           ISSMASSERT(vyaverage_input);
     971        Input* thickness_input=tria->inputs->GetInput(ThicknessObsEnum); ISSMASSERT(thickness_input);
     972        GetNormal(&normal[0],xyz_list);
     973
     974        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
     975        index1=tria->GetNodeIndex(nodes[0]);
     976        index2=tria->GetNodeIndex(nodes[1]);
     977
     978        gauss=new GaussTria();
     979        gauss->GaussEdgeCenter(index1,index2);
     980        vxaverage_input->GetParameterValue(&mean_vx,gauss);
     981        vyaverage_input->GetParameterValue(&mean_vy,gauss);
     982        delete gauss;
     983        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
     984        if (UdotN>0){
     985                /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
     986                return NULL;
     987        }
     988
     989        /* Start  looping on the number of gaussian points: */
     990        gauss=new GaussTria(index1,index2,2);
     991        for(ig=gauss->begin();ig<gauss->end();ig++){
     992
     993                gauss->GaussPoint(ig);
     994
     995                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
     996
     997                vxaverage_input->GetParameterValue(&vx,gauss);
     998                vyaverage_input->GetParameterValue(&vy,gauss);
     999                thickness_input->GetParameterValue(&thickness,gauss);
     1000                UdotN=vx*normal[0]+vy*normal[1];
     1001                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
     1002                DL= - gauss->weight*Jdet*UdotN*thickness;
     1003
     1004                for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
     1005        }
     1006
     1007        /*Clean up and return*/
     1008        delete gauss;
     1009        return pe;
     1010}
     1011/*}}}*/
    6641012/*FUNCTION Numericalflux::GetNormal {{{1*/
    6651013void Numericalflux:: GetNormal(double* normal,double xyz_list[4][3]){

  • TabularUnified issm/trunk/src/c/objects/Loads/Numericalflux.h

    r5911 r6026  
    7474                /*Numericalflux management:{{{1*/
    7575                void  GetNormal(double* normal,double xyz_list[4][3]);
    76                 ElementMatrix* CreateKMatrixInternal(void);
    77                 ElementMatrix* CreateKMatrixBoundary(void);
    78                 ElementVector* CreatePVectorInternal(void);
    79                 ElementVector* CreatePVectorBoundary(void);
     76                ElementMatrix* CreateKMatrixPrognostic(void);
     77                ElementMatrix* CreateKMatrixPrognosticInternal(void);
     78                ElementMatrix* CreateKMatrixPrognosticBoundary(void);
     79                ElementMatrix* CreateKMatrixBalancedthickness(void);
     80                ElementMatrix* CreateKMatrixBalancedthicknessInternal(void);
     81                ElementMatrix* CreateKMatrixBalancedthicknessBoundary(void);
     82                ElementMatrix* CreateKMatrixAdjointBalancedthickness(void);
     83                ElementMatrix* CreateKMatrixAdjointBalancedthicknessInternal(void);
     84                ElementMatrix* CreateKMatrixAdjointBalancedthicknessBoundary(void);
     85                ElementVector* CreatePVectorPrognostic(void);
     86                ElementVector* CreatePVectorPrognosticInternal(void);
     87                ElementVector* CreatePVectorPrognosticBoundary(void);
     88                ElementVector* CreatePVectorBalancedthickness(void);
     89                ElementVector* CreatePVectorBalancedthicknessInternal(void);
     90                ElementVector* CreatePVectorBalancedthicknessBoundary(void);
     91                ElementVector* CreatePVectorAdjointBalancedthickness(void);
     92                ElementVector* CreatePVectorAdjointBalancedthicknessInternal(void);
     93                ElementVector* CreatePVectorAdjointBalancedthicknessBoundary(void);
    8094                /*}}}*/
    8195
Note: See TracChangeset for help on using the changeset viewer.