Changeset 358


Ignore:
Timestamp:
05/12/09 14:01:34 (16 years ago)
Author:
Eric.Larour
Message:

New slope extrude module. New slope computation in Penta and Tria

Location:
issm/trunk/src
Files:
7 added
28 edited

Legend:

Unmodified
Added
Removed

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

    r343 r358  
    11531153}
    11541154
     1155void  DataSet::SlopeExtrude(Vec sg,double* sg_serial){
     1156
     1157        vector<Object*>::iterator object;
     1158        Penta* penta=NULL;
     1159
     1160        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1161               
     1162                if((*object)->Enum()==PentaEnum()){
     1163
     1164                        penta=(Penta*)(*object);
     1165                        penta->SlopeExtrude(sg,sg_serial);
     1166
     1167                }
     1168        }
     1169
     1170}
    11551171
    11561172void DataSet::ComputePressure(Vec p_g){

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

    r304 r358  
    7575                void  Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type);
    7676                void  VelocityExtrude(Vec ug,double* ug_serial);
     77                void  SlopeExtrude(Vec sg,double* sg_serial);
    7778                int   DeleteObject(Object* object);
    7879                void  ComputePressure(Vec p_g);

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

    r304 r358  
    4242                return DiagnosticHutterAnalysisEnum();
    4343        }
    44         else if (strcmp(analysis_type,"surface_slope_compute")==0){
    45                 return SurfaceSlopeComputeAnalysisEnum();
     44        else if (strcmp(analysis_type,"slope_compute")==0){
     45                return SlopeComputeAnalysisEnum();
    4646        }
    47         else if (strcmp(analysis_type,"bed_slope_compute")==0){
    48                 return BedSlopeComputeAnalysisEnum();
     47        else if (strcmp(analysis_type,"surface_slope_compute_x")==0){
     48                return SurfaceSlopeComputeXAnalysisEnum();
     49        }
     50        else if (strcmp(analysis_type,"surface_slope_compute_y")==0){
     51                return SurfaceSlopeComputeYAnalysisEnum();
     52        }
     53        else if (strcmp(analysis_type,"bed_slope_compute_x")==0){
     54                return BedSlopeComputeXAnalysisEnum();
     55        }
     56        else if (strcmp(analysis_type,"bed_slope_compute_y")==0){
     57                return BedSlopeComputeYAnalysisEnum();
    4958        }
    5059        else throw ErrorException(__FUNCT__," could not recognized analysis type!");

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

    r308 r358  
    2626int DiagnosticStokesAnalysisEnum(void){ return          30; }
    2727int DiagnosticHutterAnalysisEnum(void){ return          31; }
    28 int SurfaceSlopeComputeAnalysisEnum(void){ return       32; }
    29 int BedSlopeComputeAnalysisEnum(void){ return           33; }
     28int SlopeComputeAnalysisEnum(void){ return              32; }
     29int SurfaceSlopeComputeXAnalysisEnum(void){ return      33; }
     30int SurfaceSlopeComputeYAnalysisEnum(void){ return      34; }
     31int BedSlopeComputeXAnalysisEnum(void){ return          35; }
     32int BedSlopeComputeYAnalysisEnum(void){ return          36; }
    3033
    3134/*datasets: */

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

    r308 r358  
    4242int DiagnosticStokesAnalysisEnum(void);
    4343int DiagnosticHutterAnalysisEnum(void);
    44 int SurfaceSlopeComputeAnalysisEnum(void);
    45 int BedSlopeComputeAnalysisEnum(void);
     44int SlopeComputeAnalysisEnum(void);
     45int SurfaceSlopeComputeXAnalysisEnum(void);
     46int SurfaceSlopeComputeYAnalysisEnum(void);
     47int BedSlopeComputeXAnalysisEnum(void);
     48int BedSlopeComputeYAnalysisEnum(void);
    4649
    4750

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

    r308 r358  
    179179                                        ./ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp \
    180180                                        ./ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp\
    181                                         ./ModelProcessorx/SurfaceSlopeCompute/CreateElementsNodesAndMaterialsSurfaceSlopeCompute.cpp\
    182                                         ./ModelProcessorx/SurfaceSlopeCompute/CreateConstraintsSurfaceSlopeCompute.cpp \
    183                                         ./ModelProcessorx/SurfaceSlopeCompute/CreateLoadsSurfaceSlopeCompute.cpp\
    184                                         ./ModelProcessorx/BedSlopeCompute/CreateElementsNodesAndMaterialsBedSlopeCompute.cpp\
    185                                         ./ModelProcessorx/BedSlopeCompute/CreateConstraintsBedSlopeCompute.cpp \
    186                                         ./ModelProcessorx/BedSlopeCompute/CreateLoadsBedSlopeCompute.cpp\
     181                                        ./ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp\
     182                                        ./ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp \
     183                                        ./ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp\
    187184                                        ./ModelProcessorx/Control/CreateParametersControl.cpp\
    188185                                        ./Dofx/Dofx.h\
    … …  
    243240                                        ./ProcessParamsx/ProcessParamsx.h\
    244241                                        ./VelocityExtrudex/VelocityExtrudex.cpp\
    245                                         ./VelocityExtrudex/VelocityExtrudex.h
     242                                        ./VelocityExtrudex/VelocityExtrudex.h\
     243                                        ./SlopeExtrudex/SlopeExtrudex.cpp\
     244                                        ./SlopeExtrudex/SlopeExtrudex.h
     245
    246246
    247247
    … …  
    420420                                        ./ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp \
    421421                                        ./ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp\
    422                                         ./ModelProcessorx/SurfaceSlopeCompute/CreateElementsNodesAndMaterialsSurfaceSlopeCompute.cpp\
    423                                         ./ModelProcessorx/SurfaceSlopeCompute/CreateConstraintsSurfaceSlopeCompute.cpp \
    424                                         ./ModelProcessorx/SurfaceSlopeCompute/CreateLoadsSurfaceSlopeCompute.cpp\
    425                                         ./ModelProcessorx/BedSlopeCompute/CreateElementsNodesAndMaterialsBedSlopeCompute.cpp\
    426                                         ./ModelProcessorx/BedSlopeCompute/CreateConstraintsBedSlopeCompute.cpp \
    427                                         ./ModelProcessorx/BedSlopeCompute/CreateLoadsBedSlopeCompute.cpp\
     422                                        ./ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp\
     423                                        ./ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp \
     424                                        ./ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp\
    428425                                        ./ModelProcessorx/Control/CreateParametersControl.cpp\
    429426                                        ./Dofx/Dofx.h\
    … …  
    485482                                        ./VelocityExtrudex/VelocityExtrudex.cpp\
    486483                                        ./VelocityExtrudex/VelocityExtrudex.h\
     484                                        ./SlopeExtrudex/SlopeExtrudex.cpp\
     485                                        ./SlopeExtrudex/SlopeExtrudex.h\
    487486                                        ./parallel/diagnostic_core_nonlinear.cpp\
    488487                                        ./parallel/CreateFemModel.cpp\

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

    r300 r358  
    3030                CreateConstraintsDiagnosticHutter(pconstraints,model,model_handle);
    3131        }
    32         else if (strcmp(model->analysis_type,"surface_slope_compute")==0){
    33                 CreateConstraintsSurfaceSlopeCompute(pconstraints,model,model_handle);
    34         }
    35         else if (strcmp(model->analysis_type,"bed_slope_compute")==0){
    36                 CreateConstraintsBedSlopeCompute(pconstraints,model,model_handle);
     32        else if (strcmp(model->analysis_type,"slope_compute")==0){
     33                CreateConstraintsSlopeCompute(pconstraints,model,model_handle);
    3734        }
    3835        /*

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

    r300 r358  
    3939       
    4040        }
    41         else if ((strcmp(model->analysis_type,"surface_slope_compute")==0)){
     41        else if ((strcmp(model->analysis_type,"slope_compute")==0)){
    4242
    43                 CreateElementsNodesAndMaterialsSurfaceSlopeCompute(pelements,pnodes,pmaterials, model,model_handle);
    44        
    45         }
    46         else if ((strcmp(model->analysis_type,"bed_slope_compute")==0)){
    47 
    48                 CreateElementsNodesAndMaterialsBedSlopeCompute(pelements,pnodes,pmaterials, model,model_handle);
     43                CreateElementsNodesAndMaterialsSlopeCompute(pelements,pnodes,pmaterials, model,model_handle);
    4944       
    5045        }

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

    r300 r358  
    3434                CreateLoadsDiagnosticHutter(ploads,model,model_handle);
    3535        }
    36         else if (strcmp(model->analysis_type,"surface_slope_compute")==0){
     36        else if (strcmp(model->analysis_type,"slope_compute")==0){
    3737
    38                 CreateLoadsSurfaceSlopeCompute(ploads,model,model_handle);
    39         }
    40         else if (strcmp(model->analysis_type,"bed_slope_compute")==0){
    41 
    42                 CreateLoadsBedSlopeCompute(ploads,model,model_handle);
     38                CreateLoadsSlopeCompute(ploads,model,model_handle);
    4339        }
    4440        /*

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

    r347 r358  
    293293        ModelFetchData((void**)&model->ishutter,NULL,NULL,model_handle,"ishutter","Integer",NULL);
    294294        ModelFetchData((void**)&model->ismacayealpattyn,NULL,NULL,model_handle,"ismacayealpattyn","Integer",NULL);
    295         ModelFetchData((void**)&model->isstokes,NULL,NULL,model_handle,"isstokes","Integer",NULL);
    296295
    297296        /*!Get drag_type, drag and p,q: */
    … …  
    340339        ModelFetchData((void**)&model->mixed_layer_capacity,NULL,NULL,model_handle,"mixed_layer_capacity","Scalar",NULL);
    341340        ModelFetchData((void**)&model->thermal_exchange_velocity,NULL,NULL,model_handle,"thermal_exchange_velocity","Scalar",NULL);
    342        
     341
    343342        /*rifts: */
    344343        ModelFetchData((void**)&model->numrifts,NULL,NULL,model_handle,"numrifts","Integer",NULL);

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

    r308 r358  
    200200        void    CreateParametersDiagnosticStokes(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
    201201
    202         /*surface slope compute*/
    203         void    CreateElementsNodesAndMaterialsSurfaceSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle);
    204         void    CreateConstraintsSurfaceSlopeCompute(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
    205         void    CreateLoadsSurfaceSlopeCompute(DataSet** ploads, Model* model, ConstDataHandle model_handle);
    206         void    CreateParametersSurfaceSlopeCompute(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
    207 
    208         /*bed slope compute*/
    209         void    CreateElementsNodesAndMaterialsBedSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle);
    210         void    CreateConstraintsBedSlopeCompute(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
    211         void    CreateLoadsBedSlopeCompute(DataSet** ploads, Model* model, ConstDataHandle model_handle);
    212         void    CreateParametersBedSlopeCompute(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
     202        /*slope compute*/
     203        void    CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle);
     204        void    CreateConstraintsSlopeCompute(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
     205        void    CreateLoadsSlopeCompute(DataSet** ploads, Model* model, ConstDataHandle model_handle);
     206        void    CreateParametersSlopeCompute(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
    213207
    214208        /*control:*/

  • issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp

    r348 r358  
    11/*
    2  * CreateConstraintsBedSlopeCompute.c:
     2 * CreateConstraintsSlopeCompute.c:
    33 */
    44
    55#undef __FUNCT__
    6 #define __FUNCT__ "CreateConstraintsBedSlopeCompute"
     6#define __FUNCT__ "CreateConstraintsSlopeCompute"
    77
    88#include "../../DataSet/DataSet.h"
    … …  
    1414
    1515
    16 void    CreateConstraintsBedSlopeCompute(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){
     16void    CreateConstraintsSlopeCompute(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){
    1717
    1818

  • issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp

    r348 r358  
    11/*
    2  * CreateElementsNodesAndMaterialsBedSlopeCompute.c:
     2 * CreateElementsNodesAndMaterialsSlopeCompute.c:
    33 */
    44
    55#undef __FUNCT__
    6 #define __FUNCT__ "CreateElementsNodesAndMaterialsBedSlopeCompute"
     6#define __FUNCT__ "CreateElementsNodesAndMaterialsSlopeCompute"
    77
    88#include "../../DataSet/DataSet.h"
    … …  
    1515
    1616
    17 void    CreateElementsNodesAndMaterialsBedSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){
    18 
     17void    CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){
     18
     19
     20                /*output: int* epart, int* my_grids, double* my_bordergrids*/
     21
     22
     23        int i,j,k,n;
     24        extern int my_rank;
     25        extern int num_procs;
    1926
    2027        /*DataSets: */
    … …  
    2330        DataSet*    materials = NULL;
    2431       
     32        /*Objects: */
     33        Node*       node   = NULL;
     34        Tria*       tria = NULL;
     35        Penta*      penta = NULL;
     36
     37        int         analysis_type;
     38       
     39        /*output: */
     40        int* epart=NULL; //element partitioning.
     41        int* npart=NULL; //node partitioning.
     42        int* my_grids=NULL;
     43        double* my_bordergrids=NULL;
     44
     45        /*intermediary: */
     46        int elements_width; //size of elements
     47                       
     48        /*tria constructor input: */
     49        int tria_id;
     50        int tria_mid;
     51        int tria_mparid;
     52        int tria_g[3];
     53        double tria_h[3];
     54        double tria_s[3];
     55        double tria_b[3];
     56        double tria_k[3];
     57        double tria_melting[3];
     58        double tria_accumulation[3];
     59        int    tria_friction_type;
     60        double tria_p;
     61        double tria_q;
     62        int    tria_shelf;
     63        double tria_meanvel;/*!scaling ratio for velocities*/
     64        double tria_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
     65        double tria_viscosity_overshoot;
     66
     67        /*penta constructor input: */
     68
     69        int penta_id;
     70        int penta_mid;
     71        int penta_mparid;
     72        int penta_g[6];
     73        double penta_h[6];
     74        double penta_s[6];
     75        double penta_b[6];
     76        double penta_k[6];
     77        int penta_friction_type;
     78        double penta_p;
     79        double penta_q;
     80        int penta_shelf;
     81        int penta_onbed;
     82        int penta_onsurface;
     83        double penta_meanvel;/*!scaling ratio for velocities*/
     84        double penta_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
     85        int penta_collapse;
     86        double penta_melting[6];
     87        double penta_accumulation[6];
     88        double penta_geothermalflux[6];
     89        int penta_artdiff;
     90        int penta_thermal_steadystate;
     91        double penta_viscosity_overshoot;
     92
     93        /* node constructor input: */
     94        int node_id;
     95        int node_partitionborder=0;
     96        double node_x[3];
     97        int node_onbed;
     98        int node_onsurface;
     99        int node_upper_node_id;
     100        int node_numdofs;
     101
     102               
     103        #ifdef _PARALLEL_
     104        /*Metis partitioning: */
     105        int  range;
     106        Vec  gridborder;
     107        int  my_numgrids;
     108        int* all_numgrids=NULL;
     109        int  gridcount;
     110        int  count;
     111        #endif
     112        int  first_grid_index;
     113
     114        /*Rifts:*/
     115        int* riftsnumpenaltypairs;
     116        double** riftspenaltypairs;
     117        int* riftsfill;
     118        double* riftsfriction;
     119        double* riftpenaltypairs=NULL;
     120        int el1,el2;
     121
    25122        /*First create the elements, nodes and material properties: */
    26123        elements  = new DataSet(ElementsEnum());
    … …  
    31128        if (!model->isstokes)goto cleanup_and_return;
    32129
     130        /*Get analysis_type: */
     131        analysis_type=AnalysisTypeAsEnum(model->analysis_type);
     132
     133        /*Width of elements: */
     134        if(strcmp(model->meshtype,"2d")==0){
     135                elements_width=3; //tria elements
     136        }
     137        else{
     138                elements_width=6; //penta elements
     139        }
     140
     141        #ifdef _PARALLEL_
     142        /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
     143        if(strcmp(model->meshtype,"2d")==0){
     144                /*load elements: */
     145                ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
     146        }
     147        else{
     148                /*load elements2d: */
     149                ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
     150        }
     151
     152
     153        MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs);
     154
     155        /*Free elements and elements2d: */
     156        xfree((void**)&model->elements);
     157        xfree((void**)&model->elements2d);
     158
     159
     160        /*Deal with rifts, they have to be included into one partition only, not several: */
     161        FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts);
     162       
     163        for(i=0;i<model->numrifts;i++){
     164                riftpenaltypairs=model->riftspenaltypairs[i];
     165                for(j=0;j<model->riftsnumpenaltypairs[i];j++){
     166                        el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing
     167                        el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing
     168                        epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids;
     169                }
     170        }
     171        /*Free rifts: */
     172        xfree((void**)&riftsnumpenaltypairs);
     173        for(i=0;i<model->numrifts;i++){
     174                double* temp=riftspenaltypairs[i];
     175                xfree((void**)&temp);
     176        }
     177        xfree((void**)&riftspenaltypairs);
     178        xfree((void**)&riftsfill);
     179        xfree((void**)&riftsfriction);
     180
     181        /*Used later on: */
     182        my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int));
     183        #endif
     184
     185
     186
     187        /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */
     188
     189        /*2d mesh: */
     190        if (strcmp(model->meshtype,"2d")==0){
     191
     192                /*Fetch data needed: */
     193                ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
     194                ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
     195                ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
     196               
     197                for (i=0;i<model->numberofelements;i++){
     198
     199                #ifdef _PARALLEL_
     200                /*!All elements have been partitioned above, only create elements for this CPU: */
     201                if(my_rank==epart[i]){
     202                #endif
     203                       
     204                       
     205                        /*ids: */
     206                        tria_id=i+1; //matlab indexing.
     207                        tria_mid=-1; //no materials
     208                        tria_mparid=-1; //no materials
     209
     210                        /*vertices offsets: */
     211                        tria_g[0]=(int)*(model->elements+elements_width*i+0);
     212                        tria_g[1]=(int)*(model->elements+elements_width*i+1);
     213                        tria_g[2]=(int)*(model->elements+elements_width*i+2);
     214
     215                        /*surface and bed:*/
     216                        tria_s[0]=*(model->surface+    ((int)*(model->elements+elements_width*i+0)-1));
     217                        tria_s[1]=*(model->surface+    ((int)*(model->elements+elements_width*i+1)-1));
     218                        tria_s[2]=*(model->surface+    ((int)*(model->elements+elements_width*i+2)-1));
     219
     220                        tria_b[0]=*(model->bed+        ((int)*(model->elements+elements_width*i+0)-1));
     221                        tria_b[1]=*(model->bed+        ((int)*(model->elements+elements_width*i+1)-1));
     222                        tria_b[2]=*(model->bed+        ((int)*(model->elements+elements_width*i+2)-1));
     223
     224                        /*Create tria element using its constructor:*/
     225                        tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot);
     226
     227                        /*Add tria element to elements dataset: */
     228                        elements->AddObject(tria);
     229
     230                        #ifdef _PARALLEL_
     231                        /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
     232                         *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
     233                         into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
     234                         will hold which grids belong to this partition*/
     235                        my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
     236                        my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
     237                        my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
     238                        #endif
     239
     240                #ifdef _PARALLEL_
     241                }//if(my_rank==epart[i])
     242                #endif
     243
     244                }//for (i=0;i<numberofelements;i++)
     245
     246       
     247                /*Free data : */
     248                xfree((void**)&model->elements);
     249                xfree((void**)&model->surface);
     250                xfree((void**)&model->bed);
     251
     252        }
     253        else{ //        if (strcmp(meshtype,"2d")==0)
     254
     255                /*Fetch data needed: */
     256                ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
     257                ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
     258                ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
     259                ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat");
     260               
     261                for (i=0;i<model->numberofelements;i++){
     262                #ifdef _PARALLEL_
     263                /*We are using our element partition to decide which elements will be created on this node: */
     264                if(my_rank==epart[i]){
     265                #endif
     266
     267                       
     268                        /*name and id: */
     269                        penta_id=i+1; //matlab indexing.
     270                        penta_mid=-1; //no materials
     271                        penta_mparid=-1; //no materials
     272                       
     273
     274                        /*vertices,thickness,surface,bed and drag: */
     275                        for(j=0;j<6;j++){
     276                                penta_g[j]=(int)*(model->elements+elements_width*i+j);
     277                                penta_s[j]=*(model->surface+    ((int)*(model->elements+elements_width*i+j)-1));
     278                                penta_b[j]=*(model->bed+    ((int)*(model->elements+elements_width*i+j)-1));
     279                        }
     280                       
     281                        /*diverse: */
     282                        penta_onbed=(int)*(model->elementonbed+i);
     283
     284                        /*Create Penta using its constructor:*/
     285                        penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
     286                                        penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
     287                                        penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
     288                                        penta_thermal_steadystate,penta_viscosity_overshoot);
     289
     290                        /*Add penta element to elements dataset: */
     291                        elements->AddObject(penta);
     292       
     293
     294                        #ifdef _PARALLEL_
     295                        /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
     296                         *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
     297                         into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
     298                         will hold which grids belong to this partition*/
     299                        my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
     300                        my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
     301                        my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
     302                        my_grids[(int)*(model->elements+elements_width*i+3)-1]=1;
     303                        my_grids[(int)*(model->elements+elements_width*i+4)-1]=1;
     304                        my_grids[(int)*(model->elements+elements_width*i+5)-1]=1;
     305                        #endif
     306
     307                #ifdef _PARALLEL_
     308                }//if(my_rank==epart[i])
     309                #endif
     310
     311                }//for (i=0;i<numberofelements;i++)
     312
     313                /*Free data: */
     314                xfree((void**)&model->elements);
     315                xfree((void**)&model->surface);
     316                xfree((void**)&model->bed);
     317                xfree((void**)&model->elementonbed);
     318
     319        } //if (strcmp(meshtype,"2d")==0)
     320
     321        #ifdef _PARALLEL_
     322                /*From the element partitioning, we can determine which grids are on the inside of this cpu's
     323                 *element partition, and which are on its border with other nodes:*/
     324                gridborder=NewVec(model->numberofnodes);
     325
     326                for (i=0;i<model->numberofnodes;i++){
     327                        if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
     328                }
     329                VecAssemblyBegin(gridborder);
     330                VecAssemblyEnd(gridborder);
     331
     332                #ifdef _DEBUG_
     333                VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
     334                #endif
     335               
     336                VecToMPISerial(&my_bordergrids,gridborder);
     337
     338                #ifdef _DEBUG_
     339                if(my_rank==0){
     340                        for (i=0;i<model->numberofnodes;i++){
     341                                printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
     342                        }
     343                }
     344                #endif
     345        #endif
     346
     347        /*Partition penalties in 3d: */
     348        if(strcmp(model->meshtype,"3d")==0){
     349       
     350                /*Get penalties: */
     351                ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
     352
     353                if(model->numpenalties){
     354
     355                        model->penaltypartitioning=(int*)xmalloc(model->numpenalties*sizeof(int));
     356                        #ifdef _SERIAL_
     357                        for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=1;
     358                        #else
     359                        for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=-1;
     360
     361                        for(i=0;i<model->numpenalties;i++){
     362                                first_grid_index=(int)(*(model->penalties+i*model->numlayers+0)-1);
     363                                if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition  grids
     364                                        /*All grids that are being penalised belong to this node's internal grid partition.:*/
     365                                        model->penaltypartitioning[i]=1;
     366                                }
     367                                if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
     368                                        model->penaltypartitioning[i]=0;
     369                                }
     370                        }
     371                        #endif
     372                }
     373
     374                /*Free penalties: */
     375                xfree((void**)&model->penalties);
     376        }
     377
     378        /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids.
     379         We can therefore determine  which grids are internal to this node's partition
     380         and which ones are shared with other nodes because they are on the border of this node's partition. Knowing
     381         that, go and create the grids*/
     382
     383        /*Create nodes from x,y,z, as well as the spc values on those grids: */
     384               
     385        /*First fetch data: */
     386        if (strcmp(model->meshtype,"3d")==0){
     387                ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat");
     388                ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat");
     389        }
     390        ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat");
     391        ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat");
     392        ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat");
     393        ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
     394        ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
     395
     396
     397        /*Get number of dofs per node: */
     398        DistributeNumDofs(&node_numdofs,analysis_type);
     399
     400        for (i=0;i<model->numberofnodes;i++){
     401        #ifdef _PARALLEL_
     402        /*keep only this partition's nodes:*/
     403        if((my_grids[i]==1)){
     404        #endif
     405
     406                node_id=i+1; //matlab indexing
     407                       
     408               
     409               
     410                #ifdef _PARALLEL_
     411                if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
     412                        node_partitionborder=1;
     413                }
     414                else{
     415                        node_partitionborder=0;
     416                }
     417                #else
     418                        node_partitionborder=0;
     419                #endif
     420
     421
     422                node_x[0]=model->x[i];
     423                node_x[1]=model->y[i];
     424                node_x[2]=model->z[i];
     425
     426               
     427                node_onbed=(int)model->gridonbed[i];
     428                node_onsurface=(int)model->gridonsurface[i];   
     429                if (strcmp(model->meshtype,"3d")==0){
     430                        if (isnan(model->uppernodes[i])){
     431                                node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
     432                        }
     433                        else{
     434                                node_upper_node_id=(int)model->uppernodes[i];
     435                        }
     436                }
     437                else{
     438                        /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
     439                        node_upper_node_id=node_id;
     440                }
     441
     442                /*Create node using its constructor: */
     443                node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
     444
     445                /*set single point constraints.: */
     446                if (strcmp(model->meshtype,"3d")==0){
     447                        /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
     448                        if (model->deadgrids[i]){
     449                                for(k=1;k<=node_numdofs;k++){
     450                                        node->FreezeDof(k);
     451                                }
     452                        }
     453                }
     454                /*Add node to nodes dataset: */
     455                nodes->AddObject(node);
     456
     457        #ifdef _PARALLEL_
     458        } //if((my_grids[i]==1))
     459        #endif
     460        }
     461
     462        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
     463         * datasets, it will not be redone: */
     464        elements->Presort();
     465        nodes->Presort();
     466        materials->Presort();
     467
     468        /*Clean fetched data: */
     469        xfree((void**)&model->deadgrids);
     470        xfree((void**)&model->x);
     471        xfree((void**)&model->y);
     472        xfree((void**)&model->z);
     473        xfree((void**)&model->gridonbed);
     474        xfree((void**)&model->gridonsurface);
     475        xfree((void**)&model->uppernodes);
     476               
    33477        cleanup_and_return:
    34478
    … …  
    37481        *pnodes=nodes;
    38482        *pmaterials=materials;
     483
     484        /*Keep partitioning information into model*/
     485        model->epart=epart;
     486        model->my_grids=my_grids;
     487        model->my_bordergrids=my_bordergrids;
     488
     489        /*Free ressources:*/
     490        #ifdef _PARALLEL_
     491        xfree((void**)&all_numgrids);
     492        xfree((void**)&npart);
     493        VecFree(&gridborder);
     494        #endif
    39495}

  • issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp

    r348 r358  
    1 /*! \file CreateLoadsBedSlopeCompute.c:
     1/*! \file CreateLoadsSlopeCompute.c:
    22 */
    33
    44#undef __FUNCT__
    5 #define __FUNCT__ "CreateLoadsBedSlopeCompute"
     5#define __FUNCT__ "CreateLoadsSlopeCompute"
    66
    77#include "../../DataSet/DataSet.h"
    … …  
    1414
    1515
    16 void    CreateLoadsBedSlopeCompute(DataSet** ploads, Model* model,ConstDataHandle model_handle){
     16void    CreateLoadsSlopeCompute(DataSet** ploads, Model* model,ConstDataHandle model_handle){
    1717
    1818        DataSet*    loads    = NULL;

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

    r304 r358  
    4747#include "./GriddataMeshToGridx/GriddataMeshToGridx.h"
    4848#include "./ComputePressurex/ComputePressurex.h"
     49#include "./SlopeExtrudex/SlopeExtrudex.h"
    4950
    5051

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

    r350 r358  
    292292
    293293        }
     294        else if (
     295                        (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) ||
     296                        (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) ||
     297                        (analysis_type==BedSlopeComputeXAnalysisEnum()) ||
     298                        (analysis_type==BedSlopeComputeYAnalysisEnum())
     299                        ){
     300               
     301                CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type);
     302        }
    294303        else{
    295304                throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
    … …  
    703712                CreatePVectorDiagnosticVert( pg,inputs,analysis_type);
    704713        }
     714        else if (
     715                        (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) ||
     716                        (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) ||
     717                        (analysis_type==BedSlopeComputeXAnalysisEnum()) ||
     718                        (analysis_type==BedSlopeComputeYAnalysisEnum())
     719                        ){
     720               
     721                CreatePVectorSlopeCompute( pg,inputs,analysis_type);
     722        }
    705723        else{
    706724                throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
    … …  
    15751593
    15761594#undef __FUNCT__
     1595#define __FUNCT__ "Penta::SlopeExtrude"
     1596void  Penta::SlopeExtrude(Vec sg,double* sg_serial){
     1597
     1598        /* node data: */
     1599        const int    numgrids=6;
     1600        const int    NDOF1=1;
     1601        const int    numdofs=NDOF1*numgrids;
     1602        int          doflist[numdofs];
     1603        int          nodedof;
     1604        int          numberofdofspernode;
     1605       
     1606        Node* node=NULL;
     1607        int   i;
     1608        double slope;
     1609               
     1610
     1611        if(onbed==1){
     1612                       
     1613                GetDofList(&doflist[0],&numberofdofspernode);
     1614
     1615                if(numberofdofspernode!=1)throw ErrorException(__FUNCT__," slope can only be extruded on 1 dof per node");
     1616
     1617                /*For each node on the base of this penta,  we grab the slope. Once we know the slope, we follow the upper nodes,
     1618                 * inserting the same slope value into sg, until we reach the surface: */
     1619               
     1620                for(i=0;i<3;i++){
     1621       
     1622                        node=nodes[i]; //base nodes
     1623       
     1624                        /*get velocity for this base node: */
     1625                        slope=sg_serial[doflist[i]];
     1626
     1627                        //go throsgn all nodes which sit on top of this node, until we reach the surface,
     1628                        //and plsg  slope in sg
     1629                        for(;;){
     1630
     1631                                node->GetDofList(&nodedof,&numberofdofspernode);
     1632                                VecSetValues(sg,1,&nodedof,&slope,INSERT_VALUES);
     1633
     1634                                if (node->IsOnSurface())break;
     1635                                /*get next node: */
     1636                                node=node->GetUpperNode();
     1637                        }
     1638                }
     1639
     1640        }
     1641
     1642}
     1643
     1644
     1645#undef __FUNCT__
    15771646#define __FUNCT__ "Penta:GetB_vert"
    15781647void Penta::GetB_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4){
    … …  
    17701839        double pressure[numgrids];
    17711840        double rho_ice,g;
    1772         double xyz_list[numgrids][3];
     1841        double       xyz_list[numgrids][3];
    17731842               
    17741843        /*Get node data: */
    … …  
    17921861
    17931862}
     1863
     1864
     1865#undef __FUNCT__
     1866#define __FUNCT__ "Penta::CreateKMatrixSlopeCompute"
     1867
     1868void  Penta::CreateKMatrixSlopeCompute(Mat Kgg,void* inputs,int analysis_type){
     1869
     1870        /*Collapsed formulation: */
     1871        Tria*  tria=NULL;
     1872       
     1873        /*Is this element on the bed? :*/
     1874        if(!onbed)return;
     1875
     1876        /*Spawn Tria element from the base of the Penta: */
     1877        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     1878        tria->CreateKMatrix(Kgg,inputs, analysis_type);
     1879        delete tria;
     1880        return;
     1881
     1882}
     1883
     1884#undef __FUNCT__
     1885#define __FUNCT__ "Penta::CreatePVectorSlopeCompute"
     1886
     1887void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type){
     1888       
     1889        /*Collapsed formulation: */
     1890        Tria*  tria=NULL;
     1891       
     1892        /*Is this element on the bed? :*/
     1893        if(!onbed)return;
     1894
     1895        /*Spawn Tria element from the base of the Penta: */
     1896        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     1897        tria->CreatePVector(pg,inputs, analysis_type);
     1898        delete tria;
     1899        return;
     1900}
     1901

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

    r301 r358  
    110110                void  GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4);
    111111                void  VelocityExtrude(Vec ug,double* ug_serial);
     112                void  SlopeExtrude(Vec sg,double* sg_serial);
    112113                void  ComputePressure(Vec p_g);
     114                void  CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type);
     115                void  CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type);
     116
    113117
    114118};

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

    r352 r358  
    274274
    275275        }
     276        else if (
     277                        (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) ||
     278                        (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) ||
     279                        (analysis_type==BedSlopeComputeXAnalysisEnum()) ||
     280                        (analysis_type==BedSlopeComputeYAnalysisEnum())
     281                        ){
     282               
     283                CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type);
     284        }
    276285        else{
    277286                throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
    … …  
    666675
    667676#undef __FUNCT__
     677#define __FUNCT__ "Tria::CreateKMatrixSlopeCompute"
     678
     679void  Tria::CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type){
     680
     681        /* local declarations */
     682        int             i,j;
     683
     684        /* node data: */
     685        const int    numgrids=3;
     686        const int    NDOF1=1;
     687        const int    numdofs=NDOF1*numgrids;
     688        double       xyz_list[numgrids][3];
     689        int          doflist[numdofs];
     690        int          numberofdofspernode;
     691       
     692        /* gaussian points: */
     693        int     num_gauss,ig;
     694        double* first_gauss_area_coord  =  NULL;
     695        double* second_gauss_area_coord =  NULL;
     696        double* third_gauss_area_coord  =  NULL;
     697        double* gauss_weights           =  NULL;
     698        double  gauss_weight;
     699        double  gauss_l1l2l3[3];
     700
     701        /* matrices: */
     702        double L[1][3];
     703        double DL_scalar;
     704
     705        /* local element matrices: */
     706        double Ke_gg[numdofs][numdofs]; //local element stiffness matrix
     707        double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix evaluated at the gaussian point.
     708       
     709        double Jdet;
     710       
     711        /* Get node coordinates and dof list: */
     712        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     713        GetDofList(&doflist[0],&numberofdofspernode);
     714
     715        /* Set Ke_gg to 0: */
     716        for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
     717
     718        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     719        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     720
     721        /* Start  looping on the number of gaussian points: */
     722        for (ig=0; ig<num_gauss; ig++){
     723                /*Pick up the gaussian point: */
     724                gauss_weight=*(gauss_weights+ig);
     725                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     726                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     727                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     728
     729               
     730                /*Get L matrix: */
     731                GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
     732
     733                /* Get Jacobian determinant: */
     734                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     735               
     736                DL_scalar=gauss_weight*Jdet;
     737
     738                /*  Do the triple producte tL*D*L: */
     739                TripleMultiply( &L[0][0],1,3,1,
     740                                        &DL_scalar,1,1,0,
     741                                        &L[0][0],1,3,0,
     742                                        &Ke_gg_gaussian[0][0],0);
     743
     744                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     745                for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     746        } //for (ig=0; ig<num_gauss; ig++
     747
     748        /*Add Ke_gg to global matrix Kgg: */
     749        MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
     750               
     751        cleanup_and_return:
     752        xfree((void**)&first_gauss_area_coord);
     753        xfree((void**)&second_gauss_area_coord);
     754        xfree((void**)&third_gauss_area_coord);
     755        xfree((void**)&gauss_weights);
     756}
     757
     758#undef __FUNCT__
    668759#define __FUNCT__ "Tria::CreatePVector"
    669760void  Tria::CreatePVector(Vec pg,void* inputs,int analysis_type){
    … …  
    672763        if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
    673764                CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type);
     765        }
     766        else if (
     767                        (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) ||
     768                        (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) ||
     769                        (analysis_type==BedSlopeComputeXAnalysisEnum()) ||
     770                        (analysis_type==BedSlopeComputeYAnalysisEnum())
     771                        ){
     772               
     773                CreatePVectorSlopeCompute( pg,inputs,analysis_type);
    674774        }
    675775        else{
    … …  
    849949}
    850950
     951#undef __FUNCT__
     952#define __FUNCT__ "Tria::CreatePVectorSlopeCompute"
     953
     954void Tria::CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type){
     955
     956        int             i,j;
     957
     958        /* node data: */
     959        const int    numgrids=3;
     960        const int    NDOF1=1;
     961        const int    numdofs=NDOF1*numgrids;
     962        double       xyz_list[numgrids][3];
     963        int          doflist[numdofs];
     964        int          numberofdofspernode;
     965       
     966        /* gaussian points: */
     967        int     num_gauss,ig;
     968        double* first_gauss_area_coord  =  NULL;
     969        double* second_gauss_area_coord =  NULL;
     970        double* third_gauss_area_coord  =  NULL;
     971        double* gauss_weights           =  NULL;
     972        double  gauss_weight;
     973        double  gauss_l1l2l3[3];
     974
     975        /* Jacobian: */
     976        double Jdet;
     977
     978        /*nodal functions: */
     979        double l1l2l3[3];
     980
     981        /*element vector at the gaussian points: */
     982        double  pe_g[numdofs];
     983        double  pe_g_gaussian[numdofs];
     984        double  param[3];
     985        double  slope[2];
     986
     987        /* Get node coordinates and dof list: */
     988        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     989        GetDofList(&doflist[0],&numberofdofspernode);
     990
     991        /* Set pe_g to 0: */
     992        for(i=0;i<numdofs;i++) pe_g[i]=0.0;
     993
     994        if ( (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || (analysis_type==SurfaceSlopeComputeYAnalysisEnum())){
     995                for(i=0;i<numdofs;i++) param[i]=s[i];
     996        }
     997        if ( (analysis_type==BedSlopeComputeXAnalysisEnum()) || (analysis_type==BedSlopeComputeYAnalysisEnum())){
     998                for(i=0;i<numdofs;i++) param[i]=b[i];
     999        }
     1000
     1001        /* Get gaussian points and weights: */
     1002        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/
     1003
     1004
     1005        /* Start  looping on the number of gaussian points: */
     1006        for (ig=0; ig<num_gauss; ig++){
     1007                /*Pick up the gaussian point: */
     1008                gauss_weight=*(gauss_weights+ig);
     1009                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     1010                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     1011                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     1012
     1013                GetParameterDerivativeValue(&slope[0], &param[0],&xyz_list[0][0], gauss_l1l2l3);
     1014               
     1015                /* Get Jacobian determinant: */
     1016                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     1017               
     1018                 /*Get nodal functions: */
     1019                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     1020
     1021                /*Build pe_g_gaussian vector: */
     1022                if ( (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || (analysis_type==BedSlopeComputeXAnalysisEnum())){
     1023                        for(i=0;i<numdofs;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[0]*l1l2l3[i];
     1024                }
     1025                if ( (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) || (analysis_type==BedSlopeComputeYAnalysisEnum())){
     1026                        for(i=0;i<numdofs;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[1]*l1l2l3[i];
     1027                }
     1028
     1029                /*Add pe_g_gaussian vector to pe_g: */
     1030                for( i=0; i<numdofs; i++)pe_g[i]+=pe_g_gaussian[i];
     1031
     1032        } //for (ig=0; ig<num_gauss; ig++)
     1033
     1034        /*Add pe_g to global vector pg: */
     1035        VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
     1036
     1037        cleanup_and_return:
     1038        xfree((void**)&first_gauss_area_coord);
     1039        xfree((void**)&second_gauss_area_coord);
     1040        xfree((void**)&third_gauss_area_coord);
     1041        xfree((void**)&gauss_weights);
     1042
     1043}
    8511044
    8521045#undef __FUNCT__
    … …  
    20542247
    20552248#undef __FUNCT__
     2249#define __FUNCT__ "Tria::CreateKMatrixDiagnosticBaseVert"
     2250void  Tria::CreateKMatrixDiagnosticBaseVert(Mat Kgg,void* vinputs,int analysis_type){
     2251
     2252        int i,j;
     2253
     2254        /* node data: */
     2255        const int    numgrids=3;
     2256        const int    NDOF1=1;
     2257        const int    numdofs=NDOF1*numgrids;
     2258        double       xyz_list[numgrids][3];
     2259        int          doflist[numdofs];
     2260        int          numberofdofspernode;
     2261       
     2262        /* gaussian points: */
     2263        int     num_gauss,ig;
     2264        double* first_gauss_area_coord  =  NULL;
     2265        double* second_gauss_area_coord =  NULL;
     2266        double* third_gauss_area_coord  =  NULL;
     2267        double* gauss_weights           =  NULL;
     2268        double  gauss_weight;
     2269        double  gauss_l1l2l3[3];
     2270
     2271
     2272        /* matrices: */
     2273        double L[1][3];
     2274        double DL_scalar;
     2275
     2276        double Ke_gg[3][3]; //stiffness matrix
     2277        double Ke_gg_gaussian[3][3]; //stiffness matrix evaluated at the gaussian point.
     2278        double Jdet;
     2279
     2280        ParameterInputs* inputs=NULL;
     2281
     2282        /*recover pointers: */
     2283        inputs=(ParameterInputs*)vinputs;
     2284
     2285        /* Get node coordinates and dof list: */
     2286        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     2287        GetDofList(&doflist[0],&numberofdofspernode);
     2288
     2289        /* Set Ke_gg to 0: */
     2290        for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
     2291
     2292        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     2293        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2294
     2295        /* Start  looping on the number of gaussian points: */
     2296        for (ig=0; ig<num_gauss; ig++){
     2297                /*Pick up the gaussian point: */
     2298                gauss_weight=*(gauss_weights+ig);
     2299                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     2300                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     2301                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     2302
     2303                /* Get Jacobian determinant: */
     2304                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     2305       
     2306                //Get L matrix if viscous basal drag present:
     2307                GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
     2308
     2309                DL_scalar=gauss_weight*Jdet;
     2310               
     2311                /*  Do the triple producte tL*D*L: */
     2312                TripleMultiply( &L[0][0],1,3,1,
     2313                                        &DL_scalar,1,1,0,
     2314                                        &L[0][0],1,3,0,
     2315                                        &Ke_gg_gaussian[0][0],0);
     2316
     2317               
     2318                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     2319                for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     2320        } //for (ig=0; ig<num_gauss; ig++
     2321
     2322        /*Add Ke_gg to global matrix Kgg: */
     2323        MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
     2324               
     2325        cleanup_and_return:
     2326        xfree((void**)&first_gauss_area_coord);
     2327        xfree((void**)&second_gauss_area_coord);
     2328        xfree((void**)&third_gauss_area_coord);
     2329        xfree((void**)&gauss_weights);
     2330}
     2331
     2332#undef __FUNCT__
    20562333#define __FUNCT__ "Tria::CreateKMatrixDiagnosticSurfaceVert"
    20572334void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* vinputs,int analysis_type){

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

    r352 r358  
    7272                void  CreateKMatrixDiagnosticHoriz(Mat Kgg,void* inputs,int analysis_type);
    7373                void  CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* inputs,int analysis_type);
     74                void  CreateKMatrixDiagnosticBaseVert(Mat Kgg,void* inputs,int analysis_type);
    7475                void  CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* inputs,int analysis_type);
     76                void  CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type);
    7577                void  GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
    7678                void  GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3);
    … …  
    9395                void  CreatePVectorDiagnosticHoriz(Vec pg,void* inputs,int analysis_type);
    9496                void  CreatePVectorDiagnosticBaseVert(Vec pg,void* inputs,int analysis_type);
     97                void  CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type);
    9598                void* GetMatPar();
    9699                int   GetShelf();

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

    r246 r358  
    5555        fid=fopen(inputfilename,"rb");
    5656        if(fid==NULL) throw ErrorException(__FUNCT__,exprintf("%s%s%s","could not open file ",inputfilename," for binary reading"));
    57        
     57
    5858        _printf_("read and create finite element model:\n");
    5959        CreateFemModel(&femmodel,fid,analysis_type);
    … …  
    6262        femmodel.parameters->FindParam((void*)&u_g_initial,"u_g");
    6363        femmodel.parameters->FindParam((void*)&numberofnodes,"numberofnodes");
    64        
     64
    6565        inputs=new ParameterInputs;
    6666        inputs->Add("velocity",u_g_initial,3,numberofnodes);

  • issm/trunk/src/c/shared/Elements/ResolvePointers.cpp

    r1 r358  
    2929        for(i=0;i<num_objects;i++){
    3030
     31                /*is this object id -1? If so, drop this search, it was not requested: */
     32                if (object_ids[i]==-1)continue;
    3133
    3234                /*Check whether existing objects are correct: */

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

    r304 r358  
    2828                numdofs=2;
    2929        }
    30         else if (analysis_type==SurfaceSlopeComputeAnalysisEnum()){
    31                 numdofs=1;
    32         }
    33         else if (analysis_type==BedSlopeComputeAnalysisEnum()){
     30        else if (analysis_type==SlopeComputeAnalysisEnum()){
    3431                numdofs=1;
    3532        }

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

    r132 r358  
    5555        end
    5656end
     57
     58if (md.ismacayealpattyn==0 && md.ishutter==0 && md.isstokes==0),
     59        disp(['no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']);
     60        bool=0;return;
     61end
     62
    5763
    5864%NO NAN

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

    r308 r358  
    2323end
    2424
    25 WriteData(fid,solutiontype,'String','analysis_type');
     25if strcmp(solutiontype,'diagnostic'),
     26        WriteData(fid,'diagnostic_horiz','String','analysis_type');
     27end
    2628WriteData(fid,md.type,'String','type');
    2729WriteData(fid,md.numberofgrids,'Integer','numberofgrids');

  • issm/trunk/src/m/solutions/cielo/diagnostic.m

    r302 r358  
    2222       
    2323        displaystring(md.debug,'\n%s',['reading surface and bed slope computation model data']);
    24         md.analysis_type='surface_slope_compute'; m_ss=CreateFemModel(md);
    25         md.analysis_type='bed_slope_compute'; m_bs=CreateFemModel(md);
     24        md.analysis_type='slope_compute'; m_sl=CreateFemModel(md);
    2625       
    2726        % figure out number of dof: just for information purposes.
    28         md.dof=m_dh.nodesets.fsize; %biggest dof number
     27        md.dof=modelsize(m_dh,m_dv,m_ds,m_dhu,m_sl);
    2928
    3029        %initialize inputs
    … …  
    3332
    3433        %compute solution
    35         [u_g,p_g]=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_ss,m_bs,inputs);
     34        [u_g,p_g]=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_sl,inputs);
    3635
    3736        %Load results onto model

  • issm/trunk/src/m/solutions/cielo/diagnostic_core.m

    r326 r358  
    1 function [u_g p_g]=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_ss,m_bs,inputs);
     1function [u_g p_g]=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_sl,inputs);
    22%DIAGNOSTIC_CORE - compute the core velocity field
    33%
    44%   Usage:
    5 %      u_g=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_ss,m_bs,inputs);
     5%      u_g=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_sl,inputs);
    66%
    77
    … …  
    1616
    1717        displaystring(debug,'\n%s',['computing surface slope (x and y derivatives)...']);
    18         slopex=diagnostic_core_linear(m_ss,inputs,'surface_slope_compute_x');
    19         slopey=diagnostic_core_linear(m_ss,inputs,'surface_slope_compute_y');
     18        slopex=diagnostic_core_linear(m_sl,inputs,'surface_slope_compute_x');
     19        slopey=diagnostic_core_linear(m_sl,inputs,'surface_slope_compute_y');
    2020
    2121        if dim==3,
    2222                displaystring(debug,'\n%s',['extruding slopes in 3d...']);
    23                 slopex=SlopeExtrude(m_ss.elements,m_ss.nodes,m_ss.loads,m_ss.materials,slopex);
    24                 slopey=SlopeExtrude(m_ss.elements,m_ss.nodes,m_ss.loads,m_ss.materials,slopey);
     23                slopex=SlopeExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopex);
     24                slopey=SlopeExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopey);
    2525        end
    2626
    2727        displaystring(debug,'\n%s',['computing slopes...']);
    28         inputs=add(inputs,'surfaceslopex',slopex,'doublevec',m_ss.parameters.numberofdofspernode,m_ss.parameters.numberofnodes);
    29         inputs=add(inputs,'surfaceslopey',slopey,'doublevec',m_ss.parameters.numberofdofspernode,m_ss.parameters.numberofnodes);
     28        inputs=add(inputs,'surfaceslopex',slopex,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes);
     29        inputs=add(inputs,'surfaceslopey',slopey,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes);
    3030
    3131        displaystring(debug,'\n%s',['computing hutter velocities...']);
    … …  
    7575
    7676                displaystring(debug,'\n%s',['computing bed slope (x and y derivatives)...']);
    77                 slopex=diagnostic_core_linear(m_bs,inputs,'bed_slope_compute_x');
    78                 slopey=diagnostic_core_linear(m_bs,inputs,'bed_slope_compute_y');
    79                 slopex=SlopeExtrude(m_bs.elements,m_bs.nodes,m_bs.loads,m_bs.materials,slopex);
    80                 slopey=SlopeExtrude(m_bs.elements,m_bs.nodes,m_bs.loads,m_bs.materials,slopey);
     77                slopex=diagnostic_core_linear(m_sl,inputs,'bed_slope_compute_x');
     78                slopey=diagnostic_core_linear(m_sl,inputs,'bed_slope_compute_y');
     79                slopex=SlopeExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopex);
     80                slopey=SlopeExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopey);
    8181
    82                 inputs=add(inputs,'bedslopex',slopex,'doublevec',m_ss.parameters.numberofdofspernode,m_ss.parameters.numberofnodes);
    83                 inputs=add(inputs,'bedslopey',slopey,'doublevec',m_ss.parameters.numberofdofspernode,m_ss.parameters.numberofnodes);
     82                inputs=add(inputs,'bedslopex',slopex,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes);
     83                inputs=add(inputs,'bedslopey',slopey,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes);
    8484               
    8585                inputs=add(inputs,'velocity',u_g,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes);

  • issm/trunk/src/m/solutions/cielo/diagnostic_core_linear.m

    r326 r358  
    1010        %Update inputs in datasets
    1111        [m.elements,m.nodes, m.loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, m.loads,m.materials,inputs);
    12 
     12       
    1313        %system matrices
    1414        [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs);
    15 
     15       
    1616        %Reduce tangent matrix from g size to f size
    1717        [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets);
    … …  
    1919        %Reduce load from g size to f size
    2020        [p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets);
    21 
     21       
    2222        %Solve 
    2323        [u_f]=Solver(K_ff,p_f,[],m.parameters);
    24 
     24       
    2525        %Merge back to g set
    2626        [u_g]= Mergesolutionfromftog( u_f, m.Gmn, m.ys, m.nodesets );
    27        
     27
    2828end %end function

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

    r304 r358  
    3333                                Reducevectorgtos\
    3434                                SetStructureField\
     35                                SlopeExtrude\
    3536                                Solver\
    3637                                SpcNodes\
    … …  
    154155                          SetStructureField/SetStructureField.h
    155156
     157SlopeExtrude_SOURCES = SlopeExtrude/SlopeExtrude.cpp\
     158                          SlopeExtrude/SlopeExtrude.h
     159
    156160Solver_SOURCES = Solver/Solver.cpp\
    157161                          Solver/Solver.h
Note: See TracChangeset for help on using the changeset viewer.