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

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

some cleaning
File size: 8.9 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,double* nodes_on_iceshelf){
78
79 int i;
80 double *element_touching_iceshelf = NULL;
81 Element *element = NULL;
82 Vec vec_element_touching_iceshelf = NULL;
83
84 /*Create vector holding all the elements IsFloating flags: */
85 vec_element_touching_iceshelf=NewVec(elements->NumberOfElements(),true);
86
87 /*Loop through elements, and fill vec_element_touching_iceshelf: */
88 for(i=0;i<elements->Size();i++){
89 element=(Element*)elements->GetObjectByOffset(i);
90 VecSetValue(vec_element_touching_iceshelf,element->Sid(),element->IsNodeOnShelfFromFlags(nodes_on_iceshelf)?1.0:0.0,INSERT_VALUES);
91 }
92
93 /*Assemble vector: */
94 VecAssemblyBegin(vec_element_touching_iceshelf);
95 VecAssemblyEnd(vec_element_touching_iceshelf);
96
97 /*Serialize vector: */
98 VecToMPISerial(&element_touching_iceshelf,vec_element_touching_iceshelf);
99
100 /*free ressouces: */
101 VecFree(&vec_element_touching_iceshelf);
102
103 return element_touching_iceshelf;
104}
105/*}}}*/
106/*FUNCTION CreateNodesOnIceShelf {{{1*/
107Vec CreateNodesOnIceShelf(Nodes* nodes,Vertices* vertices,int configuration_type){
108
109 int i,numberofvertices;
110 Vec nodes_on_iceshelf = NULL;
111 Node* node = NULL;
112
113 /*First, initialize nodes_on_iceshelf, which will track which nodes have changed status: */
114 numberofvertices=vertices->NumberOfVertices();
115 nodes_on_iceshelf=NewVec(numberofvertices);
116
117 /*Loop through nodes, and fill nodes_on_iceshelf: */
118 for(i=0;i<nodes->Size();i++){
119 node=(Node*)nodes->GetObjectByOffset(i);
120 if(node->InAnalysis(configuration_type)){
121 if(node->IsFloating()){
122 VecSetValue(nodes_on_iceshelf,node->Sid(),1.0,INSERT_VALUES);
123 }
124 }
125 }
126
127 /*Assemble vector: */
128 VecAssemblyBegin(nodes_on_iceshelf);
129 VecAssemblyEnd(nodes_on_iceshelf);
130
131 return nodes_on_iceshelf;
132}
133/*}}}*/
134/*FUNCTION PotentialSheetUngrounding {{{1*/
135double* PotentialSheetUngrounding(Elements* elements,Vertices* vertices,Parameters* parameters){
136
137 int i,numberofvertices;
138 double* vertices_potentially_ungrounding = NULL;
139 Vec vec_vertices_potentially_ungrounding = NULL;
140 Element* element = NULL;
141
142 /*Initialize vector with number of vertices*/
143 numberofvertices=vertices->NumberOfVertices();
144 vec_vertices_potentially_ungrounding=NewVec(numberofvertices); //grounded vertex that could start floating
145
146 /*Fill vector vertices_potentially_floating: */
147 for(i=0;i<elements->Size();i++){
148 element=(Element*)elements->GetObjectByOffset(i);
149 element->PotentialSheetUngrounding(vec_vertices_potentially_ungrounding);
150 }
151
152 /*Assemble vector: */
153 VecAssemblyBegin(vec_vertices_potentially_ungrounding);
154 VecAssemblyEnd(vec_vertices_potentially_ungrounding);
155
156 /*Serialize vector: */
157 VecToMPISerial(&vertices_potentially_ungrounding,vec_vertices_potentially_ungrounding);
158
159 /*free ressouces and return: */
160 VecFree(&vec_vertices_potentially_ungrounding);
161 return vertices_potentially_ungrounding;
162}
163/*}}}*/
164/*FUNCTION PropagateShelfIntoConnexIceSheet {{{1*/
165double* PropagateShelfIntoConnexIceSheet(Elements* elements,Nodes* nodes,Vertices* vertices,Parameters* parameters,double* vertices_potentially_ungrounding){
166
167 int i,analysis_type;
168 int nflipped,local_nflipped;
169 double* nodes_on_iceshelf = NULL;
170 double* elements_touching_iceshelf = NULL;
171 Vec vec_nodes_on_iceshelf = NULL;
172 Element* element = NULL;
173
174 /*recover parameters: */
175 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
176
177 /*recover vec_nodes_on_iceshelf and serialize: */
178 vec_nodes_on_iceshelf=CreateNodesOnIceShelf(nodes,vertices,analysis_type);
179 VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);
180
181 nflipped=1; //bootstrap
182 while(nflipped){
183
184 /*get a list of potential elements that have nodes on ice shelf: */
185 elements_touching_iceshelf=CreateElementTouchingIceShelf(elements,nodes_on_iceshelf);
186
187 /*Go through elements_touching_iceshelf, and flag their nodes that can unground inside potential_sheet_ungrounding*/
188 local_nflipped=0;
189 for(i=0;i<elements->Size();i++){
190 element=(Element*)elements->GetObjectByOffset(i);
191 if(elements_touching_iceshelf[element->Sid()]){
192 local_nflipped+=element->UpdatePotentialSheetUngrounding(vertices_potentially_ungrounding,vec_nodes_on_iceshelf,nodes_on_iceshelf);
193 }
194 }
195
196 MPI_Allreduce(&local_nflipped,&nflipped,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
197 _printf_(VerboseConvergence()," number of grounded vertices connected to grounding line: %i\n",nflipped);
198
199 /*Avoid leaks: */
200 xfree((void**)&elements_touching_iceshelf);
201 xfree((void**)&nodes_on_iceshelf);
202
203 /*Assemble and serialize:*/
204 VecAssemblyBegin(vec_nodes_on_iceshelf);
205 VecAssemblyEnd(vec_nodes_on_iceshelf);
206 VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);
207 }
208
209 /*Free ressources:*/
210 VecFree(&vec_nodes_on_iceshelf);
211 xfree((void**)&elements_touching_iceshelf);
212
213 return nodes_on_iceshelf;
214}
215/*}}}*/
216
217/*FUNCTION UpdateShelfStatus {{{1*/
218int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf){
219
220 int i,numnods,configuration_type;
221 int nflipped=0;
222 int local_nflipped=0;
223 Vec vec_new_shelf_nodes = NULL;
224 double* new_shelf_nodes = NULL;
225 Element* element = NULL;
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/*}}}*/
Note: See TracBrowser for help on using the repository browser.