Changeset 2722


Ignore:
Timestamp:
12/11/09 09:03:16 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added balanced velocities new solution (not working yet, but who knows... one day??)

Location:
issm/trunk/src
Files:
10 added
16 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp

    r2714 r2722  
    4242int PrognosticAnalysisEnum(void){           return          250; }
    4343int BalancedthicknessAnalysisEnum(void){    return          251; }
     44int BalancedvelocitiesAnalysisEnum(void){   return          252; }
    4445//melting
    4546int MeltingAnalysisEnum(void){              return          260; }

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

    r2714 r2722  
    4343int PrognosticAnalysisEnum(void);
    4444int BalancedthicknessAnalysisEnum(void);
     45int BalancedvelocitiesAnalysisEnum(void);
    4546//melting
    4647int MeltingAnalysisEnum(void);

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

    r2714 r2722  
    222222                                        ./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\
    223223                                        ./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\
     224                                        ./ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp\
     225                                        ./ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp\
     226                                        ./ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp\
     227                                        ./ModelProcessorx/Balancedvelocities/CreateParametersBalancedvelocities.cpp\
    224228                                        ./ModelProcessorx/Qmu/CreateParametersQmu.cpp\
    225229                                        ./Dofx/Dofx.h\
    … …  
    525529                                        ./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\
    526530                                        ./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\
     531                                        ./ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp\
     532                                        ./ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp\
     533                                        ./ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp\
     534                                        ./ModelProcessorx/Balancedvelocities/CreateParametersBalancedvelocities.cpp\
    527535                                        ./ModelProcessorx/Qmu/CreateParametersQmu.cpp\
    528536                                        ./Dofx/Dofx.h\
    … …  
    625633                                        ./parallel/prognostic_core.cpp\
    626634                                        ./parallel/balancedthickness_core.cpp\
     635                                        ./parallel/balancedvelocities_core.cpp\
    627636                                        ./parallel/transient_core.cpp\
    628637                                        ./parallel/transient_core_2d.cpp\
    … …  
    636645bin_PROGRAMS =
    637646else
    638 bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe balancedthickness.exe transient.exe steadystate.exe
     647bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe balancedthickness.exe balancedvelocities.exe transient.exe steadystate.exe
    639648endif
    640649
    … …  
    656665balancedthickness_exe_CXXFLAGS= -fPIC -D_PARALLEL_
    657666
     667balancedvelocities_exe_SOURCES = parallel/balancedvelocities.cpp
     668balancedvelocities_exe_CXXFLAGS= -fPIC -D_PARALLEL_
     669
    658670transient_exe_SOURCES = parallel/transient.cpp
    659671transient_exe_CXXFLAGS= -fPIC -D_PARALLEL_

  • issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp

    r2714 r2722  
    9494
    9595        }
     96        else if (iomodel->analysis_type==BalancedvelocitiesAnalysisEnum()){
     97
     98                CreateElementsNodesAndMaterialsBalancedvelocities(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
     99                CreateConstraintsBalancedvelocities(pconstraints,iomodel,iomodel_handle);
     100                CreateLoadsBalancedvelocities(ploads,iomodel,iomodel_handle);
     101                CreateParametersBalancedvelocities(pparameters,iomodel,iomodel_handle);
     102
     103        }
    96104        else{
    97105

  • issm/trunk/src/c/ModelProcessorx/IoModel.h

    r2714 r2722  
    177177        int* my_grids; /*! grids that belong to this cpu*/
    178178        double* my_bordergrids; /*! grids that belong to this cpu, and some other cpu also*/
    179         int*     penaltypartitioning;
     179        int*    penaltypartitioning;
    180180
    181181};
    182 
    183182
    184183
    185184        /*constructor and destructor: */
    186185        IoModel*        NewIoModel(void);
    187         void    DeleteIoModel( IoModel** pthis);
     186        void     DeleteIoModel( IoModel** pthis);
    188187
    189188        /*Echo: */
    190         void    IoModelEcho(IoModel* iomodel,int which_part,int rank);
     189        void  IoModelEcho(IoModel* iomodel,int which_part,int rank);
    191190
    192191        /*Initialization with matlab workspace data, or marshall binary data: */
    193         int         IoModelInit(IoModel** piomodel,ConstDataHandle iomodel_handle);
     192        int     IoModelInit(IoModel** piomodel,ConstDataHandle iomodel_handle);
    194193
    195194        /*Creation of fem datasets: general drivers*/
    196         void    CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
    197         void    CreateParameters(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
     195        void  CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
     196        void  CreateParameters(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
    198197
    199198       
    … …  
    203202        void    CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
    204203        void    CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
    205         void    CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
    206         void    CreateParametersDiagnosticHoriz(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
     204        void  CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
     205        void  CreateParametersDiagnosticHoriz(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
    207206
    208207        /*diagnostic vertical*/
    209208        void    CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
    210209        void    CreateConstraintsDiagnosticVert(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
    211         void    CreateLoadsDiagnosticVert(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
     210        void  CreateLoadsDiagnosticVert(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
    212211
    213212        /*diagnostic hutter*/
    214213        void    CreateElementsNodesAndMaterialsDiagnosticHutter(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
    215214        void    CreateConstraintsDiagnosticHutter(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
    216         void    CreateLoadsDiagnosticHutter(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
     215        void  CreateLoadsDiagnosticHutter(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
    217216
    218217        /*diagnostic stokes*/
    219218        void    CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
    220219        void    CreateConstraintsDiagnosticStokes(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
    221         void    CreateLoadsDiagnosticStokes(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
     220        void  CreateLoadsDiagnosticStokes(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
    222221
    223222        /*slope compute*/
    224223        void    CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
    225224        void    CreateConstraintsSlopeCompute(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
    226         void    CreateLoadsSlopeCompute(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
     225        void  CreateLoadsSlopeCompute(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
    227226
    228227        /*control:*/
    229         void    CreateParametersControl(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
     228        void  CreateParametersControl(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
    230229
    231230        /*thermal:*/
    232231        void    CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
    233232        void    CreateConstraintsThermal(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
    234         void    CreateLoadsThermal(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
    235         void    CreateParametersThermal(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
     233        void  CreateLoadsThermal(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
     234        void  CreateParametersThermal(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
    236235
    237236        /*melting:*/
    238237        void    CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
    239238        void    CreateConstraintsMelting(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
    240         void    CreateLoadsMelting(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
    241         void    CreateParametersMelting(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
     239        void  CreateLoadsMelting(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
     240        void  CreateParametersMelting(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
    242241
    243242        /*prognostic:*/
    244243        void    CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
    245244        void    CreateConstraintsPrognostic(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
    246         void    CreateLoadsPrognostic(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
    247         void    CreateParametersPrognostic(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
     245        void  CreateLoadsPrognostic(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
     246        void  CreateParametersPrognostic(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
    248247
    249248        /*balancedthickness:*/
    250249        void    CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
    251250        void    CreateConstraintsBalancedthickness(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
    252         void    CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
    253         void    CreateParametersBalancedthickness(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
     251        void  CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
     252        void  CreateParametersBalancedthickness(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
     253
     254        /*balancedvelocities:*/
     255        void    CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
     256        void    CreateConstraintsBalancedvelocities(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
     257        void  CreateLoadsBalancedvelocities(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
     258        void  CreateParametersBalancedvelocities(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
    254259
    255260        /*qmu: */
    256261        void CreateParametersQmu(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
    257262
    258 
    259263#endif  /* _IOMODEL_H */

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

    r2714 r2722  
    325325                CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type);
    326326        }
     327        else if (analysis_type==BalancedvelocitiesAnalysisEnum()){
     328
     329                CreateKMatrixBalancedvelocities( Kgg,inputs,analysis_type,sub_analysis_type);
     330        }
    327331        else if (analysis_type==ThermalAnalysisEnum()){
    328332
    … …  
    344348
    345349void  Penta::CreateKMatrixBalancedthickness(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
     350
     351        /*Collapsed formulation: */
     352        Tria*  tria=NULL;
     353
     354        /*If on water, skip: */
     355        if(onwater)return;
     356
     357        /*Is this element on the bed? :*/
     358        if(!onbed)return;
     359
     360        /*Spawn Tria element from the base of the Penta: */
     361        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     362        tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
     363        delete tria;
     364        return;
     365
     366}
     367/*}}}*/
     368/*FUNCTION CreateKMatrixBalancedvelocities {{{1*/
     369#undef __FUNCT__
     370#define __FUNCT__ "Penta::CreateKMatrixBalancedvelocities"
     371
     372void  Penta::CreateKMatrixBalancedvelocities(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    346373
    347374        /*Collapsed formulation: */
    … …  
    14101437                CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type);
    14111438        }
     1439        else if (analysis_type==BalancedvelocitiesAnalysisEnum()){
     1440
     1441                CreatePVectorBalancedvelocities( pg,inputs,analysis_type,sub_analysis_type);
     1442        }
    14121443        else if (analysis_type==ThermalAnalysisEnum()){
    14131444
    … …  
    14291460
    14301461void Penta::CreatePVectorBalancedthickness( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
     1462
     1463        /*Collapsed formulation: */
     1464        Tria*  tria=NULL;
     1465
     1466        /*If on water, skip: */
     1467        if(onwater)return;
     1468
     1469        /*Is this element on the bed? :*/
     1470        if(!onbed)return;
     1471
     1472        /*Spawn Tria element from the base of the Penta: */
     1473        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     1474        tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
     1475        delete tria;
     1476        return;
     1477}
     1478/*}}}*/
     1479/*FUNCTION CreatePVectorBalancedvelocities {{{1*/
     1480#undef __FUNCT__
     1481#define __FUNCT__ "Penta::CreatePVectorBalancedvelocities"
     1482
     1483void Penta::CreatePVectorBalancedvelocities( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
    14311484
    14321485        /*Collapsed formulation: */

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

    r2714 r2722  
    118118                void  CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    119119                void  CreatePVectorBalancedthickness( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
     120                void  CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
     121                void  CreatePVectorBalancedvelocities( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
    120122
    121123                void  CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type);

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

    r2714 r2722  
    303303
    304304        }
     305        else if (analysis_type==BalancedvelocitiesAnalysisEnum()){
     306
     307                CreateKMatrixBalancedvelocities( Kgg,inputs,analysis_type,sub_analysis_type);
     308
     309        }
    305310        else{
    306311                throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
    … …  
    447452                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness1[i][j];
    448453                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
     454
     455                if(numpar->artdiff){
     456
     457                        /* Compute artificial diffusivity */
     458                        KDL[0][0]=DL_scalar*K[0][0];
     459                        KDL[1][1]=DL_scalar*K[1][1];
     460
     461                        TripleMultiply( &Bprime[0][0],2,numdof,1,
     462                                                &KDL[0][0],2,2,0,
     463                                                &Bprime[0][0],2,numdof,0,
     464                                                &Ke_gg_gaussian[0][0],0);
     465
     466                        /* Add artificial diffusivity matrix */
     467                        for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     468
     469                }
     470
     471#ifdef _DEBUGELEMENTS_
     472                if(my_rank==RANK && id==ELID){
     473                        printf("      B:\n");
     474                        for(i=0;i<3;i++){
     475                                for(j=0;j<numdof;j++){
     476                                        printf("%g ",B[i][j]);
     477                                }
     478                                printf("\n");
     479                        }
     480                        printf("      Bprime:\n");
     481                        for(i=0;i<3;i++){
     482                                for(j=0;j<numdof;j++){
     483                                        printf("%g ",Bprime[i][j]);
     484                                }
     485                                printf("\n");
     486                        }
     487                }
     488#endif
     489        } // for (ig=0; ig<num_gauss; ig++)
     490
     491        /*Add Ke_gg to global matrix Kgg: */
     492        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     493
     494#ifdef _DEBUGELEMENTS_
     495        if(my_rank==RANK && id==ELID){
     496                printf("      Ke_gg erms:\n");
     497                for( i=0; i<numdof; i++){
     498                        for (j=0;j<numdof;j++){
     499                                printf("%g ",Ke_gg[i][j]);
     500                        }
     501                        printf("\n");
     502                }
     503                printf("      Ke_gg row_indices:\n");
     504                for( i=0; i<numdof; i++){
     505                        printf("%i ",doflist[i]);
     506                }
     507
     508        }
     509#endif
     510
     511cleanup_and_return:
     512        xfree((void**)&first_gauss_area_coord);
     513        xfree((void**)&second_gauss_area_coord);
     514        xfree((void**)&third_gauss_area_coord);
     515        xfree((void**)&gauss_weights);
     516
     517}
     518/*}}}*/
     519/*FUNCTION CreateKMatrixBalancedvelocities {{{1*/
     520#undef __FUNCT__
     521#define __FUNCT__ "Tria::CreateKMatrixBalancedvelocities"
     522void  Tria::CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
     523
     524        /* local declarations */
     525        int             i,j;
     526
     527        /* node data: */
     528        const int    numgrids=3;
     529        const int    NDOF1=1;
     530        const int    numdof=NDOF1*numgrids;
     531        double       xyz_list[numgrids][3];
     532        int          doflist[numdof];
     533        int          numberofdofspernode;
     534
     535        /* gaussian points: */
     536        int     num_gauss,ig;
     537        double* first_gauss_area_coord  =  NULL;
     538        double* second_gauss_area_coord =  NULL;
     539        double* third_gauss_area_coord  =  NULL;
     540        double* gauss_weights           =  NULL;
     541        double  gauss_weight;
     542        double  gauss_l1l2l3[3];
     543
     544        /* matrices: */
     545        double L[numgrids];
     546        double B[2][numgrids];
     547        double Bprime[2][numgrids];
     548        double DL[2][2]={0.0};
     549        double DLprime[2][2]={0.0};
     550        double DL_scalar;
     551        double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix
     552        double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
     553        double Ke_gg_velocities1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
     554        double Ke_gg_velocities2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
     555        double Jdettria;
     556
     557        /*input parameters for structural analysis (diagnostic): */
     558        double  surface_normal[3];
     559        double  nx,ny,norm;
     560        double  vxvy_list[numgrids][2]={0.0};
     561        double  vx_list[numgrids]={0.0};
     562        double  vy_list[numgrids]={0.0};
     563        double  dvx[2];
     564        double  dvy[2];
     565        double  vx,vy;
     566        double  dvxdx,dvydy;
     567        double  v_gauss[2]={0.0};
     568        double  K[2][2]={0.0};
     569        double  KDL[2][2]={0.0};
     570        int     dofs[2]={0,1};
     571        int     found=0;
     572
     573        ParameterInputs* inputs=NULL;
     574
     575        /*recover pointers: */
     576        inputs=(ParameterInputs*)vinputs;
     577
     578        /*recover extra inputs from users, at current convergence iteration: */
     579        found=inputs->Recover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
     580        if(!found)throw ErrorException(__FUNCT__," could not find velocity_average  in inputs!");
     581
     582        for(i=0;i<numgrids;i++){
     583                vx_list[i]=vxvy_list[i][0];
     584                vy_list[i]=vxvy_list[i][1];
     585        }
     586
     587        /* Get node coordinates and dof list: */
     588        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     589        GetDofList(&doflist[0],&numberofdofspernode);
     590
     591        /*Modify z so that it reflects the surface*/
     592        for(i=0;i<numgrids;i++) xyz_list[i][2]=s[i];
     593
     594        /*Get normal vector to the surface*/
     595        nx=(vx_list[0]+vx_list[1]+vx_list[2])/3;
     596        ny=(vy_list[0]+vy_list[1]+vy_list[2])/3;
     597        if(nx==0 && ny==0){
     598                SurfaceNormal(&surface_normal[0],xyz_list);
     599                nx=surface_normal[0];
     600                ny=surface_normal[1];
     601        }
     602        if(nx==0 && ny==0){
     603                nx=0;
     604                ny=1;
     605        }
     606        norm=pow( pow(nx,2)+pow(ny,2) , (double).5);
     607        nx=nx/norm;
     608        ny=ny/norm;
     609
     610        //Create Artificial diffusivity once for all if requested
     611        if(numpar->artdiff){
     612                //Get the Jacobian determinant
     613                gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0;
     614                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     615
     616                //Build K matrix (artificial diffusivity matrix)
     617                v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
     618                v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
     619
     620                K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!!
     621                K[1][1]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]);
     622        }
     623
     624        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     625        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     626
     627        /* Start  looping on the number of gaussian points: */
     628        for (ig=0; ig<num_gauss; ig++){
     629                /*Pick up the gaussian point: */
     630                gauss_weight=*(gauss_weights+ig);
     631                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     632                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     633                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     634
     635                /* Get Jacobian determinant: */
     636                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     637
     638                /*Get B  and B prime matrix: */
     639                GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
     640                GetBPrime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
     641
     642                //Get vx, vy and their derivatives at gauss point
     643                GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
     644                GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
     645
     646                GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3);
     647                GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3);
     648
     649                dvxdx=dvx[0];
     650                dvydy=dvy[1];
     651
     652                DL_scalar=gauss_weight*Jdettria;
     653
     654                DLprime[0][0]=DL_scalar*nx;
     655                DLprime[1][1]=DL_scalar*ny;
     656
     657                //Do the triple product tL*D*L.
     658                //Ke_gg_velocities=B'*DLprime*Bprime;
     659                TripleMultiply( &B[0][0],2,numdof,1,
     660                                        &DLprime[0][0],2,2,0,
     661                                        &Bprime[0][0],2,numdof,0,
     662                                        &Ke_gg_velocities2[0][0],0);
     663
     664                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     665                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_velocities2[i][j];
    449666
    450667                if(numpar->artdiff){
    … …  
    15501767                CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type);
    15511768        }
     1769        else if (analysis_type==BalancedvelocitiesAnalysisEnum()){
     1770
     1771                CreatePVectorBalancedvelocities( pg,inputs,analysis_type,sub_analysis_type);
     1772        }
    15521773        else{
    15531774                throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis ",analysis_type," not supported yet"));
    … …  
    16301851
    16311852                /* Get accumulation, melting and thickness at gauss point */
     1853                GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3);
     1854                GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3);
     1855
     1856                /* Add value into pe_g: */
     1857                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(accumulation_g-melting_g)*L[i];
     1858
     1859        } // for (ig=0; ig<num_gauss; ig++)
     1860
     1861        /*Add pe_g to global matrix Kgg: */
     1862        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
     1863
     1864cleanup_and_return:
     1865        xfree((void**)&first_gauss_area_coord);
     1866        xfree((void**)&second_gauss_area_coord);
     1867        xfree((void**)&third_gauss_area_coord);
     1868        xfree((void**)&gauss_weights);
     1869
     1870}
     1871/*}}}*/
     1872/*FUNCTION CreatePVectorBalancedvelocities {{{1*/
     1873#undef __FUNCT__
     1874#define __FUNCT__ "Tria::CreatePVectorBalancedvelocities"
     1875void  Tria::CreatePVectorBalancedvelocities(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){
     1876
     1877
     1878        /* local declarations */
     1879        int             i,j;
     1880
     1881        /* node data: */
     1882        const int    numgrids=3;
     1883        const int    NDOF1=1;
     1884        const int    numdof=NDOF1*numgrids;
     1885        double       xyz_list[numgrids][3];
     1886        int          doflist[numdof];
     1887        int          numberofdofspernode;
     1888
     1889        /* gaussian points: */
     1890        int     num_gauss,ig;
     1891        double* first_gauss_area_coord  =  NULL;
     1892        double* second_gauss_area_coord =  NULL;
     1893        double* third_gauss_area_coord  =  NULL;
     1894        double* gauss_weights           =  NULL;
     1895        double  gauss_weight;
     1896        double  gauss_l1l2l3[3];
     1897
     1898        /* matrix */
     1899        double pe_g[numgrids]={0.0};
     1900        double L[numgrids];
     1901        double Jdettria;
     1902
     1903        /*input parameters for structural analysis (diagnostic): */
     1904        double  accumulation_list[numgrids]={0.0};
     1905        double  accumulation_g;
     1906        double  melting_list[numgrids]={0.0};
     1907        double  melting_g;
     1908        int     dofs[1]={0};
     1909        int     found=0;
     1910
     1911        ParameterInputs* inputs=NULL;
     1912
     1913        /*recover pointers: */
     1914        inputs=(ParameterInputs*)vinputs;
     1915
     1916        /*recover extra inputs from users, at current convergence iteration: */
     1917        found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes);
     1918        if(!found)throw ErrorException(__FUNCT__," could not find accumulation in inputs!");
     1919        found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes);
     1920        if(!found)throw ErrorException(__FUNCT__," could not find melting in inputs!");
     1921
     1922        /* Get node coordinates and dof list: */
     1923        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     1924        GetDofList(&doflist[0],&numberofdofspernode);
     1925
     1926        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     1927        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     1928
     1929        /* Start  looping on the number of gaussian points: */
     1930        for (ig=0; ig<num_gauss; ig++){
     1931                /*Pick up the gaussian point: */
     1932                gauss_weight=*(gauss_weights+ig);
     1933                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     1934                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     1935                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     1936
     1937                /* Get Jacobian determinant: */
     1938                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     1939
     1940                /*Get L matrix: */
     1941                GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     1942
     1943                /* Get accumulation, melting at gauss point */
    16321944                GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3);
    16331945                GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3);

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

    r2714 r2722  
    109109                int   GetOnBed();
    110110
    111                 void          GetThicknessList(double* thickness_list);
    112                 void          GetBedList(double* bed_list);
     111                void  GetThicknessList(double* thickness_list);
     112                void  GetBedList(double* bed_list);
    113113                Object* copy();
    114114                void  NodeConfiguration(int* tria_node_ids,Node* tria_nodes[3],int* tria_node_offsets);
    … …  
    125125                void  CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    126126                void  CreatePVectorBalancedthickness(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);
     127                void  CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
     128                void  CreatePVectorBalancedvelocities(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);
    127129                double MassFlux(double* segment,double* ug);
    128130                double GetArea(void);

  • issm/trunk/src/c/parallel/ProcessResults.cpp

    r2718 r2722  
    9090        double* m_g_serial=NULL;
    9191        double* melting=NULL;
     92
     93        Vec     v_g=NULL;
     94        double* v_g_serial=NULL;
    9295
    9396        int numberofnodes;
    … …  
    117120        if(analysis_type==BalancedthicknessAnalysisEnum()){
    118121                fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum());
     122        }
     123        if(analysis_type==BalancedvelocitiesAnalysisEnum()){
     124                fem_p=model->GetFormulation(BalancedvelocitiesAnalysisEnum());
    119125        }
    120126        fem_t=model->GetFormulation(ThermalAnalysisEnum());
    … …  
    327333                        VecFree(&h_g);
    328334                }
     335                else if(strcmp(result->GetFieldName(),"v_g")==0){
     336                        /*easy, v_g is of size numberofnodes, on 1 dof, just repartition: */
     337                        result->GetField(&v_g);
     338                        VecToMPISerial(&v_g_serial,v_g);
     339                        fem_p->parameters->FindParam(&numberofnodes,"numberofnodes");
     340                        VecToMPISerial(&partition,fem_p->partition->vector);
     341
     342                        vel=(double*)xmalloc(numberofnodes*sizeof(double));
     343
     344                        for(i=0;i<numberofnodes;i++){
     345                                vel[i]=v_g_serial[(int)partition[i]];
     346                        }
     347
     348                        /*Ok, add pressure to newresults: */
     349                        newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vel",vel,numberofnodes);
     350                        newresults->AddObject(newresult);
     351
     352                        /*do some cleanup: */
     353                        xfree((void**)&v_g_serial);
     354                        xfree((void**)&vel);
     355                        xfree((void**)&partition);
     356                        VecFree(&v_g);
     357                }
    329358                else if(strcmp(result->GetFieldName(),"s_g")==0){
    330359                        /*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */

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

    r2714 r2722  
    1818void prognostic_core(DataSet* results,Model* model, ParameterInputs* inputs);
    1919void balancedthickness_core(DataSet* results,Model* model, ParameterInputs* inputs);
     20void balancedvelocities_core(DataSet* results,Model* model, ParameterInputs* inputs);
    2021void control_core(DataSet* results,Model* model, ParameterInputs* inputs);
    2122

  • issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp

    r2714 r2722  
    5959                numdofs=1;
    6060        }
     61        else if (analysis_type==BalancedvelocitiesAnalysisEnum()){
     62                numdofs=1;
     63        }
    6164        else throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis type: ",analysis_type,"  not implemented yet!"));
    6265

  • issm/trunk/src/m/classes/public/process_solve_options.m

    r2714 r2722  
    1717
    1818%check solution type is supported
    19 if ~ismemberi(analysis_type,{'diagnostic','prognostic','thermal','steadystate','parameters','transient','balancedthickness'}),
     19if ~ismemberi(analysis_type,{'diagnostic','prognostic','thermal','steadystate','parameters','transient',...
     20                'balancedthickness','balancedvelocities'}),
    2021        error(['process_solve_options error message: analysis_type ' analysis_type ' not supported yet!']);
    2122else
    2223        %convert to enum
    23         outoptions.analysis_type=eval([upper(analysis_type(1)) lower(analysis_type(2:end)) 'AnalysisEnum']);
     24        outoptions.analysis_type=eval([upper(analysis_type(1)) lower(analysis_type(2:end)) 'AnalysisEnum()']);
    2425end
    2526if ~ismemberi(sub_analysis_type,{'steady','transient','none','horiz','adjoint','gradient','inverse','vert',''}),
    … …  
    2728else
    2829        %convert to enum
    29         outoptions.sub_analysis_type=eval([upper(sub_analysis_type(1)) lower(sub_analysis_type(2:end)) 'AnalysisEnum']);
     30        outoptions.sub_analysis_type=eval([upper(sub_analysis_type(1)) lower(sub_analysis_type(2:end)) 'AnalysisEnum()']);
    3031end
    3132

  • issm/trunk/src/m/classes/public/solve.m

    r2714 r2722  
    7373elseif md.analysis_type==BalancedthicknessAnalysisEnum,
    7474        md=balancedthickness(md);
     75
     76elseif md.analysis_type==BalancedvelocitiesAnalysisEnum,
     77        md=balancedvelocities(md);
     78
    7579else
    7680        error('solution type not supported for this model configuration.');

  • issm/trunk/src/m/enum/AnalysisTypeFromEnum.m

    r2714 r2722  
    8989end
    9090
     91if enum==BalancedvelocitiesAnalysisEnum(),
     92        string='balancedvelocities';
     93end
     94
    9195if enum==MeltingAnalysisEnum(),
    9296        string='melting';

  • issm/trunk/src/m/solutions/jpl/processresults.m

    r1871 r2722  
    126126                        newresults(i).parameter=param_g;
    127127
     128                elseif strcmpi(resultsname{j},'v_g'),
     129
     130                        v_g=results(i).v_g;
     131                        newresults(i).step=results(i).step;
     132                        newresults(i).time=results(i).time;
     133                        newresults(i).vel=v_g;
     134
    128135                else
    129136
Note: See TracChangeset for help on using the changeset viewer.