egconvert.cpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file egconvert.cpp

/*  
   ElmerGrid - A simple mesh generation and manipulation utility  
   Copyright (C) 1995- , CSC - IT Center for Science Ltd.   
 
   Author: Peter Råback
   Email: Peter.Raback@csc.fi
   Address: CSC - IT Center for Science Ltd.
            Keilaranta 14
            02101 Espoo, Finland

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

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

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

/* --------------------:  egconvert.c  :-------------------------- */

#include <stdio.h>
#include <math.h>
#include <stdarg.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>

#include "egutils.h"
#include "egdef.h"
#include "egtypes.h"
#include "egmesh.h"
#include "egnative.h"
#include "egconvert.h"

#define getline fgets(line,MAXLINESIZE,in) 


static int Getrow(char *line1,FILE *io,int upper) 
{
  int i,isend;
  char line0[MAXLINESIZE],*charend;

  for(i=0;i<MAXLINESIZE;i++) 
    line0[i] = ' ';

 newline:

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

  if(isend) return(1);

  if(line0[0] == '#' || line0[0] == '!') goto newline;
  if(strstr(line0,"#")) goto newline;

  if(upper) {
    for(i=0;i<MAXLINESIZE;i++) 
      line1[i] = toupper(line0[i]);
  }
  else {
    for(i=0;i<MAXLINESIZE;i++) 
      line1[i] = line0[i];    
  }

  return(0);
}




static int GetrowDouble(char *line1,FILE *io)
{
  int i,isend;
  char line0[MAXLINESIZE],*charend;

  for(i=0;i<MAXLINESIZE;i++) 
    line0[i] = ' ';

 newline:

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

  if(isend) return(1);

  if(line0[0] == '#' || line0[0] == '!') goto newline;
  if(strstr(line0,"#")) goto newline;

  for(i=0;i<MAXLINESIZE;i++) { 

    /* The fortran double is not recognized by C string operators */
    if( line0[i] == 'd' || line0[i] == 'D' ) {
      line1[i] = 'e';
    } else {
      line1[i] = line0[i];    
    }
  }

  return(0);
}


static int Comsolrow(char *line1,FILE *io) 
{
  int i,isend;
  char line0[MAXLINESIZE],*charend;

  for(i=0;i<MAXLINESIZE;i++) 
    line0[i] = ' ';

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

  if(isend) return(1);

  for(i=0;i<MAXLINESIZE;i++) line1[i] = line0[i];    

  return(0);
}



static void FindPointParents(struct FemType *data,struct BoundaryType *bound,
          int boundarynodes,int *nodeindx,int *boundindx,int info)
{
  int i,j,k,sideelemtype,elemind,*indx;
  int boundarytype,minboundary,maxboundary,minnode,maxnode,sideelem,elemtype;
  int sideind[MAXNODESD1],elemsides,side,sidenodes,hit,nohits;
  int *elemhits;

  info = TRUE;

  sideelem = 0;
  maxboundary = minboundary = boundindx[1];
  minnode = maxnode = nodeindx[1]; 

  for(i=1;i<=boundarynodes;i++) {
    if(maxboundary < boundindx[i]) maxboundary = boundindx[i];
    if(minboundary > boundindx[i]) minboundary = boundindx[i];
    if(maxnode < nodeindx[i]) maxnode = nodeindx[i];
    if(minnode > nodeindx[i]) minnode = nodeindx[i];
  }

  if(info) {
    printf("Boundary types are in interval [%d, %d]\n",minboundary,maxboundary);
    printf("Boundary nodes are in interval [%d, %d]\n",minnode,maxnode);
  }
  indx = Ivector(1,data->noknots);

  elemhits = Ivector(1,data->noknots);
  for(i=1;i<=data->noknots;i++) elemhits[i] = 0;

  for(elemind=1;elemind<=data->noelements;elemind++) {
    elemtype = data->elementtypes[elemind];
    elemsides = elemtype % 100;

    for(i=0;i<elemsides;i++) {
      elemhits[data->topology[elemind][i]] += 1;
    }
  }


  for(boundarytype=minboundary;boundarytype <= maxboundary;boundarytype++) {
    int boundfirst,bchits,bcsame,sideelemtype2;
    int sideind2[MAXNODESD1];

    boundfirst = 0;
  
    for(i=1;i<=data->noknots;i++) 
      indx[i] = 0;

    for(i=1;i<=boundarynodes;i++) {
      if(boundindx[i] == boundarytype) 
  indx[nodeindx[i]] = TRUE; 
    }


    for(elemind=1;elemind<=data->noelements;elemind++) {
      elemtype = data->elementtypes[elemind];
      elemsides = elemtype / 100;
      if(elemsides == 8) elemsides = 6;
      else if(elemsides == 5) elemsides = 4;
      else if(elemsides == 6) elemsides = 5;
      
      if(0) printf("ind=%d  type=%d  sides=%d\n",elemind,elemtype,elemsides);
 
      /* Check whether the bc nodes occupy every node in the selected side */
      for(side=0;side<elemsides;side++) {
  GetElementSide(elemind,side,1,data,&sideind[0],&sideelemtype);
  sidenodes = sideelemtype%100;
  
  if(0) printf("sidenodes=%d  side=%d\n",sidenodes,side);
  
  hit = TRUE;
  nohits = 0;
  for(i=0;i<sidenodes;i++) {
    if(sideind[i] <= 0) { 
      if(0) printf("sideind[%d] = %d\n",i,sideind[i]);
      hit = FALSE;
    }
    else if(sideind[i] > data->noknots) { 
      if(0) printf("sideind[%d] = %d (noknots=%d)\n",i,sideind[i],data->noknots);
      hit = FALSE;
    }
    else if(!indx[sideind[i]]) {
      hit = FALSE;
    }
    else {
      nohits++;
    }
  }
  

  if(hit == TRUE) {

    /* If all the points belong to more than one element there may be 
       another parent for the side element */
    bcsame = FALSE;
    bchits = TRUE;
    for(i=0;i<sidenodes;i++) 
      if(elemhits[sideind[i]] < 2) bchits = FALSE;

    if(bchits && boundfirst) {
      for(j=boundfirst;j<=sideelem;j++) {
        if(bound->parent2[j]) continue;
        GetElementSide(bound->parent[j],bound->side[j],1,data,&sideind2[0],&sideelemtype2);
        if(sideelemtype != sideelemtype2) continue;
        bcsame = 0;
        for(i=0;i<sidenodes;i++) 
    for(k=0;k<sidenodes;k++)
      if(sideind[i] == sideind2[k]) bcsame++;

        if(bcsame == sidenodes) {
    if(data->material[bound->parent[j]] > data->material[elemind]) {
      bound->parent2[j] = bound->parent[j];
      bound->side2[j] = bound->side[j];
      bound->parent[j] = elemind;
      bound->side[j] = side;
    }
    else {      
      bound->parent2[j] = elemind;
      bound->side2[j] = side;
    }
    goto bcset;
        }
      }
    }
    
    sideelem += 1;
    bound->parent[sideelem] = elemind;
    bound->side[sideelem] = side;
    bound->parent2[sideelem] = 0;
    bound->side2[sideelem] = 0;
    bound->types[sideelem] = boundarytype;

    if(!boundfirst) boundfirst = sideelem;

  bcset:
    continue;
  }
      }
    }    
  }

  
  free_Ivector(indx,1,data->noknots);

  if(info) printf("Found %d side elements formed by %d points.\n",
      sideelem,boundarynodes);

  bound->nosides = sideelem;

  return;
}




int LoadAbaqusInput(struct FemType *data,struct BoundaryType *bound,
        char *prefix,int info)
/* Load the grid from a format that can be read by ABAQUS 
   program designed for sructural mechanics. The commands
   understood are only those that IDEAS creates when saving
   results in ABAQUS format.
   */
{
  int noknots,noelements,elemcode,maxnodes,material;
  int mode,allocated,nvalue,nvalue2,maxknot,nosides;
  int boundarytype,boundarynodes,elsetactive;
  int *nodeindx = NULL,*boundindx = NULL;
  
  char filename[MAXFILESIZE];
  char line[MAXLINESIZE];
  int i,j,*ind = NULL;
  FILE *in;
  Real rvalues[MAXDOFS];
  int ivalues[MAXDOFS],ivalues0[MAXDOFS];


  strcpy(filename,prefix);
  if ((in = fopen(filename,"r")) == NULL) {
    AddExtension(prefix,filename,"inp");
    if ((in = fopen(filename,"r")) == NULL) {
      printf("LoadAbaqusInput: opening of the ABAQUS-file '%s' wasn't succesfull !\n",
       filename);
      return(1);
    }
  }

  printf("Reading input from ABAQUS input file %s.\n",filename);
  InitializeKnots(data);

  allocated = FALSE;
  maxknot = 0;
  elsetactive = FALSE;

  /* Because the file format doesn't provide the number of elements
     or nodes the results are read twice but registered only in the
     second time. */

omstart:

  mode = 0;
  maxnodes = 0;
  noknots = 0;
  noelements = 0;
  elemcode = 0;
  boundarytype = 0;
  boundarynodes = 0;
  material = 0;
  ivalues0[0] = ivalues0[1] = 0;


  for(;;) {
    /* getline; */

    if (Getrow(line,in,TRUE)) goto end;

    /* if(!line) goto end; */
    /* if(strstr(line,"END")) goto end; */

    if(strstr(line,"**")) {
      if(info && !allocated) printf("comment: %s",line);
    }
    else if(strrchr(line,'*')) {
      if(strstr(line,"HEAD")) {
  mode = 1;
      }
      else if(strstr(line,"NODE")) {
  if(strstr(line,"SYSTEM=R")) data->coordsystem = COORD_CART2;
  if(strstr(line,"SYSTEM=C")) data->coordsystem = COORD_AXIS;      
  if(strstr(line,"SYSTEM=P")) data->coordsystem = COORD_POLAR;      
  mode = 2;
      }
      else if(strstr(line,"ELEMENT")) {
  if(!elsetactive) material++;
  if(strstr(line,"S3R") || strstr(line,"STRI3"))
    elemcode = 303;
  else if(strstr(line,"2D4") || strstr(line,"SP4") || strstr(line,"AX4") 
    || strstr(line,"S4") || strstr(line,"CPE4")) 
    elemcode = 404;
  else if(strstr(line,"2D8") || strstr(line,"AX8"))
    elemcode = 408;
  else if(strstr(line,"3D4"))
    elemcode = 504;
  else if(strstr(line,"3D5"))
    elemcode = 605;
  else if(strstr(line,"3D6"))
    elemcode = 706;
  else if(strstr(line,"3D8"))
    elemcode = 808;
  else if(strstr(line,"3D20"))
    elemcode = 820;
  else 
    printf("Unknown element code: %s\n",line);
  
  if(maxnodes < elemcode%100) maxnodes = elemcode%100;
  mode = 3;
  if(1) printf("Loading elements of type %d starting from element %d.\n",
      elemcode,noelements);
      }
      else if(strstr(line,"BOUNDARY") || strstr(line,"CLOAD")) {
  boundarytype++;
  mode = 4;
      }
      else if(strstr(line,"NSET")) {
  boundarytype++;
  mode = 5;
      }
      else if(strstr(line,"ELSET")) {
  elsetactive = TRUE;
  material += 1;
  mode = 6;
      }
      else {
  if(!allocated) printf("unknown command: %s",line);
  mode = 0;
      }
    }

    else if(mode) {  
      switch (mode) {
  
      case 1:
  if(info) printf("Loading Abacus input file:\n%s",line);
  break;
  
      case 2:
  nvalue = StringToReal(line,rvalues,MAXNODESD2+1,',');
  i = (int)(rvalues[0]+0.5);
  if(i == 0) continue;
  noknots++;
  if(allocated) {
    ind[i] = noknots;
    data->x[noknots] = rvalues[1];
    data->y[noknots] = rvalues[2];
    data->z[noknots] = rvalues[3];
  }
  else {
    if(maxknot < i) maxknot = i;
  }
  break;
  
      case 3:
  noelements++;

  nvalue = StringToInteger(line,ivalues,MAXNODESD2+1,',');
  
  if(allocated) {
    data->elementtypes[noelements] = elemcode;
    data->material[noelements] = material;
    for(i=0;i<nvalue-1;i++) 
      data->topology[noelements][i] = ivalues[i+1];   
  }
  
  if(nvalue < elemcode % 100) {
    Getrow(line,in,TRUE);
    if(allocated) {
      if(ivalues[nvalue-1] == 0) nvalue--;        
      nvalue2 = StringToInteger(line,ivalues,MAXNODESD2+1,',');
      for(i=0;i<nvalue2;i++) 
        data->topology[noelements][nvalue-1+i] = ivalues[i];          
    }
  }
  break;

      case 4:
  nvalue = StringToInteger(line,ivalues,MAXNODESD2+1,',');

  if(ivalues[0] == ivalues0[0] & ivalues[1] != ivalues0[1]) continue;
  ivalues0[0] = ivalues[0];
  ivalues0[1] = ivalues[1];

  boundarynodes++;
  if(allocated) {
    nodeindx[boundarynodes] = ivalues[0];
    boundindx[boundarynodes] = boundarytype;
  }
  break;

      case 5:
  nvalue = StringToInteger(line,ivalues,10,',');

  if(allocated) {
    for(i=0;i<nvalue;i++) {
      boundarynodes += 1;
      nodeindx[boundarynodes] = ivalues[i];
      boundindx[boundarynodes] = boundarytype;
    }
  }
  else
    boundarynodes += nvalue;
  break;

      case 6:
  nvalue = StringToInteger(line,ivalues,10,',');

  if(allocated) {
    for(i=0;i<nvalue;i++) {
      j = ivalues[i];
      data->material[j] = material;
    }
  }
  break;

      default:
  printf("Unknown case: %d\n",mode);
      }      
    }

  }    
  end:


  if(allocated == TRUE) {
    if(info) printf("The mesh was loaded from file %s.\n",filename);

    FindPointParents(data,bound,boundarynodes,nodeindx,boundindx,info);

    /* ABAQUS format does not expect that all numbers are used
       when numbering the elements. Therefore the nodes must
       be renumberred from 1 to noknots. */

    if(noknots != maxknot) {
      if(info) printf("There are %d nodes but maximum index is %d.\n",
          noknots,maxknot);
      if(info) printf("Renumbering elements\n");
      for(j=1;j<=noelements;j++) 
  for(i=0;i < data->elementtypes[j]%100;i++) 
    data->topology[j][i] = ind[data->topology[j][i]];
    }

    free_ivector(ind,1,maxknot);
    free_Ivector(nodeindx,1,boundarynodes);
    free_Ivector(boundindx,1,boundarynodes);

    fclose(in);
    return(0);
  }

  rewind(in);
  data->noknots = noknots;
  data->noelements = noelements;
  data->maxnodes = maxnodes;
  data->dim = 3;
  
  if(info) printf("Allocating for %d knots and %d %d-node elements.\n",
      noknots,noelements,maxnodes);
  AllocateKnots(data);

  nosides = 2*boundarynodes;
  printf("There are %d boundary nodes, thus allocating %d elements\n",
   boundarynodes,nosides);
  AllocateBoundary(bound,nosides);
  nodeindx = Ivector(1,boundarynodes);
  boundindx = Ivector(1,boundarynodes);

  ind = ivector(1,maxknot);
  for(i=1;i<=maxknot;i++)
    ind[i] = 0;
    
  allocated = TRUE;
  goto omstart;
}




int LoadNastranInput(struct FemType *data,struct BoundaryType *bound,
        char *prefix,int info)
/* Load the grid from a format that in Nastran format 
   */
{
  int noknots,noelements,elemcode,maxnodes,material;
  int mode,allocated,maxknot,minknot;
  int boundarytype,boundarynodes,nodes;
  
  char filename[MAXFILESIZE];
  char line[MAXLINESIZE],*cp;
  int j,k=0;
  FILE *in;
  int ivalues0[MAXDOFS];


  strcpy(filename,prefix);
  if ((in = fopen(filename,"r")) == NULL) {
    AddExtension(prefix,filename,"nas");
    if ((in = fopen(filename,"r")) == NULL) {
      printf("LoadNastranInput: opening of the Nastran file '%s' wasn't succesfull !\n",
       filename);
      return(1);
    }
  }

  if(info) printf("Reading mesh from Nastran file %s.\n",filename);
  InitializeKnots(data);

  allocated = FALSE;
  maxknot = 0;
  minknot = 10000;

  /* Because the file format doesn't provide the number of elements
     or nodes the results are read twice but registered only in the
     second time. */
omstart:

  mode = 0;
  maxnodes = 0;
  noknots = 0;
  noelements = 0;
  elemcode = 0;
  boundarytype = 0;
  boundarynodes = 0;
  material = 0;
  ivalues0[0] = ivalues0[1] = 0;


  for(;;) {
    /* getline; */

    if (Getrow(line,in,TRUE)) goto end;

    if(line[0] == '$') {
      if(!allocated) printf("comment: %s",line);
    }
    else if(strstr(line,"GRID")) {
      noknots++;
      if(0) printf("line=%s\n",line);
      cp = &line[5];
      j = next_int(&cp);
      if(0) printf("j=%d\n",j);

      if(allocated) {
  k = next_int(&cp);
  data->x[noknots] = next_real(&cp);
  data->y[noknots] = next_real(&cp);
      }
      else {
  if(j > maxknot) maxknot = j;
  if(j < minknot) minknot = j;
      }

      if(strstr(line,"*")) 
  Getrow(line,in,TRUE);

      if(allocated) {
  cp = &line[4];
  data->z[noknots] = next_real(&cp);
      }

    }
    else if(strstr(line,"TETRA")) {
      noelements++;
      nodes = 4;
      if(nodes > maxnodes) maxnodes = nodes;
      if(allocated) {
  data->elementtypes[noelements] = 504;
  cp = &line[6];
  k = next_int(&cp);
  data->material[noelements] = next_int(&cp) + 1; 
  for(j=0;j<nodes;j++)
    data->topology[noelements][j] = next_int(&cp);            
      }      
    }
    else if(strstr(line,"PYRAM")) {
      noelements++;
      nodes = 5;
      if(nodes > maxnodes) maxnodes = nodes;
      if(allocated) {
  data->elementtypes[noelements] = 605;
  cp = &line[6];
  k = next_int(&cp);
  data->material[noelements] = next_int(&cp) + 1; 
  for(j=0;j<nodes;j++)
    data->topology[noelements][j] = next_int(&cp);            
      }      
    }
    else if(strstr(line,"PENTA")) {
      noelements++;
      nodes = 6;
      if(nodes > maxnodes) maxnodes = nodes;
      if(allocated) {
  data->elementtypes[noelements] = 706;
  cp = &line[6];
  k = next_int(&cp);
  data->material[noelements] = next_int(&cp) + 1; 
  for(j=0;j<nodes;j++)
    data->topology[noelements][j] = next_int(&cp);            
      }      
    }
     else if(strstr(line,"CHEXA")) {
      noelements++;
      nodes = 8;
      if(nodes > maxnodes) maxnodes = nodes;
      if(allocated) {
  data->elementtypes[noelements] = 808;
  cp = &line[5];
  k = next_int(&cp);
  data->material[noelements] = next_int(&cp) + 1; 
  for(j=0;j<6;j++)
    data->topology[noelements][j] = next_int(&cp);            
      }           
      Getrow(line,in,TRUE);
      if(allocated) {
  cp = &line[1];       
  for(j=6;j<8;j++) 
    data->topology[noelements][j] = k;
      }
    }
    else if(strstr(line,"ENDDAT")) {
      goto end;
    }
    else {
      printf("unknown command: %s",line);
    }
  }

  end:

  if(allocated == TRUE) {
    if(info) printf("The mesh was loaded from file %s.\n",filename);
    fclose(in);
    return(0);
  }

  rewind(in);
  data->noknots = noknots;
  data->noelements = noelements;
  data->maxnodes = maxnodes;
  data->dim = 3;
  
  printf("maxknot = %d  minknot = %d\n",maxknot,minknot);
  
  if(info) printf("Allocating for %d knots and %d %d-node elements.\n",
      noknots,noelements,maxnodes);
  AllocateKnots(data);
    
  allocated = TRUE;
  goto omstart;
}




static void ReorderFidapNodes(struct FemType *data,int element,int nodes,int typeflag) 
{
  int i,oldtopology[MAXNODESD2],*topology,dim;
  int order203[]={1,3,2};
  int order306[]={1,3,5,2,4,6};
  int order408[]={1,3,5,7,2,4,6,8};
  int order409[]={1,3,5,7,2,4,6,8,9};
  int order510[]={1,3,6,10,2,5,4,7,8,9};
  int order605[]={1,2,4,3,5};
  int order808[]={1,2,4,3,5,6,8,7};
  int order827[]={1,3,9,7,19,21,27,25,2,6,8,4,10,12,18,16,20,24,26,22,11,15,17,13,5,23,14};

  dim = data->dim;
  if(typeflag > 10) dim -= 1;

  data->elementtypes[element] = 101*nodes;
  topology = data->topology[element];

  if(dim == 1) {
    if(nodes == 3) {
      data->elementtypes[element] = 203;
      for(i=0;i<nodes;i++) 
  oldtopology[i] = topology[i];
      for(i=0;i<nodes;i++) 
  topology[i] = oldtopology[order203[i]-1];            
    }
  }
  else if(dim == 2) {
    if(nodes == 6) {
      data->elementtypes[element] = 306;
      for(i=0;i<nodes;i++) 
  oldtopology[i] = topology[i];
      for(i=0;i<nodes;i++) 
  topology[i] = oldtopology[order306[i]-1];      
    }
    else if(nodes == 8) {
      data->elementtypes[element] = 408;
      for(i=0;i<nodes;i++) 
  oldtopology[i] = topology[i];
      for(i=0;i<nodes;i++) 
  topology[i] = oldtopology[order408[i]-1];      
    }
    else if(nodes == 9) {
      data->elementtypes[element] = 409;
      for(i=0;i<nodes;i++) 
  oldtopology[i] = topology[i];
      for(i=0;i<nodes;i++) 
  topology[i] = oldtopology[order409[i]-1];      
    }
  }
  else if(dim == 3) {
    if(nodes == 4) {
      data->elementtypes[element] = 504;      
    }
    else if(nodes == 10) {
      data->elementtypes[element] = 510;      
      for(i=0;i<nodes;i++) 
  oldtopology[i] = topology[i];
      for(i=0;i<nodes;i++) 
  topology[i] = oldtopology[order510[i]-1];      
    }
    else if(nodes == 5) {
      data->elementtypes[element] = 605;      
      for(i=0;i<nodes;i++) 
  oldtopology[i] = topology[i];
      for(i=0;i<nodes;i++) 
  topology[i] = oldtopology[order605[i]-1];
    }
    else if(nodes == 6) {
      data->elementtypes[element] = 706;      
    }
    else if(nodes == 8) {
      for(i=0;i<nodes;i++) 
  oldtopology[i] = topology[i];
      for(i=0;i<nodes;i++) 
  topology[i] = oldtopology[order808[i]-1];
      data->elementtypes[element] = 808;            
    }
    else if(nodes == 27) {
      for(i=0;i<nodes;i++) 
  oldtopology[i] = topology[i];
      for(i=0;i<nodes;i++) 
  topology[i] = oldtopology[order827[i]-1];
      data->elementtypes[element] = 827;                  
    }
    else {
      printf("Unknown Fidap elementtype with %d nodes.\n",nodes);
    }

  }
  else printf("ReorderFidapNodes: unknown dimension (%d)\n",data->dim);
  if(0) printf("dim = %d element = %d elemtype = %d\n",dim,element,data->elementtypes[element]);
}




int LoadFidapInput(struct FemType *data,struct BoundaryType *boundaries,char *prefix,int info)
/* Load the grid from a format that can be read by FIDAP 
   program designed for fluid mechanics. 

   Still under implementation
   */
{
  int noknots,noelements,dim,novel,elemcode,maxnodes;
  int mode,maxknot,totelems,entity,maxentity;
  char filename[MAXFILESIZE];
  char line[MAXLINESIZE],entityname[MAXNAMESIZE];
  int i,j,k,*ind,geoflag,typeflag;
  int **topology;
  FILE *in;
  Real *vel,*temp;
  int nogroups,knotno;
  char *isio;

  AddExtension(prefix,filename,"fidap");

  if ((in = fopen(filename,"r")) == NULL) {
    AddExtension(prefix,filename,"FDNEUT");
    if ((in = fopen(filename,"r")) == NULL) {
      printf("LoadFidapInput: opening of the Fidap-file '%s' wasn't succesfull !\n",
       filename);
      return(1);
    }
  }
 
  InitializeKnots(data);

  entity = 0;
  mode = 0;
  noknots = 0;
  noelements = 0;
  dim = 0;
  elemcode = 0;
  maxnodes = 4;
  totelems = 0;
  maxentity = 0;

  for(;;) {
    isio = getline;

    if(!isio) goto end;
    if(!line) goto end;
    if(line=="") goto end;
    if(strstr(line,"END")) goto end;

    /* Control information */
    if(strstr(line,"FIDAP NEUTRAL FILE")) mode = 1;
    else if(strstr(line,"NO. OF NODES")) mode = 2;
    else if(strstr(line,"TEMPERATURE/SPECIES FLAGS")) mode = 3;
    else if(strstr(line,"PRESSURE FLAGS")) mode = 4;
    else if(strstr(line,"NODAL COORDINATES")) mode = 5;
    else if(strstr(line,"BOUNDARY CONDITIONS")) mode = 6;
    else if(strstr(line,"ELEMENT GROUPS")) mode = 7;
    else if(strstr(line,"GROUP:")) mode = 8;
    else if(strstr(line,"VELOCITY")) mode = 10;
    else if(strstr(line,"TEMPERATURE")) mode = 11;
    else if(0) printf("unknown: %s",line);
    
    switch (mode) {

    case 1: 
      if(info) printf("Loading FIDAP input file %s\n",filename);
      getline;
      if(info) printf("Name of the case: %s",line);
      mode = 0;
      break;

    case 2:
      getline;   
      if(0) printf("Reading the header info\n");
      sscanf(line,"%d%d%d%d%d",&noknots,&noelements,
       &nogroups,&dim,&novel);
      data->noknots = noknots;
      data->noelements = noelements;
      data->maxnodes = maxnodes;
      data->dim = dim;

      mode = 0;
      break;
      
    case 5:
      if(info) printf("Allocating for %d knots and %d %d-node elements.\n",
          noknots,noelements,maxnodes);
      AllocateKnots(data);
      if(info) printf("Reading the nodes\n");
      for(i=1;i<=noknots;i++) {
  getline;
  if (dim == 2)
    sscanf(line,"%d%le%le",&knotno,
     &(data->x[i]),&(data->y[i]));
  else if(dim==3) 
    sscanf(line,"%d%le%le%le",&knotno,
     &(data->x[i]),&(data->y[i]),&(data->z[i]));
      }
      break;
      
    case 8: 
      {
  int val,group,elems,nodes;
  char *cp;

  i=0;
  do val=line[i++];
  while(val!=':');i++;
  sscanf(&line[i],"%d",&group);
  
  do val=line[i++];
  while(val!=':');i++;
  sscanf(&line[i],"%d",&elems);
  
  do val=line[i++];
  while(val!=':');i++;
  sscanf(&line[i],"%d",&nodes);

  do val=line[i++];
  while(val!=':');i++;
  sscanf(&line[i],"%d",&geoflag);

  do val=line[i++];
  while(val!=':');i++;
  sscanf(&line[i],"%d",&typeflag);
  
  getline;
  i=0;
  do val=line[i++];
  while(val!=':');i++;
  sscanf(&line[i],"%s",entityname);

  if(nodes > maxnodes) {
    if(info) printf("Allocating a %d-node topology matrix\n",nodes);
    topology = Imatrix(1,noelements,0,nodes-1);
    if(totelems > 0) {
      for(j=1;j<=totelems;j++) 
        for(i=0;i<data->elementtypes[j] % 100;i++)
    topology[j][i] = data->topology[j][i];
    }
    free_Imatrix(data->topology,1,data->noelements,0,data->maxnodes-1);
    data->maxnodes = maxnodes = nodes;
    data->topology = topology;
  }

  if(0) printf("Reading %d element topologies with %d nodes for %s\n",
      elems,nodes,entityname);

  for(entity=1;entity<=maxentity;entity++) {
    k = strcmp(entityname,data->bodyname[entity]);
    if(k == 0) break;
  }

  if(entity > maxentity) {
    maxentity++;
    strcpy(data->bodyname[entity],entityname);
    if(info) printf("Found new entity: %s\n",entityname);
  }

  for(i=totelems+1;i<=totelems+elems;i++) {
    getline;

    cp = line;
    j = next_int(&cp);
          for(j=0;j<nodes;j++)
      data->topology[i][j] = next_int(&cp);

    ReorderFidapNodes(data,i,nodes,typeflag);

    if(data->elementtypes[i] == 0) {
      printf("Elementtype is zero!\n");
    }

    if(entity) data->material[i] = entity;
    else data->material[i] = group;
  }
  totelems += elems;
      }
    mode = 0;
    break;
      
    case 10:
      dim = 3;
      if(info) printf("Reading the velocity field\n");
      CreateVariable(data,2,dim,0.0,"Velocity",FALSE);
      vel = data->dofs[2];
      for(j=1;j<=noknots;j++) {
  getline;
  if(dim==2) 
    sscanf(line,"%le%le",&(vel[2*j-1]),&(vel[2*j]));
  if(dim==3) 
    sscanf(line,"%le%le%le",&(vel[3*j-2]),&(vel[3*j-1]),&(vel[3*j]));
      }
      mode = 0;
      break;
      
    case 11:
      if(info) printf("Reading the temperature field\n");
      CreateVariable(data,1,1,0.0,"Temperature",FALSE);
      temp = data->dofs[1];
      for(j=1;j<=noknots;j++) {
  getline;
  sscanf(line,"%le",&(temp[j]));
      }      
      mode = 0;
      break;
      
    default:      
      break;
    }
  }    

end:

  /* Renumber the nodes */
  maxknot = 0;
  for(i=1;i<=noelements;i++) 
    for(j=0;j < data->elementtypes[i]%100;j++) 
      if(data->topology[i][j] > maxknot) maxknot = data->topology[i][j];

  if(maxknot > noknots) {
    if(info) printf("Renumbering the nodes from 1 to %d\n",noknots);

    ind = ivector(1,maxknot);
    for(i=1;i<=maxknot;i++)
      ind[i] = 0;

    for(i=1;i<=noelements;i++) 
      for(j=0;j < data->elementtypes[i]%100;j++) 
  ind[data->topology[i][j]] = data->topology[i][j];
    i=0;
    for(j=1;j<=noknots;j++) {
      i++;
      while(ind[i]==0) i++;
      ind[i] = j;
    } 
    for(i=1;i<=noelements;i++) 
      for(j=0;j < data->elementtypes[i]%100;j++) 
  data->topology[i][j] = ind[data->topology[i][j]];
  }


  if(maxentity > 0) data->bodynamesexist = TRUE;

  fclose(in);
  
  if(info) printf("Finished reading the Fidap neutral file\n");


  ElementsToBoundaryConditions(data,boundaries,FALSE,TRUE);
  RenumberBoundaryTypes(data,boundaries,TRUE,0,info);  

  return(0);
}



static void ReorderAnsysNodes(struct FemType *data,int *oldtopology,
            int element,int dim,int nodes) 
{
  int i,*topology,elementtype = 0;
  int order820[]={1,2,3,4,5,6,7,8,9,10,11,12,17,18,19,20,13,14,15,16};
  int order504[]={1,2,3,5};
  int order306[]={1,2,3,5,6,8};
  int order510[]={1,2,3,5,9,10,12,17,18,19};
  int order613[]={1,2,3,4,5,9,10,11,12,17,18,19,20};
    
  if(dim == 3) {
    if(nodes == 20) {
      if(oldtopology[2] == oldtopology[3]) elementtype = 510;
      else if(oldtopology[4] == oldtopology[5]) elementtype = 613;
      else elementtype = 820;
    }
    if(nodes == 8) {
      if(oldtopology[2] == oldtopology[3]) elementtype = 504;
      else if(oldtopology[4] == oldtopology[5]) elementtype = 605;
      else elementtype = 808;
    }
    if(nodes == 4) elementtype = 504;
    if(nodes == 10) elementtype = 510;
  }
  else if(dim == 2) {
    if(nodes == 9) elementtype = 408;
    if(nodes == 8) {
      if(oldtopology[3] == oldtopology[6]) 
  elementtype = 306;
      else
  elementtype = 408;
    }
    if(nodes == 4) elementtype = 404;
    if(nodes == 10) elementtype = 310;
    if(nodes == 6) elementtype = 306;
    if(nodes == 3) elementtype = 303;
  }
  else if(dim == 1) {
    if(nodes == 4) elementtype = 204;
    if(nodes == 3) elementtype = 203;
    if(nodes == 2) elementtype = 202;
  }

  if(!elementtype) {
    printf("Unknown elementtype in element %d (%d nodes, %d dim).\n",
     element,nodes,dim);
  }

  data->elementtypes[element] = elementtype;
  topology = data->topology[element];

  switch (elementtype) {
 
  case 820:
    for(i=0;i<elementtype%100;i++) 
      topology[i] = oldtopology[order820[i]-1];
    break;

  case 504:
    for(i=0;i<elementtype%100;i++) 
      topology[i] = oldtopology[order504[i]-1];
    break;


  case 510:
    for(i=0;i<elementtype%100;i++) {
      if(oldtopology[2] == oldtopology[3]) 
  topology[i] = oldtopology[order510[i]-1];
      else 
  topology[i] = oldtopology[i];
    } 
    break;

  case 605:
    for(i=0;i<elementtype%100;i++) {
      topology[i] = oldtopology[i];
    } 
    break;

  case 613:
    for(i=0;i<elementtype%100;i++) {
      if(oldtopology[4] == oldtopology[5]) 
  topology[i] = oldtopology[order613[i]-1];
      else 
  topology[i] = oldtopology[i];
    } 
    break;

  case 306:
    for(i=0;i<elementtype%100;i++) {
      if(oldtopology[3] == oldtopology[6]) 
  topology[i] = oldtopology[order306[i]-1];    
      else 
  topology[i] = oldtopology[i];
    }
    break;

  default:
    for(i=0;i<elementtype%100;i++) 
      topology[i] = oldtopology[i];

  }  
}


int LoadAnsysInput(struct FemType *data,struct BoundaryType *bound,
          char *prefix,int info)
/* This procedure reads the FEM mesh as written by Ansys. */
{
  int noknots = 0,noelements = 0,nosides,sidetype,currenttype;
  int maxindx,*indx,*revindx,topology[100],ind;
  int i,j,k,l,imax,*nodeindx,*boundindx,boundarynodes;
  int noansystypes,*ansysdim,*ansysnodes,*ansystypes,boundarytypes = 0;
  int namesexist,maxside,sides;
  Real x,y,z = 0;
  FILE *in;
  char *cp,line[MAXLINESIZE],filename[MAXFILESIZE],
    text[MAXNAMESIZE],text2[MAXNAMESIZE];


  /* ExportMesh.header */

  sprintf(filename,"%s.header",prefix);
  if ((in = fopen(filename,"r")) == NULL) {
    printf("LoadAnsysInput: The opening of the header-file %s failed!\n",
     filename);
    return(1);
  }

  if(info) printf("Calculating Ansys elementtypes from %s\n",filename);
  for(i=0;getline;i++);

  noansystypes = i-1;
  printf("There seems to be %d elementytypes in file %s.\n",noansystypes,filename);
  
  ansysdim = Ivector(1,noansystypes);
  ansysnodes = Ivector(1,noansystypes);
  ansystypes = Ivector(1,noansystypes);

  rewind(in);
  for(i=0;i<=noansystypes;i++) {
    Real dummy1,dummy2,dummy3;
    getline;

    /* Ansys writes decimal points also for integers and therefore these 
       values are read in as real numbers. */

    sscanf(line,"%le %le %le",&dummy1,&dummy2,&dummy3);

    if(i==0) {
      noelements = (int) (dummy1+0.5);
      noknots = (int) (dummy2+0.5);
      boundarytypes = (int) (dummy3+0.5);
    }
    else {
      ansysdim[i] = (int) (dummy1+0.5);
      ansysnodes[i] = (int) (dummy2+0.5);
      ansystypes[i] = (int) (dummy3+0.5);
    }
  }
  fclose(in);

  printf("Ansys file has %d elements, %d nodes and %d boundary types.\n",
   noelements,noknots,boundarytypes);
  
 /* ExportMesh.names */

  sprintf(filename,"%s.names",prefix);
  in = fopen(filename,"r");
  if(in == NULL) 
    namesexist = FALSE;
  else
    namesexist = TRUE;

  if(namesexist) printf("Using names of bodies and boundaries\n");


  /* ExportMesh.node */

  sprintf(filename,"%s.node",prefix);
  if ((in = fopen(filename,"r")) == NULL) {
    printf("LoadAnsysInput: The opening of the nodes-file %s failed!\n",
     filename);
    return(2);
  }

  if(info) printf("Calculating Ansys nodes from %s\n",filename);
  for(i=0;getline;i++);

  if(info) printf("There seems to be %d nodes in file %s.\n",i,filename);
  if(i != noknots) printf("Conflicting number of nodes %d vs %d!\n",i,noknots);

  /* Make room and initialize the mesh */
  InitializeKnots(data);

  /* Even 2D elements may form a 3D object! */
  data->dim = 3;
  data->maxnodes = 1;

  for(i=1;i<=noansystypes;i++) {
    if(ansysdim[i] > data->dim) data->dim = ansysdim[i];
    if(ansysnodes[i] > data->maxnodes) data->maxnodes = ansysnodes[i];
  }    
  if(data->maxnodes < 8) data->maxnodes = 8;

  data->noknots = noknots;
  data->noelements = noelements;

  if(info) printf("Allocating for %d nodes and %d elements with max. %d nodes in %d-dim.\n",
      noknots,noelements,data->maxnodes,data->dim);
  AllocateKnots(data);

  indx = Ivector(1,noknots);
  for(i=1;i<=noknots;i++) indx[i] = 0;

  if(info) printf("Loading %d Ansys nodes from %s\n",noknots,filename);
  rewind(in);
  for(i=1;i<=noknots;i++) {
    getline; cp=line;

    indx[i] = next_int(&cp);
    if(cp[0] == '.') cp++;

    x = next_real(&cp);
    if(!cp) x = y = z = 0.0;
    else {
      y = next_real(&cp);
      if(!cp) y = z = 0.0;
      else if(data->dim == 3) {
  z = next_real(&cp);
  if(!cp) z = 0.0;
      }
    }
    data->x[i] = x;
    data->y[i] = y;
    if(data->dim == 3) data->z[i] = z;
  }
  fclose(in);
  
  /* reorder the indexes */
  maxindx = indx[1];
  for(i=1;i<=noknots;i++) 
    if(indx[i] > maxindx) maxindx = indx[i];
  revindx = Ivector(0,maxindx);

  if(maxindx > noknots) {
    printf("There are %d nodes with indexes up to %d.\n",
     noknots,maxindx);
    for(i=1;i<=maxindx;i++) 
      revindx[i] = 0;
    for(i=1;i<=noknots;i++) 
      revindx[indx[i]] = i;
  }
  else {
    for(i=1;i<=noknots;i++) 
      revindx[i] = i;
  }

  /* ExportMesh.elem */

  sprintf(filename,"%s.elem",prefix);
  if ((in = fopen(filename,"r")) == NULL) {
    printf("LoadAnsysInput: The opening of the element-file %s failed!\n",
     filename);
    return(4);
  }
   
  if(info) printf("Loading %d Ansys elements from %s\n",noelements,filename);

  for(j=1;j<=noelements;j++) {

    getline; 
    cp=line;

    for(i=0;i<8;i++) {
      ind = next_int(&cp);
      if(cp[0] == '.') cp++;
      topology[i] = revindx[ind];
    }

    data->material[j] = next_int(&cp);
    currenttype = next_int(&cp);
    
    for(k=1;k<=noansystypes;k++) 
      if(ansystypes[k] == currenttype) break;
    if(ansystypes[k] != currenttype) k=1;

    if(ansysnodes[k] > 8) {
      getline; 
      cp=line;

      if(ansysnodes[k] == 10 && topology[2] != topology[3])
  imax = 10;
      else 
  imax = 20;

      for(i=8;i<imax;i++) {
  ind = next_int(&cp);
  if(cp[0] == '.') cp++;
  topology[i] = revindx[ind];
      }
    }

    ReorderAnsysNodes(data,&topology[0],j,ansysdim[k],ansysnodes[k]);
  }      
  fclose(in);


  /* ExportMesh.boundary */

  sprintf(filename,"%s.boundary",prefix);
  printf("Calculating nodes in file %s\n",filename);
  if ((in = fopen(filename,"r")) == NULL) {
    printf("LoadAnsysInput: The opening of the boundary-file %s failed!\n",
     filename);
    return(5);
  }

  j = 0;
  for(i=0;getline;i++)
    if(!strstr(line,"Boundary")) j++;

  boundarynodes = j;
  nosides = 2*boundarynodes;

  if(info) printf("There are %d boundary nodes, allocating %d elements\n",
      boundarynodes,nosides);
  AllocateBoundary(bound,nosides);
  nodeindx = Ivector(1,boundarynodes);
  boundindx = Ivector(1,boundarynodes);

  if(info) printf("Loading Ansys boundary from %s\n",filename);

  for(i=1;i<=boundarynodes;i++) nodeindx[i] = boundindx[i] = 0;

  rewind(in);
  maxside = 0;
  j = 0;
  for(i=0;getline;i++) {
    if(strstr(line,"Boundary")) {
      sscanf(line,"%d",&sidetype);
      maxside = MAX(sidetype,maxside);
    }
    else {
      j++;
      sscanf(line,"%d",&k);
      nodeindx[j] = revindx[k];
      if(!nodeindx[j]) printf("The boundary %dth node %d is not in index list\n",j,k);
      boundindx[j] = sidetype;
    }
  }
  printf("Found %d boundary nodes with %d as maximum side.\n",j,maxside);
  fclose(in);

  FindPointParents(data,bound,boundarynodes,nodeindx,boundindx,info);

  if(namesexist) {
    int bcind,*bctypes = NULL,*bctypeused = NULL,*bcused = NULL,newsides = 0;

    data->bodynamesexist = TRUE;
    if(bound[0].nosides) {
      newsides = 0;
      bctypes = Ivector(1,maxside);
      bctypeused = Ivector(1,maxside);
      bcused = Ivector(1,bound[0].nosides);      
      for(i=1;i<=bound[0].nosides;i++) bcused[i] = FALSE;      
      data->boundarynamesexist = TRUE;
      for(i=1;i<=maxside;i++) 
  bctypeused[i] = FALSE;      
    }

    sprintf(filename,"%s.names",prefix);
    in = fopen(filename,"r");
            
    for(;;) {
      if(Getrow(line,in,TRUE)) break;
      sscanf(line,"%d%s%s%d",&bcind,&text[0],&text2[0],&sides);
      
      if(strstr(text2,"BODY")) {      
  getline;
  sscanf(line,"%d%d",&j,&bcind);
  strcpy(data->bodyname[bcind],text); 
      }
      else if(strstr(text2,"BOUNDARY")) {      
  /* Read the boundary groups belonging to a particular name */
  for(i=1;i<=maxside;i++) 
    bctypes[i] = 0;
  for(i=1;i<=sides;i++) {
    getline;
    sscanf(line,"%d%d",&j,&bcind);
    bctypes[bcind] = TRUE;
  }
  
  /* Find 1st unsed boundarytype */
  for(i=1;i<=maxside;i++) 
    if(bctypes[i] && !bctypeused[i]) break;
  
  bcind = i;
  bctypeused[bcind] = TRUE;
  if(0) printf("First unused boundary is of type %d\n",bcind);
  strcpy(data->boundaryname[bcind],text);
  
  /* Check which of the BCs have already been named */
  k = l = 0;
  for(i=1;i<=bound[0].nosides;i++) {
    j = bound[0].types[i];
    
    /* The bc is not given any name, hence it can't be a duplicate */
    if(!bctypes[j]) continue;
    
    if(!bcused[i]) {
      k++;
      bcused[i] = bcind;
    }
    else {
      l++;
      if(newsides == 0) AllocateBoundary(&bound[1],bound[0].nosides);
      newsides++;
      bound[1].types[newsides] = bcind;
      bound[1].parent[newsides] =  bound[0].parent[i];
      bound[1].side[newsides] = bound[0].side[i];
      bound[1].parent2[newsides] = bound[0].parent2[i];
      bound[1].side2[newsides] = bound[0].side2[i];
      bound[1].normal[newsides] = bound[0].normal[i];
    }
  }
  if(info) printf("There are %d boundary elements with name %s.\n",k+l,data->boundaryname[bcind]);
      }
    }

    fclose(in);

    /* Put the indexes of all conditions with the same name to be same */
    if(bound[0].nosides) {
      for(i=1;i<=bound[0].nosides;i++) 
  if(bcused[i]) bound[0].types[i] = bcused[i];
      if(newsides) {
  bound[1].nosides = newsides;
  if(info) printf("Created %d additional boundary elements to achieve unique naming.\n",newsides);
      }
      free_Ivector(bctypes,1,maxside);    
      free_Ivector(bctypeused,1,maxside);    
      free_Ivector(bcused,1,bound[0].nosides);
    }
  }

  free_Ivector(boundindx,1,boundarynodes);
  free_Ivector(nodeindx,1,boundarynodes);

  return(0);
}




static void ReorderFieldviewNodes(struct FemType *data,int *oldtopology,
          int element,int dim,int nodes) 
{
  int i,*topology,elementtype = 0;
  int order808[]={1,2,4,3,5,6,8,7};
  int order706[]={1,4,6,2,3,5};
  int order404[]={1,2,3,4};
    
  if(dim == 3) {
    if(nodes == 8) elementtype = 808;
    if(nodes == 6) elementtype = 706;
    if(nodes == 4) elementtype = 504;
  }
  else if(dim == 2) {
    if(nodes == 4) elementtype = 404;
    if(nodes == 3) elementtype = 303;
  }

  if(!elementtype) {
    printf("Unknown elementtype in element %d (%d nodes, %d dim).\n",
     element,nodes,dim);
  }

  data->elementtypes[element] = elementtype;
  topology = data->topology[element];

  for(i=0;i<elementtype%100;i++) {
    if(elementtype == 808) topology[i] = oldtopology[order808[i]-1];
    else if(elementtype == 706) topology[i] = oldtopology[order706[i]-1];
    else if(elementtype == 404) topology[i] = oldtopology[order404[i]-1];
    else topology[i] = oldtopology[i];
  }
}




int LoadFieldviewInput(struct FemType *data,struct BoundaryType *bound,char *prefix,int info)
/* Load the grid from a format that can be read by FieldView
   program by PointWise. This is a suitable format to read files created
   by GridGen. */
{
  int noknots,noelements,elemcode,maxnodes;
  int mode,totelems,entity;
  char filename[MAXFILESIZE];
  char line[MAXLINESIZE],*cp;
  int i,j,k;
  FILE *in;
  Real x,y,z;
  int maxindx,sidenodes;
  char *isio;
  int nobound,nobulk = 0,maxsidenodes,*boundtypes = NULL,**boundtopos = NULL,*boundnodes = NULL,*origtopology;

  if ((in = fopen(prefix,"r")) == NULL) {
    AddExtension(prefix,filename,"dat");
    if ((in = fopen(filename,"r")) == NULL) {
      printf("LoadFieldviewInput: opening of the Fieldview-file '%s' wasn't succesfull !\n",
       filename);
      return(1);
    }
  }
 
  InitializeKnots(data);
  data->dim = 3;
  data->created = TRUE;

  entity = 0;
  mode = 0;
  noknots = 0;
  noelements = 0;
  elemcode = 0;

  maxnodes = 8;
  maxsidenodes = 4;
  maxindx = 0;
  sidenodes = 0;

  data->maxnodes = maxnodes;

  totelems = 0;

  for(;;) {
    
    if(mode == 0) {
      isio = getline;
      
      if(!isio) goto end;
      if(!line) goto end;
      if(line=="") goto end;
      if(strstr(line,"END")) goto end;
      
      /* Control information */
      if(strstr(line,"FIELDVIEW")) mode = 1;
      else if(strstr(line,"CONSTANTS")) mode = 2;
      else if(strstr(line,"GRIDS")) mode = 3;
      else if(strstr(line,"Boundary Table")) mode = 4;
      else if(strstr(line,"Variable Names")) mode = 5;
      else if(strstr(line,"Nodes")) mode = 6;
      else if(strstr(line,"Boundary Faces")) mode = 7;
      else if(strstr(line,"Elements")) mode = 8;
      else if(strstr(line,"Variables")) mode = 9;
      else if(0) printf("unknown: %s",line);
    }      

    switch (mode) {

    case 1: 
      printf("This is indeed a Fieldview input file.\n");
      mode = 0;
      break;
      
    case 2:
      if(0) printf("Constants block\n");
      mode = 0;
      break;

    case 3:
      if(0) printf("Grids block\n");
      mode = 0;
      break;

    case 4:
      if(0) printf("Boundary Table\n");
      mode = 0;
      break;

    case 5:
      if(0) printf("Variable names\n");
      mode = 0;
      break;

    case 6:

      getline;   
      sscanf(line,"%d",&noknots);
      data->noknots = noknots;

      if(info) printf("Loading %d node coordinates\n",noknots);
      
      data->x = Rvector(1,noknots);
      data->y = Rvector(1,noknots);
      data->z = Rvector(1,noknots);
      
      for(i=1;i<=noknots;i++) {
  getline;
  sscanf(line,"%le%le%le",&x,&y,&z);
  data->x[i] = x;
  data->y[i] = y;
  data->z[i] = z;
      }
      mode = 0;
      break;
      
    case 7:
      
      getline;   
      sscanf(line,"%d",&nobound);

      if(info) printf("Loading %d boundary element definitions\n",nobound);

      boundtypes = Ivector(1,nobound);
      boundtopos = Imatrix(1,nobound,0,maxsidenodes-1);
      boundnodes = Ivector(1,nobound);
      
      for(i=1;i<=nobound;i++) {
  getline; cp=line;
  
  boundtypes[i]= next_int(&cp);
  maxsidenodes = next_int(&cp);

  for(j=0;j<maxsidenodes && cp;j++) 
    boundtopos[i][j] = next_int(&cp);

  boundnodes[i] = j;
      }
      mode = 0;
      break;
      
      
    case 8: 

      if(info) printf("Loading bulk element definitions\n");

      if(maxsidenodes == 4) noelements = noknots + nobound;
      else noelements = 6*noknots + nobound;

      origtopology = Ivector(0,maxnodes-1);
      data->topology = Imatrix(1,noelements,0,maxnodes-1);
      data->material = Ivector(1,noelements);
      data->elementtypes = Ivector(1,noelements);
      
      for(i=0;;) {
  getline; cp=line;

  if(strstr(line,"Variables")) mode = 9;
  if(mode != 8) break;
  
  i++;

  k = next_int(&cp);

  if(k == 2) maxnodes = 8;
  else if(k == 1) maxnodes = 4;

  data->material[i]= next_int(&cp);

  for(j=0;j<maxnodes && cp;j++) {
    origtopology[j] = next_int(&cp);
    if(maxindx < origtopology[j]) maxindx = origtopology[j];
  } 

  ReorderFieldviewNodes(data,origtopology,i,3,j);

  if(nobulk+nobound == noelements+1) 
    printf("Too few elements (%d) were allocated!!\n",noelements);

      }
      nobulk = i;

      printf("Found %d bulk elements\n",nobulk);
      if(nobulk+nobound > noelements) printf("Too few elements (%d) were allocated!!\n",noelements);
      printf("Allocated %.4lg %% too many elements\n",
       noelements*100.0/(nobulk+nobound)-100.0);


      for(i=1;i<=nobound;i++) {
  ReorderFieldviewNodes(data,boundtopos[i],i+nobulk,2,boundnodes[i]);
  data->material[i+nobulk] = boundtypes[i];
      }

      data->noelements = noelements = nobulk + nobound;

      mode = 0;
      break;


    case 9: 

      printf("Variables\n");
      
      mode = 0;
      break;
      
    default:      
      break;
    }
  }    
  
end:

  if(maxindx != noknots) 
    printf("The maximum index %d differs from the number of nodes %d !\n",maxindx,noknots);

  ElementsToBoundaryConditions(data,bound,FALSE,TRUE);
  
  return(0);
}




int LoadTriangleInput(struct FemType *data,struct BoundaryType *bound,
          char *prefix,int info)
/* This procedure reads the mesh assuming Triangle format
   */
{
  int noknots,noelements,maxnodes,elematts,nodeatts,dim;
  int elementtype,bcmarkers,sideelemtype;
  int i,j,k,*boundnodes;
  FILE *in;
  char *cp,line[MAXLINESIZE],elemfile[MAXFILESIZE],nodefile[MAXFILESIZE], 
    polyfile[MAXLINESIZE];


  if(info) printf("Loading mesh in Triangle format from file %s.*\n",prefix);

  sprintf(nodefile,"%s.node",prefix);
  if ((in = fopen(nodefile,"r")) == NULL) {
    printf("LoadElmerInput: The opening of the nodes file %s failed!\n",nodefile);
    return(1);
  }
  else 
    printf("Loading nodes from file %s\n",nodefile);

  getline;
  sscanf(line,"%d %d %d %d",&noknots,&dim,&nodeatts,&bcmarkers);
  fclose(in);

  if(dim != 2) {
    printf("LoadTriangleInput assumes that the space dimension is 2, not %d.\n",dim);
    return(2);
  }

  sprintf(elemfile,"%s.ele",prefix);
  if ((in = fopen(elemfile,"r")) == NULL) {
    printf("LoadElmerInput: The opening of the element file %s failed!\n",elemfile);
    return(3);
  }
  else 
    printf("Loading elements from file %s\n",elemfile);

  getline;
  sscanf(line,"%d %d %d",&noelements,&maxnodes,&elematts);
  fclose(in);


  InitializeKnots(data);
  data->dim = dim;
  data->maxnodes = maxnodes;
  data->noelements = noelements;
  data->noknots = noknots;
  elementtype = 300 + maxnodes;

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

  boundnodes = Ivector(1,noknots);
  for(i=1;i<=noknots;i++) 
    boundnodes[i] = 0;

  in = fopen(nodefile,"r");
  getline;
  for(i=1; i <= noknots; i++) {
    getline;
    cp = line;
    j = next_int(&cp);
    if(j != i) printf("LoadTriangleInput: nodes i=%d j=%d\n",i,j);
    data->x[i] = next_real(&cp);
    data->y[i] = next_real(&cp);
    for(j=0;j<nodeatts;j++) {
      next_real(&cp);
    }
    if(bcmarkers > 0) 
      boundnodes[i] = next_int(&cp);
  }
  fclose(in);

  in = fopen(elemfile,"r");
  getline;
  for(i=1; i <= noelements; i++) {
    getline;
    cp = line;
    data->elementtypes[i] = elementtype;
    j = next_int(&cp);
    if(j != i) printf("LoadTriangleInput: elem i=%d j=%d\n",i,j);
    for(j=0;j<3;j++)
      data->topology[i][j] = next_int(&cp);
    if(maxnodes == 6) {
      data->topology[i][4] = next_int(&cp);
      data->topology[i][5] = next_int(&cp);
      data->topology[i][3] = next_int(&cp);
    }
    data->material[i] = 1;
  }
  fclose(in);


  sprintf(polyfile,"%s.poly",prefix);
  if ((in = fopen(polyfile,"r")) == NULL) {
    printf("LoadElmerInput: The opening of the poly file %s failed!\n",polyfile);
    return(1);
  }
  else 
    printf("Loading nodes from file %s\n",polyfile);

  {
    int bcelems,markers,ind1,ind2,bctype,j2,k2,hit = 0;
    int elemsides,sideind[2],side,elemind = 0;

    bctype = 1;
    elemsides = 3;

    getline;
    getline;
    sscanf(line,"%d %d",&bcelems,&markers);

    CreateInverseTopology(data,info);

    AllocateBoundary(bound,bcelems);

    for(i=1;i<=bcelems;i++) {
      
      getline;
      if(markers)
  sscanf(line,"%d %d %d %d",&j,&ind1,&ind2,&bctype);
      else 
  sscanf(line,"%d %d %d",&j,&ind1,&ind2);
     
      /* find an element which owns both the nodes */
      for(j=1;j<=data->maxinvtopo;j++) {
  hit = FALSE;
  k = data->invtopo[j][ind1];
  if(!k) break;

  for(j2=1;j2<=data->maxinvtopo;j2++) { 
    k2 = data->invtopo[j2][ind2];
    if(!k2) break;
    if(k == k2) {
      hit = TRUE;
      elemind = k;
      break;
    }
  }
  if(hit) break;
      }
      if(!hit) return(1);


      /* Find the correct side of the triangular element */
      for(side=0;side<elemsides;side++) {
  GetElementSide(elemind,side,1,data,&sideind[0],&sideelemtype);
  
  hit = FALSE;
  if(sideind[0] == ind1 && sideind[1] == ind2) hit = TRUE;
  if(sideind[0] == ind2 && sideind[1] == ind1) hit = TRUE;

  if(hit) {
    bound->parent[i] = elemind;
    bound->side[i] = side;
    bound->parent2[i] = 0;
    bound->side2[i] = 0;
    bound->types[i] = bctype;
  }
      }
    }
  } 

  printf("Successfully read the mesh from the Triangle input file.\n");

  return(0);
}




int LoadMeditInput(struct FemType *data,struct BoundaryType *bound,
       char *prefix,int info)
/* This procedure reads the mesh assuming Medit format
   */
{
  int noknots = 0,noelements = 0,maxnodes,dim = 0,elementtype;
  int i,j,allocated;
  FILE *in;
  char *cp,line[MAXLINESIZE],nodefile[MAXFILESIZE];


  sprintf(nodefile,"%s.mesh",prefix);
  if(info) printf("Loading mesh in Medit format from file %s\n",prefix);

  if ((in = fopen(nodefile,"r")) == NULL) {
    printf("LoadElmerInput: The opening of the mesh file %s failed!\n",nodefile);
    return(1);
  }

  allocated = FALSE;
  maxnodes = 0;
 
allocate:

  if(allocated) {
    InitializeKnots(data);
    data->dim = dim;
    data->maxnodes = maxnodes;
    data->noelements = noelements;
    data->noknots = noknots;
    elementtype = 300 + maxnodes;

    if(info) printf("Allocating for %d knots and %d elements.\n",noknots,noelements);
    AllocateKnots(data);
    in = fopen(nodefile,"r");
  }


  for(;;) {
    if(Getrow(line,in,TRUE)) goto end;
    if(!line) goto end;
    if(strstr(line,"END")) goto end;
    
    if(strstr(line,"DIMENSION")) {
      if(Getrow(line,in,TRUE)) goto end;
      cp = line;
      dim = next_int(&cp);
      printf("dim = %d %s",dim,line);
    }
    else if(strstr(line,"VERTICES")) {
      printf("verts: %s",line);

      if(Getrow(line,in,TRUE)) goto end;
      cp = line;
      noknots = next_int(&cp);      

      printf("noknots = %d %s",noknots,line);

      for(i=1; i <= noknots; i++) {
  getline;
#if 0
  printf("i=%d line=%s",i,line);
#endif
  if(allocated) {
    cp = line;
#if 0
    printf("cp = %s",cp);
#endif

    data->x[i] = next_real(&cp);
    data->y[i] = next_real(&cp);
    if(dim > 2) data->z[i] = next_real(&cp);
  }
      }
    }
    else if(strstr(line,"TRIANGLES")) {
      if(Getrow(line,in,TRUE)) goto end;
      cp = line;
      noelements = next_int(&cp);      

      printf("noelements = %d %s",noelements,line);

      elementtype = 303;
      if(maxnodes < 3) maxnodes = 3;

      for(i=1; i <= noelements; i++) {
  getline;
  if(allocated) {
    cp = line;
    data->elementtypes[i] = elementtype;
    for(j=0;j<3;j++)
      data->topology[i][j] = next_int(&cp);
    data->material[i] = next_int(&cp);
  }
      }
    }
#if 0
    else printf("unknown command: %s",line);
#endif
  }    

end:
  fclose(in);

printf("ALLOCATED=%d\n",allocated);

  if(!allocated) {
    allocated = TRUE;
    goto allocate;
  }

  printf("Successfully read the mesh from the Medit input file.\n");

  return(0);
}




int LoadGidInput(struct FemType *data,struct BoundaryType *bound,
     char *prefix,int info)
/* Load the grid from GID mesh format */
{
  int noknots,noelements,elemcode,maxnodes,material,foundsame;
  int mode,allocated,maxknot,nosides,sideelemtype;
  int boundarytype,materialtype,boundarynodes,side,parent,elemsides;
  int dim = 0, elemnodes = 0, elembasis = 0, elemtype = 0, bulkdone, usedmax = 0,hits;
  int minbulk,maxbulk,minbound,maxbound,label,debug;
  int *usedno = NULL, **usedelem = NULL;  
  char filename[MAXFILESIZE],line[MAXLINESIZE],*cp;
  int i,j,k,n,ind,inds[MAXNODESD2],sideind[MAXNODESD1];
  FILE *in;
  Real x,y,z = 0;

  debug = FALSE;

  strcpy(filename,prefix);
  if ((in = fopen(filename,"r")) == NULL) {
    AddExtension(prefix,filename,"msh");
    if ((in = fopen(filename,"r")) == NULL) {
      printf("LoadAbaqusInput: opening of the GID-file '%s' wasn't succesfull !\n",
       filename);
      return(1);
    }
  }

  printf("Reading mesh from GID file %s.\n",filename);
  InitializeKnots(data);

  allocated = FALSE;
  maxknot = 0;

  /* Because the file format doesn't provide the number of elements
     or nodes the results are read twice but registered only in the
     second time. */

  minbulk = minbound = 1000;
  maxbulk = maxbound = 0;

omstart:

  mode = 0;
  maxnodes = 0;
  noknots = 0;
  noelements = 0;
  elemcode = 0;
  boundarytype = 0;
  boundarynodes = 0;
  material = 0;
  nosides = 0;
  bulkdone = FALSE;

  for(;;) {

    if(Getrow(line,in,FALSE)) goto end;
    if(!line) goto end;

    if(strstr(line,"MESH")) {
      if(debug) printf("MESH\n");

      if(strstr(line,"dimension 3")) 
  dim = 3;
      else if(strstr(line,"dimension 2")) 
  dim = 2;
      else printf("Unknown dimension\n");
      
      if(strstr(line,"ElemType Tetrahedra"))
  elembasis = 500;
      else if(strstr(line,"ElemType Triangle"))
  elembasis = 300;
      else if(strstr(line,"ElemType Linear"))
  elembasis = 200;
      else printf("Unknown elementtype: %s\n",line);

      if(strstr(line,"Nnode 4"))
  elemnodes = 4;
      else if(strstr(line,"Nnode 3"))
  elemnodes = 3;
      else if(strstr(line,"Nnode 2"))
  elemnodes = 2;
      else printf("Unknown elementnode: %s\n",line);

      if(elemnodes > maxnodes) maxnodes = elemnodes;      
      elemtype = elembasis + elemnodes;
      mode = 0;
      
      if(debug) printf("dim=%d elemtype=%d\n",dim,elemtype);
    }
    else if(strstr(line,"Coordinates")) {
      if(debug) printf("Start coords\n");
      mode = 1;
    }
    else if(strstr(line,"end coordinates")) {
      if(debug) printf("End coords\n");
      mode = 0;
    }
    else if(strstr(line,"Elements")) {
      if(bulkdone) { 
  if(debug) printf("Start boundary elems\n");
  mode = 3;
  boundarytype++;
      }
      else {
  if(debug) printf("Start bulk elems\n"); 
  mode = 2;
      }
    }
    else if(strstr(line,"end elements")) {
      if(debug) printf("End elems\n");

      mode = 0;
      if(!bulkdone && allocated) {
  usedno = Ivector(1,data->noknots);

  for(j=1;j<=data->noknots;j++)
    usedno[j] = 0;

  for(j=1;j<=data->noelements;j++) {
    n = data->elementtypes[j] % 100; 
    for(i=0;i<n;i++) {
      ind = data->topology[j][i];
      usedno[data->topology[j][i]] += 1;
    }
  }

  usedmax = 0;
  for(i=1;i<=data->noknots;i++)
    if(usedno[i] > usedmax) usedmax = usedno[i];

  for(j=1;j<=data->noknots;j++)
    usedno[j] = 0;

  usedelem = Imatrix(1,data->noknots,1,usedmax);
  for(j=1;j<=data->noknots;j++)
    for(i=1;i<=usedmax;i++)
      usedelem[j][i] = 0;

  for(j=1;j<=data->noelements;j++) {
    n = data->elementtypes[j] % 100;
    for(i=0;i<n;i++) {
      ind = data->topology[j][i];
      usedno[ind] += 1;
      k = usedno[ind];
      usedelem[ind][k] = j;
    }
  }
      }
      bulkdone = TRUE;
    }
    else if(!mode) {
      if(debug) printf("mode: %d %s\n",mode,line);
    }
    else if(mode) {  
      switch (mode) {
  
      case 1:
  cp = line;
  ind = next_int(&cp);
  if(ind > noknots) noknots = ind;

  x = next_real(&cp);
  y = next_real(&cp);
  if(dim == 3) z = next_real(&cp);
  
  if(allocated) {
    data->x[ind] = x;
    data->y[ind] = y;
    if(dim == 3) data->z[ind] = z;
  }
  break;

      case 2:
  cp = line;
  ind = next_int(&cp);
  if(ind > noelements) noelements = ind;
  
  for(i=0;i<elemnodes;i++) {
    k = next_int(&cp);
    if(allocated) {
      data->topology[ind][i] = k;
      data->elementtypes[ind] = elemtype;
      data->material[ind] = 1;
    }
  }
  label = next_int(&cp);
  
  if(allocated) {
    if(label) { 
      materialtype = label-minbulk+1;
    }
    else if(maxbound) {
      materialtype = maxbulk-minbulk + 2;
    }
    else { 
      materialtype = 1;
    }
    data->material[ind] = materialtype;
  }
  else {
    if(label > maxbulk) maxbulk = label;
    if(label < minbulk) minbulk = label;    
  }

  break;
  
      case 3:
  cp = line;
  ind = next_int(&cp);
  nosides++;

  for(i=0;i<elemnodes;i++) {
    k = next_int(&cp);
    inds[i] = k;
  }
  label = next_int(&cp);
  
  if(!allocated) {
    if(label) {
      if(label > maxbound) maxbound = label;
      if(label < minbound) minbound = label;
    }
  }

  if(allocated) {
    
    if(label) { 
      boundarytype = label-minbound+1;
    }
    else if(maxbound) {
      boundarytype = maxbound-minbound + 2;
    }
    else { 
      boundarytype = 1;
    }

    foundsame = FALSE;
    for(i=1;i<=usedno[inds[0]];i++) {
      parent = usedelem[inds[0]][i];

      elemsides = data->elementtypes[parent] % 100;
      
      for(side=0;side<elemsides;side++) {

        GetElementSide(parent,side,1,data,&sideind[0],&sideelemtype);

        if(elemnodes != sideelemtype%100) printf("LoadGidMesh: bug?\n");
        
        hits = 0;
        for(j=0;j<elemnodes;j++) 
    for(k=0;k<elemnodes;k++) 
      if(sideind[j] == inds[k]) hits++;

        if(hits < elemnodes) continue;
        
        if(!foundsame) {
    foundsame++;
    bound->parent[nosides] = parent;
    bound->side[nosides] = side;
    bound->parent2[nosides] = 0;
    bound->side2[nosides] = 0;
    bound->types[nosides] = boundarytype;   
        }
        else if(foundsame == 1) {
    if(parent == bound->parent[nosides]) continue;
    foundsame++;
    bound->parent2[nosides] = parent;
    bound->side2[nosides] = side;     
        }
        else if(foundsame > 1) {
    printf("Boundary %d has more than 2 parents\n",nosides);
        }
      }
    }
    if(!foundsame) {
      printf("Did not find parent for side %d\n",nosides);
      nosides--;
    }
    else {
      if(debug) printf("Parent of side %d is %d\n",nosides,bound->parent[nosides]);
    }
  }
      }
    }
  }

end:

  if(!allocated) {
    rewind(in);
    data->noknots = noknots;
    data->noelements = noelements;
    data->maxnodes = maxnodes;
    data->dim = dim;
    
    if(info) {
      printf("Allocating for %d knots and %d %d-node elements.\n",
       noknots,noelements,maxnodes);
      printf("Initial material indexes are at interval %d to %d.\n",minbulk,maxbulk);
    }  
    AllocateKnots(data);

    if(info) {
      printf("Allocating %d boundary elements\n",nosides);
      printf("Initial boundary indexes are at interval %d to %d.\n",minbound,maxbound);
    }
    AllocateBoundary(bound,nosides);

    bound->nosides = nosides;
    bound->created = TRUE;
    
    nosides = 0;
    bulkdone = FALSE;
    boundarytype = 0;

    allocated = TRUE;
    goto omstart;
  }

  bound->nosides = nosides;
  free_Ivector(usedno,1,data->noknots);
  free_Imatrix(usedelem,1,data->noknots,1,usedmax);

  if(info) printf("The mesh was loaded from file %s.\n",filename);
  return(0);
}



static void ReorderComsolNodes(int elementtype,int *topo) 
{
  int i,tmptopo[MAXNODESD2];
  int order404[]={1,2,4,3};
  int order808[]={1,2,4,3,5,6,8,7};
 
   
  switch (elementtype) {
 
  case 404:
    for(i=0;i<elementtype%100;i++) 
      tmptopo[i] = topo[i];
    for(i=0;i<elementtype%100;i++) 
      topo[i] = tmptopo[order404[i]-1];
    break;

  case 808:
    for(i=0;i<elementtype%100;i++) 
      tmptopo[i] = topo[i];
    for(i=0;i<elementtype%100;i++) 
      topo[i] = tmptopo[order808[i]-1];
    break;

  default:
    break;

  }  
}



int LoadComsolMesh(struct FemType *data,struct BoundaryType *bound,char *prefix,int info)
/* Load the grid in Comsol Multiphysics mesh format */
{
  int noknots,noelements,elemcode,maxnodes,material,allocated;

  int dim = 0, elemnodes = 0, elembasis = 0, elemtype;
  int debug,offset,domains,mindom,minbc,elemdim = 0;
  char filename[MAXFILESIZE],line[MAXLINESIZE],*cp;
  int i,j,k;
  FILE *in;

  strcpy(filename,prefix);
  if ((in = fopen(filename,"r")) == NULL) {
    AddExtension(prefix,filename,"mphtxt");
    if ((in = fopen(filename,"r")) == NULL) {
      printf("LoadComsolMesh: opening of the Comsol mesh file '%s' wasn't succesfull !\n",
       filename);
      return(1);
    }
  }

  printf("Reading mesh from Comsol mesh file %s.\n",filename);
  InitializeKnots(data);

  debug = FALSE;
  allocated = FALSE;

  mindom = 1000;
  minbc = 1000;
  offset = 1;

omstart:

  maxnodes = 0;
  noknots = 0;
  noelements = 0;
  elemcode = 0;
  material = 0;
  domains = 0;


  for(;;) {

    if(Comsolrow(line,in)) goto end;
    if(!line) goto end;

    if(strstr(line,"# sdim")) {
      cp = line;
      dim = next_int(&cp);
      if(debug) printf("dim=%d\n",dim);
    }

    else if(strstr(line,"# number of mesh points")) {
      cp = line;
      noknots = next_int(&cp);
      if(debug) printf("noknots=%d\n",noknots);
    }

    else if(strstr(line,"# lowest mesh point index")) {
      cp = line;
      offset = 1 - next_int(&cp);
      if(debug) printf("offset=%d\n",offset);
    }

    else if(strstr(line,"# type name")) {
      if(strstr(line,"vtx")) elembasis = 100;
      else if(strstr(line,"edg")) elembasis = 200;
      else if(strstr(line,"tri")) elembasis = 300;
      else if(strstr(line,"quad")) elembasis = 400;
      else if(strstr(line,"tet")) elembasis = 500;
      else if(strstr(line,"prism")) elembasis = 700;
      else if(strstr(line,"hex")) elembasis = 800;
      else printf("unknown element type = %s",line);
    }

    else if(strstr(line,"# number of nodes per element")) {
      cp = line;
      elemnodes = next_int(&cp);
      if(elemnodes > maxnodes) maxnodes = elemnodes;      
      if(debug) printf("elemnodes=%d\n",elemnodes);           
    }

    else if(strstr(line,"# Mesh point coordinates")) {
      printf("Loading %d coordinates\n",noknots);

      for(i=1;i<=noknots;i++) {
  Comsolrow(line,in); 

  if(allocated) {
    cp = line;
    data->x[i] = next_real(&cp);
    data->y[i] = next_real(&cp);
    if(dim == 3) data->z[i] = next_real(&cp);
  }
      }
    }

    else if(strstr(line,"# number of elements")) {

      cp = line;
      k = next_int(&cp);

      Comsolrow(line,in);             
      elemtype = elemnodes + elembasis;
      elemdim = GetElementDimension(elemtype);

      if(debug) printf("Loading %d elements of type %d\n",k,elemtype);
      domains = noelements;

      for(i=1;i<=k;i++) {
  Comsolrow(line,in); 

  if(dim == 3 && elembasis < 300) continue;
  if(dim == 2 && elembasis < 200) continue;

  noelements = noelements + 1;
  if(allocated) {
    data->elementtypes[noelements] = elemtype;
    data->material[noelements] = 1;   
    cp = line;
    for(j=0;j<elemnodes;j++)
      data->topology[noelements][j] = next_int(&cp) + offset;

    if(1) ReorderComsolNodes(elemtype,data->topology[noelements]);

  }
      }
    }

    else if(strstr(line,"# number of geometric entity indices") || 
      strstr(line,"# number of domains")) {

      cp = line;
      k = next_int(&cp);

      Comsolrow(line,in);             
      if(debug) printf("Loading %d domains for the elements\n",k);

      for(i=1;i<=k;i++) {
  Comsolrow(line,in); 

  if(dim == 3 && elembasis < 300) continue;
  if(dim == 2 && elembasis < 200) continue;

  domains = domains + 1;
  cp = line;
  material = next_int(&cp);

  if(allocated) {
    if(elemdim < dim) 
      material = material - minbc + 1;
    else 
      material = material - mindom + 1;
    data->material[domains] = material;   
  }
  else {
    if(elemdim < dim) {
      if(minbc > material) minbc = material;
    }
    else { 
      if(mindom > material) mindom = material;    
    }
  }

      }
    }

    else if(strstr(line,"#")) {
      if(debug) printf("Unused command:  %s",line);
    }

  }

end:

  if(!allocated) {

    if(noknots == 0 || noelements == 0 || maxnodes == 0) {
       printf("Invalid mesh consits of %d knots and %d %d-node elements.\n",
       noknots,noelements,maxnodes);     
       fclose(in);
       return(2);
    }

    rewind(in);
    data->noknots = noknots;
    data->noelements = noelements;
    data->maxnodes = maxnodes;
    data->dim = dim;
    
    if(info) {
      printf("Allocating for %d knots and %d %d-node elements.\n",
       noknots,noelements,maxnodes);
    }  
    AllocateKnots(data);
    allocated = TRUE;

    goto omstart;    
  }
  fclose(in);

  if(info) printf("The Comsol mesh was loaded from file %s.\n\n",filename);
  ElementsToBoundaryConditions(data,bound,FALSE,TRUE);

  return(0);
}


static int GmshToElmerType(int gmshtype)
{
  int elmertype = 0;

  switch (gmshtype) {
      
  case 1:        
    elmertype = 202;
    break;
  case 2:        
    elmertype = 303;
    break;
  case 3:        
    elmertype = 404;
    break;
  case 4:        
    elmertype = 504;
    break;
  case 5:        
    elmertype = 808;
    break;
  case 6:        
    elmertype = 706;
    break;
  case 7:        
    elmertype = 605;
    break;
  case 8:        
    elmertype = 203;
    break;
  case 9:        
    elmertype = 306;
    break;
  case 10:        
    elmertype = 409;
    break;
  case 11:        
    elmertype = 510;
    break;
  case 12:        
    elmertype = 827;
    break;
  case 15:        
    elmertype = 101;
    break;

  case 16:
    elmertype = 408;
    break;
  case 17:
    elmertype = 820;
    break;
  case 18:
    elmertype = 715;
    break;
  case 19:
    elmertype = 613;
    break;
  case 21:
    elmertype = 310;
    break;

  default:
    printf("Gmsh element %d does not have an Elmer counterpart!\n",gmshtype);
  }

  return(elmertype);
}


static void GmshToElmerIndx(int elemtype,int *topology)
{
  int i=0,nodes=0,oldtopology[MAXNODESD2];
  int reorder, *porder;

  int order510[]={0,1,2,3,4,5,6,7,9,8};
  int order613[]={0,1,2,3,4,5,8,10,6,7,9,11,12};
  int order715[]={0,1,2,3,4,5,6,9,7,8,10,11,12,14,13};
  int order820[]={0,1,2,3,4,5,6,7,8,11,12,9,10,12,14,15,16,18,19,17};


  reorder = FALSE;

  switch (elemtype) {
      
  case 510:        
    reorder = TRUE;
    porder = &order510[0];
    break;

  case 613:        
    reorder = TRUE;
    porder = &order613[0];
    break;

  case 715:        
    reorder = TRUE;
    porder = &order715[0];
    break;

  case 820:        
    reorder = TRUE;
    porder = &order820[0];
    break;

  }

  if( reorder ) {
    nodes = elemtype % 100;
    for(i=0;i<nodes;i++) 
      oldtopology[i] = topology[i];
    for(i=0;i<nodes;i++) 
      topology[i] = oldtopology[porder[i]];
  }
}


static int LoadGmshInput1(struct FemType *data,struct BoundaryType *bound,
        char *filename,int info)
{
  int noknots = 0,noelements = 0,maxnodes,elematts,nodeatts,nosides,dim;
  int sideind[MAXNODESD1],elemind[MAXNODESD2],tottypes,elementtype,bcmarkers;
  int i,j,k,dummyint,*boundnodes,allocated,*revindx,maxindx;
  int elemno, gmshtype, regphys, regelem, elemnodes,maxelemtype,elemdim;
  FILE *in;
  char *cp,line[MAXLINESIZE];


  if ((in = fopen(filename,"r")) == NULL) {
    printf("LoadElmerInput: The opening of the mesh file %s failed!\n",filename);
    return(1);
  }
  if(info) printf("Loading mesh in Gmsh format 1.0 from file %s\n",filename);

  allocated = FALSE;
  dim = 3;
  maxnodes = 0;
  maxindx = 0;
  maxelemtype = 0;

allocate:

  if(allocated) {
    InitializeKnots(data);
    data->dim = dim;
    data->maxnodes = maxnodes;
    data->noelements = noelements;
    data->noknots = noknots;

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

    if(maxindx > noknots) {
      revindx = Ivector(1,maxindx);
      for(i=1;i<=maxindx;i++) revindx[i] = 0;
    }
    in = fopen(filename,"r");
  }


  for(;;) {
    if(Getrow(line,in,TRUE)) goto end;
    if(!line) goto end;
    if(strstr(line,"END")) goto end;

    if(strstr(line,"$NOD")) {
      
      getline;
      cp = line;
      noknots = next_int(&cp);

      for(i=1; i <= noknots; i++) {
  getline;
  cp = line;

  j = next_int(&cp);
  if(allocated) {
    if(maxindx > noknots) revindx[j] = i;
    data->x[i] = next_real(&cp);
    data->y[i] = next_real(&cp);
    if(dim > 2) data->z[i] = next_real(&cp);
  }
  else {
    maxindx = MAX(j,maxindx);
  }
      }
      getline;
      if(!strstr(line,"$ENDNOD")) {
  printf("NOD section should end to string ENDNOD\n");
  printf("%s\n",line);
      }
    }
    
    if(strstr(line,"$ELM")) {
      getline;
      cp = line;
      noelements = next_int(&cp);

      for(i=1; i <= noelements; i++) {
  getline;
  
  cp = line;
  elemno = next_int(&cp);
  gmshtype = next_int(&cp);
  regphys = next_int(&cp);
  regelem = next_int(&cp);
  elemnodes = next_int(&cp);

  if(allocated) {
    elementtype = GmshToElmerType(gmshtype);
    maxelemtype = MAX(maxelemtype,elementtype);
    data->elementtypes[i] = elementtype;
    data->material[i] = regphys;

    if(elemnodes != elementtype%100) {
      printf("Conflict in elementtypes %d and number of nodes %d!\n",
       elementtype,elemnodes);
    }   
    for(j=0;j<elemnodes;j++)
      elemind[j] = next_int(&cp);
    
    GmshToElmerIndx(elementtype,elemind);   

    for(j=0;j<elemnodes;j++)
      data->topology[i][j] = elemind[j];
  }
  else {
    maxnodes = MAX(elemnodes,maxnodes);
  }

      }
      getline;
      if(!strstr(line,"$ENDELM")) 
  printf("ELM section should end to string ENDELM\n");
    }

  }

 end:

  fclose(in);

  if(!allocated) {
    allocated = TRUE;
    goto allocate;
  }


  if(maxindx > noknots) {
    printf("Renumbering the Gmsh nodes from %d to %d\n",maxindx,noknots);

    for(i=1; i <= noelements; i++) {
      elementtype = data->elementtypes[i];
      elemnodes = elementtype % 100; 

      for(j=0;j<elemnodes;j++) {
  k = data->topology[i][j];
  if(k <= 0 || k > maxindx) 
    printf("index out of bounds %d\n",k);
  else if(revindx[k] <= 0) 
    printf("unkonwn node %d %d in element %d\n",k,revindx[k],i);
  else 
    data->topology[i][j] = revindx[k];
      }      
    }
    free_Ivector(revindx,1,maxindx);
  }

  ElementsToBoundaryConditions(data,bound,FALSE,info);

  printf("Successfully read the mesh from the Gmsh input file.\n");

  return(0);
}


static int LoadGmshInput2(struct FemType *data,struct BoundaryType *bound,
        char *filename,int info)
{
  int noknots = 0,noelements = 0,maxnodes,elematts,nodeatts,nosides,dim,notags;
  int sideind[MAXNODESD1],elemind[MAXNODESD2],tottypes,elementtype,bcmarkers;
  int i,j,k,dummyint,*boundnodes,allocated,*revindx,maxindx;
  int elemno, gmshtype, tagphys, taggeom, tagpart, elemnodes,maxelemtype,elemdim;
  int usetaggeom,tagmat,verno;
  FILE *in;
  char *cp,line[MAXLINESIZE];


  if ((in = fopen(filename,"r")) == NULL) {
    printf("LoadElmerInput: The opening of the mesh file %s failed!\n",filename);
    return(1);
  }
  if(info) printf("Loading mesh in Gmsh format 2.0 from file %s\n",filename);

  allocated = FALSE;
  dim = 3;
  maxnodes = 0;
  maxindx = 0;
  maxelemtype = 0;
  usetaggeom = FALSE;

omstart:

  for(;;) {
    if(Getrow(line,in,FALSE)) goto end;
    if(!line) goto end;
    if(strstr(line,"$End")) continue;

    if(strstr(line,"$MeshFormat")) {
      getline;
      cp = line;
      verno = next_int(&cp);

      if(verno != 2) {
  printf("Version number is not compatible with the parser: %d\n",verno);
      }

      getline;
      if(!strstr(line,"$EndMeshFormat")) {
  printf("$MeshFormat section should end to string $EndMeshFormat:\n%s\n",line);
      }      
    }
      
    else if(strstr(line,"$Nodes")) {
      getline;
      cp = line;
      noknots = next_int(&cp);

      for(i=1; i <= noknots; i++) {
  getline;
  cp = line;

  j = next_int(&cp);
  if(allocated) {
    if(maxindx > noknots) revindx[j] = i;
    data->x[i] = next_real(&cp);
    data->y[i] = next_real(&cp);
    if(dim > 2) data->z[i] = next_real(&cp);
  }
  else {
    maxindx = MAX(j,maxindx);
  }
      }
      getline;
      if(!strstr(line,"$EndNodes")) {
  printf("$Nodes section should end to string $EndNodes:\n%s\n",line);
      }           
    }
    
    else if(strstr(line,"$Elements")) {
      getline;
      cp = line;
      noelements = next_int(&cp);

      for(i=1; i <= noelements; i++) {
  getline;
  
  cp = line;
  elemno = next_int(&cp);
  gmshtype = next_int(&cp);
  elementtype = GmshToElmerType(gmshtype);

  if(allocated) {
    elemnodes = elementtype % 100;
    data->elementtypes[i] = elementtype;

    /* Point does not seem to have physical properties */
    notags = next_int(&cp);
    if(notags > 0) tagphys = next_int(&cp);
    if(notags > 1) taggeom = next_int(&cp);
    if(notags > 2) tagpart = next_int(&cp);
    for(j=4;j<=notags;j++)
      next_int(&cp);

    if(tagphys) {
      tagmat = tagphys;
    }
    else {
      tagmat = taggeom;
      usetaggeom = TRUE;
    }

    data->material[i] = tagmat;
    for(j=0;j<elemnodes;j++)
      elemind[j] = next_int(&cp);

    GmshToElmerIndx(elementtype,elemind);   

    for(j=0;j<elemnodes;j++)
      data->topology[i][j] = elemind[j];
  }
  else {
    maxelemtype = MAX(maxelemtype,elementtype);
  }
  
      }
      getline;
      if(!strstr(line,"$EndElements")) {
  printf("$Elements section should end to string $EndElements:\n%s\n",line);
      }   
    }
    else if(strstr(line,"$PhysicalNames")) {
      if(info) printf("Physical names are not accounted for\n");
      getline;
      cp = line;
      i = next_int(&cp);
      for(;i>0;i--) getline;

      getline;
      if(!strstr(line,"$EndPhysicalNames")) {
  printf("$PhysicalNames section should end to string $EndPhysicalNames:\n%s\n",line);
      }   
    }
    else {
      if(!allocated) printf("Untreated command: %s",line);
    }

  }

 end:


  if(!allocated) {
    maxnodes = maxelemtype % 100;
    InitializeKnots(data);
    data->dim = dim;
    data->maxnodes = maxnodes;
    data->noelements = noelements;
    data->noknots = noknots;

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

    if(maxindx > noknots) {
      revindx = Ivector(1,maxindx);
      for(i=1;i<=maxindx;i++) revindx[i] = 0;
    }
    rewind(in);
    allocated = TRUE;
    goto omstart;
  }

  if(maxindx > noknots) {
    printf("Renumbering the Gmsh nodes from %d to %d\n",maxindx,noknots);

    for(i=1; i <= noelements; i++) {
      elementtype = data->elementtypes[i];
      elemnodes = elementtype % 100; 

      for(j=0;j<elemnodes;j++) {
  k = data->topology[i][j];
  if(k <= 0 || k > maxindx) 
    printf("index out of bounds %d\n",k);
  else if(revindx[k] <= 0) 
    printf("unkonwn node %d %d in element %d\n",k,revindx[k],i);
  else 
    data->topology[i][j] = revindx[k];
      }      
    }
    free_Ivector(revindx,1,maxindx);
  }

  ElementsToBoundaryConditions(data,bound,FALSE,info);

  /* The geometric entities are rather randomly numbered */
  if( usetaggeom ) {
    RenumberBoundaryTypes(data,bound,TRUE,0,info);
    RenumberMaterialTypes(data,bound,info);
  }

  if(info) printf("Successfully read the mesh from the Gmsh input file.\n");

  return(0);
}



int LoadGmshInput(struct FemType *data,struct BoundaryType *bound,
       char *prefix,int info)
{
  FILE *in;
  char line[MAXLINESIZE],filename[MAXFILESIZE];
  int errno;

  sprintf(filename,"%s",prefix);
  if ((in = fopen(filename,"r")) == NULL) {
    sprintf(filename,"%s.msh",prefix);
    if ((in = fopen(filename,"r")) == NULL) {
      printf("LoadElmerInput: The opening of the mesh file %s failed!\n",filename);
      return(1);
    }
  }

  Getrow(line,in,FALSE);
  fclose(in);

  if(info) {
    printf("Format chosen using the first line: %s",line);
  }

  if(strstr(line,"$MeshFormat")) 
    errno = LoadGmshInput2(data,bound,filename,info);
  else {
    printf("*****************************************************\n");
    printf("The $MeshFormat was not given, assuming Gmsh 1 format\n");
    printf("This version of Gmsh format is no longer supported\n");
    printf("Please use Gsmh 2 version for output\n");
    printf("*****************************************************\n");
    
    errno = LoadGmshInput1(data,bound,filename,info);
  }     

  return(errno);
}




static int UnvToElmerType(int unvtype)
{
  int elmertype;

  switch (unvtype) {

  case 11: 
  case 21:
    elmertype = 202;
    break;

  case 22:
  case 23:
    elmertype = 203;
    break;

  case 41:
  case 51:
  case 61:
  case 74:
  case 81:
  case 91:
    elmertype = 303;
    break;

  case 42:
  case 52:
  case 62:
  case 72:
  case 82:
  case 92:
    elmertype = 306;
    break;

  case 43:
  case 53:
  case 63:
  case 73:
  case 93:
    elmertype = 310;
    break;

  case 44:
  case 54:
  case 64:
  case 71:
  case 84:
  case 94:
    elmertype = 404;
    break;

  case 45:
  case 46:
  case 56:
  case 66:
  case 76:
  case 96:
    elmertype = 408;
    break;

  case 111:
    elmertype = 504;
    break;

  case 118:
    elmertype = 510;
    break;

  case 112:
    elmertype = 706;
    break;

  case 113:
    elmertype = 715;
    break;

  case 115:
    elmertype = 808;
    break;

  case 116:
    elmertype = 820;
    break;

  default:
    elmertype = 0;
    printf("Unknown elementtype in universal mesh format: %d\n",unvtype);
  }

  return(elmertype);
}


static int UnvRedundantIndexes(int nonodes,int *ind)
{
  int i,j,redundant;
  
  redundant = FALSE;
  for(i=0;i<nonodes;i++) {
    if( ind[i] == 0 ) redundant = TRUE;
    for(j=i+1;j<nonodes;j++) 
      if(ind[i] == ind[j]) redundant = TRUE;
  }
  if( redundant ) {
    printf("Redundant element %d: ",nonodes);
    for(i=0;i<nonodes;i++)
      printf(" %d ",ind[i]);
    printf("\n");
  }

  return(redundant);
}


static void UnvToElmerIndx(int elemtype,int *topology)
{
  int i=0,nodes=0,oldtopology[MAXNODESD2];
  int reorder, *porder;

  int order510[]={1,3,5,10,2,4,6,7,8,9};
  int order408[]={1,3,5,7,2,4,6,8};
  int order820[]={1,3,5,7,13,15,17,19,2,4,6,8,9,10,11,12,14,16,18,20};


  reorder = FALSE;

  switch (elemtype) {
      
  case 510:        
    reorder = TRUE;
    porder = &order510[0];
    break;

  case 408:        
    reorder = TRUE;
    porder = &order408[0];
    break;

  case 820:        
    reorder = TRUE;
    porder = &order820[0];
    break;
   
  }

  if( reorder ) {
    nodes = elemtype % 100;
    for(i=0;i<nodes;i++) 
      oldtopology[i] = topology[i];
    for(i=0;i<nodes;i++) 
      topology[i] = oldtopology[porder[i]-1];
  }
}




int LoadUniversalMesh(struct FemType *data,struct BoundaryType *bound,
          char *prefix,int info)
     /* Load the grid in universal file format */
{
  int noknots,totknots,noelements,elemcode,maxnodes;
  int allocated,maxknot,dim,ind,lines;
  int reordernodes,reorderelements,nogroups,maxnodeind,maxelem,elid,unvtype,elmertype;
  int nonodes,group,grouptype,mode,nopoints,nodeind,matind,physind,colorind;
  int debug,mingroup,maxgroup,nogroup,noentities,dummy;
  int *u2eind,*u2eelem;
  int *elementtypes;
  char filename[MAXFILESIZE],line[MAXLINESIZE],*cp;
  int i,j,k,l,n;
  char entityname[MAXNAMESIZE];
  FILE *in;


  strcpy(filename,prefix);
  if ((in = fopen(filename,"r")) == NULL) {
    AddExtension(prefix,filename,"unv");
    if ((in = fopen(filename,"r")) == NULL) {
      printf("LoadUniversalMesh: opening of the universal mesh file '%s' wasn't succesfull !\n",
       filename);
      return(1);
    }
  }
 
  printf("Reading mesh from universal mesh file %s.\n",filename);
  InitializeKnots(data);

  dim = 3;
  allocated = FALSE;
  reordernodes = FALSE;
  reorderelements = FALSE;

  debug = FALSE;
  if( debug ){
    elementtypes = Ivector(0,820);
    for(i=0;i<=820;i++) elementtypes[i] = FALSE;
  }

  maxnodeind = 0;
  maxnodes = 0;
  maxelem = 0;


omstart:

  if(info) {
    if(allocated) 
      printf("Second round for reading data\n");
    else 
      printf("First round for allocating data\n");
  }

  noknots = 0;
  noelements = 0;
  nogroups = 0;
  nopoints = 0;
  group = 0;


  for(;;) { 

    if(0) printf("line: %d  %s\n",mode,line);

  nextline:
    if( !strncmp(line,"    -1",6)) mode = 0;
    if( Getrow(line,in,FALSE)) goto end;
    if(!line) goto end;

    if( !strncmp(line,"    -1",6)) mode = 0;
    else if( !strncmp(line,"  2411",6)) mode = 2411;
    else if( !strncmp(line,"  2412",6)) mode = 2412;
    else if( !strncmp(line,"  2467",6)) mode = 2467;
    else if( !strncmp(line,"  2435",6)) mode = 2435;
    else if( !strncmp(line,"   781",6)) mode = 781;
    else if( !strncmp(line,"   780",6)) mode = 780;
    else if( allocated && strncmp(line,"      ",6)) printf("Unknown mode: %s",line);


    if(debug && mode) printf("Current mode is %d\n",mode);

    /* node definition */
    if( mode == 2411 || mode == 781 ) {
      if(debug) printf("Reading nodes in mode %d\n",mode);
      for(;;) {
  GetrowDouble(line,in);
  if( !strncmp(line,"    -1",6)) goto nextline;

  cp = line;
  nodeind = next_int(&cp);
  /* Three other fields omitted: two coordinate systems and color */
  noknots += 1;
  GetrowDouble(line,in);
  
  if(allocated) {
    if(reordernodes) {
      if(u2eind[nodeind]) 
        printf("Reordering node %d already set (%d vs. %d)\n",
         nodeind,u2eind[nodeind],noknots);
      else
        u2eind[nodeind] = noknots;
    }

    cp = line;
    data->x[noknots] = next_real(&cp);
    data->y[noknots] = next_real(&cp);
    data->z[noknots] = next_real(&cp);
  }
  else {
    if(nodeind != noknots) reordernodes = TRUE;
    maxnodeind = MAX(maxnodeind,nodeind);
  }
      }
    }

    if( mode == 2412 ) {
      if(debug) printf("Reading elements from field %d\n",mode);
      for(;;) {
  Getrow(line,in,FALSE);
  if( !strncmp(line,"    -1",6)) goto nextline;
  
  noelements += 1;
  cp = line;
  elid = next_int(&cp);
  unvtype = next_int(&cp);
  physind = next_int(&cp);
  matind = next_int(&cp);
  colorind = next_int(&cp);
  nonodes = next_int(&cp);
  
  if (!allocated) {
    maxnodes = MAX(maxnodes, nonodes);
    if(elid != noelements) reorderelements = TRUE;
    maxelem = MAX(maxelem, elid);
  }
  
  if(unvtype == 11 || unvtype == 21 || unvtype == 22 ) Getrow(line,in,FALSE);
  Getrow(line,in,FALSE);
  cp = line;

  elmertype = UnvToElmerType(unvtype); 
  if(!elmertype) {
    printf("Unknown elementtype %d %d %d %d %d %d\n",
     elid,unvtype,physind,matind,colorind,nonodes);
    printf("line: %s\n",line);
    bigerror("done");
  }

  if(elmertype == 510 )      
    lines = 1;
  else if(elmertype == 820 ) 
    lines = 2;
  else
    lines = 0;

  if(allocated) {
    if(reorderelements) u2eelem[elid] = noelements;

    if(debug && !elementtypes[elmertype]) {
      elementtypes[elmertype] = TRUE;
      printf("new elementtype in elmer: %d (unv: %d)\n",elmertype,unvtype);
    }

    if(elmertype % 100 != nonodes) {
      printf("nonodes = %d elemtype = %d elid = %d\n",nonodes,elmertype,elid);
      nonodes = elmertype % 100;
    }
    
      
    data->elementtypes[noelements] = elmertype;
    for(i=0;i<nonodes;i++) {
      if( lines > 0 && i >= 8 ) {
        if( i%8 == 0 ) {
    Getrow(line,in,FALSE);
    cp = line;
        }
      }
      data->topology[noelements][i] = next_int(&cp);
    }

    UnvRedundantIndexes(nonodes,data->topology[noelements]);

    UnvToElmerIndx(elmertype,data->topology[noelements]);   

    /* should this be physical property or material property? */
    data->material[noelements] = physind;
  }
  else {
    for(i=1;i<=lines;i++) 
      Getrow(line,in,FALSE);    
  }
      }
    }


    if( mode == 780 ) {
      int physind2,matind2;

      if(debug) printf("Reading elements from field %d\n",mode);
      for(;;) {
  Getrow(line,in,FALSE);
  if( !strncmp(line,"    -1",6)) goto nextline;
  
  noelements += 1;
  cp = line;
  elid = next_int(&cp);
  unvtype = next_int(&cp);

  physind = next_int(&cp);
  physind2 = next_int(&cp);
  matind = next_int(&cp);
  matind2 = next_int(&cp);
  colorind = next_int(&cp);
  nonodes = next_int(&cp);
  
  if (!allocated) {
    maxnodes = MAX(maxnodes, nonodes);
    if(elid != noelements) reorderelements = TRUE;
    maxelem = MAX(maxelem, elid);
  }
  
  if(unvtype == 11 || unvtype == 21) Getrow(line,in,FALSE);
  Getrow(line,in,FALSE);
  cp = line;
  if(allocated) {
    if(reorderelements) u2eelem[elid] = noelements;

    elmertype = UnvToElmerType(unvtype); 

    if(debug && !elementtypes[elmertype]) {
      elementtypes[elmertype] = TRUE;
      printf("new elementtype in elmer: %d (unv: %d)\n",elmertype,unvtype);
    }

    if(elmertype % 100 != nonodes) {
      printf("nonodes = %d elemtype = %d elid = %d\n",nonodes,elmertype,elid);
      nonodes = elmertype % 100;
    }

    data->elementtypes[noelements] = elmertype;
    for(i=0;i<nonodes;i++)
      data->topology[noelements][i] = next_int(&cp);

    UnvRedundantIndexes(nonodes,data->topology[noelements]);

    UnvToElmerIndx(elmertype,data->topology[noelements]);   

    /* should this be physical property or material property? */
    data->material[noelements] = physind;
  }
      }    
    }  

    if( mode == 2467 || mode == 2435) {
      if(debug) printf("Reading groups in mode %d\n",mode);
      
      for(;;) {
  Getrow(line,in,FALSE);
  if( !strncmp(line,"    -1",6)) goto nextline;
  
  cp = line;
  nogroup = next_int(&cp);
  for(i=1;i<=6;i++)
    dummy = next_int(&cp);
  noentities = next_int(&cp);

  Getrow(line,in,FALSE);  
  if( !strncmp(line,"    -1",6)) goto nextline;
  
  /* Used for the empty group created by salome */
  /* if( mode == 2467 && !strncmp(line,"            ",12)) continue; */
  
  group++;
  k = 0;
  if(allocated) {
    sscanf(line,"%s",entityname);
    strcpy(data->bodyname[group],entityname);
    data->bodynamesexist = TRUE;
    data->boundarynamesexist = TRUE;

    if(info) printf("Reading group %d with %d entities: %s\n",
        nogroup,noentities,entityname);
  }
  if(noentities == 0) Getrow(line,in,FALSE);

  for(i=0;i<noentities;i++) {
    if(i%2 == 0) {
      Getrow(line,in,FALSE);
      if( !strncmp(line,"    -1",6)) goto nextline;
      cp = line;
    }   

    grouptype = next_int(&cp);
    ind = next_int(&cp);
    dummy = next_int(&cp);
    dummy = next_int(&cp);

    if(ind == 0) continue;

    if( grouptype == 8 ) {
      if(allocated) {
        if(reorderelements) ind = u2eelem[ind];
        elemcode = data->elementtypes[ind];
        data->material[ind] = group;
      }
    }
    else if(grouptype == 7) {
      nopoints += 1;
      if(allocated) {
        elemcode = 101;
        data->material[noelements+nopoints] = group;
        data->elementtypes[noelements+nopoints] = elemcode;
        data->topology[noelements+nopoints][0] = ind;
      }
    }
    else {
    }
  }
  if(k && allocated && info)
    printf("Found new group %d with elements %d: %s\n",group,elemcode,entityname);

      }
    }

  }

end:

  exit;
  if(info) printf("Done reading\n");


  if(!allocated) {

    if(reordernodes) {
      if(info) printf("Reordering %d nodes with indexes up to %d\n",noknots,maxnodeind);
      u2eind = Ivector(1,maxnodeind);
      for(i=1;i<=maxnodeind;i++) u2eind[i] = 0;
    }
    if(reorderelements) {
      if(info) printf("Reordering %d elements with indexes up to %d\n",noelements,maxelem);
      u2eelem = Ivector(1,maxelem);
      for(i=1;i<=maxelem;i++) u2eelem[i] = 0;
    }

    if(noknots == 0 || noelements == 0 || maxnodes == 0) {
      printf("Invalid mesh consits of %d knots and %d %d-node elements.\n",
       noknots,noelements,maxnodes);     
      fclose(in);
      return(2);
    }

    rewind(in);
    totknots = noknots;
    data->noknots = noknots;
    data->noelements = noelements + nopoints;
    data->maxnodes = maxnodes;
    data->dim = dim;
    
    if(info) {
      printf("Allocating for %d knots and %d %d-node elements in %d dims.\n",
       noknots,noelements,maxnodes,dim);
    }  
    AllocateKnots(data);
    allocated = TRUE;

    goto omstart;    
  }
  fclose(in);


  if(reordernodes) {
    for(j=1;j<=noelements;j++)
      for(i=0;i<data->elementtypes[j]%100;i++)
  data->topology[j][i] = u2eind[data->topology[j][i]];
    free_Ivector(u2eind,1,maxnodeind);
  }
  if(reorderelements) {
    free_Ivector(u2eelem,1,maxelem);
  }


  mingroup = maxgroup = data->material[1];
  for(i=1;i<=data->noelements;i++) {
    mingroup = MIN( mingroup, data->material[i]);
    maxgroup = MAX( maxgroup, data->material[i]);
  }
  if(info) printf("The group interval is [%d,%d]\n",mingroup,maxgroup);
  if(mingroup == 0) {   
    if(info) {
      if(!maxgroup) printf("No material groups were successfully applied\n");
      printf("Unset elements were given material index %d\n",maxgroup+1);    
    }
    for(i=1;i<=data->noelements;i++) 
      if(data->material[i] == 0) data->material[i] = maxgroup + 1;
  }
    
  ElementsToBoundaryConditions(data,bound,TRUE,info);
 
  if(info) printf("The Universal mesh was loaded from file %s.\n\n",filename);

  return(0);
}




int LoadCGsimMesh(struct FemType *data,char *prefix,int info)
/* Load the mesh from postprocessing format of CGsim */
{
  int noknots,noelements,maxnodes,material,allocated,dim,debug,thismat,thisknots,thiselems;
  char filename[MAXFILESIZE],line[MAXLINESIZE],*cp;
  int i,j,inds[MAXNODESD2],savedofs;
  Real dummyreal;
  FILE *in;


  strcpy(filename,prefix);
  if ((in = fopen(filename,"r")) == NULL) {
    AddExtension(prefix,filename,"plt");
    if ((in = fopen(filename,"r")) == NULL) {
      printf("LoadCGsimMesh: opening of the CGsim mesh file '%s' wasn't succesfull !\n",
       filename);
      return(1);
    }
  }

  printf("Reading mesh from CGsim mesh file %s.\n",filename);
  InitializeKnots(data);

  debug = FALSE;
  allocated = FALSE;
  savedofs = FALSE;

omstart:

  maxnodes = 4;
  noknots = 0;
  noelements = 0;
  material = 0;
  dim = 2;
  thismat = 0;


  for(;;) {
    

    if(Getrow(line,in,FALSE)) goto end;
    if(!line) goto end;
    
    cp = strstr(line,"ZONE");    
    if(!cp) continue;
    

    thismat += 1;
    cp = strstr(line," N=");
    cp += 3;
    thisknots = next_int(&cp);

    cp = strstr(line,",E=");
    cp += 3;
    thiselems = next_int(&cp);

    if(debug) {
      printf("%s",line);
      printf("thismat = %d knots = %d elems = %d\n",thismat,thisknots,thiselems);
    }

    for(i=1;i<=thisknots;i++) {
      getline;

      if(allocated) {
  cp = line;
  data->x[noknots+i] = next_real(&cp);
  data->y[noknots+i] = next_real(&cp);
  data->z[noknots+i] = 0.0;

  if(savedofs == 1) {
    for(j=1;j<=4;j++)
      dummyreal = next_real(&cp);     
    data->dofs[1][noknots+i] = next_real(&cp);
  }
  else if(savedofs == 5) {
    for(j=1;j<=5;j++)   
      data->dofs[j][noknots+i] = next_real(&cp);    
  }

      }
    }

    for(i=1;i<=thiselems;i++) {
      getline;

      if(allocated) {
  cp = line;
  for(j=0;j<4;j++) 
    inds[j] = next_int(&cp);
  for(j=0;j<4;j++) 
    data->topology[noelements+i][j] = inds[j]+noknots;
  if(inds[2] == inds[3]) 
    data->elementtypes[noelements+i] = 303;
  else
    data->elementtypes[noelements+i] = 404;   
  data->material[noelements+i] = thismat;
      }
   }

    noknots += thisknots;
    noelements += thiselems;
  }

 end:

  if(!allocated) {
    if(noknots == 0 || noelements == 0 || maxnodes == 0) {
       printf("Invalid mesh consits of %d knots and %d %d-node elements.\n",
       noknots,noelements,maxnodes);     
       fclose(in);
       return(2);
    }

    rewind(in);
    data->noknots = noknots;
    data->noelements = noelements;
    data->maxnodes = maxnodes;
    data->dim = dim;

    
    if(info) {
      printf("Allocating for %d knots and %d %d-node elements.\n",
       noknots,noelements,maxnodes);
    }  
    AllocateKnots(data);

    if(savedofs == 1) {
      CreateVariable(data,1,1,0.0,"Temperature",FALSE);
    }
    else if(savedofs == 5) {
      CreateVariable(data,1,1,0.0,"dTdX",FALSE);
      CreateVariable(data,2,1,0.0,"dTdY",FALSE);
      CreateVariable(data,3,1,0.0,"Qx",FALSE);
      CreateVariable(data,4,1,0.0,"Qy",FALSE);
      CreateVariable(data,5,1,0.0,"Temperature",FALSE);
    }

    allocated = TRUE;
    goto omstart;    
  }
  fclose(in);

  if(info) printf("The CGsim mesh was loaded from file %s.\n\n",filename);
  return(0);
}