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("****************** Elmergrid ************************\n");
106
  printf("This program can create simple 2D structured meshes consisting of\n");
107
  printf("linear, quadratic or cubic rectangles or triangles. The meshes may\n");
108
  printf("also be extruded and revolved to create 3D forms. In addition many\n");
109
  printf("mesh formats may be imported into Elmer software. Some options have\n");
110
  printf("not been properly tested. Contact the author if you face problems.\n\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("B) Inline mode which expects at least three input parameters\n");
117
  printf("   'ElmerGrid 1 3 test'\n\n");
118
  printf("The first parameter defines the input file format:\n");
119
  printf("1)  .grd      : ElmerGrid file format\n");
120
  printf("2)  .mesh.*   : Elmer input format\n");
121
  printf("3)  .ep       : Elmer output format\n");
122
  printf("4)  .ansys    : Ansys input format\n");
123
  printf("5)  .inp      : Abaqus input format by Ideas\n");
124
  printf("6)  .fil      : Abaqus output format\n");
125
  printf("7)  .FDNEUT   : Gambit (Fidap) neutral file\n");
126
  printf("8)  .unv      : Universal mesh file format\n");
127
  printf("9)  .mphtxt   : Comsol Multiphysics mesh format\n");
128
  printf("10) .dat      : Fieldview format\n");
129
  printf("11) .node,.ele: Triangle 2D mesh format\n");
130
  printf("12) .mesh     : Medit mesh format\n");
131
  printf("13) .msh      : GID mesh format\n");
132
  printf("14) .msh      : Gmsh mesh format\n");
133
  printf("15) .ep.i     : Partitioned ElmerPost format\n");
134
  printf("16) .2dm      : 2D triangular FVCOM format\n");
135
#if 0
136
  printf("17) .msh      : Nastran format\n");
137
  printf("18) .msh      : CGsim format\n");
138
  printf("19) .geo      : Geo format\n");
139
  printf("20) .tra      : Cedrat Flux 2D format\n");
140
  printf("21) .pf3      : Cedrat Flux 3D format\n");
141
#endif 
142

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

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

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

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

209
  printf("\nKeywords are related to mesh partitioning for parallel ElmerSolver runs:\n");
210
  printf("-partition int[3]    : the mesh will be partitioned in cartesian main directions\n");
211
  printf("-partorder real[3]   : in the 'partition' method set the direction of the ordering\n");
212
  printf("-parttol real        : in the 'partition' method set the tolerance for ordering\n");
213
  printf("-partcell int[3]     : the mesh will be partitioned in cells of fixed sizes\n");
214
  printf("-partcyl int[3]      : the mesh will be partitioned in cylindrical main directions\n");
215
#if USE_METIS
216
  printf("-metiskway int       : mesh will be partitioned with Metis using graph Kway routine\n");
217
  printf("-metisrec int        : mesh will be partitioned with Metis using graph Recursive routine\n");
218
  printf("-metiscontig         : enforce that the metis partitions are contiguous\n");
219
  printf("-metisvol            : minimize total communication volume in Metis\n");
220
  printf("-metisminconn        : minimize the maximum connectivity count in Metis\n");
221
  printf("-metisseed int       : random number generator seed for Metis algorithms\n");
222
  printf("-metisncuts int      : number of competing partitions to generate\n");
223
#endif
224
  printf("-partdual            : use the dual graph in partition method (when available)\n");
225
  printf("-halo                : create halo for the partitioning for DG\n");
226
  printf("-halobc              : create halo for the partitioning at boundaries only\n");
227
  printf("-haloz / -halor      : create halo for the the special z- or r-partitioning\n");
228
  printf("-halogreedy          : create halo being greedy over the partition interfaces\n");
229
  printf("-indirect            : create indirect connections (102 elements) in the partitioning\n");
230
  printf("-periodic int[3]     : periodic coordinate directions for parallel & conforming meshes\n");
231
  printf("-partoptim           : apply aggressive optimization to node sharing\n");
232
  printf("-partnobcoptim       : do not apply optimization to bc ownership sharing\n");
233
  printf("-partbw              : minimize the bandwidth of partition-partition couplings\n");
234
  printf("-parthypre           : number the nodes continuously partitionwise\n");
235
  printf("-partzbc             : partition connected BCs separately to partitions in Z-direction\n");
236
  printf("-partrbc             : partition connected BCs separately to partitions in R-direction\n");
237
#if USE_METIS
238
  printf("-metisbc             : partition connected BCs separately to partitions by Metis\n");
239
#endif
240
  printf("-partlayers int      : extend boundary partitioning by element layers\n");
241

242
  printf("\nKeywords are related to (nearly obsolete) ElmerPost format:\n");
243
  printf("-partjoin int        : number of ElmerPost partitions in the data to be joined\n");
244
  printf("-saveinterval int[3] : the first, last and step for fusing parallel data\n");
245
  printf("-nobound             : disable saving of boundary elements in ElmerPost format\n");
246

247
  if(0) printf("-names               : conserve name information where applicable\n");
248
}
249
 
250

251
void Goodbye()
252
{
253
  printf("\nThank you for using Elmergrid!\n");
254
  printf("Send bug reports and feature wishes to [email protected]\n");
255
  exit(0);
256
}
257

258

259

260
#if USE_MATC
261
char *mtc_domath(const char *);
262
void mtc_init(FILE * input, FILE * output, FILE * error);
263
#endif
264

265

266
static int Getline(char *line1,FILE *io) 
267
{
268
  int i,isend;
269
  char line0[MAXLINESIZE],*charend,*matcpntr,*matcpntr0;
270

271
  for(i=0;i<MAXLINESIZE;i++)
272
    line0[i] = 0x00;
273

274
 newline:
275

276
  charend = fgets(line0,MAXLINESIZE,io);
277
  isend = (charend == NULL);
278

279
  if(isend) return(1);
280

281
  if(line0[0] == '#' || line0[0] == '%' || line0[0] == '!') goto newline;
282
  if(!matcactive && line0[0] == '*') goto newline;
283

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

295
#if USE_MATC
296
  if(matcactive) {
297
    matcpntr0 = strchr(line0,'$');
298
    if(matcpntr0) {
299
      matcpntr = mtc_domath(&matcpntr0[1]);
300
      if(matcpntr) {
301
	strcpy(matcpntr0, matcpntr);
302
	if(0) printf("A: %s\n%s\n",matcpntr0,matcpntr);
303
      }
304
    }
305
  }
306
#endif 
307

308
  if(strstr(line0,"subcell boundaries")) goto newline;
309
  if(strstr(line0,"material structure")) goto newline;
310
  if(strstr(line0,"mode")) goto newline;
311
  if(strstr(line0,"type")) goto newline;
312

313
  for(i=0;i<MAXLINESIZE;i++) 
314
    line1[i] = toupper(line0[i]);
315

316
  if(iodebug) printf("line: %s\n",line1);
317

318
  return(0);
319
}
320

321

322

323
void InitGrid(struct GridType *grid)
324
/* Initializes the grid of a specific mesh. A grid can differ 
325
   between various differential equations. 
326
   */
327
#define MAXBODYID 1000
328
{
329
  int i,j;
330

331
  grid->layered = FALSE;
332
  grid->layeredbc = TRUE;
333
  grid->layerbcoffset = 0;
334
  grid->triangles = FALSE;
335
  grid->triangleangle = 0.0;
336
  grid->partitions = FALSE;
337
  grid->wantedelems = 0;
338
  grid->wantedelems3d = 0;
339
  grid->wantednodes3d = 0;
340
  grid->limitdx = 0.0;
341
  grid->limitdxverify = FALSE;
342

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

369
  grid->rotatecurve = FALSE;
370
  grid->curverad = 0.5;
371
  grid->curveangle = 90.0;
372
  grid->curvezet = 1.0;
373

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

402
  grid->zmaterialmapexists = FALSE;
403
  grid->zhelicityexists = FALSE;
404
  grid->zhelicity = 0.0;
405

406
  /* Initializes the numbering of the cells. */
407
  for(j=0;j<=MAXCELLS+1;j++)
408
    for(i=0;i<=MAXCELLS+1;i++) 
409
      grid->structure[j][i] = MAT_NOTHING; 
410

411
  for(j=0;j<=MAXCELLS+1;j++)
412
    for(i=0;i<=MAXCELLS+1;i++) 
413
      grid->numbered[j][i] = 0; 
414

415
  grid->noboundaries = 0;
416
  for(i=0;i<MAXBOUNDARIES;i++) {
417
    grid->boundint[i] = 0;
418
    grid->boundext[i] = 0;
419
    grid->boundsolid[i] = 0;
420
    grid->boundtype[i] = 0;
421
  }
422

423
  grid->mappings = 0;
424
  for(i=0;i<MAXMAPPINGS;i++) {
425
    grid->mappingtype[i] = 0;
426
    grid->mappingline[i] = 0;    
427
    grid->mappinglimits[2*i] = 0;    
428
    grid->mappinglimits[2*i+1] = 0;    
429
    grid->mappingpoints[i] = 0;
430
  }
431
}
432

433

434
static void ExampleGrid1D(struct GridType **grids,int *nogrids,int info) 
435
/* Creates an example grid that might be used to analyze 
436
   flow through a step. */
437
{
438
  int j;
439
  struct GridType *grid;
440

441
  (*nogrids) = 1;
442
  (*grids) = (struct GridType*)
443
    malloc((size_t) (*nogrids)*sizeof(struct GridType)); 
444

445
  grid = grids[0];
446

447
  InitGrid(grid);
448

449
  grid->nonodes = 2;
450
  grid->xcells = 3;
451
  grid->ycells = 1;
452
  grid->wantedelems = 40;
453
  grid->firstmaterial = 1;
454
  grid->lastmaterial = 1;
455
  grid->autoratio = 1;
456
  grid->coordsystem = COORD_CART1;
457
  grid->numbering = NUMBER_1D;
458

459
  grid->x[0] = 0.0;
460
  grid->x[1] = 1.0;
461
  grid->x[2] = 3.0;
462
  grid->x[3] = 4.0;
463

464
  grid->xexpand[1] = 2.0;
465
  grid->xexpand[3] = 0.5;
466
  grid->xdens[2] = 0.5;
467

468
  for(j=1;j<=3;j++)
469
    grid->structure[1][j] = 1;
470

471
  grid->noboundaries = 2;
472
  grid->boundint[0] = 1;
473
  grid->boundint[1] = 1;
474

475
  grid->boundext[0] = -1;
476
  grid->boundext[1] = -2;
477

478
  grid->boundtype[0] = 1;
479
  grid->boundtype[1] = 2;
480

481
  grid->boundsolid[0] = 1;
482
  grid->boundsolid[1] = 1;
483

484
  grid->mappings = 0;
485
  
486
  if(info) printf("A simple 1D example mesh was created\n");
487
}
488

489

490
static void ExampleGrid2D(struct GridType **grids,int *nogrids,int info) 
491
/* Creates an example grid that might be used to analyze 
492
   flow through a step. */
493
{
494
  int j;
495
  struct GridType *grid;
496

497
  (*nogrids) = 1;
498
  (*grids) = (struct GridType*)
499
    malloc((size_t) (*nogrids)*sizeof(struct GridType)); 
500

501
  grid = grids[0];
502

503
  InitGrid(grid);
504

505
  grid->xcells = 4;
506
  grid->ycells = 3;
507
  grid->wantedelems = 100;
508
  grid->totxelems = grid->totyelems = grid->totzelems = 10;
509
  grid->firstmaterial = 1;
510
  grid->lastmaterial = 1;
511
  grid->autoratio = 1;
512
  grid->coordsystem = COORD_CART2;
513
  grid->numbering = NUMBER_XY;
514

515
  grid->x[0] = -1.0;
516
  grid->x[1] = 0.0;
517
  grid->x[2] = 2.0;
518
  grid->x[3] = 6.0;
519
  grid->x[4] = 7.0;
520

521
  grid->y[0] = 0.0;
522
  grid->y[1] = 0.5;
523
  grid->y[2] = 0.75;
524
  grid->y[3] = 1.0;
525

526
  grid->xexpand[2] = 0.4;
527
  grid->xexpand[3] = 5.0;
528

529
  grid->yexpand[2] = 2.0;
530
  grid->yexpand[3] = 0.5;
531

532
  for(j=1;j<=3;j++)
533
    grid->structure[j][1] = 2;
534
  for(j=1;j<=3;j++)
535
    grid->structure[j][2] = 1;
536
  grid->structure[1][2] = 0;
537
  for(j=1;j<=3;j++)
538
    grid->structure[j][3] = 1;
539
  for(j=1;j<=3;j++)
540
    grid->structure[j][4] = 3;
541

542
  grid->noboundaries = 3;
543
  grid->boundint[0] = 1;
544
  grid->boundint[1] = 1;
545
  grid->boundint[2] = 1;
546

547
  grid->boundext[0] = 2;
548
  grid->boundext[1] = 3;
549
  grid->boundext[2] = 0;
550

551
  grid->boundtype[0] = 1;
552
  grid->boundtype[1] = 2;
553
  grid->boundtype[2] = 3;
554

555
  grid->boundsolid[0] = 1;
556
  grid->boundsolid[1] = 1;
557
  grid->boundsolid[2] = 1;
558

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

581
  
582
static void ExampleGrid3D(struct GridType **grids,int *nogrids,int info) 
583
/* Creates an example grid that might be used to analyze 
584
   a simple acceleration sensor. */
585
{
586
  int i,j;
587
  struct GridType *grid;
588

589
  (*nogrids) = 1;
590
  (*grids) = (struct GridType*)
591
    malloc((size_t) (*nogrids)*sizeof(struct GridType)); 
592

593
  grid = grids[0];
594

595
  InitGrid(grid);
596

597
  grid->dimension = 3;
598
  grid->coordsystem = COORD_CART3;
599
  grid->layered = TRUE;
600

601
  grid->xcells = 4;
602
  grid->ycells = 5;
603
  grid->zcells = 5;
604

605
  grid->wantedelems = 1000;
606
  grid->firstmaterial = 1;
607
  grid->lastmaterial = 3;
608
  grid->autoratio = 1;
609
  grid->coordsystem = COORD_CART3;
610
  grid->numbering = NUMBER_XY;
611

612
  grid->x[0] = 0.0;
613
  grid->x[1] = 0.1;
614
  grid->x[2] = 0.4;
615
  grid->x[3] = 0.5;
616
  grid->x[4] = 0.6;
617

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

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

632
  for(j=1;j<=5;j++)
633
    for(i=1;i<=4;i++)
634
      grid->structure[j][i] = 1;
635

636
  grid->structure[2][2] = 3;
637
  grid->structure[2][3] = 3;
638
  grid->structure[3][3] = 3;
639
  grid->structure[4][3] = 3;
640
  grid->structure[4][2] = 3;
641

642
  grid->structure[3][2] = 2;
643

644
  grid->zlastmaterial[5] = 3;
645
  grid->zlastmaterial[2] = 1;
646
  grid->zlastmaterial[3] = 2;
647
  grid->zlastmaterial[4] = 1;
648
  grid->zlastmaterial[5] = 3;
649

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

656
  grid->noboundaries = 4;
657
  grid->boundint[0] = 1;
658
  grid->boundint[1] = 1;
659
  grid->boundint[2] = 1;
660
  grid->boundint[2] = 2;
661

662
  grid->boundext[0] = 0;
663
  grid->boundext[1] = 2;
664
  grid->boundext[2] = 3;
665
  grid->boundext[3] = 3;
666

667
  grid->boundtype[0] = 1;
668
  grid->boundtype[1] = 2;
669
  grid->boundtype[2] = 3;
670
  grid->boundtype[3] = 4;
671

672
  grid->boundsolid[0] = 1;
673
  grid->boundsolid[1] = 1;
674
  grid->boundsolid[2] = 1;
675
  grid->boundsolid[3] = 1;
676

677
  if(info) printf("A simple 3D example mesh was created\n");
678
}
679

680

681
void CreateExampleGrid(int dim,struct GridType **grids,int *nogrids,int info) 
682
{
683
  switch (dim) {
684
    
685
  case 1: 
686
    ExampleGrid1D(grids,nogrids,info);
687
    break;
688

689
  case 2: 
690
    ExampleGrid2D(grids,nogrids,info);
691
    break;
692

693
  case 3: 
694
    ExampleGrid3D(grids,nogrids,info);
695
    break;
696
  }
697
}
698

699

700

701

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

720
  if(0) printf("SetElementDivision\n");
721

722
  if(grid->dimension == 1) 
723
    grid->numbering = NUMBER_1D;
724

725
  nx = grid->xcells;
726
  ny = grid->ycells;
727
  linearlimit = 0.001;
728

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

749
  for(i=1; i <= ny ;i++) {
750
    if(fabs(1.-grid->yexpand[i]) < linearlimit) 
751
      grid->ylinear[i] = TRUE;
752
    else if(fabs(1.+grid->yexpand[i]) < linearlimit) 
753
      grid->ylinear[i] = TRUE;
754
    else 
755
      grid->ylinear[i] = FALSE;
756
  }
757

758

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

776
  sumxelems = 0;
777
  for(i=1;i<=nx;i++) {
778
    if(grid->dimension == 1 && !grid->numbered[1][i]) continue;
779
    sumxelems += grid->xelems[i];
780
  }
781
  sumyelems = 0;
782
  for(i=1;i<=ny;i++)  
783
    sumyelems += grid->yelems[i];
784

785
  if(grid->dimension == 1) {
786
    grid->autoratio = 2;
787
    grid->totxelems = grid->wantedelems;
788
    grid->totyelems = 1;
789
  }
790

791
  /* Allocate elements for both axis separately */
792
  if(grid->autoratio == 2) {
793

794
    wantedelemsx = (int) (1.0*grid->totxelems / relh);
795
    wantedelemsy = (int) (1.0*grid->totyelems / relh);
796

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

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

878

879
  /* Both axis dependently */
880
  if(grid->autoratio == 1)  {
881
    
882
    wantedelems = (int) (1.0*grid->wantedelems / (relh*relh));
883

884
    sumxyelems = 0;
885
    for(;;) {
886
      dxmax = 0.0;
887

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

917
	dx *= grid->xdens[i];
918
	if(dx > dxmax) {
919
	  dxmax = dx;
920
	  nxmax = i;
921
	}
922
      }      
923

924

925
      dymax = 0.0;
926
      for(i=1;i<=ny;i++) {
927

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

962
      if(dxmax > dymax) {
963
	grid->xelems[nxmax] += 1;
964
	sumxelems++;
965
	dxlimit = dxmax;
966
      }
967
      else {
968
	grid->yelems[nymax] += 1;
969
	sumyelems++;
970
	dxlimit = dymax;
971
      }
972
      
973
      sumxyelems = 0;
974
      for(j=1;j<=grid->ycells;j++)
975
	for(i=1;i<=grid->xcells;i++) 
976
	  if(grid->numbered[j][i]) sumxyelems += grid->xelems[i] * grid->yelems[j];
977

978
      if(wantedelems <= sumxyelems) break;
979
    }
980
  }    
981

982

983

984
  if(grid->autoratio == 3 || grid->limitdxverify)  {
985
    
986
    if(grid->autoratio == 3) {
987
      dxlimit = relh * grid->limitdx;
988
      dxmax = dymax = dxlimit;
989
    }
990

991
    for(i=1;i<=nx;i++) {
992

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

1030

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

1076

1077
  /* Put the linearity flags if there is only one division. */
1078
  for(i=1; i<= nx ;i++) 
1079
    if (grid->xelems[i] == 1) 
1080
      grid->xlinear[i] = TRUE; 
1081
  for(i=1; i <= ny ;i++) 
1082
    if(grid->yelems[i] == 1) 
1083
      grid->ylinear[i] = TRUE;
1084
 
1085

1086
  /* Set the size of the first element within each subcell */
1087
  for(i=1;i<=nx;i++) {
1088

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

1117
  for(i=1;i<=ny;i++) {
1118

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

1147
  /* Calculate the total number of elements */
1148
  sumxyelems = 0;
1149
  for(j=1;j<=grid->ycells;j++)
1150
    for(i=1;i<=grid->xcells;i++) 
1151
      if(grid->numbered[j][i]) sumxyelems += grid->xelems[i] * grid->yelems[j];
1152

1153
  grid->noelements = sumxyelems;
1154

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

1157
  grid->dx0 = dxmax;
1158
  grid->dy0 = dymax;
1159
}
1160

1161

1162

1163

1164
void SetCellData(struct GridType *grid,struct CellType *cell,int info)
1165
/* Sets the data that can directly be derived from type GridType. 
1166
 */
1167
{
1168
  int i,j,cnew=0;
1169

1170
  for(j=1;j<= grid->ycells ;j++)                   /* cells direction up    */
1171
    for(i=1;i<= grid->xcells; i++)                 /* cells direction right */      
1172
      if( (cnew = grid->numbered[j][i]) ) {          /* if cell is occupied   */
1173

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

1235
	/* Set the neighbour cell indices for sides and corners */ 
1236
	cell[cnew].neighbour[LEFT] = grid->numbered[j][i-1];
1237
	cell[cnew].neighbour[DOWN] = grid->numbered[j-1][i];
1238
	cell[cnew].neighbour[RIGHT]= grid->numbered[j][i+1];  
1239
	cell[cnew].neighbour[UP]   = grid->numbered[j+1][i];
1240
	cell[cnew].neighbour[4] = grid->numbered[j-1][i-1];
1241
	cell[cnew].neighbour[5] = grid->numbered[j-1][i+1];
1242
	cell[cnew].neighbour[6] = grid->numbered[j+1][i+1];
1243
	cell[cnew].neighbour[7] = grid->numbered[j+1][i-1];
1244
      }
1245

1246
  if(info) printf("%d cells were created.\n",grid->nocells);
1247
}
1248

1249

1250

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

1273
  numbering = grid->numbering;
1274
  nonodes = grid->nonodes;
1275
  knotno  = 0;
1276
  elemno  = 0;
1277

1278
  if(numbering == NUMBER_XY) {
1279
    xcells = grid->xcells;
1280
    ycells = grid->ycells;
1281
    xelems = grid->xelems;
1282
    yelems = grid->yelems;
1283
  }
1284
  else if(numbering == NUMBER_YX) {
1285
    xcells = grid->ycells;
1286
    ycells = grid->xcells;
1287
    xelems = grid->yelems;
1288
    yelems = grid->xelems;
1289
  }
1290
  else { 
1291
    printf("No %d numbering scheme exists!\n",numbering); 
1292
    return(1);
1293
  }
1294

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

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

1401
	    /* Number the elements */
1402
	    if(cnew && level > 0) {
1403
	      if(level==1)  
1404
		cell[cnew].elem1st = elemno+1;	      
1405
	      if(level==2) 
1406
		cell[cnew].elemwidth  = (elemno+1) - cell[cnew].elem1st;
1407
	      elemno += xelems[i];	    
1408
	    }
1409
	  }
1410

1411
	  if(center && level < yelems[j]) {
1412
	    if (cnew) { 
1413
	      if(!cleft && sidenodes) 
1414
		knotno++;
1415
	      if(level==0 && center==1) 
1416
		cell[cnew].leftcenter = knotno;
1417
	      if(level==0 && center==2)
1418
		cell[cnew].left2center = knotno;		
1419
	      if(level==1 && center==1)
1420
		cell[cnew].levelwidthcenter = knotno - cell[cnew].leftcenter;
1421
	      if(centernodes && sidenodes)
1422
		knotno += degree * xelems[i];
1423
	      else 
1424
		knotno += xelems[i];
1425
	    } 
1426
	  }
1427
	
1428
	} /* x-cell loop */
1429
  
1430

1431
  maxwidth = 0;
1432
  for(j=1;j<= ycells ;j++) 
1433
    for(i=1;i<= xcells; i++) {
1434
      if(numbering == NUMBER_XY)
1435
	cnew = grid->numbered[j][i];
1436
      else if(numbering == NUMBER_YX)
1437
	cnew = grid->numbered[i][j];      
1438
      if(cnew) {
1439
	width = cell[cnew].left2nd - cell[cnew].left1st;
1440
	if(width > maxwidth) 
1441
	  maxwidth = width;
1442
	width = cell[cnew].levelwidth;
1443
	if(width > maxwidth) 
1444
	  maxwidth = width;
1445
	width = cell[cnew].leftlast-cell[cnew].left2nd
1446
	  -(yelems[j]-2)*cell[cnew].levelwidth;
1447
	if(width > maxwidth) 
1448
	  maxwidth = width;
1449
      }  
1450
    }
1451
  maxwidth += 2;
1452

1453
  grid->maxwidth = maxwidth;
1454
  grid->noknots = knotno;
1455
  grid->noelements = elemno;
1456

1457
  if(info) {
1458
    printf("Numbered %d knots in %d %d-node elements.\n",knotno,elemno,nonodes);
1459
    if(numbering == NUMBER_XY) 
1460
      printf("Numbering order was <x><y> and max levelwidth %d.\n",
1461
	     maxwidth);
1462
    else if(numbering == NUMBER_YX) 
1463
      printf("Numbering order was <y><x> and max levelwidth %d.\n",
1464
	     maxwidth);
1465
  }
1466

1467
  return(0);
1468
}
1469

1470

1471

1472
static int SetCellKnots1D(struct GridType *grid, struct CellType *cell,int info)
1473
{
1474
  int i;
1475
  int degree,nonodes;
1476
  int cnew,cleft;
1477
  int elemno,knotno;
1478
  int maxwidth;
1479
  int xcells,*xelems;
1480

1481
  nonodes = grid->nonodes;
1482
  knotno  = 0;
1483
  elemno  = 0;
1484

1485
  xcells = grid->xcells;
1486
  xelems = grid->xelems;
1487

1488
  switch (nonodes) {
1489
  case 2:
1490
    degree = 1;
1491
    break;
1492
  case 3:
1493
    degree = 2;
1494
    break;
1495
  case 4:
1496
    degree = 3;
1497
    break;
1498

1499
  default:
1500
    printf("SetCellKnots1D: No numbering scheme for %d-node elements.\n",
1501
	   grid->nonodes);
1502
    return(2);
1503
  }
1504

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

1533
  return(0);
1534
}
1535

1536

1537

1538
void CreateCells(struct GridType *grid,struct CellType **cell,int info)
1539
{
1540
  (*cell) = (struct CellType*)
1541
    malloc((size_t) (grid->nocells+1)*sizeof(struct CellType)); 
1542

1543
  SetCellData(grid,*cell,info);
1544

1545
  if(grid->dimension == 1) 
1546
    SetCellKnots1D(grid,*cell,info);
1547
  else
1548
    SetCellKnots(grid,*cell,info);
1549
}
1550

1551

1552
void DestroyCells(struct CellType **cell)
1553
{
1554
  free(cell);
1555
}
1556

1557

1558

1559
int GetKnotIndex(struct CellType *cell,int i,int j)
1560
/* Given the cell and knot indices gives the corresponding 
1561
   global knot index. The indices are expected to be in the 
1562
   range [0..n] and [0..m]. Requires only the structure CellType. 
1563
   */
1564
{
1565
  int ind=0,aid=0,maxj=0;
1566

1567
  if(cell->numbering == NUMBER_1D) {
1568
    ind = cell->left1st;
1569
    if(cell->nonodes == 2) 
1570
      ind += i;
1571
    else if(cell->nonodes == 3) 
1572
      ind += 2*i;
1573
    else if(cell->nonodes == 4) 
1574
      ind += 3*i;
1575
    return(ind);
1576
  }
1577

1578
  if(cell->numbering == NUMBER_XY) 
1579
    maxj = cell->yelem;
1580
  else if(cell->numbering == NUMBER_YX) {
1581
    aid = j; j = i; i = aid;
1582
    maxj = cell->xelem;
1583
  }
1584
  else {
1585
    maxj = 0;
1586
    bigerror("GetKnotIndex: Unknown numbering scheme!");
1587
  }
1588

1589
  if(j == 0) 
1590
    ind = cell->left1st;
1591
  else if (j == maxj)
1592
    ind = cell->leftlast;
1593
  else
1594
    ind = cell->left2nd + (j-1) * cell->levelwidth;
1595
  
1596
  if(cell->nonodes == 4) 
1597
    ind += i;
1598
  else if(cell->nonodes == 5)
1599
    ind += i;
1600
  else if(cell->nonodes == 8 || cell->nonodes == 9) 
1601
    ind += 2*i;
1602
  else if(cell->nonodes == 12 || cell->nonodes == 16) 
1603
    ind += 3*i;
1604

1605
  return(ind);
1606
}
1607

1608

1609
int GetKnotCoordinate(struct CellType *cell,int i,int j,Real *x,Real *y)
1610
/* Given the cell and element indices inside the cell gives the 
1611
   corresponding global knot numbers. The indices are expected 
1612
   to be in the range [0..n] and [0..m]. Requires only the 
1613
   structure CellType.   
1614
   */
1615
{
1616
  int ind;
1617

1618
  if(cell->xlinear == 1)
1619
    (*x) = cell->xcorner[BOTLEFT] + cell->dx1 * i;
1620
  else if(cell->xlinear == 0)
1621
    (*x) = cell->xcorner[BOTLEFT] + cell->dx1 * 
1622
    (1.- pow(cell->xratio,(Real)(i))) / (1.-cell->xratio);
1623
  else if(cell->xlinear == 2) {
1624
    if(i<=cell->xelem/2) {
1625
      (*x) = cell->xcorner[BOTLEFT] + cell->dx1 * 
1626
	(1.- pow(cell->xratio,(Real)(i))) / (1.-cell->xratio);
1627
    }
1628
    else {
1629
      (*x) = cell->xcorner[BOTRIGHT] - cell->dx1 * 
1630
	(1.- pow(cell->xratio,(Real)(cell->xelem-i))) / (1.-cell->xratio);
1631
    }
1632
  }
1633

1634
  if(cell->ylinear == 1)
1635
    (*y) = cell->ycorner[BOTLEFT] + cell->dy1 * j;
1636
  else if(cell->ylinear == 0) 
1637
    (*y) = cell->ycorner[BOTLEFT] + cell->dy1 * 
1638
           (1.- pow(cell->yratio,(Real)(j))) / (1.-cell->yratio);
1639
  else if(cell->ylinear == 2) {
1640
    if(j<=cell->yelem/2) {
1641
      (*y) = cell->ycorner[BOTLEFT] + cell->dy1 * 
1642
	(1.- pow(cell->yratio,(Real)(j))) / (1.-cell->yratio);
1643
    }
1644
    else {
1645
      (*y) = cell->ycorner[TOPLEFT] - cell->dy1 * 
1646
	(1.- pow(cell->yratio,(Real)(cell->yelem-j))) / (1.-cell->yratio);
1647
    }
1648
  }
1649

1650
  ind = GetKnotIndex(cell,i,j);
1651

1652
  return(ind);
1653
}
1654

1655

1656

1657
int GetElementIndices(struct CellType *cell,int i,int j,int *ind)
1658
/* For given element gives the coordinates and index for each knot.
1659
   The indices i and j are expected to range from [1..n] and [1..m]. 
1660
   requires only the structure CellType.
1661
   */
1662
{
1663
  int nonodes,numbering,elemind=0;
1664
  
1665
  nonodes = cell->nonodes;
1666
  numbering = cell->numbering;
1667

1668
  if(numbering != NUMBER_1D) {
1669
    ind[TOPRIGHT] = GetKnotIndex(cell,i,j);
1670
    ind[BOTRIGHT] = GetKnotIndex(cell,i,j-1);
1671
    ind[TOPLEFT]  = GetKnotIndex(cell,i-1,j);
1672
    ind[BOTLEFT]  = GetKnotIndex(cell,i-1,j-1);
1673
  }    
1674

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

1723
  else if(numbering == NUMBER_YX) { 
1724
    elemind = cell->elem1st+(j-1) + (i-1)*cell->elemwidth;
1725
    
1726
    if(nonodes == 4) return(elemind);
1727

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

1772
  else if(numbering == NUMBER_1D) {
1773
    elemind = cell->elem1st+(i-1);
1774
    ind[0] = GetKnotIndex(cell,i-1,1);
1775
    if(nonodes == 2) {
1776
      ind[1] = ind[0] + 1;
1777
    }
1778
    else if(nonodes == 3) {
1779
      ind[2] = ind[0] + 1;
1780
      ind[1] = ind[0] + 2;
1781
    }
1782
    else if(nonodes == 4) {
1783
      ind[2] = ind[0] + 1;
1784
      ind[3] = ind[0] + 2;
1785
      ind[1] = ind[0] + 3;
1786
    }
1787
  }  
1788
  return(elemind);
1789
}
1790

1791

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

1805
  return(elemind);
1806
}
1807

1808

1809

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

1822
  k = nonodes = cell->nonodes;
1823
  numbering = cell->numbering;
1824

1825
  if(numbering == NUMBER_1D) {
1826
    elemind = cell->elem1st+(i-1);    
1827
    ind[0] = GetKnotCoordinate(cell,i-1,j-1,&globalcoord[0],
1828
			       &globalcoord[nonodes]);
1829
    ind[1] = GetKnotCoordinate(cell,i,j-1,&globalcoord[1],
1830
			       &globalcoord[nonodes+1]);    
1831
    
1832
    if(nonodes == 3) {
1833
      globalcoord[2] = (globalcoord[0]+globalcoord[1])/2.0;
1834
      globalcoord[5] = (globalcoord[3]+globalcoord[4])/2.0;
1835
      ind[2] = ind[0] + 1;
1836
    }
1837
    else if(nonodes == 4) {
1838
      globalcoord[2] = (2.0*globalcoord[0]+globalcoord[1])/3.0;
1839
      globalcoord[6] = (2.0*globalcoord[4]+globalcoord[5])/3.0;
1840
      ind[2] = ind[0] + 1;
1841
      globalcoord[3] = (globalcoord[0]+2.0*globalcoord[1])/3.0;
1842
      globalcoord[7] = (globalcoord[4]+2.0*globalcoord[5])/3.0;
1843
      ind[3] = ind[0] + 2;
1844
    }
1845

1846
    return(elemind);
1847
  }
1848
  else if(numbering == NUMBER_XY) 
1849
    elemind = cell->elem1st+(i-1) + (j-1)*cell->elemwidth;
1850
  else if(numbering == NUMBER_YX)
1851
    elemind = cell->elem1st+(j-1) + (i-1)*cell->elemwidth;
1852

1853

1854
  ind[TOPRIGHT] = GetKnotCoordinate(cell,i,j,&globalcoord[TOPRIGHT],
1855
				    &globalcoord[TOPRIGHT+nonodes]);
1856
  ind[BOTRIGHT] = GetKnotCoordinate(cell,i,j-1,&globalcoord[BOTRIGHT],
1857
				    &globalcoord[BOTRIGHT+nonodes]);
1858
  ind[TOPLEFT] = GetKnotCoordinate(cell,i-1,j,&globalcoord[TOPLEFT],
1859
				    &globalcoord[TOPLEFT+nonodes]);
1860
  ind[BOTLEFT] = GetKnotCoordinate(cell,i-1,j-1,&globalcoord[BOTLEFT],
1861
				    &globalcoord[BOTLEFT+nonodes]);
1862
  if(nonodes == 4) return(elemind);
1863

1864
  GetElementIndices(cell,i,j,ind);
1865

1866
#if 0  
1867
  if(cell->xlinear) 
1868
    xrat = 1.0;
1869
  else
1870
    xrat = sqrt(cell->xratio);
1871

1872
  if(cell->ylinear) 
1873
    yrat = 1.0;
1874
  else 
1875
    yrat = sqrt(cell->yratio);
1876
#endif
1877
  xrat = yrat = 1.0;
1878

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

1927
#if 0
1928
  for(i=0;i<k;i++) 
1929
    printf("ind[%d]=%d  x=%.4lg  y=%.4lg\n",i,ind[i],
1930
	   globalcoord[i],globalcoord[i+k]);
1931
#endif
1932

1933
  return(elemind);
1934
}
1935

1936

1937
int GetSideInfo(struct CellType *cell,int cellno,int side,int element,
1938
		int *elemind)
1939
/* Given the side and element numbers returns the indices of 
1940
   the side element. When the end of the side is reached the function 
1941
   returns FALSE, otherwise TRUE.
1942
   */
1943
{
1944
  int more,sideno;                 
1945
  int ind[MAXNODESD2];
1946

1947
  more = TRUE;
1948

1949
  if(cell[cellno].numbering == NUMBER_1D) {
1950

1951
    switch(side) {
1952
    case 0:
1953
      more = FALSE;
1954
      elemind[0] = GetElementIndices(&(cell[cellno]),1,element,&(ind[0]));
1955
      if((sideno = cell[cellno].neighbour[LEFT]))
1956
	elemind[1] = GetElementIndices(&(cell[sideno]),cell[sideno].xelem,element,&(ind[0]));
1957
      else 
1958
	elemind[1] = 0;
1959
      break;
1960
     
1961
      
1962
    case 1:
1963
      more = FALSE; 
1964
      elemind[0] = GetElementIndices(&(cell)[cellno],cell[cellno].xelem,element,&(ind[0]));
1965
      if((sideno = cell[cellno].neighbour[RIGHT]))
1966
	elemind[1] = GetElementIndices(&(cell[sideno]),1,element,&(ind[0]));
1967
      else 
1968
	elemind[1] = 0;
1969
      break;
1970
    }
1971

1972
    return(more);
1973
  }
1974

1975

1976
  switch(side) {
1977

1978
  case LEFT:
1979
    if(element == cell[cellno].yelem) more = FALSE;
1980
    elemind[0] = GetElementIndices(&(cell[cellno]),1,element,&(ind[0]));
1981
    if((sideno = cell[cellno].neighbour[LEFT]))
1982
      elemind[1] = GetElementIndices(&(cell[sideno]),cell[sideno].xelem,element,&(ind[0]));
1983
    else 
1984
      elemind[1] = 0;
1985
    break;
1986

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

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

2005
  case UP:
2006
    if(element == cell[cellno].xelem) more = FALSE; 
2007
    elemind[0] = GetElementIndices(&(cell)[cellno],element,cell[cellno].yelem,&(ind[0]));
2008
    if((sideno = cell[cellno].neighbour[UP]))
2009
      elemind[1] = GetElementIndices(&(cell[sideno]),element,1,&(ind[0]));
2010
    else 
2011
      elemind[1] = 0;
2012
    break;
2013
    
2014
  default:
2015
    printf("Impossible option in GetSideInfo.\n");
2016
  }
2017

2018
  return(more);
2019
}
2020

2021

2022

2023
void SetElementDivisionExtruded(struct GridType *grid,int info)
2024
{
2025
  int i,nzmax = 0,sumzelems;
2026
  Real ratio,linearlimit;
2027
  Real dzmax = 0,dz = 0;
2028
  
2029
  linearlimit = 0.001;
2030

2031
  if(grid->autoratio) {
2032
    for(i=1;i<=grid->zcells;i++)  
2033
      grid->zelems[i] = grid->minzelems;
2034
  }
2035
  else {
2036
    for(i=1;i<=grid->zcells;i++)  
2037
      if(grid->zelems[i] < grid->minzelems) grid->zelems[i] = grid->minzelems;
2038
  }
2039

2040
  sumzelems = grid->zcells * grid->minzelems;
2041
  if(sumzelems > grid->totzelems) {
2042
#if 0
2043
    printf("SetElementDivision: %d is too few elements in z-direction (min. %d)\n",
2044
	   grid->totzelems,sumzelems);
2045
#endif
2046
    grid->totzelems = sumzelems;
2047
  }      
2048

2049
  /* Put the linearity flags. */
2050
  for(i=1; i<= grid->zcells ;i++) {
2051
    if (fabs(1.-grid->zexpand[i]) < linearlimit) 
2052
      grid->zlinear[i] = TRUE; 
2053
    else 
2054
      grid->zlinear[i] = FALSE;
2055
  }
2056

2057
  if(grid->autoratio) {
2058
    int active;
2059
    for(;;) {
2060
      dzmax = 0.0;
2061
      active = FALSE;
2062

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

2092
	if(dz > dzmax) {
2093
	  dzmax = dz;
2094
	  nzmax = i;
2095
	}
2096

2097
	if(grid->autoratio) {
2098
	  if(fabs(dz - grid->dx0) < fabs( dz*(1.0/(1.0-1.0/grid->zelems[i])) - grid->dx0) ) {
2099
	    grid->zelems[i] += 1;
2100
	    sumzelems++;
2101
	    active = TRUE;
2102
	  }
2103
	}
2104
      }
2105

2106
      if(grid->autoratio) {
2107
	if(!active) break;
2108
      }
2109
      else {
2110
	grid->zelems[nzmax] += 1;
2111
	sumzelems++;
2112
	if(sumzelems >= grid->totzelems) break;
2113
      }
2114
    }
2115
  }
2116

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

2147
  if(info) printf("Created %d extruded divisions.\n",
2148
		  grid->totzelems);
2149
  grid->dz0 = dzmax;
2150
}
2151

2152

2153

2154
void SetElementDivisionCylinder(struct GridType *grid,int info)
2155
{
2156
  int i,k;
2157
  Real ratio,eps;
2158
  Real dzmax;
2159
  
2160
  eps = 1.0e-8;
2161

2162
  grid->zcells = 2*grid->rotateblocks;
2163
  grid->z[0] = 0.0;
2164
  for(i=0;grid->x[i]<eps;i++); k=i;  
2165
  grid->rotateradius1 = grid->x[k];
2166

2167
  if(grid->rotateradius2 < 0.0) {
2168
    grid->rotatecartesian = TRUE;
2169
    grid->rotateradius2 = 1.0e6;
2170
  }
2171
  if(grid->rotateradius2 <= sqrt(2.0)*grid->rotateradius1) {
2172
    if(k+1 <= grid->xcells) 
2173
      grid->rotateradius2 = grid->x[k+1];
2174
    grid->rotateradius2 = MAX(sqrt(2.0)*grid->rotateradius1,
2175
			      grid->rotateradius2);
2176
  }
2177

2178
  if(!grid->xlinear[k]) 
2179
    printf("SetElementDivisionCylinder: The division must be linear!\n");
2180

2181
  for(i=1;i<=grid->zcells;i++) {
2182
    grid->z[i] = i*(2.0*FM_PI)/8.0; 
2183
    grid->zelems[i]  = grid->xelems[k];
2184
    grid->zlinear[i] = grid->xlinear[k];
2185
    grid->zratios[i] = grid->xratios[k];
2186
    grid->zexpand[i] = grid->xexpand[k];
2187
    grid->zmaterial[i] = 0;
2188
  }
2189

2190
  grid->totzelems =  grid->zcells * grid->xelems[k];
2191

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

2219
    if(dzmax < grid->dz[i]) dzmax = grid->dz[i];
2220
  }
2221

2222
  if(info) printf("Created %d divisions in %d cells for rotation [%.2lg  %.2lg].\n",
2223
		  grid->totzelems,grid->zcells,
2224
		  grid->rotateradius1,grid->rotateradius2);
2225
  grid->dz0 = dzmax;
2226
}
2227

2228

2229

2230
static int GetCommand(char *line1,char *line2,FILE *io) 
2231
/* Line1 for commands and line2 for arguments. */
2232
{
2233
  int i,j,isend;
2234
  char line0[MAXLINESIZE],*charend,*matcpntr0,*matcpntr;
2235
  int gotlinefeed;
2236

2237
 newline:
2238

2239
  for(i=0;i<MAXLINESIZE;i++)
2240
    line2[i] = line1[i] = line0[i] = 0x00;
2241

2242
  charend = fgets(line0,MAXLINESIZE,io);
2243
  isend = (charend == NULL);
2244

2245
  if(isend) return(1);
2246

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

2249
  if(line0[0] == '#' || line0[0] == '%' || line0[0] == '!' || line0[0] == '\n') goto newline;
2250
  if(!matcactive && line0[0] == '*') goto newline;
2251

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

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

2294
  /* After these commands there will be no nextline even though there is no equality sign */
2295
  if(strstr(line1,"END")) return(0);
2296
  if(strstr(line1,"NEW MESH")) return(0);
2297

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

2320
    gotlinefeed = FALSE;
2321
    for(i=0;i<MAXLINESIZE;i++) {
2322
      if(line2[i] == '\n' || line2[i] == '\0' ) {
2323
	gotlinefeed = TRUE;
2324
	break;
2325
      }
2326
    }
2327
    if(!gotlinefeed) {
2328
      printf("There is a risk that somethings was missed in line:\n");
2329
      printf("%s\n",line2);
2330
      smallerror("Check your output line length!\n");
2331
    }
2332

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

2344
#if USE_MATC
2345
    if(matcactive) {
2346
      matcpntr0 = strchr(line2,'$');
2347
      if(matcpntr0) {
2348
	matcpntr = mtc_domath(&matcpntr0[1]);
2349
	if(matcpntr) {
2350
	  strcpy(matcpntr0, matcpntr);
2351
	  if(0) printf("C: %s\n%s\n",matcpntr0,matcpntr);
2352
	}
2353
      }
2354
    }
2355
#endif 
2356
  }
2357

2358
  if(iodebug) {
2359
    printf("command: %s\n",line1);
2360
    printf("params:  %s\n",line2);
2361
  }
2362
  
2363
  return(0);
2364
}
2365

2366

2367

2368
int SaveElmergrid(struct GridType *grid,int nogrids,char *prefix,int info)
2369
{
2370
  int sameline,maxsameline;
2371
  int i,j,dim;
2372
  FILE *out;
2373
  char filename[MAXFILESIZE];
2374

2375
  AddExtension(prefix,filename,"grd");
2376
  out = fopen(filename,"w");
2377
  dim = grid->dimension;
2378
  if(grid->coordsystem == COORD_CART1) dim = 1;
2379

2380
  j = 0;
2381
  sameline = TRUE;
2382
  maxsameline = 6;
2383
  if(grid->xcells > maxsameline) sameline = FALSE;
2384
  if(dim >= 2 && grid->ycells > maxsameline) sameline = FALSE;
2385
  if(dim >= 3 && grid->zcells > maxsameline) sameline = FALSE;
2386
  
2387
  fprintf(out,"#####  ElmerGrid input file for structured grid generation  ######\n");
2388
  fprintf(out,"Version = 210903\n");
2389

2390
  fprintf(out,"Coordinate System = ");
2391
  if(grid->coordsystem == COORD_AXIS)
2392
    fprintf(out,"2D Axisymmetric\n");
2393
  else if(grid->coordsystem == COORD_POLAR)
2394
    fprintf(out,"2D Polar\n");
2395
  else 
2396
    fprintf(out,"Cartesian %dD\n",dim);
2397
 
2398
  fprintf(out,"Subcell Divisions in %dD = ",dim);
2399
  if(dim >= 1) fprintf(out,"%d ",grid->xcells);
2400
  if(dim >= 2) fprintf(out,"%d ",grid->ycells);
2401
  if(dim >= 3) fprintf(out,"%d ",grid->zcells);
2402
  fprintf(out,"\n");
2403

2404
  fprintf(out,"Subcell Limits 1 %s",sameline ? "= ":"\n  ");
2405
  for(i=0;i <= grid->xcells;i++) 
2406
    fprintf(out,"%.5lg ",grid->x[i]); 
2407
  fprintf(out,"\n");
2408
    
2409
  if(dim >= 2) {
2410
    fprintf(out,"Subcell Limits 2 %s",sameline ? "= ":"\n  ");
2411
    for(i=0;i <= grid->ycells;i++) 
2412
      fprintf(out,"%.5lg ",grid->y[i]); 
2413
    fprintf(out,"\n");
2414
  }
2415
  
2416
  if(dim >= 3) {
2417
    fprintf(out,"Subcell Limits 3 %s",sameline ? "= ":"\n  ");
2418
    for(i=0;i <= grid->zcells;i++) 
2419
      fprintf(out,"%.5lg ",grid->z[i]); 
2420
    fprintf(out,"\n");
2421
  }  
2422

2423
  fprintf(out,"Material Structure in %dD\n",dim==1 ? 1:2);  
2424
  for(j=grid->ycells;j>=1;j--) {
2425
    fprintf(out,"  ");
2426
    for(i=1;i<=grid->xcells;i++) 
2427
      fprintf(out,"%-5d",grid->structure[j][i]);
2428
    fprintf(out,"\n");
2429
  }
2430
  fprintf(out,"End\n");
2431

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

2447
  j = 0;
2448
  if(grid[j].rotate) {
2449
    fprintf(out,"Revolve Blocks = %d\n",grid[j].rotateblocks);
2450
    fprintf(out,"Revolve Radius = %-8.3lg\n",grid[j].rotateradius2);
2451
    if(fabs(grid[j].rotateimprove-1.0) > 1.0e-10)
2452
      fprintf(out,"Revolve Improve = %-8.3lg\n",grid[j].rotateimprove);
2453
    
2454
  }
2455
  if(grid[j].rotatecurve) {
2456
    fprintf(out,"Revolve Curve Direct = %-8.3lg\n",grid[j].curvezet);
2457
    fprintf(out,"Revolve Curve Radius = %-8.3lg\n",grid[j].curverad);
2458
    fprintf(out,"Revolve Curve Angle = %-8.3lg\n",grid[j].curveangle);
2459
  }
2460

2461
  if(grid[j].coordsystem == COORD_POLAR) {
2462
    fprintf(out,"Polar Radius = %.3lg\n",grid[j].polarradius);
2463
  } 
2464

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

2482
    if(grid[j].noboundaries > 0) {
2483
      fprintf(out,"Boundary Definitions\n");
2484
      fprintf(out,"# %-8s %-8s %-8s\n","type","out","int"); 
2485
      for(i=0;i<grid[j].noboundaries;i++)
2486
	fprintf(out,"  %-8d %-8d %-8d %-8d\n",
2487
		grid[j].boundtype[i],grid[j].boundext[i],
2488
		grid[j].boundint[i], grid[j].boundsolid[i]);
2489
      fprintf(out,"End\n");
2490
    }
2491

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

2507
    if(dim == 2)
2508
      fprintf(out,"Coordinate Ratios = %-8.3lg\n",grid[j].xyratio);
2509
    if(dim == 3)
2510
      fprintf(out,"Coordinate Ratios = %-8.3lg %-8.3lg\n",
2511
	      grid[j].xyratio,grid[j].xzratio);
2512
 
2513
    fprintf(out,"Minimum Element Divisions = %d",grid[j].minxelems);
2514
    if(dim >= 2) fprintf(out," %d",grid[j].minyelems);
2515
    if(dim >= 3) fprintf(out," %d",grid[j].minzelems);
2516
    fprintf(out,"\n");
2517

2518
    fprintf(out,"Element Ratios 1 %s",sameline ? "= ":"\n  ");
2519
    for(i=1;i<=grid[j].xcells;i++) 
2520
      fprintf(out,"%.3lg ",grid[j].xexpand[i]); 
2521
    fprintf(out,"\n");
2522
    if(dim >= 2) {
2523
      fprintf(out,"Element Ratios 2 %s",sameline ? "= ":"\n  ");
2524
      for(i=1;i<=grid[j].ycells;i++) 
2525
	fprintf(out,"%.3lg ",grid[j].yexpand[i]); 
2526
      fprintf(out,"\n");
2527
    }
2528
    if(dim >= 3) {
2529
      fprintf(out,"Element Ratios 3 %s",sameline ? "= ":"\n  ");
2530
      for(i=1;i<=grid[j].zcells;i++) 
2531
	fprintf(out,"%.3lg ",grid[j].zexpand[i]); 
2532
      fprintf(out,"\n");
2533
    }
2534

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

2573
  }
2574

2575
  if(info) printf("The Elmergrid input was saved to file %s.\n",filename);
2576
  fclose(out);
2577

2578
  return(0);
2579
}
2580

2581

2582

2583

2584

2585
static int LoadElmergridOld(struct GridType **grid,int *nogrids,char *prefix,int info) 
2586
{
2587
  char filename[MAXFILESIZE];
2588
  FILE *in;
2589
  int i,j,k,l,error=0;
2590
  Real scaling;
2591
  char *cp;
2592
  int mode,noknots,noelements,dim,axisymmetric;
2593
  int elemcode,maxnodes,totelems,nogrids0;
2594
  int minmat,maxmat;
2595
  char line[MAXLINESIZE];
2596

2597

2598
  AddExtension(prefix,filename,"grd");
2599
  if ((in = fopen(filename,"r")) == NULL) {
2600
    printf("LoadElmergrid: opening of the file '%s' wasn't successful !\n",filename);
2601
    return(1);
2602
  }
2603

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

2606
  InitGrid(grid[*nogrids]);
2607
  k = *nogrids;
2608
  nogrids0 = *nogrids;
2609

2610
  mode = 0;
2611
  noknots = 0;
2612
  noelements = 0;
2613
  dim = 0;
2614
  axisymmetric = FALSE;
2615
  elemcode = 0;
2616
  maxnodes = 4;
2617
  totelems = 0;
2618
  scaling = 1.0;
2619

2620

2621

2622
  Getline(line,in);
2623
  for(;;) {
2624
    if(Getline(line,in)) goto end;
2625
    if(line[0]=='\0') goto end;
2626

2627
    if(strstr(line,"END")) goto end;
2628
    if(strstr(line,"RESULTS")) goto end;
2629

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

2661

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

2691
      Getline(line,in);
2692

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

2719

2720
      for(j=1;j<=grid[k]->dimension;j++) {
2721
	Getline(line,in);
2722
	cp=line;
2723

2724
	if(j==1) for(i=0;i<=grid[k]->xcells;i++) grid[k]->x[i] = next_real(&cp);
2725
	if(j==2) for(i=0;i<=grid[k]->ycells;i++) grid[k]->y[i] = next_real(&cp);
2726
	if(j==3) for(i=0;i<=grid[k]->zcells;i++) grid[k]->z[i] = next_real(&cp);
2727
      }
2728

2729
      printf("Loading material structure\n");
2730

2731
      for(j=grid[k]->ycells;j>=1;j--) {
2732
	
2733
	Getline(line,in);
2734
	cp=line;
2735
	
2736
	for(i=1;i<=grid[k]->xcells;i++) 
2737
	  grid[k]->structure[j][i] = next_int(&cp);
2738
      }
2739

2740
      minmat = maxmat = grid[k]->structure[1][1];
2741
      for(j=grid[k]->ycells;j>=1;j--) 
2742
	for(i=1;i<=grid[k]->xcells;i++) {
2743
	  if(minmat > grid[k]->structure[j][i])
2744
	    minmat = grid[k]->structure[j][i];
2745
	  if(maxmat < grid[k]->structure[j][i])
2746
	    maxmat = grid[k]->structure[j][i];
2747
	}      
2748
      if(minmat < 0) 
2749
	printf("LoadElmergrid: please use positive material indices.\n");
2750

2751
      mode = 0;
2752
      break;
2753

2754
    case 3:
2755
    case 31:
2756
    case 32:
2757

2758
      /* I don't know how to set this, luckily this piece of code should be obsolete */
2759
      l = 1;
2760
      for(i=grid[k]->mappings;i<grid[k]->mappings+l;i++) {
2761
	Getline(line,in);
2762
	cp=line; 
2763

2764
	grid[k]->mappingtype[i] = next_int(&cp);
2765
	if(mode == 32) grid[k]->mappingtype[i] += 50*SGN(grid[k]->mappingtype[i]);
2766

2767
	grid[k]->mappingline[i] = next_int(&cp);
2768
	grid[k]->mappinglimits[2*i] = next_real(&cp);
2769
	grid[k]->mappinglimits[2*i+1] = next_real(&cp);
2770
	grid[k]->mappingpoints[i] = next_int(&cp);
2771
	grid[k]->mappingparams[i] = Rvector(0,grid[k]->mappingpoints[i]);
2772
	for(j=0;j<grid[k]->mappingpoints[i];j++) 
2773
	  grid[k]->mappingparams[i][j] = next_real(&cp);
2774
      }
2775
      
2776
      printf("Loaded %d geometry mappings\n",l);
2777
      grid[k]->mappings += l;
2778

2779
      mode = 0;
2780
      break;
2781
      
2782
    case 4: /* NUMBERING */
2783
      if(strstr(line,"HORIZ")) grid[k]->numbering = NUMBER_XY;
2784
      if(strstr(line,"VERTI")) grid[k]->numbering = NUMBER_YX;
2785
      mode = 0;
2786
      break;
2787

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

2801
    case 6: /* ELEMENTS */
2802
      sscanf(line,"%d",&(*grid)[k].wantedelems);
2803
      mode = 0;
2804
      break;
2805

2806
    case 7: /* NODES */
2807
      sscanf(line,"%d",&(*grid)[k].nonodes);      
2808
      
2809
      (*grid)[k].elemmidpoints = FALSE;
2810
      if((*grid)[k].nonodes == 4) 
2811
	(*grid)[k].elemorder = 1;
2812
      if((*grid)[k].nonodes == 8) 
2813
	(*grid)[k].elemorder = 2;
2814
      if((*grid)[k].nonodes == 16) 
2815
	(*grid)[k].elemorder = 3;
2816

2817
      if((*grid)[k].nonodes == 9) { 
2818
	(*grid)[k].elemorder = 2;
2819
	(*grid)[k].elemmidpoints = TRUE;
2820
      }
2821
      if((*grid)[k].nonodes == 12) { 
2822
	(*grid)[k].elemorder = 3;
2823
	(*grid)[k].elemmidpoints = TRUE;
2824
      }
2825

2826

2827
      mode = 0;
2828
      break;
2829

2830
    case 8: /* TRIANGLES */
2831
      (*grid)[k].triangles = TRUE;
2832
      mode = 0;
2833
      break;
2834

2835
    case 17: /* SQUARES */
2836
      (*grid)[k].triangles = FALSE;
2837
      mode = 0;
2838
      break;
2839

2840
    case 16: /* ELEMENTTYPE and ELEMENTCODE */
2841
      sscanf(line,"%d",&elemcode);
2842
      if(elemcode/100 == 2) {
2843
	(*grid)[k].triangles = FALSE;      
2844
	(*grid)[k].nonodes = elemcode%100;
2845
      }
2846
      else if(elemcode/100 == 4) {
2847
	(*grid)[k].triangles = FALSE;      
2848
	(*grid)[k].nonodes = elemcode%100;
2849
      }
2850
      else if(elemcode/100 == 3) {  
2851
	(*grid)[k].triangles = TRUE;      
2852
	if(elemcode%100 == 3)       (*grid)[k].nonodes = 4;
2853
	else if(elemcode%100 == 6)  (*grid)[k].nonodes = 9;
2854
	else if(elemcode%100 == 10) (*grid)[k].nonodes = 16;	
2855
      }
2856

2857
      (*grid)[k].elemmidpoints = FALSE;
2858
      if((*grid)[k].nonodes == 4) 
2859
	(*grid)[k].elemorder = 1;
2860
      if((*grid)[k].nonodes == 8) 
2861
	(*grid)[k].elemorder = 2;
2862
      if((*grid)[k].nonodes == 16) 
2863
	(*grid)[k].elemorder = 3;
2864

2865
      if((*grid)[k].nonodes == 9) { 
2866
	(*grid)[k].elemorder = 2;
2867
	(*grid)[k].elemmidpoints = TRUE;
2868
      }
2869
      if((*grid)[k].nonodes == 12) { 
2870
	(*grid)[k].elemorder = 3;
2871
	(*grid)[k].elemmidpoints = TRUE;
2872
      }
2873

2874
      mode = 0;
2875
      break;
2876

2877
    case 10: /* COORDINATE RATIO */
2878
      if((*grid)[k].dimension == 2) 
2879
	sscanf(line,"%le",&(*grid)[k].xyratio);
2880
      if((*grid)[k].dimension == 3) 
2881
	sscanf(line,"%le %le",&(*grid)[k].xyratio,&(*grid)[k].xzratio);      
2882
      mode = 0;
2883
      break;
2884

2885
    case 11: /* MATERIALS */
2886
      sscanf(line,"%d %d",&(*grid)[k].firstmaterial,&(*grid)[k].lastmaterial);      
2887
      mode = 0;
2888
      break;
2889

2890
    case 12: /* LAYERES */
2891
      for(i=1;i<=(*grid)[k].zcells;i++) {
2892
	Getline(line,in);
2893
	sscanf(line,"%d %d %d\n",
2894
		&(*grid)[k].zfirstmaterial[i],&(*grid)[k].zlastmaterial[i],&(*grid)[k].zmaterial[i]); 
2895
      }
2896
      mode = 0;
2897
      break;
2898

2899
    case 13: /* ELEMENT RATIOS */
2900
      printf("Loading element ratios\n");
2901

2902
      for (j=1;j<=(*grid)[k].dimension;j++) {
2903
	Getline(line,in);
2904
	cp = line;
2905

2906
	if(j==1) for(i=1;i<=(*grid)[k].xcells;i++) (*grid)[k].xexpand[i] = next_real(&cp);
2907
	if(j==2) for(i=1;i<=(*grid)[k].ycells;i++) (*grid)[k].yexpand[i] = next_real(&cp);
2908
	if(j==3) for(i=1;i<=(*grid)[k].zcells;i++) (*grid)[k].zexpand[i] = next_real(&cp);
2909
      }
2910
      mode = 0;
2911
      break;
2912

2913
    case 29: /* ELEMENT NUMBER */
2914
      printf("Loading element numbers\n");
2915

2916
      for (j=1;j<=(*grid)[k].dimension;j++) {
2917
	Getline(line,in);
2918
	cp = line;
2919
	if(j==1) for(i=1;i<=(*grid)[k].xcells;i++) (*grid)[k].xelems[i] = next_int(&cp);
2920
	if(j==2) for(i=1;i<=(*grid)[k].ycells;i++) (*grid)[k].yelems[i] = next_int(&cp);
2921
	if(j==3) for(i=1;i<=(*grid)[k].zcells;i++) (*grid)[k].zelems[i] = next_int(&cp);
2922
      }
2923
      (*grid)[k].autoratio = 0;
2924
      mode = 0;
2925
      break;
2926

2927
    case 27: /* ELEMENT MINIMUM */
2928
      printf("Loading minimum number of elements\n");
2929
      if((*grid)[k].dimension == 1) 
2930
	sscanf(line,"%d",&(*grid)[k].minxelems);
2931
      if((*grid)[k].dimension == 2) 
2932
	sscanf(line,"%d %d",&(*grid)[k].minxelems,&(*grid)[k].minyelems);
2933
      if((*grid)[k].dimension == 3) 
2934
	sscanf(line,"%d %d %d",&(*grid)[k].minxelems,&(*grid)[k].minyelems,&(*grid)[k].minzelems);
2935
      mode = 0;
2936
      break;
2937

2938
    case 14: /* ELEMENT DENSITIES */
2939
      printf("Loading element densities\n");
2940
      for (j=1;j<=(*grid)[k].dimension;j++) {
2941
	Getline(line,in);
2942
	cp = line;
2943

2944
	if(j==1) for(i=1;i<=(*grid)[k].xcells;i++) (*grid)[k].xdens[i] = next_real(&cp);
2945
	if(j==2) for(i=1;i<=(*grid)[k].ycells;i++) (*grid)[k].ydens[i] = next_real(&cp);
2946
	if(j==3) for(i=1;i<=(*grid)[k].zcells;i++) (*grid)[k].zdens[i] = next_real(&cp);
2947
      }
2948
      mode = 0;
2949
      break;
2950

2951
    case 15: /* BOUNDARY CONDITIONS */
2952
      sscanf(line,"%d",&(*grid)[k].noboundaries);
2953
      printf("Loading %d boundary conditions\n",(*grid)[k].noboundaries);
2954

2955
      for(i=0;i<(*grid)[k].noboundaries;i++) {
2956
	Getline(line,in);
2957
	sscanf(line,"%d %d %d %d",
2958
	       &(*grid)[k].boundtype[i],&(*grid)[k].boundext[i],
2959
	       &(*grid)[k].boundint[i],&(*grid)[k].boundsolid[i]);
2960
      }  
2961
      mode = 0;
2962
      break;
2963

2964
    case 20: /* ROTATE */
2965
      (*grid)[k].rotate = TRUE;
2966
      mode = 0;
2967
      break;
2968

2969
    case 21: /* ROTRAD */
2970
      sscanf(line,"%le",&(*grid)[k].rotateradius2);
2971
      mode = 0;
2972
      break;
2973

2974
    case 22: /* ROTBLOCK */
2975
      sscanf(line,"%d",&(*grid)[k].rotateblocks);
2976
      if(0) printf("Reading blocks %d\n",(*grid)[k].rotateblocks);
2977
      mode = 0;
2978
      break;
2979

2980
    case 24: /* ROTIMP */
2981
      sscanf(line,"%le",&(*grid)[k].rotateimprove);
2982
      mode = 0;
2983
      break;
2984

2985
    case 30: /* POLAR RADIUS */
2986
      sscanf(line,"%le",&(*grid)[k].polarradius);
2987
      mode = 0;
2988
      break;
2989

2990
    case 25: /* ROTCURVE */
2991
      (*grid)[k].rotatecurve = TRUE;
2992
      sscanf(line,"%le%le%le",&(*grid)[k].curvezet,
2993
	     &(*grid)[k].curverad,&(*grid)[k].curveangle);
2994
      mode = 0;
2995
      break;
2996

2997
    case 26: /* REDUCE ELEMENT */
2998
      sscanf(line,"%d%d",&(*grid)[k].reduceordermatmin,
2999
	     &(*grid)[k].reduceordermatmax);
3000
      mode = 0;
3001
      break;
3002

3003
    case 28: /* LAYERED BO */
3004
      sscanf(line,"%d",&(*grid)[k].layeredbc);
3005
      mode = 0;
3006
      break;
3007

3008
    case 23: /* SCALING */
3009
      sscanf(line,"%le",&scaling);
3010
      for(i=0;i<=(*grid)[k].xcells;i++) (*grid)[k].x[i] *= scaling;
3011
      if((*grid)[k].dimension > 1) 
3012
	for(i=0;i<=(*grid)[k].ycells;i++) (*grid)[k].y[i] *= scaling;
3013
      if((*grid)[k].dimension == 3) 
3014
	for(i=0;i<=(*grid)[k].ycells;i++) (*grid)[k].z[i] *= scaling;
3015

3016
      (*grid)[k].rotateradius2 *= scaling;
3017
      (*grid)[k].curverad *= scaling;
3018
      (*grid)[k].curvezet *= scaling;
3019
      mode = 0;
3020
      break;
3021

3022
    default:
3023
      if(0) printf("Unknown case: %s",line);
3024
    }
3025

3026
  }
3027

3028
end:
3029

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

3032
  for(k=nogrids0;k < (*nogrids) && k<MAXCASES;k++) {
3033
    SetElementDivision(&(*grid)[k],1.0,info);
3034
  }
3035

3036

3037
  fclose(in);
3038
  return(error);
3039
}
3040

3041

3042

3043
int LoadElmergrid(struct GridType **grid,int *nogrids,char *prefix,int info) 
3044
{
3045
  char filename[MAXFILESIZE];
3046
  char command[MAXLINESIZE],params[MAXLINESIZE];
3047
  FILE *in;
3048
  int i,j,k;
3049
  char *cp;
3050
  int noknots,noelements,dim,axisymmetric;
3051
  int elemcode,maxnodes,totelems,nogrids0,minmat,maxmat;
3052
  long code;
3053
  Real raid;
3054

3055
  AddExtension(prefix,filename,"grd");
3056
  if ((in = fopen(filename,"r")) == NULL) {
3057
    printf("LoadElmergrid: opening of the file '%s' wasn't successful !\n",filename);
3058
    return(1);
3059
  }
3060

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

3063
  InitGrid(grid[*nogrids]);
3064
  k = *nogrids;
3065
  nogrids0 = *nogrids;
3066

3067
  iodebug = FALSE;
3068

3069
  noknots = 0;
3070
  noelements = 0;
3071
  dim = 0;
3072
  axisymmetric = FALSE;
3073
  elemcode = 0;
3074
  maxnodes = 4;
3075
  totelems = 0;
3076

3077
  matcactive = FALSE;
3078

3079
  for(;;) {
3080
    if(GetCommand(command,params,in)) {
3081
      if(info) printf("Reached the end of command file\n");
3082
      goto end;
3083
    }    
3084

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

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

3164
      if(grid[k]->xcells >= MAXCELLS || grid[k]->ycells >= MAXCELLS || grid[k]->zcells >= MAXCELLS) {
3165
	printf("LoadElmergrid: Too many subcells [%d %d %d] vs. %d:\n",
3166
	       grid[k]->xcells,grid[k]->ycells,grid[k]->zcells,MAXCELLS);
3167
     }
3168

3169
      /* Initialize the default structure with ones */
3170
      for(j=grid[k]->ycells;j>=1;j--) 
3171
	for(i=1;i<=grid[k]->xcells;i++) 
3172
	  grid[k]->structure[j][i] = 1;
3173
    }
3174
    
3175
    else if(strstr(command,"MINIMUM ELEMENT DIVISION")) {
3176
      if(info) printf("Loading minimum number of elements\n");
3177

3178
      if((*grid)[k].dimension == 1) 
3179
	sscanf(params,"%d",&(*grid)[k].minxelems); 
3180

3181
      if((*grid)[k].dimension == 2)  
3182
	sscanf(params,"%d %d",&(*grid)[k].minxelems,&(*grid)[k].minyelems);
3183

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

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

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

3294
    else if(strstr(command,"MATERIAL STRUCTURE")) {
3295
      if(info) printf("Loading material structure\n");
3296

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

3323
      grid[k]->maxmaterial = maxmat;
3324
    }
3325
    else if(strstr(command,"MATERIALS INTERVAL")) { 
3326
      sscanf(params,"%d %d",&(*grid)[k].firstmaterial,&(*grid)[k].lastmaterial);      
3327
    }
3328
     
3329
    else if(strstr(command,"REVOLVE")) {
3330
      if(0) printf("revolve: %s %s\n",command,params);
3331

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

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

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

3497
      for(i=0;i<MAXLINESIZE;i++) params[i] = toupper(params[i]);
3498
      for(i=(*grid)[k].mappings;i<MAXMAPPINGS;i++) {
3499
	if(i>(*grid)[k].mappings) Getline(params,in);
3500

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

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

3535
    else {
3536
      if(0) printf("Unknown command: %s",command);
3537
    }
3538
  }
3539

3540
end:
3541

3542
  if(info) printf("Found %d divisions for grid\n",*nogrids);
3543
  
3544
  for(k=nogrids0;k < (*nogrids) && k<MAXCASES;k++) {
3545

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

3574
  
3575
  /* for(k=nogrids0;k < (*nogrids) && k<MAXCASES;k++) {
3576
    SetElementDivision(&(*grid)[k],relh,info);
3577
    } */
3578

3579
  fclose(in);
3580
  return(0);
3581
}
3582

3583

3584

3585
void InitParameters(struct ElmergridType *eg)
3586
{
3587
  int i;
3588

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

3681

3682

3683

3684
int InlineParameters(struct ElmergridType *eg,int argc,char *argv[],int first,int info)
3685
{
3686
  int arg,i,dim;
3687
  char command[MAXLINESIZE];
3688
  
3689
  dim = eg->dim;
3690

3691
  printf("Elmergrid reading in-line arguments\n");
3692

3693
  /* Type of input file */
3694
  if(first > 3) {
3695
    for(i=0;i<MAXLINESIZE;i++) command[i] = ' ';
3696

3697
    strcpy(command,argv[1]);
3698
    for(i=0;i<MAXLINESIZE;i++) command[i] = toupper(command[i]);
3699
    for(i=0;i<=MAXINMETHODS;i++) {
3700
      if(strstr(command,InMethods[i])) {
3701
	eg->inmethod = i;
3702
	break;
3703
      }
3704
    }
3705
    if(i>MAXINMETHODS) eg->inmethod = atoi(argv[1]);
3706
    
3707

3708
    /* Type of output file (fewer options) */
3709
    strcpy(command,argv[2]);
3710
    for(i=0;i<MAXLINESIZE;i++) command[i] = toupper(command[i]);
3711
    for(i=1;i<=MAXOUTMETHODS;i++) {
3712
      if(strstr(command,OutMethods[i])) {
3713
	eg->outmethod = i;
3714
	break;
3715
      }
3716
    }
3717
    if(i>MAXOUTMETHODS) eg->outmethod = atoi(argv[2]);
3718
 
3719
    /* Name of output file */
3720
    strcpy(eg->filesin[0],argv[3]);
3721
    strcpy(eg->filesout[0],eg->filesin[0]);
3722
    strcpy(eg->infofile,eg->filesin[0]);
3723
  }
3724

3725

3726
  /* The optional inline parameters */
3727

3728
  for(arg=first;arg <argc; arg++) {
3729

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

4373
  {
4374
    int badpoint;
4375
    char *ptr1,*ptr2;
4376
    ptr1 = strrchr(eg->filesout[0], '.');
4377
    if(ptr1) {
4378
      badpoint=FALSE;
4379
      ptr2 = strrchr(eg->filesout[0], '/');
4380
      if(ptr2 && ptr2 > ptr1) badpoint = TRUE;
4381
      if(!badpoint) *ptr1 = '\0';
4382
    }
4383
  }
4384

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

4387
  return(0);
4388
}
4389

4390

4391

4392

4393
int LoadCommands(char *prefix,struct ElmergridType *eg,
4394
		 struct GridType *grid, int mode,int info) 
4395
{
4396
  char filename[MAXFILESIZE];
4397
  char command[MAXLINESIZE],params[MAXLINESIZE],*cp;
4398

4399
  FILE *in=NULL;
4400
  int i,j,iostat;
4401

4402
  iodebug = FALSE;
4403

4404
  if( mode == 0) {  
4405
    if ((in = fopen("ELMERGRID_STARTINFO","r"))) {
4406
      iostat = fscanf(in,"%s",filename);
4407
      fclose(in);
4408
      printf("Using the file %s defined in ELMERGRID_STARTINFO\n",filename);
4409
      if ((in = fopen(filename,"r")) == NULL) {
4410
	printf("LoadCommands: opening of the file '%s' wasn't successful !\n",filename);
4411
	return(1);
4412
      }    
4413
      else printf("Loading ElmerGrid commands from file '%s'.\n",filename);    
4414
    }    
4415
    else 
4416
      return(2);
4417
  }
4418
  else if(mode == 1) { 
4419
    AddExtension(prefix,filename,"eg");
4420
    if ((in = fopen(filename,"r")) == NULL) {
4421
      printf("LoadCommands: opening of the file '%s' wasn't successful !\n",filename);
4422
      return(3);
4423
    }    
4424
    if(info) printf("Loading ElmerGrid commands from file '%s'.\n",filename);    
4425
  }
4426
  else if(mode == 2) {
4427
    AddExtension(prefix,filename,"grd");
4428
    if ((in = fopen(filename,"r")) == NULL) {
4429
      printf("LoadCommands: opening of the file '%s' wasn't successful !\n",filename);
4430
      return(4);
4431
    }    
4432
    if(info) printf("\nLoading ElmerGrid commands from file '%s'.\n",filename);
4433
  }
4434

4435

4436

4437
  for(;;) {
4438

4439
    if(GetCommand(command,params,in)) {
4440
      printf("Reached the end of command file\n");
4441
      goto end;
4442
    }    
4443

4444
    /* If the mode is the command file mode read also the file information from the command file. */
4445

4446
    if(mode <= 1) {
4447
      if(strstr(command,"INPUT FILE")) {
4448
	sscanf(params,"%s",eg->filesin[0]);
4449
      }
4450

4451
      else if(strstr(command,"OUTPUT FILE")) {
4452
	sscanf(params,"%s",eg->filesout[0]);
4453
	eg->filerenamed = TRUE;
4454
      }
4455

4456
      else if(strstr(command,"INPUT MODE")) {
4457
	for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4458
	
4459
	for(i=0;i<=MAXINMETHODS;i++) {
4460
	  if(strstr(params,InMethods[i])) {
4461
	    eg->inmethod = i;
4462
	    break;
4463
	  }
4464
	}
4465
	if(i>MAXINMETHODS) sscanf(params,"%d",&eg->inmethod);
4466
      }
4467

4468
      else if(strstr(command,"OUTPUT MODE")) {
4469
	for(j=0;j<MAXLINESIZE;j++) params[j] = toupper(params[j]);
4470
	
4471
	/* Type of output file (fewer options) */
4472
	for(i=1;i<=MAXOUTMETHODS;i++) {
4473
	  if(strstr(params,OutMethods[i])) {
4474
	    eg->outmethod = i;
4475
	    break;
4476
	  }
4477
	}
4478
	if(i>MAXOUTMETHODS) sscanf(params,"%d",&eg->outmethod);	
4479
      }
4480
    }    
4481
    /* End of command file specific part */
4482

4483

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

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

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

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

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

4875
    }
4876
  }
4877

4878
end:
4879
  printf("Read commands from a file\n");
4880
  fclose(in);
4881
  return(0);
4882
}
4883

4884

4885

4886

4887

4888
static int FindParentSide(struct FemType *data,struct BoundaryType *bound,
4889
			  int sideelem,int sideelemtype,int *sideind)
4890
{
4891
  int i,j,sideelemtype2,elemind,parent,normal,elemtype;
4892
  int elemsides,side,sidenodes,nohits,hit,hit1,hit2;
4893
  int sideind2[MAXNODESD1];
4894
  int debug;
4895
  
4896
  hit = FALSE;
4897
  elemsides = 0;
4898
  elemtype = 0;
4899
  hit1 = FALSE;
4900
  hit2 = FALSE;
4901

4902
  debug = FALSE;
4903
  
4904
  for(parent=1;parent<=2;parent++) {
4905
    if(parent == 1) 
4906
      elemind = bound->parent[sideelem];
4907
    else
4908
      elemind = bound->parent2[sideelem];
4909
    
4910
    if(elemind > 0) {
4911
      elemtype = data->elementtypes[elemind];
4912
      elemsides = elemtype / 100;
4913

4914
      if(elemsides == 8) elemsides = 6;
4915
      else if(elemsides == 7) elemsides = 5;
4916
      else if(elemsides == 6) elemsides = 5;
4917
      else if(elemsides == 5) elemsides = 4;
4918
      
4919
      for(normal=1;normal >= -1;normal -= 2) {
4920

4921
	for(side=0;side<elemsides;side++) {
4922

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

4925
	  GetElementSide(elemind,side,normal,data,&sideind2[0],&sideelemtype2);
4926
	  
4927
	  if(sideelemtype2 < 300 && sideelemtype > 300) break;	
4928
	  if(sideelemtype2 < 200 && sideelemtype > 200) break;		
4929
	  if(sideelemtype != sideelemtype2) continue;
4930

4931
	  sidenodes = sideelemtype / 100;
4932

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

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

4956
      
4957
      /* this finding of sides does not guarantee that normals are oriented correctly */
4958
      normal = 1;
4959
      hit = FALSE;
4960
 
4961
      for(side=0;;side++) {
4962

4963
	if(0) printf("side = %d\n",side);
4964

4965
	GetElementSide(elemind,side,normal,data,&sideind2[0],&sideelemtype2);
4966

4967
	if(sideelemtype2 == 0 ) break;
4968
	if(sideelemtype2 < 300 && sideelemtype > 300) break;	
4969
	if(sideelemtype2 < 200 && sideelemtype > 200) break;		
4970
	if(sideelemtype != sideelemtype2) continue;
4971
		
4972
	sidenodes = sideelemtype % 100;
4973

4974
	nohits = 0;
4975
	for(j=0;j<sidenodes;j++) 
4976
	  for(i=0;i<sidenodes;i++) 
4977
	    if(sideind[i] == sideind2[j]) nohits++;
4978
	if(nohits == sidenodes) {
4979
	  hit = TRUE;
4980
	  if(parent == 1) {
4981
	    hit1 = TRUE;
4982
	    bound->side[sideelem] = side;
4983
	  }
4984
	  else {
4985
	    hit2 = TRUE;
4986
	    bound->side2[sideelem] = side;	      
4987
	  }
4988
	  goto skip;
4989
	}
4990
	
4991
      }
4992

4993
    skip:  
4994
      if(!hit) {
4995
	printf("FindParentSide: cannot locate BC element in parent %d: %d\n",parent,elemind);
4996
	printf("BC elem %d of type %d with corner indexes: ",sideelem,sideelemtype);
4997
	for(i=0;i<sideelemtype/100;i++)
4998
	  printf(" %d ",sideind[i]);
4999
	printf("\n");
5000

5001
	printf("Bulk elem %d of type %d with corner indexes: ",elemind,elemtype);
5002
	j = elemtype/100;
5003
	if( j >= 5 && j<=7 ) j = j-1;
5004
	for(i=0;i<j;i++)
5005
	  printf(" %d ",data->topology[elemind][i]);
5006
	printf("\n");             
5007
      }
5008
    }
5009
  }
5010

5011
  if(hit1 || hit2) 
5012
    return(0);
5013
  else
5014
    return(1);
5015
}
5016

5017

5018
static int Getnamerow(char *line,FILE *io,int upper) 
5019
{
5020
  int i,isend;
5021
  char *charend;
5022

5023

5024
  charend = fgets(line,MAXLINESIZE,io);
5025
  isend = (charend == NULL);
5026

5027
  if(isend) 
5028
    return(1);
5029
  else
5030
    return(0);
5031
}
5032

5033

5034

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

5051

5052
  sprintf(directoryname,"%s",prefix);
5053
  cdstat = chdir(directoryname);
5054

5055
  if(info) {
5056
    if(cdstat) 
5057
      printf("Loading mesh in ElmerSolver format from root directory.\n");
5058
    else
5059
      printf("Loading mesh in ElmerSolver format from directory %s.\n",directoryname);
5060
  }
5061

5062
  InitializeKnots(data);
5063

5064

5065
  sprintf(filename,"%s","mesh.header");
5066
  if ((in = fopen(filename,"r")) == NULL) {
5067
    printf("LoadElmerInput: The opening of the header-file %s failed!\n",
5068
	   filename);
5069
    return(1);
5070
  }
5071
  else 
5072
    printf("Loading header from %s\n",filename);
5073

5074
  GETLINE;
5075
  sscanf(line,"%d %d %d",&noknots,&noelements,&nosides);
5076
  GETLINE;
5077
  sscanf(line,"%d",&tottypes);
5078

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

5092
  data->dim = GetElementDimension(maxelemtype);
5093

5094
  data->maxnodes = maxnodes;
5095
  data->noknots = noknots;
5096
  data->noelements = noelements0 = noelements;
5097

5098

5099
  if(info) printf("Allocating for %d knots and %d elements.\n",
5100
		  noknots,noelements);
5101
  AllocateKnots(data);
5102

5103

5104
  sprintf(filename,"%s","mesh.nodes");
5105
  if ((in = fopen(filename,"r")) == NULL) {
5106
    if(info) printf("LoadElmerInput: The opening of the nodes-file %s failed!\n",
5107
		    filename);
5108
    bigerror("Cannot continue without nodes file!\n");
5109
  }
5110
  else 
5111
    printf("Loading %d Elmer nodes from %s\n",noknots,filename);
5112

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

5128

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

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

5244

5245

5246
  falseparents = 0;
5247
  noparents = 0;
5248
  bctopocreated = FALSE;
5249

5250
  sprintf(filename,"%s","mesh.boundary");
5251
  if ((in = fopen(filename,"r")) == NULL) {
5252
    printf("LoadElmerInput: The opening of the boundary-file %s failed!\n",
5253
	   filename);
5254
    return(4);
5255
  }
5256
  else {
5257
    if(info) printf("Loading %d boundary elements from %s\n",nosides,filename);
5258
  }
5259

5260
  if( nosides > 0 ) {
5261
    AllocateBoundary(bound,nosides);
5262
    data->noboundaries = 1;
5263
  };
5264

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

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

5344
      printf("elementtype = %d %d %d\n",i,elementtype,sideind[0]);
5345
      noparents++;
5346
    }
5347
  }
5348
  
5349
  if( falseparents ) {
5350
    printf("There seems to be %d false parents in the mesh\n",falseparents);
5351
  }
5352
  if( noparents ) {
5353
    printf("There seems to be %d bc elements without parents in the mesh\n",noparents);
5354
  }
5355

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

5373

5374
  sprintf(filename,"%s","mesh.names");
5375
  if ((in = fopen(filename,"r") )) {
5376
    int isbody,started,nameproblem;
5377
    
5378
    isbody = TRUE;
5379
    nameproblem = FALSE;
5380

5381
    if( nonames ) {
5382
      printf("Ignoring > mesh.names < because it was explicitly requested!\n");
5383
      goto namesend;
5384
    }
5385

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

5388
    for(;;) {
5389
      if(Getnamerow(line,in,FALSE)) goto namesend;
5390

5391
      if(strstr(line,"names for boundaries")) {
5392
	if(info) printf("Reading names for mesh boundaries\n");
5393
	isbody = FALSE;
5394
	continue;
5395
      }
5396
      else if(strstr(line,"names for bodies")) {
5397
	if(info) printf("Reading names for mesh bodies\n");	
5398
	isbody = TRUE;
5399
	continue;
5400
      }
5401

5402
      /* get position for entity name */
5403
      ptr1 = strchr( line,'$');
5404
      if(!ptr1) continue;
5405
      ptr1++;
5406

5407
      /* get position for entity index and read it */
5408
      ptr2 = strchr( line,'=');
5409
      if(!ptr2) continue;
5410
      ptr2++;      
5411
      j = next_int(&ptr2);
5412

5413
      /* Initialize the entity name by white spaces */
5414
      for(i=0;i<MAXLINESIZE;i++) 
5415
	line2[i] = ' ';
5416

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

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

5461
    if( nameproblem ) {
5462
      data->boundarynamesexist = FALSE;
5463
      data->bodynamesexist = FALSE;
5464
      printf("Number of indexes beyond range: %d\n",nameproblem);
5465
      smallerror("Omitting use of names because some indexes are beyond range, code some more...");
5466
    }
5467
  }
5468

5469

5470
  if(!cdstat) cdstat = chdir("..");
5471

5472
  if(info) printf("Elmer mesh loaded successfully\n");
5473

5474
  return(0);
5475
}
5476

5477

5478
int SaveElmerInput(struct FemType *data,struct BoundaryType *bound,
5479
		   char *prefix,int decimals,int nooverwrite, int info)
5480
/* Saves the mesh in a form that may be used as input 
5481
   in Elmer calculations. 
5482
   */
5483
{
5484
  int noknots,noelements,material,sumsides,elemtype,fail,cdstat,bcdim;
5485
  int sideelemtype,conelemtype,nodesd1,nodesd2,newtype;
5486
  int i,j,k,l,bulktypes[MAXELEMENTTYPE+1],sidetypes[MAXELEMENTTYPE+1];
5487
  int alltypes[MAXELEMENTTYPE+1],tottypes;
5488
  int ind[MAXNODESD1],ind2[MAXNODESD1],usedbody[MAXBODIES],usedbc[MAXBCS];
5489
  FILE *out;
5490
  char filename[MAXFILESIZE], outstyle[MAXFILESIZE];
5491
  char directoryname[MAXFILESIZE];
5492

5493
  if(!data->created) {
5494
    printf("You tried to save points that were never created.\n");
5495
    return(1);
5496
  }
5497

5498
  noelements = data->noelements;
5499
  noknots = data->noknots;
5500
  sumsides = 0;
5501

5502
  for(i=0;i<=MAXELEMENTTYPE;i++)
5503
    alltypes[i] = bulktypes[i] = sidetypes[i] = 0;
5504

5505
  for(i=0;i<MAXBODIES;i++)
5506
    usedbody[i] = 0;
5507
  for(i=0;i<MAXBCS;i++)
5508
    usedbc[i] = 0;
5509

5510
  sprintf(directoryname,"%s",prefix);
5511

5512
  if(info) printf("Saving mesh in ElmerSolver format to directory %s.\n",
5513
		  directoryname);
5514

5515
  fail = chdir(directoryname);
5516
  if(fail) {
5517
#ifdef MINGW32
5518
    fail = mkdir(directoryname);
5519
#else
5520
    fail = mkdir(directoryname,0750);
5521
#endif
5522
    if(fail) {
5523
      printf("Could not create a result directory!\n");
5524
      return(1);
5525
    }
5526
    else {
5527
      cdstat = chdir(directoryname);
5528
    }
5529
  }
5530
  else {
5531
    if(info) printf("Reusing an existing directory\n");
5532
    if(nooverwrite) {
5533
      if ((out = fopen("mesh.header", "r"))) {
5534
	printf("Mesh seems to already exist, writing is cancelled!\n"); 
5535
	return(1);
5536
      }
5537
    }
5538
  }
5539

5540
  
5541
  sprintf(filename,"%s","mesh.nodes");
5542
  out = fopen(filename,"w");
5543

5544
  if(info) printf("Saving %d coordinates to %s.\n",noknots,filename);  
5545
  if(out == NULL) {
5546
    printf("opening of file was not successful\n");
5547
    return(2);
5548
  }
5549

5550
  sprintf(outstyle,"%%d %%d %%.%dg %%.%dg %%.%dg\n",decimals,decimals,decimals);
5551
  for(i=1; i <= noknots; i++) 
5552
    fprintf(out,outstyle,i,-1,data->x[i],data->y[i],data->z[i]);    
5553

5554
  fclose(out);
5555

5556
  sprintf(filename,"%s","mesh.elements");
5557
  out = fopen(filename,"w");
5558
  if(info) printf("Saving %d element topologies to %s.\n",noelements,filename);
5559
  if(out == NULL) {
5560
    printf("opening of file was not successful\n");
5561
    return(3);
5562
  }
5563

5564
  for(i=1;i<=noelements;i++) {
5565
    elemtype = data->elementtypes[i];
5566
    material = data->material[i];
5567

5568
    if(material < MAXBODIES) usedbody[material] += 1;
5569
    fprintf(out,"%d %d %d",i,material,elemtype);
5570

5571
    bulktypes[elemtype] += 1;
5572
    nodesd2 = elemtype%100;
5573
    for(j=0;j < nodesd2;j++) 
5574
      fprintf(out," %d",data->topology[i][j]);
5575
    fprintf(out,"\n");          
5576
  }
5577
  fclose(out);
5578

5579

5580
  sprintf(filename,"%s","mesh.boundary");
5581
  out = fopen(filename,"w");
5582
  if(info) printf("Saving boundary elements to %s.\n",filename);
5583
  if(out == NULL) {
5584
    printf("opening of file was not successful\n");
5585
    return(4);
5586
  }
5587

5588
  sumsides = 0;
5589

5590

5591
  /* Save normal boundaries */
5592
  for(j=0;j < MAXBOUNDARIES;j++) {
5593
    if(bound[j].created == FALSE) continue;
5594
    if(bound[j].nosides == 0) continue;
5595
    
5596
    for(i=1; i <= bound[j].nosides; i++) {
5597
      GetBoundaryElement(i,&bound[j],data,ind,&sideelemtype); 
5598
      sumsides++;
5599
      
5600
      fprintf(out,"%d %d %d %d ",
5601
	      sumsides,bound[j].types[i],bound[j].parent[i],bound[j].parent2[i]);
5602
      fprintf(out,"%d",sideelemtype);
5603

5604
      k = bound[j].types[i];
5605
      if(k < MAXBCS) {	
5606
	bcdim = GetElementDimension(sideelemtype);
5607
	usedbc[k] = MAX(usedbc[k],bcdim+1);
5608
      }
5609
	
5610
      sidetypes[sideelemtype] += 1;
5611
      nodesd1 = sideelemtype%100;
5612
      for(l=0;l<nodesd1;l++)
5613
	fprintf(out," %d",ind[l]);
5614
      fprintf(out,"\n");
5615
    }
5616
  }
5617

5618
  newtype = 0;
5619
  for(j=0;j < MAXBOUNDARIES;j++) {
5620
    if(bound[j].created == FALSE) continue;
5621
    for(i=1; i <= bound[j].nosides; i++) 
5622
      newtype = MAX(newtype, bound[j].types[i]);
5623
  }  
5624
  fclose(out);
5625

5626
  tottypes = 0;
5627
  for(i=0;i<=MAXELEMENTTYPE;i++) {
5628
    alltypes[i] = bulktypes[i] + sidetypes[i];
5629
    if(alltypes[i]) tottypes++;
5630
  }
5631

5632
  sprintf(filename,"%s","mesh.header");
5633
  out = fopen(filename,"w");
5634
  if(info) printf("Saving header info to %s.\n",filename);  
5635
  if(out == NULL) {
5636
    printf("opening of file was not successful\n");
5637
    return(4);
5638
  }
5639
  fprintf(out,"%-6d %-6d %-6d\n",
5640
	  noknots,noelements,sumsides);
5641
  fprintf(out,"%-6d\n",tottypes);
5642
  for(i=0;i<=MAXELEMENTTYPE;i++) {
5643
    if(alltypes[i]) 
5644
      fprintf(out,"%-6d %-6d\n",i,bulktypes[i]+sidetypes[i]);
5645
  }
5646
  fclose(out);
5647

5648

5649
  if(data->boundarynamesexist || data->bodynamesexist) {
5650
    sprintf(filename,"%s","mesh.names");
5651
    out = fopen(filename,"w");
5652
    if(info) printf("Saving names info to %s.\n",filename);  
5653
    if(out == NULL) {
5654
      printf("opening of file was not successful\n");
5655
      return(5);
5656
    }
5657
    
5658
    if(data->bodynamesexist) {
5659
      fprintf(out,"! ----- names for bodies -----\n");
5660
      for(i=1;i<MAXBODIES;i++)
5661
	if(usedbody[i]) {
5662
	  if(data->bodyname[i])
5663
	    fprintf(out,"$ %s = %d\n",data->bodyname[i],i);
5664
	  else
5665
	    fprintf(out,"$ body%d = %d\n",i,i);	    
5666
	}
5667
    }     
5668
    if(data->boundarynamesexist) {
5669
      fprintf(out,"! ----- names for boundaries -----\n");
5670
      for(i=1;i<MAXBCS;i++) 	
5671
	if(usedbc[i]) {
5672
	  bcdim = usedbc[i]-1;
5673
	  if(data->boundaryname[i]) 
5674
	    fprintf(out,"$ %s = %d\n",data->boundaryname[i],i);
5675
	  else if(bcdim == 2) { 
5676
	    fprintf(out,"$ surf_bc%d = %d\n",i,i);
5677
	  }
5678
	  else if(bcdim == 1) {
5679
	    fprintf(out,"$ line_bc%d = %d\n",i,i);
5680
	  }
5681
	  else if(bcdim == 0) {
5682
	    fprintf(out,"$ node_bc%d = %d\n",i,i);
5683
	  }
5684
	}
5685
    }
5686
    fclose(out);
5687

5688
    sprintf(filename,"%s","entities.sif");
5689
    out = fopen(filename,"w");
5690
    if(info) printf("Saving entities info to %s.\n",filename);  
5691
    if(out == NULL) {
5692
      printf("opening of file was not successful\n");
5693
      return(5);
5694
    }
5695

5696
    if(data->bodynamesexist) {
5697
      fprintf(out,"!------ Skeleton for body section -----\n");
5698
      j = 0;
5699
      for(i=1;i<MAXBODIES;i++) {
5700
	if(usedbody[i]) {
5701
	  j = j + 1;
5702
	  fprintf(out,"Body %d\n",j);
5703
	  if(data->bodyname[i])
5704
	    fprintf(out,"  Name = %s\n",data->bodyname[i]);
5705
	  else
5706
	    fprintf(out,"  Name = body%d\n",i);	    
5707
	  fprintf(out,"End\n\n");
5708
	}
5709
      }
5710
    }
5711

5712
    if(data->boundarynamesexist) {
5713
      fprintf(out,"!------ Skeleton for boundary section -----\n");
5714
      j = 0;
5715
      for(i=1;i<MAXBCS;i++) {
5716
	if(usedbc[i]) {
5717
	  j = j + 1;
5718
	  fprintf(out,"Boundary Condition %d\n",j);
5719
	  if(data->boundaryname[i])
5720
	    fprintf(out,"  Name = %s\n",data->boundaryname[i]);
5721
	  else
5722
	    fprintf(out,"  Name = bc%d\n",i);	    
5723
	  fprintf(out,"End\n\n");
5724
	}
5725
      }
5726
    }
5727
    fclose(out);
5728
  }
5729
  
5730
  if(data->nodepermexist) {
5731
    sprintf(filename,"%s","mesh.nodeperm");
5732
    out = fopen(filename,"w");
5733

5734
    if(info) printf("Saving initial node permutation to %s.\n",filename);  
5735
    if(out == NULL) {
5736
      printf("opening of file was not successful\n");
5737
      return(3);
5738
    }
5739
    for(i=1; i <= noknots; i++) 
5740
      fprintf(out,"%d %d\n",i,data->nodeperm[i]);
5741
  }
5742

5743

5744
  cdstat = chdir("..");
5745
  
5746
  return(0);
5747
}
5748

5749

5750
int CreateElmerGridMesh(struct GridType *grid,
5751
			struct FemType *data,struct BoundaryType *boundaries,
5752
			Real relh,int info) {
5753
  int j;  
5754
  struct CellType *cell;
5755
  
5756
  for(j=0;j<MAXBOUNDARIES;j++) {
5757
    boundaries[j].created = FALSE;
5758
    boundaries[j].nosides = FALSE;
5759
  }
5760

5761
  SetElementDivision(grid,relh,info);
5762
  
5763
  cell = (struct CellType*)
5764
    malloc((size_t) (grid->nocells+1)*sizeof(struct CellType)); 
5765
  SetCellData(grid,cell,info);
5766

5767
  if(grid->dimension == 1) 
5768
    SetCellKnots1D(grid,cell,info);
5769
  else
5770
    SetCellKnots(grid,cell,info);
5771

5772
  CreateKnots(grid,cell,data,0,0);
5773

5774
  if(grid->noboundaries > 0) {
5775
    for(j=0;j<grid->noboundaries;j++) {
5776
      if(grid->boundsolid[j] < 4) {
5777
	CreateBoundary(cell,data,&(boundaries[j]),grid->boundext[j],grid->boundint[j],
5778
		       1,grid->boundtype[j],info);  
5779
      } 
5780
      else { 
5781
	CreatePoints(cell,data,&(boundaries[j]),grid->boundext[j],grid->boundint[j],
5782
		     grid->boundsolid[j],grid->boundtype[j],info); 	    
5783
      }
5784
    }
5785
  }
5786

5787
  free(cell);
5788

5789
  return 0;
5790
}
5791

5792

5793
Product

Resources

Company
