Changeset 4280


Ignore:
Timestamp:
06/28/10 16:17:30 (15 years ago)
Author:
Mathieu Morlighem
Message:

fixed diagnosticvert

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

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/UpdateElementsDiagnosticVert.cpp

    r4242 r4280  
    3232        IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface");
    3333        IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
    34         IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B");
    35         IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n");
     34        IoModelFetchData(&iomodel->melting_rate,NULL,NULL,iomodel_handle,"melting_rate");
     35        IoModelFetchData(&iomodel->accumulation_rate,NULL,NULL,iomodel_handle,"accumulation_rate");
    3636        IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
    3737        IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
     
    5959        xfree((void**)&iomodel->elementonsurface);
    6060        xfree((void**)&iomodel->elementonwater);
     61        xfree((void**)&iomodel->melting_rate);
     62        xfree((void**)&iomodel->accumulation_rate);
    6163        xfree((void**)&iomodel->vx);
    6264        xfree((void**)&iomodel->vy);
  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r4274 r4280  
    318318                InputUpdateFromSolutionDiagnosticHoriz( solution);
    319319        }
     320        else if (analysis_type==DiagnosticVertAnalysisEnum){
     321                InputUpdateFromSolutionDiagnosticVert( solution);
     322        }
    320323        else if (analysis_type==DiagnosticStokesAnalysisEnum){
    321324                InputUpdateFromSolutionDiagnosticStokes( solution);
     
    15751578
    15761579        /*Recover nodes ids needed to initialize the node hook.*/
    1577         for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook.
    1578                 penta_node_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
     1580        for(i=0;i<6;i++){
     1581                //go recover node ids, needed to initialize the node hook.
     1582                //WARNING: We assume P1 elements here!!!!!
     1583                penta_node_ids[i]=iomodel->nodecounter+(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
    15791584        }
    15801585
     
    19431948        inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum);
    19441949       
     1950        for(i=0;i<numvertices;i++){
     1951                pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
     1952        }
     1953        /*Now, we have to move the previous Vx and Vy inputs  to old
     1954         * status, otherwise, we'll wipe them off: */
     1955        this->inputs->ChangeEnum(VxEnum,VxOldEnum);
     1956        this->inputs->ChangeEnum(VyEnum,VyOldEnum);
     1957        this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
     1958
     1959        /*Add vx and vy as inputs to the tria element: */
     1960        this->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
     1961        this->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
     1962        this->inputs->AddInput(new TriaVertexInput(VelEnum,vel));
     1963        this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
     1964}
     1965
     1966/*}}}*/
     1967/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticVert {{{1*/
     1968void  Penta::InputUpdateFromSolutionDiagnosticVert(double* solution){
     1969
     1970        int i;
     1971
     1972        const int    numvertices=6;
     1973        const int    numdofpervertex=1;
     1974        const int    numdof=numdofpervertex*numvertices;
     1975
     1976        int          doflist[numdof];
     1977        double       values[numdof];
     1978        double       vx[numvertices];
     1979        double       vy[numvertices];
     1980        double       vz[numvertices];
     1981        double       vel[numvertices];
     1982        int          dummy;
     1983        double       pressure[numvertices];
     1984        double       surface[numvertices];
     1985        double       rho_ice,g;
     1986        double       xyz_list[numvertices][3];
     1987        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}};
     1988
     1989        Input*       VxInput=NULL;
     1990        double*      VxPtr=NULL;
     1991        Input*       VyInput=NULL;
     1992        double*      VyPtr=NULL;
     1993
     1994        /*Get dof list: */
     1995        GetDofList(&doflist[0],&dummy);
     1996
     1997        /*Get node data: */
     1998        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
     1999
     2000        /*Use the dof list to index into the solution vector: */
     2001        for(i=0;i<numdof;i++){
     2002                values[i]=solution[doflist[i]];
     2003        }
     2004
     2005        /*Ok, we have vz in values, fill in vz array: */
     2006        for(i=0;i<numvertices;i++){
     2007                vz[i]=values[i*numdofpervertex+0];
     2008        }
     2009
     2010        /*Get Vx and Vy*/
     2011        VxInput=inputs->GetInput(VxEnum);
     2012        if (VxInput){
     2013                if (VxInput->Enum()!=PentaVertexInputEnum) ISSMERROR("Cannot compute Vel as Vx is of type %s",EnumAsString(VxInput->Enum()));
     2014                VxInput->GetValuesPtr(&VxPtr,&dummy);
     2015                for(i=0;i<numvertices;i++) vx[i]=VxPtr[i];
     2016        }
     2017        else for(i=0;i<numvertices;i++) vx[i]=0.0;
     2018
     2019        VyInput=inputs->GetInput(VyEnum);
     2020        if (VyInput){
     2021                if (VyInput->Enum()!=PentaVertexInputEnum) ISSMERROR("Cannot compute Vel as Vy is of type %s",EnumAsString(VyInput->Enum()));
     2022                VyInput->GetValuesPtr(&VyPtr,&dummy);
     2023                for(i=0;i<numvertices;i++) vy[i]=VyPtr[i];
     2024        }
     2025        else for(i=0;i<numvertices;i++) vy[i]=0.0;
     2026
     2027        /*Now Compute vel*/
     2028        for(i=0;i<numvertices;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
     2029
     2030        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D,
     2031         *so the pressure is just the pressure at the z elevation: */
     2032        rho_ice=matpar->GetRhoIce();
     2033        g=matpar->GetG();
     2034        inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum);
     2035
    19452036        for(i=0;i<numvertices;i++){
    19462037                pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
  • issm/trunk/src/c/objects/Elements/Penta.h

    r4267 r4280  
    164164                void      GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);
    165165                void      GetStrainRateStokes(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);
    166                 Penta*    GetUpperElement(void);
     166                Penta*  GetUpperElement(void);
    167167                void      InputExtrude(int enum_type);
    168                 void      InputUpdateFromSolutionBalancedthickness( double* solutiong);
    169                 void      InputUpdateFromSolutionBalancedthickness2( double* solutiong);
    170                 void      InputUpdateFromSolutionBalancedvelocities( double* solutiong);
    171                 void      InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
    172                 void      InputUpdateFromSolutionDiagnosticStokes( double* solutiong);
    173                 void      InputUpdateFromSolutionPrognostic( double* solutiong);
    174                 void      InputUpdateFromSolutionPrognostic2(double* solutiong);
    175                 void      InputUpdateFromSolutionSlopeCompute( double* solutiong);
     168                void    InputUpdateFromSolutionBalancedthickness( double* solutiong);
     169                void    InputUpdateFromSolutionBalancedthickness2( double* solutiong);
     170                void    InputUpdateFromSolutionBalancedvelocities( double* solutiong);
     171                void    InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
     172                void    InputUpdateFromSolutionDiagnosticVert( double* solutiong);
     173                void    InputUpdateFromSolutionDiagnosticStokes( double* solutiong);
     174                void    InputUpdateFromSolutionPrognostic( double* solutiong);
     175                void    InputUpdateFromSolutionPrognostic2(double* solutiong);
     176                void    InputUpdateFromSolutionSlopeCompute( double* solutiong);
    176177                bool      IsInput(int name);
    177178                bool      IsOnSurface(void);
Note: See TracChangeset for help on using the changeset viewer.