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

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

NEW: adding Archive/23390-24306
File size: 8.0 KB
RevLine 
[24307]1Index: ../trunk-jpl/src/c/classes/Loads/Channel.cpp
2===================================================================
3--- ../trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 23959)
4+++ ../trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 23960)
5@@ -27,101 +27,39 @@
6 /*}}}*/
7 Channel::Channel(int channel_id,int i,int index,IoModel* iomodel){/*{{{*/
8
9- /* Intermediary */
10- int j;
11- int pos1,pos2,pos3,pos4;
12- int num_nodes;
13+ this->id=channel_id;
14+ this->parameters = NULL;
15+ this->element = NULL;
16+ this->nodes = NULL;
17
18- /*channel constructor data: */
19- int channel_elem_ids[2];
20- int channel_vertex_ids[2];
21- int channel_node_ids[4];
22- int channel_type;
23-
24- /*Get edge*/
25+ /*Get edge info*/
26 int i1 = iomodel->faces[4*index+0];
27 int i2 = iomodel->faces[4*index+1];
28 int e1 = iomodel->faces[4*index+2];
29 int e2 = iomodel->faces[4*index+3];
30
31- /*First, see wether this is an internal or boundary edge (if e2=-1)*/
32- if(e2==-1){
33- /* Boundary edge, only one element */
34- num_nodes=2;
35- channel_type=BoundaryEnum;
36- channel_elem_ids[0]=e1;
37- }
38- else{
39- /* internal edge: connected to 2 elements */
40- num_nodes=4;
41- channel_type=InternalEnum;
42- channel_elem_ids[0]=e1;
43- channel_elem_ids[1]=e2;
44- }
45+ /*Set Element hook (4th column may be -1 for boundary edges)*/
46+ this->helement = new Hook(&e1,1);
47
48- /*1: Get vertices ids*/
49+ /*Set Vertices hooks (4th column may be -1 for boundary edges)*/
50+ int channel_vertex_ids[2];
51 channel_vertex_ids[0]=i1;
52 channel_vertex_ids[1]=i2;
53+ this->hvertices =new Hook(&channel_vertex_ids[0],2);
54
55- /*2: Get node ids*/
56- if (channel_type==InternalEnum){
57+ /*Set Nodes hooks (!! Assumes P1 CG)*/
58+ int channel_node_ids[2];
59+ channel_node_ids[0]=i1;
60+ channel_node_ids[1]=i2;
61+ this->hnodes=new Hook(&channel_node_ids[0],2);
62
63- /*Now, we must get the nodes of the 4 nodes located on the edge*/
64-
65- /*2: Get the column where these ids are located in the index*/
66- pos1=pos2=pos3=pos4=UNDEF;
67- for(j=0;j<3;j++){
68- if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
69- if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
70- if(iomodel->elements[3*(e2-1)+j]==i1) pos3=j+1;
71- if(iomodel->elements[3*(e2-1)+j]==i2) pos4=j+1;
72- }
73- _assert_(pos1!=UNDEF && pos2!=UNDEF && pos3!=UNDEF && pos4!=UNDEF);
74-
75- /*3: We have the id of the elements and the position of the vertices in the index
76- * we can compute their dofs!*/
77- channel_node_ids[0]=3*(e1-1)+pos1;
78- channel_node_ids[1]=3*(e1-1)+pos2;
79- channel_node_ids[2]=3*(e2-1)+pos3;
80- channel_node_ids[3]=3*(e2-1)+pos4;
81- }
82- else{
83-
84- /*2: Get the column where these ids are located in the index*/
85- pos1=pos2=UNDEF;
86- for(j=0;j<3;j++){
87- if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
88- if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
89- }
90- _assert_(pos1!=UNDEF && pos2!=UNDEF);
91-
92- /*3: We have the id of the elements and the position of the vertices in the index
93- * we can compute their dofs!*/
94- channel_node_ids[0]=3*(e1-1)+pos1;
95- channel_node_ids[1]=3*(e1-1)+pos2;
96- }
97-
98- /*Ok, we have everything to build the object: */
99- this->id=channel_id;
100-
101- /*Hooks: */
102- this->hnodes =new Hook(channel_node_ids,num_nodes);
103- this->hvertices =new Hook(&channel_vertex_ids[0],2);
104- this->helement =new Hook(channel_elem_ids,1); // take only the first element for now
105-
106- //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
107- this->parameters=NULL;
108- this->element=NULL;
109- this->nodes=NULL;
110-}
111-/*}}}*/
112+}/*}}}*/
113 Channel::~Channel(){/*{{{*/
114 this->parameters=NULL;
115 delete helement;
116 delete hnodes;
117 delete hvertices;
118-}
119-/*}}}*/
120+}/*}}}*/
121
122 /*Object virtual functions definitions:*/
123 Object* Channel::copy() {/*{{{*/
124@@ -233,8 +171,13 @@
125 int analysis_type;
126 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
127
128- if(analysis_type != HydrologyGlaDSAnalysisEnum) _error_("Analysis not supported!!");
129- _error_("STOP");
130+ switch(analysis_type){
131+ case HydrologyGlaDSAnalysisEnum:
132+ Ke = this->CreateKMatrixHydrologyGlaDS();
133+ break;
134+ default:
135+ _error_("Don't know why we should be here");
136+ }
137
138 /*Add to global matrix*/
139 if(Ke){
140@@ -250,9 +193,15 @@
141 ElementVector* pe=NULL;
142 int analysis_type;
143 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
144- if(analysis_type != HydrologyGlaDSAnalysisEnum) _error_("Analysis not supported!!");
145- _error_("STOP");
146
147+ switch(analysis_type){
148+ case HydrologyGlaDSAnalysisEnum:
149+ pe = this->CreatePVectorHydrologyGlaDS();
150+ break;
151+ default:
152+ _error_("Don't know why we should be here");
153+ }
154+
155 /*Add to global matrix*/
156 if(pe){
157 pe->AddToGlobal(pf);
158@@ -371,3 +320,31 @@
159 *po_nz=o_nz;
160 }
161 /*}}}*/
162+
163+/*Channel specific functions*/
164+ElementVector* Channel::CreatePVectorHydrologyGlaDS(void){/*{{{*/
165+
166+ _error_("not implemented :( ");
167+
168+ /*Initialize Element matrix and return if necessary*/
169+ Tria* tria=(Tria*)element;
170+ if(!tria->IsIceInElement()) return NULL;
171+ ElementVector* pe=new ElementVector(nodes,NUMNODES,this->parameters);
172+
173+ /*Clean up and return*/
174+ return pe;
175+}
176+/*}}}*/
177+ElementMatrix* Channel::CreateKMatrixHydrologyGlaDS(void){/*{{{*/
178+
179+ _error_("not implemented :( ");
180+
181+ /*Initialize Element matrix and return if necessary*/
182+ Tria* tria=(Tria*)element;
183+ if(!tria->IsIceInElement()) return NULL;
184+ ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES,this->parameters);
185+
186+ /*Clean up and return*/
187+ return Ke;
188+}
189+/*}}}*/
190Index: ../trunk-jpl/src/c/classes/Loads/Channel.h
191===================================================================
192--- ../trunk-jpl/src/c/classes/Loads/Channel.h (revision 23959)
193+++ ../trunk-jpl/src/c/classes/Loads/Channel.h (revision 23960)
194@@ -71,6 +71,8 @@
195 /*}}}*/
196 /*Channel management:{{{*/
197 void GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
198+ ElementVector* CreatePVectorHydrologyGlaDS(void);
199+ ElementMatrix* CreateKMatrixHydrologyGlaDS(void);
200 /*}}}*/
201
202 };
203Index: ../trunk-jpl/src/c/classes/Loads/Moulin.cpp
204===================================================================
205--- ../trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23959)
206+++ ../trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23960)
207@@ -158,20 +158,19 @@
208 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
209
210 switch(analysis_type){
211-
212- case HydrologyShaktiAnalysisEnum:
213- pe = this->CreatePVectorHydrologyShakti();
214- break;
215- case HydrologyDCInefficientAnalysisEnum:
216- pe = CreatePVectorHydrologyDCInefficient();
217- break;
218- case HydrologyDCEfficientAnalysisEnum:
219- pe = CreatePVectorHydrologyDCEfficient();
220- break;
221- default:
222- _error_("Don't know why we should be here");
223- /*No loads applied, do nothing: */
224- return;
225+ case HydrologyShaktiAnalysisEnum:
226+ pe = this->CreatePVectorHydrologyShakti();
227+ break;
228+ case HydrologyDCInefficientAnalysisEnum:
229+ pe = CreatePVectorHydrologyDCInefficient();
230+ break;
231+ case HydrologyDCEfficientAnalysisEnum:
232+ pe = CreatePVectorHydrologyDCEfficient();
233+ break;
234+ default:
235+ _error_("Don't know why we should be here");
236+ /*No loads applied, do nothing: */
237+ return;
238 }
239 if(pe){
240 pe->AddToGlobal(pf);
241Index: ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
242===================================================================
243--- ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23959)
244+++ ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23960)
245@@ -25,7 +25,25 @@
246 /*Now, do we really want GlaDS?*/
247 if(hydrology_model!=HydrologyGlaDSEnum) return;
248
249- /*No extra loads*/
250+ /*Add channels?*/
251+ bool ischannels;
252+ iomodel->FindConstant(&ischannels,"md.hydrology.ischannels");
253+ if(ischannels){
254+ /*Get faces (edges in 2d)*/
255+ CreateFaces(iomodel);
256+ for(int i=0;i<iomodel->numberoffaces;i++){
257+ /*Get left and right elements*/
258+ int element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
259+
260+ /*Now, if this element is not in the partition*/
261+ if(!iomodel->my_elements[element]) continue;
262+
263+ /* Add load */
264+ loads->AddObject(new Channel(i+1,i,i,iomodel));
265+ }
266+
267+ /*Free data: */
268+ }
269 }/*}}}*/
270 void HydrologyGlaDSAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
271
Note: See TracBrowser for help on using the repository browser.