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

Last change on this file since 3360 was 3359, checked in by Mathieu Morlighem, 15 years ago

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