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 |
|
---|
16 | void 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]=3*i+1;
|
---|
194 | tria_g[1]=3*i+2;
|
---|
195 | tria_g[2]=3*i+3;
|
---|
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 | ISSMASSERT(k>=0 && k<iomodel->numberofnodes);
|
---|
336 |
|
---|
337 | //Get node properties
|
---|
338 | node_id=i*3+j+1;
|
---|
339 | node_partitionborder=0;
|
---|
340 | node_x[0]=iomodel->x[k];
|
---|
341 | node_x[1]=iomodel->y[k];
|
---|
342 | node_x[2]=iomodel->z[k];
|
---|
343 | node_sigma=(iomodel->z[k]-iomodel->bed[k])/(iomodel->thickness[k]);
|
---|
344 | node_onbed=(int)iomodel->gridonbed[k];
|
---|
345 | node_onsurface=(int)iomodel->gridonsurface[k];
|
---|
346 | node_onshelf=(int)iomodel->gridoniceshelf[k];
|
---|
347 | node_onsheet=(int)iomodel->gridonicesheet[k];
|
---|
348 | if (strcmp(iomodel->meshtype,"3d")==0){
|
---|
349 | if (isnan(iomodel->uppernodes[k])){
|
---|
350 | node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves.
|
---|
351 | }
|
---|
352 | else{
|
---|
353 | node_upper_node_id=(int)iomodel->uppernodes[k];
|
---|
354 | }
|
---|
355 | }
|
---|
356 | else{
|
---|
357 | /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
|
---|
358 | node_upper_node_id=node_id;
|
---|
359 | }
|
---|
360 |
|
---|
361 | /*Create node using its constructor: */
|
---|
362 | 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);
|
---|
363 |
|
---|
364 | /*Add node to nodes dataset: */
|
---|
365 | nodes->AddObject(node);
|
---|
366 |
|
---|
367 | #ifdef _PARALLEL_
|
---|
368 | }
|
---|
369 | #endif
|
---|
370 | }
|
---|
371 | }
|
---|
372 |
|
---|
373 | /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
|
---|
374 | * datasets, it will not be redone: */
|
---|
375 | elements->Presort();
|
---|
376 | nodes->Presort();
|
---|
377 | materials->Presort();
|
---|
378 |
|
---|
379 | /*Clean fetched data: */
|
---|
380 | xfree((void**)&iomodel->deadgrids);
|
---|
381 | xfree((void**)&iomodel->elements);
|
---|
382 | xfree((void**)&iomodel->x);
|
---|
383 | xfree((void**)&iomodel->y);
|
---|
384 | xfree((void**)&iomodel->z);
|
---|
385 | xfree((void**)&iomodel->thickness);
|
---|
386 | xfree((void**)&iomodel->bed);
|
---|
387 | xfree((void**)&iomodel->gridonbed);
|
---|
388 | xfree((void**)&iomodel->gridonsurface);
|
---|
389 | xfree((void**)&iomodel->uppernodes);
|
---|
390 | xfree((void**)&iomodel->gridonicesheet);
|
---|
391 | xfree((void**)&iomodel->gridoniceshelf);
|
---|
392 |
|
---|
393 | /*Keep partitioning information into iomodel*/
|
---|
394 | iomodel->epart=epart;
|
---|
395 | iomodel->my_grids=my_grids;
|
---|
396 | iomodel->my_bordergrids=my_bordergrids;
|
---|
397 |
|
---|
398 | /*Free ressources:*/
|
---|
399 | #ifdef _PARALLEL_
|
---|
400 | xfree((void**)&all_numgrids);
|
---|
401 | xfree((void**)&npart);
|
---|
402 | VecFree(&gridborder);
|
---|
403 | #endif
|
---|
404 |
|
---|
405 | cleanup_and_return:
|
---|
406 |
|
---|
407 | /*Assign output pointer: */
|
---|
408 | *pelements=elements;
|
---|
409 | *pnodes=nodes;
|
---|
410 | *pmaterials=materials;
|
---|
411 | }
|
---|