source: issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp@ 72

Last change on this file since 72 was 72, checked in by seroussi, 16 years ago

forgot file change for overshooting
File size: 20.7 KB
Line 
1/*
2 * CreateElementsNodesAndMaterialsDiagnosticHoriz.c:
3 */
4
5#undef __FUNCT__
6#define __FUNCT__ "CreateElementsNodesAndMaterialsDiagnosticHoriz"
7
8#include "../DataSet/DataSet.h"
9#include "../toolkits/toolkits.h"
10#include "../EnumDefinitions/EnumDefinitions.h"
11#include "../objects/objects.h"
12#include "../shared/shared.h"
13#include "./ModelProcessorx.h"
14#include "./MeshPartitionx/MeshPartitionx.h"
15
16
17int CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){
18
19
20 /*output: int* epart, int* my_grids, double* my_bordergrids*/
21
22
23 int noerr=1;
24 int i,j,k,n;
25 extern int my_rank;
26 extern int num_procs;
27
28 /*DataSets: */
29 DataSet* elements = NULL;
30 DataSet* nodes = NULL;
31 DataSet* materials = NULL;
32
33 /*Objects: */
34 Node* node = NULL;
35 Tria* tria = NULL;
36 Penta* penta = NULL;
37 Matice* matice = NULL;
38 Matpar* matpar = NULL;
39
40 int analysis_type;
41
42 /*output: */
43 int* epart=NULL; //element partitioning.
44 int* npart=NULL; //node partitioning.
45 int* my_grids=NULL;
46 double* my_bordergrids=NULL;
47
48
49 /*intermediary: */
50 int elements_width; //size of elements
51 double B_avg;
52
53 /*tria constructor input: */
54 int tria_id;
55 int tria_mid;
56 int tria_mparid;
57 int tria_g[3];
58 double tria_h[3];
59 double tria_s[3];
60 double tria_b[3];
61 double tria_k[3];
62 int tria_friction_type;
63 double tria_p;
64 double tria_q;
65 int tria_shelf;
66 double tria_meanvel;/*!scaling ratio for velocities*/
67 double tria_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
68 double tria_viscosity_overshoot;
69
70 /*matice constructor input: */
71 int matice_mid;
72 double matice_B;
73 double matice_n;
74
75 /*penta constructor input: */
76
77 int penta_id;
78 int penta_mid;
79 int penta_mparid;
80 int penta_g[6];
81 double penta_h[6];
82 double penta_s[6];
83 double penta_b[6];
84 double penta_k[6];
85 int penta_friction_type;
86 double penta_p;
87 double penta_q;
88 int penta_shelf;
89 int penta_onbed;
90 int penta_onsurface;
91 double penta_meanvel;/*!scaling ratio for velocities*/
92 double penta_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
93 int penta_collapse;
94 double penta_melting[6];
95 double penta_accumulation[6];
96 double penta_geothermalflux[6];
97 int penta_artdiff;
98 int penta_thermal_steadystate;
99
100 /*matpar constructor input: */
101 int matpar_mid;
102 double matpar_rho_ice;
103 double matpar_rho_water;
104 double matpar_heatcapacity;
105 double matpar_thermalconductivity;
106 double matpar_latentheat;
107 double matpar_beta;
108 double matpar_meltingpoint;
109 double matpar_mixed_layer_capacity;
110 double matpar_thermal_exchange_velocity;
111 double matpar_g;
112 double matpar_viscosity_overshoot;
113
114 /* node constructor input: */
115 int node_id;
116 int node_partitionborder=0;
117 double node_x[3];
118 int node_onbed;
119 int node_onsurface;
120 int node_numdofs;
121
122
123 #ifdef _PARALLEL_
124 /*Metis partitioning: */
125 int range;
126 Vec gridborder;
127 int my_numgrids;
128 int* all_numgrids=NULL;
129 int gridcount;
130 int count;
131 #endif
132 int first_grid_index;
133
134 /*Rifts:*/
135 int* riftsnumpenaltypairs;
136 double** riftspenaltypairs;
137 int* riftsfill;
138 double* riftsfriction;
139 double* riftpenaltypairs=NULL;
140 int el1,el2;
141
142 /*Penalty partitioning: */
143 int num_grids3d_collapsed;
144 double* double_penalties_grids3d_collapsed=NULL;
145 double* double_penalties_grids3d_noncollapsed=NULL;
146 int grid_ids[6];
147 int num_grid_ids;
148 int grid_id;
149
150 /*Get analysis_type: */
151 analysis_type=AnalysisTypeAsEnum(model->analysis_type);
152
153 /*First create the elements, nodes and material properties: */
154 elements = new DataSet(ElementsEnum());
155 nodes = new DataSet(NodesEnum());
156 materials = new DataSet(MaterialsEnum());
157
158 /*Width of elements: */
159 if(strcmp(model->meshtype,"2d")==0){
160 elements_width=3; //tria elements
161 }
162 else{
163 elements_width=6; //penta elements
164 }
165
166
167 #ifdef _PARALLEL_
168 /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
169 if(strcmp(model->meshtype,"2d")==0){
170 /*load elements: */
171 ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
172 }
173 else{
174 /*load elements2d: */
175 ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
176 }
177
178 MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs);
179
180 /*Free elements and elements2d: */
181 xfree((void**)&model->elements);
182 xfree((void**)&model->elements2d);
183
184
185 /*Deal with rifts, they have to be included into one partition only, not several: */
186 FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts);
187
188 for(i=0;i<model->numrifts;i++){
189 riftpenaltypairs=model->riftspenaltypairs[i];
190 for(j=0;j<model->riftsnumpenaltypairs[i];j++){
191 el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing
192 el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing
193 epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids;
194 }
195 }
196 /*Free rifts: */
197 xfree((void**)&riftsnumpenaltypairs);
198 for(i=0;i<model->numrifts;i++){
199 double* temp=riftspenaltypairs[i];
200 xfree((void**)&temp);
201 }
202 xfree((void**)&riftspenaltypairs);
203 xfree((void**)&riftsfill);
204 xfree((void**)&riftsfriction);
205
206 /*Used later on: */
207 my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int));
208 #endif
209
210
211
212 /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */
213
214 /*2d mesh: */
215 if (strcmp(model->meshtype,"2d")==0){
216
217 /*Fetch data needed: */
218 ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
219 ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
220 ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
221 ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
222 ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat");
223 ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat");
224 ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat");
225 ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
226 ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat");
227 ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat");
228
229 for (i=0;i<model->numberofelements;i++){
230
231 #ifdef _PARALLEL_
232 /*!All elements have been partitioned above, only create elements for this CPU: */
233 if(my_rank==epart[i]){
234 #endif
235
236
237 /*ids: */
238 tria_id=i+1; //matlab indexing.
239 tria_mid=i+1; //refers to the corresponding material property card
240 tria_mparid=model->numberofelements+1;//refers to the corresponding parmat property card
241
242 /*vertices offsets: */
243 tria_g[0]=(int)*(model->elements+elements_width*i+0);
244 tria_g[1]=(int)*(model->elements+elements_width*i+1);
245 tria_g[2]=(int)*(model->elements+elements_width*i+2);
246
247 /*thickness,surface and bed:*/
248 tria_h[0]= *(model->thickness+ ((int)*(model->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
249 tria_h[1]=*(model->thickness+ ((int)*(model->elements+elements_width*i+1)-1));
250 tria_h[2]=*(model->thickness+ ((int)*(model->elements+elements_width*i+2)-1)) ;
251
252 tria_s[0]=*(model->surface+ ((int)*(model->elements+elements_width*i+0)-1));
253 tria_s[1]=*(model->surface+ ((int)*(model->elements+elements_width*i+1)-1));
254 tria_s[2]=*(model->surface+ ((int)*(model->elements+elements_width*i+2)-1));
255
256 tria_b[0]=*(model->bed+ ((int)*(model->elements+elements_width*i+0)-1));
257 tria_b[1]=*(model->bed+ ((int)*(model->elements+elements_width*i+1)-1));
258 tria_b[2]=*(model->bed+ ((int)*(model->elements+elements_width*i+2)-1));
259
260 /*basal drag:*/
261 tria_friction_type=(int)model->drag_type;
262
263 tria_k[0]=*(model->drag+ ((int)*(model->elements+elements_width*i+0)-1));
264 tria_k[1]=*(model->drag+ ((int)*(model->elements+elements_width*i+1)-1));
265 tria_k[2]=*(model->drag+ ((int)*(model->elements+elements_width*i+2)-1));
266
267 tria_p=model->p[i];
268 tria_q=model->q[i];
269
270 /*element on iceshelf?:*/
271 tria_shelf=(int)*(model->elementoniceshelf+i);
272
273 tria_meanvel=model->meanvel;
274 tria_epsvel=model->epsvel;
275
276 /*viscosity_overshoot*/
277 tria_viscosity_overshoot=model->viscosity_overshoot;
278
279 /*Create tria element using its constructor:*/
280 tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot);
281
282 /*Add tria element to elements dataset: */
283 elements->AddObject(tria);
284
285 /*Deal with material property card: */
286 matice_mid=i+1; //same as the material id from the geom2 elements.
287
288 /*Average B over 3 grid elements: */
289 B_avg=0;
290 for(j=0;j<3;j++){
291 B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1));
292 }
293 B_avg=B_avg/3;
294 matice_B=B_avg;
295 matice_n=(double)*(model->n+i);
296
297 /*Create matice using its constructor:*/
298 matice= new Matice(matice_mid,matice_B,matice_n);
299
300 /*Add matice element to materials dataset: */
301 materials->AddObject(matice);
302
303 #ifdef _PARALLEL_
304 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
305 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
306 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
307 will hold which grids belong to this partition*/
308 my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
309 my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
310 my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
311 #endif
312
313 #ifdef _PARALLEL_
314 }//if(my_rank==epart[i])
315 #endif
316
317 }//for (i=0;i<numberofelements;i++)
318
319
320 /*Free data : */
321 /*xfree((void**)&model->elements);
322 xfree((void**)&model->thickness);
323 xfree((void**)&model->surface);
324 xfree((void**)&model->bed);
325 xfree((void**)&model->drag);
326 xfree((void**)&model->p);
327 xfree((void**)&model->q);
328 xfree((void**)&model->elementoniceshelf);
329 xfree((void**)&model->B);
330 xfree((void**)&model->n);*/
331 }
332 else{ // if (strcmp(meshtype,"2d")==0)
333
334 /*Fetch data needed: */
335 ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
336 ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
337 ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
338 ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
339 ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat");
340 ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat");
341 ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat");
342 ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
343 ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat");
344 ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat");
345 ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
346 ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat");
347 ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat");
348
349 for (i=0;i<model->numberofelements;i++){
350 #ifdef _PARALLEL_
351 /*We are using our element partition to decide which elements will be created on this node: */
352 if(my_rank==epart[i]){
353 #endif
354
355
356 /*name and id: */
357 penta_id=i+1; //matlab indexing.
358 penta_mid=i+1; //refers to the corresponding material property card
359 penta_mparid=model->numberofelements+1;//refers to the corresponding parmat property card
360
361 /*vertices,thickness,surface,bed and drag: */
362 for(j=0;j<6;j++){
363 penta_g[j]=(int)*(model->elements+elements_width*i+j);
364 penta_h[j]=*(model->thickness+ ((int)*(model->elements+elements_width*i+j)-1));
365 penta_s[j]=*(model->surface+ ((int)*(model->elements+elements_width*i+j)-1));
366 penta_b[j]=*(model->bed+ ((int)*(model->elements+elements_width*i+j)-1));
367 penta_k[j]=*(model->drag+ ((int)*(model->elements+elements_width*i+j)-1));
368 }
369
370 /*basal drag:*/
371 penta_friction_type=(int)model->drag_type;
372
373 penta_p=model->p[i];
374 penta_q=model->q[i];
375
376 /*diverse: */
377 penta_shelf=(int)*(model->elementoniceshelf+i);
378 penta_onbed=(int)*(model->elementonbed+i);
379 penta_onsurface=(int)*(model->elementonsurface+i);
380 penta_meanvel=model->meanvel;
381 penta_epsvel=model->epsvel;
382
383 if (*(model->elements_type+2*i+0)==MacAyealEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base.
384 penta_collapse=1;
385 }
386 else{
387 penta_collapse=0;
388 }
389
390
391 /*Create Penta using its constructor:*/
392 penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
393 penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
394 penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
395 penta_thermal_steadystate);
396
397 /*Add penta element to elements dataset: */
398 elements->AddObject(penta);
399
400
401 /*Deal with material:*/
402 matice_mid=i+1; //same as the material id from the geom2 elements.
403 /*Average B over 6 element grids: */
404 B_avg=0;
405 for(j=0;j<6;j++){
406 B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1));
407 }
408 B_avg=B_avg/6;
409 matice_B= B_avg;
410 matice_n=(double)*(model->n+i);
411
412 /*Create matice using its constructor:*/
413 matice= new Matice(matice_mid,matice_B,matice_n);
414
415 /*Add matice element to materials dataset: */
416 materials->AddObject(matice);
417
418 #ifdef _PARALLEL_
419 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
420 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
421 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
422 will hold which grids belong to this partition*/
423 my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
424 my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
425 my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
426 my_grids[(int)*(model->elements+elements_width*i+3)-1]=1;
427 my_grids[(int)*(model->elements+elements_width*i+4)-1]=1;
428 my_grids[(int)*(model->elements+elements_width*i+5)-1]=1;
429 #endif
430
431 #ifdef _PARALLEL_
432 }//if(my_rank==epart[i])
433 #endif
434
435 }//for (i=0;i<numberofelements;i++)
436
437 /*Free data: */
438 xfree((void**)&model->elements);
439 xfree((void**)&model->thickness);
440 xfree((void**)&model->surface);
441 xfree((void**)&model->bed);
442 xfree((void**)&model->drag);
443 xfree((void**)&model->p);
444 xfree((void**)&model->q);
445 xfree((void**)&model->elementoniceshelf);
446 xfree((void**)&model->elementonbed);
447 xfree((void**)&model->elementonsurface);
448 xfree((void**)&model->elements_type);
449 xfree((void**)&model->n);
450 xfree((void**)&model->B);
451
452 } //if (strcmp(meshtype,"2d")==0)
453
454 /*Add one constant material property to materials: */
455 matpar_mid=model->numberofelements+1; //put it at the end of the materials
456 matpar_g=model->g;
457 matpar_viscosity_overshoot=model->viscosity_overshoot;
458 matpar_rho_ice=model->rho_ice;
459 matpar_rho_water=model->rho_water;
460 matpar_thermalconductivity=model->thermalconductivity;
461 matpar_heatcapacity=model->heatcapacity;
462 matpar_latentheat=model->latentheat;
463 matpar_beta=model->beta;
464 matpar_meltingpoint=model->meltingpoint;
465 matpar_mixed_layer_capacity=model->mixed_layer_capacity;
466 matpar_thermal_exchange_velocity=model->thermal_exchange_velocity;
467
468 /*Create matpar object using its constructor: */
469 matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,
470 matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,
471 matpar_thermal_exchange_velocity,matpar_g,matpar_viscosity_overshoot);
472
473 /*Add to materials datset: */
474 materials->AddObject(matpar);
475
476 #ifdef _PARALLEL_
477 /*From the element partitioning, we can determine which grids are on the inside of this cpu's
478 *element partition, and which are on its border with other nodes:*/
479 gridborder=NewVec(model->numberofnodes);
480
481 for (i=0;i<model->numberofnodes;i++){
482 if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
483 }
484 VecAssemblyBegin(gridborder);
485 VecAssemblyEnd(gridborder);
486
487 #ifdef _DEBUG_
488 VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
489 #endif
490
491 VecToMPISerial(&my_bordergrids,gridborder);
492
493 #ifdef _DEBUG_
494 if(my_rank==0){
495 for (i=0;i<model->numberofnodes;i++){
496 printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
497 }
498 }
499 #endif
500 #endif
501
502 /*Partition penalties in 3d: */
503 if(strcmp(model->meshtype,"3d")==0){
504
505 /*Get penalties: */
506 ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
507
508 if(model->numpenalties){
509
510 model->penaltypartitioning=(int*)xmalloc(model->numpenalties*sizeof(int));
511 for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=-1;
512
513 for(i=0;i<model->numpenalties;i++){
514 first_grid_index=(int)(*(model->penalties+i*model->numlayers+0)-1);
515 if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition grids
516 /*All grids that are being penalised belong to this node's internal grid partition.:*/
517 model->penaltypartitioning[i]=1;
518 }
519 if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
520 model->penaltypartitioning[i]=0;
521 }
522 }
523 }
524
525 /*Free penalties: */
526 xfree((void**)&model->penalties);
527 }
528
529 /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids.
530 We can therefore determine which grids are internal to this node's partition
531 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing
532 that, go and create the grids*/
533
534 /*Create nodes from x,y,z, as well as the spc values on those grids: */
535
536 /*First fetch data: */
537 if (strcmp(model->meshtype,"3d")==0){
538 ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat");
539 }
540 ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat");
541 ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat");
542 ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat");
543 ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
544 ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
545
546
547 /*Get number of dofs per node: */
548 DistributeNumDofs(&node_numdofs,analysis_type);
549
550 for (i=0;i<model->numberofnodes;i++){
551 #ifdef _PARALLEL_
552 /*keep only this partition's nodes:*/
553 if((my_grids[i]==1)){
554 #endif
555
556 node_id=i+1; //matlab indexing
557
558
559
560 #ifdef _PARALLEL_
561 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
562 node_partitionborder=1;
563 }
564 else{
565 node_partitionborder=0;
566 }
567 #else
568 node_partitionborder=0;
569 #endif
570
571
572 node_x[0]=model->x[i];
573 node_x[1]=model->y[i];
574 node_x[2]=model->z[i];
575
576 node_onbed=(int)model->gridonbed[i];
577 node_onsurface=(int)model->gridonsurface[i];
578
579 /*Create node using its constructor: */
580 node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface);
581
582 /*set single point constraints.: */
583 if (strcmp(model->meshtype,"3d")==0){
584 /*On a 3d mesh, we may have collapsed grids. Spc all their dofs: */
585 if (model->deadgrids[i]){
586 for(k=0;k<node_numdofs;k++){
587 node->DofInSSet(k);
588 }
589 }
590 }
591 /*Add node to nodes dataset: */
592 nodes->AddObject(node);
593
594 #ifdef _PARALLEL_
595 } //if((my_grids[i]==1))
596 #endif
597 }
598
599 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
600 * datasets, it will not be redone: */
601 elements->Presort();
602 nodes->Presort();
603 materials->Presort();
604
605 /*Clean fetched data: */
606 xfree((void**)&model->deadgrids);
607 xfree((void**)&model->x);
608 xfree((void**)&model->y);
609 xfree((void**)&model->z);
610 xfree((void**)&model->gridonbed);
611 xfree((void**)&model->gridonsurface);
612
613 cleanup_and_return:
614
615 /*Assign output pointer: */
616 *pelements=elements;
617 *pnodes=nodes;
618 *pmaterials=materials;
619
620 /*Keep partitioning information into model*/
621 model->epart=epart;
622 model->my_grids=my_grids;
623 model->my_bordergrids=my_bordergrids;
624
625 /*Free ressources:*/
626 #ifdef _PARALLEL_
627 xfree((void**)&all_numgrids);
628 VecFree(&gridborder);
629 #endif
630
631 return noerr;
632
633}
Note: See TracBrowser for help on using the repository browser.