source: issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp@ 1310

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

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