source: issm/oecreview/Archive/23390-24306/ISSM-23959-23960.diff@ 28275

Last change on this file since 28275 was 24307, checked in by Mathieu Morlighem, 5 years ago

NEW: adding Archive/23390-24306
File size: 8.0 KB

  • ../trunk-jpl/src/c/classes/Loads/Channel.cpp

     
    2727/*}}}*/
    2828Channel::Channel(int channel_id,int i,int index,IoModel* iomodel){/*{{{*/
    2929
    30         /* Intermediary */
    31         int  j;
    32         int  pos1,pos2,pos3,pos4;
    33         int  num_nodes;
     30        this->id=channel_id;
     31        this->parameters = NULL;
     32        this->element    = NULL;
     33        this->nodes      = NULL;
    3434
    35         /*channel constructor data: */
    36         int   channel_elem_ids[2];
    37         int   channel_vertex_ids[2];
    38         int   channel_node_ids[4];
    39         int   channel_type;
    40 
    41         /*Get edge*/
     35        /*Get edge info*/
    4236        int i1 = iomodel->faces[4*index+0];
    4337        int i2 = iomodel->faces[4*index+1];
    4438        int e1 = iomodel->faces[4*index+2];
    4539        int e2 = iomodel->faces[4*index+3];
    4640
    47         /*First, see wether this is an internal or boundary edge (if e2=-1)*/
    48         if(e2==-1){
    49                 /* Boundary edge, only one element */
    50                 num_nodes=2;
    51                 channel_type=BoundaryEnum;
    52                 channel_elem_ids[0]=e1;
    53         }
    54         else{
    55                 /* internal edge: connected to 2 elements */
    56                  num_nodes=4;
    57                 channel_type=InternalEnum;
    58                 channel_elem_ids[0]=e1;
    59                 channel_elem_ids[1]=e2;
    60         }
     41        /*Set Element hook (4th column may be -1 for boundary edges)*/
     42        this->helement  = new Hook(&e1,1);
    6143
    62         /*1: Get vertices ids*/
     44        /*Set Vertices hooks (4th column may be -1 for boundary edges)*/
     45        int channel_vertex_ids[2];
    6346        channel_vertex_ids[0]=i1;
    6447        channel_vertex_ids[1]=i2;
     48        this->hvertices =new Hook(&channel_vertex_ids[0],2);
    6549
    66         /*2: Get node ids*/
    67         if (channel_type==InternalEnum){
     50        /*Set Nodes hooks (!! Assumes P1 CG)*/
     51        int channel_node_ids[2];
     52        channel_node_ids[0]=i1;
     53        channel_node_ids[1]=i2;
     54        this->hnodes=new Hook(&channel_node_ids[0],2);
    6855
    69                 /*Now, we must get the nodes of the 4 nodes located on the edge*/
    70 
    71                 /*2: Get the column where these ids are located in the index*/
    72                 pos1=pos2=pos3=pos4=UNDEF;
    73                 for(j=0;j<3;j++){
    74                         if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
    75                         if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
    76                         if(iomodel->elements[3*(e2-1)+j]==i1) pos3=j+1;
    77                         if(iomodel->elements[3*(e2-1)+j]==i2) pos4=j+1;
    78                 }
    79                 _assert_(pos1!=UNDEF && pos2!=UNDEF && pos3!=UNDEF && pos4!=UNDEF);
    80 
    81                 /*3: We have the id of the elements and the position of the vertices in the index
    82                  * we can compute their dofs!*/
    83                 channel_node_ids[0]=3*(e1-1)+pos1;
    84                 channel_node_ids[1]=3*(e1-1)+pos2;
    85                 channel_node_ids[2]=3*(e2-1)+pos3;
    86                 channel_node_ids[3]=3*(e2-1)+pos4;
    87         }
    88         else{
    89 
    90                 /*2: Get the column where these ids are located in the index*/
    91                 pos1=pos2=UNDEF;
    92                 for(j=0;j<3;j++){
    93                         if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
    94                         if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
    95                 }
    96                 _assert_(pos1!=UNDEF && pos2!=UNDEF);
    97 
    98                 /*3: We have the id of the elements and the position of the vertices in the index
    99                  * we can compute their dofs!*/
    100                 channel_node_ids[0]=3*(e1-1)+pos1;
    101                 channel_node_ids[1]=3*(e1-1)+pos2;
    102         }
    103 
    104         /*Ok, we have everything to build the object: */
    105         this->id=channel_id;
    106 
    107         /*Hooks: */
    108         this->hnodes    =new Hook(channel_node_ids,num_nodes);
    109         this->hvertices =new Hook(&channel_vertex_ids[0],2);
    110         this->helement  =new Hook(channel_elem_ids,1); // take only the first element for now
    111 
    112         //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
    113         this->parameters=NULL;
    114         this->element=NULL;
    115         this->nodes=NULL;
    116 }
    117 /*}}}*/
     56}/*}}}*/
    11857Channel::~Channel(){/*{{{*/
    11958        this->parameters=NULL;
    12059        delete helement;
    12160        delete hnodes;
    12261        delete hvertices;
    123 }
    124 /*}}}*/
     62}/*}}}*/
    12563
    12664/*Object virtual functions definitions:*/
    12765Object* Channel::copy() {/*{{{*/
    … …  
    233171        int analysis_type;
    234172        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    235173
    236         if(analysis_type != HydrologyGlaDSAnalysisEnum) _error_("Analysis not supported!!");
    237         _error_("STOP");
     174        switch(analysis_type){
     175                case HydrologyGlaDSAnalysisEnum:
     176                        Ke = this->CreateKMatrixHydrologyGlaDS();
     177                        break;
     178                default:
     179                        _error_("Don't know why we should be here");
     180        }
    238181
    239182        /*Add to global matrix*/
    240183        if(Ke){
    … …  
    250193        ElementVector* pe=NULL;
    251194        int analysis_type;
    252195        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    253         if(analysis_type != HydrologyGlaDSAnalysisEnum) _error_("Analysis not supported!!");
    254         _error_("STOP");
    255196
     197        switch(analysis_type){
     198                case HydrologyGlaDSAnalysisEnum:
     199                        pe = this->CreatePVectorHydrologyGlaDS();
     200                        break;
     201                default:
     202                        _error_("Don't know why we should be here");
     203        }
     204
    256205        /*Add to global matrix*/
    257206        if(pe){
    258207                pe->AddToGlobal(pf);
    … …  
    371320        *po_nz=o_nz;
    372321}
    373322/*}}}*/
     323
     324/*Channel specific functions*/
     325ElementVector* Channel::CreatePVectorHydrologyGlaDS(void){/*{{{*/
     326
     327        _error_("not implemented :( ");
     328
     329        /*Initialize Element matrix and return if necessary*/
     330        Tria*  tria=(Tria*)element;
     331        if(!tria->IsIceInElement()) return NULL;
     332        ElementVector* pe=new ElementVector(nodes,NUMNODES,this->parameters);
     333
     334        /*Clean up and return*/
     335        return pe;
     336}
     337/*}}}*/
     338ElementMatrix* Channel::CreateKMatrixHydrologyGlaDS(void){/*{{{*/
     339
     340        _error_("not implemented :( ");
     341
     342        /*Initialize Element matrix and return if necessary*/
     343        Tria*  tria=(Tria*)element;
     344        if(!tria->IsIceInElement()) return NULL;
     345        ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES,this->parameters);
     346
     347        /*Clean up and return*/
     348        return Ke;
     349}
     350/*}}}*/

  • ../trunk-jpl/src/c/classes/Loads/Channel.h

     
    7171                /*}}}*/
    7272                /*Channel management:{{{*/
    7373                void           GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
     74                ElementVector* CreatePVectorHydrologyGlaDS(void);
     75                ElementMatrix* CreateKMatrixHydrologyGlaDS(void);
    7476                /*}}}*/
    7577
    7678};

  • ../trunk-jpl/src/c/classes/Loads/Moulin.cpp

     
    158158        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    159159
    160160        switch(analysis_type){
    161 
    162         case HydrologyShaktiAnalysisEnum:
    163                 pe = this->CreatePVectorHydrologyShakti();
    164                 break;
    165         case HydrologyDCInefficientAnalysisEnum:
    166                 pe = CreatePVectorHydrologyDCInefficient();
    167                 break;
    168         case HydrologyDCEfficientAnalysisEnum:
    169                 pe = CreatePVectorHydrologyDCEfficient();
    170                 break;
    171         default:
    172                 _error_("Don't know why we should be here");
    173                 /*No loads applied, do nothing: */
    174                 return;
     161                case HydrologyShaktiAnalysisEnum:
     162                        pe = this->CreatePVectorHydrologyShakti();
     163                        break;
     164                case HydrologyDCInefficientAnalysisEnum:
     165                        pe = CreatePVectorHydrologyDCInefficient();
     166                        break;
     167                case HydrologyDCEfficientAnalysisEnum:
     168                        pe = CreatePVectorHydrologyDCEfficient();
     169                        break;
     170                default:
     171                        _error_("Don't know why we should be here");
     172                        /*No loads applied, do nothing: */
     173                        return;
    175174        }
    176175        if(pe){
    177176                pe->AddToGlobal(pf);

  • ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp

     
    2525        /*Now, do we really want GlaDS?*/
    2626        if(hydrology_model!=HydrologyGlaDSEnum) return;
    2727
    28         /*No extra loads*/
     28        /*Add channels?*/
     29        bool ischannels;
     30        iomodel->FindConstant(&ischannels,"md.hydrology.ischannels");
     31        if(ischannels){
     32                /*Get faces (edges in 2d)*/
     33                CreateFaces(iomodel);
     34                for(int i=0;i<iomodel->numberoffaces;i++){
     35                        /*Get left and right elements*/
     36                        int element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
     37
     38                        /*Now, if this element is not in the partition*/
     39                        if(!iomodel->my_elements[element]) continue;
     40
     41                        /* Add load */
     42                        loads->AddObject(new Channel(i+1,i,i,iomodel));
     43                }
     44
     45                /*Free data: */
     46        }
    2947}/*}}}*/
    3048void HydrologyGlaDSAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
    3149
Note: See TracBrowser for help on using the repository browser.