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

Last change on this file since 10356 was 10356, checked in by seroussi, 14 years ago

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