Changeset 23960

Timestamp:

05/29/19 22:26:27 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on channels

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyGlaDSAnalysis.cpp (modified) (1 diff)
• classes/Loads/Channel.cpp (modified) (6 diffs)
• classes/Loads/Channel.h (modified) (1 diff)
• classes/Loads/Moulin.cpp (modified) (1 diff)

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

@@ -26 +26 @@
         if(hydrology_model!=HydrologyGlaDSEnum) return;
 
-        /*No extra loads*/
+        /*Add channels?*/
+        bool ischannels;
+        iomodel->FindConstant(&ischannels,"md.hydrology.ischannels");
+        if(ischannels){
+                /*Get faces (edges in 2d)*/
+                CreateFaces(iomodel);
+                for(int i=0;i<iomodel->numberoffaces;i++){
+                        /*Get left and right elements*/
+                        int element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
+
+                        /*Now, if this element is not in the partition*/
+                        if(!iomodel->my_elements[element]) continue;
+
+                        /* Add load */
+                        loads->AddObject(new Channel(i+1,i,i,iomodel));
+                }
+
+                /*Free data: */
+        }
 }/*}}}*/
 void HydrologyGlaDSAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/

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

@@ -28 +28 @@
 Channel::Channel(int channel_id,int i,int index,IoModel* iomodel){/*{{{*/
 
-        /* Intermediary */
-        int  j;
-        int  pos1,pos2,pos3,pos4;
-        int  num_nodes;
-
-        /*channel constructor data: */
-        int   channel_elem_ids[2];
-        int   channel_vertex_ids[2];
-        int   channel_node_ids[4];
-        int   channel_type;
-
-        /*Get edge*/
+        this->id=channel_id;
+        this->parameters = NULL;
+        this->element    = NULL;
+        this->nodes      = NULL;
+
+        /*Get edge info*/
         int i1 = iomodel->faces[4*index+0];
         int i2 = iomodel->faces[4*index+1];
@@ -45 +39 @@
         int e2 = iomodel->faces[4*index+3];
 
-        /*First, see wether this is an internal or boundary edge (if e2=-1)*/
-        if(e2==-1){
-                /* Boundary edge, only one element */
-                num_nodes=2;
-                channel_type=BoundaryEnum;
-                channel_elem_ids[0]=e1;
-        }
-        else{
-                /* internal edge: connected to 2 elements */
-                 num_nodes=4;
-                channel_type=InternalEnum;
-                channel_elem_ids[0]=e1;
-                channel_elem_ids[1]=e2;
-        }
-
-        /*1: Get vertices ids*/
+        /*Set Element hook (4th column may be -1 for boundary edges)*/
+        this->helement  = new Hook(&e1,1);
+
+        /*Set Vertices hooks (4th column may be -1 for boundary edges)*/
+        int channel_vertex_ids[2];
         channel_vertex_ids[0]=i1;
         channel_vertex_ids[1]=i2;
-
-        /*2: Get node ids*/
-        if (channel_type==InternalEnum){
-
-                /*Now, we must get the nodes of the 4 nodes located on the edge*/
-
-                /*2: Get the column where these ids are located in the index*/
-                pos1=pos2=pos3=pos4=UNDEF;
-                for(j=0;j<3;j++){
-                        if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
-                        if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
-                        if(iomodel->elements[3*(e2-1)+j]==i1) pos3=j+1;
-                        if(iomodel->elements[3*(e2-1)+j]==i2) pos4=j+1;
-                }
-                _assert_(pos1!=UNDEF && pos2!=UNDEF && pos3!=UNDEF && pos4!=UNDEF);
-
-                /*3: We have the id of the elements and the position of the vertices in the index
-                 * we can compute their dofs!*/
-                channel_node_ids[0]=3*(e1-1)+pos1;
-                channel_node_ids[1]=3*(e1-1)+pos2;
-                channel_node_ids[2]=3*(e2-1)+pos3;
-                channel_node_ids[3]=3*(e2-1)+pos4;
-        }
-        else{
-
-                /*2: Get the column where these ids are located in the index*/
-                pos1=pos2=UNDEF;
-                for(j=0;j<3;j++){
-                        if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
-                        if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
-                }
-                _assert_(pos1!=UNDEF && pos2!=UNDEF);
-
-                /*3: We have the id of the elements and the position of the vertices in the index
-                 * we can compute their dofs!*/
-                channel_node_ids[0]=3*(e1-1)+pos1;
-                channel_node_ids[1]=3*(e1-1)+pos2;
-        }
-
-        /*Ok, we have everything to build the object: */
-        this->id=channel_id;
-
-        /*Hooks: */
-        this->hnodes    =new Hook(channel_node_ids,num_nodes);
         this->hvertices =new Hook(&channel_vertex_ids[0],2);
-        this->helement  =new Hook(channel_elem_ids,1); // take only the first element for now
-
-        //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
-        this->parameters=NULL;
-        this->element=NULL;
-        this->nodes=NULL;
-}
-/*}}}*/
+
+        /*Set Nodes hooks (!! Assumes P1 CG)*/
+        int channel_node_ids[2];
+        channel_node_ids[0]=i1;
+        channel_node_ids[1]=i2;
+        this->hnodes=new Hook(&channel_node_ids[0],2);
+
+}/*}}}*/
 Channel::~Channel(){/*{{{*/
         this->parameters=NULL;
@@ -121 +60 @@
         delete hnodes;
         delete hvertices;
-}
-/*}}}*/
+}/*}}}*/
 
 /*Object virtual functions definitions:*/
@@ -234 +172 @@
         this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 
-        if(analysis_type != HydrologyGlaDSAnalysisEnum) _error_("Analysis not supported!!");
-        _error_("STOP");
+        switch(analysis_type){
+                case HydrologyGlaDSAnalysisEnum:
+                        Ke = this->CreateKMatrixHydrologyGlaDS();
+                        break;
+                default:
+                        _error_("Don't know why we should be here");
+        }
 
         /*Add to global matrix*/
@@ -251 +194 @@
         int analysis_type;
         this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-        if(analysis_type != HydrologyGlaDSAnalysisEnum) _error_("Analysis not supported!!");
-        _error_("STOP");
+
+        switch(analysis_type){
+                case HydrologyGlaDSAnalysisEnum:
+                        pe = this->CreatePVectorHydrologyGlaDS();
+                        break;
+                default:
+                        _error_("Don't know why we should be here");
+        }
 
         /*Add to global matrix*/
@@ -372 +321 @@
 }
 /*}}}*/
+
+/*Channel specific functions*/
+ElementVector* Channel::CreatePVectorHydrologyGlaDS(void){/*{{{*/
+
+        _error_("not implemented :( ");
+
+        /*Initialize Element matrix and return if necessary*/
+        Tria*  tria=(Tria*)element;
+        if(!tria->IsIceInElement()) return NULL;
+        ElementVector* pe=new ElementVector(nodes,NUMNODES,this->parameters);
+
+        /*Clean up and return*/
+        return pe;
+}
+/*}}}*/
+ElementMatrix* Channel::CreateKMatrixHydrologyGlaDS(void){/*{{{*/
+
+        _error_("not implemented :( ");
+
+        /*Initialize Element matrix and return if necessary*/
+        Tria*  tria=(Tria*)element;
+        if(!tria->IsIceInElement()) return NULL;
+        ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES,this->parameters);
+
+        /*Clean up and return*/
+        return Ke;
+}
+/*}}}*/

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

@@ -72 +72 @@
                 /*Channel management:{{{*/
                 void           GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
+                ElementVector* CreatePVectorHydrologyGlaDS(void);
+                ElementMatrix* CreateKMatrixHydrologyGlaDS(void);
                 /*}}}*/
 

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

@@ -159 +159 @@
 
         switch(analysis_type){
-
-        case HydrologyShaktiAnalysisEnum:
-                pe = this->CreatePVectorHydrologyShakti();
-                break;
-        case HydrologyDCInefficientAnalysisEnum:
-                pe = CreatePVectorHydrologyDCInefficient();
-                break;
-        case HydrologyDCEfficientAnalysisEnum:
-                pe = CreatePVectorHydrologyDCEfficient();
-                break;
-        default:
-                _error_("Don't know why we should be here");
-                /*No loads applied, do nothing: */
-                return;
+                case HydrologyShaktiAnalysisEnum:
+                        pe = this->CreatePVectorHydrologyShakti();
+                        break;
+                case HydrologyDCInefficientAnalysisEnum:
+                        pe = CreatePVectorHydrologyDCInefficient();
+                        break;
+                case HydrologyDCEfficientAnalysisEnum:
+                        pe = CreatePVectorHydrologyDCEfficient();
+                        break;
+                default:
+                        _error_("Don't know why we should be here");
+                        /*No loads applied, do nothing: */
+                        return;
         }
         if(pe){