Changeset 3359


Ignore:
Timestamp:
03/31/10 16:43:20 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added discontinuous galerkin elements and prognostic2 solution, to be continued...

Location:
issm/trunk/src
Files:
3 added
26 edited

Legend:

Unmodified
Added
Removed

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

    r3332 r3359  
    2121
    2222
     23        _printf_("   Configuring elements...\n");
    2324        elements->Configure(elements,loads,nodes,materials,parameters);
     25        _printf_("   Configuring loads...\n");
    2426        loads->Configure(elements,loads,nodes,materials,parameters);
     27        _printf_("   Configuring nodes...\n");
    2528        nodes->Configure(elements,loads,nodes,materials,parameters);
     29        _printf_("   Configuring parameters...\n");
    2630        parameters->Configure(elements,loads, nodes, materials,parameters);
    2731

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

    r3355 r3359  
    266266                        dataset->AddObject(icefront);
    267267                }
     268                else if(enum_type==NumericalfluxEnum()){
     269                        Numericalflux* numericalflux=NULL;
     270                        numericalflux=new Numericalflux();
     271                        numericalflux->Demarshall(&marshalled_dataset);
     272                        dataset->AddObject(numericalflux);
     273                }
    268274                else if(enum_type==RgbEnum()){
    269275                        Rgb* rgb=NULL;
    … …  
    816822                }
    817823                if(EnumIsLoad((*object)->Enum())){
    818 
    819824                        load=(Load*)(*object);
    820 
    821825                        load->Configure(elements,nodes,materials);
    822826                }

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

    r3354 r3359  
    8080int PenpairEnum(void){                      return          433; }
    8181int PengridEnum(void){                      return          434; }
     82int NumericalfluxEnum(void){                return          435; }
    8283/*Materials: */
    8384int MaterialEnum(void){                     return          440; }

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

    r3354 r3359  
    8181int PenpairEnum(void);
    8282int PengridEnum(void);
     83int NumericalfluxEnum(void);
    8384/*Materials: */
    8485int MaterialEnum(void);

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

    r3354 r3359  
    8484                                        ./objects/Riftfront.cpp\
    8585                                        ./objects/Riftfront.h\
     86                                        ./objects/Numericalflux.cpp\
     87                                        ./objects/Numericalflux.h\
    8688                                        ./objects/Param.cpp\
    8789                                        ./objects/Param.h\
    … …  
    468470                                        ./objects/Riftfront.cpp\
    469471                                        ./objects/Riftfront.h\
     472                                        ./objects/Numericalflux.cpp\
     473                                        ./objects/Numericalflux.h\
    470474                                        ./objects/Param.cpp\
    471475                                        ./objects/Param.h\

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

    r3332 r3359  
    241241        count++;
    242242        param= new Param(count,"numberofnodes",DOUBLE);
    243         param->SetDouble(iomodel->numberofnodes);
     243        if (iomodel->analysis_type==Prognostic2AnalysisEnum()) param->SetDouble(3*iomodel->numberofelements);
     244        else param->SetDouble(iomodel->numberofnodes);
    244245        parameters->AddObject(param);
    245246

  • issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp

    r3332 r3359  
    9595                }
    9696
    97 
    9897                element=(int)(*(iomodel->pressureload+segment_width*i+segment_width-2)-1); //element is in the penultimate column (grid1 grid2 ... elem fill)
    9998

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

    r3332 r3359  
    9595        iomodel-> spctemperature=NULL;
    9696        iomodel-> spcthickness=NULL;
     97        iomodel->numberofedges=0;
     98        iomodel->edges=NULL;
    9799       
    98100        /*!materials: */
    … …  
    247249        xfree((void**)&iomodel->spcthickness);
    248250        xfree((void**)&iomodel->spctemperature);
     251        xfree((void**)&iomodel->edges);
    249252        xfree((void**)&iomodel->geothermalflux);
    250253        xfree((void**)&iomodel->melting);

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

    r3354 r3359  
    9595        double* spcthickness;
    9696        double* geothermalflux;
     97        int     numberofedges;
     98        double* edges;
    9799       
    98100        /*materials: */

  • issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateConstraintsPrognostic2.cpp

    r3354 r3359  
    22 * CreateConstraintsPrognostic2.c:
    33 */
    4 
    54
    65#include "../../DataSet/DataSet.h"
    … …  
    1110#include "../IoModel.h"
    1211
    13 
    1412void    CreateConstraintsPrognostic2(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
    15 
    16 
    17         int i;
    18         int count;
    1913       
    2014        DataSet* constraints = NULL;
    21         Spc*    spc  = NULL;
    22 
    23         /*spc intermediary data: */
    24         int spc_sid;
    25         int spc_node;
    26         int spc_dof;
    27         double spc_value;
    28        
    29         double* spcthickness=NULL;
    3015       
    3116        /*Create constraints: */
    3217        constraints = new DataSet(ConstraintsEnum());
    33 
    34         /*Fetch data: */
    35         IoModelFetchData(&spcthickness,NULL,NULL,iomodel_handle,"spcthickness");
    36 
    37         count=0;
    38 
    39         /*Create spcs from x,y,z, as well as the spc values on those spcs: */
    40         for (i=0;i<iomodel->numberofnodes;i++){
    41         #ifdef _PARALLEL_
    42         /*keep only this partition's nodes:*/
    43         if((iomodel->my_grids[i]==1)){
    44         #endif
    45 
    46                 if ((int)spcthickness[2*i]){
    47        
    48                         /*This grid needs to be spc'd: */
    49 
    50                         spc_sid=count;
    51                         spc_node=i+1;
    52                         spc_dof=1; //we enforce first translation degree of freedom, for temperature
    53                         spc_value=*(spcthickness+2*i+1);
    54 
    55                         spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
    56                         constraints->AddObject(spc);
    57                         count++;
    58                 }
    59 
    60         #ifdef _PARALLEL_
    61         } //if((my_grids[i]==1))
    62         #endif
    63         }
    64 
    65         /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
    66          * datasets, it will not be redone: */
    67         constraints->Presort();
    68 
    69         /*Free data: */
    70         xfree((void**)&spcthickness);
    71        
    72         cleanup_and_return:
    7318       
    7419        /*Assign output pointer: */

  • issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp

    r3354 r3359  
    246246        else{ //        if (strcmp(meshtype,"2d")==0)
    247247                ISSMERROR(exprintf("not implemented yet"));
    248 
    249                 /*Fetch data needed: */
    250                 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
    251                 IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
    252                 IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface");
    253                 IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
    254                 IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf");
    255                 IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed");
    256                 IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface");
    257                 IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
    258        
    259                 for (i=0;i<iomodel->numberofelements;i++){
    260                 #ifdef _PARALLEL_
    261                 /*We are using our element partition to decide which elements will be created on this node: */
    262                 if(my_rank==epart[i]){
    263                 #endif
    264 
    265                        
    266                         /*name and id: */
    267                         penta_id=i+1; //matlab indexing.
    268                         penta_mid=-1;
    269                         penta_mparid=-1; //no need for materials
    270                         penta_numparid=1;
    271 
    272                         /*vertices,thickness,surface,bed and drag: */
    273                         for(j=0;j<6;j++){
    274                                 penta_g[j]=(int)*(iomodel->elements+elements_width*i+j);
    275                                 penta_h[j]=*(iomodel->thickness+    ((int)*(iomodel->elements+elements_width*i+j)-1));
    276                                 penta_s[j]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+j)-1));
    277                                 penta_b[j]=*(iomodel->bed+    ((int)*(iomodel->elements+elements_width*i+j)-1));
    278                         }
    279 
    280                         /*diverse: */
    281                         penta_shelf=(int)*(iomodel->elementoniceshelf+i);
    282                         penta_onbed=(int)*(iomodel->elementonbed+i);
    283                         penta_onsurface=(int)*(iomodel->elementonsurface+i);
    284                         penta_collapse=1;
    285                         penta_onwater=(bool)*(iomodel->elementonwater+i);
    286        
    287 
    288                         /*Create Penta using its constructor:*/
    289                         penta= new Penta( penta_id,penta_mid,penta_mparid,penta_numparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
    290                                         penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,
    291                                         penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,
    292                                         penta_thermal_steadystate,penta_onwater);
    293 
    294                         /*Add penta element to elements dataset: */
    295                         elements->AddObject(penta);
    296        
    297                         #ifdef _PARALLEL_
    298                         /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
    299                          *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
    300                          into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
    301                          will hold which grids belong to this partition*/
    302                         my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
    303                         my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
    304                         my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
    305                         my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;
    306                         my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;
    307                         my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;
    308                         #endif
    309 
    310                 #ifdef _PARALLEL_
    311                 }//if(my_rank==epart[i])
    312                 #endif
    313 
    314                 }//for (i=0;i<numberofelements;i++)
    315 
    316                 /*Free data: */
    317                 xfree((void**)&iomodel->elements);
    318                 xfree((void**)&iomodel->thickness);
    319                 xfree((void**)&iomodel->surface);
    320                 xfree((void**)&iomodel->bed);
    321                 xfree((void**)&iomodel->elementoniceshelf);
    322                 xfree((void**)&iomodel->elementonbed);
    323                 xfree((void**)&iomodel->elementonsurface);
    324                 xfree((void**)&iomodel->elementonwater);
    325 
    326248        } //if (strcmp(meshtype,"2d")==0)
    327249
    … …  
    373295        materials->AddObject(matpar);
    374296
    375         /*Partition penalties in 3d: */
    376         if(strcmp(iomodel->meshtype,"3d")==0){
    377        
    378                 /*Get penalties: */
    379                 IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");
    380 
    381                 if(iomodel->numpenalties){
    382 
    383                         iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int));
    384                         #ifdef _SERIAL_
    385                         for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1;
    386                         #else
    387                         for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1;
    388 
    389                         for(i=0;i<iomodel->numpenalties;i++){
    390                                 first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1);
    391                                 if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition  grids
    392                                         /*All grids that are being penalised belong to this node's internal grid partition.:*/
    393                                         iomodel->penaltypartitioning[i]=1;
    394                                 }
    395                                 if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
    396                                         iomodel->penaltypartitioning[i]=0;
    397                                 }
    398                         }
    399                         #endif
    400                 }
    401 
    402                 /*Free penalties: */
    403                 xfree((void**)&iomodel->penalties);
    404         }
    405 
    406297        /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids.
    407298         We can therefore determine  which grids are internal to this node's partition
    … …  
    416307                IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
    417308        }
     309        IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
    418310        IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x");
    419311        IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y");
    … …  
    430322        DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);
    431323
    432         for (i=0;i<iomodel->numberofnodes;i++){
    433         #ifdef _PARALLEL_
    434         /*keep only this partition's nodes:*/
    435         if((my_grids[i]==1)){
    436         #endif
    437 
    438                 node_id=i+1; //matlab indexing
    439                        
    440                
    441                
    442                 #ifdef _PARALLEL_
    443                 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
    444                         node_partitionborder=1;
    445                 }
    446                 else{
     324        /*Build Nodes dataset -> 3 for each element*/
     325        for (i=0;i<iomodel->numberofelements;i++){
     326                for (j=0;j<3;j++){
     327
     328                        #ifdef _PARALLEL_
     329                        /*!All elements have been partitioned above, only create elements for this CPU: */
     330                        if(my_rank==epart[i]){
     331                        #endif
     332
     333                        //Get id of the vertex on which the current node is located
     334                        k=(int)*(iomodel->elements+elements_width*i+j)-1; //(Matlab to C indexing)
     335
     336                        //Get node properties
     337                        node_id=i*3+j+1;
    447338                        node_partitionborder=0;
    448                 }
    449                 #else
    450                         node_partitionborder=0;
    451                 #endif
    452 
    453                 node_x[0]=iomodel->x[i];
    454                 node_x[1]=iomodel->y[i];
    455                 node_x[2]=iomodel->z[i];
    456                 node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
    457                
    458                 node_onbed=(int)iomodel->gridonbed[i];
    459                 node_onsurface=(int)iomodel->gridonsurface[i]; 
    460                 node_onshelf=(int)iomodel->gridoniceshelf[i];   
    461                 node_onsheet=(int)iomodel->gridonicesheet[i];   
    462 
    463                 if (strcmp(iomodel->meshtype,"3d")==0){
    464                         if (isnan(iomodel->uppernodes[i])){
    465                                 node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
     339                        node_x[0]=iomodel->x[k];
     340                        node_x[1]=iomodel->y[k];
     341                        node_x[2]=iomodel->z[k];
     342                        node_sigma=(iomodel->z[k]-iomodel->bed[k])/(iomodel->thickness[k]);
     343                        node_onbed=(int)iomodel->gridonbed[k];
     344                        node_onsurface=(int)iomodel->gridonsurface[k]; 
     345                        node_onshelf=(int)iomodel->gridoniceshelf[k];   
     346                        node_onsheet=(int)iomodel->gridonicesheet[k];   
     347                        if (strcmp(iomodel->meshtype,"3d")==0){
     348                                if (isnan(iomodel->uppernodes[k])){
     349                                        node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
     350                                }
     351                                else{
     352                                        node_upper_node_id=(int)iomodel->uppernodes[k];
     353                                }
    466354                        }
    467355                        else{
    468                                 node_upper_node_id=(int)iomodel->uppernodes[i];
     356                                /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
     357                                node_upper_node_id=node_id;
    469358                        }
     359
     360                        /*Create node using its constructor: */
     361                        node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
     362
     363                        /*Add node to nodes dataset: */
     364                        nodes->AddObject(node);
     365
     366                        #ifdef _PARALLEL_
     367                        }
     368                        #endif
    470369                }
    471                 else{
    472                         /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
    473                         node_upper_node_id=node_id;
    474                 }
    475 
    476                 /*Create node using its constructor: */
    477                 node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
    478 
    479                 /*set single point constraints.: */
    480                 if (strcmp(iomodel->meshtype,"3d")==0){
    481                         /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
    482                         if (!iomodel->gridonbed[i]){
    483                                 for(k=1;k<=node_numdofs;k++){
    484                                         node->FreezeDof(k);
    485                                 }
    486                         }
    487                 }
    488                 /*Add node to nodes dataset: */
    489                 nodes->AddObject(node);
    490 
    491         #ifdef _PARALLEL_
    492         } //if((my_grids[i]==1))
    493         #endif
    494370        }
    495371
    … …  
    502378        /*Clean fetched data: */
    503379        xfree((void**)&iomodel->deadgrids);
     380        xfree((void**)&iomodel->elements);
    504381        xfree((void**)&iomodel->x);
    505382        xfree((void**)&iomodel->y);
    … …  
    513390        xfree((void**)&iomodel->gridoniceshelf);
    514391       
    515 
    516392        /*Keep partitioning information into iomodel*/
    517393        iomodel->epart=epart;
    … …  
    532408        *pnodes=nodes;
    533409        *pmaterials=materials;
    534 
    535410}

  • issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp

    r3354 r3359  
    11/*! \file CreateLoadsPrognostic2.c:
    22 */
    3 
    43
    54#include "../../DataSet/DataSet.h"
    … …  
    98#include "../../shared/shared.h"
    109#include "../../include/macros.h"
     10#include "../../include/typedefs.h"
    1111#include "../IoModel.h"
    12 
    1312
    1413void    CreateLoadsPrognostic2(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1514
    16         DataSet*    loads    = NULL;
     15        int i,j;
     16        int i1,i2;
     17        int pos1,pos2;
     18        double e1,e2;
     19
     20        extern int my_rank;
     21        extern int num_procs;
     22
     23        DataSet*       loads        = NULL;
     24        Numericalflux* numericalflux= NULL;
     25
     26        /*numericalflux intermediary data: */
     27        char numericalflux_type[NUMERICALFLUXSTRING];
     28        int  numericalflux_id;
     29        int  numericalflux_node_ids[MAX_NUMERICALFLUX_NODES];
     30        int  numericalflux_elem_ids[MAX_NUMERICALFLUX_ELEMS];
    1731
    1832        /*Create loads: */
    1933        loads   = new DataSet(LoadsEnum());
    20        
    21        
     34
     35        /*Get edges and elements*/
     36        IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges");
     37        IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
     38
     39        /*First load data:*/
     40        for (i=0;i<iomodel->numberofedges;i++){
     41
     42                /*Get left and right elements*/
     43                e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2]
     44                e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2]
     45
     46                #ifdef _PARALLEL_
     47                if (iomodel->epart[e1]!=my_rank){
     48                        /*This load does not belong to this cluster node, as it references an element
     49                         *that does not belong to this node's partition. Drop this 'i':*/
     50                        continue;
     51                }
     52                #endif
     53
     54                /*Create load*/
     55                numericalflux_id=i+1; //Matlab indexing
     56
     57                if (!isnan(e2)){
     58                        strcpy(numericalflux_type,"internal");
     59
     60                        numericalflux_elem_ids[0]=(int)e1+1;//id is in matlab index
     61                        numericalflux_elem_ids[1]=(int)e2+1;//id is in matlab index
     62
     63                        /*Now, we must get the nodes of the 4 nodes located on the edge*/
     64
     65                        /*1: Get vertices ids*/
     66                        i1=(int)iomodel->edges[4*i+0];
     67                        i2=(int)iomodel->edges[4*i+1];
     68
     69                        /*2: Get the column where these ids are located in the index*/
     70                        pos1=pos2=UNDEF;
     71                        for(j=0;j<3;j++){
     72                                if (iomodel->elements[3*(int)e1+j]==i1) pos1=j;
     73                                if (iomodel->elements[3*(int)e2+j]==i1) pos2=j;
     74                        }
     75                        ISSMASSERT(pos1!=UNDEF && pos2!=UNDEF);
     76
     77                        /*3: We have the id of the elements and the position of the vertices in the index
     78                         * we can compute their dofs!*/
     79                        numericalflux_node_ids[0]=3*(int)e1+pos1+1;      //id is in matlab index
     80                        numericalflux_node_ids[1]=3*(int)e1+(pos1+1)%3+1;
     81                        numericalflux_node_ids[2]=3*(int)e2+pos2+1;
     82                        numericalflux_node_ids[3]=3*(int)e2+(pos2-1)%3+1;
     83                }
     84                else{
     85                        strcpy(numericalflux_type,"boundary");
     86
     87                        numericalflux_elem_ids[0]=(int)e1+1;
     88
     89                        /*1: Get vertices ids*/
     90                        i1=(int)iomodel->edges[4*i+0];
     91                        i2=(int)iomodel->edges[4*i+1];
     92
     93                        /*2: Get the column where these ids are located in the index*/
     94                        pos1==UNDEF;
     95                        for(j=0;j<3;j++){
     96                                if (iomodel->elements[3*(int)e1+j]==i1) pos1=j;
     97                        }
     98                        ISSMASSERT(pos1!=UNDEF);
     99
     100                        /*3: We have the id of the elements and the position of the vertices in the index
     101                         * we can compute their dofs!*/
     102                        numericalflux_node_ids[0]=3*(int)e1+pos1+1; //id is in matlab index
     103                        numericalflux_node_ids[1]=3*(int)e1+(pos1+1)%3+1;
     104                }
     105
     106                numericalflux = new Numericalflux(numericalflux_id,numericalflux_type,numericalflux_node_ids,numericalflux_elem_ids);
     107
     108                loads->AddObject(numericalflux);
     109
     110        }
     111
     112        /*Free data: */
     113        xfree((void**)&iomodel->edges);
     114        xfree((void**)&iomodel->elements);
     115
    22116        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
    23117         * datasets, it will not be redone: */
    24118        loads->Presort();
    25 
    26         cleanup_and_return:
    27119
    28120        /*Assign output pointer: */
    … …  
    30122
    31123}
    32 
    33 

  • issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp

    r3354 r3359  
    99#include "../../objects/objects.h"
    1010#include "../../shared/shared.h"
     11#include "../../include/macros.h"
    1112#include "../IoModel.h"
    1213
    … …  
    1617        DataSet* parameters=NULL;
    1718        int      count;
    18         int      i;
     19        int      i,j;
    1920        int      dim;
     21        double*  elements;
     22        int*     part;
    2023
    2124        double* vx=NULL;
     25        double* vx_g=NULL;
    2226        double* vy=NULL;
    23         double* vz=NULL;
    24         double* u_g=NULL;
    25         double* pressure=NULL;
    26         double* temperature=NULL;
     27        double* vy_g=NULL;
    2728        double* thickness=NULL;
    28         double* surface=NULL;
    29         double* bed=NULL;
     29        double* h_g=NULL;
    3030        double* accumulation=NULL;
     31        double* a_g=NULL;
    3132        double* melting=NULL;
     33        double* m_g=NULL;
    3234
    3335        /*recover parameters : */
    … …  
    3638        count=parameters->Size();
    3739
    38         /*Get vx and vy: */
     40        /*Create transformatiob vector for DG inputs*/
     41        IoModelFetchData(&elements,NULL,NULL,iomodel_handle,"elements");
     42        part=(int*)xcalloc(iomodel->numberofelements*3,sizeof(int));
     43        for(i=0;i<iomodel->numberofelements;i++){
     44                for(j=0;j<3;j++){
     45                        part[3*i+j]=(int)elements[3*i+j]-1; //Matlab to C indexing
     46                        ISSMASSERT(part[3*i+j]<iomodel->numberofnodes);
     47                }
     48        }
     49        xfree((void**)&elements);
     50
     51        /*Get vx: */
    3952        IoModelFetchData(&vx,NULL,NULL,iomodel_handle,"vx");
    40         IoModelFetchData(&vy,NULL,NULL,iomodel_handle,"vy");
    41         IoModelFetchData(&vz,NULL,NULL,iomodel_handle,"vz");
    42 
    43         u_g=(double*)xcalloc(iomodel->numberofnodes*3,sizeof(double));
    44 
    45         if(vx) for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+0]=vx[i]/iomodel->yts;
    46         if(vy) for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+1]=vy[i]/iomodel->yts;
    47         if(vz) for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+2]=vz[i]/iomodel->yts;
     53        vx_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
     54        if(vx) for(i=0;i<3*iomodel->numberofelements;i++) vx_g[i]=vx[part[i]]/iomodel->yts;
    4855
    4956        count++;
    50         param= new Param(count,"u_g",DOUBLEVEC);
    51         param->SetDoubleVec(u_g,3*iomodel->numberofnodes,3);
     57        param= new Param(count,"vx_g",DOUBLEVEC);
     58        param->SetDoubleVec(vx_g,3*iomodel->numberofelements,1);
    5259        parameters->AddObject(param);
     60        xfree((void**)&vx);
    5361
     62        /*Get vy: */
     63        IoModelFetchData(&vy,NULL,NULL,iomodel_handle,"vy");
     64        vy_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
     65        if(vy) for(i=0;i<3*iomodel->numberofelements;i++) vy_g[i]=vy[part[i]]/iomodel->yts;
    5466
    55         xfree((void**)&vx);
     67        count++;
     68        param= new Param(count,"vy_g",DOUBLEVEC);
     69        param->SetDoubleVec(vy_g,3*iomodel->numberofelements,1);
     70        parameters->AddObject(param);
    5671        xfree((void**)&vy);
    57         xfree((void**)&vz);
    58         xfree((void**)&u_g);
    59 
    60         /*Get pressure if 3d iomodel: */
    61         parameters->FindParam(&dim,"dim");
    62         if (dim==3){
    63                 IoModelFetchData(&pressure,NULL,NULL,iomodel_handle,"pressure");
    64                
    65                 count++;
    66                 param= new Param(count,"p_g",DOUBLEVEC);
    67                 if(pressure) param->SetDoubleVec(pressure,iomodel->numberofnodes,1);
    68                 else param->SetDoubleVec(pressure,0,0);
    69                 parameters->AddObject(param);
    70 
    71                 /*Free pressure: */
    72                 xfree((void**)&pressure);
    73         }
    74 
    75         /*Get temperature if 3d iomodel: */
    76         parameters->FindParam(&dim,"dim");
    77         if (dim==3){
    78                 IoModelFetchData(&temperature,NULL,NULL,iomodel_handle,"temperature");
    79                
    80                 count++;
    81                 param= new Param(count,"t_g",DOUBLEVEC);
    82                 if(temperature) param->SetDoubleVec(temperature,iomodel->numberofnodes,1);
    83                 else param->SetDoubleVec(temperature,0,0);
    84                 parameters->AddObject(param);
    85 
    86                 /*Free temperature: */
    87                 xfree((void**)&temperature);
    88         }
    8972
    9073        /*Get thickness: */
    9174        IoModelFetchData(&thickness,NULL,NULL,iomodel_handle,"thickness");
    92        
     75        h_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
     76        if(thickness) for(i=0;i<3*iomodel->numberofelements;i++) h_g[i]=thickness[part[i]]/iomodel->yts;
     77
    9378        count++;
    9479        param= new Param(count,"h_g",DOUBLEVEC);
    95         if(thickness) param->SetDoubleVec(thickness,iomodel->numberofnodes,1);
    96         else param->SetDoubleVec(thickness,0,0);
     80        param->SetDoubleVec(h_g,3*iomodel->numberofelements,1);
    9781        parameters->AddObject(param);
    98 
    99         /*Free thickness: */
    10082        xfree((void**)&thickness);
    10183
    102         /*Get surface: */
    103         IoModelFetchData(&surface,NULL,NULL,iomodel_handle,"surface");
    104        
     84        /*Get accumulation: */
     85        IoModelFetchData(&accumulation,NULL,NULL,iomodel_handle,"accumulation");
     86        a_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
     87        if(accumulation) for(i=0;i<3*iomodel->numberofelements;i++) a_g[i]=accumulation[part[i]]/iomodel->yts;
     88
    10589        count++;
    106         param= new Param(count,"s_g",DOUBLEVEC);
    107         if(surface) param->SetDoubleVec(surface,iomodel->numberofnodes,1);
    108         else param->SetDoubleVec(surface,0,0);
     90        param= new Param(count,"a_g",DOUBLEVEC);
     91        param->SetDoubleVec(a_g,3*iomodel->numberofelements,1);
    10992        parameters->AddObject(param);
    110 
    111         /*Free surface: */
    112         xfree((void**)&surface);
    113 
    114         /*Get bed: */
    115         IoModelFetchData(&bed,NULL,NULL,iomodel_handle,"bed");
    116        
    117         count++;
    118         param= new Param(count,"b_g",DOUBLEVEC);
    119         if(bed) param->SetDoubleVec(bed,iomodel->numberofnodes,1);
    120         else param->SetDoubleVec(bed,0,0);
    121         parameters->AddObject(param);
    122 
    123         /*Free bed: */
    124         xfree((void**)&bed);
     93        xfree((void**)&accumulation);
    12594
    12695        /*Get melting: */
    12796        IoModelFetchData(&melting,NULL,NULL,iomodel_handle,"melting");
    128         if(melting) for(i=0;i<iomodel->numberofnodes;i++)melting[i]=melting[i]/iomodel->yts;
    129        
     97        m_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
     98        if(melting) for(i=0;i<3*iomodel->numberofelements;i++) m_g[i]=melting[part[i]]/iomodel->yts;
     99
    130100        count++;
    131101        param= new Param(count,"m_g",DOUBLEVEC);
    132         if(melting) param->SetDoubleVec(melting,iomodel->numberofnodes,1);
    133         else param->SetDoubleVec(melting,0,1);
     102        param->SetDoubleVec(m_g,3*iomodel->numberofelements,1);
    134103        parameters->AddObject(param);
    135 
    136         /*Free melting: */
    137104        xfree((void**)&melting);
    138105
    139         /*Get accumulation: */
    140         IoModelFetchData(&accumulation,NULL,NULL,iomodel_handle,"accumulation");
    141         if(accumulation) for(i=0;i<iomodel->numberofnodes;i++)accumulation[i]=accumulation[i]/iomodel->yts;
    142        
    143         count++;
    144         param= new Param(count,"a_g",DOUBLEVEC);
    145         if(accumulation) param->SetDoubleVec(accumulation,iomodel->numberofnodes,1);
    146         else param->SetDoubleVec(accumulation,0,0);
    147         parameters->AddObject(param);
    148 
    149         /*Free accumulation: */
    150         xfree((void**)&accumulation);
    151 
     106        /*Free partitioning vector*/
     107        xfree((void**)&part);
    152108
    153109        /*Assign output pointer: */

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

    r3332 r3359  
    1414void    CreateConstraintsSlopeCompute(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1515
    16 
    1716        DataSet* constraints = NULL;
    1817
    1918        /*Create constraints: */
    2019        constraints = new DataSet(ConstraintsEnum());
    21 
    22         /*Now, is the flag isstokes on? otherwise, do nothing: */
    23         if (!iomodel->isstokes)goto cleanup_and_return;
    24        
    25         cleanup_and_return:     
    2620       
    2721        /*Assign output pointer: */

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

    r3335 r3359  
    2020/*Assertion macro*/
    2121#define ISSMASSERT(statement)\
    22   if (!(statement)) ISSMERROR(exprintf("Assertion %s failed, please report bug to an ISSM developer",#statement))
     22  if (!(statement)) ISSMERROR(exprintf("Assertion \"%s\" failed, please report bug to an ISSM developer",#statement))
    2323
    2424/*The following macros hide the error exception handling in a matlab module. Just put

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

    r3332 r3359  
    173173        /*Link this load with its nodes: */
    174174        if (strcmp(type,"segment")==0){
    175                
    176175                ResolvePointers((Object**)nodes,node_ids,node_offsets,2,nodesin);
    177        
    178176        }
    179177        else{

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

    r3332 r3359  
    600600                }
    601601                else if (
     602                                (strcmp(field_name,"vx")==0) ||
     603                                (strcmp(field_name,"vy")==0) ||
    602604                                (strcmp(field_name,"thickness")==0) ||
    603605                                (strcmp(field_name,"temperature")==0) ||

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

    r3332 r3359  
    389389
    390390                CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type);
    391 
    392391        }
    393392        else if (analysis_type==PrognosticAnalysisEnum()){
    394393
    395394                CreateKMatrixPrognostic( Kgg,inputs,analysis_type,sub_analysis_type);
    396 
     395        }
     396        else if (analysis_type==Prognostic2AnalysisEnum()){
     397
     398                CreateKMatrixPrognostic2(Kgg,inputs,analysis_type,sub_analysis_type);
    397399        }
    398400        else if (analysis_type==BalancedthicknessAnalysisEnum()){
    399401
    400402                CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type);
    401 
    402403        }
    403404        else if (analysis_type==BalancedvelocitiesAnalysisEnum()){
    404405
    405406                CreateKMatrixBalancedvelocities( Kgg,inputs,analysis_type,sub_analysis_type);
    406 
    407407        }
    408408        else{
     409
    409410                ISSMERROR(exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
    410411        }
    … …  
    16041605                }
    16051606#endif
     1607        } // for (ig=0; ig<num_gauss; ig++)
     1608
     1609        /*Add Ke_gg to global matrix Kgg: */
     1610        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     1611
     1612#ifdef _DEBUGELEMENTS_
     1613        if(my_rank==RANK && id==ELID){
     1614                printf("      Ke_gg erms:\n");
     1615                for( i=0; i<numdof; i++){
     1616                        for (j=0;j<numdof;j++){
     1617                                printf("%g ",Ke_gg[i][j]);
     1618                        }
     1619                        printf("\n");
     1620                }
     1621                printf("      Ke_gg row_indices:\n");
     1622                for( i=0; i<numdof; i++){
     1623                        printf("%i ",doflist[i]);
     1624                }
     1625
     1626        }
     1627#endif
     1628
     1629cleanup_and_return:
     1630        xfree((void**)&first_gauss_area_coord);
     1631        xfree((void**)&second_gauss_area_coord);
     1632        xfree((void**)&third_gauss_area_coord);
     1633        xfree((void**)&gauss_weights);
     1634
     1635}
     1636/*}}}*/
     1637/*FUNCTION CreateKMatrixPrognostic2 {{{1*/
     1638void  Tria::CreateKMatrixPrognostic2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
     1639
     1640        /* local declarations */
     1641        int             i,j;
     1642
     1643        /* node data: */
     1644        const int    numgrids=3;
     1645        const int    NDOF1=1;
     1646        const int    numdof=NDOF1*numgrids;
     1647        double       xyz_list[numgrids][3];
     1648        int          doflist[numdof];
     1649        int          numberofdofspernode;
     1650
     1651        /* gaussian points: */
     1652        int     num_gauss,ig;
     1653        double* first_gauss_area_coord  =  NULL;
     1654        double* second_gauss_area_coord =  NULL;
     1655        double* third_gauss_area_coord  =  NULL;
     1656        double* gauss_weights           =  NULL;
     1657        double  gauss_weight;
     1658        double  gauss_l1l2l3[3];
     1659
     1660        /* matrices: */
     1661        double L[numgrids];
     1662        double B[2][numgrids];
     1663        double Bprime[2][numgrids];
     1664        double DL[2][2]={0.0};
     1665        double DLprime[2][2]={0.0};
     1666        double DL_scalar;
     1667        double Ke_gg[numdof][numdof]={0.0};
     1668        double Ke_gg1[numdof][numdof]={0.0};
     1669        double Ke_gg2[numdof][numdof]={0.0};
     1670        double Jdettria;
     1671
     1672        /*input parameters for structural analysis (diagnostic): */
     1673        double  vx_list[numgrids]={0.0};
     1674        double  vy_list[numgrids]={0.0};
     1675        double  vx,vy;
     1676        double  dt;
     1677        int     dofs[1]={0};
     1678        int     found;
     1679
     1680        ParameterInputs* inputs=NULL;
     1681
     1682        /*recover pointers: */
     1683        inputs=(ParameterInputs*)vinputs;
     1684
     1685        /*recover extra inputs from users, at current convergence iteration: */
     1686        found=inputs->Recover("vx_average",&vx_list[0],1,dofs,numgrids,(void**)nodes);
     1687        if(!found)ISSMERROR(" could not find vx_average in inputs!");
     1688        found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes);
     1689        if(!found)ISSMERROR(" could not find vy_average in inputs!");
     1690        found=inputs->Recover("dt",&dt);
     1691        if(!found)ISSMERROR(" could not find dt in inputs!");
     1692
     1693        /* Get node coordinates and dof list: */
     1694        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     1695        GetDofList(&doflist[0],&numberofdofspernode);
     1696
     1697        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     1698        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     1699
     1700        /* Start  looping on the number of gaussian points: */
     1701        for (ig=0; ig<num_gauss; ig++){
     1702
     1703                /*Pick up the gaussian point: */
     1704                gauss_weight=*(gauss_weights+ig);
     1705                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     1706                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     1707                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     1708
     1709                /* Get Jacobian determinant: */
     1710                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     1711
     1712                /*Get L matrix: */
     1713                GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     1714
     1715                DL_scalar=gauss_weight*Jdettria;
     1716
     1717                /*  Do the triple product tL*D*L: */
     1718                TripleMultiply( &L[0],1,numdof,1,
     1719                                        &DL_scalar,1,1,0,
     1720                                        &L[0],1,numdof,0,
     1721                                        &Ke_gg1[0][0],0);
     1722
     1723                /*Get B  and B prime matrix: */
     1724                /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
     1725                GetB_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
     1726                GetBPrime_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
     1727
     1728                //Get vx, vy and their derivatives at gauss point
     1729                GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
     1730                GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
     1731
     1732                DL_scalar=-dt*gauss_weight*Jdettria;
     1733
     1734                DLprime[0][0]=DL_scalar*vx;
     1735                DLprime[1][1]=DL_scalar*vy;
     1736
     1737                //Do the triple product tL*D*L.
     1738                TripleMultiply( &B[0][0],2,numdof,1,
     1739                                        &DLprime[0][0],2,2,0,
     1740                                        &Bprime[0][0],2,numdof,0,
     1741                                        &Ke_gg2[0][0],0);
     1742
     1743                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     1744                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg1[i][j];
     1745                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg2[i][j];
     1746
    16061747        } // for (ig=0; ig<num_gauss; ig++)
    16071748
    … …  
    18281969               
    18291970                        CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
    1830                
    18311971                }
    18321972                else ISSMERROR(exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
    … …  
    18391979
    18401980                CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
     1981        }
     1982        else if (analysis_type==Prognostic2AnalysisEnum()){
     1983
     1984                CreatePVectorPrognostic2( pg,inputs,analysis_type,sub_analysis_type);
    18411985        }
    18421986        else if (analysis_type==BalancedthicknessAnalysisEnum()){
    … …  
    23212465void  Tria::CreatePVectorPrognostic(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){
    23222466
     2467
     2468        /* local declarations */
     2469        int             i,j;
     2470
     2471        /* node data: */
     2472        const int    numgrids=3;
     2473        const int    NDOF1=1;
     2474        const int    numdof=NDOF1*numgrids;
     2475        double       xyz_list[numgrids][3];
     2476        int          doflist[numdof];
     2477        int          numberofdofspernode;
     2478
     2479        /* gaussian points: */
     2480        int     num_gauss,ig;
     2481        double* first_gauss_area_coord  =  NULL;
     2482        double* second_gauss_area_coord =  NULL;
     2483        double* third_gauss_area_coord  =  NULL;
     2484        double* gauss_weights           =  NULL;
     2485        double  gauss_weight;
     2486        double  gauss_l1l2l3[3];
     2487
     2488        /* matrix */
     2489        double pe_g[numgrids]={0.0};
     2490        double L[numgrids];
     2491        double Jdettria;
     2492
     2493        /*input parameters for structural analysis (diagnostic): */
     2494        double  accumulation_list[numgrids]={0.0};
     2495        double  accumulation_g;
     2496        double  melting_list[numgrids]={0.0};
     2497        double  melting_g;
     2498        double  thickness_list[numgrids]={0.0};
     2499        double  thickness_g;
     2500        double  dt;
     2501        int     dofs[1]={0};
     2502        int     found=0;
     2503
     2504        ParameterInputs* inputs=NULL;
     2505
     2506        /*recover pointers: */
     2507        inputs=(ParameterInputs*)vinputs;
     2508
     2509        /*recover extra inputs from users, at current convergence iteration: */
     2510        found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes);
     2511        if(!found)ISSMERROR(" could not find accumulation in inputs!");
     2512        found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes);
     2513        if(!found)ISSMERROR(" could not find melting in inputs!");
     2514        found=inputs->Recover("thickness",&thickness_list[0],1,dofs,numgrids,(void**)nodes);
     2515        if(!found)ISSMERROR(" could not find thickness in inputs!");
     2516        found=inputs->Recover("dt",&dt);
     2517        if(!found)ISSMERROR(" could not find dt in inputs!");
     2518
     2519        /* Get node coordinates and dof list: */
     2520        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     2521        GetDofList(&doflist[0],&numberofdofspernode);
     2522
     2523        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     2524        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2525
     2526        /* Start  looping on the number of gaussian points: */
     2527        for (ig=0; ig<num_gauss; ig++){
     2528                /*Pick up the gaussian point: */
     2529                gauss_weight=*(gauss_weights+ig);
     2530                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     2531                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     2532                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     2533
     2534                /* Get Jacobian determinant: */
     2535                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     2536
     2537                /*Get L matrix: */
     2538                GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     2539
     2540                /* Get accumulation, melting and thickness at gauss point */
     2541                GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3);
     2542                GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3);
     2543                GetParameterValue(&thickness_g, &thickness_list[0],gauss_l1l2l3);
     2544
     2545                /* Add value into pe_g: */
     2546                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i];
     2547
     2548        } // for (ig=0; ig<num_gauss; ig++)
     2549
     2550        /*Add pe_g to global matrix Kgg: */
     2551        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
     2552
     2553cleanup_and_return:
     2554        xfree((void**)&first_gauss_area_coord);
     2555        xfree((void**)&second_gauss_area_coord);
     2556        xfree((void**)&third_gauss_area_coord);
     2557        xfree((void**)&gauss_weights);
     2558
     2559}
     2560/*}}}*/
     2561/*FUNCTION CreatePVectorPrognostic2 {{{1*/
     2562void  Tria::CreatePVectorPrognostic2(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){
    23232563
    23242564        /* local declarations */

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

    r3180 r3359  
    125125                void  CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    126126                void  CreatePVectorPrognostic(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);
     127                void  CreateKMatrixPrognostic2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
     128                void  CreatePVectorPrognostic2(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);
    127129                void  CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    128130                void  CreatePVectorBalancedthickness(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);

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

    r2790 r3359  
    2525#include "./Penpair.h"
    2626#include "./Pengrid.h"
     27#include "./Numericalflux.h"
    2728#include "./Param.h"
    2829#include "./Element.h"

  • issm/trunk/src/m/classes/@model/model.m

    r3273 r3359  
    3737        md.nodeconnectivity=NaN;
    3838        md.elementconnectivity=NaN;
     39        md.edges=NaN;
    3940
    4041        %Initial 2d mesh

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

    r3354 r3359  
    280280md.y=bamgmesh_out.Vertices(:,2);
    281281md.elements=bamgmesh_out.Triangles(:,1:3);
     282md.edges=bamgmesh_out.ElementEdges;
    282283md.segments=bamgmesh_out.Segments(:,1:3);
    283284md.segmentmarkers=bamgmesh_out.Segments(:,4);

  • issm/trunk/src/m/classes/public/display/displaymesh.m

    r3273 r3359  
    3333fielddisplay(md,'y','nodes y coordinate');
    3434fielddisplay(md,'z','nodes z coordinate');
     35fielddisplay(md,'edges','edges of the 2d mesh (node1 node2 element1 element2)');
     36
     37disp(sprintf('\n      Properties:'));
     38fielddisplay(md,'type','mesh type');
    3539fielddisplay(md,'numlayers','number of extrusion layers');
    3640fielddisplay(md,'extrusionexponent','exponent for extrusion');
    3741fielddisplay(md,'dof','maximum number of dofs solved');
    38 
    39 disp(sprintf('\n      Properties:'));
    40 fielddisplay(md,'type','mesh type');
    4142fielddisplay(md,'bamg','Geometry and 2d mesh properties (if generated by Bamg)');
    4243fielddisplay(md,'penalties','penalties list');

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

    r3203 r3359  
    8282
    8383WriteData(fid,md.pressureload,'Mat','pressureload');
     84WriteData(fid,md.edges,'Mat','edges');
    8485
    8586WriteData(fid,md.geothermalflux,'Mat','geothermalflux');

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

    r3354 r3359  
    1313        displaystring(md.verbose,'%s',['reading prognostic2 model data']);
    1414        md.analysis_type=Prognostic2AnalysisEnum; models.p=CreateFemModel(md);
    15         error('STOP here');
    1615
    1716        % figure out number of dof: just for information purposes.
    … …  
    2120        displaystring(md.verbose,'\n%s',['setup inputs...']);
    2221        inputs=inputlist;
    23         inputs=add(inputs,'velocity',models.p.parameters.u_g,'doublevec',3,models.p.parameters.numberofnodes);
     22        inputs=add(inputs,'vx',models.p.parameters.vx_g,'doublevec',1,models.p.parameters.numberofnodes);
     23        inputs=add(inputs,'vy',models.p.parameters.vy_g,'doublevec',1,models.p.parameters.numberofnodes);
    2424        inputs=add(inputs,'thickness',models.p.parameters.h_g,'doublevec',1,models.p.parameters.numberofnodes);
    2525        inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.numberofnodes);
    … …  
    2828
    2929        displaystring(md.verbose,'\n%s',['call computational core:']);
    30         results=prognostic2_core(models,inputs,PrognosticAnalysisEnum(),NoneAnalysisEnum());
     30        results=prognostic2_core(models,inputs,Prognostic2AnalysisEnum(),NoneAnalysisEnum());
    3131
    3232        displaystring(md.verbose,'\n%s',['load results...']);

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

    r3354 r3359  
    1212        displaystring(m.parameters.verbose,'\n%s',['depth averaging velocity...']);
    1313        %Take only the first two dofs of m.parameters.u_g
    14         u_g=get(inputs,'velocity',[1 1 0 0]);
    15         u_g=FieldDepthAverage(m.elements,m.nodes,m.loads,m.materials,m.parameters,u_g,'velocity');
    16         inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m.parameters.numberofnodes);
     14        vx_g=get(inputs,'vx',1);
     15        vy_g=get(inputs,'vy',1);
     16        vx_g=FieldDepthAverage(m.elements,m.nodes,m.loads,m.materials,m.parameters,vx_g,'vx');
     17        vy_g=FieldDepthAverage(m.elements,m.nodes,m.loads,m.materials,m.parameters,vy_g,'vy');
     18        inputs=add(inputs,'vx_average',vx_g,'doublevec',1,m.parameters.numberofnodes);
     19        inputs=add(inputs,'vy_average',vy_g,'doublevec',1,m.parameters.numberofnodes);
    1720
    1821        displaystring(m.parameters.verbose,'\n%s',['call computational core:']);
    1922        results.h_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type);
     23        error('STOP here');
    2024
    2125        displaystring(m.parameters.verbose,'\n%s',['extrude computed thickness on all layers:']);
Note: See TracChangeset for help on using the changeset viewer.