Changeset 3429


Ignore:
Timestamp:
04/07/10 16:12:51 (15 years ago)
Author:
Eric.Larour
Message:

New simplified ModelProcessor.

Location:
issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r3332 r3429  
    1515
    1616        int i,j;
    17         int count;
     17        int count=0;
    1818       
    1919        DataSet* constraints = NULL;
    2020        Spc*    spc  = NULL;
    2121        Rgb*    rgb  = NULL;
     22        int     node1,node2;
    2223
    23         /*spc intermediary data: */
    24         int spc_sid;
    25         int spc_node;
    26         int spc_dof;
    27         double spc_value;
    28 
    29         /*rgb constructor data: */
    30         int rgb_id;
    31         int rgb_dof;
    32         int rgb_nodeid1;
    33         int rgb_nodeid2;
    34         int grid1,grid2;
    35        
    36         double* spcvelocity=NULL;
    37         double* gridonhutter=NULL;
     24        /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
     25        if (!iomodel->ismacayealpattyn)goto cleanup_and_return;
    3826       
    3927        /*Create constraints: */
    4028        constraints = new DataSet(ConstraintsEnum());
    4129
    42         /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
    43         if (!iomodel->ismacayealpattyn)goto cleanup_and_return;
     30        /*Spcs: fetch data: */
     31        IoModelFetchData(&iomodel->spcvelocity,NULL,NULL,iomodel_handle,"spcvelocity");
     32        IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");
    4433
    45         /*Fetch data: */
    46         IoModelFetchData(&spcvelocity,NULL,NULL,iomodel_handle,"spcvelocity");
    47         IoModelFetchData(&gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");
     34        count=1; //matlab indexing
     35        /*Create spcs from x,y,z, as well as the spc values on those spcs: */
     36        for (i=0;i<iomodel->numberofvertices;i++){
     37                if(iomodel->my_vertices[i]){
    4838
    49         count=0;
     39                        if ((int)iomodel->spcvelocity[6*i+0] | (int)iomodel->gridonhutter[i])
     40                                constraints->AddObject(new Spc(count,i+1,1,*(iomodel->spcvelocity+6*i+3)/iomodel->yts)); //add count'th spc, on node i+1, setting dof 1 to vx.
    5041
    51         /*Create spcs from x,y,z, as well as the spc values on those spcs: */
    52         for (i=0;i<iomodel->numberofnodes;i++){
    53         #ifdef _PARALLEL_
    54         /*keep only this partition's nodes:*/
    55         if((iomodel->my_grids[i]==1)){
    56         #endif
     42                        count++;
     43                       
     44                        if ((int)iomodel->spcvelocity[6*i+1] | (int)iomodel->gridonhutter[i])
     45                                constraints->AddObject(new Spc(count,i+1,2,*(iomodel->spcvelocity+6*i+4)/iomodel->yts)); //add count'th spc, on node i+1, setting dof 2 to vy
    5746
    58                 if ((int)spcvelocity[6*i+0] | (int)gridonhutter[i]){
    59        
    60                         /*This grid needs to be spc'd to vx_obs:*/
    61 
    62                         spc_sid=count;
    63                         spc_node=i+1;
    64                         spc_dof=1; //we enforce first x translation degree of freedom
    65                         spc_value=*(spcvelocity+6*i+3)/iomodel->yts;
    66 
    67                         spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
    68                         constraints->AddObject(spc);
    6947                        count++;
    7048                }
     49        }
     50         
     51        /*Free data: */
     52        xfree((void**)&iomodel->spcvelocity);
     53        xfree((void**)&iomodel->gridonhutter);
    7154
    72                 if ((int)spcvelocity[6*i+1] | (int)gridonhutter[i]){
    73 
    74                         /*This grid needs to be spc'd to vy_obs:*/
    75 
    76                         spc_sid=count;
    77                         spc_node=i+1;
    78                         spc_dof=2; //we enforce first y translation degree of freedom
    79                         spc_value=*(spcvelocity+6*i+4)/iomodel->yts;
    80                        
    81                         spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
    82                         constraints->AddObject(spc);
    83                         count++;
    84                 }
    85 
    86         #ifdef _PARALLEL_
    87         } //if((my_grids[i]==1))
    88         #endif
    89         }
    90 
    91         /*Create penalties loads for linking of collapsed formulation grids with non collapsed grids: */
    92 
    93         /*First fetch penalties: */
    94         if (strcmp(iomodel->meshtype,"3d")==0){
    95                 IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");
    96         }
    97 
    98         count=0;
     55       
     56        /*penalty loads: */
     57        IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");
     58       
     59        count=1; //matlab indexing
    9960        if(iomodel->numpenalties){
    10061
    101                 /*First deal with internal grids: */
    10262                for (i=0;i<iomodel->numpenalties;i++){
    103                         if (iomodel->penaltypartitioning[i]>=0){ //this penalty belongs to at least this CPU
    10463
    105                                 for(j=1;j<iomodel->numlayers;j++){
    106                                         /*We are pairing grids along a vertical profile.*/
    107                                         grid1=(int)*(iomodel->penalties+iomodel->numlayers*i);
    108                                         grid2=(int)*(iomodel->penalties+iomodel->numlayers*i+j);
     64                        for(j=1;j<iomodel->numlayers;j++){
     65       
     66                                /*We are pairing nodes along a vertical profile.*/
     67                                node1=(int)*(iomodel->penalties+iomodel->numlayers*i);
     68                                node2=(int)*(iomodel->penalties+iomodel->numlayers*i+j);
    10969
    110                                         rgb_id=count+1; //matlab indexing
    111                                         rgb_dof=1;
    112                                         rgb_nodeid1=grid1;
    113                                         rgb_nodeid2=grid2;
    114                                         rgb= new Rgb(rgb_id,rgb_nodeid1,rgb_nodeid2,rgb_dof);
    115                                         constraints->AddObject(rgb);
    116                                         count++;
     70                                constraints->AddObject(new Rgb(count,node1,node2,1));  //add count'th Rgb on dof 1 between node1 and node2
     71                               
     72                                count++;
     73                               
     74                                constraints->AddObject(new Rgb(count,node1,node2,2));  //add count'th Rgb on dof 1 between node1 and node2
    11775
    118                                         rgb_id=count+1; //matlab indexing
    119                                         rgb_dof=2;
    120                                         rgb_nodeid1=grid1;
    121                                         rgb_nodeid2=grid2;
    122                                         rgb= new Rgb(rgb_id,rgb_nodeid1,rgb_nodeid2,rgb_dof);
    123                                         constraints->AddObject(rgb);
    124                                         count++;
    125                                 } //for ( i=0; i<numpenalties; i++)
    12676                        }
    12777                }
    … …  
    13585        constraints->Presort();
    13686       
    137         /*Free data: */
    138         xfree((void**)&spcvelocity);
    139         xfree((void**)&gridonhutter);
    140        
    14187        cleanup_and_return:
    14288        /*Assign output pointer: */

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

    r3423 r3429  
    2323        DataSet*    materials = NULL;
    2424       
     25        /*Objects: */
     26        Node*       node   = NULL;
     27        Vertex*     vertex = NULL;
     28        Tria*       tria = NULL;
     29        Penta*      penta = NULL;
     30        Matice*     matice  = NULL;
     31        Matpar*     matpar  = NULL;
     32
    2533        /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
    2634        if (!iomodel->ismacayealpattyn)goto cleanup_and_return;

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

    r3359 r3429  
    1313void    CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1414
    15         int i;
    16         int element;
    17 
    18         extern int my_rank;
    19         extern int num_procs;
    20        
    2115        DataSet*    loads    = NULL;
    2216        Icefront*   icefront = NULL;
    … …  
    2418
    2519        int segment_width;
    26         int element_type;
    27         int i1,i2,i3,i4;
    28         int i0,range;
     20        int element;
    2921
    30                
    31         /*icefront intermediary data: */
    32         char icefront_type[ICEFRONTSTRING];
    33         int  icefront_fill;
    34         int  icefront_element_type;
    35         int  icefront_sid;
    36         int  icefront_eid;
    37         int  icefront_mparid;
    38         int  icefront_node_ids[MAX_ICEFRONT_GRIDS];
    39         double icefront_h[MAX_ICEFRONT_GRIDS];
    40         double  icefront_b[MAX_ICEFRONT_GRIDS];
    41 
    42         /*rifts: */
    43         char riftfront_type[RIFTFRONTSTRING];
    44         int  riftfront_id;
    45         int  riftfront_node_ids[2];
    46         int  riftfront_mparid;
    47         double riftfront_h[2];
    48         double riftfront_b[2];
    49         double riftfront_s[2];
    50         double riftfront_normal[2];
    51         double riftfront_length;
    52         int    riftfront_fill;
    53         double riftfront_friction;
    54         double riftfront_fraction;
    55         double riftfront_fractionincrement;
    56         bool   riftfront_shelf;
    57         double riftfront_penalty_offset;
    58         int riftfront_penalty_lock;
    59         bool riftfront_active;
    60         bool riftfront_frozen;
    61         int  riftfront_counter;
    62         bool riftfront_prestable;
    63         int el1,el2;
    64         int grid1,grid2;
    65         double normal[2];
    66         double length;
    67         int    fill;
    68         double friction;
    69         double fraction;
    70         double fractionincrement;
    71        
    7222        /*Create loads: */
    7323        loads   = new DataSet(LoadsEnum());
    … …  
    8636        for (i=0;i<iomodel->numberofpressureloads;i++){
    8737               
    88                 if (strcmp(iomodel->meshtype,"2d")==0){
    89                         segment_width=4;
    90                         element_type=TriaEnum();
    91                 }
    92                 else{
    93                         segment_width=6;
    94                         element_type=PentaEnum();
    95                 }
    96 
     38                /*Retrieve element to which this icefront belongs: */
     39                if (strcmp(iomodel->meshtype,"2d")==0) segment_width=4;
     40                else segment_width=6;
    9741                element=(int)(*(iomodel->pressureload+segment_width*i+segment_width-2)-1); //element is in the penultimate column (grid1 grid2 ... elem fill)
    9842
    99                 #ifdef _PARALLEL_
    100                 if (iomodel->epart[element]!=my_rank){
    101                         /*This load does not belong to this cluster node, as it references an element
    102                          *that does not belong to this node's partition. Drop this 'i':*/
    103                         continue;
    104                 }
     43                /*Now, if this element is not in the partition, pass: */
     44                if(!iomodel->my_elements[i])continue;
     45               
     46                /*Do not create ice front if Hutter or Stokes elements*/
     47                if ((int)*(iomodel->elements_type+2*element+0)==(HutterFormulationEnum() || StokesFormulationEnum()))continue;
    10548
    106                 #endif
    107        
    108                 /*Do not create ice front if Hutter or Stokes elements*/
    109                 if ((int)*(iomodel->elements_type+2*element+0)==(HutterFormulationEnum() || StokesFormulationEnum())){
    110                         continue;
    111                 }
    112 
    113                 icefront_mparid=iomodel->numberofelements+1; //matlab indexing
    114                 icefront_sid=i+1; //matlab indexing
    115                 icefront_fill=(int)*(iomodel->pressureload+segment_width*i+segment_width-1);
    116                 icefront_eid=(int) *(iomodel->pressureload+segment_width*i+segment_width-2); //matlab indexing
    117                 icefront_element_type=element_type;
    118 
    119                 i1=(int)*(iomodel->pressureload+segment_width*i+0);
    120                 i2=(int)*(iomodel->pressureload+segment_width*i+1);
    121                        
    122                 icefront_node_ids[0]=i1;
    123                 icefront_node_ids[1]=i2;
    124                
    125                 icefront_h[0]=iomodel->thickness[i1-1];
    126                 icefront_h[1]=iomodel->thickness[i2-1];
    127 
    128                 icefront_b[0]=iomodel->bed[i1-1];
    129                 icefront_b[1]=iomodel->bed[i2-1];
    130                
    131                 if ((int)*(iomodel->elements_type+2*element+0)==MacAyealFormulationEnum()){ //this is a collapsed 3d element, icefront will be 2d
    132                         strcpy(icefront_type,"segment");
    133                 }
    134                 else if ((int)*(iomodel->elements_type+2*element+0)==PattynFormulationEnum()){ //this is a real 3d element, icefront will be 3d.
    135                         strcpy(icefront_type,"quad");
    136                         i3=(int)*(iomodel->pressureload+segment_width*i+2);
    137                         i4=(int)*(iomodel->pressureload+segment_width*i+3);
    138                         icefront_node_ids[2]=i3;
    139                         icefront_node_ids[3]=i4;
    140                        
    141                         icefront_h[2]=iomodel->thickness[i3-1];
    142                         icefront_h[3]=iomodel->thickness[i4-1];
    143 
    144                         icefront_b[2]=iomodel->bed[i3-1];
    145                         icefront_b[3]=iomodel->bed[i4-1];
    146                 }
    147                 else{
    148                         ISSMERROR(exprintf(" element type %i not supported yet",(int)*(iomodel->elements_type+2*element+0)));
    149                 }
    150 
    151                 icefront = new Icefront(icefront_type,icefront_fill,icefront_sid,icefront_mparid,icefront_eid,icefront_element_type,icefront_node_ids,icefront_h,icefront_b);
    152                
    153                 loads->AddObject(icefront);
     49                /*Create and  add load: */
     50                loads->AddObject(new Icefront(i,iomodel));
    15451
    15552        }
    … …  
    16057        xfree((void**)&iomodel->bed);
    16158
    162         /*Create penpair loads also for rift grids, so that they don't penetrate one another, if needed: */
    163         /*First fetch rifts: */
     59        /*Create Riffront loads for rifts: */
    16460        IoModelFetchData(&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo");
    16561        IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
    … …  
    17167                for(i=0;i<iomodel->numrifts;i++){
    17268                               
    173                         el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2);
    174                         #ifdef _PARALLEL_
    175                         if (iomodel->epart[el1-1]!=my_rank){
    176                                 /*This load does not belong to this cluster node, as it references an element
    177                                  *that does not belong to this node's partition. Drop this 'j':*/
    178                                 continue;
     69                        if(iomodel->my_elements[(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2)]){
     70
     71                                loads->AddObject(new Riftfront(i,iomodel));
    17972                        }
    180                         #endif
    181 
    182                         /*Ok, retrieve all the data needed to add a penalty between the two grids: */
    183                         el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+3);
    184 
    185                         grid1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+0);
    186                         grid2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+1);
    187                        
    188                         normal[0]=*(iomodel->riftinfo+RIFTINFOSIZE*i+4);
    189                         normal[1]=*(iomodel->riftinfo+RIFTINFOSIZE*i+5);
    190                         length=*(iomodel->riftinfo+RIFTINFOSIZE*i+6);
    191                        
    192                         fill = (int)*(iomodel->riftinfo+RIFTINFOSIZE*i+7);
    193                         friction=*(iomodel->riftinfo+RIFTINFOSIZE*i+8);
    194                         fraction=*(iomodel->riftinfo+RIFTINFOSIZE*i+9);
    195                         fractionincrement=*(iomodel->riftinfo+RIFTINFOSIZE*i+10);
    196        
    197                         strcpy(riftfront_type,"2d");
    198                         riftfront_id=i+1; //matlab indexing
    199                         riftfront_node_ids[0]=grid1;
    200                         riftfront_node_ids[1]=grid2;
    201                         riftfront_mparid=iomodel->numberofelements+1; //matlab indexing
    202 
    203                         riftfront_h[0]=iomodel->thickness[grid1-1];
    204                         riftfront_h[1]=iomodel->thickness[grid2-1];
    205 
    206                         riftfront_b[0]=iomodel->bed[grid1-1];
    207                         riftfront_b[1]=iomodel->bed[grid2-1];
    208 
    209                         riftfront_s[0]=iomodel->surface[grid1-1];
    210                         riftfront_s[1]=iomodel->surface[grid2-1];
    211 
    212                         riftfront_normal[0]=normal[0];
    213                         riftfront_normal[1]=normal[1];
    214                         riftfront_length=length;
    215                        
    216                         riftfront_fill=fill;
    217                         riftfront_friction=friction;
    218                         riftfront_fraction=fraction;
    219                         riftfront_fractionincrement=fractionincrement;
    220                         riftfront_shelf=(bool)iomodel->gridoniceshelf[grid1-1];
    221 
    222                         riftfront_penalty_offset=iomodel->penalty_offset;
    223                         riftfront_penalty_lock=iomodel->penalty_lock;
    224 
    225                         riftfront_active=0;
    226                         riftfront_frozen=0;
    227                         riftfront_counter=0;
    228                         riftfront_prestable=0;
    229                        
    230                         riftfront=new Riftfront(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction, riftfront_fractionincrement, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_frozen,riftfront_counter,riftfront_prestable,riftfront_shelf);
    231 
    232                         loads->AddObject(riftfront);
    23373                }
    23474        }
    235 
    236 
     75                               
     76       
    23777        /*free ressources: */
    23878        xfree((void**)&iomodel->riftinfo);

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

    r3332 r3429  
    3535        IoModelFetchData(&vz,NULL,NULL,iomodel_handle,"vz");
    3636
    37         ug=(double*)xcalloc(iomodel->numberofnodes*3,sizeof(double));
     37        ug=(double*)xcalloc(iomodel->numberofvertices*3,sizeof(double));
    3838
    39         if(vx) for(i=0;i<iomodel->numberofnodes;i++)ug[3*i+0]=vx[i]/iomodel->yts;
    40         if(vy) for(i=0;i<iomodel->numberofnodes;i++)ug[3*i+1]=vy[i]/iomodel->yts;
    41         if(vz) for(i=0;i<iomodel->numberofnodes;i++)ug[3*i+2]=vz[i]/iomodel->yts;
     39        if(vx) for(i=0;i<iomodel->numberofvertices;i++)ug[3*i+0]=vx[i]/iomodel->yts;
     40        if(vy) for(i=0;i<iomodel->numberofvertices;i++)ug[3*i+1]=vy[i]/iomodel->yts;
     41        if(vz) for(i=0;i<iomodel->numberofvertices;i++)ug[3*i+2]=vz[i]/iomodel->yts;
    4242
    4343        count++;
    4444        param= new Param(count,"u_g",DOUBLEVEC);
    45         param->SetDoubleVec(ug,3*iomodel->numberofnodes,3);
     45        param->SetDoubleVec(ug,3*iomodel->numberofvertices,3);
    4646        parameters->AddObject(param);
    4747
Note: See TracChangeset for help on using the changeset viewer.