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

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

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