CoCalc -- egnative.c

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/elmergrid/src/egnative.c
⁵²²⁸ views
1
/*   
2
   ElmerGrid - A simple mesh generation and manipulation utility  
3
   Copyright (C) 1995- , CSC - IT Center for Science Ltd.   
4

5
   Author: Peter Raback
6
   Email: [email protected]
7
   Address: CSC - IT Center for Science Ltd.
8
            Keilaranta 14
9
            02101 Espoo, Finland
10

11
   This program is free software; you can redistribute it and/or
12
   modify it under the terms of the GNU General Public License
13
   as published by the Free Software Foundation; either version 2
14
   of the License, or (at your option) any later version.
15

16
   This program is distributed in the hope that it will be useful,
17
   but WITHOUT ANY WARRANTY; without even the implied warranty of
18
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19
   GNU General Public License for more details.
20

21
   You should have received a copy of the GNU General Public License
22
   along with this program; if not, write to the Free Software
23
   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
24
*/
25

26
/* -------------------------------:  egnative.c  :----------------------------
27
   This module includes routines for I/O of native formats of Elmer. 
28
*/
29

30
#include <stdio.h>
31
#include <unistd.h>
32
#include <stddef.h>
33
#include <stdlib.h>
34
#include <string.h>
35
#include <math.h>
36

37
#if HAVE_UNISTD_H
38
#include <unistd.h>
39
#endif
40

41
#include <stdarg.h>
42
#include <stdio.h>
43
#include <limits.h>
44
#include <ctype.h>
45
#include <sys/stat.h>
46
#include <sys/types.h>
47

48
#include "egutils.h"
49
#include "egdef.h"
50
#include "egtypes.h"
51
#include "egmesh.h"
52
/* #include "egparallel.h" */
53
#include "egnative.h"
54
/*#include "../config.h"*/
55

56
#define GETLINE ioptr=fgets(line,MAXLINESIZE,in) 
57
static char *ioptr;
58

59

60
#define DEBUG 0
61

62

63
int matcactive=FALSE, iodebug=FALSE;
64

65
#define MAXINMETHODS 21
66
const char *InMethods[] = {
67
  /*0*/ "EG",
68
  /*1*/ "ELMERGRID",
69
  /*2*/ "ELMERSOLVER",
70
  /*3*/ "ELMERPOST",
71
  /*4*/ "ANSYS",
72
  /*5*/ "IDEAS",
73
  /*6*/ "ABAQUS",
74
  /*7*/ "FIDAP",
75
  /*8*/ "UNV",
76
  /*9*/ "COMSOL",
77
  /*10*/ "FIELDVIEW",
78
  /*11*/ "TRIANGLE",
79
  /*12*/ "MEDIT",
80
  /*13*/ "GID",
81
  /*14*/ "GMSH",
82
  /*15*/ "PARTITIONED",
83
  /*16*/ "FVCOM",
84
  /*17*/ "NASTRAN",
85
  /*18*/ "CGSIM",
86
  /*19*/ "GEO",
87
  /*20*/ "FLUX2D",
88
  /*21*/ "FLUX3D",
89
};
90

91

92
#define MAXOUTMETHODS 5
93
const char *OutMethods[] = {
94
  /*0*/ "EG",
95
  /*1*/ "ELMERGRID",
96
  /*2*/ "ELMERSOLVER",
97
  /*3*/ "ELMERPOST",
98
  /*4*/ "GMSH",
99
  /*5*/ "VTU",
100
};
101

102

103
void Instructions()
104
{
105
  printf("This program can create simple 2D structured meshes consisting of\n");
106
  printf("linear, quadratic or cubic rectangles or triangles. The meshes may\n");
107
  printf("also be extruded and revolved to create 3D forms. In addition many\n");
108
  printf("mesh formats may be imported into Elmer software. Some options have\n");
109
  printf("not been properly tested. Contact the author if you face problems.\n");
110
  printf("------------------------------------------------------------------\n");
111

112
  printf("The program has two operation modes\n");
113
  printf("A) Command file mode which has the command file as the only argument\n");
114
  printf("   'ElmerGrid commandfile.eg'\n\n");
115

116
  printf("------------------------------------------------------------------\n");
117
  printf("B) Inline mode which expects at least three input parameters\n");
118
  printf("   'ElmerGrid 1 3 test'\n\n");
119
  printf("The first parameter defines the input file format:\n");
120
  printf("1)  .grd      : ElmerGrid file format\n");
121
  printf("2)  .mesh.*   : Elmer input format\n");
122
  printf("3)  .ep       : Elmer output format\n");
123
  printf("4)  .ansys    : Ansys input format\n");
124
  printf("5)  .inp      : Abaqus input format by Ideas\n");
125
  printf("6)  .fil      : Abaqus output format\n");
126
  printf("7)  .FDNEUT   : Gambit (Fidap) neutral file\n");
127
  printf("8)  .unv      : Universal mesh file format\n");
128
  printf("9)  .mphtxt   : Comsol Multiphysics mesh format\n");
129
  printf("10) .dat      : Fieldview format\n");
130
  printf("11) .node,.ele: Triangle 2D mesh format\n");
131
  printf("12) .mesh     : Medit mesh format\n");
132
  printf("13) .msh      : GID mesh format\n");
133
  printf("14) .msh      : Gmsh mesh format\n");
134
  printf("15) .ep.i     : Partitioned ElmerPost format\n");
135
  printf("16) .2dm      : 2D triangular FVCOM format\n");
136
#if 0
137
  printf("17) .msh      : Nastran format\n");
138
  printf("18) .msh      : CGsim format\n");
139
  printf("19) .geo      : Geo format\n");
140
  printf("20) .tra      : Cedrat Flux 2D format\n");
141
  printf("21) .pf3      : Cedrat Flux 3D format\n");
142
#endif 
143

144
  printf("\nThe second parameter defines the output file format:\n");
145
  printf("1)  .grd      : ElmerGrid file format\n");
146
  printf("2)  .mesh.*   : ElmerSolver format (also partitioned .part format)\n");
147
  printf("3)  .ep       : ElmerPost format\n");
148
  printf("4)  .msh      : Gmsh mesh format\n");
149
  printf("5)  .vtu      : VTK ascii XML format\n");
150
#if 0
151
  printf("5)  .inp      : Abaqus input format\n");
152
  printf("7)  .fidap    : Fidap format\n");
153
  printf("18) .ep       : Fastcap input format.\n");
154
#endif
155
  printf("------------------------------------------------------------------\n");
156

157
  printf("\nThe third parameter is the name of the input file.\n");
158
  printf("If the file does not exist, an example with the same name is created.\n");
159
  printf("The default output file name is the same with a different suffix.\n\n");
160

161
  printf("There are several additional in-line parameters that are\n");
162
  printf("taken into account only when applicable to the given format.\n");
163

164
  printf("-out str             : name of the output file\n");
165
  printf("-in str              : name of a secondary input file\n");
166
  printf("-decimals            : number of decimals in the saved mesh (eg. 8)\n");
167
  printf("-bin                 : save mesh for ElmerSolver in binary format\n");
168
  printf("-sbin                : save mesh for ElmerSolver in binary & single precision format\n");
169
  printf("-relh real           : give relative mesh density parameter for ElmerGrid meshing\n");
170
  printf("-triangles           : rectangles will be divided to triangles\n");
171
  printf("-merge real          : merges nodes that are close to each other\n");
172
  printf("-order real[3]       : reorder elements and nodes using c1*x+c2*y+c3*z\n");
173
  printf("-centralize          : set the center of the mesh to origin\n");
174
  printf("-scale real[3]       : scale the coordinates with vector real[3]\n");
175
  printf("-translate real[3]   : translate the nodes with vector real[3]\n");
176
  printf("-rotate real[3]      : rotate around the main axis with angles real[3]\n");
177
  printf("-clone int[3]        : make identical copies of the mesh\n");
178
  printf("-clonesize real[3]   : the size of the mesh to be cloned if larger to the original\n");
179
  printf("-mirror int[3]       : copy the mesh around the origin in coordinate directions\n");
180
  printf("-cloneinds           : when performing cloning should cloned entities be given new indexes\n");
181
  printf("-unite               : the meshes will be united\n");
182
  printf("-unitenooverlap      : the meshes will be united without overlap in entity numbering\n");
183
  printf("-polar real          : map 2D mesh to a cylindrical shell with given radius\n");
184
  printf("-cylinder            : map 2D/3D cylindrical mesh to a cartesian mesh\n");
185
  printf("-reduce int[2]       : reduce element order at material interval [int1 int2]\n");
186
  printf("-increase            : increase element order from linear to quadratic\n");
187
  printf("-bcoffset int        : add an offset to the boundary conditions\n");
188
  printf("-discont int         : make the boundary to have secondary nodes\n");
189
  printf("-connect int         : make the boundary to have internal connection among its elements\n");
190
  printf("-removeintbcs        : remove internal boundaries if they are not needed\n");
191
  printf("-removelowdim        : remove boundaries that are two ranks lower than highest dim\n");
192
  printf("-removeunused        : remove nodes that are not used in any element\n");
193
  printf("-bulkorder           : renumber materials types from 1 so that every number is used\n");
194
  printf("-boundorder          : renumber boundary types from 1 so that every number is used\n");
195
  printf("-autoclean           : this performs the united action of the four above\n");
196
  printf("-multidim            : keep lower order entities even if they are not boundaries\n");
197
  printf("-bulkbound int[3]    : set the intersection of materials [int1 int2] to be boundary int3\n");
198
  printf("-boundbound int[3]   : set the intersection of boundaries [int1 int2] to be boundary int3\n");
199
  printf("-bulktype int[3]     : set material types in interval [int1 int2] to type int3\n");
200
  printf("-boundtype int[3]    : set sidetypes in interval [int1 int2] to type int3\n");
201
  printf("-layer int[2] real[2]: make a boundary layer for given boundary\n");
202
  printf("-layermove int       : apply Jacobi filter int times to move the layered mesh\n");
203
  printf("-divlayer int[2] real[2]: make a boundary layer for given boundary\n");
204
  printf("-3d / -2d / -1d      : mesh is 3, 2 or 1-dimensional (applies to examples)\n");
205
  printf("-isoparam            : ensure that higher order elements are convex\n");
206
  printf("-nonames             : disable use of mesh.names even if it would be supported by the format\n");
207
  printf("-nosave              : disable saving part altogether\n");
208
  printf("-nooverwrite         : if mesh already exists don't overwrite it\n");
209
  printf("-vtuone              : start real node indexes in vtu file from one\n");
210
  printf("-timer               : show timer information\n");
211
  printf("-infofile str        : file for saving the timer and size information\n");
212

213
  printf("\nKeywords are related to mesh partitioning for parallel ElmerSolver runs:\n");
214
  printf("-partition int[3]    : the mesh will be partitioned in cartesian main directions\n");
215
  printf("-partorder real[3]   : in the 'partition' method set the direction of the ordering\n");
216
  printf("-parttol real        : in the 'partition' method set the tolerance for ordering\n");
217
  printf("-partcell int[3]     : the mesh will be partitioned in cells of fixed sizes\n");
218
  printf("-partcyl int[3]      : the mesh will be partitioned in cylindrical main directions\n");
219
#if USE_METIS
220
  printf("-metiskway int       : mesh will be partitioned with Metis using graph Kway routine\n");
221
  printf("-metisrec int        : mesh will be partitioned with Metis using graph Recursive routine\n");
222
  printf("-metiscontig         : enforce that the metis partitions are contiguous\n");
223
  printf("-metisvol            : minimize total communication volume in Metis\n");
224
  printf("-metisminconn        : minimize the maximum connectivity count in Metis\n");
225
  printf("-metisseed int       : random number generator seed for Metis algorithms\n");
226
  printf("-metisncuts int      : number of competing partitions to generate\n");
227
#endif
228
  printf("-partdual            : use the dual graph in partition method (when available)\n");
229
  printf("-halo                : create halo for the partitioning for DG\n");
230
  printf("-halobc              : create halo for the partitioning at boundaries only\n");
231
  printf("-haloz / -halor      : create halo for the the special z- or r-partitioning\n");
232
  printf("-halogreedy          : create halo being greedy over the partition interfaces\n");
233
  printf("-periodic int[3]     : periodic coordinate directions for parallel & conforming meshes\n");
234
  printf("-partoptim           : apply aggressive optimization to node sharing\n");
235
  printf("-partnobcoptim       : do not apply optimization to bc ownership sharing\n");
236
  printf("-partbw              : minimize the bandwidth of partition-partition couplings\n");
237
  printf("-parthypre           : number the nodes continuously partitionwise\n");
238
  printf("-partzbc             : partition connected BCs separately to partitions in Z-direction\n");
239
  printf("-partrbc             : partition connected BCs separately to partitions in R-direction\n");
240
#if USE_METIS
241
  printf("-metisbc             : partition connected BCs separately to partitions by Metis\n");
242
#endif
243
  printf("-partlayers int      : extend boundary partitioning by element layers\n");
244
  
245
  
246
#if 0
247
  printf("\nKeywords are related to (nearly obsolete) ElmerPost format:\n");
248
  printf("-partjoin int        : number of ElmerPost partitions in the data to be joined\n");
249
  printf("-saveinterval int[3] : the first, last and step for fusing parallel data\n");
250
  printf("-nobound             : disable saving of boundary elements in ElmerPost format\n");
251
#endif
252
  
253
  if(0) printf("-names               : conserve name information where applicable\n");
254
  printf("------------------------------------------------------------------\n");
255
}
256
 
257

258
void Goodbye()
259
{
260
  printf("\nThank you for using Elmergrid!\n");
261
  printf("Send bug reports and feature wishes to [email protected]\n");
262
  exit(0);
263
}
264

265

266

267
#if USE_MATC
268
char *mtc_domath(const char *);
269
void mtc_init(FILE * input, FILE * output, FILE * error);
270
#endif
271

272

273
static int Getline(char *line1,FILE *io) 
274
{
275
  int i,isend;
276
  char line0[MAXLINESIZE],*charend,*matcpntr,*matcpntr0;
277

278
  for(i=0;i<MAXLINESIZE;i++)
279
    line0[i] = 0x00;
280

281
 newline:
282

283
  charend = fgets(line0,MAXLINESIZE,io);
284
  isend = (charend == NULL);
285

286
  if(isend) return(1);
287

288
  if(line0[0] == '#' || line0[0] == '%' || line0[0] == '!') goto newline;
289
  if(!matcactive && line0[0] == '*') goto newline;
290

291
  if(!matcactive && strchr(line0,'$') ) {
292
    #if USE_MATC
293
      printf("\n a $ found and MATC has been compiled but not activated,\n");
294
      printf("either remove the $ or add 'MATC = True' to grd input file.\n");
295
    #else
296
      printf("\n a $ found and MATC has not been compiled into ElmerGrid,\n");
297
      printf("either remove the $ or compile ElmerGrid with MATC\n");
298
    #endif // USE_MATC
299
      bigerror("Error: $ found in command and MATC is not active");
300
  }
301

302
#if USE_MATC
303
  if(matcactive) {
304
    matcpntr0 = strchr(line0,'$');
305
    if(matcpntr0) {
306
      matcpntr = mtc_domath(&matcpntr0[1]);
307
      if(matcpntr) {
308
	strcpy(matcpntr0, matcpntr);
309
	if(0) printf("A: %s\n%s\n",matcpntr0,matcpntr);
310
      }
311
    }
312
  }
313
#endif 
314

315
  if(strstr(line0,"subcell boundaries")) goto newline;
316
  if(strstr(line0,"material structure")) goto newline;
317
  if(strstr(line0,"mode")) goto newline;
318
  if(strstr(line0,"type")) goto newline;
319

320
  for(i=0;i<MAXLINESIZE;i++) 
321
    line1[i] = toupper(line0[i]);
322

323
  if(iodebug) printf("line: %s\n",line1);
324

325
  return(0);
326
}
327

328

329

330
void InitGrid(struct GridType *grid)
331
/* Initializes the grid of a specific mesh. A grid can differ 
332
   between various differential equations. 
333
   */
334
#define MAXBODYID 1000
335
{
336
  int i,j;
337

338
  grid->layered = FALSE;
339
  grid->layeredbc = TRUE;
340
  grid->layerbcoffset = 0;
341
  grid->triangles = FALSE;
342
  grid->triangleangle = 0.0;
343
  grid->partitions = FALSE;
344
  grid->wantedelems = 0;
345
  grid->wantedelems3d = 0;
346
  grid->wantednodes3d = 0;
347
  grid->limitdx = 0.0;
348
  grid->limitdxverify = FALSE;
349

350
  grid->dimension = 2;
351
  grid->xcells = grid->ycells = grid->zcells = 0;
352
  grid->nocells = 0;
353
  grid->noknots = grid->noelements = 0;
354
  grid->totxelems = grid->totyelems = grid->totzelems = 0;
355
  grid->minxelems = grid->minyelems = 1;
356
  grid->minzelems = 2;
357
  grid->firstmaterial = 1;
358
  grid->lastmaterial = MAXBODYID;
359
  grid->autoratio = 1;
360
  grid->xyratio = 1.0;
361
  grid->xzratio = 1.0;
362
  grid->nonodes = 4;
363
  grid->numbering = NUMBER_XY;
364
  grid->rotate = FALSE;
365
  grid->rotateradius1 = 0.0;
366
  grid->rotateradius2 = 0.0;
367
  grid->rotateimprove = 1.0;
368
  grid->rotateblocks = 4;
369
  grid->rotatecartesian = FALSE;
370
  grid->reduceordermatmin = 0;
371
  grid->reduceordermatmax = 0;
372
  grid->polarradius = 1.0;
373
  grid->elemorder = 1;
374
  grid->elemmidpoints = FALSE;
375

376
  grid->rotatecurve = FALSE;
377
  grid->curverad = 0.5;
378
  grid->curveangle = 90.0;
379
  grid->curvezet = 1.0;
380

381
  for(i=0;i<=MAXCELLS;i++) {
382
    grid->xelems[i] = 0;
383
    grid->xlinear[i] = TRUE;
384
    grid->xratios[i] = 1.;
385
    grid->xexpand[i] = 1.;
386
    grid->xdens[i] = 1.;
387
    grid->x[i] = 0.; 
388
  }
389
  for(i=0;i<=MAXCELLS;i++) {
390
    grid->yelems[i] = 0;
391
    grid->ylinear[i] = TRUE;
392
    grid->yratios[i] = 1.0;
393
    grid->yexpand[i] = 1.0;
394
    grid->ydens[i] = 1.0;
395
    grid->y[i] = 0.;
396
  }
397
  for(i=0;i<=MAXCELLS;i++) {
398
    grid->zelems[i] = 0;
399
    grid->zlinear[i] = TRUE;
400
    grid->zratios[i] = 1.0;
401
    grid->zexpand[i] = 1.0;
402
    grid->zdens[i] = 1.0;
403
    grid->z[i] = 0.;
404
    grid->zfirstmaterial[i] = 1;
405
    grid->zlastmaterial[i] = MAXBODYID;
406
    grid->zmaterial[i] = 0; 
407
  }
408

409
  grid->zmaterialmapexists = FALSE;
410
  grid->zhelicityexists = FALSE;
411
  grid->zhelicity = 0.0;
412

413
  /* Initializes the numbering of the cells. */
414
  for(j=0;j<=MAXCELLS+1;j++)
415
    for(i=0;i<=MAXCELLS+1;i++) 
416
      grid->structure[j][i] = MAT_NOTHING; 
417

418
  for(j=0;j<=MAXCELLS+1;j++)
419
    for(i=0;i<=MAXCELLS+1;i++) 
420
      grid->numbered[j][i] = 0; 
421

422
  grid->noboundaries = 0;
423
  for(i=0;i<MAXBOUNDARIES;i++) {
424
    grid->boundint[i] = 0;
425
    grid->boundext[i] = 0;
426
    grid->boundsolid[i] = 0;
427
    grid->boundtype[i] = 0;
428
  }
429

430
  grid->mappings = 0;
431
  for(i=0;i<MAXMAPPINGS;i++) {
432
    grid->mappingtype[i] = 0;
433
    grid->mappingline[i] = 0;    
434
    grid->mappinglimits[2*i] = 0;    
435
    grid->mappinglimits[2*i+1] = 0;    
436
    grid->mappingpoints[i] = 0;
437
  }
438
}
439

440

441
static void ExampleGrid1D(struct GridType **grids,int *nogrids,int info) 
442
/* Creates an example grid that might be used to analyze 
443
   flow through a step. */
444
{
445
  int j;
446
  struct GridType *grid;
447

448
  (*nogrids) = 1;
449
  (*grids) = (struct GridType*)
450
    malloc((size_t) (*nogrids)*sizeof(struct GridType)); 
451

452
  grid = grids[0];
453

454
  InitGrid(grid);
455

456
  grid->nonodes = 2;
457
  grid->xcells = 3;
458
  grid->ycells = 1;
459
  grid->wantedelems = 40;
460
  grid->firstmaterial = 1;
461
  grid->lastmaterial = 1;
462
  grid->autoratio = 1;
463
  grid->coordsystem = COORD_CART1;
464
  grid->numbering = NUMBER_1D;
465

466
  grid->x[0] = 0.0;
467
  grid->x[1] = 1.0;
468
  grid->x[2] = 3.0;
469
  grid->x[3] = 4.0;
470

471
  grid->xexpand[1] = 2.0;
472
  grid->xexpand[3] = 0.5;
473
  grid->xdens[2] = 0.5;
474

475
  for(j=1;j<=3;j++)
476
    grid->structure[1][j] = 1;
477

478
  grid->noboundaries = 2;
479
  grid->boundint[0] = 1;
480
  grid->boundint[1] = 1;
481

482
  grid->boundext[0] = -1;
483
  grid->boundext[1] = -2;
484

485
  grid->boundtype[0] = 1;
486
  grid->boundtype[1] = 2;
487

488
  grid->boundsolid[0] = 1;
489
  grid->boundsolid[1] = 1;
490

491
  grid->mappings = 0;
492
  
493
  if(info) printf("A simple 1D example mesh was created\n");
494
}
495

496

497
static void ExampleGrid2D(struct GridType **grids,int *nogrids,int info) 
498
/* Creates an example grid that might be used to analyze 
499
   flow through a step. */
500
{
501
  int j;
502
  struct GridType *grid;
503

504
  (*nogrids) = 1;
505
  (*grids) = (struct GridType*)
506
    malloc((size_t) (*nogrids)*sizeof(struct GridType)); 
507

508
  grid = grids[0];
509

510
  InitGrid(grid);
511

512
  grid->xcells = 4;
513
  grid->ycells = 3;
514
  grid->wantedelems = 100;
515
  grid->totxelems = grid->totyelems = grid->totzelems = 10;
516
  grid->firstmaterial = 1;
517
  grid->lastmaterial = 1;
518
  grid->autoratio = 1;
519
  grid->coordsystem = COORD_CART2;
520
  grid->numbering = NUMBER_XY;
521

522
  grid->x[0] = -1.0;
523
  grid->x[1] = 0.0;
524
  grid->x[2] = 2.0;
525
  grid->x[3] = 6.0;
526
  grid->x[4] = 7.0;
527

528
  grid->y[0] = 0.0;
529
  grid->y[1] = 0.5;
530
  grid->y[2] = 0.75;
531
  grid->y[3] = 1.0;
532

533
  grid->xexpand[2] = 0.4;
534
  grid->xexpand[3] = 5.0;
535

536
  grid->yexpand[2] = 2.0;
537
  grid->yexpand[3] = 0.5;
538

539
  for(j=1;j<=3;j++)
540
    grid->structure[j][1] = 2;
541
  for(j=1;j<=3;j++)
542
    grid->structure[j][2] = 1;
543
  grid->structure[1][2] = 0;
544
  for(j=1;j<=3;j++)
545
    grid->structure[j][3] = 1;
546
  for(j=1;j<=3;j++)
547
    grid->structure[j][4] = 3;
548

549
  grid->noboundaries = 3;
550
  grid->boundint[0] = 1;
551
  grid->boundint[1] = 1;
552
  grid->boundint[2] = 1;
553

554
  grid->boundext[0] = 2;
555
  grid->boundext[1] = 3;
556
  grid->boundext[2] = 0;
557

558
  grid->boundtype[0] = 1;
559
  grid->boundtype[1] = 2;
560
  grid->boundtype[2] = 3;
561

562
  grid->boundsolid[0] = 1;
563
  grid->boundsolid[1] = 1;
564
  grid->boundsolid[2] = 1;
565

566
  grid->mappings = 2;
567
  grid->mappingtype[0] = 1;
568
  grid->mappingtype[1] = 1;
569
  grid->mappingline[0] = 0;
570
  grid->mappingline[1] = 3;
571
  grid->mappinglimits[0] = grid->mappinglimits[1] = 0.5;
572
  grid->mappinglimits[2] = grid->mappinglimits[3] = 0.5;
573
  grid->mappingpoints[0] = grid->mappingpoints[1] = 4;
574
  grid->mappingparams[0] = Rvector(0,grid->mappingpoints[0]);
575
  grid->mappingparams[1] = Rvector(0,grid->mappingpoints[1]);
576
  grid->mappingparams[0][0] = 2.0;
577
  grid->mappingparams[0][1] = 0.0;
578
  grid->mappingparams[0][2] = 6.0;
579
  grid->mappingparams[0][3] = -0.2;
580
  grid->mappingparams[1][0] = 0.0;
581
  grid->mappingparams[1][1] = 0.0;
582
  grid->mappingparams[1][2] = 6.0;
583
  grid->mappingparams[1][3] = 0.3;
584
  
585
  if(info) printf("A simple 2D example mesh was created\n");
586
}
587

588
  
589
static void ExampleGrid3D(struct GridType **grids,int *nogrids,int info) 
590
/* Creates an example grid that might be used to analyze 
591
   a simple acceleration sensor. */
592
{
593
  int i,j;
594
  struct GridType *grid;
595

596
  (*nogrids) = 1;
597
  (*grids) = (struct GridType*)
598
    malloc((size_t) (*nogrids)*sizeof(struct GridType)); 
599

600
  grid = grids[0];
601

602
  InitGrid(grid);
603

604
  grid->dimension = 3;
605
  grid->coordsystem = COORD_CART3;
606
  grid->layered = TRUE;
607

608
  grid->xcells = 4;
609
  grid->ycells = 5;
610
  grid->zcells = 5;
611

612
  grid->wantedelems = 1000;
613
  grid->firstmaterial = 1;
614
  grid->lastmaterial = 3;
615
  grid->autoratio = 1;
616
  grid->coordsystem = COORD_CART3;
617
  grid->numbering = NUMBER_XY;
618

619
  grid->x[0] = 0.0;
620
  grid->x[1] = 0.1;
621
  grid->x[2] = 0.4;
622
  grid->x[3] = 0.5;
623
  grid->x[4] = 0.6;
624

625
  grid->y[0] = 0.0;
626
  grid->y[1] = 0.1;
627
  grid->y[2] = 0.2;
628
  grid->y[3] = 0.4;
629
  grid->y[4] = 0.5;
630
  grid->y[5] = 0.6;
631

632
  grid->z[0] = 0.0;
633
  grid->z[1] = 0.1;
634
  grid->z[2] = 0.2;
635
  grid->z[3] = 0.4;
636
  grid->z[4] = 0.5;
637
  grid->z[5] = 0.6;
638

639
  for(j=1;j<=5;j++)
640
    for(i=1;i<=4;i++)
641
      grid->structure[j][i] = 1;
642

643
  grid->structure[2][2] = 3;
644
  grid->structure[2][3] = 3;
645
  grid->structure[3][3] = 3;
646
  grid->structure[4][3] = 3;
647
  grid->structure[4][2] = 3;
648

649
  grid->structure[3][2] = 2;
650

651
  grid->zlastmaterial[5] = 3;
652
  grid->zlastmaterial[2] = 1;
653
  grid->zlastmaterial[3] = 2;
654
  grid->zlastmaterial[4] = 1;
655
  grid->zlastmaterial[5] = 3;
656

657
  grid->zmaterial[1] = 1;
658
  grid->zmaterial[2] = 2;
659
  grid->zmaterial[3] = 3;
660
  grid->zmaterial[4] = 4;
661
  grid->zmaterial[5] = 5;
662

663
  grid->noboundaries = 4;
664
  grid->boundint[0] = 1;
665
  grid->boundint[1] = 1;
666
  grid->boundint[2] = 1;
667
  grid->boundint[2] = 2;
668

669
  grid->boundext[0] = 0;
670
  grid->boundext[1] = 2;
671
  grid->boundext[2] = 3;
672
  grid->boundext[3] = 3;
673

674
  grid->boundtype[0] = 1;
675
  grid->boundtype[1] = 2;
676
  grid->boundtype[2] = 3;
677
  grid->boundtype[3] = 4;
678

679
  grid->boundsolid[0] = 1;
680
  grid->boundsolid[1] = 1;
681
  grid->boundsolid[2] = 1;
682
  grid->boundsolid[3] = 1;
683

684
  if(info) printf("A simple 3D example mesh was created\n");
685
}
686

687

688
void CreateExampleGrid(int dim,struct GridType **grids,int *nogrids,int info) 
689
{
690
  switch (dim) {
691
    
692
  case 1: 
693
    ExampleGrid1D(grids,nogrids,info);
694
    break;
695

696
  case 2: 
697
    ExampleGrid2D(grids,nogrids,info);
698
    break;
699

700
  case 3: 
701
    ExampleGrid3D(grids,nogrids,info);
702
    break;
703
  }
704
}
705

706

707

708

709
void SetElementDivision(struct GridType *grid,Real relh,int info)
710
/* Given the densities and ratios in each cell finds the 
711
   optimum way to divide the mesh into elements. 
712
   The procedure is the following:
713
   For each subcell set the minimum number of elements
714
   then add one element at a time till the number of 
715
   elements is the desired one. The added element
716
   is always set to the subcell having the sparsest mesh.
717
   Also numbers the cells taking into consideration only 
718
   materials that have indices in interval [firstmat,lastmat].
719
   */
720
{
721
  int i,j,nx,ny,nxmax = 0,nymax = 0;
722
  int sumxelems,sumyelems,sumxyelems;
723
  int wantedelems,wantedelemsx,wantedelemsy;
724
  Real ratio,linearlimit;
725
  Real dxmax = 0,dymax = 0,dx = 0,dy = 0,dxlimit = 0;
726

727
  if(0) printf("SetElementDivision\n");
728

729
  if(grid->dimension == 1) 
730
    grid->numbering = NUMBER_1D;
731

732
  nx = grid->xcells;
733
  ny = grid->ycells;
734
  linearlimit = 0.001;
735

736
  /* Lets number the cells from left to right and up to down. */
737
  grid->nocells = 0;
738
  for(j=1;j<=ny;j++)
739
    for(i=1;i<=nx;i++) {
740
      if (grid->structure[j][i] >= grid->firstmaterial  &&  
741
	  grid->structure[j][i] <= grid->lastmaterial)   
742
        grid->numbered[j][i] = ++grid->nocells; 
743
    }
744
  if(0) printf("The mesh is divided into %d separate subcells.\n",grid->nocells);
745
  
746
  /* Put the linearity flags. */
747
  for(i=1; i<= nx ;i++) {
748
    if (fabs(1.-grid->xexpand[i]) < linearlimit) 
749
      grid->xlinear[i] = TRUE; 
750
    else if (fabs(1.+grid->xexpand[i]) < linearlimit) 
751
      grid->xlinear[i] = TRUE; 
752
    else 
753
      grid->xlinear[i] = FALSE;
754
  }
755

756
  for(i=1; i <= ny ;i++) {
757
    if(fabs(1.-grid->yexpand[i]) < linearlimit) 
758
      grid->ylinear[i] = TRUE;
759
    else if(fabs(1.+grid->yexpand[i]) < linearlimit) 
760
      grid->ylinear[i] = TRUE;
761
    else 
762
      grid->ylinear[i] = FALSE;
763
  }
764

765

766
 /* If there are no materials no elements either.
767
     Otherwise at least grid->minelems elements 
768
     If you want to number elements even if they are 
769
     not there change this initialization to minelems. */
770
  if(grid->autoratio) {
771
    for(i=1;i<=nx;i++)  
772
      grid->xelems[i] = grid->minxelems;
773
    for(i=1;i<=ny;i++)  
774
      grid->yelems[i] = grid->minyelems;
775
  }    
776
  else {
777
    for(i=1;i<=nx;i++)  
778
      if(grid->xelems[i] < grid->minxelems) grid->xelems[i] = grid->minxelems;
779
    for(i=1;i<=ny;i++)  
780
      if(grid->yelems[i] < grid->minyelems) grid->yelems[i] = grid->minyelems;
781
  }
782

783
  sumxelems = 0;
784
  for(i=1;i<=nx;i++) {
785
    if(grid->dimension == 1 && !grid->numbered[1][i]) continue;
786
    sumxelems += grid->xelems[i];
787
  }
788
  sumyelems = 0;
789
  for(i=1;i<=ny;i++)  
790
    sumyelems += grid->yelems[i];
791

792
  if(grid->dimension == 1) {
793
    grid->autoratio = 2;
794
    grid->totxelems = grid->wantedelems;
795
    grid->totyelems = 1;
796
  }
797

798
  /* Allocate elements for both axis separately */
799
  if(grid->autoratio == 2) {
800

801
    wantedelemsx = (int) (1.0*grid->totxelems / relh);
802
    wantedelemsy = (int) (1.0*grid->totyelems / relh);
803

804
    if(sumxelems > wantedelemsx) {
805
      printf("SetElementDivision: %d is too few elements in x-direction\n",grid->totxelems);
806
      wantedelemsx = sumxelems+1;
807
    }      
808
    if(sumyelems > wantedelemsy ) {
809
      printf("SetElementDivision: %d is too few elements in y-direction\n",grid->totyelems);
810
      wantedelemsy = sumyelems+1;
811
    }   
812
    
813
    for(;sumxelems < wantedelemsx;) {
814
      dxmax = 0.0;
815
      for(i=1;i<=nx;i++) {
816
	if(grid->xelems[i] == 0) continue;
817
	
818
	/* Don't put elements to passive subcells */
819
	if(grid->dimension == 1 && !grid->numbered[1][i]) continue; 
820
	
821
	if(grid->xlinear[i] == TRUE || grid->xelems[i]==1) 
822
	  dx = (grid->x[i] - grid->x[i-1])/grid->xelems[i];
823
	else {
824
	  if(grid->xexpand[i] > 0.0) {
825
	    ratio = pow(grid->xexpand[i],1./(grid->xelems[i]-1.));
826
	    dx = (grid->x[i] - grid->x[i-1]) * 
827
	      (1.-ratio) / (1.-pow(ratio,(Real)(grid->xelems[i])));
828
	    if(ratio < 1.)   
829
	      dx *= grid->xexpand[i];
830
	  }
831
	  else {	    
832
	    if(grid->xelems[i]==2) {
833
	      dx = (grid->x[i] - grid->x[i-1])/grid->xelems[i];
834
	    } 
835
	    else if(grid->xelems[i]%2 == 0) {
836
	      ratio = pow(-grid->xexpand[i],1./(grid->xelems[i]/2-1.));
837
	      dx = 0.5 * (grid->x[i] - grid->x[i-1]) * 
838
		(1.-ratio) / (1.-pow(ratio,(Real)(grid->xelems[i]/2)));
839
	    }
840
	    else if(grid->xelems[i]%2 == 1) {
841
	      ratio = pow(-grid->xexpand[i],1./(grid->xelems[i]/2));
842
	      dx = (grid->x[i] - grid->x[i-1]) / 
843
		(2.0*(1.-pow(ratio,(Real)(grid->xelems[i]/2)))/
844
		 (1-ratio) + pow(ratio,(Real)(grid->xelems[i]/2+0.5)));
845
	    }
846
	  }
847
	}
848
	dx *= grid->xdens[i];
849
	if(dx > dxmax) {
850
	  dxmax = dx;
851
	  nxmax = i;
852
	}
853
      }      
854
      grid->xelems[nxmax] += 1;
855
      sumxelems++;
856
    }
857
  
858

859
    if(grid->dimension > 1) {
860
      for(;sumyelems < wantedelemsy;) {
861
	dymax = 0.0;
862
	for(i=1;i<=ny;i++) {
863
	  if(grid->yelems[i] == 0) continue;
864
	  if(grid->ylinear[i] || grid->yelems[i]==1) 
865
	    dy = (grid->y[i] - grid->y[i-1])/grid->yelems[i];
866
	  else {
867
	    ratio = pow(grid->yexpand[i],1./(grid->yelems[i]-1));
868
	    dy = (grid->y[i] - grid->y[i-1]) * 
869
	      (1.-ratio)/(1.-pow(ratio,(Real)(grid->yelems[i])));
870
	    if(ratio < 1.)   
871
	      dy *= grid->yexpand[i];
872
	  }
873
	  dy *= grid->ydens[i];
874
	  if(dy > dymax) {
875
	    dymax = dy;
876
	    nymax = i;
877
	  }
878
	}      
879
	grid->yelems[nymax] += 1;
880
	sumyelems++;
881
      }
882
    }
883
  }   
884

885

886
  /* Both axis dependently */
887
  if(grid->autoratio == 1)  {
888
    
889
    wantedelems = (int) (1.0*grid->wantedelems / (relh*relh));
890

891
    sumxyelems = 0;
892
    for(;;) {
893
      dxmax = 0.0;
894

895
      for(i=1;i<=nx;i++) {
896
	
897
	if(grid->xelems[i] == 0) continue;
898
	if(grid->xlinear[i] || grid->xelems[i]==1) 
899
	  dx = (grid->x[i] - grid->x[i-1])/grid->xelems[i];
900
	else {
901
	  if(grid->xexpand[i] > 0.0) {
902
	    ratio = pow(grid->xexpand[i],1./(grid->xelems[i]-1.));
903
	    dx = (grid->x[i] - grid->x[i-1]) * 
904
	      (1.-ratio) / (1.-pow(ratio,(Real)(grid->xelems[i])));
905
	    if(ratio < 1.)   
906
	      dx *= grid->xexpand[i];
907
	  }
908
	  else if(grid->xelems[i]==2) {
909
	    dx = (grid->x[i] - grid->x[i-1])/grid->xelems[i];
910
	  } 
911
	  else if(grid->xelems[i]%2 == 0) {
912
	    ratio = pow(-grid->xexpand[i],1./(grid->xelems[i]/2-1.));
913
	    dx = 0.5 * (grid->x[i] - grid->x[i-1]) * 
914
	      (1.-ratio) / (1.-pow(ratio,(Real)(grid->xelems[i]/2)));
915
	  }
916
	  else if(grid->xelems[i]%2 == 1) {
917
	    ratio = pow(-grid->xexpand[i],1./(grid->xelems[i]/2));
918
	    dx = (grid->x[i] - grid->x[i-1]) / 
919
	      (2.0*(1.-pow(ratio,(Real)(grid->xelems[i]/2)))/
920
	       (1-ratio) + pow(ratio,(Real)(grid->xelems[i]/2+0.5)));
921
	  }
922
	}
923

924
	dx *= grid->xdens[i];
925
	if(dx > dxmax) {
926
	  dxmax = dx;
927
	  nxmax = i;
928
	}
929
      }      
930

931

932
      dymax = 0.0;
933
      for(i=1;i<=ny;i++) {
934

935
	if(grid->yelems[i] == 0) continue;
936
	if(grid->ylinear[i] || grid->yelems[i]==1) 
937
	  dy = (grid->y[i] - grid->y[i-1])/grid->yelems[i];
938
	else {
939
	  if(grid->yexpand[i] > 0.0) {
940
	    ratio = pow(grid->yexpand[i],1./(grid->yelems[i]-1));
941
	    dy = (grid->y[i] - grid->y[i-1]) * 
942
	      (1.-ratio)/(1.-pow(ratio,(Real)(grid->yelems[i])));
943
	    if(ratio < 1.)   
944
	      dy *= grid->yexpand[i];
945
	  }
946
	  else if(grid->yelems[i]==2) {
947
	    dy = (grid->y[i] - grid->y[i-1])/grid->yelems[i];
948
	  } 
949
	  else if(grid->yelems[i]%2 == 0) {
950
	    ratio = pow(-grid->yexpand[i],1./(grid->yelems[i]/2-1.));
951
	    dy = 0.5 * (grid->y[i] - grid->y[i-1]) * 
952
	      (1.-ratio) / (1.-pow(ratio,(Real)(grid->yelems[i]/2)));
953
	  }
954
	  else if(grid->yelems[i]%2 == 1) {
955
	    ratio = pow(-grid->yexpand[i],1./(grid->yelems[i]/2));
956
	    dy = (grid->y[i] - grid->y[i-1]) / 
957
	      (2.0*(1.-pow(ratio,(Real)(grid->yelems[i]/2)))/
958
	       (1-ratio) + pow(ratio,(Real)(grid->yelems[i]/2+0.5)));
959
	  }
960
	}
961
	dy *= grid->ydens[i];
962
	if(dy > dymax) {
963
	  dymax = dy;
964
	  nymax = i;
965
	}
966
      }      
967
      dymax /= grid->xyratio;
968

969
      if(dxmax > dymax) {
970
	grid->xelems[nxmax] += 1;
971
	sumxelems++;
972
	dxlimit = dxmax;
973
      }
974
      else {
975
	grid->yelems[nymax] += 1;
976
	sumyelems++;
977
	dxlimit = dymax;
978
      }
979
      
980
      sumxyelems = 0;
981
      for(j=1;j<=grid->ycells;j++)
982
	for(i=1;i<=grid->xcells;i++) 
983
	  if(grid->numbered[j][i]) sumxyelems += grid->xelems[i] * grid->yelems[j];
984

985
      if(wantedelems <= sumxyelems) break;
986
    }
987
  }    
988

989

990

991
  if(grid->autoratio == 3 || grid->limitdxverify)  {
992
    
993
    if(grid->autoratio == 3) {
994
      dxlimit = relh * grid->limitdx;
995
      dxmax = dymax = dxlimit;
996
    }
997

998
    for(i=1;i<=nx;i++) {
999

1000
      for(;;) {       
1001
	if(grid->xlinear[i] || grid->xelems[i]==0) 
1002
	  dx = (grid->x[i] - grid->x[i-1])/(grid->xelems[i]+1);
1003
	else {
1004
	  if(grid->xexpand[i] > 0.0) {
1005
	    ratio = pow(grid->xexpand[i],1./grid->xelems[i]);
1006
	    dx = (grid->x[i] - grid->x[i-1]) * 
1007
	      (1.-ratio) / (1.-pow(ratio,(Real)(grid->xelems[i]+1)));
1008
	    if(ratio < 1.)   
1009
	      dx *= grid->xexpand[i];
1010
	  }
1011
	  else if(grid->xelems[i]==1) {
1012
	    dx = (grid->x[i] - grid->x[i-1])/(grid->xelems[i]+1);
1013
	  } 
1014
	  else if((grid->xelems[i]+1)%2 == 0) {
1015
	    ratio = pow(-grid->xexpand[i],1./(grid->xelems[i]/2));
1016
	    dx = 0.5 * (grid->x[i] - grid->x[i-1]) * 
1017
	      (1.-ratio) / (1.-pow(ratio,(Real)(grid->xelems[i]/2+1)));
1018
	  }
1019
	  else if((grid->xelems[i]+1)%2 == 1) {
1020
	    ratio = pow(-grid->xexpand[i],1./((grid->xelems[i]+1)/2));
1021
	    dx = (grid->x[i] - grid->x[i-1]) / 
1022
	      (2.0*(1.-pow(ratio,(Real)((grid->xelems[i]+1)/2)))/
1023
	       (1-ratio) + pow(ratio,(Real)((grid->xelems[i]+1)/2+0.5)));
1024
	  }
1025
	}
1026
	dx *= grid->xdens[i];
1027
	
1028
	/* choose the best fit for desired density */
1029
	if(fabs(dx-dxlimit) < fabs( dx*(1+1.0/grid->xelems[i]) -dxlimit) )
1030
	  grid->xelems[i] += 1;
1031
	else
1032
	  break;
1033
      }
1034
      sumxelems += grid->xelems[i];
1035
    }      
1036

1037

1038
    for(i=1;i<=ny;i++) {
1039
      
1040
      for(;;) {
1041
	if(grid->ylinear[i] || grid->yelems[i]==0) 
1042
	  dy = (grid->y[i] - grid->y[i-1])/(grid->yelems[i]+1);
1043
	else {
1044
	  if(grid->yexpand[i] > 0.0) {
1045
	    ratio = pow(grid->yexpand[i],1./grid->yelems[i]);
1046
	    dy = (grid->y[i] - grid->y[i-1]) * 
1047
	      (1.-ratio) / (1.-pow(ratio,(Real)(grid->yelems[i]+1)));
1048
	    if(ratio < 1.)   
1049
	      dy *= grid->yexpand[i];
1050
	  }
1051
	  else if(grid->yelems[i]==1) {
1052
	    dy = (grid->y[i] - grid->y[i-1])/(grid->yelems[i]+1);
1053
	  } 
1054
	  else if((grid->yelems[i]+1)%2 == 0) {
1055
	    ratio = pow(-grid->yexpand[i],1./(grid->yelems[i]/2));
1056
	    dy = 0.5 * (grid->y[i] - grid->y[i-1]) * 
1057
	      (1.-ratio) / (1.-pow(ratio,(Real)(grid->yelems[i]/2+1)));
1058
	  }
1059
	  else if((grid->yelems[i]+1)%2 == 1) {
1060
	    ratio = pow(-grid->yexpand[i],1./((grid->yelems[i]+1)/2));
1061
	    dy = (grid->y[i] - grid->y[i-1]) / 
1062
	      (2.0*(1.-pow(ratio,(Real)((grid->yelems[i]+1)/2)))/
1063
	       (1-ratio) + pow(ratio,(Real)((grid->yelems[i]+1)/2+0.5)));
1064
	  }
1065
	}
1066
	
1067
	dy *= grid->ydens[i] / grid->xyratio;
1068
	/* choose the best fit for desired density */
1069
	if(fabs(dy-dxlimit) < fabs( dy*(1+1.0/grid->yelems[i]) -dxlimit) )
1070
	  grid->yelems[i] += 1;
1071
	else
1072
	  break; 
1073
      }
1074
      sumyelems += grid->yelems[i];
1075
    }      
1076
    
1077
    sumxyelems = 0;
1078
    for(j=1;j<=grid->ycells;j++)
1079
      for(i=1;i<=grid->xcells;i++) 
1080
	if(grid->numbered[j][i]) sumxyelems += grid->xelems[i] * grid->yelems[j];
1081
  }
1082

1083

1084
  /* Put the linearity flags if there is only one division. */
1085
  for(i=1; i<= nx ;i++) 
1086
    if (grid->xelems[i] == 1) 
1087
      grid->xlinear[i] = TRUE; 
1088
  for(i=1; i <= ny ;i++) 
1089
    if(grid->yelems[i] == 1) 
1090
      grid->ylinear[i] = TRUE;
1091
 
1092

1093
  /* Set the size of the first element within each subcell */
1094
  for(i=1;i<=nx;i++) {
1095

1096
    if(grid->xlinear[i] == TRUE) 
1097
      grid->dx[i] = (grid->x[i] - grid->x[i-1])/grid->xelems[i];
1098
    else {
1099
      if(grid->xexpand[i] > 0.0) {
1100
	ratio = pow(grid->xexpand[i],1./(grid->xelems[i]-1.));
1101
	grid->xratios[i] = ratio;
1102
	grid->dx[i] = (grid->x[i] - grid->x[i-1]) * 
1103
	  (1.-ratio) / (1.-pow(ratio,(Real)(grid->xelems[i])));
1104
      }
1105
      else if(grid->xelems[i] == 2) {
1106
	grid->dx[i] = (grid->x[i] - grid->x[i-1])/grid->xelems[i];
1107
      }
1108
      else if(grid->xelems[i]%2 == 0) {
1109
	ratio = pow(-grid->xexpand[i],1./(grid->xelems[i]/2-1.));
1110
	grid->xratios[i] = ratio;
1111
	grid->dx[i] = 0.5 * (grid->x[i] - grid->x[i-1]) * 
1112
	  (1.-ratio) / (1.-pow(ratio,(Real)(grid->xelems[i]/2)));
1113
      }
1114
      else if(grid->xelems[i]%2 == 1) {
1115
	ratio = pow(-grid->xexpand[i],1./(grid->xelems[i]/2));
1116
	grid->xratios[i] = ratio;
1117
	grid->dx[i] = (grid->x[i] - grid->x[i-1]) / 
1118
	  (2.0*(1.-pow(ratio,(Real)(grid->xelems[i]/2)))/
1119
	   (1-ratio) + pow(ratio,(Real)(grid->xelems[i]/2+0.5)));
1120
      }
1121
    }	
1122
  }
1123

1124
  for(i=1;i<=ny;i++) {
1125

1126
    if(grid->ylinear[i] == TRUE) 
1127
      grid->dy[i] = (grid->y[i] - grid->y[i-1])/grid->yelems[i];
1128
    else {
1129
      if(grid->yexpand[i] > 0.0) {
1130
      ratio = pow(grid->yexpand[i],1./(grid->yelems[i]-1.));
1131
      grid->yratios[i] = ratio;
1132
      grid->dy[i] = (grid->y[i] - grid->y[i-1]) * 
1133
	(1.-ratio)/(1.-pow(ratio,(Real)(grid->yelems[i])));
1134
      }
1135
      else if(grid->yelems[i] == 2) {
1136
	grid->dy[i] = (grid->y[i] - grid->y[i-1])/grid->yelems[i];
1137
      }
1138
      else if(grid->yelems[i]%2 == 0) {
1139
	ratio = pow(-grid->yexpand[i],1./(grid->yelems[i]/2-1.));
1140
	grid->yratios[i] = ratio;
1141
	grid->dy[i] = 0.5 * (grid->y[i] - grid->y[i-1]) * 
1142
	  (1.-ratio) / (1.-pow(ratio,(Real)(grid->yelems[i]/2)));
1143
      }
1144
      else if(grid->yelems[i]%2 == 1) {
1145
	ratio = pow(-grid->yexpand[i],1./(grid->yelems[i]/2));
1146
	grid->yratios[i] = ratio;
1147
	grid->dy[i] = (grid->y[i] - grid->y[i-1]) / 
1148
	  (2.0*(1.-pow(ratio,(Real)(grid->yelems[i]/2)))/
1149
	   (1-ratio) + pow(ratio,(Real)(grid->yelems[i]/2+0.5)));
1150
      }
1151
    }
1152
  }
1153

1154
  /* Calculate the total number of elements */
1155
  sumxyelems = 0;
1156
  for(j=1;j<=grid->ycells;j++)
1157
    for(i=1;i<=grid->xcells;i++) 
1158
      if(grid->numbered[j][i]) sumxyelems += grid->xelems[i] * grid->yelems[j];
1159

1160
  grid->noelements = sumxyelems;
1161

1162
  if(0) printf("Created a total of %d elements\n",sumxyelems);
1163

1164
  grid->dx0 = dxmax;
1165
  grid->dy0 = dymax;
1166
}
1167

1168

1169

1170

1171
void SetCellData(struct GridType *grid,struct CellType *cell,int info)
1172
/* Sets the data that can directly be derived from type GridType. 
1173
 */
1174
{
1175
  int i,j,cnew=0;
1176

1177
  for(j=1;j<= grid->ycells ;j++)                   /* cells direction up    */
1178
    for(i=1;i<= grid->xcells; i++)                 /* cells direction right */      
1179
      if( (cnew = grid->numbered[j][i]) ) {          /* if cell is occupied   */
1180

1181
        /* Initialize the numbering to zero  */
1182
        cell[cnew].left1st = 0;
1183
        cell[cnew].left2nd = 0;
1184
        cell[cnew].leftlast = 0;
1185
        cell[cnew].levelwidth = 0;
1186
        cell[cnew].levelwidthcenter = 0;
1187
	cell[cnew].leftcenter = 0;
1188
	cell[cnew].elemwidth = 0;
1189
	cell[cnew].elem1st = 0;
1190
	
1191
	/* Set the element type */
1192
	cell[cnew].nonodes = grid->nonodes;
1193
	cell[cnew].numbering = grid->numbering;
1194
	
1195
	/* Set the cell coordinates */
1196
        cell[cnew].xcorner[TOPLEFT] = cell[cnew].xcorner[BOTLEFT] = grid->x[i-1];
1197
	cell[cnew].xcorner[TOPRIGHT] = cell[cnew].xcorner[BOTRIGHT] = grid->x[i];
1198
        cell[cnew].ycorner[BOTRIGHT] = cell[cnew].ycorner[BOTLEFT] = grid->y[j-1];
1199
        cell[cnew].ycorner[TOPRIGHT] = cell[cnew].ycorner[TOPLEFT] = grid->y[j];
1200
	
1201
	/* Set the cell width in both directions */	
1202
	cell[cnew].xwidth = grid->x[i] - grid->x[i-1];
1203
	cell[cnew].ywidth = grid->y[j] - grid->y[j-1];
1204
	
1205
	/* Set the number of elements in a subcell */
1206
	cell[cnew].xelem = grid->xelems[i];
1207
	cell[cnew].yelem = grid->yelems[j];
1208
	
1209
	/* Set the the ratio of elements when going left to right and down to up */
1210
	cell[cnew].xratio = grid->xratios[i];
1211
	cell[cnew].yratio = grid->yratios[j];
1212
	cell[cnew].xlinear = grid->xlinear[i];
1213
	cell[cnew].ylinear = grid->ylinear[j];
1214
	cell[cnew].xlinear = grid->xlinear[i];
1215
	cell[cnew].ylinear = grid->ylinear[j];
1216
	if(grid->xexpand[i] < 0.0 && grid->xlinear[i] == FALSE) {
1217
	  if(grid->xelems[i] == 2) cell[cnew].xlinear = 1;
1218
	  else cell[cnew].xlinear = 2;
1219
	}
1220
	if(grid->yexpand[j] < 0.0 && grid->ylinear[j] == FALSE) {
1221
	  if(grid->yelems[j] == 2) cell[cnew].ylinear = 1;
1222
	  else cell[cnew].ylinear = 2;
1223
	}
1224
	
1225
	/* Set the length of the first element in both directions  */
1226
	cell[cnew].dx1 = grid->dx[i];
1227
	cell[cnew].dy1 = grid->dy[j];
1228
	
1229
	/* Set the material in question */
1230
	cell[cnew].material = grid->structure[j][i];
1231
	
1232
	/* Set the boundary types for sides and corners */
1233
	cell[cnew].boundary[LEFT] = grid->structure[j][i-1];
1234
	cell[cnew].boundary[DOWN] = grid->structure[j-1][i];
1235
	cell[cnew].boundary[RIGHT]= grid->structure[j][i+1];  
1236
	cell[cnew].boundary[UP]   = grid->structure[j+1][i];
1237
	cell[cnew].boundary[4] = grid->structure[j-1][i-1];
1238
	cell[cnew].boundary[5] = grid->structure[j-1][i+1];
1239
	cell[cnew].boundary[6] = grid->structure[j+1][i+1];
1240
	cell[cnew].boundary[7] = grid->structure[j+1][i-1];
1241

1242
	/* Set the neighbour cell indices for sides and corners */ 
1243
	cell[cnew].neighbour[LEFT] = grid->numbered[j][i-1];
1244
	cell[cnew].neighbour[DOWN] = grid->numbered[j-1][i];
1245
	cell[cnew].neighbour[RIGHT]= grid->numbered[j][i+1];  
1246
	cell[cnew].neighbour[UP]   = grid->numbered[j+1][i];
1247
	cell[cnew].neighbour[4] = grid->numbered[j-1][i-1];
1248
	cell[cnew].neighbour[5] = grid->numbered[j-1][i+1];
1249
	cell[cnew].neighbour[6] = grid->numbered[j+1][i+1];
1250
	cell[cnew].neighbour[7] = grid->numbered[j+1][i-1];
1251
      }
1252

1253
  if(info) printf("%d cells were created.\n",grid->nocells);
1254
}
1255

1256

1257

1258
static int SetCellKnots(struct GridType *grid, struct CellType *cell,int info)
1259
/* Uses given mesh to number the knots present in the cells.
1260
   The knots are numbered independently of the cells from left 
1261
   to right and up to down. Only the numbers of four knots at the 
1262
   left side of each cell are saved for later use. The number of 
1263
   each knot can later be easily recalled using simple functions. 
1264
   The subroutine was initially written for ordering the knots
1265
   from left to right and down to up. However, this does not always
1266
   produce reasonably small bandwidth and therefore a numbering 
1267
   scheme for the other order was later added. Therefore the algorithm 
1268
   may not be that clear for this other scheme.
1269
   Numbers the knots in rectangular elements. There may be 
1270
   4, 5, 8, 9, 12 or 16 nodes in each elements.
1271
   */
1272
{
1273
  int i,j,level,center;
1274
  int degree,centernodes,sidenodes,nonodes;
1275
  int cnew=0,cup=0,cleft=0,cleftup=0;
1276
  int elemno,knotno;
1277
  int maxwidth,width,numbering;
1278
  int xcells,ycells,*yelems=NULL,*xelems=NULL;
1279

1280
  numbering = grid->numbering;
1281
  nonodes = grid->nonodes;
1282
  knotno  = 0;
1283
  elemno  = 0;
1284

1285
  if(numbering == NUMBER_XY) {
1286
    xcells = grid->xcells;
1287
    ycells = grid->ycells;
1288
    xelems = grid->xelems;
1289
    yelems = grid->yelems;
1290
  }
1291
  else if(numbering == NUMBER_YX) {
1292
    xcells = grid->ycells;
1293
    ycells = grid->xcells;
1294
    xelems = grid->yelems;
1295
    yelems = grid->xelems;
1296
  }
1297
  else { 
1298
    printf("No %d numbering scheme exists!\n",numbering); 
1299
    return(1);
1300
  }
1301

1302
  switch (nonodes) {
1303
  case 4:
1304
    centernodes = FALSE;
1305
    sidenodes = FALSE;
1306
    degree = 1;
1307
    break;
1308
  case 5:
1309
    centernodes = TRUE;
1310
    sidenodes = FALSE;
1311
    degree = 2;
1312
    break;
1313
  case 8:
1314
    centernodes = FALSE;
1315
    sidenodes = TRUE;
1316
    degree = 2;
1317
    break;
1318
  case 9:
1319
    centernodes = TRUE;
1320
    sidenodes = TRUE;
1321
    degree = 2;
1322
    break;
1323
  case 12:
1324
    centernodes = FALSE;
1325
    sidenodes = TRUE;
1326
    degree = 3;
1327
    break;
1328
  case 16:
1329
    centernodes = TRUE;
1330
    sidenodes = TRUE;
1331
    degree = 3;
1332
    break;
1333
  default:
1334
    printf("SetCellKnots: No numbering scheme for %d-node elements.\n",
1335
	   grid->nonodes);
1336
    return(2);
1337
  }
1338

1339
  for(j=1;j<= ycells ;j++)                     /* cells direction up     */
1340
    for(level=0; level <= yelems[j] ;level++)  /* lines inside cells     */
1341
      for(center=0; center < degree; center++) 
1342
      for(i=1;i<= xcells; i++)                 /* cells direction right  */
1343
	{        
1344
	  if(numbering == NUMBER_XY) {	  
1345
	    cnew = grid->numbered[j][i];
1346
	    cleft= grid->numbered[j][i-1];
1347
	    cup  = grid->numbered[j+1][i];
1348
	    cleftup = grid->numbered[j+1][i-1];
1349
	    if(cnew) {
1350
	      cell[cnew].xind = i;
1351
	      cell[cnew].yind = j;	      
1352
	    }
1353
	  }	  
1354
	  else if(numbering == NUMBER_YX) {
1355
	    cnew = grid->numbered[i][j];
1356
	    cleft= grid->numbered[i-1][j];
1357
	    cup  = grid->numbered[i][j+1];
1358
	    cleftup = grid->numbered[i-1][j+1];
1359
	    if(cnew) {
1360
	      cell[cnew].xind = j;
1361
	      cell[cnew].yind = i;	      
1362
	    }
1363
	  }
1364
	  
1365
	  if(center == 0) {
1366
	    /* the first line of an occupied cell is unnumbered, 
1367
	       this happens only in the bottom */
1368
	    if (cnew && level==0 && !cell[cnew].left1st) { 
1369
	      if(!cleft) knotno++;
1370
	      cell[cnew].left1st = knotno;
1371
	      if(sidenodes) 
1372
		knotno += degree * xelems[i];
1373
	      else
1374
		knotno += xelems[i];
1375
	    } /* level=0 */
1376
	    
1377
	    /* the n:th line of an occupied cell */
1378
	    else if (cnew && level > 0 && level < yelems[j]) {   
1379
	      if(!cleft) knotno++;
1380
	      if(level==1) 
1381
		cell[cnew].left2nd = knotno;
1382
	      if(level==2) 
1383
		cell[cnew].levelwidth = knotno - cell[cnew].left2nd; 	      
1384
	      if(sidenodes) 
1385
		knotno += degree * xelems[i];
1386
	      else 
1387
		knotno += xelems[i];
1388
	    } /* 1 <= level < n */	  
1389
	  
1390
	    /* last line of this cell or first line of upper cell */
1391
	    else if ((cnew || cup) && level == yelems[j]) { 
1392
	      if(!cleft && !cleftup) 
1393
		knotno++;
1394
	      if(cnew)
1395
		cell[cnew].leftlast = knotno;
1396
	      if(level==2) 
1397
		cell[cnew].levelwidth = knotno - cell[cnew].left2nd; 	      
1398
	      if(yelems[j] == 1) 
1399
		cell[cnew].left2nd = cell[cnew].leftlast;
1400
	      if(cup)
1401
		cell[cup].left1st = knotno;
1402
	      if(sidenodes) 
1403
		knotno += degree * xelems[i];
1404
	      else 
1405
		knotno += xelems[i];
1406
	    } /* level=n */
1407

1408
	    /* Number the elements */
1409
	    if(cnew && level > 0) {
1410
	      if(level==1)  
1411
		cell[cnew].elem1st = elemno+1;	      
1412
	      if(level==2) 
1413
		cell[cnew].elemwidth  = (elemno+1) - cell[cnew].elem1st;
1414
	      elemno += xelems[i];	    
1415
	    }
1416
	  }
1417

1418
	  if(center && level < yelems[j]) {
1419
	    if (cnew) { 
1420
	      if(!cleft && sidenodes) 
1421
		knotno++;
1422
	      if(level==0 && center==1) 
1423
		cell[cnew].leftcenter = knotno;
1424
	      if(level==0 && center==2)
1425
		cell[cnew].left2center = knotno;		
1426
	      if(level==1 && center==1)
1427
		cell[cnew].levelwidthcenter = knotno - cell[cnew].leftcenter;
1428
	      if(centernodes && sidenodes)
1429
		knotno += degree * xelems[i];
1430
	      else 
1431
		knotno += xelems[i];
1432
	    } 
1433
	  }
1434
	
1435
	} /* x-cell loop */
1436
  
1437

1438
  maxwidth = 0;
1439
  for(j=1;j<= ycells ;j++) 
1440
    for(i=1;i<= xcells; i++) {
1441
      if(numbering == NUMBER_XY)
1442
	cnew = grid->numbered[j][i];
1443
      else if(numbering == NUMBER_YX)
1444
	cnew = grid->numbered[i][j];      
1445
      if(cnew) {
1446
	width = cell[cnew].left2nd - cell[cnew].left1st;
1447
	if(width > maxwidth) 
1448
	  maxwidth = width;
1449
	width = cell[cnew].levelwidth;
1450
	if(width > maxwidth) 
1451
	  maxwidth = width;
1452
	width = cell[cnew].leftlast-cell[cnew].left2nd
1453
	  -(yelems[j]-2)*cell[cnew].levelwidth;
1454
	if(width > maxwidth) 
1455
	  maxwidth = width;
1456
      }  
1457
    }
1458
  maxwidth += 2;
1459

1460
  grid->maxwidth = maxwidth;
1461
  grid->noknots = knotno;
1462
  grid->noelements = elemno;
1463

1464
  if(info) {
1465
    printf("Numbered %d knots in %d %d-node elements.\n",knotno,elemno,nonodes);
1466
    if(numbering == NUMBER_XY) 
1467
      printf("Numbering order was <x><y> and max levelwidth %d.\n",
1468
	     maxwidth);
1469
    else if(numbering == NUMBER_YX) 
1470
      printf("Numbering order was <y><x> and max levelwidth %d.\n",
1471
	     maxwidth);
1472
  }
1473

1474
  return(0);
1475
}
1476

1477

1478

1479
static int SetCellKnots1D(struct GridType *grid, struct CellType *cell,int info)
1480
{
1481
  int i;
1482
  int degree,nonodes;
1483
  int cnew,cleft;
1484
  int elemno,knotno;
1485
  int maxwidth;
1486
  int xcells,*xelems;
1487

1488
  nonodes = grid->nonodes;
1489
  knotno  = 0;
1490
  elemno  = 0;
1491

1492
  xcells = grid->xcells;
1493
  xelems = grid->xelems;
1494

1495
  switch (nonodes) {
1496
  case 2:
1497
    degree = 1;
1498
    break;
1499
  case 3:
1500
    degree = 2;
1501
    break;
1502
  case 4:
1503
    degree = 3;
1504
    break;
1505

1506
  default:
1507
    printf("SetCellKnots1D: No numbering scheme for %d-node elements.\n",
1508
	   grid->nonodes);
1509
    return(2);
1510
  }
1511

1512
  for(i=1;i<= xcells; i++) {        
1513
    
1514
    cnew = grid->numbered[1][i];
1515
    cleft= grid->numbered[1][i-1];
1516
    if(cnew) cell[cnew].xind = i;
1517
    
1518
    if (cnew) { 
1519
      /* Number the nodes */
1520
      if(!cleft) knotno++;
1521
      cell[cnew].left1st = knotno;
1522
      knotno += degree * xelems[i];
1523
      
1524
      /* Number the elements */
1525
      cell[cnew].elem1st = elemno+1;	      
1526
      elemno += xelems[i];	    
1527
    }
1528
  } /* x-cell loop */
1529
  
1530
  maxwidth = degree;
1531
  
1532
  grid->maxwidth = maxwidth;
1533
  grid->noknots = knotno;
1534
  grid->noelements = elemno;
1535
  
1536
  if(info) {
1537
    printf("Numbered %d knots in %d %d-node elements.\n",knotno,elemno,nonodes);
1538
  }
1539

1540
  return(0);
1541
}
1542

1543

1544

1545
void CreateCells(struct GridType *grid,struct CellType **cell,int info)
1546
{
1547
  (*cell) = (struct CellType*)
1548
    malloc((size_t) (grid->nocells+1)*sizeof(struct CellType)); 
1549

1550
  SetCellData(grid,*cell,info);
1551

1552
  if(grid->dimension == 1) 
1553
    SetCellKnots1D(grid,*cell,info);
1554
  else
1555
    SetCellKnots(grid,*cell,info);
1556
}
1557

1558

1559
void DestroyCells(struct CellType **cell)
1560
{
1561
  free(cell);
1562
}
1563

1564

1565

1566
int GetKnotIndex(struct CellType *cell,int i,int j)
1567
/* Given the cell and knot indices gives the corresponding 
1568
   global knot index. The indices are expected to be in the 
1569
   range [0..n] and [0..m]. Requires only the structure CellType. 
1570
   */
1571
{
1572
  int ind=0,aid=0,maxj=0;
1573

1574
  if(cell->numbering == NUMBER_1D) {
1575
    ind = cell->left1st;
1576
    if(cell->nonodes == 2) 
1577
      ind += i;
1578
    else if(cell->nonodes == 3) 
1579
      ind += 2*i;
1580
    else if(cell->nonodes == 4) 
1581
      ind += 3*i;
1582
    return(ind);
1583
  }
1584

1585
  if(cell->numbering == NUMBER_XY) 
1586
    maxj = cell->yelem;
1587
  else if(cell->numbering == NUMBER_YX) {
1588
    aid = j; j = i; i = aid;
1589
    maxj = cell->xelem;
1590
  }
1591
  else {
1592
    maxj = 0;
1593
    bigerror("GetKnotIndex: Unknown numbering scheme!");
1594
  }
1595

1596
  if(j == 0) 
1597
    ind = cell->left1st;
1598
  else if (j == maxj)
1599
    ind = cell->leftlast;
1600
  else
1601
    ind = cell->left2nd + (j-1) * cell->levelwidth;
1602
  
1603
  if(cell->nonodes == 4) 
1604
    ind += i;
1605
  else if(cell->nonodes == 5)
1606
    ind += i;
1607
  else if(cell->nonodes == 8 || cell->nonodes == 9) 
1608
    ind += 2*i;
1609
  else if(cell->nonodes == 12 || cell->nonodes == 16) 
1610
    ind += 3*i;
1611

1612
  return(ind);
1613
}
1614

1615

1616
int GetKnotCoordinate(struct CellType *cell,int i,int j,Real *x,Real *y)
1617
/* Given the cell and element indices inside the cell gives the 
1618
   corresponding global knot numbers. The indices are expected 
1619
   to be in the range [0..n] and [0..m]. Requires only the 
1620
   structure CellType.   
1621
   */
1622
{
1623
  int ind;
1624

1625
  if(cell->xlinear == 1)
1626
    (*x) = cell->xcorner[BOTLEFT] + cell->dx1 * i;
1627
  else if(cell->xlinear == 0)
1628
    (*x) = cell->xcorner[BOTLEFT] + cell->dx1 * 
1629
    (1.- pow(cell->xratio,(Real)(i))) / (1.-cell->xratio);
1630
  else if(cell->xlinear == 2) {
1631
    if(i<=cell->xelem/2) {
1632
      (*x) = cell->xcorner[BOTLEFT] + cell->dx1 * 
1633
	(1.- pow(cell->xratio,(Real)(i))) / (1.-cell->xratio);
1634
    }
1635
    else {
1636
      (*x) = cell->xcorner[BOTRIGHT] - cell->dx1 * 
1637
	(1.- pow(cell->xratio,(Real)(cell->xelem-i))) / (1.-cell->xratio);
1638
    }
1639
  }
1640

1641
  if(cell->ylinear == 1)
1642
    (*y) = cell->ycorner[BOTLEFT] + cell->dy1 * j;
1643
  else if(cell->ylinear == 0) 
1644
    (*y) = cell->ycorner[BOTLEFT] + cell->dy1 * 
1645
           (1.- pow(cell->yratio,(Real)(j))) / (1.-cell->yratio);
1646
  else if(cell->ylinear == 2) {
1647
    if(j<=cell->yelem/2) {
1648
      (*y) = cell->ycorner[BOTLEFT] + cell->dy1 * 
1649
	(1.- pow(cell->yratio,(Real)(j))) / (1.-cell->yratio);
1650
    }
1651
    else {
1652
      (*y) = cell->ycorner[TOPLEFT] - cell->dy1 * 
1653
	(1.- pow(cell->yratio,(Real)(cell->yelem-j))) / (1.-cell->yratio);
1654
    }
1655
  }
1656

1657
  ind = GetKnotIndex(cell,i,j);
1658

1659
  return(ind);
1660
}
1661

1662

1663

1664
int GetElementIndices(struct CellType *cell,int i,int j,int *ind)
1665
/* For given element gives the coordinates and index for each knot.
1666
   The indices i and j are expected to range from [1..n] and [1..m]. 
1667
   requires only the structure CellType.
1668
   */
1669
{
1670
  int nonodes,numbering,elemind=0;
1671
  
1672
  nonodes = cell->nonodes;
1673
  numbering = cell->numbering;
1674

1675
  if(numbering != NUMBER_1D) {
1676
    ind[TOPRIGHT] = GetKnotIndex(cell,i,j);
1677
    ind[BOTRIGHT] = GetKnotIndex(cell,i,j-1);
1678
    ind[TOPLEFT]  = GetKnotIndex(cell,i-1,j);
1679
    ind[BOTLEFT]  = GetKnotIndex(cell,i-1,j-1);
1680
  }    
1681

1682
  if(numbering == NUMBER_XY) {
1683
    elemind = cell->elem1st+(i-1) + (j-1)*cell->elemwidth;
1684
    if(nonodes == 4) return(elemind);
1685
        
1686
    if(nonodes == 5) {
1687
      ind[4] = cell->leftcenter + (j-1) * cell->levelwidthcenter + i;
1688
    }
1689
    else if(nonodes == 8) {
1690
      ind[4] = ind[0]+1;
1691
      ind[6] = ind[3]+1;
1692
      ind[5] = cell->leftcenter + (j-1) * cell->levelwidthcenter + i;
1693
      ind[7] = ind[5]-1;
1694
    }
1695
    else if(nonodes == 9) {
1696
      ind[4] = ind[0]+1;
1697
      ind[6] = ind[3]+1;
1698
      ind[5] = cell->leftcenter + (j-1) * cell->levelwidthcenter + 2*i;
1699
      ind[7] = ind[5]-2;
1700
      ind[8] = ind[5]-1;
1701
    }
1702
    else if(nonodes == 12) {
1703
      ind[4]  = ind[0]+1;
1704
      ind[5]  = ind[0]+2;
1705
      ind[9]  = ind[3]+1;
1706
      ind[8]  = ind[3]+2;
1707
      ind[6]  = cell->leftcenter + (j-1) * cell->levelwidthcenter + i;
1708
      ind[11] = ind[6]-1;
1709
      ind[7]  = cell->left2center + (j-1) * cell->levelwidthcenter + i;
1710
      ind[10] = ind[7]-1;
1711
    }
1712
    else if(nonodes == 16) {
1713
      ind[4]  = ind[0]+1;
1714
      ind[5]  = ind[0]+2;
1715
      ind[9]  = ind[3]+1;
1716
      ind[8]  = ind[3]+2;
1717
      ind[6]  = cell->leftcenter + (j-1) * cell->levelwidthcenter + 3*i;
1718
      ind[11] = ind[6]-3;
1719
      ind[7]  = cell->left2center + (j-1) * cell->levelwidthcenter + 3*i;
1720
      ind[10] = ind[7]-3;
1721
      ind[13] = ind[6]-1;
1722
      ind[12] = ind[6]-2;
1723
      ind[14] = ind[7]-1;
1724
      ind[15] = ind[7]-2;
1725
    }
1726
    else 
1727
      printf("GetElementIndices: not implemented for %d nodes.\n",nonodes);    
1728
  }
1729

1730
  else if(numbering == NUMBER_YX) { 
1731
    elemind = cell->elem1st+(j-1) + (i-1)*cell->elemwidth;
1732
    
1733
    if(nonodes == 4) return(elemind);
1734

1735
    if(nonodes == 5) {
1736
      ind[4] = cell->leftcenter + (i-1) * cell->levelwidthcenter + j;
1737
    }
1738
    else if (nonodes==8) {
1739
      ind[7] = ind[0]+1;
1740
      ind[5] = ind[1]+1;
1741
      ind[6] = cell->leftcenter + (i-1) * cell->levelwidthcenter + j;
1742
      ind[4] = ind[6]-1;
1743
    }
1744
    else if(nonodes == 9) {
1745
      ind[7] = ind[0]+1;
1746
      ind[5] = ind[1]+1;
1747
      ind[6] = cell->leftcenter + (i-1) * cell->levelwidthcenter + 2*j;
1748
      ind[4] = ind[6]-2;
1749
      ind[8] = ind[6]-1;
1750
    }
1751
    else if(nonodes == 12) {
1752
      ind[11]= ind[0]+1;
1753
      ind[10]= ind[0]+2;
1754
      ind[6] = ind[1]+1;
1755
      ind[7] = ind[1]+2;
1756
      ind[9] = cell->leftcenter + (i-1) * cell->levelwidthcenter + j;
1757
      ind[4] = ind[9]-1;
1758
      ind[8] = cell->left2center + (i-1) * cell->levelwidthcenter + j;
1759
      ind[5] = ind[8]-1;
1760
    }
1761
    else if(nonodes == 16) {
1762
      ind[11]= ind[0]+1;
1763
      ind[10]= ind[0]+2;
1764
      ind[6] = ind[1]+1;
1765
      ind[7] = ind[1]+2;
1766
      ind[9] = cell->leftcenter + (i-1) * cell->levelwidthcenter + 3*j;
1767
      ind[4] = ind[9]-3;
1768
      ind[15]= ind[9]-1;
1769
      ind[12]= ind[9]-2;
1770
      ind[8] = cell->left2center + (i-1) * cell->levelwidthcenter + 3*j;
1771
      ind[5] = ind[8]-3;
1772
      ind[14]= ind[8]-1;
1773
      ind[13]= ind[8]-2;
1774
    }
1775
    else 
1776
      printf("GetElementIndices: not implemented for %d nodes.\n",nonodes);
1777
  }
1778

1779
  else if(numbering == NUMBER_1D) {
1780
    elemind = cell->elem1st+(i-1);
1781
    ind[0] = GetKnotIndex(cell,i-1,1);
1782
    if(nonodes == 2) {
1783
      ind[1] = ind[0] + 1;
1784
    }
1785
    else if(nonodes == 3) {
1786
      ind[2] = ind[0] + 1;
1787
      ind[1] = ind[0] + 2;
1788
    }
1789
    else if(nonodes == 4) {
1790
      ind[2] = ind[0] + 1;
1791
      ind[3] = ind[0] + 2;
1792
      ind[1] = ind[0] + 3;
1793
    }
1794
  }  
1795
  return(elemind);
1796
}
1797

1798

1799
int GetElementIndex(struct CellType *cell,int i,int j)
1800
/* For given element gives the element index. 
1801
   The indices i and j are expected to range from [1..n] and [1..m]. 
1802
   requires only the structure CellType.
1803
   */
1804
{
1805
  int elemind=0;
1806
 
1807
  if(cell->numbering == NUMBER_XY) 
1808
    elemind = cell->elem1st+(i-1) + (j-1)*cell->elemwidth;
1809
  else if(cell->numbering == NUMBER_YX) 
1810
    elemind = cell->elem1st+(j-1) + (i-1)*cell->elemwidth;
1811

1812
  return(elemind);
1813
}
1814

1815

1816

1817
int GetElementCoordinates(struct CellType *cell,int i,int j,
1818
			  Real *globalcoord,int *ind)
1819
/* For given element gives the coordinates and index for each knot.
1820
   The indices i and j are expected to range from [1..n] and [1..m]. 
1821
   requires only the structure CellType
1822
   Note that this subroutine assumes that the elements are
1823
   rectangular.
1824
   */
1825
{
1826
  int k,nonodes,numbering,elemind=0;
1827
  Real xrat,yrat;
1828

1829
  k = nonodes = cell->nonodes;
1830
  numbering = cell->numbering;
1831

1832
  if(numbering == NUMBER_1D) {
1833
    elemind = cell->elem1st+(i-1);    
1834
    ind[0] = GetKnotCoordinate(cell,i-1,j-1,&globalcoord[0],
1835
			       &globalcoord[nonodes]);
1836
    ind[1] = GetKnotCoordinate(cell,i,j-1,&globalcoord[1],
1837
			       &globalcoord[nonodes+1]);    
1838
    
1839
    if(nonodes == 3) {
1840
      globalcoord[2] = (globalcoord[0]+globalcoord[1])/2.0;
1841
      globalcoord[5] = (globalcoord[3]+globalcoord[4])/2.0;
1842
      ind[2] = ind[0] + 1;
1843
    }
1844
    else if(nonodes == 4) {
1845
      globalcoord[2] = (2.0*globalcoord[0]+globalcoord[1])/3.0;
1846
      globalcoord[6] = (2.0*globalcoord[4]+globalcoord[5])/3.0;
1847
      ind[2] = ind[0] + 1;
1848
      globalcoord[3] = (globalcoord[0]+2.0*globalcoord[1])/3.0;
1849
      globalcoord[7] = (globalcoord[4]+2.0*globalcoord[5])/3.0;
1850
      ind[3] = ind[0] + 2;
1851
    }
1852

1853
    return(elemind);
1854
  }
1855
  else if(numbering == NUMBER_XY) 
1856
    elemind = cell->elem1st+(i-1) + (j-1)*cell->elemwidth;
1857
  else if(numbering == NUMBER_YX)
1858
    elemind = cell->elem1st+(j-1) + (i-1)*cell->elemwidth;
1859

1860

1861
  ind[TOPRIGHT] = GetKnotCoordinate(cell,i,j,&globalcoord[TOPRIGHT],
1862
				    &globalcoord[TOPRIGHT+nonodes]);
1863
  ind[BOTRIGHT] = GetKnotCoordinate(cell,i,j-1,&globalcoord[BOTRIGHT],
1864
				    &globalcoord[BOTRIGHT+nonodes]);
1865
  ind[TOPLEFT] = GetKnotCoordinate(cell,i-1,j,&globalcoord[TOPLEFT],
1866
				    &globalcoord[TOPLEFT+nonodes]);
1867
  ind[BOTLEFT] = GetKnotCoordinate(cell,i-1,j-1,&globalcoord[BOTLEFT],
1868
				    &globalcoord[BOTLEFT+nonodes]);
1869
  if(nonodes == 4) return(elemind);
1870

1871
  GetElementIndices(cell,i,j,ind);
1872

1873
#if 0  
1874
  if(cell->xlinear) 
1875
    xrat = 1.0;
1876
  else
1877
    xrat = sqrt(cell->xratio);
1878

1879
  if(cell->ylinear) 
1880
    yrat = 1.0;
1881
  else 
1882
    yrat = sqrt(cell->yratio);
1883
#endif
1884
  xrat = yrat = 1.0;
1885

1886
  if(nonodes == 5) {
1887
    globalcoord[4] = 0.5*(globalcoord[0]+globalcoord[2]);
1888
    globalcoord[4+k] = 0.5*(globalcoord[k]+globalcoord[2+k]);
1889
  }
1890
  else if(nonodes == 8 || nonodes == 9) {
1891
    globalcoord[4] = (xrat*globalcoord[0]+globalcoord[1])/(1+xrat);
1892
    globalcoord[4+k] = globalcoord[k];    
1893
    globalcoord[5] = globalcoord[1];
1894
    globalcoord[5+k] = 
1895
      (yrat*globalcoord[1+k]+globalcoord[2+k])/(1+yrat);    
1896
    globalcoord[6] = globalcoord[4];
1897
    globalcoord[6+k] = globalcoord[2+k];
1898
    globalcoord[7] = globalcoord[3];
1899
    globalcoord[7+k] = globalcoord[5+k];
1900
    if(nonodes == 9) {
1901
      globalcoord[8] = globalcoord[4];
1902
      globalcoord[8+k] = globalcoord[5+k];
1903
    }
1904
  }
1905
  else if(nonodes == 12 || nonodes == 16) {
1906
    globalcoord[4] = (2.*globalcoord[0]+globalcoord[1])/3.0;
1907
    globalcoord[4+k] = (2.*globalcoord[k]+globalcoord[1+k])/3.0;
1908
    globalcoord[5] = (2.*globalcoord[1]+globalcoord[0])/3.0;
1909
    globalcoord[5+k] = (2.*globalcoord[1+k]+globalcoord[k])/3.0;
1910
    globalcoord[6] = (2.*globalcoord[1]+globalcoord[2])/3.0;
1911
    globalcoord[6+k] = (2.*globalcoord[1+k]+globalcoord[2+k])/3.0;
1912
    globalcoord[7] = (2.*globalcoord[2]+globalcoord[1])/3.0;
1913
    globalcoord[7+k] = (2.*globalcoord[2+k]+globalcoord[1+k])/3.0;
1914
    globalcoord[8] = (2.*globalcoord[2]+globalcoord[3])/3.0;
1915
    globalcoord[8+k] = (2.*globalcoord[2+k]+globalcoord[3+k])/3.0;
1916
    globalcoord[9] = (2.*globalcoord[3]+globalcoord[2])/3.0;
1917
    globalcoord[9+k] = (2.*globalcoord[3+k]+globalcoord[2+k])/3.0;
1918
    globalcoord[10] = (2.*globalcoord[3]+globalcoord[0])/3.0;
1919
    globalcoord[10+k] = (2.*globalcoord[3+k]+globalcoord[k])/3.0;
1920
    globalcoord[11] = (2.*globalcoord[0]+globalcoord[3])/3.0;
1921
    globalcoord[11+k] = (2.*globalcoord[k]+globalcoord[3+k])/3.0;
1922
    if(nonodes == 16) {
1923
      globalcoord[12] = (2.*globalcoord[11]+globalcoord[6])/3.0;
1924
      globalcoord[12+k] = (2.*globalcoord[11+k]+globalcoord[6+k])/3.0;
1925
      globalcoord[13] = (2.*globalcoord[6]+globalcoord[11])/3.0;
1926
      globalcoord[13+k] = (2.*globalcoord[6+k]+globalcoord[11+k])/3.0;
1927
      globalcoord[14] = (2.*globalcoord[7]+globalcoord[10])/3.0;
1928
      globalcoord[14+k] = (2.*globalcoord[7+k]+globalcoord[10+k])/3.0;
1929
      globalcoord[15] = (2.*globalcoord[10]+globalcoord[7])/3.0;
1930
      globalcoord[15+k] = (2.*globalcoord[10+k]+globalcoord[7+k])/3.0;      
1931
    }
1932
  }
1933

1934
#if 0
1935
  for(i=0;i<k;i++) 
1936
    printf("ind[%d]=%d  x=%.4lg  y=%.4lg\n",i,ind[i],
1937
	   globalcoord[i],globalcoord[i+k]);
1938
#endif
1939

1940
  return(elemind);
1941
}
1942

1943

1944
int GetSideInfo(struct CellType *cell,int cellno,int side,int element,
1945
		int *elemind)
1946
/* Given the side and element numbers returns the indices of 
1947
   the side element. When the end of the side is reached the function 
1948
   returns FALSE, otherwise TRUE.
1949
   */
1950
{
1951
  int more,sideno;                 
1952
  int ind[MAXNODESD2];
1953

1954
  more = TRUE;
1955

1956
  if(cell[cellno].numbering == NUMBER_1D) {
1957

1958
    switch(side) {
1959
    case 0:
1960
      more = FALSE;
1961
      elemind[0] = GetElementIndices(&(cell[cellno]),1,element,&(ind[0]));
1962
      if((sideno = cell[cellno].neighbour[LEFT]))
1963
	elemind[1] = GetElementIndices(&(cell[sideno]),cell[sideno].xelem,element,&(ind[0]));
1964
      else 
1965
	elemind[1] = 0;
1966
      break;
1967
     
1968
      
1969
    case 1:
1970
      more = FALSE; 
1971
      elemind[0] = GetElementIndices(&(cell)[cellno],cell[cellno].xelem,element,&(ind[0]));
1972
      if((sideno = cell[cellno].neighbour[RIGHT]))
1973
	elemind[1] = GetElementIndices(&(cell[sideno]),1,element,&(ind[0]));
1974
      else 
1975
	elemind[1] = 0;
1976
      break;
1977
    }
1978

1979
    return(more);
1980
  }
1981

1982

1983
  switch(side) {
1984

1985
  case LEFT:
1986
    if(element == cell[cellno].yelem) more = FALSE;
1987
    elemind[0] = GetElementIndices(&(cell[cellno]),1,element,&(ind[0]));
1988
    if((sideno = cell[cellno].neighbour[LEFT]))
1989
      elemind[1] = GetElementIndices(&(cell[sideno]),cell[sideno].xelem,element,&(ind[0]));
1990
    else 
1991
      elemind[1] = 0;
1992
    break;
1993

1994
  case RIGHT:
1995
    if(element == cell[cellno].yelem) more = FALSE; 
1996
    elemind[0] = GetElementIndices(&(cell)[cellno],cell[cellno].xelem,element,&(ind[0]));
1997
    if((sideno = cell[cellno].neighbour[RIGHT]))
1998
      elemind[1] = GetElementIndices(&(cell[sideno]),1,element,&(ind[0]));
1999
    else 
2000
      elemind[1] = 0;
2001
    break;
2002

2003
  case DOWN:
2004
    if(element == cell[cellno].xelem) more = FALSE;
2005
    elemind[0] = GetElementIndices(&(cell)[cellno],element,1,&(ind[0]));
2006
    if((sideno = cell[cellno].neighbour[DOWN]))
2007
      elemind[1] = GetElementIndices(&(cell[sideno]),element,cell[sideno].yelem,&(ind[0]));
2008
    else 
2009
      elemind[1] = 0;
2010
    break;
2011

2012
  case UP:
2013
    if(element == cell[cellno].xelem) more = FALSE; 
2014
    elemind[0] = GetElementIndices(&(cell)[cellno],element,cell[cellno].yelem,&(ind[0]));
2015
    if((sideno = cell[cellno].neighbour[UP]))
2016
      elemind[1] = GetElementIndices(&(cell[sideno]),element,1,&(ind[0]));
2017
    else 
2018
      elemind[1] = 0;
2019
    break;
2020
    
2021
  default:
2022
    printf("Impossible option in GetSideInfo.\n");
2023
  }
2024

2025
  return(more);
2026
}
2027

2028

2029

2030
void SetElementDivisionExtruded(struct GridType *grid,int info)
2031
{
2032
  int i,nzmax = 0,sumzelems;
2033
  Real ratio,linearlimit;
2034
  Real dzmax = 0,dz = 0;
2035
  
2036
  linearlimit = 0.001;
2037

2038
  if(grid->autoratio) {
2039
    for(i=1;i<=grid->zcells;i++)  
2040
      grid->zelems[i] = grid->minzelems;
2041
  }
2042
  else {
2043
    for(i=1;i<=grid->zcells;i++)  
2044
      if(grid->zelems[i] < grid->minzelems) grid->zelems[i] = grid->minzelems;
2045
  }
2046

2047
  sumzelems = grid->zcells * grid->minzelems;
2048
  if(sumzelems > grid->totzelems) {
2049
#if 0
2050
    printf("SetElementDivision: %d is too few elements in z-direction (min. %d)\n",
2051
	   grid->totzelems,sumzelems);
2052
#endif
2053
    grid->totzelems = sumzelems;
2054
  }      
2055

2056
  /* Put the linearity flags. */
2057
  for(i=1; i<= grid->zcells ;i++) {
2058
    if (fabs(1.-grid->zexpand[i]) < linearlimit) 
2059
      grid->zlinear[i] = TRUE; 
2060
    else 
2061
      grid->zlinear[i] = FALSE;
2062
  }
2063

2064
  if(grid->autoratio) {
2065
    int active;
2066
    for(;;) {
2067
      dzmax = 0.0;
2068
      active = FALSE;
2069

2070
      for(i=1;i<=grid->zcells;i++) {
2071
	if(grid->zelems[i] == 0) continue;
2072
	if(grid->zlinear[i] == TRUE)
2073
	  dz = (grid->z[i] - grid->z[i-1])/(grid->zelems[i]+1);
2074
	else {
2075
	  if(grid->zexpand[i] > 0.0) {
2076
	    ratio = pow(grid->zexpand[i],1./(1.*grid->zelems[i]));
2077
	    dz = (grid->z[i] - grid->z[i-1]) * 
2078
	      (1.-ratio) / (1.-pow(ratio,(Real)(grid->zelems[i]+1)));
2079
	    if(ratio < 1.)   
2080
	      dz *= grid->zexpand[i];
2081
	  }
2082
	  else if(grid->zelems[i]==1) {
2083
	    dz = (grid->z[i] - grid->z[i-1])/(grid->zelems[i]+1);
2084
	  } 
2085
	  else if((grid->zelems[i]+1)%2 == 0) {
2086
	    ratio = pow(-grid->zexpand[i],1./((grid->zelems[i]+1)/2-1.));
2087
	    dz = 0.5 * (grid->z[i] - grid->z[i-1]) * 
2088
	      (1.-ratio) / (1.-pow(ratio,(Real)((grid->zelems[i]+1)/2)));
2089
	  }
2090
	  else if((grid->zelems[i]+1)%2 == 1) {
2091
	    ratio = pow(-grid->zexpand[i],1./((grid->zelems[i]+1)/2));
2092
	    dz = (grid->z[i] - grid->z[i-1]) / 
2093
	      (2.0*(1.-pow(ratio,(Real)((grid->zelems[i]+1)/2)))/
2094
	       (1-ratio) + pow(ratio,(Real)((grid->zelems[i]+1)/2+0.5)));
2095
	  }
2096
	}
2097
	dz *= grid->zdens[i] / grid->xzratio;
2098

2099
	if(dz > dzmax) {
2100
	  dzmax = dz;
2101
	  nzmax = i;
2102
	}
2103

2104
	if(grid->autoratio) {
2105
	  if(fabs(dz - grid->dx0) < fabs( dz*(1.0/(1.0-1.0/grid->zelems[i])) - grid->dx0) ) {
2106
	    grid->zelems[i] += 1;
2107
	    sumzelems++;
2108
	    active = TRUE;
2109
	  }
2110
	}
2111
      }
2112

2113
      if(grid->autoratio) {
2114
	if(!active) break;
2115
      }
2116
      else {
2117
	grid->zelems[nzmax] += 1;
2118
	sumzelems++;
2119
	if(sumzelems >= grid->totzelems) break;
2120
      }
2121
    }
2122
  }
2123

2124
  /* Set the size of the first element within each subcell */
2125
  grid->totzelems = 0;
2126
  for(i=1;i<=grid->zcells;i++) {
2127
    grid->totzelems += grid->zelems[i];
2128
    if(grid->zlinear[i] == TRUE || grid->zelems[i]==1) 
2129
      grid->dz[i] = (grid->z[i] - grid->z[i-1])/grid->zelems[i];
2130
    else if(grid->zexpand[i] > 0.0) {
2131
      ratio = pow(grid->zexpand[i],1./(grid->zelems[i]-1.));
2132
      grid->zratios[i] = ratio;
2133
      grid->dz[i] = (grid->z[i] - grid->z[i-1]) * 
2134
	(1.-ratio) / (1.-pow(ratio,(Real)(grid->zelems[i])));
2135
    }	
2136
    else if(grid->zelems[i] == 2) {
2137
      grid->dz[i] = (grid->z[i] - grid->z[i-1])/grid->zelems[i];
2138
    }
2139
    else if(grid->zelems[i]%2 == 0) {
2140
      ratio = pow(-grid->zexpand[i],1./(grid->zelems[i]/2-1.));
2141
      grid->zratios[i] = ratio;
2142
	grid->dz[i] = 0.5 * (grid->z[i] - grid->z[i-1]) * 
2143
	  (1.-ratio) / (1.-pow(ratio,(Real)(grid->zelems[i]/2)));
2144
    }
2145
    else if(grid->zelems[i]%2 == 1) {
2146
      ratio = pow(-grid->zexpand[i],1./(grid->zelems[i]/2));
2147
      grid->zratios[i] = ratio;
2148
      grid->dz[i] = (grid->z[i] - grid->z[i-1]) / 
2149
	(2.0*(1.-pow(ratio,(Real)(grid->zelems[i]/2)))/
2150
	 (1-ratio) + pow(ratio,(Real)(grid->zelems[i]/2+0.5)));
2151
    }
2152
  }
2153

2154
  if(info) printf("Created %d extruded divisions.\n",
2155
		  grid->totzelems);
2156
  grid->dz0 = dzmax;
2157
}
2158

2159

2160

2161
void SetElementDivisionCylinder(struct GridType *grid,int info)
2162
{
2163
  int i,k;
2164
  Real ratio,eps;
2165
  Real dzmax;
2166
  
2167
  eps = 1.0e-8;
2168

2169
  grid->zcells = 2*grid->rotateblocks;
2170
  grid->z[0] = 0.0;
2171
  for(i=0;grid->x[i]<eps;i++); k=i;  
2172
  grid->rotateradius1 = grid->x[k];
2173

2174
  if(grid->rotateradius2 < 0.0) {
2175
    grid->rotatecartesian = TRUE;
2176
    grid->rotateradius2 = 1.0e6;
2177
  }
2178
  if(grid->rotateradius2 <= sqrt(2.0)*grid->rotateradius1) {
2179
    if(k+1 <= grid->xcells) 
2180
      grid->rotateradius2 = grid->x[k+1];
2181
    grid->rotateradius2 = MAX(sqrt(2.0)*grid->rotateradius1,
2182
			      grid->rotateradius2);
2183
  }
2184

2185
  if(!grid->xlinear[k]) 
2186
    printf("SetElementDivisionCylinder: The division must be linear!\n");
2187

2188
  for(i=1;i<=grid->zcells;i++) {
2189
    grid->z[i] = i*(2.0*FM_PI)/8.0; 
2190
    grid->zelems[i]  = grid->xelems[k];
2191
    grid->zlinear[i] = grid->xlinear[k];
2192
    grid->zratios[i] = grid->xratios[k];
2193
    grid->zexpand[i] = grid->xexpand[k];
2194
    grid->zmaterial[i] = 0;
2195
  }
2196

2197
  grid->totzelems =  grid->zcells * grid->xelems[k];
2198

2199
  dzmax = 0.0;
2200
  for(i=1;i<=grid->zcells;i++) {
2201
    if(grid->zlinear[i] == TRUE || grid->zelems[i]==1) 
2202
      grid->dz[i] = (grid->z[i] - grid->z[i-1])/grid->zelems[i];
2203
    else if(grid->zexpand[i] > 0.0) {
2204
      ratio = pow(grid->zexpand[i],1./(grid->zelems[i]-1.));
2205
      grid->zratios[i] = ratio;
2206
      grid->dz[i] = (grid->z[i] - grid->z[i-1]) * 
2207
	(1.-ratio) / (1.-pow(ratio,(Real)(grid->zelems[i])));
2208
    }	
2209
    else if(grid->zelems[i] == 2) {
2210
      grid->dz[i] = (grid->z[i] - grid->z[i-1])/grid->zelems[i];
2211
    }
2212
    else if(grid->zelems[i]%2 == 0) {
2213
      ratio = pow(-grid->zexpand[i],1./(grid->zelems[i]/2-1.));
2214
      grid->zratios[i] = ratio;
2215
	grid->dz[i] = 0.5 * (grid->z[i] - grid->z[i-1]) * 
2216
	  (1.-ratio) / (1.-pow(ratio,(Real)(grid->zelems[i]/2)));
2217
    }
2218
    else if(grid->zelems[i]%2 == 1) {
2219
      ratio = pow(-grid->zexpand[i],1./(grid->zelems[i]/2));
2220
      grid->zratios[i] = ratio;
2221
      grid->dz[i] = (grid->z[i] - grid->z[i-1]) / 
2222
	(2.0*(1.-pow(ratio,(Real)(grid->zelems[i]/2)))/
2223
	 (1-ratio) + pow(ratio,(Real)(grid->zelems[i]/2+0.5)));
2224
    }
2225

2226
    if(dzmax < grid->dz[i]) dzmax = grid->dz[i];
2227
  }
2228

2229
  if(info) printf("Created %d divisions in %d cells for rotation [%.2lg  %.2lg].\n",
2230
		  grid->totzelems,grid->zcells,
2231
		  grid->rotateradius1,grid->rotateradius2);
2232
  grid->dz0 = dzmax;
2233
}
2234

2235

2236

2237
static int GetCommand(char *line1,char *line2,FILE *io) 
2238
/* Line1 for commands and line2 for arguments. */
2239
{
2240
  int i,j,isend;
2241
  char line0[MAXLINESIZE],*charend,*matcpntr0,*matcpntr;
2242
  int gotlinefeed;
2243

2244
 newline:
2245

2246
  for(i=0;i<MAXLINESIZE;i++)
2247
    line2[i] = line1[i] = line0[i] = 0x00;
2248

2249
  charend = fgets(line0,MAXLINESIZE,io);
2250
  isend = (charend == NULL);
2251

2252
  if(isend) return(1);
2253

2254
  if(strlen(line0)<=strspn(line0," \t\r\n")) goto newline;
2255

2256
  if(line0[0] == '#' || line0[0] == '%' || line0[0] == '!' || line0[0] == '\n') goto newline;
2257
  if(!matcactive && line0[0] == '*') goto newline;
2258

2259
  if(!matcactive && strchr(line0,'$') ) {
2260
    #if USE_MATC
2261
      printf("\n a $ found and MATC has been compiled but not activated,\n");
2262
      printf("either remove the $ or add 'MATC = True' to grd input file.\n");
2263
    #else
2264
      printf("\n a $ found and MATC has not been compiled into ElmerGrid,\n");
2265
      printf("either remove the $ or compile ElmerGrid with MATC\n");
2266
    #endif // USE_MATC
2267
      bigerror("Error: $ found in command and MATC is not active");
2268
  }
2269

2270
#if USE_MATC
2271
  if(matcactive) {
2272
    matcpntr0 = strchr(line0,'$');
2273
    if(matcpntr0) {
2274
      matcpntr = mtc_domath(&matcpntr0[1]);
2275
      if(matcpntr) {
2276
	strcpy(matcpntr0, matcpntr);
2277
	if(0) printf("B: %s\n%s\n",matcpntr0,matcpntr);
2278
      }
2279
      else {
2280
	if(0) printf("B0: %s\n",matcpntr0);
2281
	goto newline;
2282
      }
2283
    }
2284
  }
2285
#endif 
2286
  
2287
  gotlinefeed = FALSE;
2288
  j = 0;
2289
  for(i=0;i<MAXLINESIZE;i++) {
2290
    if(line0[i] == '\n' || line0[i] == '\0' ) {
2291
      gotlinefeed = TRUE;
2292
      break;
2293
    }
2294
    if(line0[i] == '=') {
2295
      j = i;
2296
      break;
2297
    }
2298
    line1[i] = toupper(line0[i]);
2299
  }
2300

2301
  /* After these commands there will be no nextline even though there is no equality sign */
2302
  if(strstr(line1,"END")) return(0);
2303
  if(strstr(line1,"NEW MESH")) return(0);
2304

2305
  if(j) { /* Arguments are actually on the same line after '=' */
2306
    for(i=j+1;i<MAXLINESIZE;i++) {
2307
      line2[i-j-1] = line0[i];    
2308
      if( line0[i] == '\n' || line0[i] == '\0' ) {
2309
	gotlinefeed = TRUE;
2310
	break;
2311
      }
2312
    }  
2313
    if(!gotlinefeed) {
2314
      printf("There is a risk that somethings was missed in line:\n");
2315
      printf("%s\n",line0);
2316
      smallerror("Check your output line length!\n");
2317
    }
2318
  }
2319
  else { /* Arguments are on the next line */
2320
  newline2:
2321
    charend = fgets(line2,MAXLINESIZE,io);
2322
    isend = (charend == NULL);
2323
    if(isend) return(2);
2324
    if(line2[0] == '#' || line2[0] == '%' || line2[0] == '!') goto newline2;
2325
    if(!matcactive && line2[0] == '*') goto newline2;
2326

2327
    gotlinefeed = FALSE;
2328
    for(i=0;i<MAXLINESIZE;i++) {
2329
      if(line2[i] == '\n' || line2[i] == '\0' ) {
2330
	gotlinefeed = TRUE;
2331
	break;
2332
      }
2333
    }
2334
    if(!gotlinefeed) {
2335
      printf("There is a risk that somethings was missed in line:\n");
2336
      printf("%s\n",line2);
2337
      smallerror("Check your output line length!\n");
2338
    }
2339

2340
    if(!matcactive && strchr(line0,'$') ) {
2341
      #if USE_MATC
2342
        printf("\n a $ found and MATC has been compiled but not activated,\n");
2343
        printf("either remove the $ or add 'MATC = True' to grd input file.\n");
2344
      #else
2345
        printf("\n a $ found and MATC has not been compiled into ElmerGrid,\n");
2346
        printf("either remove the $ or compile ElmerGrid with MATC\n");
2347
      #endif // USE_MATC
2348
        bigerror("Error: $ found in command and MATC is not active");
2349
    }
2350

2351
#if USE_MATC
2352
    if(matcactive) {
2353
      matcpntr0 = strchr(line2,'$');
2354
      if(matcpntr0) {
2355
	matcpntr = mtc_domath(&matcpntr0[1]);
2356
	if(matcpntr) {
2357
	  strcpy(matcpntr0, matcpntr);
2358
	  if(0) printf("C: %s\n%s\n",matcpntr0,matcpntr);
2359
	}
2360
      }
2361
    }
2362
#endif 
2363
  }
2364

2365
  if(iodebug) {
2366
    printf("command: %s\n",line1);
2367
    printf("params:  %s\n",line2);
2368
  }
2369
  
2370
  return(0);
2371
}
2372

2373

2374

2375
int SaveElmergrid(struct GridType *grid,int nogrids,char *prefix,int info)
2376
{
2377
  int sameline,maxsameline;
2378
  int i,j,dim;
2379
  FILE *out;
2380
  char filename[MAXFILESIZE];
2381

2382
  AddExtension(prefix,filename,"grd");
2383
  out = fopen(filename,"w");
2384
  dim = grid->dimension;
2385
  if(grid->coordsystem == COORD_CART1) dim = 1;
2386

2387
  j = 0;
2388
  sameline = TRUE;
2389
  maxsameline = 6;
2390
  if(grid->xcells > maxsameline) sameline = FALSE;
2391
  if(dim >= 2 && grid->ycells > maxsameline) sameline = FALSE;
2392
  if(dim >= 3 && grid->zcells > maxsameline) sameline = FALSE;
2393
  
2394
  fprintf(out,"#####  ElmerGrid input file for structured grid generation  ######\n");
2395
  fprintf(out,"Version = 210903\n");
2396

2397
  fprintf(out,"Coordinate System = ");
2398
  if(grid->coordsystem == COORD_AXIS)
2399
    fprintf(out,"2D Axisymmetric\n");
2400
  else if(grid->coordsystem == COORD_POLAR)
2401
    fprintf(out,"2D Polar\n");
2402
  else 
2403
    fprintf(out,"Cartesian %dD\n",dim);
2404
 
2405
  fprintf(out,"Subcell Divisions in %dD = ",dim);
2406
  if(dim >= 1) fprintf(out,"%d ",grid->xcells);
2407
  if(dim >= 2) fprintf(out,"%d ",grid->ycells);
2408
  if(dim >= 3) fprintf(out,"%d ",grid->zcells);
2409
  fprintf(out,"\n");
2410

2411
  fprintf(out,"Subcell Limits 1 %s",sameline ? "= ":"\n  ");
2412
  for(i=0;i <= grid->xcells;i++) 
2413
    fprintf(out,"%.5lg ",grid->x[i]); 
2414
  fprintf(out,"\n");
2415
    
2416
  if(dim >= 2) {
2417
    fprintf(out,"Subcell Limits 2 %s",sameline ? "= ":"\n  ");
2418
    for(i=0;i <= grid->ycells;i++) 
2419
      fprintf(out,"%.5lg ",grid->y[i]); 
2420
    fprintf(out,"\n");
2421
  }
2422
  
2423
  if(dim >= 3) {
2424
    fprintf(out,"Subcell Limits 3 %s",sameline ? "= ":"\n  ");
2425
    for(i=0;i <= grid->zcells;i++) 
2426
      fprintf(out,"%.5lg ",grid->z[i]); 
2427
    fprintf(out,"\n");
2428
  }  
2429

2430
  fprintf(out,"Material Structure in %dD\n",dim==1 ? 1:2);  
2431
  for(j=grid->ycells;j>=1;j--) {
2432
    fprintf(out,"  ");
2433
    for(i=1;i<=grid->xcells;i++) 
2434
      fprintf(out,"%-5d",grid->structure[j][i]);
2435
    fprintf(out,"\n");
2436
  }
2437
  fprintf(out,"End\n");
2438

2439
  if(grid->mappings > 0) {
2440
    fprintf(out,"Geometry Mappings\n");
2441
    fprintf(out,"# mode  line  limits(2)   Np  params(Np)\n");
2442
    for(i=0;i<grid->mappings;i++) {
2443
      fprintf(out,"  %-5d %-5d %-7.5lg %-7.5lg %-3d ",
2444
	      grid->mappingtype[i],grid->mappingline[i],
2445
	      grid->mappinglimits[2*i],grid->mappinglimits[2*i+1],
2446
	      grid->mappingpoints[i]);
2447
      for(j=0;j<grid->mappingpoints[i];j++) 
2448
	fprintf(out,"%.4lg ",grid->mappingparams[i][j]);
2449
      fprintf(out,"\n");
2450
    }
2451
    fprintf(out,"End\n");
2452
  }
2453

2454
  j = 0;
2455
  if(grid[j].rotate) {
2456
    fprintf(out,"Revolve Blocks = %d\n",grid[j].rotateblocks);
2457
    fprintf(out,"Revolve Radius = %-8.3lg\n",grid[j].rotateradius2);
2458
    if(fabs(grid[j].rotateimprove-1.0) > 1.0e-10)
2459
      fprintf(out,"Revolve Improve = %-8.3lg\n",grid[j].rotateimprove);
2460
    
2461
  }
2462
  if(grid[j].rotatecurve) {
2463
    fprintf(out,"Revolve Curve Direct = %-8.3lg\n",grid[j].curvezet);
2464
    fprintf(out,"Revolve Curve Radius = %-8.3lg\n",grid[j].curverad);
2465
    fprintf(out,"Revolve Curve Angle = %-8.3lg\n",grid[j].curveangle);
2466
  }
2467

2468
  if(grid[j].coordsystem == COORD_POLAR) {
2469
    fprintf(out,"Polar Radius = %.3lg\n",grid[j].polarradius);
2470
  } 
2471

2472
  for(j=0;j<nogrids;j++) {
2473
    
2474
    if(j>0) fprintf(out,"\nStart New Mesh\n");
2475
  
2476
    fprintf(out,"Materials Interval = %d %d\n",
2477
	    grid[j].firstmaterial,grid[j].lastmaterial);
2478
  
2479
    if(dim == 3) {
2480
      fprintf(out,"Extruded Structure\n");
2481
      fprintf(out,"# %-8s %-8s %-8s\n","1stmat", "lastmat","newmat");
2482
      for(i=1;i<=grid[j].zcells;i++) 
2483
	fprintf(out,"  %-8d %-8d %-8d\n",
2484
		grid[j].zfirstmaterial[i],grid[j].zlastmaterial[i],
2485
		grid[j].zmaterial[i]); 
2486
      fprintf(out,"End\n");    
2487
    }
2488

2489
    if(grid[j].noboundaries > 0) {
2490
      fprintf(out,"Boundary Definitions\n");
2491
      fprintf(out,"# %-8s %-8s %-8s\n","type","out","int"); 
2492
      for(i=0;i<grid[j].noboundaries;i++)
2493
	fprintf(out,"  %-8d %-8d %-8d %-8d\n",
2494
		grid[j].boundtype[i],grid[j].boundext[i],
2495
		grid[j].boundint[i], grid[j].boundsolid[i]);
2496
      fprintf(out,"End\n");
2497
    }
2498

2499
    if(grid->numbering == NUMBER_XY)
2500
      fprintf(out,"Numbering = Horizontal\n");
2501
    if(grid->numbering == NUMBER_YX)
2502
      fprintf(out,"Numbering = Vertical\n");
2503
        
2504
    fprintf(out,"Element Degree = %d\n",grid[j].elemorder);
2505
    fprintf(out,"Element Innernodes = %s\n",grid[j].elemmidpoints ? "True" : "False");
2506
    fprintf(out,"Triangles = %s\n",grid[j].triangles ? "True" : "False");
2507
    if(grid[j].autoratio) 
2508
      fprintf(out,"Surface Elements = %d\n",grid[j].wantedelems);
2509
    if(dim == 3 && grid[j].wantedelems3d) 
2510
      fprintf(out,"Volume Elements = %d\n",grid[j].wantedelems3d);
2511
    if(dim == 3 && grid[j].wantednodes3d) 
2512
      fprintf(out,"Volume Nodes = %d\n",grid[j].wantednodes3d);
2513

2514
    if(dim == 2)
2515
      fprintf(out,"Coordinate Ratios = %-8.3lg\n",grid[j].xyratio);
2516
    if(dim == 3)
2517
      fprintf(out,"Coordinate Ratios = %-8.3lg %-8.3lg\n",
2518
	      grid[j].xyratio,grid[j].xzratio);
2519
 
2520
    fprintf(out,"Minimum Element Divisions = %d",grid[j].minxelems);
2521
    if(dim >= 2) fprintf(out," %d",grid[j].minyelems);
2522
    if(dim >= 3) fprintf(out," %d",grid[j].minzelems);
2523
    fprintf(out,"\n");
2524

2525
    fprintf(out,"Element Ratios 1 %s",sameline ? "= ":"\n  ");
2526
    for(i=1;i<=grid[j].xcells;i++) 
2527
      fprintf(out,"%.3lg ",grid[j].xexpand[i]); 
2528
    fprintf(out,"\n");
2529
    if(dim >= 2) {
2530
      fprintf(out,"Element Ratios 2 %s",sameline ? "= ":"\n  ");
2531
      for(i=1;i<=grid[j].ycells;i++) 
2532
	fprintf(out,"%.3lg ",grid[j].yexpand[i]); 
2533
      fprintf(out,"\n");
2534
    }
2535
    if(dim >= 3) {
2536
      fprintf(out,"Element Ratios 3 %s",sameline ? "= ":"\n  ");
2537
      for(i=1;i<=grid[j].zcells;i++) 
2538
	fprintf(out,"%.3lg ",grid[j].zexpand[i]); 
2539
      fprintf(out,"\n");
2540
    }
2541

2542
    if(grid[j].autoratio) {
2543
      fprintf(out,"Element Densities 1 %s",sameline ? "= ":"\n  ");
2544
      for(i=1;i<=grid[j].xcells;i++) 
2545
	fprintf(out,"%.3lg ",grid[j].xdens[i]); 
2546
      fprintf(out,"\n");
2547
      if(dim >= 2) {
2548
	fprintf(out,"Element Densities 2 %s",sameline ? "= ":"\n  ");
2549
	for(i=1;i<=grid[j].ycells;i++) 
2550
	  fprintf(out,"%.3lg ",grid[j].ydens[i]); 
2551
	fprintf(out,"\n");
2552
      }
2553
      if(dim >= 3) {       
2554
	fprintf(out,"Element Densities 3 %s",sameline ? "= ":"\n  ");
2555
	for(i=1;i<=grid[j].zcells;i++) 
2556
	  fprintf(out,"%.3lg ",grid[j].zdens[i]); 
2557
	fprintf(out,"\n");
2558
      }
2559
    }
2560
    else {
2561
      fprintf(out,"Element Divisions 1 %s",sameline ? "= ":"\n  ");
2562
      for(i=1;i<=grid[j].xcells;i++) 
2563
	fprintf(out,"%d ",grid[j].xelems[i]); 
2564
      fprintf(out,"\n");
2565
      if(dim >= 2) {
2566
	fprintf(out,"Element Divisions 2 %s",sameline ? "= ":"\n  ");
2567
	for(i=1;i<=grid[j].ycells;i++) 
2568
	  fprintf(out,"%d ",grid[j].yelems[i]); 
2569
	fprintf(out,"\n");
2570
      }
2571
      if(dim >= 3) {       
2572
	fprintf(out,"Element Divisions 3 %s",sameline ? "= ":"\n  ");
2573
	for(i=1;i<=grid[j].zcells;i++) 
2574
	  fprintf(out,"%d ",grid[j].zelems[i]); 
2575
	fprintf(out,"\n");
2576
      }
2577
    }
2578
    
2579

2580
  }
2581

2582
  if(info) printf("The Elmergrid input was saved to file %s.\n",filename);
2583
  fclose(out);
2584

2585
  return(0);
2586
}
2587

2588

2589

2590

2591

2592
static int LoadElmergridOld(struct GridType **grid,int *nogrids,char *prefix,int info) 
2593
{
2594
  char filename[MAXFILESIZE];
2595
  FILE *in;
2596
  int i,j,k,l,error=0;
2597
  Real scaling;
2598
  char *cp;
2599
  int mode,noknots,noelements,dim,axisymmetric;
2600
  int elemcode,maxnodes,totelems,nogrids0;
2601
  int minmat,maxmat;
2602
  char line[MAXLINESIZE];
2603

2604

2605
  AddExtension(prefix,filename,"grd");
2606
  if ((in = fopen(filename,"r")) == NULL) {
2607
    printf("LoadElmergrid: opening of the file '%s' wasn't successful !\n",filename);
2608
    return(1);
2609
  }
2610

2611
  if(info) printf("Loading the geometry from file '%s'.\n",filename);
2612

2613
  InitGrid(grid[*nogrids]);
2614
  k = *nogrids;
2615
  nogrids0 = *nogrids;
2616

2617
  mode = 0;
2618
  noknots = 0;
2619
  noelements = 0;
2620
  dim = 0;
2621
  axisymmetric = FALSE;
2622
  elemcode = 0;
2623
  maxnodes = 4;
2624
  totelems = 0;
2625
  scaling = 1.0;
2626

2627

2628

2629
  Getline(line,in);
2630
  for(;;) {
2631
    if(Getline(line,in)) goto end;
2632
    if(line[0]=='\0') goto end;
2633

2634
    if(strstr(line,"END")) goto end;
2635
    if(strstr(line,"RESULTS")) goto end;
2636

2637
    /* Control information */
2638
    if(strstr(line,"VERSION")) mode = 1;
2639
    else if(strstr(line,"GEOMETRY")) mode = 2;
2640
    else if(strstr(line,"MAPPINGS IN")) mode = 31;
2641
    else if(strstr(line,"MAPPINGS OUT")) mode = 32;
2642
    else if(strstr(line,"MAPPINGS")) mode = 3;
2643
    else if(strstr(line,"NUMBERING")) mode = 4;
2644
    else if(strstr(line,"MESHING")) mode = 5;
2645
    else if(strstr(line,"ELEMENTS")) mode = 6;
2646
    else if(strstr(line,"ELEMENT NUMBER")) mode = 29;
2647
    else if(strstr(line,"NODES")) mode = 7;
2648
    else if(strstr(line,"TRIANGLE")) mode = 8;
2649
    else if(strstr(line,"SQUARE")) mode = 17;
2650
    else if(strstr(line,"COORDINATE RATIO"))  mode = 10;
2651
    else if(strstr(line,"MATERIALS")) mode = 11;
2652
    else if(strstr(line,"LAYERED ST")) mode = 12;
2653
    else if(strstr(line,"ELEMENT RAT")) mode = 13;
2654
    else if(strstr(line,"ELEMENT DENS")) mode = 14;
2655
    else if(strstr(line,"ELEMENT MINIMUM")) mode = 27;
2656
    else if(strstr(line,"BOUNDARY COND")) mode = 15;
2657
    else if(strstr(line,"ELEMENTTYPE") || strstr(line,"ELEMENTCODE")) mode = 16;
2658
    else if(strstr(line,"ROTATE")) mode = 20;
2659
    else if(strstr(line,"ROTRAD")) mode = 21;
2660
    else if(strstr(line,"ROTBLOCK")) mode = 22;
2661
    else if(strstr(line,"ROTIMP")) mode = 24;
2662
    else if(strstr(line,"ROTCURVE")) mode = 25;
2663
    else if(strstr(line,"REDUCE ELEMENT")) mode = 26;
2664
    else if(strstr(line,"SCALING")) mode = 23;
2665
    else if(strstr(line,"LAYERED BO")) mode = 28;
2666
    else if(strstr(line,"POLAR RADIUS")) mode = 30;
2667

2668

2669
    switch (mode) {
2670
    case 1: 
2671
      printf("Loading Elmergrid file: %s\n",line);
2672
      mode = 0;
2673
      break;
2674
      
2675
    case 2:
2676
      grid[k]->dimension = 2;
2677
      if(strstr(line,"CARTES") && strstr(line,"1D")) {
2678
	grid[k]->coordsystem = COORD_CART1;
2679
	grid[k]->dimension = 1;
2680
      }
2681
      else if(strstr(line,"CARTES") && strstr(line,"2D")) 
2682
	grid[k]->coordsystem = COORD_CART2;
2683
      else if(strstr(line,"AXIS") && strstr(line,"2D")) 
2684
	grid[k]->coordsystem = COORD_AXIS;
2685
      else if(strstr(line,"POLAR") && strstr(line,"2D")) 
2686
	grid[k]->coordsystem = COORD_POLAR;
2687
      else if(strstr(line,"CARTES") && strstr(line,"3D")) {
2688
	grid[k]->coordsystem = COORD_CART3;
2689
	grid[k]->dimension = 3;
2690
      }
2691
      else if(strstr(line,"CYLINDRICAL")) {
2692
	grid[k]->coordsystem = COORD_CYL;
2693
	grid[k]->dimension = 3;
2694
      }
2695
      else printf("Unknown coordinate system: %s\n",line);
2696
      printf("Defining the coordinate system (%d-DIM).\n",grid[k]->dimension);
2697

2698
      Getline(line,in);
2699

2700
      if(grid[k]->dimension == 1) {
2701
	sscanf(line,"%d",&(*grid)[k].xcells);
2702
	grid[k]->ycells = 1;	
2703
      }
2704
      if(grid[k]->dimension == 2) 
2705
	sscanf(line,"%d %d",&(*grid)[k].xcells,&(*grid)[k].ycells);
2706
      if(grid[k]->dimension == 3) 
2707
	sscanf(line,"%d %d %d",&(*grid)[k].xcells,&(*grid)[k].ycells,&(*grid)[k].zcells);      
2708
      if(grid[k]->xcells >= MAXCELLS || grid[k]->ycells >= MAXCELLS || 
2709
	 grid[k]->zcells >= MAXCELLS) {
2710
	printf("LoadGrid: Too many subcells [%d %d %d] vs. %d:\n",
2711
	       grid[k]->xcells,grid[k]->ycells,grid[k]->zcells,MAXCELLS);
2712
      }
2713
      
2714
      if(grid[k]->dimension == 1) {
2715
	printf("Loading [%d] subcell intervals in 1D\n",
2716
	       grid[k]->xcells);
2717
      }
2718
      else if(grid[k]->dimension == 2) {
2719
	printf("Loading [%d %d] subcell intervals in 2D\n",
2720
	       grid[k]->xcells,grid[k]->ycells);   
2721
      } else {
2722
	printf("Loading [%d %d %d] subcell intervals in 3D\n",
2723
	       grid[k]->xcells,grid[k]->ycells,grid[k]->zcells);   
2724
      }
2725

2726

2727
      for(j=1;j<=grid[k]->dimension;j++) {
2728
	Getline(line,in);
2729
	cp=line;
2730

2731
	if(j==1) for(i=0;i<=grid[k]->xcells;i++) grid[k]->x[i] = next_real(&cp);
2732
	if(j==2) for(i=0;i<=grid[k]->ycells;i++) grid[k]->y[i] = next_real(&cp);
2733
	if(j==3) for(i=0;i<=grid[k]->zcells;i++) grid[k]->z[i] = next_real(&cp);
2734
      }
2735

2736
      printf("Loading material structure\n");
2737

2738
      for(j=grid[k]->ycells;j>=1;j--) {
2739
	
2740
	Getline(line,in);
2741
	cp=line;
2742
	
2743
	for(i=1;i<=grid[k]->xcells;i++) 
2744
	  grid[k]->structure[j][i] = next_int(&cp);
2745
      }
2746

2747
      minmat = maxmat = grid[k]->structure[1][1];
2748
      for(j=grid[k]->ycells;j>=1;j--) 
2749
	for(i=1;i<=grid[k]->xcells;i++) {
2750
	  if(minmat > grid[k]->structure[j][i])
2751
	    minmat = grid[k]->structure[j][i];
2752
	  if(maxmat < grid[k]->structure[j][i])
2753
	    maxmat = grid[k]->structure[j][i];
2754
	}      
2755
      if(minmat < 0) 
2756
	printf("LoadElmergrid: please use positive material indices.\n");
2757

2758
      mode = 0;
2759
      break;
2760

2761
    case 3:
2762
    case 31:
2763
    case 32:
2764

2765
      /* I don't know how to set this, luckily this piece of code should be obsolete */
2766
      l = 1;
2767
      for(i=grid[k]->mappings;i<grid[k]->mappings+l;i++) {
2768
	Getline(line,in);
2769
	cp=line; 
2770

2771
	grid[k]->mappingtype[i] = next_int(&cp);
2772
	if(mode == 32) grid[k]->mappingtype[i] += 50*SGN(grid[k]->mappingtype[i]);
2773

2774
	grid[k]->mappingline[i] = next_int(&cp);
2775
	grid[k]->mappinglimits[2*i] = next_real(&cp);
2776
	grid[k]->mappinglimits[2*i+1] = next_real(&cp);
2777
	grid[k]->mappingpoints[i] = next_int(&cp);
2778
	grid[k]->mappingparams[i] = Rvector(0,grid[k]->mappingpoints[i]);
2779
	for(j=0;j<grid[k]->mappingpoints[i];j++) 
2780
	  grid[k]->mappingparams[i][j] = next_real(&cp);
2781
      }
2782
      
2783
      printf("Loaded %d geometry mappings\n",l);
2784
      grid[k]->mappings += l;
2785

2786
      mode = 0;
2787
      break;
2788
      
2789
    case 4: /* NUMBERING */
2790
      if(strstr(line,"HORIZ")) grid[k]->numbering = NUMBER_XY;
2791
      if(strstr(line,"VERTI")) grid[k]->numbering = NUMBER_YX;
2792
      mode = 0;
2793
      break;
2794

2795
    case 5: /* MESHING */
2796
      if((*nogrids) >= MAXCASES) {
2797
	printf("There are more grids than was allocated for!\n"); 
2798
	printf("Ignoring meshes starting from %d\n.",(*nogrids)+1);
2799
	goto end;
2800
      }
2801
      (*nogrids)++;
2802
      printf("Loading element meshing no %d\n",*nogrids);
2803
      k = *nogrids - 1;	           
2804
      if(k > nogrids0) (*grid)[k] = (*grid)[k-1];	 
2805
      mode = 0;
2806
      break;
2807

2808
    case 6: /* ELEMENTS */
2809
      sscanf(line,"%d",&(*grid)[k].wantedelems);
2810
      mode = 0;
2811
      break;
2812

2813
    case 7: /* NODES */
2814
      sscanf(line,"%d",&(*grid)[k].nonodes);      
2815
      
2816
      (*grid)[k].elemmidpoints = FALSE;
2817
      if((*grid)[k].nonodes == 4) 
2818
	(*grid)[k].elemorder = 1;
2819
      if((*grid)[k].nonodes == 8) 
2820
	(*grid)[k].elemorder = 2;
2821
      if((*grid)[k].nonodes == 16) 
2822
	(*grid)[k].elemorder = 3;
2823

2824
      if((*grid)[k].nonodes == 9) { 
2825
	(*grid)[k].elemorder = 2;
2826
	(*grid)[k].elemmidpoints = TRUE;
2827
      }
2828
      if((*grid)[k].nonodes == 12) { 
2829
	(*grid)[k].elemorder = 3;
2830
	(*grid)[k].elemmidpoints = TRUE;
2831
      }
2832

2833

2834
      mode = 0;
2835
      break;
2836

2837
    case 8: /* TRIANGLES */
2838
      (*grid)[k].triangles = TRUE;
2839
      mode = 0;
2840
      break;
2841

2842
    case 17: /* SQUARES */
2843
      (*grid)[k].triangles = FALSE;
2844
      mode = 0;
2845
      break;
2846

2847
    case 16: /* ELEMENTTYPE and ELEMENTCODE */
2848
      sscanf(line,"%d",&elemcode);
2849
      if(elemcode/100 == 2) {
2850
	(*grid)[k].triangles = FALSE;      
2851
	(*grid)[k].nonodes = elemcode%100;
2852
      }
2853
      else if(elemcode/100 == 4) {
2854
	(*grid)[k].triangles = FALSE;      
2855
	(*grid)[k].nonodes = elemcode%100;
2856
      }
2857
      else if(elemcode/100 == 3) {  
2858
	(*grid)[k].triangles = TRUE;      
2859
	if(elemcode%100 == 3)       (*grid)[k].nonodes = 4;
2860
	else if(elemcode%100 == 6)  (*grid)[k].nonodes = 9;
2861
	else if(elemcode%100 == 10) (*grid)[k].nonodes = 16;	
2862
      }
2863

2864
      (*grid)[k].elemmidpoints = FALSE;
2865
      if((*grid)[k].nonodes == 4) 
2866
	(*grid)[k].elemorder = 1;
2867
      if((*grid)[k].nonodes == 8) 
2868
	(*grid)[k].elemorder = 2;
2869
      if((*grid)[k].nonodes == 16) 
2870
	(*grid)[k].elemorder = 3;
2871

2872
      if((*grid)[k].nonodes == 9) { 
2873
	(*grid)[k].elemorder = 2;
2874
	(*grid)[k].elemmidpoints = TRUE;
2875
      }
2876
      if((*grid)[k].nonodes == 12) { 
2877
	(*grid)[k].elemorder = 3;
2878
	(*grid)[k].elemmidpoints = TRUE;
2879
      }
2880

2881
      mode = 0;
2882
      break;
2883

2884
    case 10: /* COORDINATE RATIO */
2885
      if((*grid)[k].dimension == 2) 
2886
	sscanf(line,"%le",&(*grid)[k].xyratio);
2887
      if((*grid)[k].dimension == 3) 
2888
	sscanf(line,"%le %le",&(*grid)[k].xyratio,&(*grid)[k].xzratio);      
2889
      mode = 0;
2890
      break;
2891

2892
    case 11: /* MATERIALS */
2893
      sscanf(line,"%d %d",&(*grid)[k].firstmaterial,&(*grid)[k].lastmaterial);      
2894
      mode = 0;
2895
      break;
2896

2897
    case 12: /* LAYERES */
2898
      for(i=1;i<=(*grid)[k].zcells;i++) {
2899
	Getline(line,in);
2900
	sscanf(line,"%d %d %d\n",
2901
		&(*grid)[k].zfirstmaterial[i],&(*grid)[k].zlastmaterial[i],&(*grid)[k].zmaterial[i]); 
2902
      }
2903
      mode = 0;
2904
      break;
2905

2906
    case 13: /* ELEMENT RATIOS */
2907
      printf("Loading element ratios\n");
2908

2909
      for (j=1;j<=(*grid)[k].dimension;j++) {
2910
	Getline(line,in);
2911
	cp = line;
2912

2913
	if(j==1) for(i=1;i<=(*grid)[k].xcells;i++) (*grid)[k].xexpand[i] = next_real(&cp);
2914
	if(j==2) for(i=1;i<=(*grid)[k].ycells;i++) (*grid)[k].yexpand[i] = next_real(&cp);
2915
	if(j==3) for(i=1;i<=(*grid)[k].zcells;i++) (*grid)[k].zexpand[i] = next_real(&cp);
2916
      }
2917
      mode = 0;
2918
      break;
2919

2920
    case 29: /* ELEMENT NUMBER */
2921
      printf("Loading element numbers\n");
2922

2923
      for (j=1;j<=(*grid)[k].dimension;j++) {
2924
	Getline(line,in);
2925
	cp = line;
2926
	if(j==1) for(i=1;i<=(*grid)[k].xcells;i++) (*grid)[k].xelems[i] = next_int(&cp);
2927
	if(j==2) for(i=1;i<=(*grid)[k].ycells;i++) (*grid)[k].yelems[i] = next_int(&cp);
2928
	if(j==3) for(i=1;i<=(*grid)[k].zcells;i++) (*grid)[k].zelems[i] = next_int(&cp);
2929
      }
2930
      (*grid)[k].autoratio = 0;
2931
      mode = 0;
2932
      break;
2933

2934
    case 27: /* ELEMENT MINIMUM */
2935
      printf("Loading minimum number of elements\n");
2936
      if((*grid)[k].dimension == 1) 
2937
	sscanf(line,"%d",&(*grid)[k].minxelems);
2938
      if((*grid)[k].dimension == 2) 
2939
	sscanf(line,"%d %d",&(*grid)[k].minxelems,&(*grid)[k].minyelems);
2940
      if((*grid)[k].dimension == 3) 
2941
	sscanf(line,"%d %d %d",&(*grid)[k].minxelems,&(*grid)[k].minyelems,&(*grid)[k].minzelems);
2942
      mode = 0;
2943
      break;
2944

2945
    case 14: /* ELEMENT DENSITIES */
2946
      printf("Loading element densities\n");
2947
      for (j=1;j<=(*grid)[k].dimension;j++) {
2948
	Getline(line,in);
2949
	cp = line;
2950

2951
	if(j==1) for(i=1;i<=(*grid)[k].xcells;i++) (*grid)[k].xdens[i] = next_real(&cp);
2952
	if(j==2) for(i=1;i<=(*grid)[k].ycells;i++) (*grid)[k].ydens[i] = next_real(&cp);
2953
	if(j==3) for(i=1;i<=(*grid)[k].zcells;i++) (*grid)[k].zdens[i] = next_real(&cp);
2954
      }
2955
      mode = 0;
2956
      break;
2957

2958
    case 15: /* BOUNDARY CONDITIONS */
2959
      sscanf(line,"%d",&(*grid)[k].noboundaries);
2960
      printf("Loading %d boundary conditions\n",(*grid)[k].noboundaries);
2961

2962
      for(i=0;i<(*grid)[k].noboundaries;i++) {
2963
	Getline(line,in);
2964
	sscanf(line,"%d %d %d %d",
2965
	       &(*grid)[k].boundtype[i],&(*grid)[k].boundext[i],
2966
	       &(*grid)[k].boundint[i],&(*grid)[k].boundsolid[i]);
2967
      }  
2968
      mode = 0;
2969
      break;
2970

2971
    case 20: /* ROTATE */
2972
      (*grid)[k].rotate = TRUE;
2973
      mode = 0;
2974
      break;
2975

2976
    case 21: /* ROTRAD */
2977
      sscanf(line,"%le",&(*grid)[k].rotateradius2);
2978
      mode = 0;
2979
      break;
2980

2981
    case 22: /* ROTBLOCK */
2982
      sscanf(line,"%d",&(*grid)[k].rotateblocks);
2983
      if(0) printf("Reading blocks %d\n",(*grid)[k].rotateblocks);
2984
      mode = 0;
2985
      break;
2986

2987
    case 24: /* ROTIMP */
2988
      sscanf(line,"%le",&(*grid)[k].rotateimprove);
2989
      mode = 0;
2990
      break;
2991

2992
    case 30: /* POLAR RADIUS */
2993
      sscanf(line,"%le",&(*grid)[k].polarradius);
2994
      mode = 0;
2995
      break;
2996

2997
    case 25: /* ROTCURVE */
2998
      (*grid)[k].rotatecurve = TRUE;
2999
      sscanf(line,"%le%le%le",&(*grid)[k].curvezet,
3000
	     &(*grid)[k].curverad,&(*grid)[k].curveangle);
3001
      mode = 0;
3002
      break;
3003

3004
    case 26: /* REDUCE ELEMENT */
3005
      sscanf(line,"%d%d",&(*grid)[k].reduceordermatmin,
3006
	     &(*grid)[k].reduceordermatmax);
3007
      mode = 0;
3008
      break;
3009

3010
    case 28: /* LAYERED BO */
3011
      sscanf(line,"%d",&(*grid)[k].layeredbc);
3012
      mode = 0;
3013
      break;
3014

3015
    case 23: /* SCALING */
3016
      sscanf(line,"%le",&scaling);
3017
      for(i=0;i<=(*grid)[k].xcells;i++) (*grid)[k].x[i] *= scaling;
3018
      if((*grid)[k].dimension > 1) 
3019
	for(i=0;i<=(*grid)[k].ycells;i++) (*grid)[k].y[i] *= scaling;
3020
      if((*grid)[k].dimension == 3) 
3021
	for(i=0;i<=(*grid)[k].ycells;i++) (*grid)[k].z[i] *= scaling;
3022

3023
      (*grid)[k].rotateradius2 *= scaling;
3024
      (*grid)[k].curverad *= scaling;
3025
      (*grid)[k].curvezet *= scaling;
3026
      mode = 0;
3027
      break;
3028

3029
    default:
3030
      if(0) printf("Unknown case: %s",line);
3031
    }
3032

3033
  }
3034

3035
end:
3036

3037
  if(info) printf("Found %d divisions for grid\n",*nogrids);
3038

3039
  for(k=nogrids0;k < (*nogrids) && k<MAXCASES;k++) {
3040
    SetElementDivision(&(*grid)[k],1.0,info);
3041
  }
3042

3043

3044
  fclose(in);
3045
  return(error);
3046
}
3047

3048

3049

3050
int LoadElmergrid(struct GridType **grid,int *nogrids,char *prefix,int info) 
3051
{
3052
  char filename[MAXFILESIZE];
3053
  char command[MAXLINESIZE],params[MAXLINESIZE];
3054
  FILE *in;
3055
  int i,j,k;
3056
  char *cp;
3057
  int noknots,noelements,dim,axisymmetric;
3058
  int elemcode,maxnodes,totelems,nogrids0,minmat,maxmat;
3059
  long code;
3060
  Real raid;
3061

3062
  AddExtension(prefix,filename,"grd");
3063
  if ((in = fopen(filename,"r")) == NULL) {
3064
    printf("LoadElmergrid: opening of the file '%s' wasn't successful !\n",filename);
3065
    return(1);
3066
  }
3067

3068
  if(info) printf("Loading the geometry from file '%s'.\n",filename);
3069

3070
  InitGrid(grid[*nogrids]);
3071
  k = *nogrids;
3072
  nogrids0 = *nogrids;
3073

3074
  iodebug = FALSE;
3075

3076
  noknots = 0;
3077
  noelements = 0;
3078
  dim = 0;
3079
  axisymmetric = FALSE;
3080
  elemcode = 0;
3081
  maxnodes = 4;
3082
  totelems = 0;
3083

3084
  matcactive = FALSE;
3085

3086
  for(;;) {
3087
    if(GetCommand(command,params,in)) {
3088
      if(info) printf("Reached the end of command file\n");
3089
      goto end;
3090
    }    
3091

3092
    /* Control information */
3093
    if(strstr(command,"VERSION")) {
3094
      if(strstr(command,"080500")) {
3095
	if(info) printf("Loading old version of Elmergrid file.\n");
3096
	i = LoadElmergridOld(grid,nogrids,prefix,info);
3097
	return(i);
3098
      }
3099
      else {
3100
	sscanf(params,"%ld",&code);
3101
	if(code == 210903) 
3102
	  if(info) printf("Loading ElmerGrid file version: %ld\n",code);
3103
	else {
3104
	  printf("Unknown ElmerGrid file version: %ld\n",code);
3105
	  return(2);
3106
	}
3107
      }
3108
      *nogrids += 1;
3109
    }      
3110

3111
    else if(strstr(command,"DEBUG IO")) {
3112
      for(i=0;i<MAXLINESIZE;i++) params[i] = toupper(params[i]);
3113
      if(strstr(params,"FALSE")) 
3114
	iodebug = FALSE;
3115
      else {
3116
	iodebug = TRUE;
3117
	printf("IO debugging activated\n");
3118
      }
3119
    }
3120
 
3121
    else if(strstr(command,"MATC")) {
3122
      for(i=0;i<MAXLINESIZE;i++) params[i] = toupper(params[i]);
3123
      if(strstr(params,"FALSE")) 
3124
	matcactive = FALSE;
3125
      else {
3126
#if USE_MATC
3127
	matcactive = TRUE;
3128
	mtc_init(NULL, stdout, stderr);
3129
	strcpy(command, "format( 12 )");	
3130
	mtc_domath(command);	 
3131
	printf("MATC language activated with 12 digit accuracy.\n");	
3132
#else
3133
        matcactive = FALSE;
3134
        printf("Unable to enable MATC, as it is not even compiled.\n");
3135
        bigerror("Unable to enable MATC, as it is not even compiled.");
3136
#endif 
3137
      }
3138
    }
3139
    
3140
    else if(strstr(command,"COORDINATE SYSTEM")) {
3141
      for(i=0;i<MAXLINESIZE;i++) params[i] = toupper(params[i]);
3142
      grid[k]->dimension = 2;
3143
      if(strstr(params,"CARTESIAN 1D")) {
3144
	grid[k]->coordsystem = COORD_CART1;
3145
	grid[k]->dimension = 1;
3146
      }
3147
      else if(strstr(params,"CARTESIAN 2D")) 
3148
	grid[k]->coordsystem = COORD_CART2;
3149
      else if(strstr(params,"AXISYMMETRIC")) 
3150
	grid[k]->coordsystem = COORD_AXIS;
3151
      else if(strstr(params,"POLAR"))
3152
	grid[k]->coordsystem = COORD_POLAR;
3153
      else if(strstr(params,"CARTESIAN 3D")) {
3154
	grid[k]->coordsystem = COORD_CART3;
3155
	grid[k]->dimension = 3;
3156
      }
3157
      else printf("Unknown coordinate system: %s\n",params);
3158
      if(info) printf("Defining the coordinate system (%d-DIM).\n",grid[k]->dimension);
3159
    }
3160
    
3161
    else if(strstr(command,"SUBCELL DIVISIONS")) {
3162
      if(grid[k]->dimension == 1) {
3163
	sscanf(params,"%d",&(*grid)[k].xcells);
3164
	grid[k]->ycells = 1;	
3165
      }
3166
      else if(grid[k]->dimension == 2) 
3167
	sscanf(params,"%d %d",&(*grid)[k].xcells,&(*grid)[k].ycells);
3168
      else if(grid[k]->dimension == 3) 
3169
	sscanf(params,"%d %d %d",&(*grid)[k].xcells,&(*grid)[k].ycells,&(*grid)[k].zcells);      
3170

3171
      if(grid[k]->xcells >= MAXCELLS || grid[k]->ycells >= MAXCELLS || grid[k]->zcells >= MAXCELLS) {
3172
	printf("LoadElmergrid: Too many subcells [%d %d %d] vs. %d:\n",
3173
	       grid[k]->xcells,grid[k]->ycells,grid[k]->zcells,MAXCELLS);
3174
     }
3175

3176
      /* Initialize the default structure with ones */
3177
      for(j=grid[k]->ycells;j>=1;j--) 
3178
	for(i=1;i<=grid[k]->xcells;i++) 
3179
	  grid[k]->structure[j][i] = 1;
3180
    }
3181
    
3182
    else if(strstr(command,"MINIMUM ELEMENT DIVISION")) {
3183
      if(info) printf("Loading minimum number of elements\n");
3184

3185
      if((*grid)[k].dimension == 1) 
3186
	sscanf(params,"%d",&(*grid)[k].minxelems); 
3187

3188
      if((*grid)[k].dimension == 2)  
3189
	sscanf(params,"%d %d",&(*grid)[k].minxelems,&(*grid)[k].minyelems);
3190

3191
      if((*grid)[k].dimension == 3) 
3192
	sscanf(params,"%d %d %d",&(*grid)[k].minxelems,
3193
	       &(*grid)[k].minyelems,&(*grid)[k].minzelems);
3194
    }      
3195
    
3196
    else if(strstr(command,"SUBCELL LIMITS 1")) {
3197
      if(info) printf("Loading %d subcell limits in X-direction\n",grid[k]->xcells+1);
3198
      cp = params;
3199
      for(i=0;i<=grid[k]->xcells;i++) {
3200
	grid[k]->x[i] = next_real(&cp);
3201
	if(i > 0 && grid[k]->x[i] < grid[k]->x[i-1]) {
3202
	  printf("Subcell limits 1(%d): %12.6le %12.6le\n",i,grid[k]->x[i],grid[k]->x[i-1]);
3203
	  bigerror("Values for limits 1 should be a growing series, existing\n");
3204
	}
3205
      }
3206
    }    
3207
    else if(strstr(command,"SUBCELL LIMITS 2")) {
3208
      if(info) printf("Loading %d subcell limits in Y-direction\n",grid[k]->ycells+1);
3209
      cp = params;
3210
      for(i=0;i<=grid[k]->ycells;i++) {
3211
	grid[k]->y[i] = next_real(&cp);
3212
	if(i > 0 && grid[k]->y[i] < grid[k]->y[i-1]) {
3213
	  printf("Subcell limits 2(%d): %12.6le %12.6le\n",i,grid[k]->y[i],grid[k]->y[i-1]);
3214
	  bigerror("Values for limits should be a growing series, existing\n");
3215
	}
3216
      }
3217
    }      
3218
    else if(strstr(command,"SUBCELL LIMITS 3")) {
3219
      if(info) printf("Loading %d subcell limits in Z-direction\n",grid[k]->zcells+1);
3220
      cp = params;
3221
      for(i=0;i<=grid[k]->zcells;i++) {
3222
	grid[k]->z[i] = next_real(&cp);
3223
	if(i > 0 && grid[k]->z[i] < grid[k]->z[i-1]) {
3224
	  printf("Subcell limits 3(%d): %12.6le %12.6le\n",i,grid[k]->z[i],grid[k]->z[i-1]);
3225
	  bigerror("Values for limits should be a growing series, existing\n");
3226
	}
3227
      }
3228
    }
3229

3230
    else if(strstr(command,"SUBCELL SIZES 1")) {
3231
      if(info) printf("Loading %d subcell sizes in X-direction\n",grid[k]->xcells);
3232
      cp = params;
3233
      for(i=1;i<=grid[k]->xcells;i++) grid[k]->x[i] = next_real(&cp);
3234
      for(i=1;i<=grid[k]->xcells;i++) grid[k]->x[i] = grid[k]->x[i-1] + grid[k]->x[i];
3235
    }      
3236
    else if(strstr(command,"SUBCELL SIZES 2")) {
3237
      if(info) printf("Loading %d subcell sizes in Y-direction\n",grid[k]->ycells);
3238
      cp = params;
3239
      for(i=1;i<=grid[k]->ycells;i++) grid[k]->y[i] = next_real(&cp);
3240
      for(i=1;i<=grid[k]->ycells;i++) grid[k]->y[i] = grid[k]->y[i-1] + grid[k]->y[i];
3241
    }      
3242
    else if(strstr(command,"SUBCELL SIZES 3")) {
3243
      if(info) printf("Loading %d subcell sizes in Z-direction\n",grid[k]->zcells);
3244
      cp = params;
3245
      for(i=1;i<=grid[k]->zcells;i++) grid[k]->z[i] = next_real(&cp);
3246
      for(i=1;i<=grid[k]->zcells;i++) grid[k]->z[i] = grid[k]->z[i-1] + grid[k]->z[i];
3247
    }
3248

3249
    else if(strstr(command,"SUBCELL ORIGIN 1")) {
3250
      for(i=0;i<MAXLINESIZE;i++) params[i] = toupper(params[i]);
3251
      if(strstr(params,"CENTER")) {
3252
	raid = 0.5 * (grid[k]->x[0] + grid[k]->x[grid[k]->xcells]);
3253
      }
3254
      else if(strstr(params,"LEFT") || strstr(params,"MIN") ) {
3255
	raid = grid[k]->x[0];
3256
      }
3257
      else if(strstr(params,"RIGHT") || strstr(params,"MAX") ) {
3258
	raid = grid[k]->x[grid[k]->xcells];
3259
      }
3260
      else {
3261
	cp = params;
3262
	raid = next_real(&cp);
3263
      }
3264
      for(i=0;i<=grid[k]->xcells;i++) grid[k]->x[i] -= raid;
3265
    }
3266
    else if(strstr(command,"SUBCELL ORIGIN 2")) {
3267
      for(i=0;i<MAXLINESIZE;i++) params[i] = toupper(params[i]);
3268
      if(strstr(params,"CENTER")) {
3269
	raid = 0.5 * (grid[k]->y[0] + grid[k]->y[grid[k]->ycells]);
3270
      }
3271
      else if(strstr(params,"LEFT")) {
3272
	raid = grid[k]->y[0];
3273
      }
3274
      else if(strstr(params,"RIGHT")) {
3275
	raid = grid[k]->y[grid[k]->ycells];
3276
      }
3277
      else {
3278
	cp = params;
3279
	raid = next_real(&cp);
3280
      }      
3281
      for(i=0;i<=grid[k]->ycells;i++) grid[k]->y[i] -= raid;
3282
    }
3283
    else if(strstr(command,"SUBCELL ORIGIN 3")) {
3284
      for(i=0;i<MAXLINESIZE;i++) params[i] = toupper(params[i]);
3285
      if(strstr(params,"CENTER")) {
3286
	raid = 0.5 * (grid[k]->z[0] + grid[k]->z[grid[k]->zcells]);
3287
      }
3288
      else if(strstr(params,"LEFT")) {
3289
	raid = grid[k]->z[0];
3290
      }
3291
      else if(strstr(params,"RIGHT")) {
3292
	raid = grid[k]->z[grid[k]->zcells];
3293
      }
3294
      else {
3295
	cp = params;
3296
	raid = next_real(&cp);
3297
      }
3298
      for(i=0;i<=grid[k]->zcells;i++) grid[k]->z[i] -= raid;      
3299
    }
3300

3301
    else if(strstr(command,"MATERIAL STRUCTURE")) {
3302
      if(info) printf("Loading material structure\n");
3303

3304
      /* Initialize the default structure with zeros */
3305
      for(j=grid[k]->ycells;j>=1;j--) 
3306
	for(i=1;i<=grid[k]->xcells;i++) 
3307
	  grid[k]->structure[j][i] = 0;
3308
     
3309
      for(j=grid[k]->ycells;j>=1;j--) {
3310
	if(j < grid[k]->ycells) Getline(params,in);
3311
	cp=params;
3312
	for(i=1;i<=grid[k]->xcells;i++) 
3313
	  grid[k]->structure[j][i] = next_int(&cp);
3314
      }      
3315
      minmat = maxmat = grid[k]->structure[1][1];
3316
      for(j=grid[k]->ycells;j>=1;j--) 
3317
	for(i=1;i<=grid[k]->xcells;i++) {
3318
	  if(minmat > grid[k]->structure[j][i])
3319
	    minmat = grid[k]->structure[j][i];
3320
	  if(maxmat < grid[k]->structure[j][i])
3321
	    maxmat = grid[k]->structure[j][i];
3322
	}      
3323
      if(minmat < 0) 
3324
	printf("LoadElmergrid: please use positive material indices.\n");
3325
      if(maxmat > MAXBODYID) 
3326
	printf("LoadElmergrid: material indices larger to %d may create problems.\n",
3327
	       MAXBODYID);
3328
      printf("LoadElmergrid: materials interval is [%d,%d]\n",minmat,maxmat);
3329

3330
      grid[k]->maxmaterial = maxmat;
3331
    }
3332
    else if(strstr(command,"MATERIALS INTERVAL")) { 
3333
      sscanf(params,"%d %d",&(*grid)[k].firstmaterial,&(*grid)[k].lastmaterial);      
3334
    }
3335
     
3336
    else if(strstr(command,"REVOLVE")) {
3337
      if(0) printf("revolve: %s %s\n",command,params);
3338

3339
      if(strstr(command,"REVOLVE RADIUS")) {
3340
	(*grid)[k].rotate = TRUE;
3341
	sscanf(params,"%le",&(*grid)[k].rotateradius2);
3342
      }
3343
      else if(strstr(command,"REVOLVE BLOCKS")) {
3344
	(*grid)[k].rotate = TRUE;
3345
	sscanf(params,"%d",&(*grid)[k].rotateblocks);
3346
      }
3347
      else if(strstr(command,"REVOLVE IMPROVE")) {
3348
	(*grid)[k].rotate = TRUE;
3349
	sscanf(params,"%le",&(*grid)[k].rotateimprove);
3350
      }
3351
      else if(strstr(command,"REVOLVE RADIUS")) {
3352
	sscanf(params,"%le",&(*grid)[k].polarradius);
3353
      }
3354
      else if(strstr(command,"REVOLVE CURVE DIRECT")) {
3355
	(*grid)[k].rotatecurve = TRUE;
3356
	sscanf(params,"%le",&(*grid)[k].curvezet);
3357
      }
3358
      else if(strstr(command,"REVOLVE CURVE RADIUS")) {
3359
	(*grid)[k].rotatecurve = TRUE;
3360
	sscanf(params,"%le",&(*grid)[k].curverad);
3361
      }
3362
      else if(strstr(command,"REVOLVE CURVE ANGLE")) {
3363
	(*grid)[k].rotatecurve = TRUE;
3364
	sscanf(params,"%le",&(*grid)[k].curveangle);
3365
      }
3366
    }
3367

3368
    else if(strstr(command,"REDUCE ORDER INTERVAL")) {
3369
      sscanf(params,"%d%d",&(*grid)[k].reduceordermatmin,
3370
	     &(*grid)[k].reduceordermatmax);
3371
    }
3372
    
3373
    else if(strstr(command,"BOUNDARY DEFINITION")) {
3374
      printf("Loading boundary conditions\n");
3375
      
3376
      for(i=0;i<MAXBOUNDARIES;i++) {
3377
	if(i>0) Getline(params,in);
3378
	for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
3379
	if(strstr(params,"END")) break;
3380
	sscanf(params,"%d %d %d %d",
3381
	       &(*grid)[k].boundtype[i],&(*grid)[k].boundext[i],
3382
	       &(*grid)[k].boundint[i],&(*grid)[k].boundsolid[i]);
3383
      }  
3384
      printf("Found %d boundaries\n",i);
3385
      (*grid)[k].noboundaries = i;
3386
    }
3387
    
3388
    /* else if(strstr(command,"LAYERED BOUNDARIES")) {
3389
      for(i=0;i<MAXLINESIZE;i++) params[i] = toupper(params[i]);
3390
      if(strstr(params,"TRUE")) (*grid)[k].layeredbc = 1;
3391
      if(strstr(params,"FALSE")) (*grid)[k].layeredbc = 0;
3392
    } */
3393
    
3394
    else if(strstr(command,"NUMBERING")) {
3395
      for(i=0;i<MAXLINESIZE;i++) params[i] = toupper(params[i]);
3396
      if(strstr(params,"HORIZONTAL")) (*grid)[k].numbering = NUMBER_XY;
3397
      if(strstr(params,"VERTICAL")) (*grid)[k].numbering = NUMBER_YX;
3398
    }
3399
    
3400
    else if(strstr(command,"ELEMENT DEGREE")) {
3401
      sscanf(params,"%d",&(*grid)[k].elemorder);
3402
    }
3403
    
3404
    else if(strstr(command,"ELEMENT INNERNODES")) {
3405
      for(i=0;i<MAXLINESIZE;i++) params[i] = toupper(params[i]);
3406
      if(strstr(params,"TRUE")) (*grid)[k].elemmidpoints = TRUE;
3407
      if(strstr(params,"FALSE")) (*grid)[k].elemmidpoints = FALSE;
3408
    }
3409
    else if(strstr(command,"ELEMENTTYPE") || strstr(command,"ELEMENTCODE")) {
3410
      sscanf(params,"%d",&elemcode);
3411
    }
3412
    
3413
    else if(strstr(command,"TRIANGLES CRITICAL ANGLE")) {
3414
      sscanf(params,"%le",&(*grid)[k].triangleangle);      
3415
    }
3416
    else if(strstr(command,"TRIANGLES")) {
3417
      for(i=0;i<MAXLINESIZE;i++) params[i] = toupper(params[i]);
3418
      if(strstr(params,"TRUE")) (*grid)[k].triangles = TRUE;
3419
      if(strstr(params,"FALSE")) (*grid)[k].triangles = FALSE;
3420
    }
3421
    
3422
    else if(strstr(command,"PLANE ELEMENTS")) {
3423
      sscanf(params,"%d",&(*grid)[k].wantedelems);
3424
    }
3425
    else if(strstr(command,"SURFACE ELEMENTS")) {
3426
      sscanf(params,"%d",&(*grid)[k].wantedelems);
3427
    }
3428
    else if(strstr(command,"REFERENCE DENSITY")) {
3429
      sscanf(params,"%le",&(*grid)[k].limitdx);
3430
      (*grid)[k].autoratio = 3;     
3431
    }
3432
    else if(strstr(command,"VERIFY DENSITY")) {
3433
      for(i=0;i<MAXLINESIZE;i++) params[i] = toupper(params[i]);
3434
      if(strstr(params,"TRUE")) (*grid)[k].limitdxverify = TRUE;
3435
      if(strstr(params,"FALSE")) (*grid)[k].limitdxverify = FALSE;
3436
    }
3437
    else if(strstr(command,"VOLUME ELEMENTS")) {
3438
      sscanf(params,"%d",&(*grid)[k].wantedelems3d);
3439
    }
3440
    else if(strstr(command,"VOLUME NODES")) {
3441
      sscanf(params,"%d",&(*grid)[k].wantednodes3d);
3442
    }
3443

3444
    else if(strstr(command,"COORDINATE RATIO")) {
3445
      if((*grid)[k].dimension == 2) 
3446
	sscanf(params,"%le",&(*grid)[k].xyratio);
3447
      if((*grid)[k].dimension == 3) 
3448
	sscanf(params,"%le %le",&(*grid)[k].xyratio,&(*grid)[k].xzratio);      
3449
    }
3450
    
3451
    else if(strstr(command,"ELEMENT RATIOS 1")) {
3452
      cp = params;
3453
      for(i=1;i<=(*grid)[k].xcells;i++) (*grid)[k].xexpand[i] = next_real(&cp);
3454
    }
3455
    else if(strstr(command,"ELEMENT RATIOS 2")) {
3456
      cp = params;
3457
      for(i=1;i<=(*grid)[k].ycells;i++) (*grid)[k].yexpand[i] = next_real(&cp);
3458
    }
3459
    else if(strstr(command,"ELEMENT RATIOS 3")) {
3460
      cp = params;
3461
      for(i=1;i<=(*grid)[k].zcells;i++) (*grid)[k].zexpand[i] = next_real(&cp);
3462
    }
3463
    
3464
    else if(strstr(command,"ELEMENT DENSITIES 1")) {
3465
      cp = params;
3466
      for(i=1;i<=(*grid)[k].xcells;i++) (*grid)[k].xdens[i] = next_real(&cp);
3467
    }
3468
    else if(strstr(command,"ELEMENT DENSITIES 2")) {
3469
      cp = params;
3470
      for(i=1;i<=(*grid)[k].ycells;i++) (*grid)[k].ydens[i] = next_real(&cp);
3471
    }
3472
    else if(strstr(command,"ELEMENT DENSITIES 3")) {
3473
      cp = params;
3474
      for(i=1;i<=(*grid)[k].zcells;i++) (*grid)[k].zdens[i] = next_real(&cp);
3475
    }
3476
    
3477
    else if(strstr(command,"ELEMENT DIVISIONS 1")) {
3478
      cp = params;
3479
      for(i=1;i<=(*grid)[k].xcells;i++) (*grid)[k].xelems[i] = next_int(&cp);
3480
      (*grid)[k].autoratio = 0;
3481
    }
3482
    else if(strstr(command,"ELEMENT DIVISIONS 2")) {
3483
      cp = params;
3484
      for(i=1;i<=(*grid)[k].ycells;i++) (*grid)[k].yelems[i] = next_int(&cp);
3485
      (*grid)[k].autoratio = 0;
3486
    }
3487
    else if(strstr(command,"ELEMENT DIVISIONS 3")) {
3488
      cp = params;
3489
      for(i=1;i<=(*grid)[k].zcells;i++) (*grid)[k].zelems[i] = next_int(&cp);
3490
      (*grid)[k].autoratio = 0;
3491
    }
3492
    
3493
    /* else if(strstr(command,"EXTRUDED STRUCTURE")) {
3494
      for(i=1;i<=(*grid)[k].zcells;i++) {
3495
	if(i>1) Getline(params,in);
3496
	sscanf(params,"%d %d %d\n",
3497
	       &(*grid)[k].zfirstmaterial[i],&(*grid)[k].zlastmaterial[i],&(*grid)[k].zmaterial[i]); 
3498
      }
3499
    } */
3500
    
3501
    else if(strstr(command,"GEOMETRY MAPPINGS")) {     
3502
      if(k > 0) (*grid)[k].mappings = 0;
3503

3504
      for(i=0;i<MAXLINESIZE;i++) params[i] = toupper(params[i]);
3505
      for(i=(*grid)[k].mappings;i<MAXMAPPINGS;i++) {
3506
	if(i>(*grid)[k].mappings) Getline(params,in);
3507

3508
	if(strstr(params,"END")) break;
3509
	cp=params; 
3510
	(*grid)[k].mappingtype[i] = next_int(&cp);	
3511
#if 0
3512
	(*grid)[k].mappingtype[i] += 50*SGN((*grid)[k].mappingtype[i]);
3513
#endif
3514
	(*grid)[k].mappingline[i] = next_int(&cp);
3515
	(*grid)[k].mappinglimits[2*i] = next_real(&cp);
3516
	(*grid)[k].mappinglimits[2*i+1] = next_real(&cp);
3517
	(*grid)[k].mappingpoints[i] = next_int(&cp);
3518
	(*grid)[k].mappingparams[i] = Rvector(0,(*grid)[k].mappingpoints[i]);
3519
	for(j=0;j<(*grid)[k].mappingpoints[i];j++) 
3520
	  (*grid)[k].mappingparams[i][j] = next_real(&cp);
3521
      }      
3522
      printf("Loaded %d geometry mappings\n",i);
3523
      (*grid)[k].mappings = i;      
3524
    }
3525

3526
    else if(strstr(command,"END") ) {      
3527
      if(0) printf("End of field\n");
3528
    }
3529
      
3530
    else if(strstr(command,"START NEW MESH")) {
3531
      if((*nogrids) >= MAXCASES) {
3532
	printf("There are more grids than was allocated for!\n"); 
3533
	printf("Ignoring meshes starting from %d\n.",(*nogrids)+1);
3534
	goto end;
3535
      }
3536
      (*nogrids)++;
3537
      printf("\nLoading element meshing no %d\n",*nogrids);
3538
      k = *nogrids - 1;	           
3539
      if(k > nogrids0) (*grid)[k] = (*grid)[k-1];	 
3540
    }
3541

3542
    else {
3543
      if(0) printf("Unknown command: %s",command);
3544
    }
3545
  }
3546

3547
end:
3548

3549
  if(info) printf("Found %d divisions for grid\n",*nogrids);
3550
  
3551
  for(k=nogrids0;k < (*nogrids) && k<MAXCASES;k++) {
3552

3553
    if(elemcode == 0) {
3554
      if((*grid)[k].dimension == 1) {
3555
	(*grid)[k].nonodes = (*grid)[k].elemorder + 1;
3556
      }
3557
      else if((*grid)[k].elemmidpoints == FALSE) {
3558
	(*grid)[k].nonodes = 4 * (*grid)[k].elemorder;
3559
      }					
3560
      else {
3561
	if((*grid)[k].elemorder == 2) (*grid)[k].nonodes = 9;
3562
	if((*grid)[k].elemorder == 3) (*grid)[k].nonodes = 16;	
3563
      }
3564
    }
3565
    else if(elemcode/100 == 2) {
3566
      (*grid)[k].triangles = FALSE;      
3567
      (*grid)[k].nonodes = elemcode%100;
3568
    }
3569
    else if(elemcode/100 == 4) {
3570
      (*grid)[k].triangles = FALSE;      
3571
      (*grid)[k].nonodes = elemcode%100;
3572
    }
3573
    else if(elemcode/100 == 3) {  
3574
      (*grid)[k].triangles = TRUE;      
3575
      if(elemcode%100 == 3)       (*grid)[k].nonodes = 4;
3576
      else if(elemcode%100 == 6)  (*grid)[k].nonodes = 9;
3577
      else if(elemcode%100 == 10) (*grid)[k].nonodes = 16;	
3578
    }    
3579
  }
3580

3581
  
3582
  /* for(k=nogrids0;k < (*nogrids) && k<MAXCASES;k++) {
3583
    SetElementDivision(&(*grid)[k],relh,info);
3584
    } */
3585

3586
  fclose(in);
3587
  return(0);
3588
}
3589

3590

3591

3592
void InitParameters(struct ElmergridType *eg)
3593
{
3594
  int i;
3595

3596
  eg->relh = 1.0;
3597
  eg->inmethod = 0;
3598
  eg->outmethod = 0;
3599
  eg->silent = FALSE;
3600
  eg->binary = 0;
3601
  eg->nofilesin = 1;
3602
  eg->unitemeshes = FALSE;
3603
  eg->unitenooverlap = FALSE;
3604
  eg->triangles = FALSE;
3605
  eg->triangleangle = 0.0;
3606
  eg->rotate = FALSE;
3607
  eg->polar = FALSE;
3608
  eg->cylinder = FALSE;
3609
  eg->usenames = TRUE;
3610
  eg->layers = 0;
3611
  eg->layereps = 0.0;
3612
  eg->layermove = 0;
3613
  eg->partitions = 0;
3614
  eg->elements3d = 0;
3615
  eg->nodes3d = 0;
3616
  eg->metis = 0;
3617
  eg->metis_contig = FALSE;
3618
  eg->metis_volcut = FALSE;
3619
  eg->metis_seed = 0;
3620
  eg->metis_ncuts = 1;
3621
  eg->metis_minconn = FALSE;
3622
  eg->partopt = 0;
3623
  eg->partoptim = FALSE;
3624
  eg->partbcoptim = TRUE;
3625
  eg->partjoin = 0;
3626
  for(i=0;i<MAXHALOMODES;i++) eg->parthalo[i] = FALSE;
3627
  eg->reduce = FALSE;
3628
  eg->increase = FALSE;
3629
  eg->translate = FALSE;
3630
  eg->isoparam = FALSE;
3631
  eg->multidim = FALSE;
3632
  eg->removelowdim = FALSE;
3633
  eg->removeintbcs = FALSE;
3634
  eg->removeunused = FALSE;
3635
  eg->dim = 3;
3636
  eg->center = FALSE;
3637
  eg->scale = FALSE;
3638
  eg->order = FALSE;
3639
  eg->boundbounds = 0;
3640
  eg->saveinterval[0] = eg->saveinterval[1] = eg->saveinterval[2] = 0;
3641
  eg->bulkbounds = 0;
3642
  eg->partorder = FALSE;
3643
  eg->findsides = FALSE;
3644
  eg->parthypre = FALSE;
3645
  eg->partdual = FALSE;
3646
  eg->partbcz = 0;
3647
  eg->partbcr = 0;
3648
  eg->partbclayers = 1;
3649
  eg->partbcmetis = 0;
3650
  eg->partbw = FALSE;
3651
  eg->saveboundaries = TRUE;
3652
  eg->vtuone = FALSE;
3653
  eg->timeron = FALSE;
3654
  eg->nosave = FALSE;
3655
  eg->nooverwrite = FALSE;
3656
  eg->merge = FALSE;
3657
  eg->bcoffset = FALSE;
3658
  eg->periodic = 0;
3659
  eg->periodicdim[0] = 0;
3660
  eg->periodicdim[1] = 0;
3661
  eg->periodicdim[2] = 0;
3662
  eg->bulkorder = FALSE;
3663
  eg->boundorder = FALSE;
3664
  eg->sidemappings = 0;
3665
  eg->bulkmappings = 0;
3666
  eg->coordinatemap[0] = eg->coordinatemap[1] = eg->coordinatemap[2] = 0;
3667
  eg->clone[0] = eg->clone[1] = eg->clone[2] = 0;
3668
  eg->clonesize[0] = eg->clonesize[1] = eg->clonesize[2] = 0.0;
3669
  eg->mirror[0] = eg->mirror[1] = eg->mirror[2] = 0;
3670
  eg->cloneinds = FALSE;
3671
  eg->mirrorbc = 0;
3672
  eg->decimals = 12;
3673
  eg->discont = 0;
3674
  eg->connect = 0;
3675
  eg->connectboundsnosets = 0;
3676
  
3677
  eg->rotatecurve = FALSE;
3678
  eg->curverad = 0.5;
3679
  eg->curveangle = 90.0;
3680
  eg->curvezet = 0.0;
3681
  eg->parttol = 0.0;
3682
  eg->filerenamed = FALSE;
3683
  
3684
  for(i=0;i<MAXSIDEBULK;i++) 
3685
    eg->sidebulk[i] = 0;
3686
}
3687

3688

3689

3690

3691
int InlineParameters(struct ElmergridType *eg,int argc,char *argv[],int first,int info)
3692
{
3693
  int arg,i,dim;
3694
  char command[MAXLINESIZE];
3695
  
3696
  dim = eg->dim;
3697

3698
  printf("Elmergrid reading in-line arguments\n");
3699

3700
  /* Type of input file */
3701
  if(first > 3) {
3702
    for(i=0;i<MAXLINESIZE;i++) command[i] = ' ';
3703

3704
    strcpy(command,argv[1]);
3705
    for(i=0;i<MAXLINESIZE;i++) command[i] = toupper(command[i]);
3706
    for(i=0;i<=MAXINMETHODS;i++) {
3707
      if(strstr(command,InMethods[i])) {
3708
	eg->inmethod = i;
3709
	break;
3710
      }
3711
    }
3712
    if(i>MAXINMETHODS) eg->inmethod = atoi(argv[1]);
3713
    
3714

3715
    /* Type of output file (fewer options) */
3716
    strcpy(command,argv[2]);
3717
    for(i=0;i<MAXLINESIZE;i++) command[i] = toupper(command[i]);
3718
    for(i=1;i<=MAXOUTMETHODS;i++) {
3719
      if(strstr(command,OutMethods[i])) {
3720
	eg->outmethod = i;
3721
	break;
3722
      }
3723
    }
3724
    if(i>MAXOUTMETHODS) eg->outmethod = atoi(argv[2]);
3725
 
3726
    /* Name of output file */
3727
    strcpy(eg->filesin[0],argv[3]);
3728
    strcpy(eg->filesout[0],eg->filesin[0]);
3729
    strcpy(eg->infofile,eg->filesin[0]);
3730
  }
3731

3732

3733
  /* The optional inline parameters */
3734

3735
  for(arg=first;arg <argc; arg++) {
3736

3737
    if(strcmp(argv[arg],"-silent") == 0) {
3738
      eg->silent = TRUE;
3739
      info = FALSE;
3740
    }
3741
    else if(strcmp(argv[arg],"-verbose") == 0) {
3742
      eg->silent = FALSE;
3743
      info = TRUE;
3744
    }   
3745
    else if(strcmp(argv[arg],"-in") ==0 ) {
3746
      if(arg+1 >= argc) {
3747
	printf("The secondary input file name is required as a parameter\n");
3748
	return(1);
3749
      }
3750
      else {
3751
	strcpy(eg->filesin[eg->nofilesin],argv[arg+1]);
3752
	printf("A secondary input file %s will be loaded.\n",eg->filesin[eg->nofilesin]);
3753
	eg->nofilesin++;
3754
      }
3755
    }
3756
    else if(strcmp(argv[arg],"-out") == 0) {
3757
      if(arg+1 >= argc) {
3758
	printf("The output name is required as a parameter\n");
3759
	return(2);
3760
      }
3761
      else {
3762
	strcpy(eg->filesout[0],argv[arg+1]);
3763
	eg->filerenamed = TRUE;
3764
      }
3765
    }
3766
    else if(strcmp(argv[arg],"-bin") == 0) {
3767
      eg->binary = 1;
3768
    }
3769
    else if(strcmp(argv[arg],"-sbin") == 0) {
3770
      eg->binary = 2;
3771
    }
3772
    else if(strcmp(argv[arg],"-decimals") == 0) {
3773
      eg->decimals = atoi(argv[arg+1]);
3774
    }
3775
    else if(strcmp(argv[arg],"-triangles") ==0) {
3776
      eg->triangles = TRUE;
3777
      printf("The rectangles will be split to triangles.\n");
3778
      if(arg+1 < argc) {
3779
	if(strcmp(argv[arg+1],"-")) {
3780
	  eg->triangleangle = atof(argv[arg+1]);
3781
	}
3782
      }
3783
    }
3784
    else if(strcmp(argv[arg],"-merge") == 0) {
3785
      if(arg+1 >= argc) {
3786
	printf("Give a parameter for critical distance.\n");
3787
	return(3);
3788
      }
3789
      else {
3790
	eg->merge = TRUE;
3791
	eg->cmerge = atof(argv[arg+1]);
3792
      }
3793
    }
3794
    else if(strcmp(argv[arg],"-relh") == 0) {
3795
      if(arg+1 >= argc) {
3796
	printf("Give a relative mesh density related to the specifications\n");
3797
	return(3);
3798
      }
3799
      else {
3800
	eg->relh = atof(argv[arg+1]);
3801
      }
3802
    }
3803
    else if(strcmp(argv[arg],"-order") == 0) {
3804
      if(arg+dim >= argc) {
3805
	printf("Give %d parameters for the order vector.\n",dim);
3806
 	return(4);
3807
      }
3808
      else {
3809
	eg->order = TRUE;
3810
	eg->corder[0] = atof(argv[arg+1]);
3811
	eg->corder[1] = atof(argv[arg+2]);
3812
	if(dim==3) eg->corder[2] = atof(argv[arg+3]);
3813
      }
3814
    }
3815
    else if(strcmp(argv[arg],"-parttol") == 0) {
3816
      if(arg+1 >= argc) {
3817
	printf("Give a tolerance for geometric partition algorithms\n");
3818
	return(3);
3819
      }
3820
      else {
3821
	eg->parttol = atof(argv[arg+1]);
3822
      }
3823
    }
3824
    else if(strcmp(argv[arg],"-autoorder") == 0) {
3825
      eg->order = 2;
3826
    }
3827
    else if(strcmp(argv[arg],"-halo") == 0) {
3828
      eg->parthalo[1] = TRUE;
3829
    }
3830
    else if(strcmp(argv[arg],"-halobc") == 0) {
3831
      eg->parthalo[2] = TRUE;
3832
    }
3833
    else if(strcmp(argv[arg],"-halodb") == 0) {
3834
      eg->parthalo[1] = TRUE;
3835
      eg->parthalo[2] = TRUE;            
3836
    }   
3837
    else if(strcmp(argv[arg],"-haloz") == 0) {
3838
      eg->parthalo[3] = TRUE;
3839
    }
3840
    else if(strcmp(argv[arg],"-halor") == 0) {
3841
      eg->parthalo[3] = TRUE;
3842
    }
3843
    else if(strcmp(argv[arg],"-halogreedy") == 0) {
3844
      eg->parthalo[4] = TRUE;
3845
    }    
3846
    else if(strcmp(argv[arg],"-metisorder") == 0) {
3847
      eg->order = 3;
3848
    }    
3849
    else if(strcmp(argv[arg],"-centralize") == 0) {
3850
      eg->center = TRUE;
3851
    }
3852
    else if(strcmp(argv[arg],"-scale") == 0) {
3853
      if(arg+dim >= argc) {
3854
	printf("Give %d parameters for the scaling.\n",dim);
3855
 	return(5);
3856
     }
3857
      else {
3858
	eg->scale = TRUE;
3859
	eg->cscale[0] = atof(argv[arg+1]);
3860
	eg->cscale[1] = atof(argv[arg+2]);
3861
	if(dim==3) eg->cscale[2] = atof(argv[arg+3]);
3862
      }
3863
    }
3864
    else if(strcmp(argv[arg],"-translate") == 0) {
3865
      if(arg+dim >= argc) {
3866
	printf("Give %d parameters for the translate vector.\n",dim);
3867
	return(6);
3868
      }
3869
      else {
3870
	eg->translate = TRUE;
3871
	eg->ctranslate[0] = atof(argv[arg+1]);
3872
	eg->ctranslate[1] = atof(argv[arg+2]);
3873
	if(dim==3) eg->ctranslate[2] = atof(argv[arg+3]);
3874
      }
3875
    }
3876
    else if(strcmp(argv[arg],"-saveinterval") == 0) {
3877
      if(arg+dim >= argc) {
3878
	printf("Give min, max and step for the interval.\n");
3879
	return(7);
3880
      }
3881
      else {
3882
	eg->saveinterval[0] = atoi(argv[arg+1]);
3883
	eg->saveinterval[1] = atoi(argv[arg+2]);
3884
	eg->saveinterval[2] = atoi(argv[arg+3]);
3885
      }
3886
    }
3887
    else if(strcmp(argv[arg],"-rotate") == 0 || strcmp(argv[arg],"-rotate") == 0) {
3888
      if(arg+dim >= argc) {
3889
	printf("Give three parameters for the rotation angles.\n");
3890
	return(8);
3891
      }
3892
      else {
3893
	eg->rotate = TRUE;
3894
	eg->crotate[0] = atof(argv[arg+1]);
3895
	eg->crotate[1] = atof(argv[arg+2]);
3896
	eg->crotate[2] = atof(argv[arg+3]);
3897
      }
3898
    }
3899
    else if(strcmp(argv[arg],"-clone") == 0) {
3900
      if(arg+dim >= argc) {
3901
	printf("Give the number of clones in each %d directions.\n",dim);
3902
 	return(9);
3903
     }
3904
      else {
3905
	eg->clone[0] = atoi(argv[arg+1]);
3906
	eg->clone[1] = atoi(argv[arg+2]);
3907
	if(dim == 3) eg->clone[2] = atoi(argv[arg+3]);
3908
      }
3909
    }
3910
    else if(strcmp(argv[arg],"-clonesize") == 0) {
3911
      if(arg+dim >= argc) {
3912
	printf("Give the clone size in each %d directions.\n",dim);
3913
 	return(10);
3914
      }
3915
      else {
3916
	eg->clonesize[0] = atof(argv[arg+1]);
3917
	eg->clonesize[1] = atof(argv[arg+2]);
3918
	if(dim == 3) eg->clonesize[2] = atof(argv[arg+3]);
3919
      }
3920
    }
3921
    else if(strcmp(argv[arg],"-cloneinds") == 0) {
3922
      eg->cloneinds = TRUE;
3923
    }
3924
    else if(strcmp(argv[arg],"-mirror") == 0) {
3925
      if(arg+dim >= argc) {
3926
	printf("Give the symmetry of the coordinate directions, eg. 1 1 0\n");
3927
      }
3928
      else {
3929
	eg->mirror[0] = atoi(argv[arg+1]);
3930
	eg->mirror[1] = atoi(argv[arg+2]);
3931
	if(dim == 3) eg->mirror[2] = atoi(argv[arg+3]);
3932
      }
3933
    }
3934
    else if(strcmp(argv[arg],"-mirrorbc") == 0) {
3935
      if(arg+1 >= argc) {
3936
	printf("Give the number of symmetry BC.\n");
3937
 	return(11);
3938
      }
3939
      else {
3940
	eg->mirrorbc = atoi(argv[arg+1]);
3941
      }
3942
    }
3943
    else if(strcmp(argv[arg],"-unite") == 0) {
3944
      eg->unitemeshes = TRUE;
3945
      printf("The meshes will be united.\n");
3946
    }   
3947
    else if(strcmp(argv[arg],"-unitenooverlap") == 0) {
3948
      eg->unitemeshes = TRUE;
3949
      eg->unitenooverlap = TRUE;
3950
      printf("The meshes will be united without overlap in BCs or bodies.\n");
3951
    }   
3952
    else if(strcmp(argv[arg],"-nonames") == 0) {
3953
      eg->usenames = FALSE;
3954
      printf("Names will be omitted even if they would exist\n");
3955
    }   
3956
    else if(strcmp(argv[arg],"-multidim") == 0) {
3957
      eg->multidim = TRUE;
3958
      printf("Lower dimensional entities may be bulk too!\n");
3959
    }   
3960
    else if(strcmp(argv[arg],"-removelowdim") == 0) {
3961
      eg->removelowdim = TRUE;
3962
      printf("Lower dimensional boundaries will be removed\n");
3963
    }   
3964
    else if(strcmp(argv[arg],"-removeintbcs") == 0) {
3965
      eg->removeintbcs = TRUE;
3966
      printf("Lower dimensional boundaries will be removed\n");
3967
    }   
3968
    else if(strcmp(argv[arg],"-removeunused") == 0) {
3969
      eg->removeunused = TRUE;
3970
      printf("Nodes that do not appear in any element will be removed\n");
3971
    }   
3972
    else if(strcmp(argv[arg],"-autoclean") == 0) {
3973
      eg->removelowdim = TRUE;
3974
      eg->bulkorder = TRUE;
3975
      eg->boundorder = TRUE;
3976
      eg->removeunused = TRUE;
3977
      printf("Lower dimensional boundaries will be removed\n");
3978
      printf("Materials and boundaries will be renumbered\n");
3979
      printf("Nodes that do not appear in any element will be removed\n");
3980
    }   
3981
    else if(strcmp(argv[arg],"-polar") == 0) {
3982
      eg->polar = TRUE;
3983
      printf("Making transformation to polar coordinates.\n");
3984
      if(arg+1 >= argc) {
3985
	printf("The preferred radius is required as a parameter\n");
3986
	eg->polarradius = 1.0;
3987
      }
3988
      else {
3989
	eg->polarradius = atoi(argv[arg+1]);
3990
      }
3991
    }
3992
    else if(strcmp(argv[arg],"-cylinder") == 0) {
3993
      eg->cylinder = TRUE;
3994
      printf("Making transformation from cylindrical to cartesian coordinates.\n");
3995
    }
3996
    else if(strcmp(argv[arg],"-reduce") == 0) {
3997
      if(arg+2 >= argc) {
3998
	printf("Give two material for the interval.\n");
3999
 	return(12);
4000
      }
4001
      else {
4002
	eg->reduce = TRUE;      
4003
	eg->reducemat1 = atoi(argv[arg+1]);
4004
	eg->reducemat2 = atoi(argv[arg+2]);
4005
      }
4006
    }
4007
    else if(strcmp(argv[arg],"-increase") == 0) {
4008
      eg->increase = TRUE;
4009
    }
4010
    else if(strcmp(argv[arg],"-bulkorder") == 0) {
4011
      eg->bulkorder = TRUE;
4012
    }
4013
    else if(strcmp(argv[arg],"-boundorder") == 0) {
4014
      eg->boundorder = TRUE;
4015
    }
4016
    else if(strcmp(argv[arg],"-partition") == 0  ||
4017
       strcmp(argv[arg],"-partcell") == 0  || 
4018
       strcmp(argv[arg],"-partcyl") == 0 ) {
4019
      if(arg+dim >= argc) {
4020
	printf("The number of partitions in %d dims is required as parameters.\n",dim);
4021
	return(13);
4022
      }
4023
      else {
4024
	eg->partitions = 1;
4025
	eg->partdim[0] = atoi(argv[arg+1]);
4026
	eg->partdim[1] = atoi(argv[arg+2]);
4027
	if(dim == 3) eg->partdim[2] = atoi(argv[arg+3]);
4028
	eg->partitions = 1;
4029
	for(i=0;i<3;i++) {
4030
	  if(eg->partdim[i] == 0) eg->partdim[i] = 1;
4031
	  eg->partitions *= eg->partdim[i];
4032
	}
4033
	eg->partopt = -1;
4034
	if( strcmp(argv[arg],"-partition") == 0 ) {
4035
	  if(arg+4 < argc) 
4036
	    if(argv[arg+4][0] != '-') eg->partopt = atoi(argv[arg+4]);
4037
	}
4038
	else if( strcmp( argv[arg],"-partcell") == 0 )  {
4039
	  eg->partopt = 2;
4040
	} else if( strcmp( argv[arg],"-partcyl") == 0 ) {
4041
	  eg->partopt = 3;
4042
	}
4043
	  
4044
	printf("The mesh will be partitioned geometrically to %d partitions.\n",
4045
	       eg->partitions);
4046
      }
4047
    }
4048
    else if(strcmp(argv[arg],"-partorder") == 0) {
4049
      if(arg+dim >= argc) {
4050
	printf("Give %d parameters for the order vector.\n",dim);
4051
 	return(14);
4052
      }
4053
      else {
4054
	eg->partorder = 1;
4055
	eg->partcorder[0] = atof(argv[arg+1]);
4056
	eg->partcorder[1] = atof(argv[arg+2]);
4057
	if(dim==3) eg->partcorder[2] = atof(argv[arg+3]);
4058
      }
4059
    }
4060
    else if(strcmp(argv[arg],"-partoptim") == 0) {
4061
      eg->partoptim = TRUE;
4062
      printf("Aggressive optimization will be applied to node sharing.\n");
4063
    }
4064
    else if(strcmp(argv[arg],"-partnobcoptim") == 0) {
4065
      eg->partbcoptim = FALSE;
4066
      printf("Aggressive optimization will not be applied to parent element sharing.\n");
4067
    }
4068
    else if(strcmp(argv[arg],"-partbw") == 0) {
4069
      eg->partbw = TRUE;
4070
      printf("Bandwidth will be optimized for partitions.\n");
4071
    }
4072
    else if(strcmp(argv[arg],"-parthypre") == 0) {
4073
      eg->parthypre = TRUE;
4074
      printf("Numbering of partitions will be made continuous.\n");
4075
    }
4076
    else if(strcmp(argv[arg],"-partdual") == 0) {
4077
      eg->partdual = TRUE;
4078
      printf("Using dual (elemental) graph in partitioning.\n");
4079
    }
4080
    else if(strcmp(argv[arg],"-metis") == 0 ||
4081
       strcmp(argv[arg],"-metisrec") == 0 ||
4082
       strcmp(argv[arg],"-metiskway") == 0 ) {
4083
#if USE_METIS
4084
      if(arg+1 >= argc) {
4085
	printf("The number of partitions is required as a parameter\n");
4086
	return(15);
4087
      }
4088
      else {
4089
	eg->metis = atoi(argv[arg+1]);
4090
	printf("The mesh will be partitioned with Metis to %d partitions.\n",eg->metis);
4091
	eg->partopt = 0;
4092
	if(strcmp(argv[arg],"-metisrec") == 0)
4093
	  eg->partopt = 2;
4094
	else if(strcmp(argv[arg],"-metiskway") == 0 )
4095
	  eg->partopt = 3;
4096
	else if(arg+2 < argc) 
4097
	  if(argv[arg+2][0] != '-') eg->partopt = atoi(argv[arg+2]);
4098
      }    
4099
#else
4100
      printf("This version of ElmerGrid was compiled without Metis library!\n");
4101
#endif     
4102
    }    
4103
    else if(strcmp(argv[arg],"-metisseed") == 0 ) {
4104
      if(arg+1 >= argc) {
4105
	printf("The random number seed is required as parameter for -metisseed!\n");
4106
	return(15);
4107
      }
4108
      else {
4109
	eg->metis_seed = atoi(argv[arg+1]);
4110
	printf("Seed for Metis partitioning routines: %d\n",eg->metis_seed);
4111
      }
4112
    }
4113
    else if(strcmp(argv[arg],"-metisncuts") == 0 ) {
4114
      if(arg+1 >= argc) {
4115
	printf("The number of parameters is required as parameter for -metisncuts!\n");
4116
	return(15);
4117
      }
4118
      else {
4119
	eg->metis_ncuts = atoi(argv[arg+1]);
4120
	printf("Number of competing partitions to generate : %d\n",eg->metis_ncuts);
4121
      }
4122
    }   
4123
    else if(strcmp(argv[arg],"-partjoin") == 0) {
4124
      if(arg+1 >= argc) {
4125
	printf("The number of partitions is required as a parameter!\n");
4126
	return(15);
4127
      }
4128
      else {
4129
	eg->partjoin = atoi(argv[arg+1]);
4130
	printf("The results will joined using %d partitions.\n",eg->partjoin);
4131
      }
4132
    }
4133
    else if(strcmp(argv[arg],"-partconnect") == 0 || strcmp(argv[arg],"-partzbc") == 0 ) {
4134
      if(arg+1 >= argc) {
4135
	printf("The number of 1D partitions is required as a parameter!\n");
4136
	return(15);
4137
      }
4138
      else {
4139
	eg->partbcz = atoi(argv[arg+1]);
4140
	printf("The connected BCs will be partitioned to %d partitions in Z.\n",eg->partbcz);
4141
      }
4142
    }
4143
    else if(strcmp(argv[arg],"-partrbc") == 0 ) {
4144
      if(arg+1 >= argc) {
4145
	printf("The number of 1D partitions is required as a parameter!\n");
4146
	return(15);
4147
      }
4148
      else {
4149
	eg->partbcr = atoi(argv[arg+1]);
4150
	printf("The connected BCs will be partitioned to %d partitions in R.\n",eg->partbcr);
4151
      }
4152
    }
4153
    else if(strcmp(argv[arg],"-partlayers") == 0) {
4154
      if(arg+1 >= argc) {
4155
	printf("The number of layers to be extended is required as a parameter\n");
4156
	return(15);
4157
      }
4158
      else {
4159
	eg->partbclayers = atoi(argv[arg+1]);
4160
	printf("The boundary partitioning will be extended by %d layers.\n",eg->partbclayers);
4161
      }
4162
    }
4163
    else if(strcmp(argv[arg],"-metiscontig") == 0 ) {
4164
      eg->metis_contig = TRUE;
4165
    }
4166
    else if(strcmp(argv[arg],"-metisvol") == 0 ) {
4167
      eg->metis_volcut = TRUE;
4168
    }
4169
    else if(strcmp(argv[arg],"-metisminconn") == 0 ) {
4170
      eg->metis_minconn = TRUE;
4171
    }    
4172
    else if(strcmp(argv[arg],"-metisconnect") == 0 || strcmp(argv[arg],"-metisbc") == 0 ) {
4173
      if(arg+1 >= argc) {
4174
	printf("The number of Metis partitions is required as a parameter\n");
4175
	return(15);
4176
      }
4177
      else {
4178
	eg->partbcmetis = atoi(argv[arg+1]);
4179
	printf("The connected BCs will be partitioned to %d partitions by Metis.\n",eg->partbcmetis);
4180
      }
4181
    }
4182
    else if(strcmp(argv[arg],"-periodic") == 0) {
4183
      if(arg+dim >= argc) {
4184
	printf("Give the periodic coordinate directions (e.g. 1 1 0)\n");
4185
 	return(16);
4186
      }
4187
      else {
4188
	eg->periodicdim[0] = atoi(argv[arg+1]);
4189
	eg->periodicdim[1] = atoi(argv[arg+2]);
4190
	if(dim == 3) eg->periodicdim[2] = atoi(argv[arg+3]);
4191
      }
4192
    }
4193
    else if(strcmp(argv[arg],"-discont") == 0) {
4194
      if(arg+1 >= argc) {
4195
	printf("Give the discontinuous boundary conditions.\n");
4196
 	return(17);
4197
      }
4198
      else {
4199
	eg->discontbounds[eg->discont] = atoi(argv[arg+1]);
4200
	eg->discont++;
4201
      }
4202
    }
4203
    else if(strcmp(argv[arg],"-connect") == 0) {
4204
      if(arg+1 >= argc) {
4205
	printf("Give the connected boundary conditions.\n");
4206
 	return(10);
4207
      }
4208
      else {
4209
	eg->connectboundsnosets += 1;
4210
	for(i=arg+1;i<argc && strncmp(argv[i],"-",1); i++) { 
4211
	  eg->connectbounds[eg->connect] = atoi(argv[i]);
4212
	  eg->connectboundsset[eg->connect] = eg->connectboundsnosets;
4213
	  eg->connect++;
4214
	}
4215
      }
4216
    } 
4217
    else if(strcmp(argv[arg],"-connectall") == 0) {
4218
      eg->connectboundsnosets += 1;
4219
      eg->connectbounds[eg->connect] = -1;
4220
      eg->connectboundsset[eg->connect] = eg->connectboundsnosets;
4221
      eg->connect++;
4222
    }
4223
    else if(strcmp(argv[arg],"-connectint") == 0) {
4224
      eg->connectboundsnosets += 1;
4225
      eg->connectbounds[eg->connect] = -2;
4226
      eg->connectboundsset[eg->connect] = eg->connectboundsnosets;
4227
      eg->connect++;
4228
    }
4229
    else if(strcmp(argv[arg],"-connectfree") == 0) {
4230
      eg->connectboundsnosets += 1;
4231
      eg->connectbounds[eg->connect] = -3;
4232
      eg->connectboundsset[eg->connect] = eg->connectboundsnosets;
4233
      eg->connect++;
4234
    }
4235
    else if(strcmp(argv[arg],"-boundbound") == 0) {
4236
      for(i=arg+1;i<=arg+3 && i<argc; i++) {
4237
	eg->boundbound[3*eg->boundbounds+i-(1+arg)] = atoi(argv[i]);
4238
	if((i-arg)%3 == 0) eg->boundbounds++;
4239
      }
4240
    } 
4241
    else if(strcmp(argv[arg],"-bulkbound") == 0) {
4242
      for(i=arg+1;i<=arg+3 && i<argc; i++) {
4243
	eg->bulkbound[3*eg->bulkbounds+i-(1+arg)] = atoi(argv[i]);
4244
	if((i-arg)%3 == 0) eg->bulkbounds++;
4245
      }
4246
    } 
4247
    else if(strcmp(argv[arg],"-boundtype") == 0) {
4248
      for(i=arg+1;i<argc && strncmp(argv[i],"-",1); i++) 
4249
	eg->sidemap[3*eg->sidemappings+i-1-arg] = atoi(argv[i]);
4250
      eg->sidemappings++;
4251
    } 
4252
    else if(strcmp(argv[arg],"-bulktype") == 0) {
4253
      for(i=arg+1;i<argc && strncmp(argv[i],"-",1); i++) 
4254
	eg->bulkmap[3*eg->bulkmappings+i-1-arg] = atoi(argv[i]);
4255
      eg->bulkmappings++;
4256
    } 
4257
    else if(strcmp(argv[arg],"-coordinatemap") == 0) {
4258
      if( arg+3 >= argc ) {
4259
	printf("Give three parameters for the index permutation\n");
4260
	return(18);
4261
      }
4262
      else {
4263
	for(i=0;i<3;i++) 
4264
	  eg->coordinatemap[i] = atoi(argv[arg+1+i]);
4265
      }
4266
    } 
4267
    else if(strcmp(argv[arg],"-layer") == 0) {
4268
      if(arg+4 >= argc) {
4269
	printf("Give four parameters for the layer: boundary, elements, thickness, ratio.\n");
4270
	return(18);
4271
      }
4272
      else if(eg->layers == MAXBOUNDARIES) {
4273
	printf("There can only be %d layers, sorry.\n",MAXBOUNDARIES);
4274
	return(19);
4275
      }
4276
      else {
4277
	eg->layerbounds[eg->layers] = atoi(argv[arg+1]);
4278
	eg->layernumber[eg->layers] = atoi(argv[arg+2]);
4279
	eg->layerthickness[eg->layers] = atof(argv[arg+3]);
4280
	eg->layerratios[eg->layers] = atof(argv[arg+4]);
4281
	eg->layerparents[eg->layers] = 0;
4282
	eg->layers++;
4283
      }
4284
    }   
4285
    else if(strcmp(argv[arg],"-layermove") == 0) {
4286
      if(arg+1 >= argc) {
4287
	printf("Give maximum number of Jacobi filters.\n");
4288
 	return(20);
4289
      }
4290
      else {
4291
	eg->layermove = atoi(argv[arg+1]);
4292
      }
4293
    }
4294
    /* This uses a very dirty trick where the variables related to argument -layer are used 
4295
       with a negative indexing */ 
4296
    else if(strcmp(argv[arg],"-divlayer") == 0) {
4297
      if(arg+4 >= argc) {
4298
	printf("Give four parameters for the layer: boundary, elements, relative thickness, ratio.\n");
4299
	return(21);
4300
      }
4301
      else if(abs(eg->layers) == MAXBOUNDARIES) {
4302
	printf("There can only be %d layers, sorry.\n",MAXBOUNDARIES);
4303
	return(22);
4304
      }
4305
      else {
4306
	eg->layerbounds[abs(eg->layers)] = atoi(argv[arg+1]);
4307
	eg->layernumber[abs(eg->layers)] = atoi(argv[arg+2]);
4308
	eg->layerthickness[abs(eg->layers)] = atof(argv[arg+3]);
4309
	eg->layerratios[abs(eg->layers)] = atof(argv[arg+4]);
4310
	eg->layerparents[abs(eg->layers)] = 0;
4311
	eg->layers--;
4312
      }
4313
    }
4314
    else if(strcmp(argv[arg],"-3d") == 0) {
4315
      eg->dim = dim = 3;
4316
    }
4317
    else if(strcmp(argv[arg],"-2d") == 0) {
4318
      eg->dim = dim = 2;
4319
    }
4320
    else if(strcmp(argv[arg],"-1d") == 0) {
4321
      eg->dim = dim = 1;
4322
    }
4323
    else if(strcmp(argv[arg],"-isoparam") == 0) {
4324
      eg->isoparam = TRUE;
4325
    }
4326
    else if(strcmp(argv[arg],"-nobound") == 0) {
4327
      eg->saveboundaries = FALSE;
4328
    }
4329
    else if(strcmp(argv[arg],"-vtuone") == 0) {
4330
      eg->vtuone = TRUE;
4331
    }
4332
    else if(strcmp(argv[arg],"-nosave") == 0) {
4333
      eg->nosave = TRUE;
4334
    }
4335
    else if(strcmp(argv[arg],"-nooverwrite") == 0) {
4336
      eg->nooverwrite = TRUE;
4337
    }
4338
    else if(strcmp(argv[arg],"-timer") == 0) {
4339
      eg->timeron = TRUE;
4340
    }
4341
    else if(strcmp(argv[arg],"-infofile") == 0) {
4342
      eg->timeron = TRUE;
4343
      if(arg+1 >= argc) {
4344
	printf("The output name is required as a parameter\n");
4345
	return(2);
4346
      }
4347
      else {
4348
	strcpy(eg->infofile,argv[arg+1]);
4349
      }
4350
    }    
4351
    /* The following keywords are not actively used */
4352
    else if(strcmp(argv[arg],"-bcoffset") == 0) {
4353
      eg->bcoffset = atoi(argv[arg+1]);
4354
    }
4355
    else if(strcmp(argv[arg],"-noelements") == 0) {
4356
      eg->elements3d = atoi(argv[arg+1]);
4357
    }
4358
    else if(strcmp(argv[arg],"-nonodes") == 0) {
4359
      eg->nodes3d = atoi(argv[arg+1]);
4360
    }
4361
    else if(strcmp(argv[arg],"-sidefind") == 0) {
4362
      eg->findsides = 0;
4363
      for(i=arg+1;i<argc && strncmp(argv[i],"-",1); i++) {
4364
	eg->sidebulk[i-1-arg] = atoi(argv[i]);
4365
	eg->findsides++;
4366
      }
4367
    } 
4368
    else if(strcmp(argv[arg],"-findbound") == 0) {
4369
      eg->findsides = 0;
4370
      for(i=arg+1;i+1<argc && strncmp(argv[i],"-",1); i += 2) {
4371
	eg->sidebulk[i-1-arg] = atoi(argv[i]);
4372
	eg->sidebulk[i-arg] = atoi(argv[i+1]);
4373
	eg->findsides++;
4374
      }
4375
    }
4376
    else if(strcmp(argv[arg],"-") == 0 ) {
4377
      printf("Unknown in-line argument: %s\n",argv[arg]);
4378
      bigerror("Cannot deal with argument, aborting!");
4379
    }
4380
    
4381
  }
4382

4383
  {
4384
    int badpoint;
4385
    char *ptr1,*ptr2;
4386
    ptr1 = strrchr(eg->filesout[0], '.');
4387
    if(ptr1) {
4388
      badpoint=FALSE;
4389
      ptr2 = strrchr(eg->filesout[0], '/');
4390
      if(ptr2 && ptr2 > ptr1) badpoint = TRUE;
4391
      if(!badpoint) *ptr1 = '\0';
4392
    }
4393
  }
4394

4395
  printf("Output will be saved to file %s.\n",eg->filesout[0]);
4396

4397
  return(0);
4398
}
4399

4400

4401

4402

4403
int LoadCommands(char *prefix,struct ElmergridType *eg,
4404
		 struct GridType *grid, int mode,int info) 
4405
{
4406
  char filename[MAXFILESIZE];
4407
  char command[MAXLINESIZE],params[MAXLINESIZE],*cp;
4408

4409
  FILE *in=NULL;
4410
  int i,j,iostat;
4411

4412
  iodebug = FALSE;
4413

4414
  if( mode == 0) {  
4415
    if ((in = fopen("ELMERGRID_STARTINFO","r"))) {
4416
      iostat = fscanf(in,"%s",filename);
4417
      fclose(in);
4418
      printf("Using the file %s defined in ELMERGRID_STARTINFO\n",filename);
4419
      if ((in = fopen(filename,"r")) == NULL) {
4420
	printf("LoadCommands: opening of the file '%s' wasn't successful !\n",filename);
4421
	return(1);
4422
      }    
4423
      else printf("Loading ElmerGrid commands from file '%s'.\n",filename);    
4424
    }    
4425
    else 
4426
      return(2);
4427
  }
4428
  else if(mode == 1) { 
4429
    AddExtension(prefix,filename,"eg");
4430
    if ((in = fopen(filename,"r")) == NULL) {
4431
      printf("LoadCommands: opening of the file '%s' wasn't successful !\n",filename);
4432
      return(3);
4433
    }    
4434
    if(info) printf("Loading ElmerGrid commands from file '%s'.\n",filename);    
4435
  }
4436
  else if(mode == 2) {
4437
    AddExtension(prefix,filename,"grd");
4438
    if ((in = fopen(filename,"r")) == NULL) {
4439
      printf("LoadCommands: opening of the file '%s' wasn't successful !\n",filename);
4440
      return(4);
4441
    }    
4442
    if(info) printf("\nLoading ElmerGrid commands from file '%s'.\n",filename);
4443
  }
4444

4445

4446

4447
  for(;;) {
4448

4449
    if(GetCommand(command,params,in)) {
4450
      printf("Reached the end of command file\n");
4451
      goto end;
4452
    }    
4453

4454
    /* If the mode is the command file mode read also the file information from the command file. */
4455

4456
    if(mode <= 1) {
4457
      if(strstr(command,"INPUT FILE")) {
4458
	sscanf(params,"%s",eg->filesin[0]);
4459
      }
4460

4461
      else if(strstr(command,"OUTPUT FILE")) {
4462
	sscanf(params,"%s",eg->filesout[0]);
4463
	eg->filerenamed = TRUE;
4464
      }
4465

4466
      else if(strstr(command,"INPUT MODE")) {
4467
	for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4468
	
4469
	for(i=0;i<=MAXINMETHODS;i++) {
4470
	  if(strstr(params,InMethods[i])) {
4471
	    eg->inmethod = i;
4472
	    break;
4473
	  }
4474
	}
4475
	if(i>MAXINMETHODS) sscanf(params,"%d",&eg->inmethod);
4476
      }
4477

4478
      else if(strstr(command,"OUTPUT MODE")) {
4479
	for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4480
	
4481
	/* Type of output file (fewer options) */
4482
	for(i=1;i<=MAXOUTMETHODS;i++) {
4483
	  if(strstr(params,OutMethods[i])) {
4484
	    eg->outmethod = i;
4485
	    break;
4486
	  }
4487
	}
4488
	if(i>MAXOUTMETHODS) sscanf(params,"%d",&eg->outmethod);	
4489
      }
4490
    }    
4491
    /* End of command file specific part */
4492

4493

4494
    if(strstr(command,"DECIMALS")) {
4495
      sscanf(params,"%d",&eg->decimals);
4496
    }
4497
    else if(strstr(command,"BOUNDARY OFFSET")) {
4498
      sscanf(params,"%d",&eg->bcoffset);
4499
    }
4500
    else if(strstr(command,"TRIANGLES CRITICAL ANGLE")) {
4501
      sscanf(params,"%le",&eg->triangleangle);
4502
    }
4503
    else if(strstr(command,"TRIANGLES")) {
4504
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4505
      if(strstr(params,"TRUE")) eg->triangles = TRUE;      
4506
    }
4507
    else if(strstr(command,"MERGE NODES")) {
4508
      eg->merge = TRUE;
4509
      sscanf(params,"%le",&eg->cmerge);
4510
    }
4511
    else if(strstr(command,"UNITE")) {
4512
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4513
      if(strstr(params,"TRUE")) eg->unitemeshes = TRUE;      
4514
    }
4515
    else if(strstr(command,"UNITENOOVERLAP")) {
4516
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4517
      if(strstr(params,"TRUE")) {
4518
	eg->unitemeshes = TRUE;      
4519
	eg->unitenooverlap = TRUE;      
4520
      }
4521
    }
4522
    
4523
    else if(strstr(command,"ORDER NODES")) {
4524
      eg->order = TRUE;
4525
      if(eg->dim == 1) 
4526
	sscanf(params,"%le",&eg->corder[0]);
4527
      else if(eg->dim == 2) 
4528
	sscanf(params,"%le%le",&eg->corder[0],&eg->corder[1]);
4529
      else if(eg->dim == 3) 
4530
	sscanf(params,"%le%le%le",&eg->corder[0],&eg->corder[1],&eg->corder[2]);
4531
    }
4532
    else if(strstr(command,"SCALE")) {
4533
      eg->scale = TRUE;
4534
      if(eg->dim == 1) 
4535
	sscanf(params,"%le",&eg->cscale[0]);
4536
      else if(eg->dim == 2) 
4537
	sscanf(params,"%le%le",&eg->cscale[0],&eg->cscale[1]);
4538
      else if(eg->dim == 3) 
4539
	sscanf(params,"%le%le%le",&eg->cscale[0],&eg->cscale[1],&eg->cscale[2]);
4540
    }
4541
    else if(strstr(command,"CENTRALIZE")) {
4542
      eg->center = TRUE;
4543
    }
4544
    else if(strstr(command,"TRANSLATE")) {
4545
      eg->translate = TRUE;
4546
      if(eg->dim == 1) 
4547
	sscanf(params,"%le",&eg->ctranslate[0]);
4548
      else if(eg->dim == 2) 
4549
	sscanf(params,"%le%le",&eg->ctranslate[0],&eg->ctranslate[1]);
4550
      else if(eg->dim == 3) 
4551
	sscanf(params,"%le%le%le",&eg->ctranslate[0],&eg->ctranslate[1],&eg->ctranslate[2]);
4552
    }
4553
    else if(strstr(command,"ROTATE MESH")) {
4554
      eg->rotate = TRUE;
4555
      sscanf(params,"%le%le%le",&eg->crotate[0],&eg->crotate[1],&eg->crotate[2]);
4556
    }
4557
    else if(strstr(command,"CLONE")) {
4558
      if(strstr(command,"CLONE SIZE")) {
4559
	if(eg->dim == 1) 
4560
	  sscanf(params,"%le",&eg->clonesize[0]);
4561
	else if(eg->dim == 2) 
4562
	  sscanf(params,"%le%le",&eg->clonesize[0],&eg->clonesize[1]);
4563
	else if(eg->dim == 3) 
4564
	  sscanf(params,"%le%le%le",&eg->clonesize[0],&eg->clonesize[1],&eg->clonesize[2]);	
4565
      }
4566
      else {
4567
	if(eg->dim == 1) 
4568
	  sscanf(params,"%d",&eg->clone[0]);
4569
	else if(eg->dim == 2) 
4570
	  sscanf(params,"%d%d",&eg->clone[0],&eg->clone[1]);
4571
	else if(eg->dim == 3) 
4572
	  sscanf(params,"%d%d%d",&eg->clone[0],&eg->clone[1],&eg->clone[2]);
4573
      }
4574
    }
4575
    else if(strstr(command,"MERGE")) {
4576
      eg->merge = TRUE;
4577
      sscanf(params,"%le",&eg->cmerge);
4578
    }
4579
    else if(strstr(command,"MIRROR")) {
4580
      if(eg->dim == 1) 
4581
	sscanf(params,"%d",&eg->mirror[0]);
4582
      else if(eg->dim == 2) 
4583
	sscanf(params,"%d%d",&eg->mirror[0],&eg->mirror[1]);
4584
      else if(eg->dim == 3) 
4585
	sscanf(params,"%d%d%d",&eg->mirror[0],&eg->mirror[1],&eg->mirror[2]);
4586
    }
4587
    else if(strstr(command,"POLAR RADIUS")) {
4588
      eg->polar = TRUE;
4589
      sscanf(params,"%le",&eg->polarradius);
4590
    }
4591
    else if(strstr(command,"CYLINDER")) {
4592
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4593
      if(strstr(params,"TRUE")) eg->cylinder = TRUE;      
4594
    }
4595
    else if(strstr(command,"REDUCE DEGREE")) {
4596
      eg->reduce = TRUE;
4597
      sscanf(params,"%d%d",&eg->reducemat1,&eg->reducemat2);
4598
    }
4599
    else if(strstr(command,"INCREASE DEGREE")) {
4600
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4601
      if(strstr(params,"TRUE")) eg->increase = TRUE;      
4602
    }
4603
    else if(strstr(command,"METIS OPTION")) {
4604
#if USE_METIS
4605
      sscanf(params,"%d",&eg->partopt);
4606
#else
4607
      printf("This version of ElmerGrid was compiled without Metis library!\n");
4608
#endif
4609
    }
4610
    else if(strstr(command,"METIS")) {
4611
#if USE_METIS
4612
      sscanf(params,"%d",&eg->metis);
4613
#else
4614
      printf("This version of ElmerGrid was compiled without Metis library!\n");
4615
#endif
4616
    }
4617
    else if(strstr(command,"METISKWAY")) {
4618
#if USE_METIS
4619
      sscanf(params,"%d",&eg->metis);
4620
      eg->partopt = 3;
4621
#else
4622
      printf("This version of ElmerGrid was compiled without Metis library!\n");
4623
#endif
4624
    }
4625
    else if(strstr(command,"METISREC")) {
4626
#if USE_METIS
4627
      sscanf(params,"%d",&eg->metis);
4628
      eg->partopt = 2;
4629
#else
4630
      printf("This version of ElmerGrid was compiled without Metis library!\n");
4631
#endif
4632
    }
4633
    else if(strstr(command,"PARTITION DUAL")) {
4634
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4635
      if(strstr(params,"TRUE")) eg->partdual = TRUE;      
4636
    }
4637
    else if(strstr(command,"PARTJOIN")) {
4638
      sscanf(params,"%d",&eg->partjoin);
4639
    }
4640
    else if(strstr(command,"PARTITION ORDER")) {
4641
      eg->partorder = 1;
4642
      if(eg->dim == 2) sscanf(params,"%le%le",&eg->partcorder[0],&eg->partcorder[1]);
4643
      if(eg->dim == 3) sscanf(params,"%le%le%le",&eg->partcorder[0],
4644
			      &eg->partcorder[1],&eg->partcorder[2]);      
4645
    }
4646
    else if(strstr(command,"PARTITION") || strstr(command,"PARTCYL") || strstr(command,"PARTCELL")) {
4647
      if(eg->dim == 2) sscanf(params,"%d%d",&eg->partdim[0],&eg->partdim[1]);
4648
      if(eg->dim == 3) sscanf(params,"%d%d%d",&eg->partdim[0],&eg->partdim[1],&eg->partdim[2]);
4649
      eg->partitions = 1;
4650
      for(i=0;i<eg->dim;i++) {
4651
	if(eg->partdim[i] < 1) eg->partdim[i] = 1;
4652
	eg->partitions *= eg->partdim[i];
4653
      }
4654
      if( strstr(command,"PARTCYL") ) eg->partopt = 3;
4655
      if( strstr(command,"PARTCCELL") ) eg->partopt = 2;
4656
    }
4657
    else if(strstr(command,"PERIODIC")) {
4658
      if(eg->dim == 2) sscanf(params,"%d%d",&eg->periodicdim[0],&eg->periodicdim[1]);
4659
      if(eg->dim == 3) sscanf(params,"%d%d%d",&eg->periodicdim[0],
4660
			      &eg->periodicdim[1],&eg->periodicdim[2]);
4661
    }
4662
    else if(strstr(command,"HALO")) {
4663
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4664
      if(strstr(params,"TRUE")) eg->parthalo[1] = TRUE;      
4665
    }
4666
    else if(strstr(command,"BOUNDARY HALO")) {
4667
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4668
      if(strstr(params,"TRUE")) eg->parthalo[2] = TRUE;
4669
    }
4670
    else if(strstr(command,"EXTRUDED HALO")) {
4671
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4672
      if(strstr(params,"TRUE")) eg->parthalo[3] = TRUE;
4673
    }
4674
    else if(strstr(command,"GREEDY HALO")) {
4675
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4676
      if(strstr(params,"TRUE")) eg->parthalo[4] = TRUE;
4677
    }    
4678
    else if(strstr(command,"PARTBW")) {
4679
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4680
      if(strstr(params,"TRUE")) eg->partbw = TRUE;      
4681
    }
4682
    else if(strstr(command,"PARTHYPRE")) {
4683
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4684
      if(strstr(params,"TRUE")) eg->parthypre = TRUE;      
4685
    }
4686
    else if(strstr(command,"BOUNDARY TYPE MAPPINGS")) {
4687
      for(i=0;i<MAXMAPPINGS;i++) {
4688
	if(i>0) Getline(params,in);
4689
	for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4690
	if(strstr(params,"END")) break;
4691
	cp = params;      
4692
	sscanf(params,"%d%d%d",&eg->sidemap[3*i],&eg->sidemap[3*i+1],&eg->sidemap[3*i+2]);
4693
      }
4694
      printf("Found %d boundary type mappings\n",i);
4695
      eg->sidemappings = i;
4696
    }
4697
    else if(strstr(command,"BULK TYPE MAPPINGS")) {
4698
      for(i=0;i<MAXMAPPINGS;i++) {
4699
	if(i>0) Getline(params,in);
4700
	for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4701
	if(strstr(params,"END")) break;
4702
	cp = params;      
4703
	sscanf(params,"%d%d%d",&eg->bulkmap[3*i],&eg->bulkmap[3*i+1],&eg->bulkmap[3*i+2]);
4704
      }
4705
      printf("Found %d bulk type mappings\n",i);
4706
      eg->bulkmappings = i;
4707
    }
4708
    else if(strstr(command,"COORDINATE MAPPING")) {
4709
      sscanf(params,"%d%d%d",&eg->coordinatemap[0],&eg->coordinatemap[1],&eg->coordinatemap[2]);
4710
    }	
4711
    else if(strstr(command,"BOUNDARY BOUNDARY")) {
4712
      for(i=0;i<MAXBOUNDARIES;i++) {
4713
	if(i>0) Getline(params,in);
4714
	for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4715
	if(strstr(params,"END")) break;
4716
	cp = params;      
4717
	sscanf(params,"%d%d%d",&eg->boundbound[3*i+2],&eg->boundbound[3*i],&eg->boundbound[3*i+1]);
4718
      }
4719
      printf("Found %d boundary boundary definitions\n",i);
4720
      eg->boundbounds = i;
4721
    }
4722
    else if(strstr(command,"MATERIAL BOUNDARY")) {
4723
      for(i=0;i<MAXBOUNDARIES;i++) {
4724
	if(i>0) Getline(params,in);
4725
	for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4726
	if(strstr(params,"END")) break;
4727
	cp = params;      
4728
	sscanf(params,"%d%d%d",&eg->bulkbound[3*i+2],&eg->bulkbound[3*i],&eg->bulkbound[3*i+1]);
4729
      }
4730
      printf("Found %d material boundary definitions\n",i);
4731
      eg->bulkbounds = i;
4732
    }
4733

4734
    else if(strstr(command,"RENUMBER BOUNDARY")) {
4735
      for(i=0;i<MAXBOUNDARIES;i++) {
4736
	for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4737
	if(strstr(params,"END")) break;
4738
	cp = params;      
4739
	sscanf(params,"%d%d%d",&eg->sidemap[3*i],&eg->sidemap[3*i+1],&eg->sidemap[3*i+2]);
4740
      }
4741
      printf("Found %d boundary mappings\n",i);
4742
      eg->sidemappings = i;
4743
    }
4744
    else if(strstr(command,"RENUMBER MATERIAL")) {
4745
      for(i=0;i<MAXBOUNDARIES;i++) {
4746
	for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4747
	if(strstr(params,"END")) break;
4748
	cp = params;      
4749
	sscanf(params,"%d%d%d",&eg->bulkmap[3*i],&eg->bulkmap[3*i+1],&eg->bulkmap[3*i+2]);
4750
      }
4751
      printf("Found %d material mappings\n",i);
4752
      eg->bulkmappings = i;
4753
    }
4754

4755
    else if(strstr(command,"BOUNDARY LAYER")) {
4756
      if(strstr(command,"BOUNDARY LAYER MOVE")) {
4757
	sscanf(params,"%d",&eg->layermove);
4758
      }
4759
      else if(strstr(command,"BOUNDARY LAYER EPSILON")) {
4760
	sscanf(params,"%le",&eg->layereps);
4761
      }
4762
      else {
4763
	for(i=0;i<MAXBOUNDARIES;i++) {
4764
	  if(i>0) Getline(params,in);
4765
	  for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4766
	  cp = params;      
4767

4768
	  if(strstr(params,"END") || strstr(params,"End") ) break;
4769
	  eg->layerbounds[i] = next_int(&cp);
4770
	  eg->layernumber[i] = next_int(&cp);
4771
	  eg->layerthickness[i] = next_real(&cp);
4772
	  eg->layerratios[i] = next_real(&cp);
4773
	  eg->layerparents[i] = next_int(&cp);	  
4774
	}
4775
	printf("Found %d boundary layers\n",i);
4776
	eg->layers = i;
4777
      }
4778
    }
4779
    else if(strstr(command,"REMOVE LOWER DIMENSIONAL BOUNDARIES")) {
4780
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4781
      if(strstr(params,"TRUE")) eg->removelowdim = TRUE; 
4782
    }
4783
    else if(strstr(command,"REMOVE INTERNAL BOUNDARIES")) {
4784
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4785
      if(strstr(params,"TRUE")) eg->removeintbcs = TRUE; 
4786
    }
4787
    else if(strstr(command,"REMOVE UNUSED NODES")) {
4788
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4789
      if(strstr(params,"TRUE")) eg->removeunused = TRUE; 
4790
    }
4791
    else if(strstr(command,"NO MESH NAMES")) {
4792
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4793
      if(strstr(params,"TRUE")) eg->usenames = FALSE; 
4794
    }
4795
    else if(strstr(command,"REORDER MATERIAL")) {
4796
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4797
      if(strstr(params,"TRUE")) eg->bulkorder = TRUE; 
4798
    }
4799
    else if(strstr(command,"REORDER BOUNDARY")) {
4800
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4801
      if(strstr(params,"TRUE")) eg->boundorder = TRUE; 
4802
    }
4803
    else if(strstr(command,"DIMENSION")) {
4804
      sscanf(params,"%d",&eg->dim);
4805
    }
4806
    else if(strstr(command,"ISOPARAMETRIC")) {
4807
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4808
      if(strstr(params,"TRUE")) eg->isoparam = TRUE;
4809
    }
4810
    else if(strstr(command,"NO BOUNDARY")) {
4811
      for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4812
      if(strstr(params,"TRUE")) eg->saveboundaries = FALSE;
4813
    }
4814
    else if(strstr(command,"LAYERED BOUNDARIES")) {
4815
      for(i=0;i<MAXLINESIZE;i++) params[i] = toupper(params[i]);
4816
      if(strstr(params,"TRUE")) grid->layeredbc = 1;
4817
      if(strstr(params,"FALSE")) grid->layeredbc = 0;
4818
    }
4819
    else if(strstr(command,"EXTRUDED")) {
4820
      grid->dimension = 3;
4821

4822
      if(strstr(command,"EXTRUDED DIVISIONS")) {
4823
	sscanf(params,"%d",&grid->zcells);		
4824
      }
4825
      if(strstr(command,"EXTRUDED BC OFFSET")) {
4826
	sscanf(params,"%d",&grid->layerbcoffset);		
4827
      }
4828
      else if(strstr(command,"EXTRUDED LIMITS")) {
4829
	cp = params;
4830
	for(i=0;i<=grid->zcells;i++ ) { 
4831
	  grid->z[i] = next_real(&cp);
4832
	  if(i > 0 && grid->z[i] < grid->z[i-1]) {
4833
	    printf("Extruded limits %d: %12.6le %12.6le\n",i,grid->z[i],grid->z[i-1]);
4834
	    bigerror("Values for limits should be a growing series, existing\n");
4835
	  }
4836
	}
4837
      }
4838
      else if(strstr(command,"EXTRUDED SIZES")) {
4839
	cp = params;
4840
	for(i=1;i<=grid->zcells;i++) grid->z[i] = next_real(&cp);
4841
	for(i=1;i<=grid->zcells;i++) grid->z[i] = grid->z[i-1] + grid->z[i];
4842
      }     
4843
      else if(strstr(command,"EXTRUDED ELEMENTS")) {
4844
	cp = params;
4845
	for(i=1;i<=grid->zcells;i++) grid->zelems[i] = next_int(&cp);
4846
	grid->autoratio = FALSE;    
4847
      } 
4848
      else if(strstr(command,"EXTRUDED RATIOS")) {
4849
	cp = params;
4850
	for(i=1;i<=grid->zcells;i++) grid->zexpand[i] = next_real(&cp);
4851
      }
4852
      else if(strstr(command,"EXTRUDED DENSITIES")) {
4853
	cp = params;
4854
	for(i=1;i<=grid->zcells;i++) grid->zdens[i] = next_real(&cp);
4855
      }
4856
      else if(strstr(command,"EXTRUDED STRUCTURE")) {
4857
	for(i=1;i<= grid->zcells;i++) {
4858
	  if(i>1) Getline(params,in);
4859
	  sscanf(params,"%d %d %d\n",
4860
		 &grid->zfirstmaterial[i],&grid->zlastmaterial[i],&grid->zmaterial[i]); 
4861
	}
4862
      }
4863
      else if(strstr(command,"EXTRUDED MAX MATERIAL")) {
4864
	sscanf(params,"%d",&grid->maxmaterial);		
4865
      }
4866
      else if(strstr(command,"EXTRUDED MATERIAL MAPPINGS")) {
4867
	grid->zmaterialmap = Imatrix(1,grid->zcells,1,grid->maxmaterial);
4868
	for(i=1;i<=grid->zcells;i++) {
4869
	  if(i>1) Getline(params,in); 
4870
	  cp = params;
4871
	  for(j=1;j<=grid->maxmaterial;j++)
4872
	    grid->zmaterialmap[i][j] = next_int(&cp);
4873
	}
4874
	grid->zmaterialmapexists = TRUE;
4875
      }
4876
      else if(strstr(command,"EXTRUDED HELICITY")) {
4877
	sscanf(params,"%le",&grid->zhelicity);			
4878
	grid->zhelicityexists = TRUE;
4879
      }
4880

4881
    }
4882
  }
4883

4884
end:
4885
  printf("Read commands from a file\n");
4886
  fclose(in);
4887
  return(0);
4888
}
4889

4890

4891

4892

4893

4894
static int FindParentSide(struct FemType *data,struct BoundaryType *bound,
4895
			  int sideelem,int sideelemtype,int *sideind)
4896
{
4897
  int i,j,sideelemtype2,elemind,parent,normal,elemtype;
4898
  int elemsides,side,sidenodes,nohits,hit,hit1,hit2;
4899
  int sideind2[MAXNODESD1];
4900
  int debug;
4901
  
4902
  hit = FALSE;
4903
  elemsides = 0;
4904
  elemtype = 0;
4905
  hit1 = FALSE;
4906
  hit2 = FALSE;
4907

4908
  debug = FALSE;
4909
  
4910
  for(parent=1;parent<=2;parent++) {
4911
    if(parent == 1) 
4912
      elemind = bound->parent[sideelem];
4913
    else
4914
      elemind = bound->parent2[sideelem];
4915
    
4916
    if(elemind > 0) {
4917
      elemtype = data->elementtypes[elemind];
4918
      elemsides = elemtype / 100;
4919

4920
      if(elemsides == 8) elemsides = 6;
4921
      else if(elemsides == 7) elemsides = 5;
4922
      else if(elemsides == 6) elemsides = 5;
4923
      else if(elemsides == 5) elemsides = 4;
4924
      
4925
      for(normal=1;normal >= -1;normal -= 2) {
4926

4927
	for(side=0;side<elemsides;side++) {
4928

4929
	  if(debug) printf("elem = %d %d %d %d\n",elemind,elemsides,normal,side);
4930

4931
	  GetElementSide(elemind,side,normal,data,&sideind2[0],&sideelemtype2);
4932
	  
4933
	  if(sideelemtype2 < 300 && sideelemtype > 300) break;	
4934
	  if(sideelemtype2 < 200 && sideelemtype > 200) break;		
4935
	  if(sideelemtype != sideelemtype2) continue;
4936

4937
	  sidenodes = sideelemtype / 100;
4938

4939
	  for(j=0;j<sidenodes;j++) {
4940
	    if(debug) printf("sidenode: %d %d %d\n",j,sideind[j],sideind2[j]);
4941

4942
	    hit = TRUE;
4943
	    for(i=0;i<sidenodes;i++) 
4944
	      if(sideind[(i+j)%sidenodes] != sideind2[i]) hit = FALSE;
4945
	    
4946
	    if(hit == TRUE) {
4947
	      if(parent == 1) {
4948
		hit1 = TRUE;
4949
		bound->side[sideelem] = side;
4950
		bound->normal[sideelem] = normal;
4951
	      }
4952
	      else {
4953
		hit2 = TRUE;
4954
		bound->side2[sideelem] = side;	      
4955
	      }
4956
	      goto skip;
4957
	    }
4958
	  }
4959
	}
4960
      }
4961

4962
      
4963
      /* this finding of sides does not guarantee that normals are oriented correctly */
4964
      normal = 1;
4965
      hit = FALSE;
4966
 
4967
      for(side=0;;side++) {
4968

4969
	if(0) printf("side = %d\n",side);
4970

4971
	GetElementSide(elemind,side,normal,data,&sideind2[0],&sideelemtype2);
4972

4973
	if(sideelemtype2 == 0 ) break;
4974
	if(sideelemtype2 < 300 && sideelemtype > 300) break;	
4975
	if(sideelemtype2 < 200 && sideelemtype > 200) break;		
4976
	if(sideelemtype != sideelemtype2) continue;
4977
		
4978
	sidenodes = sideelemtype % 100;
4979

4980
	nohits = 0;
4981
	for(j=0;j<sidenodes;j++) 
4982
	  for(i=0;i<sidenodes;i++) 
4983
	    if(sideind[i] == sideind2[j]) nohits++;
4984
	if(nohits == sidenodes) {
4985
	  hit = TRUE;
4986
	  if(parent == 1) {
4987
	    hit1 = TRUE;
4988
	    bound->side[sideelem] = side;
4989
	  }
4990
	  else {
4991
	    hit2 = TRUE;
4992
	    bound->side2[sideelem] = side;	      
4993
	  }
4994
	  goto skip;
4995
	}
4996
	
4997
      }
4998

4999
    skip:  
5000
      if(!hit) {
5001
	printf("FindParentSide: cannot locate BC element in parent %d: %d\n",parent,elemind);
5002
	printf("BC elem %d of type %d with corner indexes: ",sideelem,sideelemtype);
5003
	for(i=0;i<sideelemtype/100;i++)
5004
	  printf(" %d ",sideind[i]);
5005
	printf("\n");
5006

5007
	printf("Bulk elem %d of type %d with corner indexes: ",elemind,elemtype);
5008
	j = elemtype/100;
5009
	if( j >= 5 && j<=7 ) j = j-1;
5010
	for(i=0;i<j;i++)
5011
	  printf(" %d ",data->topology[elemind][i]);
5012
	printf("\n");             
5013
      }
5014
    }
5015
  }
5016

5017
  if(hit1 || hit2) 
5018
    return(0);
5019
  else
5020
    return(1);
5021
}
5022

5023

5024
static int Getnamerow(char *line,FILE *io,int upper) 
5025
{
5026
  int i,isend;
5027
  char *charend;
5028

5029

5030
  charend = fgets(line,MAXLINESIZE,io);
5031
  isend = (charend == NULL);
5032

5033
  if(isend) 
5034
    return(1);
5035
  else
5036
    return(0);
5037
}
5038

5039

5040

5041
int LoadElmerInput(struct FemType *data,struct BoundaryType *bound,
5042
		   char *prefix,int nonames, int info)
5043
/* This procedure reads the mesh assuming ElmerSolver format.
5044
   */
5045
{
5046
  int noknots,noelements,nosides,maxelemtype,maxnodes,nonodes;
5047
  int sideind[MAXNODESD1],tottypes,elementtype;
5048
  int i,j,k,l,dummyint,cdstat,fail;
5049
  int falseparents,noparents,bctopocreated;
5050
  int activeperm,activeelemperm,mini,maxi,minelem,maxelem,p1,p2;
5051
  int *nodeperm,*elemperm,*invperm,*invelemperm;
5052
  int iostat,noelements0;
5053
  FILE *in;
5054
  char line[MAXLINESIZE],line2[MAXLINESIZE],filename[MAXFILESIZE],directoryname[MAXFILESIZE];
5055
  char *ptr1,*ptr2;
5056

5057

5058
  sprintf(directoryname,"%s",prefix);
5059
  cdstat = chdir(directoryname);
5060

5061
  if(info) {
5062
    if(cdstat) 
5063
      printf("Loading mesh in ElmerSolver format from root directory.\n");
5064
    else
5065
      printf("Loading mesh in ElmerSolver format from directory %s.\n",directoryname);
5066
  }
5067

5068
  InitializeKnots(data);
5069

5070

5071
  sprintf(filename,"%s","mesh.header");
5072
  if ((in = fopen(filename,"r")) == NULL) {
5073
    printf("LoadElmerInput: The opening of the header-file %s failed!\n",
5074
	   filename);
5075
    return(1);
5076
  }
5077
  else 
5078
    printf("Loading header from %s\n",filename);
5079

5080
  GETLINE;
5081
  sscanf(line,"%d %d %d",&noknots,&noelements,&nosides);
5082
  GETLINE;
5083
  sscanf(line,"%d",&tottypes);
5084

5085
  maxelemtype = 0;
5086
  maxnodes = 0;
5087
  for(i=1;i<=tottypes;i++) {   
5088
    GETLINE;
5089
    sscanf(line,"%d",&dummyint);
5090
    maxelemtype = MAX( dummyint, maxelemtype );
5091
    j = maxelemtype % 100;
5092
    maxnodes = MAX( j, maxnodes );
5093
  }
5094
  printf("Maximum elementtype index is: %d\n",maxelemtype);
5095
  printf("Maximum number of nodes in element is: %d\n",maxnodes);
5096
  fclose(in);
5097

5098
  data->dim = GetElementDimension(maxelemtype);
5099

5100
  data->maxnodes = maxnodes;
5101
  data->noknots = noknots;
5102
  data->noelements = noelements0 = noelements;
5103

5104

5105
  if(info) printf("Allocating for %d knots and %d elements.\n",
5106
		  noknots,noelements);
5107
  AllocateKnots(data);
5108

5109

5110
  sprintf(filename,"%s","mesh.nodes");
5111
  if ((in = fopen(filename,"r")) == NULL) {
5112
    if(info) printf("LoadElmerInput: The opening of the nodes-file %s failed!\n",
5113
		    filename);
5114
    bigerror("Cannot continue without nodes file!\n");
5115
  }
5116
  else 
5117
    printf("Loading %d Elmer nodes from %s\n",noknots,filename);
5118

5119
  activeperm = FALSE;
5120
  for(i=1; i <= noknots; i++) {
5121
    GETLINE;
5122
    sscanf(line,"%d %d %le %le %le",
5123
	   &j, &dummyint, &(data->x[i]),&(data->y[i]),&(data->z[i]));
5124
    if(j != i && !activeperm) {
5125
      printf("LoadElmerInput: The node number (%d) at node %d is not compact, creating permutation\n",j,i);
5126
      activeperm = TRUE;
5127
      nodeperm = Ivector(1,noknots);
5128
      for(k=1;k<i;k++) nodeperm[k] = k;
5129
    }
5130
    if(activeperm) nodeperm[i] = j;
5131
  }
5132
  fclose(in);
5133

5134

5135
  /* Create inverse permutation for nodes */
5136
  if(activeperm) {
5137
    for(i=1;i<=noknots;i++) {
5138
      if(i==1) {
5139
	mini = nodeperm[i];
5140
	maxi = nodeperm[i];
5141
      }
5142
      else {
5143
	mini = MIN(nodeperm[i],mini);
5144
	maxi = MAX(nodeperm[i],maxi);
5145
      }
5146
    }
5147
    if(info) printf("LoadElmerInput: Node index range is: [%d %d]\n",mini,maxi);
5148
    invperm = Ivector(mini,maxi);
5149
    for(i=mini;i<=maxi;i++)
5150
      invperm[i] = -1;
5151
    for(i=1;i<=noknots;i++) {
5152
      j = nodeperm[i];
5153
      if( invperm[j] > 0 ) 
5154
	printf("LoadElmerInput: Node %d is redundant which may be problematic!\n",j);      
5155
      else
5156
	invperm[j] = i;
5157
    }
5158
  }
5159
  else {
5160
    mini = 1;
5161
    maxi = noknots;
5162
  }
5163
  
5164
  
5165
  activeelemperm = FALSE;
5166
  sprintf(filename,"%s","mesh.elements");
5167
  if ((in = fopen(filename,"r")) == NULL) {
5168
    printf("LoadElmerInput: The opening of the element-file %s failed!\n",
5169
	   filename);
5170
    bigerror("Cannot continue without element file!\n");
5171
  }
5172
  else 
5173
    if(info) printf("Loading %d bulk elements from %s\n",noelements,filename);
5174
  
5175
  for(i=1; i <= noelements; i++) {
5176
    iostat = fscanf(in,"%d",&j);
5177
    if(iostat <= 0 ) {
5178
      printf("LoadElmerInput: Failed reading element line %d, reducing size of element table to %d!\n",i,i-1);
5179
      data->noelements = noelements = i-1;
5180
      break;
5181
    }
5182
    
5183
    if(i != j && !activeelemperm) {
5184
      printf("LoadElmerInput: The element numbering (%d) at element %d is not compact, creating permutation\n",j,i);
5185
      activeelemperm = TRUE;
5186
      elemperm = Ivector(1,noelements0);
5187
      for(k=1; k < i; k++)
5188
	elemperm[k] = k;
5189
    }
5190
    if( activeelemperm ) elemperm[i] = j;
5191
    iostat = fscanf(in,"%d %d",&(data->material[i]),&elementtype);
5192
    if( iostat < 2 ) {
5193
      printf("LoadElmerInput: Failed reading definitions for bulk element %d\n",j);
5194
      bigerror("Cannot continue without this data!\n");
5195
    }
5196
    if(elementtype > maxelemtype ) {
5197
      printf("Invalid bulk elementtype: %d\n",elementtype);
5198
      bigerror("Cannot continue with invalid elements");
5199
    }
5200
    data->elementtypes[i] = elementtype;
5201
    nonodes = elementtype % 100;
5202
    if( nonodes > maxnodes ) {
5203
      printf("Number of nodes %d in element %d is greater than allocated maximum %d\n",nonodes,j,maxnodes);
5204
      bigerror("Cannot continue with invalid elements");
5205
    }
5206
    for(k=0;k<nonodes;k++) {
5207
      iostat = fscanf(in,"%d",&l);
5208
      if( l < mini || l > maxi ) {
5209
	printf("Node %d in element %d is out of range: %d\n",k+1,j,l);
5210
	bigerror("Cannot continue with this node numbering");
5211
      }
5212
      if( activeperm )
5213
	data->topology[i][k] = invperm[l];
5214
      else
5215
	data->topology[i][k] = l;
5216
    }
5217
  }
5218
  fclose(in);
5219
 
5220

5221
  /* Create inverse permutation for bulk elements */
5222
  if(activeelemperm) {
5223
    for(i=1;i<=noelements;i++) {
5224
      if(i==1) {
5225
	minelem = elemperm[i];
5226
	maxelem = elemperm[i];
5227
      }
5228
      else {
5229
	minelem = MIN(elemperm[i],minelem);
5230
	maxelem = MAX(elemperm[i],maxelem);
5231
      }
5232
    }
5233
    if(info) printf("LoadElmerInput: Element index range is: [%d %d]\n",minelem,maxelem);
5234
    invelemperm = Ivector(minelem,maxelem);
5235
    for(i=minelem;i<=maxelem;i++)
5236
      invelemperm[i] = -1;
5237
    for(i=1;i<=noelements;i++) {
5238
      j = elemperm[i];
5239
      if( invelemperm[j] > 0 )
5240
	printf("LoadElmerInput: Element %d is redundant which may be problematic!\n",j);      
5241
      else	
5242
	invelemperm[j] = i;
5243
    }
5244
  }
5245
  else {
5246
    minelem = 1;
5247
    maxelem = noelements;
5248
  }
5249

5250

5251

5252
  falseparents = 0;
5253
  noparents = 0;
5254
  bctopocreated = FALSE;
5255

5256
  sprintf(filename,"%s","mesh.boundary");
5257
  if ((in = fopen(filename,"r")) == NULL) {
5258
    printf("LoadElmerInput: The opening of the boundary-file %s failed!\n",
5259
	   filename);
5260
    return(4);
5261
  }
5262
  else {
5263
    if(info) printf("Loading %d boundary elements from %s\n",nosides,filename);
5264
  }
5265

5266
  if( nosides > 0 ) {
5267
    AllocateBoundary(bound,nosides);
5268
    data->noboundaries = 1;
5269
  };
5270

5271
  i = 0;
5272
  for(k=1; k <= nosides; k++) {
5273
    
5274
    iostat = fscanf(in,"%d",&dummyint);
5275
    if( iostat < 1 ) {
5276
      printf("LoadElmerInput: Failed reading boundary element line %d, reducing size of element table to %d!\n",k,i);
5277
      bound->nosides = nosides = i;
5278
      break;
5279
    }      
5280
    i++;
5281

5282
    iostat = fscanf(in,"%d %d %d %d",&(bound->types[i]),&p1,&p2,&elementtype);
5283
    if(iostat < 4 ) {
5284
      printf("LoadElmerInput: Failed reading definitions for boundary element %d\n",k);
5285
      bigerror("Cannot continue without this data!\n"); 
5286
    }    
5287
    if( p1 > 0 && (p1 < minelem || p1 > maxelem ) ) {
5288
      printf("Parent in boundary element %d out of range: %d\n",k,p1);    
5289
      bigerror("Cannot continue with bad parents");
5290
    }
5291
    if( p2 > 0 && (p2 < minelem || p2 > maxelem ) ) {
5292
      printf("Parent in boundary element %d out of range: %d\n",k,p2);
5293
      bigerror("Cannot continue with bad parents");
5294
    }
5295
      
5296
    if(activeelemperm) {
5297
      if( p1 > 0 ) p1 = invelemperm[p1];
5298
      if( p2 > 0 ) p2 = invelemperm[p2];
5299
    }
5300
    
5301
    if(elementtype > maxelemtype ) {
5302
      printf("Invalid boundary elementtype: %d\n",elementtype);
5303
      bigerror("Cannot continue with invalid elements");
5304
    }
5305
    nonodes = elementtype % 100;
5306
    if( nonodes > maxnodes ) {
5307
      printf("Number of nodes %d in side element %d is greater than allocated maximum %d\n",nonodes,dummyint,maxnodes);
5308
      bigerror("Cannot continue with invalid elements");
5309
    }
5310
    
5311
    for(j=0;j< nonodes ;j++) { 
5312
      iostat = fscanf(in,"%d",&l);
5313
      if(activeperm) 
5314
	sideind[j] = invperm[l];
5315
      else
5316
	sideind[j] = l;
5317
    }
5318
          
5319
    if( p1 == 0 && p2 != 0 ) {
5320
      bound->parent[i] = p2;
5321
      bound->parent2[i] = p1;
5322
    }
5323
    else {
5324
      bound->parent[i] = p1;
5325
      bound->parent2[i] = p2;
5326
    }
5327
    
5328
    if(bound->parent[i] > 0) {
5329
      fail = FindParentSide(data,bound,i,elementtype,sideind);
5330
      if(fail) falseparents++;      
5331
    }
5332
    else {
5333
#if 0
5334
      printf("Parents not specified for side %d with inds: ",dummyint);
5335
      for(j=0;j< elementtype%100 ;j++) 
5336
	printf("%d ",sideind[j]);
5337
      printf("and type: %d\n",bound->types[i]);   
5338
#endif
5339
      if( !bctopocreated ) {
5340
	bound->elementtypes = Ivector(1,nosides);
5341
	for(j=1;j<=nosides;j++)
5342
	  bound->elementtypes[j] = 0;
5343
	bound->topology = Imatrix(1,nosides,0,data->maxnodes-1);
5344
	bctopocreated = TRUE;
5345
      }
5346
      for(j=0;j< elementtype%100 ;j++) 
5347
	bound->topology[i][j] = sideind[j];
5348
      bound->elementtypes[i] = elementtype;
5349

5350
      printf("elementtype = %d %d %d\n",i,elementtype,sideind[0]);
5351
      noparents++;
5352
    }
5353
  }
5354
  
5355
  if( falseparents ) {
5356
    printf("There seems to be %d false parents in the mesh\n",falseparents);
5357
  }
5358
  if( noparents ) {
5359
    printf("There seems to be %d bc elements without parents in the mesh\n",noparents);
5360
  }
5361

5362
  bound->nosides = i;
5363
  fclose(in); 
5364
  
5365
  /* Save node permutation for later use */
5366
  data->nodepermexist = activeperm;
5367
  if(activeperm) {
5368
    data->nodeperm = nodeperm;
5369
    free_Ivector(invperm,mini,maxi);
5370
  }
5371
  
5372
  /* Element permutation is irrelevant probably for practical purposes (?) and hence it is forgotten. */
5373
  if(activeelemperm) {
5374
    free_Ivector(invelemperm,minelem,maxelem);
5375
    free_Ivector(elemperm,1,noelements0);
5376
  }
5377
  
5378

5379

5380
  sprintf(filename,"%s","mesh.names");
5381
  if ((in = fopen(filename,"r") )) {
5382
    int isbody,started,nameproblem;
5383
    
5384
    isbody = TRUE;
5385
    nameproblem = FALSE;
5386

5387
    if( nonames ) {
5388
      printf("Ignoring > mesh.names < because it was explicitly requested!\n");
5389
      goto namesend;
5390
    }
5391

5392
    if(info) printf("Loading names for mesh parts from file %s\n",filename);
5393

5394
    for(;;) {
5395
      if(Getnamerow(line,in,FALSE)) goto namesend;
5396

5397
      if(strstr(line,"names for boundaries")) {
5398
	if(info) printf("Reading names for mesh boundaries\n");
5399
	isbody = FALSE;
5400
	continue;
5401
      }
5402
      else if(strstr(line,"names for bodies")) {
5403
	if(info) printf("Reading names for mesh bodies\n");	
5404
	isbody = TRUE;
5405
	continue;
5406
      }
5407

5408
      /* get position for entity name */
5409
      ptr1 = strchr( line,'$');
5410
      if(!ptr1) continue;
5411
      ptr1++;
5412

5413
      /* get position for entity index and read it */
5414
      ptr2 = strchr( line,'=');
5415
      if(!ptr2) continue;
5416
      ptr2++;      
5417
      j = next_int(&ptr2);
5418

5419
      /* Initialize the entity name by white spaces */
5420
      for(i=0;i<MAXLINESIZE;i++) 
5421
	line2[i] = ' ';
5422

5423
      started = FALSE;
5424
      k = 0;
5425
      for(i=0;i<MAXLINESIZE;i++) {
5426
	if( ptr1[0] == '=' ) {
5427
	  /* remove possible trailing white space */
5428
	  if(line2[k-1] == ' ') k--;	    
5429
	  line2[k] = '\0';
5430
	  break;
5431
	}
5432
	if(started || ptr1[0] != ' ') {
5433
	  /* remove possible leading white space */
5434
	  line2[k] = ptr1[0];
5435
	  started = TRUE;
5436
	  k++;
5437
	}
5438
	ptr1++;
5439
      }
5440

5441
      /* Copy the entityname to mesh structure */
5442
      if( isbody ) {
5443
	if(j < 0 || j >= MAXBODIES ) {
5444
	  if(!nameproblem) printf("Cannot treat names for bodies beyond %d\n",j);
5445
	  nameproblem++;
5446
	}
5447
	else {
5448
	  if(!data->bodyname[j]) data->bodyname[j] = Cvector(0,MAXNAMESIZE);
5449
	  strcpy(data->bodyname[j],line2);	
5450
	  data->bodynamesexist = TRUE;
5451
	}
5452
      }
5453
      else {
5454
	if(j < 0 || j >= MAXBCS ) {
5455
	  if(!nameproblem) printf("Cannot treat names for boundaries beyond %d\n",j);
5456
	  nameproblem++;
5457
	}
5458
	else {
5459
	  if(!data->boundaryname[j]) data->boundaryname[j] = Cvector(0,MAXNAMESIZE);
5460
	  strcpy(data->boundaryname[j],line2);	
5461
	  data->boundarynamesexist = TRUE;
5462
	}
5463
      }
5464
    }
5465
    namesend:
5466

5467
    if( nameproblem ) {
5468
      data->boundarynamesexist = FALSE;
5469
      data->bodynamesexist = FALSE;
5470
      printf("Number of indexes beyond range: %d\n",nameproblem);
5471
      smallerror("Omitting use of names because some indexes are beyond range, code some more...");
5472
    }
5473
  }
5474

5475

5476
  if(!cdstat) cdstat = chdir("..");
5477

5478
  if(info) printf("Elmer mesh loaded successfully\n");
5479

5480
  return(0);
5481
}
5482

5483

5484
int SaveElmerInput(struct FemType *data,struct BoundaryType *bound,
5485
		   char *prefix,int decimals,int binary,int nooverwrite, int info)
5486
/* Saves the mesh in a form that may be used as input 
5487
   in Elmer calculations. 
5488
   */
5489
{
5490
  int noknots,noelements,material,sumsides,elemtype,fail,cdstat,bcdim;
5491
  int sideelemtype,conelemtype,nodesd1,nodesd2,newtype,singleprec;
5492
  int i,j,k,l,bulktypes[MAXELEMENTTYPE+1],sidetypes[MAXELEMENTTYPE+1];
5493
  int alltypes[MAXELEMENTTYPE+1],tottypes;
5494
  int ind[MAXNODESD1],ind2[MAXNODESD1],usedbody[MAXBODIES],usedbc[MAXBCS];
5495
  FILE *out;
5496
  char filename[MAXFILESIZE], outstyle[MAXFILESIZE];
5497
  char directoryname[MAXFILESIZE];
5498
  double coords[3];
5499
  float scoords[3];
5500
  
5501
  if(!data->created) {
5502
    printf("You tried to save points that were never created.\n");
5503
    return(1);
5504
  }
5505

5506
  singleprec = FALSE;
5507
  if(binary == 2) singleprec = TRUE;
5508
  
5509
  noelements = data->noelements;
5510
  noknots = data->noknots;
5511
  sumsides = 0;
5512

5513
  for(i=0;i<=MAXELEMENTTYPE;i++)
5514
    alltypes[i] = bulktypes[i] = sidetypes[i] = 0;
5515

5516
  for(i=0;i<MAXBODIES;i++)
5517
    usedbody[i] = 0;
5518
  for(i=0;i<MAXBCS;i++)
5519
    usedbc[i] = 0;
5520

5521
  sprintf(directoryname,"%s",prefix);
5522

5523
  if(info) printf("Saving mesh in ElmerSolver format to directory %s.\n",
5524
		  directoryname);
5525

5526
  fail = chdir(directoryname);
5527
  if(fail) {
5528
#ifdef MINGW32
5529
    fail = mkdir(directoryname);
5530
#else
5531
    fail = mkdir(directoryname,0750);
5532
#endif
5533
    if(fail) {
5534
      printf("Could not create a result directory!\n");
5535
      return(1);
5536
    }
5537
    else {
5538
      cdstat = chdir(directoryname);
5539
    }
5540
  }
5541
  else {
5542
    if(info) printf("Reusing an existing directory\n");
5543
    if(nooverwrite) {
5544
      if ((out = fopen("mesh.header", "r"))) {
5545
	printf("Mesh seems to already exist, writing is cancelled!\n"); 
5546
	return(1);
5547
      }
5548
    }
5549
  }
5550

5551
  /* Save mesh nodes in binary or ascii format */  
5552
  if(binary) {  
5553
    if(singleprec) 
5554
      sprintf(filename,"%s","mesh.nodes.sbin");
5555
    else 
5556
      sprintf(filename,"%s","mesh.nodes.bin");
5557
    out = fopen(filename,"wb");
5558
  } else {
5559
    sprintf(filename,"%s","mesh.nodes");
5560
    out = fopen(filename,"w");
5561
  }
5562
    
5563
  if(info) printf("Saving %d coordinates to %s.\n",noknots,filename);  
5564
  if(out == NULL) {
5565
    printf("opening of file was not successful\n");
5566
    return(2);
5567
  }
5568
  
5569
  if(binary) {
5570
    for(i=1; i <= noknots; i++) {
5571
      fwrite(&i,sizeof(int),1,out); 
5572
      /* Note that in binary format we don't save the obsolite "-1". */
5573
      if(singleprec) {
5574
	/* We save the coordinates in single precidion format */
5575
	scoords[0] = data->x[i];
5576
	scoords[1] = data->y[i];
5577
	scoords[2] = data->z[i];
5578
	fwrite(scoords,sizeof(float),3,out); 
5579
      } else {
5580
	coords[0] = data->x[i];
5581
	coords[1] = data->y[i];
5582
	coords[2] = data->z[i];
5583
	fwrite(coords,sizeof(double),3,out);
5584
      }
5585
    }
5586
  }
5587
  else {
5588
    sprintf(outstyle,"%%d %%d %%.%dg %%.%dg %%.%dg\n",decimals,decimals,decimals);
5589
    for(i=1; i <= noknots; i++) 
5590
      fprintf(out,outstyle,i,-1,data->x[i],data->y[i],data->z[i]);    
5591
  }    
5592
  fclose(out);
5593

5594

5595
  /* Save bulk elements in binary or ascii format */  
5596
  if(binary) { 
5597
    sprintf(filename,"%s","mesh.elements.bin");
5598
    out = fopen(filename,"wb");
5599
  } else {
5600
    sprintf(filename,"%s","mesh.elements");
5601
    out = fopen(filename,"w");
5602
  }
5603
  
5604
  if(info) printf("Saving %d element topologies to %s.\n",noelements,filename);
5605
  if(out == NULL) {
5606
    printf("opening of file was not successful\n");
5607
    return(3);
5608
  }
5609
  
5610
  for(i=1;i<=noelements;i++) {
5611
    elemtype = data->elementtypes[i];
5612
    material = data->material[i];
5613

5614
    if(material < MAXBODIES) usedbody[material] += 1;
5615
    bulktypes[elemtype] += 1;
5616
    nodesd2 = elemtype%100;
5617

5618
    if(binary) {      
5619
      fwrite(&i,sizeof(int),1,out);
5620
      k=-1; fwrite(&k,sizeof(int),1,out);  /* partition is saved in binary also in serial as "-1". */
5621
      fwrite(&material,sizeof(int),1,out);
5622
      fwrite(&elemtype,sizeof(int),1,out);
5623
      fwrite(&data->topology[i][0],sizeof(int),nodesd2,out);
5624
    } else {      
5625
      fprintf(out,"%d %d %d",i,material,elemtype);
5626
      for(j=0;j < nodesd2;j++) 
5627
	fprintf(out," %d",data->topology[i][j]);
5628
      fprintf(out,"\n");
5629
    }
5630
  }
5631
  fclose(out);
5632

5633

5634
  /* Save bulk elements in binary or ascii format */  
5635
  if(binary) { 
5636
    sprintf(filename,"%s","mesh.boundary.bin");
5637
    out = fopen(filename,"wb");
5638
  } else {
5639
    sprintf(filename,"%s","mesh.boundary");
5640
    out = fopen(filename,"w");
5641
  }
5642
  
5643
  if(info) printf("Saving boundary elements to %s.\n",filename);
5644
  if(out == NULL) {
5645
    printf("opening of file was not successful\n");
5646
    return(4);
5647
  }
5648

5649
  sumsides = 0;
5650

5651

5652
  /* Save normal boundaries */
5653
  for(j=0;j < MAXBOUNDARIES;j++) {
5654
    if(bound[j].created == FALSE) continue;
5655
    if(bound[j].nosides == 0) continue;
5656
    
5657
    for(i=1; i <= bound[j].nosides; i++) {
5658
      GetBoundaryElement(i,&bound[j],data,ind,&sideelemtype); 
5659
      sumsides++;
5660

5661
      k = bound[j].types[i];
5662
      if(k < MAXBCS) {	
5663
	bcdim = GetElementDimension(sideelemtype);
5664
	usedbc[k] = MAX(usedbc[k],bcdim+1);
5665
      }
5666
	
5667
      sidetypes[sideelemtype] += 1;
5668
      nodesd1 = sideelemtype%100;
5669
      
5670
      if(binary) {
5671
	fwrite(&sumsides,sizeof(int),1,out);	
5672
	k=-1; fwrite(&k,sizeof(int),1,out);
5673
	fwrite(&bound[j].types[i],sizeof(int),1,out);
5674
	fwrite(&bound[j].parent[i],sizeof(int),1,out);
5675
	fwrite(&bound[j].parent2[i],sizeof(int),1,out);
5676
	fwrite(&sideelemtype,sizeof(int),1,out);
5677
	fwrite(ind,sizeof(int),nodesd1,out);
5678
      } else {
5679
	fprintf(out,"%d %d %d %d %d",
5680
		sumsides,bound[j].types[i],bound[j].parent[i],bound[j].parent2[i],sideelemtype);
5681
	for(l=0;l<nodesd1;l++)
5682
	  fprintf(out," %d",ind[l]);
5683
	fprintf(out,"\n");
5684
      }
5685
    }
5686
  }
5687

5688
  newtype = 0;
5689
  for(j=0;j < MAXBOUNDARIES;j++) {
5690
    if(bound[j].created == FALSE) continue;
5691
    for(i=1; i <= bound[j].nosides; i++) 
5692
      newtype = MAX(newtype, bound[j].types[i]);
5693
  }  
5694
  fclose(out);
5695

5696
  tottypes = 0;
5697
  for(i=0;i<=MAXELEMENTTYPE;i++) {
5698
    alltypes[i] = bulktypes[i] + sidetypes[i];
5699
    if(alltypes[i]) tottypes++;
5700
  }
5701

5702
  sprintf(filename,"%s","mesh.header");
5703
  out = fopen(filename,"w");
5704
  if(info) printf("Saving header info to %s.\n",filename);  
5705
  if(out == NULL) {
5706
    printf("opening of file was not successful\n");
5707
    return(4);
5708
  }
5709
  fprintf(out,"%-6d %-6d %-6d\n",
5710
	  noknots,noelements,sumsides);
5711
  fprintf(out,"%-6d\n",tottypes);
5712
  for(i=0;i<=MAXELEMENTTYPE;i++) {
5713
    if(alltypes[i]) 
5714
      fprintf(out,"%-6d %-6d\n",i,bulktypes[i]+sidetypes[i]);
5715
  }
5716
  fclose(out);
5717

5718

5719
  if(data->boundarynamesexist || data->bodynamesexist) {
5720
    sprintf(filename,"%s","mesh.names");
5721
    out = fopen(filename,"w");
5722
    if(info) printf("Saving names info to %s.\n",filename);  
5723
    if(out == NULL) {
5724
      printf("opening of file was not successful\n");
5725
      return(5);
5726
    }
5727
    
5728
    if(data->bodynamesexist) {
5729
      fprintf(out,"! ----- names for bodies -----\n");
5730
      for(i=1;i<MAXBODIES;i++)
5731
	if(usedbody[i]) {
5732
	  if(data->bodyname[i])
5733
	    fprintf(out,"$ %s = %d\n",data->bodyname[i],i);
5734
	  else
5735
	    fprintf(out,"$ body%d = %d\n",i,i);	    
5736
	}
5737
    }     
5738
    if(data->boundarynamesexist) {
5739
      fprintf(out,"! ----- names for boundaries -----\n");
5740
      for(i=1;i<MAXBCS;i++) 	
5741
	if(usedbc[i]) {
5742
	  bcdim = usedbc[i]-1;
5743
	  if(data->boundaryname[i]) 
5744
	    fprintf(out,"$ %s = %d\n",data->boundaryname[i],i);
5745
	  else if(bcdim == 2) { 
5746
	    fprintf(out,"$ surf_bc%d = %d\n",i,i);
5747
	  }
5748
	  else if(bcdim == 1) {
5749
	    fprintf(out,"$ line_bc%d = %d\n",i,i);
5750
	  }
5751
	  else if(bcdim == 0) {
5752
	    fprintf(out,"$ node_bc%d = %d\n",i,i);
5753
	  }
5754
	}
5755
    }
5756
    fclose(out);
5757

5758
    sprintf(filename,"%s","entities.sif");
5759
    out = fopen(filename,"w");
5760
    if(info) printf("Saving entities info to %s.\n",filename);  
5761
    if(out == NULL) {
5762
      printf("opening of file was not successful\n");
5763
      return(5);
5764
    }
5765

5766
    if(data->bodynamesexist) {
5767
      fprintf(out,"!------ Skeleton for body section -----\n");
5768
      j = 0;
5769
      for(i=1;i<MAXBODIES;i++) {
5770
	if(usedbody[i]) {
5771
	  j = j + 1;
5772
	  fprintf(out,"Body %d\n",j);
5773
	  if(data->bodyname[i])
5774
	    fprintf(out,"  Name = %s\n",data->bodyname[i]);
5775
	  else
5776
	    fprintf(out,"  Name = body%d\n",i);	    
5777
	  fprintf(out,"End\n\n");
5778
	}
5779
      }
5780
    }
5781

5782
    if(data->boundarynamesexist) {
5783
      fprintf(out,"!------ Skeleton for boundary section -----\n");
5784
      j = 0;
5785
      for(i=1;i<MAXBCS;i++) {
5786
	if(usedbc[i]) {
5787
	  j = j + 1;
5788
	  fprintf(out,"Boundary Condition %d\n",j);
5789
	  if(data->boundaryname[i])
5790
	    fprintf(out,"  Name = %s\n",data->boundaryname[i]);
5791
	  else
5792
	    fprintf(out,"  Name = bc%d\n",i);	    
5793
	  fprintf(out,"End\n\n");
5794
	}
5795
      }
5796
    }
5797
    fclose(out);
5798
  }
5799
  
5800
  if(data->nodepermexist) {
5801
    sprintf(filename,"%s","mesh.nodeperm");
5802
    out = fopen(filename,"w");
5803

5804
    if(info) printf("Saving initial node permutation to %s.\n",filename);  
5805
    if(out == NULL) {
5806
      printf("opening of file was not successful\n");
5807
      return(3);
5808
    }
5809
    for(i=1; i <= noknots; i++) 
5810
      fprintf(out,"%d %d\n",i,data->nodeperm[i]);
5811
  }
5812

5813

5814
  cdstat = chdir("..");
5815
  
5816
  return(0);
5817
}
5818

5819

5820
int CreateElmerGridMesh(struct GridType *grid,
5821
			struct FemType *data,struct BoundaryType *boundaries,
5822
			Real relh,int info) {
5823
  int j;  
5824
  struct CellType *cell;
5825
  
5826
  for(j=0;j<MAXBOUNDARIES;j++) {
5827
    boundaries[j].created = FALSE;
5828
    boundaries[j].nosides = FALSE;
5829
  }
5830

5831
  SetElementDivision(grid,relh,info);
5832
  
5833
  cell = (struct CellType*)
5834
    malloc((size_t) (grid->nocells+1)*sizeof(struct CellType)); 
5835
  SetCellData(grid,cell,info);
5836

5837
  if(grid->dimension == 1) 
5838
    SetCellKnots1D(grid,cell,info);
5839
  else
5840
    SetCellKnots(grid,cell,info);
5841

5842
  CreateKnots(grid,cell,data,0,0);
5843

5844
  if(grid->noboundaries > 0) {
5845
    for(j=0;j<grid->noboundaries;j++) {
5846
      if(grid->boundsolid[j] < 4) {
5847
	CreateBoundary(cell,data,&(boundaries[j]),grid->boundext[j],grid->boundint[j],
5848
		       1,grid->boundtype[j],info);  
5849
      } 
5850
      else { 
5851
	CreatePoints(cell,data,&(boundaries[j]),grid->boundext[j],grid->boundint[j],
5852
		     grid->boundsolid[j],grid->boundtype[j],info); 	    
5853
      }
5854
    }
5855
  }
5856

5857
  free(cell);
5858

5859
  return 0;
5860
}
5861

5862

5863
Product

Resources

Company
