/* ElmerGrid - A simple mesh generation and manipulation utility Copyright (C) 1995- , CSC - Scientific Computing Ltd. Author: Peter Råback Email: Peter.Raback@csc.fi Address: CSC - Scientific Computing 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. */ /* --------------------: femfilein.c :-------------------------- */ #include #include #include #include #include #include #include #include #include "nrutil.h" #include "common.h" #include "femdef.h" #include "femtools.h" #include "femtypes.h" #include "femknot.h" #include "femfilein.h" #define getline fgets(line,MAXLINESIZE,in) char line[MAXLINESIZE]; static int Getrow(char *line1,FILE *io,int upper) { int i,isend; char line0[MAXLINESIZE],*charend; for(i=0;i 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;itopology[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 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(0) { printf("******\n"); printf("hits=%d ind=%d type=%d sides=%d\n",nohits,elemind,elemtype,elemsides); printf("sidenodes=%d side=%d\n",sidenodes,side); for(i=0;itopology[elemind][i]); printf("\n"); hit = TRUE; } 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;iparent2[j]) continue; GetElementSide(bound->parent[j],bound->side[j],1,data,&sideind2[0],&sideelemtype2); if(sideelemtype != sideelemtype2) continue; bcsame = 0; for(i=0;imaterial[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,*boundindx; char filename[MAXFILESIZE]; char line[MAXLINESIZE]; int i,j,*ind; 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;itopology[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;itopology[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;imaterial[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; } static int ReadAbaqusField(FILE *in,char *buffer,int *argtype,int *argno) /* This subroutine reads the Abaqus file format and tries to make sence out of it. */ { int i,val,digits; static int maxargno=0,mode=0; val = fgetc(in); if(val==EOF) return(-1); if(val=='\n') val = fgetc(in); if(val=='*') { if(0) printf("start field\n"); if((*argno) != maxargno) printf("The previous field was of wrong length, debugging time!\n"); (*argno) = 0; mode = 0; val = fgetc(in); if(val=='\n') val = fgetc(in); } if(val=='I') { for(i=0;i<2;i++) { val = fgetc(in); if(val=='\n') val = fgetc(in); buffer[i] = val; } buffer[2] = '\0'; digits = atoi(buffer); for(i=0;imaxnodes = 9; data->dim = 3; for(;;) { mode = ReadAbaqusField(in,buffer,&argtype,&argno); if(0) printf("%d %d: buffer: %s\n",argtype,argno,buffer); switch (mode) { case -1: goto jump; case 0: break; case 1921: /* General info */ if(argno == 3) printf("Reading output file for Abaqus %s\n",buffer); else if(argno == 4) printf("Created on %s",buffer); else if(argno == 5) printf("%s",buffer); else if(argno == 6) printf("%s\n",buffer); else if(argno == 7) data->noelements = atoi(buffer); else if(argno == 8 && allocated == FALSE) { data->noknots = atoi(buffer); allocated = TRUE; AllocateKnots(data); indx = Ivector(0,2 * data->noknots); for(i=1;i<=2*data->noknots;i++) indx[i] = 0; } break; case 1900: /* Element definition */ if(argno == 3) elemno = atoi(buffer); else if(argno == 4) { if(strstr(buffer,"2D4") || strstr(buffer,"SP4") || strstr(buffer,"AX4")) data->elementtypes[elemno] = 404; else if(strstr(buffer,"2D8") || strstr(buffer,"AX8") || strstr(buffer,"S8R5")) data->elementtypes[elemno] = 408; else if(strstr(buffer,"3D8")) data->elementtypes[elemno] = 808; else printf("Unknown element code: %s\n",buffer); } else if(argno >= 5) data->topology[elemno][argno-5] = atoi(buffer); break; case 1901: /* Node definition */ if(argno == 3) { knotno++; if(atoi(buffer) > 2*data->noknots) printf("LoadAbaqusOutput: allocate more space for indx.\n"); else indx[atoi(buffer)] = knotno; } if(argno == 4) sscanf(buffer,"%le",&(data->x[knotno])); if(argno == 5) sscanf(buffer,"%le",&(data->y[knotno])); if(argno == 6) sscanf(buffer,"%le",&(data->z[knotno])); break; case 1933: /* Element set */ if(argno == 3) { elset++; strcpy(data->bodyname[elset],buffer); } case 1934: /* Element set continuation */ if(argno > 3) { elemno = atoi(buffer); data->material[elemno] = elset; } break; case 2001: /* Just ignore */ break; case 1: if(argno == 3) knotno = indx[atoi(buffer)]; if(argno == 5) secno = atoi(buffer); break; case 2: if(prevdog != mode) { prevdog = mode; nodogs++; CreateVariable(data,nodogs,1,0.0,"Temperature",FALSE); } break; /* Read vectors in nodes in elements */ case 11: if(prevdog != mode) { prevdog = mode; nodogs++; CreateVariable(data,nodogs,3,0.0,"Stress",FALSE); } case 12: if(prevdog != mode) { prevdog = mode; nodogs++; CreateVariable(data,nodogs,3,0.0,"Invariants",FALSE); } if(secno==1 && argno == 3) sscanf(buffer,"%le",&(data->dofs[nodogs][3*knotno-2])); if(secno==1 && argno == 4) sscanf(buffer,"%le",&(data->dofs[nodogs][3*knotno-1])); if(secno==1 && argno == 5) sscanf(buffer,"%le",&(data->dofs[nodogs][3*knotno])); break; /* Read vectors in nodes. */ case 101: if(prevdog != mode) { prevdog = mode; nodogs++; CreateVariable(data,nodogs,3,0.0,"Displacement",FALSE); } case 102: if(prevdog != mode) { prevdog = mode; nodogs++; CreateVariable(data,nodogs,3,0.0,"Velocity",FALSE); } if(argno == 3) knotno = indx[atoi(buffer)]; if(argno == 4) sscanf(buffer,"%le",&(data->dofs[nodogs][3*knotno-2])); if(argno == 5) sscanf(buffer,"%le",&(data->dofs[nodogs][3*knotno-1])); if(argno == 6) sscanf(buffer,"%le",&(data->dofs[nodogs][3*knotno])); break; default: if(ignored != mode) { printf("Record %d was ignored!\n",mode); ignored = mode; } break; } } jump: if(info) printf("Renumbering elements\n"); for(j=1;j<=data->noelements;j++) for(i=0;i < data->elementtypes[j]%100;i++) data->topology[j][i] = indx[data->topology[j][i]]; free_ivector(indx,0,2*data->noknots); fclose(in); if(info) printf("LoadAbacusInput: results were loaded from file %s.\n",filename); return(0); } 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,nvalue,nvalue2,maxknot,minknot,nosides; int boundarytype,boundarynodes,nodes,maxnode; int *nodeindx,*boundindx; char filename[MAXFILESIZE]; char line[MAXLINESIZE],*cp; int i,j,k,l,*ind; FILE *in; Real rvalues[MAXDOFS]; int ivalues[MAXDOFS],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;jtopology[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;jtopology[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;jtopology[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,j,oldtopology[MAXNODESD2],*topology,elementtype,dim; int order808[]={1,2,4,3,5,6,8,7}; int order408[]={1,3,5,7,2,4,6,8}; int order306[]={1,3,5,2,4,6}; int order203[]={1,3,2}; int order605[]={1,2,4,3,5}; 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;ielementtypes[element] = 306; for(i=0;ielementtypes[element] = 408; for(i=0;ielementtypes[element] = 504; } else if(nodes == 5) { data->elementtypes[element] = 605; for(i=0;ielementtypes[element] = 706; } else if(nodes == 8) { for(i=0;idim); } int LoadFidapInput(struct FemType *data,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,allocated,nvalue,maxknot,totelems,entity,maxentity; char filename[MAXFILESIZE]; char line[MAXLINESIZE],*cp,entityname[MAXNAMESIZE],entitylist[100][MAXNAMESIZE]; int i,j,k,*ind,geoflag,typeflag; int **topology; FILE *in; Real rvalues[MAXDOFS]; Real *vel,*temp; int ivalues[MAXDOFS]; int nogroups,knotno,nonodes,elemno; 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;ielementtypes[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++) { #if 0 k = strcmp(entityname,entitylist[entity]); #else k = strcmp(entityname,data->bodyname[entity]); #endif if(k == 0) break; } if(entity > maxentity) { maxentity++; #if 0 strcpy(entitylist[entity],entityname); #else strcpy(data->bodyname[entity],entityname); #endif 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;jtopology[i][j] = next_int(&cp); ReorderFidapNodes(data,i,nodes,typeflag); if(data->elementtypes[i] == 0) { printf("******** nolla\n"); } if(entity) data->material[i] = entity; else data->material[i] = group; } totelems += elems; } mode = 0; break; case 10: 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"); return(0); } static void ReorderAnsysNodes(struct FemType *data,int *oldtopology, int element,int dim,int nodes) { int i,j,*topology,elementtype; 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;idim = 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;ibodynamesexist = 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 int FindFemlabParents(struct FemType *data,struct BoundaryType *bound, int nosides,int sidenodes,int *boundindx,int **nodeindx) { int i,j,k,l,sideelemtype,elemind,parent,*indx,*twosided; int boundarytype,minboundary,maxboundary,sideelem; int sideind[MAXNODESD1],elemsides,side,hit,sameelem,same,nosame; if(nosides < 1) return(0); maxboundary = minboundary = boundindx[1]; for(i=1;i<=nosides;i++) { if(maxboundary < boundindx[i]) maxboundary = boundindx[i]; if(minboundary > boundindx[i]) minboundary = boundindx[i]; } printf("Boundary types are in interval [%d, %d]\n",minboundary,maxboundary); twosided = Ivector(1,nosides); indx = Ivector(1,data->noknots); sideelem = 0; sameelem = 0; for(boundarytype=minboundary;boundarytype <= maxboundary;boundarytype++) { for(i=1;i<=data->noknots;i++) indx[i] = 0; for(i=1;i<=nosides;i++) { if(boundindx[i] == boundarytype) for(j=1;j<=sidenodes;j++) indx[nodeindx[i][j]] = TRUE; } for(i=1;i<=nosides;i++) twosided[i] = 0; for(elemind=1;elemind<=data->noelements;elemind++) { elemsides = data->elementtypes[elemind] / 100; if(elemsides == 8) elemsides = 6; else if(elemsides == 6) elemsides = 5; else if(elemsides == 5) elemsides = 4; if(0) printf("elem=%d\n",elemind); for(side=0;side data->noknots) { printf("sideind[%d] = %d (noknots=%d)\n", i,sideind[i],data->noknots); hit = FALSE; } else if(!indx[sideind[i]]) hit = FALSE; } if(hit == TRUE) { same = FALSE; for(i=1;i<=nosides;i++) { if(boundindx[i] == boundarytype) { nosame = 0; for(k=0;kparent2[twosided[i]] = elemind; bound->side2[twosided[i]] = side; } else { sideelem += 1; if(same) twosided[i] = sideelem; bound->parent[sideelem] = elemind; bound->side[sideelem] = side; bound->parent2[sideelem] = 0; bound->side2[sideelem] = 0; bound->types[sideelem] = boundarytype; } same = FALSE; } } } } if(sameelem) printf("Found %d side elements that have two parents.\n",sameelem); if(sideelem == nosides) printf("Found correctly %d side elements.\n",sideelem); else { printf("Found %d side elements, should have found %d\n",sideelem,nosides); bound->nosides = sideelem; } free_Ivector(indx,1,data->noknots); free_Ivector(twosided,1,nosides); return(0); } int LoadFemlabMesh(struct FemType *data,struct BoundaryType *bound, char *prefix,int info) /* This procedure reads the FEM mesh as written by Femlab and a special Matlab routine . */ { int noknots,noelements,nosides,elemcode,sideelem,sidetype,currenttype; int elementtypes,sideind[MAXNODESD1],tottypes,elementtype; int i,j,k,l,imax,grp,dummyint,**nodeindx,*boundindx,boundarynodes,maxnodes; int dim,nonodes,sidenodes; Real x,y,z,dummy; FILE *in; char *cp,line[MAXLINESIZE],filename[MAXFILESIZE], text[MAXNAMESIZE],directoryname[MAXFILESIZE]; sprintf(filename,"%s.header",prefix); if ((in = fopen(filename,"r")) == NULL) { printf("LoadFemlabInput: The opening of the header-file %s failed!\n", filename); return(1); } for(i=0;getline;i++); fclose(in); in = fopen(filename,"r"); getline; sscanf(line,"%le",&dummy); dim = dummy+0.5; getline; sscanf(line,"%le",&dummy); noknots = dummy+0.5; getline; sscanf(line,"%le",&dummy); noelements = dummy+0.5; getline; sscanf(line,"%le",&dummy); nosides = dummy+0.5; getline; sscanf(line,"%le",&dummy); nonodes = dummy+0.5; fclose(in); printf("Femlab %dD file has %d %d-node elements, %d nodes and %d sides.\n", dim,noelements,nonodes,noknots,nosides); sprintf(filename,"%s.node",prefix); if ((in = fopen(filename,"r")) == NULL) { printf("LoadFemlabInput: The opening of the nodes-file %s failed!\n", filename); return(1); } InitializeKnots(data); /* Even 2D elements may form a 3D object! */ data->dim = 3; data->maxnodes = nonodes; data->noknots = noknots; data->noelements = noelements; AllocateKnots(data); sprintf(filename,"%s.node",prefix); if ((in = fopen(filename,"r")) == NULL) { printf("LoadFemlabInput: The opening of the nodes-file %s failed!\n", filename); return(3); } else printf("Loading %d Femlab nodes from %s\n",noknots,filename); for(i=1;i<=noknots;i++) { getline; cp=line; x = next_real(&cp); y = next_real(&cp); z = next_real(&cp); data->x[i] = x; data->y[i] = y; data->z[i] = z; } fclose(in); sprintf(filename,"%s.elem",prefix); if ((in = fopen(filename,"r")) == NULL) { printf("LoadFemlabInput: The opening of the element-file %s failed!\n", filename); return(4); } else printf("Loading %d Femlab elements from %s\n",noelements,filename); for(j=1;j<=noelements;j++) { getline; cp=line; for(i=0;itopology[j][i] = (int)(dummy); } dummy = next_real(&cp)+0.5; data->material[j] = (int)(dummy); if(nonodes == 3) data->elementtypes[j] = 303; else if(nonodes == 4) data->elementtypes[j] = 504; } fclose(in); sprintf(filename,"%s.boundary",prefix); if ((in = fopen(filename,"r")) == NULL) { printf("LoadFemlabInput: The opening of the boundary-file %s failed!\n", filename); return(5); } else printf("Reading Femlab boundaries from file %s\n",filename); AllocateBoundary(bound,nosides); if(nonodes == 3) sidenodes = 2; else if(nonodes == 4) sidenodes = 3; nodeindx = Imatrix(1,nosides,1,sidenodes); boundindx = Ivector(1,nosides); j = 0; for(j=1;j<=nosides;j++) { getline; cp=line; for(i=1;i<=sidenodes;i++) { dummy = next_real(&cp)+0.5; nodeindx[j][i] = (int)(dummy); } dummy = next_real(&cp)+0.5; boundindx[j] = (int)(dummy); } fclose(in); FindFemlabParents(data,bound,nosides,sidenodes,boundindx,nodeindx); free_Imatrix(nodeindx,1,nosides,1,sidenodes); free_Ivector(boundindx,1,nosides); return(0); } static void ReorderFieldviewNodes(struct FemType *data,int *oldtopology, int element,int dim,int nodes) { int i,j,*topology,elementtype; 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;idim = 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;jtopology = 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 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); 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,nosides,dim; int sideind[MAXNODESD1],elemind[MAXNODESD2],tottypes,elementtype,bcmarkers; int i,j,k,dummyint,*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 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; int elemsides,sidelemtype,sideind[2],side,sideelemtype,elemind; 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 %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;sideparent[i] = elemind; bound->side[i] = side; bound->parent2[i] = 0; bound->side2[i] = 0; bound->types[i] = bctype; } } } } printf("Succesfully 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,noelements,maxnodes,elematts,nodeatts,nosides,dim; int sideind[MAXNODESD1],elemind[MAXNODESD2],tottypes,elementtype,bcmarkers; int i,j,k,dummyint,*boundnodes,allocated; FILE *in; char *cp,line[MAXLINESIZE],elemfile[MAXFILESIZE],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("Succesfully 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,nvalue,maxknot,nosides,sideelemtype; int boundarytype,materialtype,boundarynodes,side,parent,elemsides; int dim, elemnodes, elembasis, elemtype, bulkdone, usedmax,hits; int minbulk,maxbulk,minbound,maxbound,label,debug; int *usedno, **usedelem; char filename[MAXFILESIZE],line[MAXLINESIZE],*cp; int i,j,k,l,n,ind,inds[MAXNODESD2],sideind[MAXNODESD1]; FILE *in; Real x,y,z; 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;itopology[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;itopology[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;itopology[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 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;sideparent[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,j,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 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;jtopology[noelements][j] = next_int(&cp) + offset; if(1) ReorderComsolNodes(elemtype,data->topology[noelements]); } } } else if(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); return(0); } static int GmshToElmerType(int gmshtype) { int elmertype; 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; default: printf("Gmsh element %d does not have an Elmer counterpart!\n",gmshtype); } return(elmertype); } static void GmshToElmerIndx(int elemtype,int elemind[]) { int tmpind[MAXNODESD2]; switch (elemtype) { case 510: tmpind[8] = elemind[8]; tmpind[9] = elemind[9]; elemind[8] = tmpind[9]; elemind[9] = tmpind[8]; break; /* There seems to be conflicting data whether this is needed or not */ case 306: tmpind[4] = elemind[4]; tmpind[5] = elemind[5]; elemind[4] = tmpind[5]; elemind[5] = tmpind[4]; break; } } static int LoadGmshInput1(struct FemType *data,struct BoundaryType *bound, char *filename,int info) { int noknots,noelements,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;jtopology[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;jtopology[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,info); printf("Succesfully 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,noelements,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; 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; omstart: for(;;) { if(Getrow(line,in,TRUE)) goto end; if(!line) goto end; if(strstr(line,"END")) continue; /* Header info is not much needed */ if(!strncasecmp(line,"$MeshFormat",11)) { Getrow(line,in,TRUE); Getrow(line,in,TRUE); if(strncasecmp(line,"$EndMeshFormat",14)) { printf("MeshFormat section should end to string $EndMeshFormat\n"); printf("%s\n",line); } } else if(!strncasecmp(line,"$Nodes",6)) { 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; } else if(!strncasecmp(line,"$Elements",9)) { 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 */ tagphys = 0; if(gmshtype != 15) { 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); } data->material[i] = tagphys; for(j=0;jtopology[i][j] = elemind[j]; } else { maxelemtype = MAX(maxelemtype,elementtype); } } getline; } 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;jtopology[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); } if(1) ElementsToBoundaryConditions(data,bound,info); if(info) printf("Succesfully 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,TRUE); fclose(in); if(strstr(line,"MESHFORMAT")) errno = LoadGmshInput2(data,bound,filename,info); else errno = LoadGmshInput1(data,bound,filename,info); return(errno); } int UnvToElmerType(int unvtype) { int elmertype; switch (unvtype) { case 11: case 21: elmertype = 202; break; case 41: case 74: case 91: elmertype = 303; break; case 42: case 51: case 62: case 72: case 92: elmertype = 306; break; case 44: case 54: case 64: case 71: case 94: elmertype = 404; break; case 111: elmertype = 504; break; case 118: elmertype = 513; 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); } int LoadUniversalMesh(struct FemType *data,char *prefix,int info) /* Load the grid in universal file format */ { int noknots,totknots,noelements,elemcode,maxnodes; int allocated,maxknot,dim,ind; int reordernodes,reorderelements,nogroups,maxnode,maxelem,elid,unvtype,elmertype; int nonodes,group,grouptype,mode,nopoints; int debug,mingroup,maxgroup; int *u2eind,*u2eelem; 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; debug = FALSE; allocated = FALSE; reordernodes = FALSE; reorderelements = FALSE; omstart: maxnodes = 0; nogroups = 0; maxnode = 0; maxelem = 0; noknots = 0; noelements = 0; nopoints = 0; group = 0; for(;;) { 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(0 && allocated && strncmp(line," ",6)) printf("Unknown command: %s",line); /* node definition */ if( mode == 2411) { if(0 && info && allocated) printf("Reading nodes\n"); for(;;) { Getrow(line,in,FALSE); if( !strncmp(line," -1",6)) goto nextline; cp = line; i = next_int(&cp); noknots += 1; if(i != noknots) reordernodes = TRUE; maxnode = MAX(maxnode,i); Getrow(line,in,FALSE); if(allocated) { if(reordernodes) { u2eind[i] = noknots; } cp = line; data->x[noknots] = next_real(&cp); data->y[noknots] = next_real(&cp); data->z[noknots] = next_real(&cp); } } } if( mode == 2412) { if(0 && info && allocated) printf("Reading elements\n"); for(;;) { Getrow(line,in,FALSE); if( !strncmp(line," -1",6)) goto nextline; noelements += 1; cp = line; elid = next_int(&cp); unvtype = next_int(&cp); i = next_int(&cp); i = next_int(&cp); i = 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) Getrow(line,in,FALSE); Getrow(line,in,FALSE); cp = line; if(allocated) { if(reorderelements) u2eelem[elid] = noelements; elmertype = UnvToElmerType(unvtype); if(elmertype%100 != nonodes) printf("nonodes = %d elemtype = %d elid = %d\n",nonodes,elmertype,elid); data->elementtypes[noelements] = elmertype; for(i=0;itopology[noelements][i] = next_int(&cp); } } } if( mode == 2467 || mode == 2435) { if(0 && info && allocated) printf("Reading groups\n"); Getrow(line,in,FALSE); if( !strncmp(line," -1",6)) goto nextline; for(;;) { 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; } for(;;) { Getrow(line,in,FALSE); if( !strncmp(line," -1",6)) goto nextline; cp = line; for(i=1;i<=2;i++) { grouptype = next_int(&cp); ind = next_int(&cp); if( ind == 0 && i==1) goto newgroup; if( ind == 0 && i==2) continue; k++; /* Temperary exception: jump over nodal element groups */ if(grouptype == 7 && mode == 2435) { if(allocated) printf("Note: in field 2435 nodal point groups are currently omitted!\n"); group--; goto newgroup; } j = next_int(&cp); j = next_int(&cp); 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 goto newgroup; if(k == 1 && allocated && info) printf("Found new group %d with elements %d: %s\n",group,elemcode,entityname); } } newgroup: continue; } } } end: exit; printf("done reading\n"); if(!allocated) { if(reordernodes) { if(info) printf("Reordering %d nodes with indexes up to %d\n",noknots,maxnode); u2eind = Ivector(1,maxnode); for(i=1;i<=maxnodes;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.\n", noknots,noelements,maxnodes); } AllocateKnots(data); allocated = TRUE; goto omstart; } fclose(in); if(reordernodes) { for(j=1;j<=noelements;j++) for(i=0;ielementtypes[j]%100;i++) data->topology[j][i] = u2eind[data->topology[j][i]]; free_Ivector(u2eind,1,maxnode); } 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(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; } if(info) printf("The Universal mesh was loaded from file %s.\n\n",filename); return(0); }