Changeset 4904


Ignore:
Timestamp:
07/30/10 11:09:56 (15 years ago)
Author:
seroussi
Message:

now compute vx and vy average before spawning tria (prog, bathick and balvel)

Location:
issm/trunk/src/c/objects
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r4899 r4904  
    19621962        if(!onbed)return;
    19631963
     1964        /*Depth Averaging Vx and Vy*/
     1965        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
     1966        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
     1967
    19641968        /*Spawn Tria element from the base of the Penta: */
    19651969        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    … …  
    19901994        /*Is this element on the bed? :*/
    19911995        if(!onbed)return;
     1996
     1997        /*Depth Averaging Vx and Vy*/
     1998        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
     1999        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
    19922000
    19932001        /*Spawn Tria element from the base of the Penta: */
    … …  
    26912699/*}}}*/
    26922700/*FUNCTION Penta::CreateKMatrixPrognostic {{{1*/
    2693 
    26942701void  Penta::CreateKMatrixPrognostic(Mat Kgg){
    26952702
    … …  
    27102717        /*Is this element on the bed? :*/
    27112718        if(!onbed)return;
     2719
     2720        /*Depth Averaging Vx and Vy*/
     2721        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
     2722        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
    27122723
    27132724        /*Spawn Tria element from the base of the Penta: */
    … …  
    27202731/*}}}*/
    27212732/*FUNCTION Penta::CreateKMatrixSlope {{{1*/
    2722 
    27232733void  Penta::CreateKMatrixSlope(Mat Kgg){
    27242734
    … …  
    36923702/*}}}*/
    36933703/*FUNCTION Penta::CreatePVectorPrognostic {{{1*/
    3694 
    36953704void Penta::CreatePVectorPrognostic( Vec pg){
    36963705

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r4900 r4904  
    23212321        double  v_gauss[2]={0.0};
    23222322
    2323 
    23242323        double  K[2][2]={0.0};
    23252324        double  KDL[2][2]={0.0};
    23262325        int     dofs[2]={0,1};
    23272326        int     found=0;
     2327        int     dim;
    23282328
    23292329        /*inputs: */
    … …  
    23382338        /*retrieve some parameters: */
    23392339        this->parameters->FindParam(&artdiff,ArtDiffEnum);
     2340        this->parameters->FindParam(&dim,DimEnum);
    23402341
    23412342        /*Retrieve all inputs we will be needed*/
    2342         vxaverage_input=inputs->GetInput(VxAverageEnum);
    2343         vyaverage_input=inputs->GetInput(VyAverageEnum);
     2343        if(dim==2){
     2344                vxaverage_input=inputs->GetInput(VxEnum);
     2345                vyaverage_input=inputs->GetInput(VyEnum);
     2346        }
     2347        else{
     2348                vxaverage_input=inputs->GetInput(VxAverageEnum);
     2349                vyaverage_input=inputs->GetInput(VyAverageEnum);
     2350        }
    23442351
    23452352        //Create Artificial diffusivity once for all if requested
    … …  
    24762483        int     dofs[1]={0};
    24772484        int     found;
     2485        int     dim;
    24782486
    24792487        /*inputs: */
    … …  
    24842492        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    24852493        GetDofList(&doflist[0],&numberofdofspernode);
     2494        this->parameters->FindParam(&dim,DimEnum);
    24862495
    24872496        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    … …  
    24892498
    24902499        /*Retrieve all inputs we will be needed*/
    2491         vx_input=inputs->GetInput(VxEnum);
    2492         vy_input=inputs->GetInput(VyEnum);
     2500        if(dim==2){
     2501                vx_input=inputs->GetInput(VxEnum);
     2502                vy_input=inputs->GetInput(VyEnum);
     2503        }
    24932504
    24942505        /* Start  looping on the number of gaussian points: */
    … …  
    25892600        int     dofs[2]={0,1};
    25902601        int     found=0;
     2602        int     dim;
    25912603
    25922604        /*inputs: */
    25932605        Input* vxaverage_input=NULL;
    25942606        Input* vyaverage_input=NULL;
     2607        Input* surface_input=NULL;
    25952608        bool artdiff;
    25962609
    25972610        /*retrieve some parameters: */
    25982611        this->parameters->FindParam(&artdiff,ArtDiffEnum);
     2612        this->parameters->FindParam(&dim,DimEnum);
    25992613
    26002614        /* Get node coordinates and dof list: */
    … …  
    26022616        GetDofList(&doflist[0],&numberofdofspernode);
    26032617
     2618        /*Retrieve all inputs we will be needed*/
     2619        surface_input=inputs->GetInput(SurfaceEnum);
     2620        if(dim==2){
     2621                vxaverage_input=inputs->GetInput(VxEnum);
     2622                vyaverage_input=inputs->GetInput(VyEnum);
     2623        }
     2624        else{
     2625                vxaverage_input=inputs->GetInput(VxAverageEnum);
     2626                vyaverage_input=inputs->GetInput(VyAverageEnum);
     2627        }
     2628
    26042629        /*Modify z so that it reflects the surface*/
    2605         inputs->GetParameterValues(&surface_list[0],&gaussgrids[0][0],3,SurfaceEnum);
     2630        surface_input->GetParameterValues(&surface_list[0],&gaussgrids[0][0],3);
    26062631        for(i=0;i<numgrids;i++) xyz_list[i][2]=surface_list[i];
    26072632
    26082633        /*Get normal vector to the surface*/
    2609         inputs->GetParameterAverage(&nx,VxAverageEnum);
    2610         inputs->GetParameterAverage(&ny,VyAverageEnum);
     2634        vxaverage_input->GetParameterAverage(&nx);
     2635        vyaverage_input->GetParameterAverage(&ny);
    26112636        if(nx==0 && ny==0){
    26122637                SurfaceNormal(&surface_normal[0],xyz_list);
    … …  
    26212646        nx=nx/norm;
    26222647        ny=ny/norm;
    2623 
    2624         /*Retrieve all inputs we will be needed*/
    2625         vxaverage_input=inputs->GetInput(VxAverageEnum);
    2626         vyaverage_input=inputs->GetInput(VyAverageEnum);
    26272648
    26282649        //Create Artificial diffusivity once for all if requested
    … …  
    27042725        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    27052726
    2706 
    2707 cleanup_and_return:
    27082727        xfree((void**)&first_gauss_area_coord);
    27092728        xfree((void**)&second_gauss_area_coord);
    27102729        xfree((void**)&third_gauss_area_coord);
    27112730        xfree((void**)&gauss_weights);
    2712 
    27132731}
    27142732/*}}}*/
    … …  
    32843302        int     dofs[2]={0,1};
    32853303        int     found;
     3304        int     dim;
    32863305
    32873306        /*inputs: */
    … …  
    32933312        /*retrieve some parameters: */
    32943313        this->parameters->FindParam(&dt,DtEnum);
     3314        this->parameters->FindParam(&dim,DimEnum);
    32953315        this->parameters->FindParam(&artdiff,ArtDiffEnum);
    32963316
    … …  
    33003320
    33013321        /*Retrieve all inputs we will be needing: */
    3302         vxaverage_input=inputs->GetInput(VxAverageEnum);
    3303         vyaverage_input=inputs->GetInput(VyAverageEnum);
     3322        if(dim==2){
     3323                vxaverage_input=inputs->GetInput(VxEnum);
     3324                vyaverage_input=inputs->GetInput(VyEnum);
     3325        }
     3326        else{
     3327                vxaverage_input=inputs->GetInput(VxAverageEnum);
     3328                vyaverage_input=inputs->GetInput(VyAverageEnum);
     3329        }
    33043330
    33053331        //Create Artificial diffusivity once for all if requested
    … …  
    34053431        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    34063432
    3407 
    3408 cleanup_and_return:
    34093433        xfree((void**)&first_gauss_area_coord);
    34103434        xfree((void**)&second_gauss_area_coord);
    … …  
    34533477        int     dofs[1]={0};
    34543478        int     found;
     3479        int     dim;
    34553480
    34563481        /*inputs: */
    … …  
    34613486        /*retrieve some parameters: */
    34623487        this->parameters->FindParam(&dt,DtEnum);
     3488        this->parameters->FindParam(&dim,DimEnum);
    34633489
    34643490        /* Get node coordinates and dof list: */
    … …  
    34703496
    34713497        /*Retrieve all inputs we will be needed*/
    3472         vxaverage_input=inputs->GetInput(VxAverageEnum);
    3473         vyaverage_input=inputs->GetInput(VyAverageEnum);
     3498        if(dim==2){
     3499                vxaverage_input=inputs->GetInput(VxEnum);
     3500                vyaverage_input=inputs->GetInput(VyEnum);
     3501        }
     3502        else{
     3503                vxaverage_input=inputs->GetInput(VxAverageEnum);
     3504                vyaverage_input=inputs->GetInput(VyAverageEnum);
     3505        }
    34743506
    34753507        /* Start  looping on the number of gaussian points: */
    … …  
    35253557        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    35263558
    3527 
    3528 cleanup_and_return:
    35293559        xfree((void**)&first_gauss_area_coord);
    35303560        xfree((void**)&second_gauss_area_coord);
    … …  
    37103740void  Tria::CreatePVectorBalancedthickness_CG(Vec pg ){
    37113741
    3712 
    37133742        /* local declarations */
    37143743        int             i,j;
    … …  
    37813810        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    37823811
    3783 cleanup_and_return:
    37843812        xfree((void**)&first_gauss_area_coord);
    37853813        xfree((void**)&second_gauss_area_coord);
    37863814        xfree((void**)&third_gauss_area_coord);
    37873815        xfree((void**)&gauss_weights);
    3788 
    37893816}
    37903817/*}}}*/
    … …  
    38663893        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    38673894
    3868 cleanup_and_return:
    38693895        xfree((void**)&first_gauss_area_coord);
    38703896        xfree((void**)&second_gauss_area_coord);
    38713897        xfree((void**)&third_gauss_area_coord);
    38723898        xfree((void**)&gauss_weights);
    3873 
    38743899}
    38753900/*}}}*/
    38763901/*FUNCTION Tria::CreatePVectorBalancedvelocities {{{1*/
    38773902void  Tria::CreatePVectorBalancedvelocities(Vec pg){
    3878 
    38793903
    38803904        /* local declarations */
    … …  
    48144838        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    48154839
    4816 cleanup_and_return:
    48174840        xfree((void**)&first_gauss_area_coord);
    48184841        xfree((void**)&second_gauss_area_coord);
    … …  
    49014924        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    49024925
    4903 cleanup_and_return:
    49044926        xfree((void**)&first_gauss_area_coord);
    49054927        xfree((void**)&second_gauss_area_coord);
    … …  
    54705492        /*OK, now we can call input method*/
    54715493        this->inputs->GetParameterValue(&value, &gauss_tria[0],enumtype);
     5494
     5495        /*Assign output pointers:*/
     5496        *pvalue=value;
     5497}
     5498/*}}}*/
     5499/*FUNCTION Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in) {{{1*/
     5500void Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in){
     5501
     5502        /*Output*/
     5503        double value;
     5504
     5505        /*Intermediaries*/
     5506        const int numnodes=3;
     5507        int       grid1=-1,grid2=-1;
     5508        int       grid3;
     5509        int       i;
     5510        double    gauss_tria[numnodes];
     5511
     5512        /*go through 3 nodes (all nodes for tria) and identify 1st and 2nd nodes: */
     5513        ISSMASSERT(nodes);
     5514        for(i=0;i<numnodes;i++){
     5515                if (node1==nodes[i]) grid1=i;
     5516                if (node2==nodes[i]) grid2=i;
     5517        }
     5518
     5519        /*Reverse grid1 and 2 if necessary*/
     5520        if (grid1>grid2){
     5521
     5522                /*Reverse grid1 and grid2*/
     5523                grid3=grid1; grid1=grid2; grid2=grid3;
     5524
     5525                /*Change segment gauss point (between -1 and +1)*/
     5526                gauss_seg=-gauss_seg;
     5527        }
     5528
     5529        /*Build Triangle Gauss point*/
     5530        if (grid1==0 && grid2==1){
     5531                /*gauss_seg is between 0 and 1*/
     5532                gauss_tria[0]=0.5*(1.-gauss_seg);
     5533                gauss_tria[1]=1.-0.5*(1.-gauss_seg);
     5534                gauss_tria[2]=0.;
     5535        }
     5536        else if (grid1==0 && grid2==2){
     5537                /*gauss_seg is between 0 and 2*/
     5538                gauss_tria[0]=0.5*(1.-gauss_seg);
     5539                gauss_tria[1]=0.;
     5540                gauss_tria[2]=1.-0.5*(1.-gauss_seg);
     5541        }
     5542        else if (grid1==1 && grid2==2){
     5543                /*gauss_seg is between 1 and 2*/
     5544                gauss_tria[0]=0.;
     5545                gauss_tria[1]=0.5*(1.-gauss_seg);
     5546                gauss_tria[2]=1.-0.5*(1.-gauss_seg);
     5547        }
     5548        else{
     5549                ISSMERROR("The 2 nodes provided are either the same or did not match current Tria nodes");
     5550        }
     5551
     5552        /*OK, now we can call input method*/
     5553        input_in->GetParameterValue(&value, &gauss_tria[0]);
    54725554
    54735555        /*Assign output pointers:*/

  • issm/trunk/src/c/objects/Elements/Tria.h

    r4899 r4904  
    145145                void    GetParameterValue(double* pvalue,Node* node,int enumtype);
    146146                void    GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype);
     147                void    GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in);
    147148                void      GetSolutionFromInputsDiagnosticHoriz(Vec solution);
    148149                void      GetSolutionFromInputsAdjointHoriz(Vec solution);

  • issm/trunk/src/c/objects/Inputs/Input.h

    r4899 r4904  
    1212class Node;
    1313class ElementResult;
    14 #include "../Node.h"
    1514/*}}}*/
    1615

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

    r4765 r4904  
    401401        double dt;
    402402        int    found;
     403        int    dim;
    403404
    404405        /*dynamic objects pointed to by hooks: */
    405406        Node**  nodes=NULL;
     407        Tria*   tria=NULL;
     408        Input*  vxaverage_input=NULL;
     409        Input*  vyaverage_input=NULL;
    406410
    407411        /*Retrieve parameters: */
    408412        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     413        parameters->FindParam(&dim,DimEnum);
    409414
    410415        /*recover objects from hooks: */
    411416        nodes=(Node**)hnodes->deliverp();
     417        tria=(Tria*)helement->delivers();
    412418
    413419        /*recover parameters: */
    … …  
    428434        GetNormal(&normal[0],xyz_list);
    429435
     436        /*Retrieve all inputs we will be needing: */
     437        if(dim==2){
     438                vxaverage_input=tria->inputs->GetInput(VxEnum);
     439                vyaverage_input=tria->inputs->GetInput(VyEnum);
     440        }
     441
    430442        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    431443        num_gauss=2;
    … …  
    443455
    444456                //Get vx, vy and v.n
    445                 this->GetParameterValue(&vx,gauss_coord,VxAverageEnum);
    446                 this->GetParameterValue(&vy,gauss_coord,VyAverageEnum);
     457                this->GetParameterValue(&vx,gauss_coord,vxaverage_input);
     458                this->GetParameterValue(&vy,gauss_coord,vyaverage_input);
    447459               
    448460                UdotN=vx*normal[0]+vy*normal[1];
    … …  
    518530        int    dofs[1]={0};
    519531        int    found;
     532        int    dim;
    520533
    521534        /*dynamic objects pointed to by hooks: */
    522535        Node**  nodes=NULL;
     536        Tria*   tria=NULL;
     537        Input*  vxaverage_input=NULL;
     538        Input*  vyaverage_input=NULL;
    523539
    524540        /*Retrieve parameters: */
    525541        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     542        parameters->FindParam(&dim,DimEnum);
    526543
    527544        /*recover objects from hooks: */
    528545        nodes=(Node**)hnodes->deliverp();
     546        tria=(Tria*)helement->delivers();
    529547
    530548        /*recover parameters: */
    … …  
    545563        GetNormal(&normal[0],xyz_list);
    546564
     565        /*Retrieve all inputs we will be needing: */
     566        if(dim==2){
     567                vxaverage_input=tria->inputs->GetInput(VxEnum);
     568                vyaverage_input=tria->inputs->GetInput(VyEnum);
     569        }
     570
    547571        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    548         this->GetParameterValue(&mean_vx,0.,VxAverageEnum);
    549         this->GetParameterValue(&mean_vy,0.,VyAverageEnum);
     572        this->GetParameterValue(&mean_vx,0.,vxaverage_input);
     573        this->GetParameterValue(&mean_vy,0.,vyaverage_input);
    550574        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    551575        if (UdotN<=0){
    … …  
    569593
    570594                //Get vx, vy and v.n
    571                 this->GetParameterValue(&vx,gauss_coord,VxAverageEnum);
    572                 this->GetParameterValue(&vy,gauss_coord,VyAverageEnum);
     595                this->GetParameterValue(&vx,gauss_coord,vxaverage_input);
     596                this->GetParameterValue(&vy,gauss_coord,vyaverage_input);
    573597
    574598                UdotN=vx*normal[0]+vy*normal[1];
    … …  
    645669        int    dofs[1]={0};
    646670        int    found;
     671        int    dim;
    647672
    648673        /*dynamic objects pointed to by hooks: */
    649674        Node**  nodes=NULL;
     675        Tria*   tria=NULL;
     676        Input*  vxaverage_input=NULL;
     677        Input*  vyaverage_input=NULL;
     678        Input*  thickness_input=NULL;
    650679
    651680        /*recover objects from hooks: */
    652681        nodes=(Node**)hnodes->deliverp();
     682        tria=(Tria*)helement->delivers();
    653683
    654684        /*Retrieve parameters: */
    655685        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     686        parameters->FindParam(&dim,DimEnum);
    656687
    657688        /*recover parameters: */
    … …  
    667698        }
    668699       
     700        /*Retrieve all inputs we will be needing: */
     701        if(dim==2){
     702                vxaverage_input=tria->inputs->GetInput(VxEnum);
     703                vyaverage_input=tria->inputs->GetInput(VyEnum);
     704        }
     705        thickness_input=tria->inputs->GetInput(ThicknessEnum);
     706
    669707        /* Get node coordinates, dof list and normal vector: */
    670708        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    … …  
    673711
    674712        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    675         this->GetParameterValue(&mean_vx,0.,VxAverageEnum);
    676         this->GetParameterValue(&mean_vy,0.,VyAverageEnum);
     713        this->GetParameterValue(&mean_vx,0.,vxaverage_input);
     714        this->GetParameterValue(&mean_vy,0.,vyaverage_input);
    677715        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    678716        if (UdotN>0){
    … …  
    696734
    697735                //Get vx, vy and v.n
    698                 this->GetParameterValue(&vx,gauss_coord,VxAverageEnum);
    699                 this->GetParameterValue(&vy,gauss_coord,VyAverageEnum);
    700                 this->GetParameterValue(&thickness,gauss_coord,ThicknessEnum);
     736                this->GetParameterValue(&vx,gauss_coord,vxaverage_input);
     737                this->GetParameterValue(&vy,gauss_coord,vyaverage_input);
     738                this->GetParameterValue(&thickness,gauss_coord,thickness_input);
    701739
    702740                UdotN=vx*normal[0]+vy*normal[1];
    … …  
    846884
    847885        /*output: */
    848         double value1,value2,value;
    849         int    type;
     886        double value;
    850887
    851888        /*dynamic objects pointed to by hooks: */
    … …  
    853890        Node** nodes=NULL;
    854891
    855         /*recover type: */
    856         inputs->GetParameterValue(&type,TypeEnum);
    857 
    858 
    859892        /*recover objects from hooks: */
    860893        tria=(Tria*)helement->delivers();
    … …  
    868901}
    869902/*}}}*/
     903/*FUNCTION Numericalflux::GetParameterValue {{{1*/
     904void Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in){
     905
     906        /*output: */
     907        double value;
     908
     909        /*dynamic objects pointed to by hooks: */
     910        Tria*  tria=NULL;
     911        Node** nodes=NULL;
     912
     913        /*recover objects from hooks: */
     914        tria=(Tria*)helement->delivers();
     915        nodes=(Node**)hnodes->deliverp();
     916
     917        /*Get value on Element 1*/
     918        tria->GetParameterValue(&value,nodes[0],nodes[1],gauss_coord,input_in);
     919
     920        /*Assign output pointers:*/
     921        *pvalue=value;
     922}
     923/*}}}*/

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

    r4575 r4904  
    7070                void  GetNormal(double* normal,double xyz_list[4][3]);
    7171                void  GetParameterValue(double* pvalue, double gauss_coord,int enumtype);
     72                void  GetParameterValue(double* pvalue, double gauss_coord,Input* input_in);
    7273                void  CreateKMatrixInternal(Mat Kgg);
    7374                void  CreateKMatrixBoundary(Mat Kgg);
Note: See TracChangeset for help on using the changeset viewer.