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

Last change on this file since 3417 was 3417, checked in by Eric.Larour, 15 years ago

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