source: issm/trunk/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationxUtils.cpp@ 10309

Last change on this file since 10309 was 10309, checked in by seroussi, 13 years ago

start to reorganize grounding lines
File size: 9.1 KB
Line 
1/*!\file GroundinglineMigrationx
2 * \brief: migration grounding line position.
3 */
4
5#include "./GroundinglineMigrationx.h"
6#include "./GroundinglineMigrationxLocal.h"
7#include "../../shared/shared.h"
8#include "../../io/io.h"
9#include "../../io/io.h"
10#include "../../include/include.h"
11#include "../../toolkits/toolkits.h"
12#include "../../EnumDefinitions/EnumDefinitions.h"
13
14/*FUNCTION PotentialSheetUngrounding {{{1*/
15double* PotentialSheetUngrounding(Elements* elements,Nodes* nodes,Parameters* parameters){
16
17 int i;
18 Element* element=NULL;
19 Vec vec_potential_sheet_ungrounding=NULL;
20 double* potential_sheet_ungrounding=NULL;
21 int numnods;
22 int configuration_type;
23
24 /*recover parameters: */
25 parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
26
27 /*First, initialize vec_new_shelf_nodes, which will track which nodes have changed status: */
28 numnods=nodes->NumberOfNodes(configuration_type);
29 vec_potential_sheet_ungrounding=NewVec(numnods);
30
31 /*Loop through elements, and fill vec_potential_sheet_ungrounding: */
32 for(i=0;i<elements->Size();i++){
33 element=(Element*)elements->GetObjectByOffset(i);
34 element->PotentialSheetUngrounding(vec_potential_sheet_ungrounding);
35 }
36
37 /*Assemble vector: */
38 VecAssemblyBegin(vec_potential_sheet_ungrounding);
39 VecAssemblyEnd(vec_potential_sheet_ungrounding);
40
41 /*Serialize vector: */
42 VecToMPISerial(&potential_sheet_ungrounding,vec_potential_sheet_ungrounding);
43
44 /*free ressouces: */
45 VecFree(&vec_potential_sheet_ungrounding);
46
47 return potential_sheet_ungrounding;
48}
49/*}}}*/
50/*FUNCTION PropagateShelfIntoConnexIceSheet {{{1*/
51double* PropagateShelfIntoConnexIceSheet(Elements* elements,Nodes* nodes,Parameters* parameters,double* potential_sheet_ungrounding){
52
53 int i;
54 Element* element=NULL;
55 int numnods;
56 int analysis_type;
57 Vec vec_nodes_on_iceshelf=NULL;
58 double* nodes_on_iceshelf=NULL;
59 int nflipped,local_nflipped;
60 double* elements_touching_iceshelf=NULL;
61
62 /*recover parameters: */
63 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
64
65 /*recover vec_nodes_on_iceshelf: */
66 vec_nodes_on_iceshelf=CreateNodesOnIceShelf(nodes,analysis_type);
67
68 nflipped=1; //bootstrap
69 while(nflipped){
70
71 /*get a list of potential elements that have nodes on ice shelf: */
72 elements_touching_iceshelf=CreateElementTouchingIceShelf(elements,vec_nodes_on_iceshelf);
73
74 /*now, go through elements_touching_iceshelf, and if they have nodes inside potential_sheet_ungrounding,
75 * flag it: */
76 local_nflipped=0;
77
78 /*serialize vec_nodes_on_iceshelf, needed by element's UpdatePotentialSheetUngrounding routine, to figure out if
79 * nodes have flipped from grounded to ungrounded: */
80 VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);
81
82 for(i=0;i<elements->Size();i++){
83 element=(Element*)elements->GetObjectByOffset(i);
84 if(elements_touching_iceshelf[element->Sid()]){
85 local_nflipped+=element->UpdatePotentialSheetUngrounding(potential_sheet_ungrounding,vec_nodes_on_iceshelf,nodes_on_iceshelf);
86 }
87 }
88
89 MPI_Allreduce(&local_nflipped,&nflipped,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
90 _printf_(VerboseConvergence()," number of grounded vertices connected to grounding line: %i\n",nflipped);
91
92 /*Avoid leaks: */
93 xfree((void**)&nodes_on_iceshelf);
94 xfree((void**)&elements_touching_iceshelf);
95
96 /*Assemble:*/
97 VecAssemblyBegin(vec_nodes_on_iceshelf);
98 VecAssemblyEnd(vec_nodes_on_iceshelf);
99 }
100
101 /*Serialize: */
102 VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);
103
104 /*Free ressources:*/
105 VecFree(&vec_nodes_on_iceshelf);
106 xfree((void**)&elements_touching_iceshelf);
107
108 return nodes_on_iceshelf;
109}
110/*}}}*/
111/*FUNCTION CreateElementOnGroundingline {{{1*/
112bool* CreateElementOnGroundingline(Elements* elements,double* element_on_iceshelf){
113
114 int i,j;
115 Element *element = NULL;
116 bool *element_on_gl = NULL;
117 bool ongl=false;
118 int *neighboorsids = NULL;
119 int sid;
120 int shelf;
121
122 /*Go through elements, and look for elements that can possibly have a grounding line migration. These
123 * are elements that touch the grounding line: */
124 element_on_gl=(bool*)xmalloc(elements->Size()*sizeof(bool));
125
126 for(i=0;i<elements->Size();i++){
127 element=(Element*)elements->GetObjectByOffset(i);
128
129 sid=element->Sid();
130 neighboorsids=element->GetHorizontalNeighboorSids();
131
132 shelf=(int)element_on_iceshelf[sid];
133 ongl=false;
134 for(j=0;j<3;j++){
135 if (neighboorsids[j]<0)continue; //this neighboor does not exist
136 if ((shelf==1) & (element_on_iceshelf[neighboorsids[j]]==1))continue; //both this element and this neighboor are on th ice shelf
137 if ((shelf==0) & (element_on_iceshelf[neighboorsids[j]]==0))continue; //both this element and this neighboor are on the ice sheet
138 ongl=true; //neighboor j is on a different location than us, ie we are touching the grounding line.
139 }
140 element_on_gl[i]=ongl;
141 }
142
143 return element_on_gl;
144}
145/*}}}*/
146/*FUNCTION CreateElementOnIceShelf {{{1*/
147double* CreateElementOnIceShelf(Elements* elements){
148
149 int i;
150 Element* element=NULL;
151 Vec vec_element_on_iceshelf=NULL;
152 double* element_on_iceshelf=NULL;
153
154 /*Create vector holding all the elements IsFloating flags: */
155 vec_element_on_iceshelf=NewVec(elements->NumberOfElements(),true);
156
157 /*Loop through elements, and fill vec_element_on_iceshelf: */
158 for(i=0;i<elements->Size();i++){
159 element=(Element*)elements->GetObjectByOffset(i);
160 VecSetValue(vec_element_on_iceshelf,element->Sid(),(int)element->IsFloating(),INSERT_VALUES);
161 }
162
163 /*Assemble vector: */
164 VecAssemblyBegin(vec_element_on_iceshelf);
165 VecAssemblyEnd(vec_element_on_iceshelf);
166
167 /*Serialize vector: */
168 VecToMPISerial(&element_on_iceshelf,vec_element_on_iceshelf);
169
170 /*free ressouces: */
171 VecFree(&vec_element_on_iceshelf);
172
173 return element_on_iceshelf;
174}
175/*}}}*/
176/*FUNCTION CreateElementTouchingIceShelf {{{1*/
177double* CreateElementTouchingIceShelf(Elements* elements,Vec vec_nodes_on_iceshelf){
178
179 int i;
180 Element* element=NULL;
181 Vec vec_element_touching_iceshelf=NULL;
182 double* nodes_on_iceshelf=NULL;
183
184 /*output: */
185 double* element_touching_iceshelf=NULL;
186
187 /*serialize: */
188 VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);
189
190 /*Create vector holding all the elements IsFloating flags: */
191 vec_element_touching_iceshelf=NewVec(elements->NumberOfElements(),true);
192
193 /*Loop through elements, and fill vec_element_touching_iceshelf: */
194 for(i=0;i<elements->Size();i++){
195 element=(Element*)elements->GetObjectByOffset(i);
196 VecSetValue(vec_element_touching_iceshelf,element->Sid(),element->IsNodeOnShelfFromFlags(nodes_on_iceshelf)?1.0:0.0,INSERT_VALUES);
197 }
198
199 /*Assemble vector: */
200 VecAssemblyBegin(vec_element_touching_iceshelf);
201 VecAssemblyEnd(vec_element_touching_iceshelf);
202
203 /*Serialize vector: */
204 VecToMPISerial(&element_touching_iceshelf,vec_element_touching_iceshelf);
205
206 /*free ressouces: */
207 VecFree(&vec_element_touching_iceshelf);
208 xfree((void**)&nodes_on_iceshelf);
209
210 return element_touching_iceshelf;
211}
212/*}}}*/
213/*FUNCTION UpdateShelfStatus {{{1*/
214int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf){
215
216 int i;
217 Element* element=NULL;
218 int configuration_type;
219 int numnods;
220 Vec vec_new_shelf_nodes=NULL;
221 double* new_shelf_nodes=NULL;
222
223 /*output: */
224 int local_nflipped=0;
225 int nflipped=0;
226
227 /*recover parameters: */
228 parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
229
230 /*First, initialize vec_new_shelf_nodes, which will track which nodes have changed status: */
231 numnods=nodes->NumberOfNodes(configuration_type);
232 vec_new_shelf_nodes=NewVec(numnods);
233
234 /*Ok, now go through the elements that have gone through grounding line migration,
235 * and update their flags: */
236 for(i=0;i<elements->Size();i++){
237 element=(Element*)elements->GetObjectByOffset(i);
238 if(element_touching_iceshelf[element->Sid()]){
239 local_nflipped+=element->UpdateShelfStatus(vec_new_shelf_nodes);
240 }
241 }
242 MPI_Allreduce(&local_nflipped,&nflipped,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
243
244 /*Serialize vec_new_shelf_nodes: */
245 VecToMPISerial(&new_shelf_nodes,vec_new_shelf_nodes);
246
247 /*Now, go through ALL elements, and update the status of the nodes, to propagate what happened at the grounding line: */
248 /*Carry out grounding line migration for those elements: */
249 for(i=0;i<elements->Size();i++){
250 element=(Element*)elements->GetObjectByOffset(i);
251 element->UpdateShelfFlags(new_shelf_nodes);
252 }
253
254 /*Free ressources: */
255 VecFree(&vec_new_shelf_nodes);
256 xfree((void**)&new_shelf_nodes);
257
258 return nflipped;
259}
260/*}}}*/
261/*FUNCTION CreateNodesOnIceShelf {{{1*/
262Vec CreateNodesOnIceShelf(Nodes* nodes,int configuration_type){
263
264 /*output: */
265 Vec nodes_on_iceshelf=NULL;
266
267 /*intermediary: */
268 int numnods;
269 int i;
270 Node* node=NULL;
271
272 /*First, initialize nodes_on_iceshelf, which will track which nodes have changed status: */
273 numnods=nodes->NumberOfNodes(configuration_type);
274 nodes_on_iceshelf=NewVec(numnods);
275
276 /*Loop through nodes, and fill nodes_on_iceshelf: */
277 for(i=0;i<nodes->Size();i++){
278 node=(Node*)nodes->GetObjectByOffset(i);
279 if(node->InAnalysis(configuration_type)){
280 if(node->IsFloating()){
281 VecSetValue(nodes_on_iceshelf,node->Sid(),1.0,INSERT_VALUES);
282 }
283 }
284 }
285
286 /*Assemble vector: */
287 VecAssemblyBegin(nodes_on_iceshelf);
288 VecAssemblyEnd(nodes_on_iceshelf);
289
290 return nodes_on_iceshelf;
291}
292/*}}}*/
Note: See TracBrowser for help on using the repository browser.