Changeset 2714


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

Added Balancedthickness solution

Location:
issm/trunk/src
Files:
19 edited

Legend:

Unmodified
Added
Removed

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

    r2333 r2714  
    77
    88/*Datasets: */
    9 int DatasetsEnum(void){                 return          100; }
    10 int ElementsEnum(void){                 return          101; }
    11 int NodesEnum(void){                    return          102; }
    12 int ConstraintsEnum(void){              return          103; }
    13 int LoadsEnum(void){                    return          104; }
    14 int MaterialsEnum(void){                return          105; }
    15 int ParametersEnum(void){               return          106; }
    16 int ResultsEnum(void){                  return          107; }
     9int DatasetsEnum(void){                     return          100; }
     10int ElementsEnum(void){                     return          101; }
     11int NodesEnum(void){                        return          102; }
     12int ConstraintsEnum(void){                  return          103; }
     13int LoadsEnum(void){                        return          104; }
     14int MaterialsEnum(void){                    return          105; }
     15int ParametersEnum(void){                   return          106; }
     16int ResultsEnum(void){                      return          107; }
    1717
    1818/*ANALYSIS TYPES: */
    19 int AnalysisEnum(void){                 return          200; }
     19int AnalysisEnum(void){                     return          200; }
    2020//diagnostic
    21 int DiagnosticAnalysisEnum(void){       return          210; }
    22 int HorizAnalysisEnum(void){            return          211; }
    23 int StokesAnalysisEnum(void){           return          212; }
    24 int HutterAnalysisEnum(void){           return          213; }
    25 int VertAnalysisEnum(void){             return          214; }
     21int DiagnosticAnalysisEnum(void){           return          210; }
     22int HorizAnalysisEnum(void){                return          211; }
     23int StokesAnalysisEnum(void){               return          212; }
     24int HutterAnalysisEnum(void){               return          213; }
     25int VertAnalysisEnum(void){                 return          214; }
    2626//control
    27 int ControlAnalysisEnum(void){          return          220; }
    28 int AdjointAnalysisEnum(void){          return          221; }
    29 int InverseAnalysisEnum(void){          return          222; }
    30 int GradientAnalysisEnum(void){         return          223; }
     27int ControlAnalysisEnum(void){              return          220; }
     28int AdjointAnalysisEnum(void){              return          221; }
     29int InverseAnalysisEnum(void){              return          222; }
     30int GradientAnalysisEnum(void){             return          223; }
    3131//thermal
    32 int ThermalAnalysisEnum(void){          return          230; }
     32int ThermalAnalysisEnum(void){              return          230; }
    3333//transient
    34 int TransientAnalysisEnum(void){        return          231; }
     34int TransientAnalysisEnum(void){            return          231; }
    3535//slope
    36 int SlopeComputeAnalysisEnum(void){     return          240; }
    37 int SurfaceXAnalysisEnum(void){         return          241; }
    38 int SurfaceYAnalysisEnum(void){         return          242; }
    39 int BedXAnalysisEnum(void){             return          243; }
    40 int BedYAnalysisEnum(void){             return          244; }
     36int SlopeComputeAnalysisEnum(void){         return          240; }
     37int SurfaceXAnalysisEnum(void){             return          241; }
     38int SurfaceYAnalysisEnum(void){             return          242; }
     39int BedXAnalysisEnum(void){                 return          243; }
     40int BedYAnalysisEnum(void){                 return          244; }
    4141//prognostic
    42 int PrognosticAnalysisEnum(void){       return          250; }
     42int PrognosticAnalysisEnum(void){           return          250; }
     43int BalancedthicknessAnalysisEnum(void){    return          251; }
    4344//melting
    44 int MeltingAnalysisEnum(void){          return          260; }
     45int MeltingAnalysisEnum(void){              return          260; }
    4546//mesh2grid
    46 int Mesh2gridAnalysisEnum(void){        return          270; }
     47int Mesh2gridAnalysisEnum(void){            return          270; }
    4748//parameters
    48 int ParametersAnalysisEnum(void){       return          280; }
     49int ParametersAnalysisEnum(void){           return          280; }
    4950//steadystate
    50 int SteadystateAnalysisEnum(void){      return          281; }
     51int SteadystateAnalysisEnum(void){          return          281; }
    5152//none
    52 int NoneAnalysisEnum(void){             return          290; }
     53int NoneAnalysisEnum(void){                 return          290; }
    5354
    5455
    5556/*Formulations: */
    56 int FormulationEnum(void){              return          300; }
    57 int NoneFormulationEnum(void){          return          301; }
    58 int HutterFormulationEnum(void){        return          302; }
    59 int MacAyealFormulationEnum(void){      return          303; }
    60 int PattynFormulationEnum(void){        return          304; }
    61 int StokesFormulationEnum(void){        return          305; }
     57int FormulationEnum(void){                  return          300; }
     58int NoneFormulationEnum(void){              return          301; }
     59int HutterFormulationEnum(void){            return          302; }
     60int MacAyealFormulationEnum(void){          return          303; }
     61int PattynFormulationEnum(void){            return          304; }
     62int StokesFormulationEnum(void){            return          305; }
    6263
    6364/*OBJECTS: */
    64 int ObjectEnum(void){                   return          400; }
     65int ObjectEnum(void){                       return          400; }
    6566/*Elements: */
    66 int ElementEnum(void){                  return          410; }
    67 int TriaEnum(void){                     return          411; }
    68 int PentaEnum(void){                    return          412; }
    69 int SingEnum(void){                     return          414; }
    70 int BeamEnum(void){                     return          415; }
     67int ElementEnum(void){                      return          410; }
     68int TriaEnum(void){                         return          411; }
     69int PentaEnum(void){                        return          412; }
     70int SingEnum(void){                         return          414; }
     71int BeamEnum(void){                         return          415; }
    7172/*Grids: */
    72 int NodeEnum(void){                     return          420; }
     73int NodeEnum(void){                         return          420; }
    7374/*Loads: */
    74 int LoadEnum(void){                     return          430; }
    75 int IcefrontEnum(void){                 return          431; }
    76 int RiftfrontEnum(void){                return          432; }
    77 int PenpairEnum(void){                  return          433; }
    78 int PengridEnum(void){                  return          434; }
     75int LoadEnum(void){                         return          430; }
     76int IcefrontEnum(void){                     return          431; }
     77int RiftfrontEnum(void){                    return          432; }
     78int PenpairEnum(void){                      return          433; }
     79int PengridEnum(void){                      return          434; }
    7980/*Materials: */
    80 int MaterialEnum(void){                 return          440; }
    81 int MaticeEnum(void){                   return          441; }
    82 int MatparEnum(void){                   return          442; }
    83 int NumparEnum(void){                   return          443; }
     81int MaterialEnum(void){                     return          440; }
     82int MaticeEnum(void){                       return          441; }
     83int MatparEnum(void){                       return          442; }
     84int NumparEnum(void){                       return          443; }
    8485/*Inputs: */
    85 int InputEnum(void){                    return          450; }
     86int InputEnum(void){                        return          450; }
    8687/*Params: */
    87 int ParamEnum(void){                    return          460; }
     88int ParamEnum(void){                        return          460; }
    8889/*Results: */
    89 int ResultEnum(void){                   return          470; }
     90int ResultEnum(void){                       return          470; }
    9091/*Rgb: */
    91 int RgbEnum(void){                      return          480; }
     92int RgbEnum(void){                          return          480; }
    9293/*Spc: */
    93 int SpcEnum(void){                      return          490; }
     94int SpcEnum(void){                          return          490; }
    9495/*DofVec: */
    95 int DofVecEnum(void){                   return          495; }
     96int DofVecEnum(void){                       return          495; }
    9697
    9798
    9899/*GEOGRAPHY:*/
    99 int GeographyEnum(void){                return          500; }
    100 int IceSheetEnum(void){                 return          502; }
    101 int IceShelfEnum(void){                 return          502; }
     100int GeographyEnum(void){                    return          500; }
     101int IceSheetEnum(void){                     return          502; }
     102int IceShelfEnum(void){                     return          502; }
    102103
    103104/*FILL:*/
    104 int WaterEnum(void){                    return          601; }
    105 int IceEnum(void){                      return          602; }
    106 int AirEnum(void){                      return          603; }
    107 int MelangeEnum(void){                  return          604; }
     105int WaterEnum(void){                        return          601; }
     106int IceEnum(void){                          return          602; }
     107int AirEnum(void){                          return          603; }
     108int MelangeEnum(void){                      return          604; }
    108109
    109110/*functions on enums: */

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

    r2333 r2714  
    4242//prognostic
    4343int PrognosticAnalysisEnum(void);
     44int BalancedthicknessAnalysisEnum(void);
    4445//melting
    4546int MeltingAnalysisEnum(void);

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

    r2591 r2714  
    218218                                        ./ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp\
    219219                                        ./ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp\
     220                                        ./ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp\
     221                                        ./ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp\
     222                                        ./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\
     223                                        ./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\
    220224                                        ./ModelProcessorx/Qmu/CreateParametersQmu.cpp\
    221225                                        ./Dofx/Dofx.h\
    … …  
    517521                                        ./ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp\
    518522                                        ./ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp\
     523                                        ./ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp\
     524                                        ./ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp\
     525                                        ./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\
     526                                        ./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\
    519527                                        ./ModelProcessorx/Qmu/CreateParametersQmu.cpp\
    520528                                        ./Dofx/Dofx.h\
    … …  
    616624                                        ./parallel/ProcessResults.cpp\
    617625                                        ./parallel/prognostic_core.cpp\
     626                                        ./parallel/balancedthickness_core.cpp\
    618627                                        ./parallel/transient_core.cpp\
    619628                                        ./parallel/transient_core_2d.cpp\
    … …  
    627636bin_PROGRAMS =
    628637else
    629 bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe transient.exe steadystate.exe
     638bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe balancedthickness.exe transient.exe steadystate.exe
    630639endif
    631640
    … …  
    644653prognostic_exe_CXXFLAGS= -fPIC -D_PARALLEL_
    645654
     655balancedthickness_exe_SOURCES = parallel/balancedthickness.cpp
     656balancedthickness_exe_CXXFLAGS= -fPIC -D_PARALLEL_
     657
    646658transient_exe_SOURCES = parallel/transient.cpp
    647659transient_exe_CXXFLAGS= -fPIC -D_PARALLEL_
    648 

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

    r2112 r2714  
    8686                                       
    8787        }
     88        else if (iomodel->analysis_type==BalancedthicknessAnalysisEnum()){
     89
     90                CreateElementsNodesAndMaterialsBalancedthickness(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
     91                CreateConstraintsBalancedthickness(pconstraints,iomodel,iomodel_handle);
     92                CreateLoadsBalancedthickness(ploads,iomodel,iomodel_handle);
     93                CreateParametersBalancedthickness(pparameters,iomodel,iomodel_handle);
     94
     95        }
    8896        else{
     97
    8998                throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s"," analysis_type: ",iomodel->analysis_type," sub_analysis_type: ",iomodel->sub_analysis_type," not supported yet!"));
     99
    90100        }
    91101        CreateParametersQmu(pparameters,iomodel,iomodel_handle);

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

    r2586 r2714  
    247247        void    CreateParametersPrognostic(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
    248248
     249        /*balancedthickness:*/
     250        void    CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
     251        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);
     254
    249255        /*qmu: */
    250256        void CreateParametersQmu(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);

  • issm/trunk/src/c/io/FetchParams.cpp

    r2337 r2714  
    6666
    6767                        if (M==0 | N==0){
    68                                 throw ErrorException(__FUNCT__,exprintf("%s%s%i%s%i%s%i%s",__FUNCT__," error message: array in parameters structure field ",count," is of size (",M,",",N,")"));
     68                                throw ErrorException(__FUNCT__,exprintf("%s%s%i (%s) %s%i%s%i%s",__FUNCT__," error message: array in parameters structure field ",count,name," is of size (",M,",",N,")"));
    6969                        }
    7070

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

    r2713 r2714  
    321321                CreateKMatrixPrognostic( Kgg,inputs,analysis_type,sub_analysis_type);
    322322        }
     323        else if (analysis_type==BalancedthicknessAnalysisEnum()){
     324
     325                CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type);
     326        }
    323327        else if (analysis_type==ThermalAnalysisEnum()){
    324328
    … …  
    332336                throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
    333337        }
     338
     339}
     340/*}}}*/
     341/*FUNCTION CreateKMatrixBalancedthickness {{{1*/
     342#undef __FUNCT__
     343#define __FUNCT__ "Penta::CreateKMatrixBalancedthickness"
     344
     345void  Penta::CreateKMatrixBalancedthickness(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
     346
     347        /*Collapsed formulation: */
     348        Tria*  tria=NULL;
     349
     350        /*If on water, skip: */
     351        if(onwater)return;
     352
     353        /*Is this element on the bed? :*/
     354        if(!onbed)return;
     355
     356        /*Spawn Tria element from the base of the Penta: */
     357        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     358        tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
     359        delete tria;
     360        return;
    334361
    335362}
    … …  
    13791406                CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
    13801407        }
     1408        else if (analysis_type==BalancedthicknessAnalysisEnum()){
     1409
     1410                CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type);
     1411        }
    13811412        else if (analysis_type==ThermalAnalysisEnum()){
    13821413
    … …  
    13911422        }       
    13921423
     1424}
     1425/*}}}*/
     1426/*FUNCTION CreatePVectorBalancedthickness {{{1*/
     1427#undef __FUNCT__
     1428#define __FUNCT__ "Penta::CreatePVectorBalancedthickness"
     1429
     1430void Penta::CreatePVectorBalancedthickness( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
     1431
     1432        /*Collapsed formulation: */
     1433        Tria*  tria=NULL;
     1434
     1435        /*If on water, skip: */
     1436        if(onwater)return;
     1437
     1438        /*Is this element on the bed? :*/
     1439        if(!onbed)return;
     1440
     1441        /*Spawn Tria element from the base of the Penta: */
     1442        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     1443        tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
     1444        delete tria;
     1445        return;
    13931446}
    13941447/*}}}*/

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

    r2333 r2714  
    116116                void  CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    117117                void  CreatePVectorPrognostic( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
     118                void  CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
     119                void  CreatePVectorBalancedthickness( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
    118120
    119121                void  CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type);

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

    r2713 r2714  
    298298
    299299        }
     300        else if (analysis_type==BalancedthicknessAnalysisEnum()){
     301
     302                CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type);
     303
     304        }
    300305        else{
    301306                throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
    302307        }
     308
     309}
     310/*}}}*/
     311/*FUNCTION CreateKMatrixBalancedthickness {{{1*/
     312#undef __FUNCT__
     313#define __FUNCT__ "Tria::CreateKMatrixBalancedthickness"
     314void  Tria::CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
     315
     316        /* local declarations */
     317        int             i,j;
     318
     319        /* node data: */
     320        const int    numgrids=3;
     321        const int    NDOF1=1;
     322        const int    numdof=NDOF1*numgrids;
     323        double       xyz_list[numgrids][3];
     324        int          doflist[numdof];
     325        int          numberofdofspernode;
     326
     327        /* gaussian points: */
     328        int     num_gauss,ig;
     329        double* first_gauss_area_coord  =  NULL;
     330        double* second_gauss_area_coord =  NULL;
     331        double* third_gauss_area_coord  =  NULL;
     332        double* gauss_weights           =  NULL;
     333        double  gauss_weight;
     334        double  gauss_l1l2l3[3];
     335
     336        /* matrices: */
     337        double L[numgrids];
     338        double B[2][numgrids];
     339        double Bprime[2][numgrids];
     340        double DL[2][2]={0.0};
     341        double DLprime[2][2]={0.0};
     342        double DL_scalar;
     343        double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix
     344        double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
     345        double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
     346        double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
     347
     348        double Jdettria;
     349
     350        /*input parameters for structural analysis (diagnostic): */
     351        double  vxvy_list[numgrids][2]={0.0};
     352        double  vx_list[numgrids]={0.0};
     353        double  vy_list[numgrids]={0.0};
     354        double  dvx[2];
     355        double  dvy[2];
     356        double  vx,vy;
     357        double  dvxdx,dvydy;
     358        double  v_gauss[2]={0.0};
     359        double  K[2][2]={0.0};
     360        double  KDL[2][2]={0.0};
     361        int     dofs[2]={0,1};
     362        int     found=0;
     363
     364        ParameterInputs* inputs=NULL;
     365
     366        /*recover pointers: */
     367        inputs=(ParameterInputs*)vinputs;
     368
     369        /*recover extra inputs from users, at current convergence iteration: */
     370        found=inputs->Recover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
     371        if(!found)throw ErrorException(__FUNCT__," could not find velocity_average  in inputs!");
     372
     373        for(i=0;i<numgrids;i++){
     374                vx_list[i]=vxvy_list[i][0];
     375                vy_list[i]=vxvy_list[i][1];
     376        }
     377
     378        /* Get node coordinates and dof list: */
     379        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     380        GetDofList(&doflist[0],&numberofdofspernode);
     381
     382        //Create Artificial diffusivity once for all if requested
     383        if(numpar->artdiff){
     384                //Get the Jacobian determinant
     385                gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0;
     386                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     387
     388                //Build K matrix (artificial diffusivity matrix)
     389                v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
     390                v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
     391
     392                K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]);
     393                K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]);
     394        }
     395
     396        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     397        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     398
     399        /* Start  looping on the number of gaussian points: */
     400        for (ig=0; ig<num_gauss; ig++){
     401                /*Pick up the gaussian point: */
     402                gauss_weight=*(gauss_weights+ig);
     403                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     404                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     405                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     406
     407                /* Get Jacobian determinant: */
     408                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     409
     410                /*Get B  and B prime matrix: */
     411                GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
     412                GetBPrime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
     413
     414                //Get vx, vy and their derivatives at gauss point
     415                GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
     416                GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
     417
     418                GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3);
     419                GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3);
     420
     421                dvxdx=dvx[0];
     422                dvydy=dvy[1];
     423
     424                DL_scalar=gauss_weight*Jdettria;
     425
     426                //Create DL and DLprime matrix
     427                DL[0][0]=DL_scalar*dvxdx;
     428                DL[1][1]=DL_scalar*dvydy;
     429
     430                DLprime[0][0]=DL_scalar*vx;
     431                DLprime[1][1]=DL_scalar*vy;
     432
     433                //Do the triple product tL*D*L.
     434                //Ke_gg_thickness=B'*DLprime*Bprime;
     435
     436                TripleMultiply( &B[0][0],2,numdof,1,
     437                                        &DL[0][0],2,2,0,
     438                                        &B[0][0],2,numdof,0,
     439                                        &Ke_gg_thickness1[0][0],0);
     440
     441                TripleMultiply( &B[0][0],2,numdof,1,
     442                                        &DLprime[0][0],2,2,0,
     443                                        &Bprime[0][0],2,numdof,0,
     444                                        &Ke_gg_thickness2[0][0],0);
     445
     446                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     447                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness1[i][j];
     448                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
     449
     450                if(numpar->artdiff){
     451
     452                        /* Compute artificial diffusivity */
     453                        KDL[0][0]=DL_scalar*K[0][0];
     454                        KDL[1][1]=DL_scalar*K[1][1];
     455
     456                        TripleMultiply( &Bprime[0][0],2,numdof,1,
     457                                                &KDL[0][0],2,2,0,
     458                                                &Bprime[0][0],2,numdof,0,
     459                                                &Ke_gg_gaussian[0][0],0);
     460
     461                        /* Add artificial diffusivity matrix */
     462                        for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     463
     464                }
     465
     466#ifdef _DEBUGELEMENTS_
     467                if(my_rank==RANK && id==ELID){
     468                        printf("      B:\n");
     469                        for(i=0;i<3;i++){
     470                                for(j=0;j<numdof;j++){
     471                                        printf("%g ",B[i][j]);
     472                                }
     473                                printf("\n");
     474                        }
     475                        printf("      Bprime:\n");
     476                        for(i=0;i<3;i++){
     477                                for(j=0;j<numdof;j++){
     478                                        printf("%g ",Bprime[i][j]);
     479                                }
     480                                printf("\n");
     481                        }
     482                }
     483#endif
     484        } // for (ig=0; ig<num_gauss; ig++)
     485
     486        /*Add Ke_gg to global matrix Kgg: */
     487        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     488
     489#ifdef _DEBUGELEMENTS_
     490        if(my_rank==RANK && id==ELID){
     491                printf("      Ke_gg erms:\n");
     492                for( i=0; i<numdof; i++){
     493                        for (j=0;j<numdof;j++){
     494                                printf("%g ",Ke_gg[i][j]);
     495                        }
     496                        printf("\n");
     497                }
     498                printf("      Ke_gg row_indices:\n");
     499                for( i=0; i<numdof; i++){
     500                        printf("%i ",doflist[i]);
     501                }
     502
     503        }
     504#endif
     505
     506cleanup_and_return:
     507        xfree((void**)&first_gauss_area_coord);
     508        xfree((void**)&second_gauss_area_coord);
     509        xfree((void**)&third_gauss_area_coord);
     510        xfree((void**)&gauss_weights);
    303511
    304512}
    … …  
    13381546                CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
    13391547        }
     1548        else if (analysis_type==BalancedthicknessAnalysisEnum()){
     1549
     1550                CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type);
     1551        }
    13401552        else{
    13411553                throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis ",analysis_type," not supported yet"));
    13421554        }
     1555
     1556}
     1557/*}}}*/
     1558/*FUNCTION CreatePVectorBalancedthickness {{{1*/
     1559#undef __FUNCT__
     1560#define __FUNCT__ "Tria::CreatePVectorBalancedthickness"
     1561void  Tria::CreatePVectorBalancedthickness(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){
     1562
     1563
     1564        /* local declarations */
     1565        int             i,j;
     1566
     1567        /* node data: */
     1568        const int    numgrids=3;
     1569        const int    NDOF1=1;
     1570        const int    numdof=NDOF1*numgrids;
     1571        double       xyz_list[numgrids][3];
     1572        int          doflist[numdof];
     1573        int          numberofdofspernode;
     1574
     1575        /* gaussian points: */
     1576        int     num_gauss,ig;
     1577        double* first_gauss_area_coord  =  NULL;
     1578        double* second_gauss_area_coord =  NULL;
     1579        double* third_gauss_area_coord  =  NULL;
     1580        double* gauss_weights           =  NULL;
     1581        double  gauss_weight;
     1582        double  gauss_l1l2l3[3];
     1583
     1584        /* matrix */
     1585        double pe_g[numgrids]={0.0};
     1586        double L[numgrids];
     1587        double Jdettria;
     1588
     1589        /*input parameters for structural analysis (diagnostic): */
     1590        double  accumulation_list[numgrids]={0.0};
     1591        double  accumulation_g;
     1592        double  melting_list[numgrids]={0.0};
     1593        double  melting_g;
     1594        double  thickness_list[numgrids]={0.0};
     1595        double  thickness_g;
     1596        int     dofs[1]={0};
     1597        int     found=0;
     1598
     1599        ParameterInputs* inputs=NULL;
     1600
     1601        /*recover pointers: */
     1602        inputs=(ParameterInputs*)vinputs;
     1603
     1604        /*recover extra inputs from users, at current convergence iteration: */
     1605        found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes);
     1606        if(!found)throw ErrorException(__FUNCT__," could not find accumulation in inputs!");
     1607        found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes);
     1608        if(!found)throw ErrorException(__FUNCT__," could not find melting in inputs!");
     1609
     1610        /* Get node coordinates and dof list: */
     1611        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     1612        GetDofList(&doflist[0],&numberofdofspernode);
     1613
     1614        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     1615        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     1616
     1617        /* Start  looping on the number of gaussian points: */
     1618        for (ig=0; ig<num_gauss; ig++){
     1619                /*Pick up the gaussian point: */
     1620                gauss_weight=*(gauss_weights+ig);
     1621                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     1622                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     1623                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     1624
     1625                /* Get Jacobian determinant: */
     1626                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     1627
     1628                /*Get L matrix: */
     1629                GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     1630
     1631                /* Get accumulation, melting and thickness at gauss point */
     1632                GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3);
     1633                GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3);
     1634
     1635                /* Add value into pe_g: */
     1636                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(accumulation_g-melting_g)*L[i];
     1637
     1638        } // for (ig=0; ig<num_gauss; ig++)
     1639
     1640        /*Add pe_g to global matrix Kgg: */
     1641        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
     1642
     1643cleanup_and_return:
     1644        xfree((void**)&first_gauss_area_coord);
     1645        xfree((void**)&second_gauss_area_coord);
     1646        xfree((void**)&third_gauss_area_coord);
     1647        xfree((void**)&gauss_weights);
    13431648
    13441649}

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

    r2334 r2714  
    123123                void  CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    124124                void  CreatePVectorPrognostic(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);
     125                void  CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
     126                void  CreatePVectorBalancedthickness(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);
    125127                double MassFlux(double* segment,double* ug);
    126128                double GetArea(void);
    127129                double GetAreaCoordinate(double x, double y, int which_one);
    128130
    129 
    130131};
    131132#endif  /* _TRIA_H */

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

    r2333 r2714  
    109109        fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum(),HutterAnalysisEnum());
    110110        fem_sl=model->GetFormulation(SlopeComputeAnalysisEnum());
    111         fem_p=model->GetFormulation(PrognosticAnalysisEnum());
     111        if(analysis_type==PrognosticAnalysisEnum()){
     112                fem_p=model->GetFormulation(PrognosticAnalysisEnum());
     113        }
     114        if(analysis_type==BalancedthicknessAnalysisEnum()){
     115                fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum());
     116        }
    112117        fem_t=model->GetFormulation(ThermalAnalysisEnum());
    113118

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

    r2383 r2714  
    1717void diagnostic_core(DataSet* results,Model* model, ParameterInputs* inputs);
    1818void prognostic_core(DataSet* results,Model* model, ParameterInputs* inputs);
     19void balancedthickness_core(DataSet* results,Model* model, ParameterInputs* inputs);
    1920void control_core(DataSet* results,Model* model, ParameterInputs* inputs);
    2021

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

    r2629 r2714  
    2828        Model* model=NULL;
    2929
    30         Vec     h_g=NULL;
    3130        Vec     u_g=NULL;
    3231        double* u_g_serial=NULL;

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

    r1120 r2714  
    5656                numdofs=1;
    5757        }
     58        else if (analysis_type==BalancedthicknessAnalysisEnum()){
     59                numdofs=1;
     60        }
    5861        else throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis type: ",analysis_type,"  not implemented yet!"));
    5962

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

    r2686 r2714  
    310310        checknan(md,fields);
    311311
    312         %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
     312        %INITIAL TEMPERATURE
    313313        fields={'temperature','spctemperature(:,2)','observed_temperature'};
    314314        checkgreater(md,fields,0)
    315315
     316end
     317
     318%BALANCEDTHICKNESS
     319if md.analysis_type==BalancedthicknessAnalysisEnum
     320
     321        %VELOCITIES MELTING AND ACCUMULATION
     322        fields={'vx','vy','accumulation','melting'};
     323        checksize(md,fields,[md.numberofgrids 1]);
     324        checknan(md,fields);
     325
     326        %SPC
     327        if any(md.spcthickness(find(md.gridonboundary))~=1),
     328                error(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcthickness']);
     329        end
    316330end
    317331

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

    r2547 r2714  
    1717
    1818%check solution type is supported
    19 if ~ismemberi(analysis_type,{'diagnostic','prognostic','thermal','steadystate','parameters','transient'}),
     19if ~ismemberi(analysis_type,{'diagnostic','prognostic','thermal','steadystate','parameters','transient','balancedthickness'}),
    2020        error(['process_solve_options error message: analysis_type ' analysis_type ' not supported yet!']);
    2121else

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

    r2540 r2714  
    7171        md=steadystate(md);
    7272
     73elseif md.analysis_type==BalancedthicknessAnalysisEnum,
     74        md=balancedthickness(md);
    7375else
    7476        error('solution type not supported for this model configuration.');

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

    r1901 r2714  
    8585end
    8686
     87if enum==BalancedthicknessAnalysisEnum(),
     88        string='balancedthickness';
     89end
     90
    8791if enum==MeltingAnalysisEnum(),
    8892        string='melting';

  • issm/trunk/src/pro/bytscl.m

    r1 r2714  
    33%Equivalent of the bytscl idl routine.
    44
    5 Min=min(min(value));
    6 Max=max(max(value));
     5Min=min(value(:));
     6Max=max(value(:));
    77Top=255;
    88
    9 value=(Top + 0.9999)*(value - Min)/(Max - Min);
     9pow=1;
     10value=(Top + 0.9999)*(value - Min).^(pow)/(Max - Min)^(pow);
Note: See TracChangeset for help on using the changeset viewer.