Changeset 4015


Ignore:
Timestamp:
06/03/10 13:15:38 (15 years ago)
Author:
Eric.Larour
Message:

Completed UpdateInputsFromSolution for diagnostic

Location:
issm/trunk/src/c
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Makefile.am

    r4013 r4015  
    963963bin_PROGRAMS =
    964964else
    965 bin_PROGRAMS = DiagnosticAnalysis.exe ThermalAnalysis.exe PrognosticAnalysis.exe Prognostic2Analysis.exe BalancedthicknessAnalysis.exe Balancedthickness2Analysis.exe BalancedvelocitiesAnalysis.exe TransientAnalysis.exe SteadystateAnalysis.exe SlopecomputeAnalysis.exe
     965bin_PROGRAMS = DiagnosticAnalysis.exe
     966dnl bin_PROGRAMS = DiagnosticAnalysis.exe ThermalAnalysis.exe PrognosticAnalysis.exe Prognostic2Analysis.exe BalancedthicknessAnalysis.exe Balancedthickness2Analysis.exe BalancedvelocitiesAnalysis.exe TransientAnalysis.exe SteadystateAnalysis.exe SlopecomputeAnalysis.exe
    966967endif
    967968

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

    r4013 r4015  
    574574                        UpdateInputsFromSolutionDiagnosticHoriz( solution,analysis_type,sub_analysis_type);
    575575                }
     576                else if (sub_analysis_type==StokesAnalysisEnum){
     577
     578                        UpdateInputsFromSolutionDiagnosticStokes( solution,analysis_type,sub_analysis_type);
     579                }
    576580                else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
    577581
    … …  
    621625        double       vx[numvertices];
    622626        double       vy[numvertices];
    623 
    624627        int          dummy;
     628        double       pressure[numvertices];
     629        double       surface[numvertices];
     630        double       rho_ice,g;
     631        double       xyz_list[numvertices][3];
     632        double       gauss[numvertices][numvertices]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
     633
    625634       
    626635        /*Get dof list: */
    627636        GetDofList(&doflist[0],&dummy);
     637
     638        /*Get node data: */
     639        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    628640
    629641        /*Use the dof list to index into the solution vector: */
    … …  
    638650        }
    639651
     652        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D,
     653         *so the pressure is just the pressure at the z elevation: */
     654        rho_ice=matpar->GetRhoIce();
     655        g=matpar->GetG();
     656        inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum);
     657       
     658        for(i=0;i<numvertices;i++){
     659                pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
     660        }
    640661        /*Now, we have to move the previous Vx and Vy inputs  to old
    641662         * status, otherwise, we'll wipe them off: */
    642663        this->inputs->ChangeEnum(VxEnum,VxOldEnum);
    643664        this->inputs->ChangeEnum(VyEnum,VyOldEnum);
     665        this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
    644666
    645667        /*Add vx and vy as inputs to the tria element: */
    646668        this->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
    647669        this->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
     670        this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
     671}
     672
     673/*}}}*/
     674/*FUNCTION Penta::UpdateInputsFromSolutionDiagnosticStokes {{{1*/
     675void  Penta::UpdateInputsFromSolutionDiagnosticStokes(double* solution, int analysis_type, int sub_analysis_type){
     676       
     677       
     678        int i;
     679
     680        const int    numvertices=6;
     681        const int    numdofpervertex=4;
     682        const int    numdof=numdofpervertex*numvertices;
     683       
     684        int          doflist[numdof];
     685        double       values[numdof];
     686        double       vx[numvertices];
     687        double       vy[numvertices];
     688        double       vz[numvertices];
     689        double       pressure[numvertices];
     690        int          dummy;
     691        double       stokesreconditioning;
     692
     693        /*Get dof list: */
     694        GetDofList(&doflist[0],&dummy);
     695
     696        /*Use the dof list to index into the solution vector: */
     697        for(i=0;i<numdof;i++){
     698                values[i]=solution[doflist[i]];
     699        }
     700
     701        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
     702        for(i=0;i<numvertices;i++){
     703                vx[i]=values[i*numdofpervertex+0];
     704                vy[i]=values[i*numdofpervertex+1];
     705                vz[i]=values[i*numdofpervertex+2];
     706                pressure[i]=values[i*numdofpervertex+3];
     707        }
     708
     709        /*Recondition pressure: */
     710        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
     711        for(i=0;i<numvertices;i++){
     712                pressure[i]=pressure[i]*stokesreconditioning;
     713        }
     714
     715       
     716        /*Now, we have to move the previous inputs  to old
     717         * status, otherwise, we'll wipe them off: */
     718        this->inputs->ChangeEnum(VxEnum,VxOldEnum);
     719        this->inputs->ChangeEnum(VyEnum,VyOldEnum);
     720        this->inputs->ChangeEnum(VzEnum,VzOldEnum);
     721        this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
     722
     723        /*Add vx and vy as inputs to the tria element: */
     724        this->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
     725        this->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
     726        this->inputs->AddInput(new PentaVertexInput(VzEnum,vz));
     727        this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
    648728}
    649729

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

    r4001 r4015  
    155155                void  UpdateInputsFromSolution(double* solution, int analysis_type, int sub_analysis_type);
    156156                void  UpdateInputsFromSolutionDiagnosticHoriz( double* solution,int analysis_type,int sub_analysis_type);
     157                void  UpdateInputsFromSolutionDiagnosticStokes( double* solution,int analysis_type,int sub_analysis_type);
    157158                void  UpdateInputsFromSolutionSlopeCompute( double* solution,int analysis_type,int sub_analysis_type);
    158159                void  UpdateInputsFromSolutionPrognostic( double* solution,int analysis_type,int sub_analysis_type);

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

    r4013 r4015  
    537537        double       vx[numvertices];
    538538        double       vy[numvertices];
     539        double       pressure[numvertices];
     540        double       thickness[numvertices];
     541        double       rho_ice,g;
     542        double       gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
     543
    539544
    540545        int          dummy;
    … …  
    554559        }
    555560
     561        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
     562         *so the pressure is just the pressure at the bedrock: */
     563        rho_ice=matpar->GetRhoIce();
     564        g=matpar->GetG();
     565        inputs->GetParameterValues(&thickness[0],&gauss[0][0],3,ThicknessEnum);
     566       
     567        for(i=0;i<numvertices;i++){
     568                pressure[i]=rho_ice*g*thickness[i];
     569        }
     570
    556571        /*Now, we have to move the previous Vx and Vy inputs  to old
    557572         * status, otherwise, we'll wipe them off: */
    558573        this->inputs->ChangeEnum(VxEnum,VxOldEnum);
    559574        this->inputs->ChangeEnum(VyEnum,VyOldEnum);
     575        this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
    560576
    561577        /*Add vx and vy as inputs to the tria element: */
    562578        this->inputs->AddInput(new TriaVertexInput(VxEnum,vx));
    563579        this->inputs->AddInput(new TriaVertexInput(VyEnum,vy));
     580        this->inputs->AddInput(new TriaVertexInput(PressureEnum,pressure));
     581
     582
    564583}
    565584

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

    r3956 r4015  
    4444                void    UpdateInputsFromConstant(int constant, int name){ISSMERROR("Not implemented yet!");}
    4545                void    UpdateInputsFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");}
    46 
    4746                void    UpdateInputsFromSolution(double* solution, int analysis_type, int sub_analysis_type){ISSMERROR("Not implemented yet!");}
    4847

  • issm/trunk/src/c/solutions/ResetBoundaryConditions.cpp

    r4013 r4015  
    77#include "../EnumDefinitions/EnumDefinitions.h"
    88
    9 void ResetBoundaryConditions(FemModel* femmodel, int analysis_type, Vec ug){
     9void ResetBoundaryConditions(FemModel* femmodel, int analysis_type, int sub_analysis_type){
    1010       
    1111        int verbose=0;
     12        Vec ug=NULL;
    1213                       
    1314        femmodel->parameters->FindParam(&verbose,VerboseEnum);
    1415        if(verbose)_printf_("%s\n"," updating boundary conditions...");
     16                       
     17        GetSolutionFromInputsx( &ug, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,  femmodel->parameters, analysis_type,sub_analysis_type);
    1518
     19        /*set current analysis: */
    1620        femmodel->SetCurrentAnalysis(analysis_type);
    1721       
    18         /*First, for this analysis_type, free existing boundary condition vectors: */
     22        /*For this analysis_type, free existing boundary condition vectors: */
    1923
    2024        //global dof set
    … …  
    2933        Reducevectorgtosx(&femmodel->ys[analysis_counter],femmodel->yg[analysis_counter]->vector,femmodel->nodesets[analysis_counter]);
    3034
     35        /*Free ressources:*/
     36        VecFree(&ug);
     37
    3138}

  • issm/trunk/src/c/solutions/diagnostic.cpp

    r4014 r4015  
    6666
    6767        _printf_("create finite element model, using analyses types statically defined above:\n");
    68         femmodel=new FemModel(iomodel,analyses,5);
     68        femmodel=new FemModel(fid,analyses,5);
    6969       
    7070        /*get parameters: */

  • issm/trunk/src/c/solutions/diagnostic_core.cpp

    r4013 r4015  
    1212
    1313void* diagnostic_core(FemModel* femmodel){
    14 
    15         /*solutions: */
    16         Vec ug=NULL;
    1714
    1815        /*parameters: */
    … …  
    3936                ReinitializeInputx(femmodel,VyEnum,QmuVyEnum);
    4037                ReinitializeInputx(femmodel,VzEnum,QmuVzEnum);
     38                ReinitializeInputx(femmodel,PressureEnum,QmuPressureEnum);
    4139        }
    4240
    … …  
    5149        }
    5250               
    53        
    5451        if(ishutter){
    5552                       
    5653                if(verbose)_printf_("%s\n"," computing hutter velocities...");
    57                 solve_linear(&ug,femmodel,DiagnosticAnalysisEnum,HutterAnalysisEnum);
     54                solve_linear(NULL,femmodel,DiagnosticAnalysisEnum,HutterAnalysisEnum);
    5855               
    59                 if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticHorizAnalysisEnum,ug);
    60                 VecFree(&ug);
     56                if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,HorizAnalysisEnum);
    6157        }
    6258
    … …  
    8076                        if(verbose)_printf_("%s\n"," update boundary conditions for stokes using velocities previously computed...");
    8177                        GetSolutionFromInputsx( &ug, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,  femmodel->parameters, DiagnosticAnalysisEnum,StokesAnalysisEnum);
    82                         ResetBoundaryConditions(femmodel,DiagnosticStokesAnalysisEnum,ug);
    83                         VecFree(&ug);
     78                        ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,StokesAnalysisEnum);
    8479
    8580                        if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ...");
    … …  
    8782                }
    8883        }
    89        
    90         /*Free ressources: */
    91         VecFree(&ug);
    9284}

  • issm/trunk/src/c/solutions/slope_core.cpp

    r4013 r4015  
    1616        bool isstokes;
    1717        bool ishutter;
    18 
    19         /*output: */
    20         Vec slope=NULL;
    2118
    2219        /*Recover some parameters: */

  • issm/trunk/src/c/solutions/solutions.h

    r4013 r4015  
    3333void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,bool conserve_loads,int analysis_type,int sub_analysis_type);
    3434void diagnostic_core_linear(Vec* ppg,FemModel* fem,int  analysis_type,int sub_analysis_type);
    35 void solver_linear(Vec* ppg,FemModel* fem,int  analysis_type,int sub_analysis_type);
     35void solver_linear(Vec* pug, FemModel* femmodel,int  analysis_type,int sub_analysis_type);
    3636
    3737
    … …  
    5656void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,int analysis_type,int sub_analysis_type);
    5757//int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femmodel,char* filename);
    58 void ResetBoundaryConditions(FemModel* femmodel, int analysis_type, Vec ug);
     58void ResetBoundaryConditions(FemModel* femmodel, int analysis_type, int sub_analysis_type);
    5959
    6060#endif

  • issm/trunk/src/c/solutions/solver_linear.cpp

    r4013 r4015  
    88#include "../modules/modules.h"
    99
    10 void solver_linear(Vec* pug,FemModel* fem,int analysis_type,int sub_analysis_type){
     10void solver_linear(Vec* pug, FemModel* fem,int analysis_type,int sub_analysis_type){
    1111
    1212        /*parameters:*/
    … …  
    6161        xfree((void**)&solver_string);
    6262
    63         /*Assign output pointers if requested:*/
    64         if (pug) *pug=ug;
     63        /*Assign output pointers:*/
     64        if(pug)*pug=ug;
    6565}
Note: See TracChangeset for help on using the changeset viewer.