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

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

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