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

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

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