source: issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp@ 3378

Last change on this file since 3378 was 3378, checked in by Mathieu Morlighem, 16 years ago

fixed DG in parallel
File size: 14.3 KB
Line 
1/*
2 * CreateElementsNodesAndMaterialsPrognostic2.c:
3 */
4
5#include "../../DataSet/DataSet.h"
6#include "../../toolkits/toolkits.h"
7#include "../../EnumDefinitions/EnumDefinitions.h"
8#include "../../objects/objects.h"
9#include "../../shared/shared.h"
10#include "../../include/macros.h"
11#include "../../include/typedefs.h"
12#include "../../MeshPartitionx/MeshPartitionx.h"
13#include "../IoModel.h"
14
15void CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
16
17 int i,j,k,n;
18 extern int my_rank;
19 extern int num_procs;
20
21 /*DataSets: */
22 DataSet* elements = NULL;
23 DataSet* nodes = NULL;
24 DataSet* materials = NULL;
25
26 /*Objects: */
27 Node* node = NULL;
28 Tria* tria = NULL;
29 Penta* penta = NULL;
30 Matice* matice = NULL;
31 Matpar* matpar = NULL;
32
33 /*output: */
34 int* epart=NULL; //element partitioning.
35 int* npart=NULL; //node partitioning.
36 int* my_grids=NULL;
37 double* my_bordergrids=NULL;
38
39 /*intermediary: */
40 int elements_width; //size of elements
41 double B_avg;
42
43 /*tria constructor input: */
44 int tria_id;
45 int tria_mid;
46 int tria_mparid;
47 int tria_numparid;
48 int tria_g[3];
49 double tria_h[3];
50 double tria_s[3];
51 double tria_b[3];
52 double tria_k[3];
53 double tria_melting[3];
54 double tria_accumulation[3];
55 double tria_geothermalflux[3];
56 int tria_friction_type;
57 double tria_p;
58 double tria_q;
59 int tria_shelf;
60 bool tria_onwater;
61
62 /*matice constructor input: */
63 int matice_mid;
64 double matice_B;
65 double matice_n;
66
67 /*penta constructor input: */
68
69 int penta_id;
70 int penta_mid;
71 int penta_mparid;
72 int penta_numparid;
73 int penta_g[6];
74 double penta_h[6];
75 double penta_s[6];
76 double penta_b[6];
77 double penta_k[6];
78 int penta_friction_type;
79 double penta_p;
80 double penta_q;
81 int penta_shelf;
82 int penta_onbed;
83 int penta_onsurface;
84 int penta_collapse;
85 double penta_melting[6];
86 double penta_accumulation[6];
87 double penta_geothermalflux[6];
88 int penta_thermal_steadystate;
89 bool penta_onwater;
90
91 /*matpar constructor input: */
92 int matpar_mid;
93 double matpar_rho_ice;
94 double matpar_rho_water;
95 double matpar_heatcapacity;
96 double matpar_thermalconductivity;
97 double matpar_latentheat;
98 double matpar_beta;
99 double matpar_meltingpoint;
100 double matpar_mixed_layer_capacity;
101 double matpar_thermal_exchange_velocity;
102 double matpar_g;
103
104 /* node constructor input: */
105 int node_id;
106 int node_partitionborder=0;
107 double node_x[3];
108 double node_sigma;
109 int node_onbed;
110 int node_onsurface;
111 int node_onshelf;
112 int node_onsheet;
113 int node_upper_node_id;
114 int node_numdofs;
115
116
117 #ifdef _PARALLEL_
118 /*Metis partitioning: */
119 int range;
120 Vec gridborder=NULL;
121 int my_numgrids;
122 int* all_numgrids=NULL;
123 int gridcount;
124 int count;
125 /*Nodes cloning*/
126 double e1,e2;
127 int i1,i2;
128 int pos;
129 #endif
130 int first_grid_index;
131
132 /*First create the elements, nodes and material properties: */
133 elements = new DataSet(ElementsEnum());
134 nodes = new DataSet(NodesEnum());
135 materials = new DataSet(MaterialsEnum());
136
137 /*Width of elements: */
138 if(strcmp(iomodel->meshtype,"2d")==0){
139 elements_width=3; //tria elements
140 }
141 else{
142 ISSMERROR("not implemented yet!");
143 elements_width=6; //penta elements
144 }
145
146 #ifdef _PARALLEL_
147 /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
148 if(strcmp(iomodel->meshtype,"2d")==0){
149 /*load elements: */
150 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
151 }
152 else{
153 /*load elements2d: */
154 IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d");
155 }
156
157
158 MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs);
159
160 /*Free elements and elements2d: */
161 xfree((void**)&iomodel->elements);
162 xfree((void**)&iomodel->elements2d);
163
164 /*Used later on: */
165 my_grids=(int*)xcalloc(3*iomodel->numberofelements,sizeof(int));
166 #endif
167
168 /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */
169
170 /*2d mesh: */
171 if (strcmp(iomodel->meshtype,"2d")==0){
172
173 /*Fetch data needed: */
174 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
175 IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
176 IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface");
177 IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
178 IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf");
179 IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
180
181 for (i=0;i<iomodel->numberofelements;i++){
182
183 #ifdef _PARALLEL_
184 /*!All elements have been partitioned above, only create elements for this CPU: */
185 if(my_rank==epart[i]){
186 #endif
187
188 /*ids: */
189 tria_id=i+1; //matlab indexing.
190 tria_mid=-1; //no need for materials
191 tria_mparid=-1; //no need for materials
192 tria_numparid=1;
193
194 /*vertices offsets: */
195 tria_g[0]=3*i+1;
196 tria_g[1]=3*i+2;
197 tria_g[2]=3*i+3;
198
199 /*thickness,surface and bed:*/
200 tria_h[0]= *(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
201 tria_h[1]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+1)-1));
202 tria_h[2]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+2)-1)) ;
203
204 tria_s[0]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+0)-1));
205 tria_s[1]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+1)-1));
206 tria_s[2]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+2)-1));
207
208 tria_b[0]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+0)-1));
209 tria_b[1]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+1)-1));
210 tria_b[2]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+2)-1));
211
212 /*element on iceshelf?:*/
213 tria_shelf=(int)*(iomodel->elementoniceshelf+i);
214 tria_onwater=(bool)*(iomodel->elementonwater+i);
215
216 /*Create tria element using its constructor:*/
217 tria=new Tria(tria_id, tria_mid, tria_mparid, tria_numparid,tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_onwater);
218
219 /*Add tria element to elements dataset: */
220 elements->AddObject(tria);
221
222 #ifdef _PARALLEL_
223 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
224 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
225 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
226 will hold which grids belong to this partition*/
227 my_grids[3*i+0]=1;
228 my_grids[3*i+1]=1;
229 my_grids[3*i+2]=1;
230 #endif
231
232 #ifdef _PARALLEL_
233 }//if(my_rank==epart[i])
234 #endif
235
236 }//for (i=0;i<numberofelements;i++)
237
238 /*Free data : */
239 xfree((void**)&iomodel->thickness);
240 xfree((void**)&iomodel->surface);
241 xfree((void**)&iomodel->bed);
242 xfree((void**)&iomodel->elementoniceshelf);
243 xfree((void**)&iomodel->elementonwater);
244
245 }
246 else{ // if (strcmp(meshtype,"2d")==0)
247 ISSMERROR(exprintf("not implemented yet"));
248 } //if (strcmp(meshtype,"2d")==0)
249
250 #ifdef _PARALLEL_
251 /*If we are in parallel, we must add the nodes that are not in this partition
252 * but share edges of elements that are in this partition. This is needed for
253 * the loads as numerical fluxes involve the dofs of all nodes on a given edge*/
254
255 /*Get edges and elements*/
256 IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges");
257
258 /*!All elements have been partitioned above, only create elements for this CPU: */
259 for (i=0;i<iomodel->numberofedges;i++){
260
261 /*Get left and right elements*/
262 e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2]
263 e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2]
264
265 /* 1) If the element e1 is in the current partition
266 * 2) and if the edge of the element is shared by another element
267 * 3) and if this element is not in the same partition:
268 * we must clone the nodes on this partition so that the loads (Numericalflux)
269 * will have access to their properties (dofs,...)*/
270 if(my_rank==epart[(int)e1] && !isnan(e2) && my_rank!=epart[(int)e2]){
271
272 /*1: Get vertices ids*/
273 i1=(int)iomodel->edges[4*i+0];
274 i2=(int)iomodel->edges[4*i+1];
275
276 /*2: Get the column where these ids are located in the index*/
277 pos==UNDEF;
278 for(j=0;j<3;j++){
279 if (iomodel->elements[3*(int)e2+j]==i1) pos=j+1;
280 }
281 ISSMASSERT(pos!=UNDEF);
282
283 /*3: We have the id of the elements and the position of the vertices in the index
284 * we can now create the corresponding nodes:*/
285 my_grids[(int)e2*3+pos-1]=1;
286 my_grids[(int)e2*3+((pos+1)%3)]=1;
287 }
288 }
289
290 /*Free data: */
291 xfree((void**)&iomodel->edges);
292
293 /*From the element partitioning, we can determine which nodes are on the inside of this cpu's
294 *element partition, and which are on its border with other nodes:*/
295 gridborder=NewVec(3*iomodel->numberofelements);
296
297 for (i=0;i<3*iomodel->numberofelements;i++){
298 if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
299 }
300 VecAssemblyBegin(gridborder);
301 VecAssemblyEnd(gridborder);
302
303 VecToMPISerial(&my_bordergrids,gridborder);
304
305 #endif
306
307 /*Add one constant material property to materials: */
308 matpar_mid=1; //put it at the end of the materials
309 matpar_g=iomodel->g;
310 matpar_rho_ice=iomodel->rho_ice;
311 matpar_rho_water=iomodel->rho_water;
312 matpar_thermalconductivity=iomodel->thermalconductivity;
313 matpar_heatcapacity=iomodel->heatcapacity;
314 matpar_latentheat=iomodel->latentheat;
315 matpar_beta=iomodel->beta;
316 matpar_meltingpoint=iomodel->meltingpoint;
317 matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity;
318 matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity;
319
320 /*Create matpar object using its constructor: */
321 matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,
322 matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,
323 matpar_thermal_exchange_velocity,matpar_g);
324
325 /*Add to materials datset: */
326 materials->AddObject(matpar);
327
328 /*Ok, let's summarise. Now, every CPU has the following array: my_grids
329 We can therefore determine which grids are internal to this node's partition
330 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing
331 that, go and create the grids*/
332
333 /*Create nodes from x,y,z, as well as the spc values on those grids: */
334
335 /*First fetch data: */
336 if (strcmp(iomodel->meshtype,"3d")==0){
337 IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids");
338 IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
339 }
340 IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x");
341 IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y");
342 IoModelFetchData(&iomodel->z,NULL,NULL,iomodel_handle,"z");
343 IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
344 IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
345 IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed");
346 IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface");
347 IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
348 IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
349
350 /*Get number of dofs per node: */
351 DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);
352
353 /*Build Nodes dataset -> 3 for each element*/
354 for (i=0;i<iomodel->numberofelements;i++){
355 for (j=0;j<3;j++){
356
357 #ifdef _PARALLEL_
358 /*!All elements have been partitioned above, only create elements for this CPU: */
359 if(my_grids[3*i+j]){
360 #endif
361
362 //Get id of the vertex on which the current node is located
363 k=(int)*(iomodel->elements+elements_width*i+j)-1; //(Matlab to C indexing)
364 ISSMASSERT(k>=0 && k<iomodel->numberofnodes);
365
366 //Get node properties
367 node_id=i*3+j+1;
368 #ifdef _PARALLEL_
369 if(my_bordergrids[node_id-1]>1.0) { //this grid belongs to a partition border
370 node_partitionborder=1;
371 }
372 else{
373 node_partitionborder=0;
374 }
375 #else
376 node_partitionborder=0;
377 #endif
378 node_x[0]=iomodel->x[k];
379 node_x[1]=iomodel->y[k];
380 node_x[2]=iomodel->z[k];
381 node_sigma=(iomodel->z[k]-iomodel->bed[k])/(iomodel->thickness[k]);
382 node_onbed=(int)iomodel->gridonbed[k];
383 node_onsurface=(int)iomodel->gridonsurface[k];
384 node_onshelf=(int)iomodel->gridoniceshelf[k];
385 node_onsheet=(int)iomodel->gridonicesheet[k];
386 if (strcmp(iomodel->meshtype,"3d")==0){
387 if (isnan(iomodel->uppernodes[k])){
388 node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves.
389 }
390 else{
391 node_upper_node_id=(int)iomodel->uppernodes[k];
392 }
393 }
394 else{
395 /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
396 node_upper_node_id=node_id;
397 }
398
399 /*Create node using its constructor: */
400 node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
401
402 /*Add node to nodes dataset: */
403 nodes->AddObject(node);
404 #ifdef _PARALLEL_
405 }
406 #endif
407 }
408 }
409
410 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
411 * datasets, it will not be redone: */
412 elements->Presort();
413 nodes->Presort();
414 materials->Presort();
415
416 /*Clean fetched data: */
417 xfree((void**)&iomodel->deadgrids);
418 xfree((void**)&iomodel->elements);
419 xfree((void**)&iomodel->x);
420 xfree((void**)&iomodel->y);
421 xfree((void**)&iomodel->z);
422 xfree((void**)&iomodel->thickness);
423 xfree((void**)&iomodel->bed);
424 xfree((void**)&iomodel->gridonbed);
425 xfree((void**)&iomodel->gridonsurface);
426 xfree((void**)&iomodel->uppernodes);
427 xfree((void**)&iomodel->gridonicesheet);
428 xfree((void**)&iomodel->gridoniceshelf);
429
430 /*Keep partitioning information into iomodel*/
431 iomodel->epart=epart;
432 iomodel->my_grids=my_grids;
433 iomodel->my_bordergrids=my_bordergrids;
434
435 /*Free ressources:*/
436 #ifdef _PARALLEL_
437 xfree((void**)&all_numgrids);
438 xfree((void**)&npart);
439 VecFree(&gridborder);
440 #endif
441
442 cleanup_and_return:
443
444 /*Assign output pointer: */
445 *pelements=elements;
446 *pnodes=nodes;
447 *pmaterials=materials;
448}
Note: See TracBrowser for help on using the repository browser.