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

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

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