/* 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. */ /* --------------------------: femknot.c :---------------------------- This module includes subroutines that formulate the mesh into structures more usuful for the user. These are the functions that should be used for assembling. The routines usually operate on structures FemType and BoundaryType. */ #include #include #include #include #include "common.h" #include "nrutil.h" #include "femdef.h" #include "femtypes.h" #include "femmesh.h" #include "femknot.h" #include "femtools.h" #include "femsolve.h" #define DEBUG 0 void GetElementInfo(int element,struct FemType *data, Real *globalcoord,int *ind,int *material) /* For a given element gives the coordinates and index for the knot. This subroutine uses the standard formulation which uses up more memory, but is easier to understand. It requires that the knots must first be stored in struct FemType. */ { int i,indi,nodesd2; nodesd2 = data->elementtypes[element]%100; for(i=0;itopology[element][i]; globalcoord[i] = data->x[indi]; globalcoord[i+nodesd2] = data->y[indi]; } (*material) = data->material[element]; } int GetElementDimension(int elementtype) { int elemdim; switch (elementtype / 100) { case 1: elemdim = 0; break; case 2: elemdim = 1; break; case 3: case 4: elemdim = 2; break; case 5: case 6: case 7: case 8: elemdim = 3; break; default: printf("GetElementDimension: unknown elementtype %d\n",elementtype); } return(elemdim); } int GetMaxElementType(struct FemType *data) { int i,maxelementtype; maxelementtype = data->elementtypes[1]; for(i=1;i<=data->noelements;i++) if(data->elementtypes[i] > maxelementtype) maxelementtype = data->elementtypes[i]; return(maxelementtype); } int GetMinElementType(struct FemType *data) { int i,minelementtype; minelementtype = data->elementtypes[1]; for(i=1;i<=data->noelements;i++) if(data->elementtypes[i] < minelementtype) minelementtype = data->elementtypes[i]; return(minelementtype); } int GetMaxElementDimension(struct FemType *data) { int maxelementtype,elemdim; maxelementtype = GetMaxElementType(data); elemdim = GetElementDimension(maxelementtype); return(elemdim); } void GetElementSide(int element,int side,int normal, struct FemType *data,int *ind,int *sideelemtype) /* Give the indices of a given side of a given element. The subroutine is valid for 4, 5, 8, 9, 12 and 16 node rectangular elements, and 3 and 6 node triangular elements. */ { int i,j,elemtype,*elemind,sides,ind2[MAXNODESD2]; elemtype = data->elementtypes[element]; elemind = data->topology[element]; sides = elemtype/100; if(side < 0 && sides > 4) side = -(side+1); switch (elemtype) { case 202: case 203: case 204: *sideelemtype = 101; ind[0] = elemind[side]; break; case 303: /* Linear triangle */ if(side < 3) { *sideelemtype = 202; ind[0] = elemind[side]; ind[1] = elemind[(side+1)%3]; } else { *sideelemtype = 101; ind[0] = elemind[side-3]; } break; case 306: /* 2nd order triangle */ if(side < 3) { *sideelemtype = 203; ind[0] = elemind[side]; ind[1] = elemind[(side+1)%3]; ind[2] = elemind[side+3]; } else { *sideelemtype = 101; ind[0] = elemind[side-3]; } break; case 310: /* 3rd order triangle */ if(side < 3) { *sideelemtype = 204; ind[0] = elemind[side]; ind[1] = elemind[(side+1)%3]; ind[2] = elemind[2*side+3]; ind[3] = elemind[2*side+4]; } else { *sideelemtype = 101; ind[0] = elemind[side-3]; } break; case 404: /* Linear quadrilateral */ case 405: if(side < 4) { *sideelemtype = 202; ind[0] = elemind[side]; ind[1] = elemind[(side+1)%4]; } else { *sideelemtype = 101; ind[0] = elemind[side-4]; } break; case 408: /* 2nd order quadrilateral */ case 409: if(side < 4) { *sideelemtype = 203; ind[0] = elemind[side]; ind[1] = elemind[(side+1)%4]; ind[2] = elemind[side+4]; } else { *sideelemtype = 101; ind[0] = elemind[side-4]; } break; case 412: /* 3rd order quadrilateral */ case 416: if(side < 4) { *sideelemtype = 204; ind[0] = elemind[side]; ind[1] = elemind[(side+1)%4]; ind[2] = elemind[2*side+4]; ind[3] = elemind[2*side+5]; } else { *sideelemtype = 101; ind[0] = elemind[side-4]; } break; case 504: /* Linear tetrahedron */ if(side < 4) { *sideelemtype = 303; if(side < 3) { ind[0] = elemind[side]; ind[1] = elemind[(side+1)%3]; ind[2] = elemind[3]; } if(side == 3) { ind[0] = elemind[0]; ind[1] = elemind[2]; ind[2] = elemind[1]; } } else if(side < 10) { *sideelemtype = 202; if(side < 7) { ind[0] = elemind[side-4]; ind[1] = elemind[3]; } else { ind[0] = elemind[side-7]; ind[1] = elemind[(side-6)%3]; } } else { *sideelemtype = 101; ind[0] = elemind[side-10]; } break; case 510: /* 2nd order tetrahedron */ if(side < 4) { *sideelemtype = 306; if(side < 3) { ind[0] = elemind[side]; ind[1] = elemind[(side+1)%3]; ind[2] = elemind[3]; ind[3] = elemind[4+side]; ind[4] = elemind[7+(side+1)%3]; ind[5] = elemind[7+side]; } else if(side == 3) { ind[0] = elemind[0]; ind[1] = elemind[1]; ind[2] = elemind[2]; ind[3] = elemind[4]; ind[4] = elemind[5]; ind[5] = elemind[6]; } } else if(side < 10) { *sideelemtype = 203; if(side < 7) { ind[0] = elemind[side-4]; ind[1] = elemind[3]; ind[2] = elemind[side+3]; } else { ind[0] = elemind[side-7]; ind[1] = elemind[(side-6)%3]; ind[2] = elemind[side-3]; } } else { *sideelemtype = 101; ind[0] = elemind[side-10]; } break; case 706: /* Linear prism or vedge element */ if(side < 3) { *sideelemtype = 404; ind[0] = elemind[side]; ind[1] = elemind[(side+1)%3]; ind[2] = elemind[(side+1)%3+3]; ind[3] = elemind[side+3]; } else if (side < 5) { *sideelemtype = 303; for(i=0;i<3;i++) ind[i] = elemind[3*(side-3)+i]; } else { *sideelemtype = 101; ind[0] = elemind[side-5]; } break; case 605: /* Linear pyramid */ if(side < 4) { *sideelemtype = 303; ind[0] = elemind[side]; ind[1] = elemind[4]; ind[2] = elemind[(side+1)%4]; } else if (side == 4) { *sideelemtype = 404; for(i=0;i<3;i++) ind[i] = elemind[i]; } else { *sideelemtype = 101; ind[0] = elemind[side-4]; } break; case 613: /* 2nd order pyramid */ if(side < 4) { *sideelemtype = 306; ind[0] = elemind[side]; ind[1] = elemind[(side+1)%4]; ind[2] = elemind[4]; ind[3] = elemind[side+5]; ind[4] = elemind[(side+1)%4+9]; ind[5] = elemind[side%4+9]; } else if (side == 4) { *sideelemtype = 408; for(i=0;i<3;i++) ind[i] = elemind[i]; for(i=0;i<3;i++) ind[i+4] = elemind[i+5]; } else { *sideelemtype = 101; ind[0] = elemind[side-4]; } break; case 808: /* Linear brick */ if(side < 6) { *sideelemtype = 404; if(side < 4) { ind[0] = elemind[side]; ind[1] = elemind[(side+1)%4]; ind[2] = elemind[(side+1)%4+4]; ind[3] = elemind[side+4]; } else if(side < 6) { for(i=0;i<4;i++) ind[i] = elemind[4*(side-4)+i]; } } else if(side < 18) { *sideelemtype = 202; if(side < 10) { ind[0] = elemind[side-6]; ind[1] = elemind[(side-5)%4]; } else if(side < 14) { ind[0] = elemind[side-6]; ind[1] = elemind[(side-9)%4+4]; } else if(side < 18) { ind[0] = elemind[side-14]; ind[1] = elemind[side-14+4]; } } else if(side < 26) { *sideelemtype = 101; ind[0] = elemind[side-18]; } break; case 820: /* 2nd order brick */ *sideelemtype = 408; if(side < 4) { ind[0] = elemind[side]; ind[1] = elemind[(side+1)%4]; ind[2] = elemind[(side+1)%4+4]; ind[3] = elemind[side+4]; ind[4] = elemind[8+side]; ind[5] = elemind[12+(side+1)%4]; ind[6] = elemind[16+side]; ind[7] = elemind[12+side]; } else if(side < 6) { for(i=0;i<4;i++) ind[i] = elemind[4*(side-4)+i]; for(i=0;i<4;i++) ind[i+4] = elemind[8*(side-4)+8+i]; } break; case 827: *sideelemtype = 409; if(side < 4) { ind[0] = elemind[side]; ind[1] = elemind[(side+1)%4]; ind[2] = elemind[(side+1)%4+4]; ind[3] = elemind[side+4]; ind[4] = elemind[8+side]; ind[5] = elemind[12+(side+1)%4]; ind[6] = elemind[16+side]; ind[7] = elemind[12+side]; ind[8] = elemind[20+side]; } else { for(i=0;i<4;i++) ind[i] = elemind[4*(side-4)+i]; for(i=0;i<4;i++) ind[i+4] = elemind[8*(side-4)+8+i]; ind[8] = elemind[20+side]; } break; default: printf("GetElementSide: unknown elementtype %d (elem=%d,side=%d)\n",elemtype,element,side); } if(normal == -1) { if(*sideelemtype == 202 || *sideelemtype == 203 || *sideelemtype == 303 || *sideelemtype == 404) { j = *sideelemtype/100-1; for(i=0;i<=j;i++) ind2[i] = ind[i]; for(i=0;i<=j;i++) ind[i] = ind2[j-i]; } #if 0 else if(normal != 1) { printf("GetElementSide: unknown option (normal=%d)\n",normal); } #endif } } int GetElementGraph(int element,int edge,struct FemType *data,int *ind) { int i,j,elemtype,basetype,elemnodes,*elemind,hit,evenodd,quadratic,side; elemtype = data->elementtypes[element]; basetype = elemtype / 100; elemnodes = elemtype % 100; quadratic = (elemnodes > basetype); elemind = data->topology[element]; ind[0] = ind[1] = 0; if(quadratic) side = edge / 2; else side = edge; switch (basetype) { case 2: if(side == 0) { ind[0] = elemind[0]; ind[1] = elemind[1]; } break; case 3: if(side < 3) { ind[0] = elemind[side]; ind[1] = elemind[(side+1)%3]; } break; case 4: if(side < 4) { ind[0] = elemind[side]; ind[1] = elemind[(side+1)%4]; } break; case 5: if(side < 3) { ind[0] = elemind[side]; ind[1] = elemind[(side+1)%3]; } else if(side < 6) { ind[0] = elemind[side-3]; ind[1] = elemind[3]; } break; case 6: if(side < 4) { ind[0] = elemind[side]; ind[1] = elemind[(side+1)%4]; } else if(side < 8) { ind[0] = elemind[side-4]; ind[1] = elemind[4]; } break; case 7: if(side < 3) { ind[0] = elemind[side]; ind[1] = elemind[(side+1)%3]; } else if(side < 6) { ind[0] = elemind[side-3]; ind[1] = elemind[side]; } else if(side < 9) { ind[0] = elemind[side-3]; ind[1] = elemind[3+(side+1)%3]; } break; case 8: if(side < 4) { ind[0] = elemind[side]; ind[1] = elemind[(side+1)%4]; } else if(side < 8) { ind[0] = elemind[side-4]; ind[1] = elemind[side]; } else if(side < 12) { ind[0] = elemind[side-4]; ind[1] = elemind[4+(side+1)%4]; } break; } hit = (ind[0] || ind[1]); if(hit && quadratic) { evenodd = edge - 2*side; switch (basetype) { case 2: ind[evenodd] = elemind[2]; break; case 3: ind[evenodd] = elemind[side+3]; break; case 4: ind[evenodd] = elemind[side+4]; break; case 5: ind[evenodd] = elemind[side+4]; break; case 6: ind[evenodd] = elemind[side+5]; break; case 7: ind[evenodd] = elemind[side+6]; break; case 8: ind[evenodd] = elemind[side+8]; break; } } return(hit); } int CalculateIndexwidth(struct FemType *data,int indxis,int *indx) { int i,ind,nonodes,indexwidth; int imax,imin,element; /* Calculate the maximum bandwidth */ indexwidth = 0; for(element=1; element <= data->noelements; element++) { imin = data->noknots; imax = 0; nonodes = data->elementtypes[element]%100; for(i=0;itopology[element][i]; if(indxis) ind = indx[ind]; if(ind == 0) continue; if(ind > imax) imax = ind; if(ind < imin) imin = ind; } if(imax-imin > indexwidth) indexwidth = imax-imin; } if(!indxis) data->indexwidth = indexwidth; return(indexwidth); } void InitializeKnots(struct FemType *data) { int i; data->timesteps = 0; data->noknots = 0; data->noelements = 0; data->coordsystem = COORD_CART2; data->numbering = NUMBER_XY; data->created = FALSE; data->variables = 0; data->maxnodes = 0; data->indexwidth = 0; data->noboundaries = 0; data->mapgeo = 1; data->nocorners = 0; data->pelems = 0; data->boundarynamesexist = FALSE; data->bodynamesexist = FALSE; data->nopartitions = 1; data->partitionexist = FALSE; data->periodicexist = FALSE; data->connectexist = FALSE; data->dualexists = FALSE; data->invtopoexists = FALSE; data->partitiontableexists = FALSE; for(i=0;iedofs[i] = 0; data->bandwidth[i] = 0; data->iterdofs[i] = 0; strcpy(data->dofname[i],""); } for(i=0;ibodyname[i],""); sprintf(data->bodyname[i],"body%d",i); } for(i=0;iboundaryname[i],""); sprintf(data->boundaryname[i],"bc%d",i); } } void AllocateKnots(struct FemType *data) { int i; data->topology = Imatrix(1,data->noelements,0,data->maxnodes-1); data->material = Ivector(1,data->noelements); data->elementtypes = Ivector(1,data->noelements); for(i=1;i<=data->noelements;i++) data->material[i] = 0; for(i=1;i<=data->noelements;i++) data->elementtypes[i] = 0; data->x = Rvector(1,data->noknots); data->y = Rvector(1,data->noknots); data->z = Rvector(1,data->noknots); for(i=1;i<=data->noknots;i++) data->x[i] = data->y[i] = data->z[i] = 0.0; data->created = TRUE; #if DEBUG printf("Allocated for %d %d-node elements resulting to %d nodes\n", data->noelements,data->maxnodes,data->noknots); #endif } static void MovePointCircle(Real *lim,int points,Real *coords, Real x,Real y,Real *dx,Real *dy) { int i; Real x0,y0,r,r1,r2,p; r2 = fabs(lim[1]-lim[0]); if(r2 > fabs(lim[2]-lim[1])) r2 = fabs(lim[2]-lim[1]); for(i=0;i= r2) continue; y0 = y-(lim[1]+coords[2*i+1]); if(y0 < 0 && lim[0] > lim[1]) continue; if(y0 > 0 && lim[2] < lim[1]) continue; if(fabs(y0) >= r2) continue; r = sqrt(x0*x0+y0*y0); if(r < 1.0e-50) continue; if(fabs(x0) > fabs(y0)) { p = fabs(x0)/r - 1.0; if(fabs(x0) <= r1) { *dx += p*x0; *dy += p*y0; } else if(fabs(x0) <= r2) { *dx += p*x0*(r2-fabs(x0))/(r2-r1); *dy += p*y0*(r2-fabs(x0))/(r2-r1); } } else { p = fabs(y0)/r - 1.0; if(fabs(y0) <= r1) { *dx += p*x0; *dy += p*y0; } else if(fabs(y0) <= r2) { *dx += p*x0*(r2-fabs(y0))/(r2-r1); *dy += p*y0*(r2-fabs(y0))/(r2-r1); } } } } static void MovePointLinear(Real *lim,int points,Real *coords, Real x,Real y,Real *dx,Real *dy) { static int i=0; Real c,d; if(y > lim[0] && y < lim[2]) { if(x <= coords[0]) { d = coords[1]; } else if(x >= coords[points-2]) { d = coords[points-1]; } else { i = 1; while(x > coords[2*i] && i < points/2-1) i++; c = (coords[2*i+1]-coords[2*i-1])/(coords[2*i]-coords[2*i-2]); d = coords[2*i-1] + c*(x-coords[2*i-2]); } if(y < lim[1]) *dy += d*(y-lim[0])/(lim[1]-lim[0]); else *dy += d*(lim[2]-y)/(lim[2]-lim[1]); } } static void MovePointAngle(Real *lim,int points,Real *coords, Real x,Real y,Real *dx,Real *dz) { static int i=0; Real c,d,x1,y1,z1,degs; degs = FM_PI/180.0; x1 = z1 = 0.0; if(y > lim[0] && y < lim[2]) { if(x <= coords[0]) { x1 = x; } else { i = 1; while(x > coords[2*i] && i <= points/2-1) { x1 = x1 + cos(degs*coords[2*i-1])*(coords[2*i]-coords[2*i-2]); z1 = z1 + sin(degs*coords[2*i-1])*(coords[2*i]-coords[2*i-2]); i++; } x1 = x1 + cos(degs*coords[2*i-1])*(x-coords[2*i-2]); z1 = z1 + sin(degs*coords[2*i-1])*(x-coords[2*i-2]); } if(y < lim[1]) { *dx += (x1-x)*(y-lim[0])/(lim[1]-lim[0]); *dz += z1*(y-lim[0])/(lim[1]-lim[0]); } else { *dx += (x1-x)*(lim[2]-y)/(lim[2]-lim[1]); *dz += z1*(lim[2]-y)/(lim[2]-lim[1]); } } } static void MovePointSinus(Real *lim,int points,Real *coords, Real x,Real y,Real *dx,Real *dy) { static int i=0; Real c,d; if(y > lim[0] && y < lim[2]) { if(x <= coords[0]) { d = 0.0; } else if(x >= coords[1]) { d = coords[3]*sin(coords[2]*2.*FM_PI); } else { c = coords[2]*2.*FM_PI/(coords[1]-coords[0]); d = coords[3]*sin(c*(x-coords[0])); } if(y < lim[1]) *dy += d*(y-lim[0])/(lim[1]-lim[0]); else *dy += d*(lim[2]-y)/(lim[2]-lim[1]); } } static void MovePointPowerSeries(Real *lim,int points,Real *coords, Real x,Real y,Real *dx,Real *dy) { int i,n; Real c,d,t,u; if(y > lim[0] && y < lim[2]) { t = x; if(coords[1] > coords[0]) { if(tcoords[1]) t = coords[1]; } else { if(t>coords[0]) t = coords[0]; if(t lim[0] && y < lim[2]) { t = x; if(coords[1] > coords[0]) { if(tcoords[1]) t = coords[1]; } else { if(t>coords[0]) t = coords[0]; if(t=1) d += u * coords[3]; for(i=2;i lim[0] && y < lim[2]) { if(x <= coords[0]) { d = coords[1]; } else if(x >= coords[points-2]) { d = coords[points-1]; } else { i = 1; while(x > coords[3*i] && i < points/3-1) i++; c = (coords[3*i+1]-coords[3*i-2])/pow((coords[3*i]-coords[3*i-3]),coords[3*i-1]); d = coords[3*i-2] + c*pow((x-coords[3*i-3]),coords[3*i-1]); } if(y < lim[1]) *dy += d*(y-lim[0])/(lim[1]-lim[0]); else *dy += d*(lim[2]-y)/(lim[2]-lim[1]); } } static void MovePointArc(Real *lim,int points,Real *coords, Real x,Real y,Real *dx,Real *dy) { static int i=0; Real sx,sy,ss,r,rat,d,x0,y0; if(y > lim[0] && y < lim[2]) { if(x <= coords[0]) { d = coords[1]; } else if(x >= coords[points-2]) { d = coords[points-1]; } else { i = 1; while(x > coords[3*i] && i < points/3-1) i++; sx = 0.5*(coords[3*i]-coords[3*i-3]); sy = 0.5*(coords[3*i+1]-coords[3*i-2]); r = coords[3*i-1]; ss = sx*sx+sy*sy; rat = sqrt(1.0/ss-1.0/(r*r))*r; x0 = coords[3*i-3] + sx - sy * rat; y0 = coords[3*i-2] + sy + sx * rat; d = y0-sqrt(r*r-(x-x0)*(x-x0))*r/fabs(r); } if(y < lim[1]) *dy += d*(y-lim[0])/(lim[1]-lim[0]); else *dy += d*(lim[2]-y)/(lim[2]-lim[1]); } } void CreateKnots(struct GridType *grid,struct CellType *cell, struct FemType *data,int noknots,int info) /* Saves information concerning the knots to a special structure to avoid repetitous calculations. This should be used unless there is a sevire lack of memory. GridType includes only rectangular 2D elements. */ { Real globalcoord[DIM*MAXNODESD2]; Real maplim[3*MAXMAPPINGS]; int material,nonodes,elemind,elemtype; int mode,level,maplevel,dim; int celli,cellj,i,j,k,l,ind[MAXNODESD2]; Real x,y,z,dx,dy,dz,size,minsize,maxsize; InitializeKnots(data); dim = data->dim = grid->dimension; nonodes = grid->nonodes; data->maxnodes = grid->nonodes; data->nocells = grid->nocells; data->noelements = grid->noelements; data->coordsystem = grid->coordsystem; data->numbering = grid->numbering; data->indexwidth = grid->maxwidth; data->noknots = MAX(noknots,grid->noknots); data->noboundaries = grid->noboundaries; for(i=0;iboundint[i] = grid->boundint[i]; data->boundext[i] = grid->boundext[i]; data->boundsolid[i] = grid->boundsolid[i]; data->boundtype[i] = grid->boundtype[i]; } AllocateKnots(data); minsize = 1.0e20; if(dim == 1) elemtype = grid->nonodes + 200; else elemtype = grid->nonodes + 400; for(i=1;i<=data->noelements;i++) data->elementtypes[i] = elemtype; /* This numbers the elements the same way the knots are numbered. */ for(cellj=1;cellj<= grid->ycells ;cellj++) { /* cells direction up */ for(j=1; j<=grid->yelems[cellj]; j++) /* lines inside cells */ for(celli=1;celli<= grid->xcells; celli++) /* cells direction right */ if(k=grid->numbered[cellj][celli]) { material = cell[k].material; for(i=1; i<=grid->xelems[celli]; i++) { elemind = GetElementCoordinates(&(cell)[k],i,j,globalcoord,ind); if(data->noknots == grid->noknots) { for(l=0;ltopology[elemind][l] = ind[l]; data->x[ind[l]] = globalcoord[l]; data->y[ind[l]] = globalcoord[l+nonodes]; } data->material[elemind] = material; } } } } /* Map the knots as defined in structures grid */ for(k=0;kmappings;k++) { j = grid->mappingline[k]; if(grid->mappingtype[k] > 0) maplim[3*k+1] = grid->y[j]; else if(grid->mappingtype[k] < 0) maplim[3*k+1] = grid->x[j]; else continue; maplim[3*k] = maplim[3*k+1] - grid->mappinglimits[2*k]; maplim[3*k+2] = maplim[3*k+1] + grid->mappinglimits[2*k+1]; } if(grid->mappings) for(level=0;level<10;level++) { maplevel = FALSE; for(k=0;kmappings;k++) if(abs(grid->mappingtype[k]/10) == level) maplevel = TRUE; if(maplevel == FALSE) continue; if(level >= 5) data->dim = 3; for(i=1;i<=data->noknots;i++) { x = data->x[i]; y = data->y[i]; dx = 0.0; dy = 0.0; dz = 0.0; for(k=0;kmappings;k++) { mode = grid->mappingtype[k]%10; switch (mode) { case 1: MovePointLinear(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], x,y,&dx,&dy); break; case 2: MovePointPower(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], x,y,&dx,&dy); break; case 3: MovePointArc(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], x,y,&dx,&dy); break; case 4: MovePointCircle(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], x,y,&dx,&dy); break; case 5: MovePointSinus(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], x,y,&dx,&dy); break; case 6: MovePointPowerSeries(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], x,y,&dx,&dy); break; case 7: MovePointPowerSeries2(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], x,y,&dx,&dy); break; case 8: MovePointAngle(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], x,y,&dx,&dz); break; case -1: MovePointLinear(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], y,x,&dy,&dx); break; case -2: MovePointPower(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], y,x,&dy,&dx); break; case -3: MovePointArc(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], y,x,&dy,&dx); break; case -4: MovePointCircle(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], y,x,&dy,&dx); break; case -5: MovePointSinus(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], y,x,&dy,&dx); break; case -6: MovePointPowerSeries(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], y,x,&dy,&dx); break; case -7: MovePointPowerSeries2(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], y,x,&dy,&dx); break; case -8: MovePointAngle(&maplim[3*k],grid->mappingpoints[k],grid->mappingparams[k], y,x,&dy,&dz); break; } } if(mode == 8 || mode == -8) { data->x[i] += dx; data->y[i] += dy; data->z[i] += dz; } else if(level >= 5) { data->z[i] += dx + dy; } else { data->x[i] += dx; data->y[i] += dy; } } } minsize = 1.0e20; maxsize = 0.0; for(i=1;i<=data->noelements;i++) { GetElementInfo(i,data,globalcoord,ind,&material); dx = globalcoord[0]-globalcoord[1]; dy = globalcoord[nonodes]-globalcoord[nonodes+1]; size = dx*dx+dy*dy; if(size < minsize) minsize = size; if(size > maxsize) maxsize = size; dx = globalcoord[0]-globalcoord[nonodes-1]; dy = globalcoord[nonodes]-globalcoord[2*nonodes-1]; size = dx*dx+dy*dy; if(size < minsize) minsize = size; if(size > maxsize) maxsize = size; } data->maxsize = sqrt(maxsize); data->minsize = sqrt(minsize); if(info) printf("Maximum elementsize is %.3le and minimum %.3le.\n", data->maxsize,data->minsize); } int CreateVariable(struct FemType *data,int variable,int unknowns, Real value,char *dofname,int eorder) /* Create variables for the given data structure */ { int i,info=FALSE; int timesteps; if(variable == 0) return(0); if(data->created == FALSE) { if(info) printf("CreateVariable: Knots must first be created!\n"); return(1); } timesteps = data->timesteps; if(timesteps < 1) timesteps = 1; if(data->edofs[variable] == 0) { data->variables += 1; data->edofs[variable] = unknowns; data->alldofs[variable] = unknowns * data->noknots; data->bandwidth[variable] = unknowns * data->indexwidth; data->dofs[variable] = Rvector(1,timesteps * data->alldofs[variable]); if(info) printf("Created variable %s with %d dofs.\n", dofname,data->alldofs[variable]); for(i=1;i<=data->alldofs[variable]*timesteps;i++) data->dofs[variable][i] = value; data->iterdofs[variable] = 1; } else if (data->edofs[variable] == unknowns) { if(info) printf("CreateVariable: Variable %d exists with correct number of dofs!\n", variable); } else { if(info) printf("CreateVariable: Variable %d exists with wrong number of dofs!\n", variable); return(2); } if(eorder) { if (data->eorder[variable] == FALSE) { data->eorder[variable] = TRUE; data->order[variable] = Ivector(1,data->alldofs[variable]); for(i=1;i<=data->alldofs[variable];i++) data->order[variable][i] = i; } if(info) printf("Created index for variable %s.\n",dofname); } strcpy(data->dofname[variable],dofname); return(0); } void DestroyKnots(struct FemType *data) { int i; if(!data->created) return; data->created = FALSE; for(i=0;iedofs[i] != 0) { if(data->edofs[i] > 0) free_Rvector(data->dofs[i],1,data->alldofs[i]); data->edofs[i] = 0; } free_Imatrix(data->topology,1,data->noelements,0,data->maxnodes-1); free_Ivector(data->material,1,data->noelements); free_Ivector(data->elementtypes,1,data->noelements); free_Rvector(data->x,1,data->noknots); free_Rvector(data->y,1,data->noknots); if(data->dim == 3) free_Rvector(data->z,1,data->noknots); if(data->nocorners > 0) free_Ivector(data->corners,1,2*data->nocorners); } void SideAreas(struct FemType *data,struct BoundaryType *bound) /* Calculate the sideares for later use into structure 'bound'. In 2D case the area means line length. */ { int i,ind[MAXNODESD1],sideelemtype; Real r1,r2,z1,z2; if(data->mapgeo == bound->maparea) return; bound->totalarea = 0.; if(data->dim != 2) { return; } for(i=1; i<=bound->nosides; i++) { GetElementSide(bound->parent[i],bound->side[i],bound->normal[i], data,ind,&sideelemtype); r1 = data->x[ind[0]]; r2 = data->x[ind[1]]; z1 = data->y[ind[0]]; z2 = data->y[ind[1]]; if(bound->coordsystem == COORD_CART2) bound->areas[i] = sqrt( (z1-z2)*(z1-z2) + (r1-r2)*(r1-r2) ); else if(bound->coordsystem == COORD_AXIS) bound->areas[i] = FM_PI * (r1+r2) * sqrt( (z1-z2)*(z1-z2) + (r1-r2)*(r1-r2) ); else if(bound->coordsystem == COORD_POLAR) bound->areas[i] = fabs((z2-z1)*(r1+r2)/2.0); bound->totalarea += bound->areas[i]; } bound->maparea = data->mapgeo; } int CreateBoundary(struct CellType *cell,struct FemType *data, struct BoundaryType *bound,int material1,int material2, int solidmat,int boundarytype) /* This subroutine makes a boundary which includes all sides that separate two materials that fullfill the conditions in the function call. If both materials are positive only the sides for which both of the materials coinside are accepted. In other cases the negative argument tells which conditions the positive argument should fullfill. Note that on a boundary where knots are created only for the other material, this material should be the latter one in the function call (material). The physical properties (emissivity) of the materials are taken from the one given by the flag 'solidmat'. */ { int i,side,more,elem,elemind[2],nosides,no,times; int sidemat,thismat,size,setpoint,dimsides,cellside; if(data->dim == 1) dimsides = 2; else dimsides = 4; if(bound->created == TRUE) { printf("CreateBoundary: You tried to recreate the boundary!\n"); return(1); } if(!data->created) { printf("CreateBoundary: You tried to create a boundary before the knots were made."); return(2); } if(material1 < 0 && material2 < 0) { printf("CreateBoundary: the material arguments are both negative"); return(3); } times = 0; bound->created = FALSE; bound->nosides = 0; if(solidmat >= 2) solidmat -= 2; startpoint: /* Go through all elements which have a boundary with the given material, but are not themself of that material. First only calculate their amount, then allocate space and tabulate them. */ nosides = 0; for(no=1; no <= data->nocells; no++) for(side=0; side < dimsides; side++) { if(data->dim == 1) cellside = 3-2*side; else cellside = side; setpoint = FALSE; sidemat = cell[no].boundary[cellside]; thismat = cell[no].material; /* The free boundary conditions are not allowed if the negative keywords are used. */ #if 0 /* This has been eliminated since it just causes confusion */ if(sidemat == MAT_ORIGO && material1 != MAT_ORIGO && (material1 <= MAT_BIGGER || material2 <= MAT_BIGGER)) continue; #endif /* Either material must be the one defined. */ if( material1 >= 0 && material1 != sidemat) continue; if( material2 >= 0 && material2 != thismat) continue; #if 0 printf("mat=[%d %d] sidemat=%d thismat=%d side=%d\n", material1,material2,sidemat,thismat,side); #endif if( material2 == -((side+2)%4+1) && sidemat == material1 && sidemat != thismat) setpoint = TRUE; if( material1 == -(side+1) && thismat == material2 && sidemat != thismat) setpoint = TRUE; if( material1 == MAT_BIGGER && sidemat > material2 ) setpoint = TRUE; if( material1 == MAT_SMALLER && sidemat < material2 ) setpoint = TRUE; if( material1 == MAT_ANYTHING && sidemat != material2 ) setpoint = TRUE; if( material2 == MAT_BIGGER && thismat > material1 ) setpoint = TRUE; if( material2 == MAT_SMALLER && thismat < material1 ) setpoint = TRUE; if( material2 == MAT_ANYTHING && thismat != material1 ) setpoint = TRUE; if( sidemat == material1 && thismat == material2 ) setpoint = TRUE; if(setpoint == TRUE) { #if 0 printf("going through boundary %d vs. %d in cell %d\n",material1,material2,no); #endif elem = 0; do { elem++; nosides++; more = GetSideInfo(cell,no,side,elem,elemind); #if 0 printf("elem=%d nosides=%d no=%d side=%d elemind=%d %d\n", elem,nosides, no, side, elemind[0], elemind[1]); #endif /* In the second round the values are tabulated. */ if(times) { /* It is assumed that the material pointed by solidmat determines the surface properties. */ bound->parent[nosides] = elemind[0]; bound->parent2[nosides] = elemind[1]; bound->side[nosides] = side; bound->side2[nosides] = (side+dimsides/2)%dimsides; bound->types[nosides] = boundarytype; /* The direction of the surface normal must be included */ if(solidmat==FIRST) { bound->material[nosides] = sidemat; bound->normal[nosides] = 1; } if(solidmat==SECOND){ bound->material[nosides] = thismat; bound->normal[nosides] = -1; } } } while(more); } } if(nosides == 0) { printf("No boundary between materials %d and %d exists.\n", material1,material2); return(0); } if(times == 0) { times++; /* Allocate space. This has sometimes led to strange errors. The allocation takes place only in the first loop. */ bound->created = TRUE; bound->nosides = size = nosides; bound->coordsystem = data->coordsystem; bound->fixedpoints = 1; bound->open = FALSE; bound->maparea = 0; bound->mapvf = 0; bound->types = Ivector(1,nosides); bound->areas = Rvector(1,nosides); bound->side = Ivector(1,nosides); bound->side2 = Ivector(1,nosides); bound->material = Ivector(1,nosides); bound->parent = Ivector(1,nosides); bound->parent2 = Ivector(1,nosides); bound->normal = Ivector(1,nosides); bound->vfcreated = FALSE; bound->gfcreated = FALSE; bound->echain = FALSE; bound->ediscont = FALSE; for(i=0;ievars[i] = FALSE; goto startpoint; } /* Calculate the areas of the side elements. */ SideAreas(data,bound); printf("%d element sides between materials %d and %d were located to type %d.\n", nosides,material1,material2,boundarytype); return(0); } int AllocateBoundary(struct BoundaryType *bound,int size) { int i; if(bound->created == TRUE) { printf("AllocateBoundary: You tried to recreate the boundary!\n"); return(1); } bound->created = TRUE; bound->nosides = size; bound->fixedpoints = 1; bound->open = FALSE; bound->maparea = 0; bound->mapvf = 0; bound->vfcreated = FALSE; bound->gfcreated = FALSE; bound->echain = FALSE; bound->ediscont = FALSE; bound->areas = Rvector(1,size); bound->material = Ivector(1,size); bound->side = Ivector(1,size); bound->side2 = Ivector(1,size); bound->parent = Ivector(1,size); bound->parent2 = Ivector(1,size); bound->types = Ivector(1,size); bound->normal = Ivector(1,size); for(i=1;i<=size;i++) { bound->areas[i] = 0.0; bound->material[i] = 0; bound->side[i] = 0; bound->side2[i] = 0; bound->parent[i] = 0; bound->parent2[i] = 0; bound->types[i] = 0; bound->normal[i] = 1; } for(i=0;ievars[i] = FALSE; return(0); } int CreateBoundaryChain(struct FemType *data,struct BoundaryType *bound,int info) { int i,j,sideind[MAXNODESD1],sideind2[MAXNODESD1]; int size,indfirst,indsecond,side,setchain,sideelemtype; if(bound->created == FALSE) { if(info) printf("CreateBoundaryChain: boundary not created!\n"); return(1); } if(bound->nosides < 2) { if(info) printf("CreateBoundaryChain: there must be at least 2 sides!\n"); return(2); } if(bound->echain == TRUE) { if(0) printf("CreateBoundaryChain: chain already exists!\n"); return(3); } setchain = FALSE; restart: /* First determine the way that the chain starts. */ GetElementSide(bound->parent[1],bound->side[1],bound->normal[1],data,sideind,&sideelemtype); GetElementSide(bound->parent[2],bound->side[2],bound->normal[2],data,sideind2,&sideelemtype); if(sideind[1] == sideind2[0]) { indfirst = sideind[0]; indsecond = sideind[1]; } else if(sideind[0] == sideind2[1]) { indfirst = sideind[1]; indsecond = sideind[0]; } else printf("CreateBoundaryChain: Impossible case!\n"); if(setchain) { bound->chain[0] = indfirst; bound->chain[1] = indsecond; } size = 1; side = 1; /* The chain is followed in positive direction while succesful, if not the direction is changed and so on... */ for(;;) { for(i=side+1;i<=bound->nosides;i++) { GetElementSide(bound->parent[i],bound->side[i],bound->normal[i],data,sideind,&sideelemtype); if(sideind[0] == indsecond) { indsecond = sideind[1]; goto jump; } else if(sideind[1] == indsecond) { indsecond = sideind[0]; goto jump; } } for(i=side-1;i>=1;i--) { GetElementSide(bound->parent[i],bound->side[i],bound->normal[i],data,sideind,&sideelemtype); if(sideind[0] == indsecond) { indsecond = sideind[1]; goto jump; } else if(sideind[1] == indsecond) { indsecond = sideind[0]; goto jump; } } goto theend; jump: size++; if(setchain) bound->chain[size] = indsecond; side = i; if(indsecond == indfirst) goto theend; } theend: if(setchain == FALSE) { if(size != bound->nosides) printf("CreateBoundaryChain: the boundary is not continuous (%d vs. %d)\n", size,bound->nosides); else printf("CreateBoundaryChain: the boundary is continuous with %d nodes.\n",size+1); bound->chain = Ivector(0,size); bound->echain = TRUE; bound->chainsize = size; setchain = TRUE; goto restart; } return(0); } int DestroyBoundary(struct BoundaryType *bound) /* Destroys boundaries of various types. */ { int i,nosides; if(!bound->created) { printf("DestroyBoundary: boundary does not exist."); return(1); } nosides = bound->nosides; if(!nosides) { bound->created = FALSE; return(2); } free_Rvector(bound->areas,1,nosides); free_Ivector(bound->material,1,nosides); free_Ivector(bound->side,1,nosides); free_Ivector(bound->side2,1,nosides); free_Ivector(bound->parent,1,nosides); free_Ivector(bound->parent2,1,nosides); if(bound->vfcreated) free_Rmatrix(bound->vf,1,nosides,1,nosides); if(bound->gfcreated) free_Rmatrix(bound->gf,1,nosides,1,nosides); for(i=0;ievars[i]) { bound->evars[i] = 0; } bound->nosides = 0; bound->created = FALSE; #if DEBUG printf("%d element sides were destroyed.\n",nosides); #endif return(0); } int CreateBoundaries(struct CellType *cell,struct FemType *data, struct BoundaryType *boundaries) { int i,j; j = 0; if(data->noboundaries > 0) for(i=0;inoboundaries;i++) { while(boundaries[j].created) { j++; if(j >= MAXBOUNDARIES) { printf("CreateBoundaries: too many boundaries %d\n",j); return(1); } } CreateBoundary(cell,data,&boundaries[j], data->boundext[i],data->boundint[i], data->boundsolid[i],data->boundtype[i]); } return(0); } int CreatePoints(struct CellType *cell,struct FemType *data, struct BoundaryType *bound, int param1,int param2,int pointmode,int pointtype) { int size,i,no,corner,times,elem,node,sideelemtype; bound->created = FALSE; bound->nosides = 0; times = 0; omstart: i = 0; /* Create nodes that are devided by the two materials specified */ if(pointmode == 4) { for(no=1; no <= data->nocells; no++) if(cell[no].material == param2) { for(corner=0; corner < 4; corner++) if(cell[no].boundary[4+corner] == param1) { i++; if(times) { bound->material[i] = param2; bound->types[i] = pointtype; bound->side[i] = 4 + corner; if(corner == BOTLEFT) elem = GetElementIndex(&cell[no],1,1); else if(corner == BOTRIGHT) elem = GetElementIndex(&cell[no],cell[no].xelem,1); else if(corner == TOPRIGHT) elem = GetElementIndex(&cell[no],cell[no].xelem,cell[no].yelem); else if(corner == TOPLEFT) elem = GetElementIndex(&cell[no],1,cell[no].yelem); bound->parent[i] = elem; } } } } if(pointmode == 5) { corner = -1; for(no=1; no <= data->nocells && corner <0; no++) { if(cell[no].xind-1 == param1 && cell[no].yind-1 == param2) corner = BOTLEFT; else if(cell[no].xind == param1 && cell[no].yind-1 == param2) corner = BOTRIGHT; else if(cell[no].xind == param1 && cell[no].yind == param2) corner = TOPRIGHT; else if(cell[no].xind-1 == param1 && cell[no].yind == param2) corner = TOPLEFT; } if(corner >= 0) { i++; no--; } if(times) { bound->types[i] = pointtype; bound->side[i] = 4 + corner; if(corner == BOTLEFT) elem = GetElementIndex(&cell[no],1,1); else if(corner == BOTRIGHT) elem = GetElementIndex(&cell[no],cell[no].xelem,1); else if(corner == TOPRIGHT) elem = GetElementIndex(&cell[no],cell[no].xelem,cell[no].yelem); else if(corner == TOPLEFT) elem = GetElementIndex(&cell[no],1,cell[no].yelem); bound->parent[i] = elem; node = data->topology[elem][corner]; printf("Found node %d at (%.3lg, %.3lg)\n",node,data->x[node],data->y[node]); } } size = i; if(times == 0 && size > 0) { AllocateBoundary(bound,size); times = 1; goto omstart; } printf("Created %d new points of type %d in the corner of materials %d and %d.\n", size,pointtype,param1,param2); return(FALSE); } static int CreateNewNodes(struct FemType *data,int *order,int material,int new) { int i,j,k,l,lmax,ind; int newsize,noknots,nonodes; int *neworder; Real *newx,*newy,*newz,*newdofs[MAXDOFS]; noknots = data->noknots; printf("Creating %d new nodes for discoutinuous boundary.\n",new); /* Allocate for the new nodes */ newsize = noknots+new; newx = Rvector(1,newsize); newy = Rvector(1,newsize); if(data->dim == 3) newz = Rvector(1,newsize); neworder = Ivector(1,newsize); for(j=1;jedofs[j]) newdofs[j] = Rvector(1,data->edofs[j]*newsize); /* Set the new coordinates and dofs */ j = 0; for(i=1;i<=noknots;i++) { j++; neworder[j] = i; newx[j] = data->x[i]; newy[j] = data->y[i]; if(data->dim == 3) newz[j] = data->z[i]; for(k=1;kedofs[k]) for(l=1;l<=lmax;l++) newdofs[k][lmax*(j-1)+l] = data->dofs[k][lmax*(i-1)+l]; } if(order[i] < 0) { j++; neworder[j] = -i; newx[j] = data->x[i]; newy[j] = data->y[i]; if(data->dim == 3) newz[j] = data->z[i]; for(k=1;kedofs[k]) for(l=1;l<=lmax;l++) newdofs[k][lmax*(j-1)+l] = data->dofs[k][lmax*(i-1)+l]; } } } /* Find the old index corresponding to the new one. */ for(i=1;i<=newsize;i++) if(neworder[i] > 0) { if(order[neworder[i]] < 0) order[neworder[i]] = -i; else order[neworder[i]] = i; } /* Set the new element topology */ for(i=1;i<=data->noelements;i++) { nonodes = data->elementtypes[i]%100; for(j=0;jtopology[i][j]; if(data->material[i] == material && order[ind] < 0) data->topology[i][j] = abs(order[ind])+1; else data->topology[i][j] = abs(order[ind]); } } /* Destroy old vectors and set the pointers to the new vectors. */ free_Rvector(data->x,1,noknots); free_Rvector(data->y,1,noknots); if(data->dim == 3) free_Rvector(data->z,1,noknots); for(k=1;kedofs[k]) free_Rvector(data->dofs[k],1,noknots); data->noknots = newsize; data->x = newx; data->y = newy; if(data->dim == 3) data->z = newz; for(k=1;kedofs[k]) { data->dofs[k] = newdofs[k]; data->alldofs[k] = data->edofs[k] * data->noknots; } } return(0); } int SetDiscontinuousBoundary(struct FemType *data,struct BoundaryType *bound, int boundtype,int endnodes,int info) /* Create secondary points for a given boundary. This feature is handy when one wants to solve problems with discontinous field variables. */ { int i,j,k,bc,ind,sideind[MAXNODESD1]; int side,parent,new,doublesides,maxtype,newbc; int newsuccess,noelements,nonodes,sideelemtype,sidenodes,disconttype; int *order; int mat1,mat2,par1,par2,mat1old,mat2old,material; static int hitsexist=FALSE,hitslength,*hits; if(boundtype < 0) { newbc = TRUE; boundtype = -boundtype; } else { newbc = FALSE; } mat1old = mat2old = 0; doublesides = 0; for(bc=0;bcnosides) continue; for(i=1;i<=bound[bc].nosides;i++) { if(bound[bc].types[i] == boundtype) { par1 = bound[bc].parent[i]; par2 = bound[bc].parent2[i]; if(par1 && par2) { doublesides++; mat1 = data->material[par1]; mat2 = data->material[par2]; if(!mat1old) mat1old = mat1; else if(mat1old != mat1) mat1old = -mat1; if(!mat2old) mat2old = mat2; else if(mat2old != mat2) mat2old = -mat2; } } } } if(!doublesides) return(1); if(mat1old > 0) material = mat1old; else if(mat2old > 0) material = mat2old; else { printf("SetDiscontinuousBoundary: impossible to make the boundary of several materials\n"); return(2); } noelements = data->noelements; order = Ivector(1,data->noknots); for(i=1;i<=data->noknots;i++) order[i] = i; if(endnodes == 1) { if(!hitsexist) { hitslength = 1.1*data->noknots; hits = Ivector(1,hitslength); hitsexist = TRUE; } else if(hitslength <= data->noknots) { free_Ivector(hits,1,hitslength); hitslength = 1.1*data->noknots; hits = Ivector(1,hitslength); } for(i=1;i<=data->noknots;i++) hits[i] = 0; for(j=1;j<=noelements;j++) { nonodes = data->elementtypes[j]%100; for(i=0;itopology[j][i]] += 1; } } if(newbc) { maxtype = 0; for(bc=0;bc 0) { new++; order[ind] = -new; } } else if(endnodes == 0) { if(order[ind] > 0) order[ind] = 0; else if(order[ind] == 0) { new++; order[ind] = -new; } } else if(endnodes == 1) { if(order[ind] > 0) { if(hits[ind] < 4) { new++; order[ind] = -new; } else order[ind] = 0; } else if(order[ind] == 0) { new++; order[ind] = -new; } } } } if(endnodes == 0 || endnodes == 1) { for(i=1;i<=data->noknots;i++) if(order[i] == 0) order[i] = i; } } if(new == 0) return(3); newsuccess = CreateNewNodes(data,order,material,new); return(newsuccess); } int FindPeriodicBoundary(struct FemType *data,struct BoundaryType *bound, int boundary1,int boundary2,int info) /* Create periodic boundary conditions for a given boundary pair boundary1, boundary2. */ { int i,j,k,l,ind,sideind[MAXNODESD1]; int side,parent,elemtype; int minp[2],maxp[2],bounds[2],dp[2],sumsides[2]; if(bound->created == FALSE) { printf("SetDiscontinuousBoundary: The boundary does not exist!\n"); return(1); } if(!bound->nosides) return(0); bounds[0] = boundary1; bounds[1] = boundary2; minp[0] = minp[1] = data->noknots+1; maxp[0] = maxp[1] = 0; sumsides[0] = sumsides[1] = 0; for(j=0;j < MAXBOUNDARIES;j++) { if(bound->created == FALSE) continue; for(i=1; i <= bound[j].nosides; i++) { for(k=0;k<=1;k++) { if(bound[j].types[i] == bounds[k]) { sumsides[k] += 1; if(bound[j].parent[i] > maxp[k]) maxp[k] = bound[j].parent[i]; if(bound[j].parent[i] < minp[k]) minp[k] = bound[j].parent[i]; } } } } for (k=0;k<=1;k++) { dp[k] = maxp[k] - minp[k]; if(info) printf("Parents of boundary %d are on the interval [%d, %d]\n", bounds[k],minp[k],maxp[k]); } if(dp[0] != dp[1] || sumsides[0] != sumsides[1]) { printf("FindPeriodicBoundary: The simple scheme cannot handle these boundaries!\n"); printf("dp=[%d %d] sumsides=[%d %d]\n",dp[0],dp[1],sumsides[0],sumsides[1]); return(1); } for(j=0;j < MAXBOUNDARIES;j++) { if(bound->created == FALSE) continue; for(i=1; i <= bound[j].nosides; i++) { for(k=0;k<=1;k++) { if(bound[j].types[i] == bounds[k]) { parent = bound[j].parent[i]; bound[j].parent2[i] = bound[j].parent[i] - minp[k] + minp[(k+1)%2]; if(!bound[j].ediscont) { bound[j].discont = Ivector(1,bound[j].nosides); for(l=1; l <= bound[j].nosides; l++) bound[j].discont[l] = 0; bound[j].ediscont = TRUE; } bound[j].discont[i] = 2+k; elemtype = data->elementtypes[parent]; if(elemtype%100 == 4) { bound[j].side2[i] = (bound[j].side[i] + 2) % 4; } else if(elemtype%100 == 8) { if(bound[j].side[i] < 4) bound[j].side2[i] = (bound[j].side[i] + 2) % 4; if(bound[j].side[i] >= 4) bound[j].side2[i] = 4 + (5 - bound[j].side[i]); } } } } } if(info) printf("Periodic boundaries were set with a simple scheme\n"); return(2); } int SetConnectedBoundary(struct FemType *data,struct BoundaryType *bound, int bctype,int connecttype,int info) /* Create connected boundary conditions for a given bctype */ { int i,j,k,l,bc,sideelemtype,sidenodes; int sideind[MAXNODESD1]; for(bc=0;bcconnectexist) { data->connect = Ivector(1,data->noknots); for(k=1;k<=data->noknots;k++) data->connect[k] = 0; data->connectexist = TRUE; } GetElementSide(bound[bc].parent[i],bound[bc].side[i],bound[bc].normal[i], data,sideind,&sideelemtype); sidenodes = sideelemtype%100; for(j=0;jconnect[k] = connecttype; } } } return(0); } int SetDiscontinuousPoints(struct FemType *data,struct PointType *point, int material) /* Create secondary point for a given point. The variable that is used to set up the boundary must not have any previously defined Dirichlet points. */ { int i,j,ind,corner,*order; int parent,new,unknowns; int newsuccess; if(point->nopoints == FALSE) { printf("SetDiscontinuousPoint: No points exists!\n"); return(0); } order = Ivector(1,data->noknots); for(i=1;i<=data->noknots;i++) order[i] = i; /* Find the number of new nodes */ new = 0; for(i=0;inopoints;i++) { parent = point->parent[i]; corner = point->corner[i]; ind = data->topology[parent][corner]; if(order[ind] >= 0) { new++; order[ind] = -new; } } if(new == 0) return(0); newsuccess = CreateNewNodes(data,order,material,new); return(newsuccess); } int FindCorners(struct GridType *grid,struct CellType *cell, struct FemType *data,int info) /* Find the nodes in the mesh that are at material corners. These nodes are often of special interest. */ { int i,j,k,ind,cellno,elem; int allocated,nocorners; nocorners = 0; allocated = FALSE; omstart: if(nocorners > 0) { data->corners = Ivector(1,2*nocorners); data->nocorners = nocorners; allocated = TRUE; } k = 0; for(i=1;i<=grid->xcells+1;i++) for(j=1;j<=grid->ycells+1;j++) { if(grid->structure[j][i] == grid->structure[j][i-1] && grid->structure[j-1][i] == grid->structure[j-1][i-1]) continue; if(grid->structure[j][i] == grid->structure[j-1][i] && grid->structure[j][i-1] == grid->structure[j-1][i-1]) continue; /* point (i,j) must now be a corner */ if(cellno = grid->numbered[j][i]) { elem = GetElementIndex(&(cell)[cellno],1,1); ind = BOTLEFT; } else if(cellno = grid->numbered[j][i-1]) { elem = GetElementIndex(&(cell)[cellno],cell[cellno].xelem,1); ind = BOTRIGHT; } else if(cellno = grid->numbered[j-1][i]) { elem = GetElementIndex(&(cell)[cellno],1,cell[cellno].yelem); ind = TOPLEFT; } else if(cellno = grid->numbered[j-1][i-1]) { elem = GetElementIndex(&(cell)[cellno],cell[cellno].xelem,cell[cellno].yelem); ind = TOPRIGHT; } else continue; /* ind is now the index of the corner knot */ k++; if(allocated == FALSE) continue; data->corners[2*k-1] = elem; data->corners[2*k] = ind; } nocorners = k; if(nocorners == 0) return(0); if(allocated == FALSE) goto omstart; if(info) printf("Found %d material corners.\n",nocorners); return(0); } int SolutionFromMeshToMesh(struct CellType *cell1, struct GridType *grid1, struct FemType *data1, struct CellType *cell2, struct GridType *grid2, struct FemType *data2, int mapgeo,int variable,int info) /* Copies variable values from data1 to data2 for two meshes that have similar geometry, but different number of elements. Its assumed that all the elements are rectangular in shape. The subroutine may, however, be applied to nonrectangular geometries also, but then the accuracy is difficult to estimate. Note that this subroutine holds only for linear elements. */ { int xcell,ycell,i1,j1,i2,j2,no1,no2; int ind1[MAXNODESD2],ind2,k; int nonodes1,nonodes2,unknowns; int elem1,elem2,mapres,fast; Real coord1[DIM*MAXNODESD2],x2,y2,rx,ry; Real epsilon=1.0e-20; Real *vector1,*vector2; if(!data1->edofs[variable] || !data2->edofs[variable]) return(1); unknowns = data1->edofs[variable]; vector1 = data1->dofs[variable]; vector2 = data2->dofs[variable]; nonodes1 = grid1->nonodes; nonodes2 = grid2->nonodes; if(nonodes1 != 4 || nonodes2 != 4) { if(info) printf("SolutionFromMeshToMesh: algorithm defined only for 4-node elements\n"); return(2); } if(mapgeo && data2->mapgeo < data1->mapgeo) data2->mapgeo = data1->mapgeo; else mapgeo = FALSE; if(data2->iterdofs[variable] < data1->iterdofs[variable]) { mapres = TRUE; data2->iterdofs[variable] = data1->iterdofs[variable]; } else mapres = FALSE; if(!mapres && !mapgeo) { if(info) printf("SolutionFromMeshToMesh: no mapping for geometry or variable %s (%d vs. %d)!\n", data1->dofname[variable],data1->iterdofs[variable],data2->iterdofs[variable]); return(0); } if(grid1->noknots == data1->noknots && grid2->noknots == data2->noknots) fast = TRUE; else fast = FALSE; for(xcell=1;xcell<=MAXCELLS;xcell++) for(ycell=1;ycell<=MAXCELLS;ycell++) /* Go through cells that are common to both grids. */ if( (no1= grid1->numbered[ycell][xcell]) && (no2= grid2->numbered[ycell][xcell]) ) { if(0) printf("xcell=%d ycell=%d no1=%d no=%d\n",xcell,ycell,no1,no2); j1 = 1; for(j2=0; j2 <= cell2[no2].yelem; j2++) { i1 = 1; for(i2=0; i2 <= cell2[no2].xelem; i2++) { /* ind2 is the original node number for the node */ ind2 = GetKnotCoordinate(&(cell2)[no2],i2,j2,&x2,&y2); GetElementCoordinates(&(cell1)[no1],i1,j1,coord1,ind1); /* Find j1 and i1 in the rectangular mesh so that the node ind2 lies in the element */ if(coord1[TOPRIGHT+nonodes1]+epsilon < y2 && j1 < cell1[no1].yelem) do { j1++; GetElementCoordinates(&(cell1)[no1],i1,j1,coord1,ind1); } while(j1 < cell1[no1].yelem && coord1[TOPRIGHT+nonodes1] < y2); if(coord1[TOPRIGHT]+epsilon < x2 && i1 < cell1[no1].xelem) do { i1++; GetElementCoordinates(&(cell1)[no1],i1,j1,coord1,ind1); } while(i1 < cell1[no1].xelem && coord1[TOPRIGHT] < x2); if(0) printf("j2=%d i2=%d j1=%d i1=%d\n",j2,i2,j1,i1); rx = (coord1[BOTRIGHT]-x2) / (coord1[BOTRIGHT]-coord1[BOTLEFT]); ry = (coord1[TOPLEFT+nonodes1]-y2) / (coord1[TOPLEFT+nonodes1]-coord1[BOTLEFT+nonodes1]); rx = MIN(rx,1.0); rx = MAX(rx,0.0); ry = MIN(ry,1.0); ry = MAX(ry,0.0); /* Find the new indices. If there are no discontinuous boundaries they are the same as the old ones. */ if(!fast) { elem1=GetElementIndex(&(cell1)[no1],i1,j1); for(k=0;k<4;k++) ind1[k] = data1->topology[elem1][k]; if(0) printf("elem1=%d ind1=[%d %d %d %d]\n",elem1,ind1[0],ind1[1],ind1[2],ind1[3]); if(i2>0 && j2>0) { elem2 = GetElementIndex(&(cell2)[no2],i2,j2); ind2 = data2->topology[elem2][TOPRIGHT]; } else if(i2==0 && j2>0) { elem2 = GetElementIndex(&(cell2)[no2],i2+1,j2); ind2 = data2->topology[elem2][TOPLEFT]; } else if(i2>0 && j2==0) { elem2 = GetElementIndex(&(cell2)[no2],i2,j2+1); ind2 = data2->topology[elem2][BOTRIGHT]; } else { elem2 = GetElementIndex(&(cell2)[no2],i2+1,j2+1); ind2 = data2->topology[elem2][BOTLEFT]; } } if(mapgeo) { data2->x[ind2] = data1->x[ind1[BOTLEFT]] * rx * ry + data1->x[ind1[BOTRIGHT]] * (1.-rx) * ry + data1->x[ind1[TOPLEFT]] * rx * (1.-ry) + data1->x[ind1[TOPRIGHT]] * (1.-rx) * (1.-ry); data2->y[ind2] = data1->y[ind1[BOTLEFT]] * rx * ry + data1->y[ind1[BOTRIGHT]] * (1.-rx) * ry + data1->y[ind1[TOPLEFT]] * rx * (1.-ry) + data1->y[ind1[TOPRIGHT]] * (1.-rx) * (1.-ry); } if(mapres) for(k=1;k<=unknowns;k++) { vector2[unknowns*(ind2-1)+k] = vector1[unknowns*(ind1[BOTLEFT]-1)+k] * rx * ry + vector1[unknowns*(ind1[BOTRIGHT]-1)+k]* (1.-rx) * ry + vector1[unknowns*(ind1[TOPLEFT]-1)+k] * rx * (1.-ry) + vector1[unknowns*(ind1[TOPRIGHT]-1)+k]* (1.-rx) * (1.-ry); } } } } if(info) { if(mapgeo) printf("Geometry was mapped from one mesh to another!\n"); if(mapres) printf("Results of %s was mapped from one mesh to another!\n", data1->dofname[variable]); } return(0); } int ElementsToTriangles(struct FemType *data,struct BoundaryType *bound, int info) /* Make triangles out of rectangular elements */ { int i,j,k,l,side,elem,i1,i2,isum,foundside,sideelemtype; int noelements,elementtype,triangletype,triangles,noknots,nonodes; int **newtopo,*newmaterial,*newelementtypes,newnodes,*indx,*needed; int sideind[MAXNODESD1], sideind2[MAXNODESD1]; Real dx1,dx2,dy1,dy2,ds1,ds2; struct FemType data2; noelements = data->noelements; noknots = data->noknots; elementtype = data->elementtypes[1]; for(i=1;i<=data->noelements;i++) { if(data->elementtypes[i] != elementtype) printf("ElementsToTriangles: Implemented only for constant elementtypes (%d,%d).\n", data->elementtypes[i],elementtype); } switch(elementtype) { case 404: triangletype = 303; newnodes = 3; triangles = 2; break; case 405: triangletype = 303; newnodes = 3; triangles = 4; break; case 409: triangletype = 306; newnodes = 6; triangles = 2; break; case 416: triangletype = 310; newnodes = 10; triangles = 2; break; default: printf("ElementsToTriangles: not implemented for elementtype %d\n",elementtype); return(1); } needed = Ivector(1,noknots); for(i=1;i<=noknots;i++) needed[i] = 0; nonodes = elementtype / 100; for(i=1;i<=noelements;i++) { for(j=0;jtopology[i][j]] += 1; } /* First divide the elements along the shorter diameter */ newtopo = Imatrix(1,triangles*noelements,0,newnodes-1); newmaterial = Ivector(1,triangles*noelements); newelementtypes = Ivector(1,triangles*noelements); data2 = *data; data2.topology = newtopo; data2.material = newmaterial; data2.elementtypes = newelementtypes; for(i=1;i<=noelements;i++) { newelementtypes[2*i-1] = triangletype; newelementtypes[2*i] = triangletype; /* compute the diagonals in order to make division along the shorter one */ dx1 = data->x[data->topology[i][0]] - data->x[data->topology[i][2]]; dy1 = data->y[data->topology[i][0]] - data->y[data->topology[i][2]]; dx2 = data->x[data->topology[i][1]] - data->x[data->topology[i][3]]; dy2 = data->y[data->topology[i][1]] - data->y[data->topology[i][3]]; ds1 = dx1*dx1+dy1*dy1; ds2 = dx2*dx2+dy2*dy2; /* In case of corner nodes favor division where corner is split */ if(needed[data->topology[i][0]] <= 2 && (needed[data->topology[i][2]]) <= 2) ds1 *= 2; if(needed[data->topology[i][1]] <= 2 && (needed[data->topology[i][3]]) <= 2) ds2 *= 2; if(elementtype == 404 || elementtype == 409 || elementtype == 416) { if(ds1 > ds2) { newtopo[2*i-1][0] = data->topology[i][0]; newtopo[2*i-1][1] = data->topology[i][1]; newtopo[2*i-1][2] = data->topology[i][3]; newmaterial[2*i-1]= data->material[i]; newtopo[2*i][0] = data->topology[i][2]; newtopo[2*i][1] = data->topology[i][3]; newtopo[2*i][2] = data->topology[i][1]; newmaterial[2*i] = data->material[i]; } else { newtopo[2*i-1][0] = data->topology[i][1]; newtopo[2*i-1][1] = data->topology[i][2]; newtopo[2*i-1][2] = data->topology[i][0]; newmaterial[2*i-1]= data->material[i]; newtopo[2*i][0] = data->topology[i][3]; newtopo[2*i][1] = data->topology[i][0]; newtopo[2*i][2] = data->topology[i][2]; newmaterial[2*i] = data->material[i]; } } if(elementtype == 409) { if(ds1 > ds2) { newtopo[2*i-1][3] = data->topology[i][4]; newtopo[2*i-1][4] = data->topology[i][8]; newtopo[2*i-1][5] = data->topology[i][7]; newtopo[2*i][3] = data->topology[i][6]; newtopo[2*i][4] = data->topology[i][8]; newtopo[2*i][5] = data->topology[i][5]; } else { newtopo[2*i-1][3] = data->topology[i][5]; newtopo[2*i-1][4] = data->topology[i][8]; newtopo[2*i-1][5] = data->topology[i][4]; newtopo[2*i][3] = data->topology[i][7]; newtopo[2*i][4] = data->topology[i][8]; newtopo[2*i][5] = data->topology[i][6]; } } if(elementtype == 416) { if(ds1 > ds2) { newtopo[2*i-1][3] = data->topology[i][4]; newtopo[2*i-1][4] = data->topology[i][5]; newtopo[2*i-1][5] = data->topology[i][13]; newtopo[2*i-1][6] = data->topology[i][15]; newtopo[2*i-1][7] = data->topology[i][10]; newtopo[2*i-1][8] = data->topology[i][11]; newtopo[2*i-1][9] = data->topology[i][12]; newtopo[2*i][3] = data->topology[i][8]; newtopo[2*i][4] = data->topology[i][9]; newtopo[2*i][5] = data->topology[i][15]; newtopo[2*i][6] = data->topology[i][13]; newtopo[2*i][7] = data->topology[i][6]; newtopo[2*i][8] = data->topology[i][7]; newtopo[2*i][9] = data->topology[i][14]; } else { newtopo[2*i-1][3] = data->topology[i][6]; newtopo[2*i-1][4] = data->topology[i][7]; newtopo[2*i-1][5] = data->topology[i][14]; newtopo[2*i-1][6] = data->topology[i][12]; newtopo[2*i-1][7] = data->topology[i][4]; newtopo[2*i-1][8] = data->topology[i][5]; newtopo[2*i-1][9] = data->topology[i][13]; newtopo[2*i][3] = data->topology[i][10]; newtopo[2*i][4] = data->topology[i][11]; newtopo[2*i][5] = data->topology[i][12]; newtopo[2*i][6] = data->topology[i][14]; newtopo[2*i][7] = data->topology[i][8]; newtopo[2*i][8] = data->topology[i][9]; newtopo[2*i][9] = data->topology[i][15]; } } if(elementtype == 405) { newtopo[4*i-3][0] = data->topology[i][0]; newtopo[4*i-3][1] = data->topology[i][1]; newtopo[4*i-3][2] = data->topology[i][4]; newmaterial[4*i-3]= data->material[i]; newtopo[4*i-2][0] = data->topology[i][1]; newtopo[4*i-2][1] = data->topology[i][2]; newtopo[4*i-2][2] = data->topology[i][4]; newmaterial[4*i-2]= data->material[i]; newtopo[4*i-1][0] = data->topology[i][2]; newtopo[4*i-1][1] = data->topology[i][3]; newtopo[4*i-1][2] = data->topology[i][4]; newmaterial[4*i-1]= data->material[i]; newtopo[4*i][0] = data->topology[i][3]; newtopo[4*i][1] = data->topology[i][0]; newtopo[4*i][2] = data->topology[i][4]; newmaterial[4*i] = data->material[i]; } } /* Then make the corresponding mapping for the BCs. This is done in a brute-force way where all the possible new elements are checked. */ for(j=0;j < MAXBOUNDARIES;j++) { if(!bound[j].created) continue; for(i=1; i <= bound[j].nosides; i++) { for(l=1;l<=2;l++) { if(l==1) k = bound[j].parent[i]; else k = bound[j].parent2[i]; if(k == 0) continue; if(l == 1) GetElementSide(bound[j].parent[i],bound[j].side[i],bound[j].normal[i], data,sideind,&sideelemtype); else GetElementSide(bound[j].parent2[i],bound[j].side2[i],bound[j].normal[i], data,sideind,&sideelemtype); if(sideelemtype/100 != 2) { printf("ElementsToTriangles: implemeted only for BCs 202 and 203\n"); continue; } /* Test for all possible elements that could be parents */ for(elem = triangles*(k-1)+1;elem <= k*triangles;elem++) { isum = 0; for(i1=0;i1<3;i1++) { if(newtopo[elem][i1] == sideind[0]) isum++; if(newtopo[elem][i1] == sideind[1]) isum++; } if(isum != 2) continue; for(side=0;side<3;side++) { GetElementSide(elem,side,bound[j].normal[i], &data2,sideind2,&sideelemtype); isum = 0; for(i1=0;i1<2;i1++) { if(sideind2[i1] == sideind[0]) isum++; if(sideind2[i1] == sideind[1]) isum++; } if(isum == 2) goto nextparent; } } nextparent: if(isum == 2) { if(l == 1) { bound[j].parent[i] = elem; bound[j].side[i] = side; } if(l == 2) { bound[j].parent2[i] = elem; bound[j].side2[i] = side; } } else { printf("Failed to find parent for side %d of %d (parent %d)\n",i,j,k); } } } } free_Imatrix(data->topology,1,noelements,0,data->maxnodes-1); free_Ivector(data->material,1,noelements); free_Ivector(data->elementtypes,1,noelements); free_Ivector(needed,1,noknots); data->topology = newtopo; data->elementtypes = newelementtypes; data->material = newmaterial; data->noelements *= triangles; data->maxnodes = newnodes; if(info) printf("The mesh was reconstructed from %d triangles\n",triangles*noelements); return(0); } int PolarCoordinates(struct FemType *data,Real rad,int info) { int i,j,k,l; Real x,y,dx,dy,fii,zet,dr; for(i=1;i<=data->noknots;i++) { zet = data->x[i]; fii = FM_PI/180.0 * data->y[i]; dr = data->z[i]; data->z[i] = zet; data->x[i] = (rad+dr) * cos(fii); data->y[i] = (rad+dr) * sin(fii); } if(info) printf("Making coordinate transformation from polar to cartesian\n"); return(0); } int CylinderCoordinates(struct FemType *data,int info) { int i,j,k,l; Real x,y,rad,fii; for(i=1;i<=data->noknots;i++) { rad = data->x[i]; fii = FM_PI/180.0 * data->y[i]; data->x[i] = rad * cos(fii); data->y[i] = rad * sin(fii); } if(info) printf("Making coordinate transformation from cylindrical to cartesian\n"); return(0); } int UniteMeshes(struct FemType *data1,struct FemType *data2, struct BoundaryType *bound1,struct BoundaryType *bound2, int info) /* Unites two meshes for one larger mesh */ { int i,j,k,l; int noelements,noknots,nonodes; int **newtopo,*newmaterial,*newelementtypes,maxnodes; Real *newx,*newy,*newz; noknots = data1->noknots + data2->noknots; noelements = data1->noelements + data2->noelements; maxnodes = MAX(data1->maxnodes,data2->maxnodes); if(data2->dim > data1->dim) data1->dim = data2->dim; if(0) printf("Uniting two meshes to %d nodes and %d elements.\n",noknots,noelements); for(j=0;j < MAXBOUNDARIES;j++) { if(!bound2[j].created) continue; for(k=j;k < MAXBOUNDARIES;k++) if(!bound1[k].created) break; #if 0 printf("k=%d j=%d\n",k,j); #endif bound1[k].created = bound2[j].created; bound1[k].nosides = bound2[j].nosides; bound1[k].coordsystem = bound2[j].coordsystem; bound1[k].side = bound2[j].side; bound1[k].side2 = bound2[j].side2; bound1[k].parent = bound2[j].parent; bound1[k].parent2 = bound2[j].parent2; bound1[k].areas = bound2[j].areas; bound1[k].material = bound2[j].material; bound1[k].vfcreated = bound2[j].vfcreated; bound1[k].gfcreated = bound2[j].gfcreated; bound1[k].echain = bound2[j].echain; bound1[k].types = bound2[j].types; bound1[k].normal = bound2[j].normal; for(i=1; i <= bound1[k].nosides; i++) { bound1[k].parent[i] += data1->noelements; if(bound1[k].parent2[i]) bound1[k].parent2[i] += data1->noelements; } } data1->maxnodes = maxnodes; newtopo = Imatrix(1,noelements,0,maxnodes-1); newmaterial = Ivector(1,noelements); newelementtypes = Ivector(1,noelements); newx = Rvector(1,noknots); newy = Rvector(1,noknots); if(data1->dim == 3) newz = Rvector(1,noknots); for(i=1;i<=data1->noknots;i++) { newx[i] = data1->x[i]; newy[i] = data1->y[i]; if(data1->dim == 3) newz[i] = data1->z[i]; } for(i=1;i<=data2->noknots;i++) { newx[i+data1->noknots] = data2->x[i]; newy[i+data1->noknots] = data2->y[i]; if(data1->dim == 3) newz[i+data1->noknots] = data2->z[i]; } for(i=1;i<=data1->noelements;i++) { newmaterial[i] = data1->material[i]; newelementtypes[i] = data1->elementtypes[i]; nonodes = newelementtypes[i]%100; for(j=0;jtopology[i][j]; } for(i=1;i<=data2->noelements;i++) { newmaterial[i+data1->noelements] = data2->material[i]; newelementtypes[i+data1->noelements] = data2->elementtypes[i]; nonodes = newelementtypes[i+data1->noelements]%100; for(j=0;jnoelements][j] = data2->topology[i][j] + data1->noknots; } free_Imatrix(data1->topology,1,data1->noelements,0,data1->maxnodes-1); free_Ivector(data1->material,1,data1->noelements); free_Rvector(data1->x,1,data1->noknots); free_Rvector(data1->y,1,data1->noknots); if(data1->dim == 3) free_Rvector(data1->z,1,data1->noknots); free_Imatrix(data2->topology,1,data2->noelements,0,data2->maxnodes-1); free_Ivector(data2->material,1,data2->noelements); free_Rvector(data2->x,1,data2->noknots); free_Rvector(data2->y,1,data2->noknots); if(data1->dim == 3) free_Rvector(data2->z,1,data2->noknots); data1->noelements = noelements; data1->noknots = noknots; data1->topology = newtopo; data1->material = newmaterial; data1->elementtypes = newelementtypes; data1->x = newx; data1->y = newy; if(data1->dim == 3) data1->z = newz; if(info) printf("Two meshes were united to one with %d nodes and %d elements.\n", noknots,noelements); return(0); } int CloneMeshes(struct FemType *data,struct BoundaryType *bound, int *ncopies,Real *meshsize,int diffmats,int info) /* Unites two meshes for one larger mesh */ { int i,j,k,l,m; int noelements,noknots,nonodes,totcopies,ind; int **newtopo,*newmaterial,*newelementtypes,maxnodes; int maxmaterial,maxtype,ncopy,bndr,nosides; Real *newx,*newy,*newz; Real maxcoord[3],mincoord[3]; int *vparent,*vparent2,*vside,*vside2,*vtypes,*vmaterial,*vnormal,*vdiscont; Real *vareas; printf("CloneMeshes: copying the mesh to a matrix\n"); if(diffmats) diffmats = 1; totcopies = 1; for(i=0;idim;i++) { if(ncopies[i] > 1) totcopies *= ncopies[i]; } if(data->dim == 2) ncopies[2] = 1; maxcoord[0] = mincoord[0] = data->x[1]; maxcoord[1] = mincoord[1] = data->y[1]; if(data->dim > 2) maxcoord[2] = mincoord[2] = data->z[1]; for(i=1;i<=data->noknots;i++) { if(data->x[i] > maxcoord[0]) maxcoord[0] = data->x[i]; if(data->x[i] < mincoord[0]) mincoord[0] = data->x[i]; if(data->y[i] > maxcoord[1]) maxcoord[1] = data->y[i]; if(data->y[i] < mincoord[1]) mincoord[1] = data->y[i]; if(data->dim > 2) { if(data->z[i] > maxcoord[2]) maxcoord[2] = data->z[i]; if(data->z[i] < mincoord[2]) mincoord[2] = data->z[i]; } } for(i=0;idim;i++) { if(maxcoord[i]-mincoord[i] > meshsize[i]) meshsize[i] = maxcoord[i]-mincoord[i]; } noknots = totcopies * data->noknots; noelements = totcopies * data->noelements; maxnodes = data->maxnodes; printf("Copying the mesh to %d identical domains.\n",totcopies); data->maxnodes = maxnodes; newtopo = Imatrix(1,noelements,0,maxnodes-1); newmaterial = Ivector(1,noelements); newelementtypes = Ivector(1,noelements); newx = Rvector(1,noknots); newy = Rvector(1,noknots); if(data->dim == 3) newz = Rvector(1,noknots); for(l=0;lnoknots;i++) { ncopy = j+k*ncopies[0]+k*l*ncopies[1]; ind = i + ncopy*data->noknots; newx[ind] = data->x[i] + j*meshsize[0]; newy[ind] = data->y[i] + k*meshsize[1]; if(data->dim == 3) newz[ind] = data->z[i] + l*meshsize[2]; } } } } maxmaterial = 0; for(i=1;i<=data->noelements;i++) if(data->material[i] > maxmaterial) maxmaterial = data->material[i]; for(l=0;lnoelements;i++) { ncopy = j+k*ncopies[0]+k*l*ncopies[1]; ind = i + ncopy*data->noelements; newmaterial[ind] = data->material[i] + diffmats*maxmaterial*ncopy; newelementtypes[ind] = data->elementtypes[i]; nonodes = newelementtypes[i]%100; for(m=0;mtopology[i][m] + ncopy*data->noknots; } } } } maxtype = 0; for(j=0;j < MAXBOUNDARIES;j++) { if(!bound[j].created) continue; for(i=1; i <= bound[j].nosides; i++) if(maxtype < bound[j].types[i]) maxtype = bound[j].types[i]; } for(bndr=0;bndr < MAXBOUNDARIES;bndr++) { if(!bound[bndr].created) continue; nosides = totcopies * bound[bndr].nosides; vparent = Ivector(1, nosides); vparent2 = Ivector(1, nosides); vside = Ivector(1, nosides); vside2 = Ivector(1, nosides); vmaterial = Ivector(1, nosides); vtypes = Ivector(1, nosides); vareas = Rvector(1, nosides); vnormal = Ivector(1, nosides); if(bound[bndr].ediscont) { vdiscont = Ivector(1, nosides); for(i=1; i <= nosides; i++) vdiscont[i] = 0; } for(l=0;lnoelements; vside[ind] = bound[bndr].side[i]; if(bound[bndr].parent2[i]) { vparent2[ind] = bound[bndr].parent2[i] + ncopy * data->noelements; vside2[ind] = bound[bndr].side2[i]; } else { vparent2[ind] = 0.0; vside2[ind] = 0.0; } vnormal[ind] = bound[bndr].normal[i]; if(bound[bndr].ediscont) vdiscont[ind] = bound[bndr].discont[i]; vtypes[ind] = bound[bndr].types[i] + diffmats * ncopy * maxtype; vareas[ind] = bound[bndr].areas[i]; vmaterial[ind] = bound[bndr].material[i] + ncopy * maxmaterial; } } } } bound[bndr].nosides = nosides; bound[bndr].side = vside; bound[bndr].side2 = vside2; bound[bndr].parent = vparent; bound[bndr].parent2 = vparent2; bound[bndr].types = vtypes; bound[bndr].areas = vareas; bound[bndr].material = vmaterial; if(bound[bndr].ediscont) bound[bndr].discont = vdiscont; bound[bndr].vfcreated = FALSE; bound[bndr].gfcreated = FALSE; } free_Imatrix(data->topology,1,data->noelements,0,data->maxnodes-1); free_Ivector(data->material,1,data->noelements); free_Rvector(data->x,1,data->noknots); free_Rvector(data->y,1,data->noknots); if(data->dim == 3) free_Rvector(data->z,1,data->noknots); data->noelements = noelements; data->noknots = noknots; data->topology = newtopo; data->material = newmaterial; data->elementtypes = newelementtypes; data->x = newx; data->y = newy; if(data->dim == 3) data->z = newz; if(info) printf("The mesh was copied to several identical meshes\n"); return(0); } int MirrorMeshes(struct FemType *data,struct BoundaryType *bound, int *symmaxis,int diffmats,Real *meshsize,int symmbound,int info) /* Makes a mirror image of a mesh and unites it with the original mesh */ { int i,j,k,l,m,left; int noelements,noknots,nonodes,totcopies,ind; int **newtopo,*newmaterial,*newelementtypes,maxnodes; int maxmaterial,maxtype,ncopy,bndr,nosides; Real *newx,*newy,*newz; Real maxcoord[3],mincoord[3]; int ind0,elem0,symm,axis1,axis2,axis3,symmcount; int *vparent,*vparent2,*vside,*vside2,*vtypes,*vmaterial,*vnormal,*vdiscont; Real *vareas; printf("MirrorMeshes: making a symmetric mapping of the mesh\n"); if(symmaxis[0]) symmaxis[0] = 1; if(symmaxis[1]) symmaxis[1] = 1; if(symmaxis[2]) symmaxis[2] = 1; if(data->dim < 3) symmaxis[2] = 0; maxcoord[0] = mincoord[0] = data->x[1]; maxcoord[1] = mincoord[1] = data->y[1]; if(data->dim > 2) maxcoord[2] = mincoord[2] = data->z[1]; for(i=1;i<=data->noknots;i++) { if(data->x[i] > maxcoord[0]) maxcoord[0] = data->x[i]; if(data->x[i] < mincoord[0]) mincoord[0] = data->x[i]; if(data->y[i] > maxcoord[1]) maxcoord[1] = data->y[i]; if(data->y[i] < mincoord[1]) mincoord[1] = data->y[i]; if(data->dim > 2) { if(data->z[i] > maxcoord[2]) maxcoord[2] = data->z[i]; if(data->z[i] < mincoord[2]) mincoord[2] = data->z[i]; } } for(i=0;idim;i++) { if(maxcoord[i]-mincoord[i] > meshsize[i]) meshsize[i] = maxcoord[i]-mincoord[i]; } if(diffmats) diffmats = 1; totcopies = 1; for(i=0;idim;i++) if(symmaxis[i]) totcopies *= 2; noknots = totcopies * data->noknots; noelements = totcopies * data->noelements; maxnodes = data->maxnodes; printf("Mirroring the mesh to %d symmetrical domains.\n",totcopies); data->maxnodes = maxnodes; newtopo = Imatrix(1,noelements,0,maxnodes-1); newmaterial = Ivector(1,noelements); newelementtypes = Ivector(1,noelements); newx = Rvector(1,noknots); newy = Rvector(1,noknots); if(data->dim == 3) newz = Rvector(1,noknots); ind0 = 0; for(axis1=0;axis1 <= symmaxis[0];axis1++) { for(axis2=0;axis2 <= symmaxis[1];axis2++) { for(axis3=0;axis3 <= symmaxis[2];axis3++) { for(i=1;i<=data->noknots;i++) { ind = i + ind0; newx[ind] = (1-2*axis1) * data->x[i]; newy[ind] = (1-2*axis2) * data->y[i]; if(data->dim == 3) newz[ind] = (1-2*axis3)*data->z[i]; newmaterial[ind] = data->material[i]; newelementtypes[ind] = data->elementtypes[i]; } ind0 += data->noknots; } } } elem0 = 0; ind0 = 0; for(axis1=0;axis1 <= symmaxis[0];axis1++) { for(axis2=0;axis2 <= symmaxis[1];axis2++) { for(axis3=0;axis3 <= symmaxis[2];axis3++) { for(i=1;i<=data->noelements;i++) { ind = i + elem0; newmaterial[ind] = data->material[i]; newelementtypes[ind] = data->elementtypes[i]; nonodes = newelementtypes[i]%100; for(m=0;mtopology[i][m] + ind0; } elem0 += data->noelements; ind0 += data->noknots; printf("elem0=%d ind0=%d\n",elem0,ind0); } } } maxtype = 0; for(j=0;j < MAXBOUNDARIES;j++) { if(!bound[j].created) continue; for(i=1; i <= bound[j].nosides; i++) if(maxtype < bound[j].types[i]) maxtype = bound[j].types[i]; } for(bndr=0;bndr < MAXBOUNDARIES;bndr++) { if(!bound[bndr].created) continue; nosides = totcopies * bound[bndr].nosides; ind = 0; vparent = Ivector(1, nosides); vparent2 = Ivector(1, nosides); vside = Ivector(1, nosides); vside2 = Ivector(1, nosides); vmaterial = Ivector(1, nosides); vtypes = Ivector(1, nosides); vareas = Rvector(1, nosides); vnormal = Ivector(1, nosides); if(bound[bndr].ediscont) { vdiscont = Ivector(1, nosides); for(i=1;i<=nosides;i++) vdiscont[i] = 0; } symmcount = 0; elem0 = 0; ind0 = 0; for(axis1=0;axis1 <= symmaxis[0];axis1++) { for(axis2=0;axis2 <= symmaxis[1];axis2++) { for(axis3=0;axis3 <= symmaxis[2];axis3++) { for(i=1; i <= bound[bndr].nosides; i++) { if(bound[bndr].types[i] == symmbound) continue; ind++; vparent[ind] = bound[bndr].parent[i] + elem0; vparent2[ind] = bound[bndr].parent2[i] + elem0; vside[ind] = bound[bndr].side[i]; vside2[ind] = bound[bndr].side2[i]; vnormal[ind] = bound[bndr].normal[i]; if(bound[bndr].ediscont) vdiscont[ind] = bound[bndr].discont[i]; vtypes[ind] = bound[bndr].types[i] + diffmats * symmcount * maxtype; vareas[ind] = bound[bndr].areas[i]; vmaterial[ind] = bound[bndr].material[i]; } symmcount++; elem0 += data->noelements; } } } nosides = ind; bound[bndr].nosides = nosides; bound[bndr].side = vside; bound[bndr].side2 = vside2; bound[bndr].parent = vparent; bound[bndr].parent2 = vparent2; bound[bndr].types = vtypes; bound[bndr].areas = vareas; bound[bndr].material = vmaterial; if(bound[bndr].ediscont) bound[bndr].discont = vdiscont; bound[bndr].vfcreated = FALSE; bound[bndr].gfcreated = FALSE; } free_Imatrix(data->topology,1,data->noelements,0,data->maxnodes-1); free_Ivector(data->material,1,data->noelements); free_Rvector(data->x,1,data->noknots); free_Rvector(data->y,1,data->noknots); if(data->dim == 3) free_Rvector(data->z,1,data->noknots); data->noelements = noelements; data->noknots = noknots; data->topology = newtopo; data->material = newmaterial; data->elementtypes = newelementtypes; data->x = newx; data->y = newy; if(data->dim == 3) data->z = newz; if(info) printf("The mesh was copied to several identical meshes\n"); return(0); } static void ReorderAutomatic(struct FemType *data,int iterations, int *origindx,Real corder[],int info) { int i,j,k,l,nonodes,maxnodes,noelements,noknots,minwidth,indexwidth; int **neighbours,*newrank,*newindx,*oldrank,*oldindx; int nocands,*cands,ind,ind2,cantdo; int elemtype,indready,iter,*localorder,*localtmp,nolocal; Real *localdist,x,y,z,dx,dy,dz,sum; iterations = 3; iter = 0; maxnodes = 8; noelements = data->noelements; noknots = data->noknots; cantdo = FALSE; for(j=1;j<=noelements;j++) { elemtype = data->elementtypes[j]; if(elemtype != 404 && elemtype != 303 && elemtype != 808) cantdo = elemtype; } if(cantdo) { printf("Automatic reordering not specified for elementtype %d\n",cantdo); return; } printf("Allocating...\n"); cands = Ivector(1,maxnodes); localorder = Ivector(1,maxnodes); localtmp = Ivector(1,maxnodes); localdist = Rvector(1,maxnodes); neighbours = Imatrix(1,noknots,1,maxnodes); newrank = Ivector(1,noknots); oldrank = Ivector(1,noknots); newindx = Ivector(1,noknots); oldindx = Ivector(1,noknots); for(i=1;i<=noknots;i++) oldindx[i] = origindx[i]; for(i=1;i<=noknots;i++) oldrank[origindx[i]] = i; minwidth = CalculateIndexwidth(data,TRUE,oldrank); if(info) printf("Indexwidth of the initial node order is %d.\n",minwidth); for(j=1;j<=noknots;j++) for(i=1;i<=maxnodes;i++) neighbours[j][i] = 0; if(info) printf("Initializing neighbours\n"); for(j=1;j<=noelements;j++) { elemtype = data->elementtypes[j]; nonodes = elemtype%100; for(i=0;itopology[j][i]; if(elemtype == 404 || elemtype == 303) { nocands = 2; cands[1] = (i+1)%nonodes; cands[2] = (i+nonodes-1)%nonodes; } else if(elemtype == 808) { nocands = 3; if(i<4) { cands[1] = (i+1)%4; cands[2] = (i+3)%4; cands[3] = i+4; } else { cands[1] = (i-4+1)%4+4; cands[2] = (i-4+3)%4+4; cands[3] = i-4; } } for(k=1;k<=nocands;k++) { ind2 = data->topology[j][cands[k]]; for(l=1;l<=maxnodes;l++) { if(neighbours[ind][l] == 0) break; if(neighbours[ind][l] == ind2) ind2 = 0; } if(ind2) neighbours[ind][l] = ind2; } } } #if 0 for(j=1;j<=noknots;j++) { printf("neigbourds[%d]= ",j); for(l=1;l<=maxnodes;l++) printf("%d ",neighbours[j][l]); printf("\n"); } #endif if(info) printf("Reordering neighbours table\n"); for(j=1;j<=noknots;j++) { nolocal = 0; dz = 0.0; for(l=1;l<=maxnodes;l++){ if(ind = neighbours[j][l]) { nolocal++; localtmp[l] = ind; dx = data->x[l] - data->x[ind]; dy = data->y[l] - data->y[ind]; if(data->dim == 3) dz = data->z[l] - data->z[ind]; localdist[l] = corder[0]*fabs(dx) + corder[1]*fabs(dy) + corder[2]*fabs(dz); } } SortIndex(nolocal,localdist,localorder); for(l=1;l<=nolocal;l++) neighbours[j][l] = localtmp[localorder[l]]; #if 0 for(l=1;l<=nolocal;l++) printf("j=%d l=%d dist=%.3le order=%d %d\n", j,l,localdist[l],localorder[l],neighbours[j][l]); #endif } for(iter=1;iter<=iterations;iter++) { if(info) printf("Optimal topology testing %d\n",iter); for(i=1;i<=noknots;i++) newrank[i] = 0; ind = 0; indready = 1; do { if(indready > ind) { for(l=noknots;l>=1;l--) if(j = oldindx[l]) break; if(info) printf("Starting over from node %d when ind=%d indready=%d\n",j,ind,indready); } else { j = newindx[indready] ; } for(l=1;ind2 = neighbours[j][l];l++) { if(ind2) { if(!newrank[ind2]) { ind++; newrank[ind2] = ind; newindx[ind] = ind2; oldindx[oldrank[ind2]] = 0; oldrank[ind2] = 0; } } } indready++; } while(ind < noknots); indexwidth = CalculateIndexwidth(data,TRUE,newrank); if(info) printf("Indexwidth of the suggested node order is %d.\n",indexwidth); for(i=1;i<=noknots;i++) oldrank[i] = newrank[i]; for(i=1;i<=noknots;i++) oldindx[i] = newindx[i]; if(indexwidth < minwidth) { for(i=1;i<=noknots;i++) origindx[i] = newindx[i]; minwidth = indexwidth; } } thened: free_Ivector(cands,1,maxnodes); free_Ivector(localorder,1,maxnodes); free_Ivector(localtmp,1,maxnodes); free_Rvector(localdist,1,maxnodes); free_Imatrix(neighbours,1,noknots,1,maxnodes); free_Ivector(newrank,1,noknots); free_Ivector(oldrank,1,noknots); free_Ivector(newindx,1,noknots); free_Ivector(oldindx,1,noknots); } void ReorderElements(struct FemType *data,struct BoundaryType *bound, int manual,Real corder[],int info) { int i,j,k,l,side,automatic; int noelements,noknots,nonodes,length; int **newtopology,*newmaterial,*newelementtypes; int *indx,*revindx,*elemindx,*revelemindx; int oldnoknots, oldnoelements; Real *newx,*newy,*newz,*arrange; Real dx,dy,dz,cx,cy,cz,cbase; noelements = oldnoelements = data->noelements; noknots = oldnoknots = data->noknots; if(info) printf("Reordering %d knots and %d elements in %d-dimensions.\n", noknots,noelements,data->dim); if(noelements > noknots) length = noelements; else length = noknots; arrange = Rvector(1,length); indx = Ivector(1,noknots); revindx = Ivector(1,noknots); elemindx = Ivector(1,noelements); revelemindx = Ivector(1,noelements); if(manual == 1) { cx = corder[0]; cy = corder[1]; cz = corder[2]; } else { Real xmin,xmax,ymin,ymax,zmin,zmax; xmin = xmax = data->x[1]; ymin = ymax = data->y[1]; if(data->dim == 3) zmin = zmax = data->z[1]; for(i=1;i<=data->noknots;i++) { if(xmin > data->x[i]) xmin = data->x[i]; if(xmax < data->x[i]) xmax = data->x[i]; if(ymin > data->y[i]) ymin = data->y[i]; if(ymax < data->y[i]) ymax = data->y[i]; if(data->dim == 3) { if(zmin > data->z[i]) zmin = data->z[i]; if(zmax < data->z[i]) zmax = data->z[i]; } } dx = xmax-xmin; dy = ymax-ymin; if(data->dim == 3) dz = zmax-zmin; else dz = 0.0; /* The second strategy seems to be better in many cases */ #if 0 cx = dx; cy = dy; cz = dz; #else cbase = 100.0; cx = pow(cbase,1.0*(dx>dy)+1.0*(dx>dz)); cy = pow(cbase,1.0*(dy>dx)+1.0*(dy>dz)); cz = pow(cbase,1.0*(dz>dx)+1.0*(dz>dx)); #endif corder[0] = cx; corder[1] = cy; corder[2] = cz; } if(info) printf("Ordering with (%.3lg*x + %.3lg*y + %.3lg*z)\n",cx,cy,cz); for(i=1;i<=noknots;i++) { arrange[i] = cx*data->x[i] + cy*data->y[i]; if(data->dim == 3) arrange[i] += cz*data->z[i]; } SortIndex(noknots,arrange,indx); if(manual == 2) ReorderAutomatic(data,0,indx,corder,TRUE); for(i=1;i<=noknots;i++) revindx[indx[i]] = i; for(j=1;j<=noelements;j++) { nonodes = data->elementtypes[j]%100; arrange[j] = 0.0; for(i=0;itopology[j][i]; arrange[j] += cx*data->x[k] + cy*data->y[k]; if(data->dim == 3) arrange[j] += cz*data->z[k]; } } SortIndex(noelements,arrange,elemindx); for(i=1;i<=noelements;i++) revelemindx[elemindx[i]] = i; #if 0 for(i=1;i<=noknots;i++) printf("i=%d indx=%d revi=%d f=%.2lg\n", i,indx[i],revindx[i],arrange[indx[i]]); #endif if(info) printf("Moving knots to new positions\n"); newx = Rvector(1,data->noknots); newy = Rvector(1,data->noknots); if(data->dim == 3) newz = Rvector(1,data->noknots); for(i=1;i<=data->noknots;i++) { newx[i] = data->x[indx[i]]; newy[i] = data->y[indx[i]]; if(data->dim == 3) newz[i] = data->z[indx[i]]; } free_Rvector(data->x,1,data->noknots); free_Rvector(data->y,1,data->noknots); if(data->dim == 3) free_Rvector(data->z,1,data->noknots); data->x = newx; data->y = newy; if(data->dim == 3) data->z = newz; if(info) printf("Moving the elements to new positions\n"); newtopology = Imatrix(1,noelements,0,data->maxnodes-1); newmaterial = Ivector(1,noelements); newelementtypes = Ivector(1,noelements); for(j=1;j<=noelements;j++) { newmaterial[j] = data->material[elemindx[j]]; newelementtypes[j] = data->elementtypes[elemindx[j]]; nonodes = newelementtypes[j]%100; for(i=0;itopology[elemindx[j]][i]; newtopology[j][i] = revindx[k]; } } data->material = newmaterial; data->elementtypes = newelementtypes; data->topology = newtopology; printf("Moving the parents of the boundary nodes.\n"); for(j=0;j < MAXBOUNDARIES;j++) { if(!bound[j].created) continue; for(i=1; i <= bound[j].nosides; i++) { bound[j].parent[i] = revelemindx[bound[j].parent[i]]; if(bound[j].parent2[i]) bound[j].parent2[i] = revelemindx[bound[j].parent2[i]]; } } i = CalculateIndexwidth(data,FALSE,indx); printf("Indexwidth of the new node order is %d.\n",i); free_Rvector(arrange,1,length); free_Ivector(indx,1,oldnoknots); free_Ivector(revindx,1,oldnoknots); free_Ivector(elemindx,1,oldnoelements); free_Ivector(revelemindx,1,oldnoelements); } int RemoveUnusedNodes(struct FemType *data,int info) { int i,j,k; int noelements,noknots,nonodes,activeknots; int *indx; noelements = data->noelements; noknots = data->noknots; indx = Ivector(1,noknots); for(i=1;i<=noknots;i++) indx[i] = 0; for(j=1;j<=noelements;j++) { nonodes = data->elementtypes[j] % 100; for(i=0;itopology[j][i] ] = 1; } activeknots = 0; for(i=1;i<=noknots;i++) { if(indx[i]) { activeknots += 1; indx[i] = activeknots; } } if( noknots == activeknots) { if(info) printf("All %d nodes were used by the mesh elements\n",noknots); return(1); } if(info) printf("Removing %d unused nodes (out of %d) from the mesh\n",noknots-activeknots,noknots); for(j=1;j<=noelements;j++) { nonodes = data->elementtypes[j] % 100; for(i=0;itopology[j][i] = indx[ data->topology[j][i] ]; } for(i=1;i<=noknots;i++) { j = indx[i]; if(!j) continue; data->x[j] = data->x[i]; data->y[j] = data->y[i]; if(data->dim == 3) data->z[j] = data->z[i]; } data->noknots = activeknots; free_Ivector(indx,1,noknots); return(0); } void RenumberBoundaryTypes(struct FemType *data,struct BoundaryType *bound, int renumber, int bcoffset, int info) { int i,j,k,l,doinit; int minbc,maxbc,*mapbc; if(renumber) { if(info) printf("Renumbering boundary types\n"); doinit = TRUE; for(j=0;j < MAXBOUNDARIES;j++) { if(!bound[j].created) continue; for(i=1;i<=bound[j].nosides;i++) { if(doinit) { maxbc = minbc = bound[j].types[i]; doinit = FALSE; } maxbc = MAX(maxbc,bound[j].types[i]); minbc = MIN(minbc,bound[j].types[i]); } } if(doinit) return; mapbc = Ivector(minbc,maxbc); for(i=minbc;i<=maxbc;i++) mapbc[i] = 0; for(j=0;j < MAXBOUNDARIES;j++) { if(!bound[j].created) continue; for(i=1;i<=bound[j].nosides;i++) mapbc[bound[j].types[i]] = TRUE; } j = 0; for(i=minbc;i<=maxbc;i++) if(mapbc[i]) { j++; mapbc[i] = j; } if(maxbc - minbc >= j || minbc != 1) { if(info) printf("Mapping boundary types from [%d %d] to [%d %d]\n",minbc,maxbc,1,j); for(j=0;j < MAXBOUNDARIES;j++) { if(!bound[j].created) continue; for(i=1;i<=bound[j].nosides;i++) bound[j].types[i] = mapbc[bound[j].types[i]]; } if(data->boundarynamesexist) { for(j=minbc;j<=MIN(maxbc,MAXBODIES-1);j++) { if(mapbc[j]) strcpy(data->boundaryname[mapbc[j]],data->boundaryname[j]); } } } free_Ivector(mapbc,minbc,maxbc); } if(bcoffset) { if(info) printf("Adding offset of %d to the BCs\n",bcoffset); for(j=0;j < MAXBOUNDARIES;j++) { if(!bound[j].created) continue; for(i=1;i<=bound[j].nosides;i++) bound[j].types[i] += bcoffset; } if(data->boundarynamesexist) { for(j=MAXBOUNDARIES-bcoffset-1;j>=0;j--) { strcpy(data->boundaryname[j+bcoffset],data->boundaryname[j]); } } } else { if(info) printf("BCs ordered continously between %d and %d\n",minbc,maxbc); } } void RenumberMaterialTypes(struct FemType *data,struct BoundaryType *bound,int info) { int i,j,k,l,noelements,doinit; int minmat,maxmat,*mapmat; if(info) printf("Setting new material types\n"); noelements = data->noelements; if(noelements < 1) { printf("There are no elements to set!\n"); return; } doinit = TRUE; for(j=1;j<=noelements;j++) { if(doinit) { maxmat = minmat = data->material[j]; doinit = FALSE; } maxmat = MAX(maxmat,data->material[j]); minmat = MIN(minmat,data->material[j]); } mapmat = Ivector(minmat,maxmat); for(i=minmat;i<=maxmat;i++) mapmat[i] = 0; for(j=1;j<=noelements;j++) mapmat[data->material[j]] = TRUE; j = 0; for(i=minmat;i<=maxmat;i++) if(mapmat[i]) mapmat[i] = ++j; if(maxmat - minmat >= j || minmat != 1) { if(info) printf("Mapping material types from [%d %d] to [%d %d]\n", minmat,maxmat,1,j); for(j=1;j<=noelements;j++) data->material[j] = mapmat[data->material[j]]; if(data->bodynamesexist) { for(j=minmat;j<=MIN(maxmat,MAXBODIES-1);j++) { if(mapmat[j]) strcpy(data->bodyname[mapmat[j]],data->bodyname[j]); } } } else { if(info) printf("Materials ordered continously between %d and %d\n",minmat,maxmat); } free_Ivector(mapmat,minmat,maxmat); } int RemoveLowerDimensionalBoundaries(struct FemType *data,struct BoundaryType *bound,int info) { int i,j,k,l,noelements; int elemtype,maxelemtype,maxelemdim,elemdim; int parent, side, sideind[MAXNODESD1],sideelemtype; int nosides, oldnosides,newnosides; if(info) printf("Removing lower dimensional boundaries\n"); noelements = data->noelements; if(noelements < 1) return(1); maxelemtype = GetMaxElementType(data); maxelemdim = GetElementDimension(maxelemtype); if(info) printf("Maximum elementtype is %d and dimension %d\n",maxelemtype,maxelemdim); if(maxelemdim < 2) return(2); oldnosides = 0; newnosides = 0; for(j=0;j < MAXBOUNDARIES;j++) { nosides = 0; if(!bound[j].created) continue; for(i=1;i<=bound[j].nosides;i++) { oldnosides++; parent = bound[j].parent[i]; side = bound[j].side[i]; GetElementSide(parent,side,1,data,sideind,&sideelemtype); if(sideelemtype > 300) elemdim = 2; else if(sideelemtype > 200) elemdim = 1; else elemdim = 0; if(maxelemdim - elemdim > 1) continue; nosides++; if(nosides == i) continue; bound[j].parent[nosides] = bound[j].parent[i]; bound[j].parent2[nosides] = bound[j].parent2[i]; bound[j].side[nosides] = bound[j].side[i]; bound[j].side2[nosides] = bound[j].side2[i]; bound[j].types[nosides] = bound[j].types[i]; } bound[j].nosides = nosides; newnosides += nosides; } if(info) printf("Removed %d (out of %d) less than %dD boundary elements\n", oldnosides-newnosides,oldnosides,maxelemdim); return(0); } static void FindEdges(struct FemType *data,struct BoundaryType *bound, int material,int sidetype,int info) { int i,j,k,l,side,identical,element; int noelements,noknots,nonodes,nomaterials,nosides,newbound; int maxelementtype,maxedgenodes,elemedges,maxelemedges,edge,dosides; int **edgetable,sideind[MAXNODESD1],sideelemtype,allocated; int *indx,*revindx; Real *arrange; nomaterials = 0; maxelementtype = 0; printf("FindEdges: Finding edges of bulk elements of type %d\n",material); maxelementtype = GetMaxElementType(data); if(maxelementtype/100 > 4) { printf("FindEdges: Implemented only for 2D elements!\n"); dosides = 0; return; } if(maxelementtype/100 <= 2) maxedgenodes = 1; else if(maxelementtype/100 <= 4) maxedgenodes = 2; maxelemedges = maxelementtype/100; edgetable = Imatrix(1,maxelemedges*nomaterials,0,maxedgenodes+1); for(i=1;i<=maxelemedges*nomaterials;i++) for(j=0;j<=maxedgenodes+1;j++) edgetable[i][j] = 0; edge = 0; for(element=1;element<=data->noelements;element++) { if(data->material[element] != material) continue; elemedges = data->elementtypes[element]/100; for(side=0;side sideind[1]) { edgetable[edge][0] = sideind[0]; edgetable[edge][1] = sideind[1]; } else { edgetable[edge][1] = sideind[0]; edgetable[edge][0] = sideind[1]; } } } } noknots = edge; arrange = Rvector(1,noknots); for(i=1;i<=noknots;i++) arrange[i] = 0.0; for(i=1;i<=noknots;i++) arrange[i] = edgetable[i][0]; indx = Ivector(1,noknots); SortIndex(noknots,arrange,indx); #if 0 printf("noknots = %d\n",noknots); for(i=1;i<=noknots;i++) printf("indx[%d]=%d edge=%d arrange[%d] = %lg arrange[indx[%d]] = %lg\n", i,indx[i],edgetable[i][0],i,arrange[i],i,arrange[indx[i]]); #endif #if 0 revindx = Ivector(1,data->noknots); for(i=1;i<=noknots;i++) revindx[indx[i]] = i; #endif allocated = FALSE; omstart: nosides = 0; for(i=1;i<=noknots;i++) { identical = FALSE; if(maxedgenodes == 1) { for(j=i+1;j<=noknots && edgetable[indx[i]][0] == edgetable[indx[j]][0];j++) identical = TRUE; for(j=i-1;j>=1 && edgetable[indx[i]][0] == edgetable[indx[j]][0];j--) identical = TRUE; } else if(maxedgenodes == 2) { for(j=i+1;j<=noknots && edgetable[indx[i]][0] == edgetable[indx[j]][0];j++) if(edgetable[indx[i]][1] == edgetable[indx[j]][1]) identical = TRUE; for(j=i-1;j>=1 && edgetable[indx[i]][0] == edgetable[indx[j]][0];j--) if(edgetable[indx[i]][1] == edgetable[indx[j]][1]) identical = TRUE; } if(identical) continue; nosides++; if(allocated) { bound[newbound].parent[nosides] = edgetable[indx[i]][maxedgenodes]; bound[newbound].parent2[nosides] = 0; bound[newbound].side[nosides] = edgetable[indx[i]][maxedgenodes+1]; bound[newbound].side2[nosides] = 0; bound[newbound].types[nosides] = sidetype; } } if(!allocated) { for(j=0;j < MAXBOUNDARIES && bound[j].created;j++); newbound = j; AllocateBoundary(&bound[newbound],nosides); allocated = TRUE; if(info) printf("Created boundary %d of type %d and size %d for material %d\n", newbound,sidetype,nosides,material); goto omstart; } free_Ivector(indx,1,noknots); free_Imatrix(edgetable,1,maxelemedges*nomaterials,0,maxedgenodes+1); } static int CompareIndexes(int elemtype,int *ind1,int *ind2) { int i,j,same,nosides,hits; hits = 0; nosides = elemtype / 100; for(i=0;idim; if(dimred) lowerdim = dim - dimred; else lowerdim = dim-1; noboundnodes = 0; for(i=1;i<=data->noknots;i++) if(boundnodes[i]) noboundnodes++; if(!noboundnodes) { printf("FindNewBoundaries: no nonzero entries in boundnodes vector!\n"); return(1); } else { if(info) printf("There are %d nonzero entries in boundnodes vector!\n",noboundnodes); } omstart: nosides = 0; for(element=1;element<=data->noelements;element++) { elemtype = data->elementtypes[element]; if(dim == 1) { minedge = 0; maxedge = elemtype/100 -1; } else if(dim == 2) { if(lowerdim == 1) { minedge = 0; maxedge = elemtype/100 -1; } else if(lowerdim == 0) { minedge = elemtype/100; maxedge = minedge + 1; } } else if(dim == 3) { if(lowerdim == 2) { minedge = 0; if(elemtype/100 == 5) maxedge = 3; else if(elemtype/100 == 6 || elemtype/100 == 7) maxedge = 4; else if(elemtype/100 == 8) maxedge = 5; } else if(lowerdim == 1) { if(elemtype/100 == 8) { minedge = 6; maxedge = 17; } else if(elemtype/100 == 5) { minedge = 4; maxedge = 9; } else printf("FindNewBoundaries: not implemented for all 3d boundaries\n"); } else if(lowerdim == 0) { if(elemtype/100 == 8) { minedge = 18; maxedge = 25; } } } for(side=minedge;side<=maxedge;side++) { GetElementSide(element,side,1,data,sideind,&sideelemtype); nonodes = sideelemtype % 100; identical = TRUE; for(i=0;i newtype) newtype = bound[j].types[i]; } newbound++; if(!suggesttype) newtype++; AllocateBoundary(&bound[newbound],nosides); allocated = TRUE; if(info) printf("Allocating for %d sides of boundary %d\n",nosides,newtype); goto omstart; } bound[newbound].nosides = nosides; if(info) printf("Found %d sides of dim %d to define boundary %d\n",nosides,lowerdim,newtype); for(i=1;i<=nosides;i++) { if(j = bound[newbound].parent2[i]) { if(bound[newbound].parent[i] > bound[newbound].parent2[i]) { bound[newbound].parent2[i] = bound[newbound].parent[i]; bound[newbound].parent[i] = j; j = bound[newbound].side2[i]; bound[newbound].side2[i] = bound[newbound].side[i]; bound[newbound].side[i] = j; } } } } else { if(lowerdim == 0) { printf("The nodes do not form a boundary!\n"); return(2); } else { lowerdim--; printf("The nodes do not form a boundary, trying with %d-dimensional elements.\n",lowerdim); goto omstart; } } return(0); } int FindBulkBoundary(struct FemType *data,int mat1,int mat2, int *boundnodes,int *noboundnodes,int info) { int i,j,k; int hits,nonodes,maxnodes,minnodes,material,elemtype; Real ds,xmin,xmax,ymin,ymax,zmin,zmax,eps; int *visited,elemind[MAXNODESD2],elemdim,*ind; Real *anglesum,dx1,dx2,dy1,dy2,dz1,dz2,ds1,ds2,dotprod; eps = 1.0e-4; *noboundnodes = 0; if(mat1 < 1 && mat2 < 1) { printf("FindBulkBoundaty: Either of the materials must be positive\n"); return(1); } else if(mat1 < 1) { i = mat1; mat1 = mat2; mat2 = i; } if(info) printf("Finding nodes between bulk elements of material %d and %d\n",mat1,mat2); visited = Ivector(1,data->noknots); for(i=1;i<=data->noknots;i++) visited[i] = 0; for(i=1;i<=data->noknots;i++) boundnodes[i] = 0; elemdim = 0; for(i=1;i<=data->noelements;i++) { material = data->material[i]; if(material == mat1) { nonodes = data->elementtypes[i] % 100; k = data->elementtypes[i]/100; if(k > elemdim) elemdim = k; for(j=0;jtopology[i][j]; visited[k] += 1; } } } maxnodes = minnodes = visited[1]; for(i=1;i<=data->noknots;i++) { if(visited[i] > maxnodes) maxnodes = visited[i]; if(visited[i] < minnodes) minnodes = visited[i]; } if(elemdim == 3 || elemdim == 4) { anglesum = Rvector(1, data->noknots); for(i=1;i<=data->noknots;i++) anglesum[i] = 0.0; for(i=1;i<=data->noelements;i++) { material = data->material[i]; if(material == mat1) { nonodes = data->elementtypes[i]/100; ind = data->topology[i]; if(nonodes == 3 || nonodes == 4) { for(k=0;kx[ind[(k+1)%nonodes]] - data->x[ind[k]]; dy1 = data->y[ind[(k+1)%nonodes]] - data->y[ind[k]]; dz1 = data->z[ind[(k+1)%nonodes]] - data->z[ind[k]]; dx2 = data->x[ind[(k+nonodes-1)%nonodes]] - data->x[ind[k]]; dy2 = data->y[ind[(k+nonodes-1)%nonodes]] - data->y[ind[k]]; dz2 = data->z[ind[(k+nonodes-1)%nonodes]] - data->z[ind[k]]; ds1 = sqrt(dx1*dx1+dy1*dy1+dz1*dz1); ds2 = sqrt(dx2*dx2+dy2*dy2+dz2*dz2); dotprod = dx1*dx2 + dy1*dy2 + dz1*dz2; anglesum[ind[k]] += acos(dotprod / (ds1*ds2)); } } } } j = 0; for(i=1;i<=data->noknots;i++) { anglesum[i] /= 2.0 * FM_PI; if(anglesum[i] > 0.99) visited[i] = 0; if(anglesum[i] > 1.01) printf("FindBulkBoundary: surpricingly large angle %.3le in node %d\n",anglesum[i],i); if(visited[i]) j++; } if(0) printf("There are %d boundary node candidates\n",j); free_Rvector(anglesum,1,data->noknots); } else { for(i=1;i<=data->noknots;i++) if(visited[i] == maxnodes) visited[i] = 0; if(maxnodes < 2) { printf("FindBulkBoundary: Nodes must belong to more than %d elements.\n",maxnodes); return(2); } } if(mat2 == 0) { for(i=1;i<=data->noelements;i++) { material = data->material[i]; if(material == mat1) continue; nonodes = data->elementtypes[i] % 100; for(j=0;jtopology[i][j]; boundnodes[k] += 1; } } for(k=1;k<=data->noknots;k++) { if(!visited[k]) boundnodes[k] = 0; else if(visited[k] < boundnodes[k]) boundnodes[k] = 0; else if(visited[k] + boundnodes[k] < maxnodes) boundnodes[k] = 1; else boundnodes[k] = 0; } } else if(mat2 == -10) { for(i=1;i<=data->noknots;i++) if(visited[i]) boundnodes[i] = 1; } else if(mat2 == -11 || mat2 == -12 || mat2 > 0) { for(i=1;i<=data->noelements;i++) { material = data->material[i]; if(material == mat1) continue; if(mat2 > 0 && material != mat2) continue; if(mat2 == -11 && material < mat1) continue; if(mat2 == -12 && material > mat1) continue; nonodes = data->elementtypes[i]%100; for(j=0;jtopology[i][j]; if(visited[k]) boundnodes[k] = 1; } } } else if(mat2 >= -2*data->dim && mat2 <= -1) { j = TRUE; for(i=1;i<=data->noknots;i++) if(visited[i]) { if(j) { xmax = xmin = data->x[k]; if(data->dim >= 2) ymax = ymin = data->y[k]; if(data->dim >= 3) zmax = zmin = data->z[k]; j = FALSE; } else { if(data->x[i] > xmax) xmax = data->x[i]; if(data->x[i] < xmin) xmin = data->x[i]; if(data->dim >= 2) { if(data->y[i] > ymax) ymax = data->y[i]; if(data->y[i] < ymin) ymin = data->y[i]; } if(data->dim >= 3) { if(data->z[i] > zmax) zmax = data->z[i]; if(data->z[i] < zmin) zmin = data->z[i]; } } } ds = (xmax-xmin)*(xmax-xmin); if(data->dim >= 2) ds += (ymax-ymin)*(ymax-ymin); if(data->dim >= 3) ds += (zmax-zmin)*(zmax-zmin); ds = sqrt(ds); eps = 1.0e-5 * ds; for(i=1;i<=data->noknots;i++) if(visited[i] < maxnodes && visited[i]) { if(data->dim == 1) { if(mat2 == -1 && fabs(data->x[i]-xmin) < eps) boundnodes[i] = 1; else if(mat2 == -2 && fabs(data->x[i]-xmax) < eps) boundnodes[i] = 1; } if(data->dim >= 2) { if(mat2 == -1 && (fabs(data->y[i]-ymin) < eps)) boundnodes[i] = 1; else if(mat2 == -3 && (fabs(data->y[i]-ymax) < eps)) boundnodes[i] = 1; else if(mat2 == -4 && (fabs(data->x[i]-xmin) < eps)) boundnodes[i] = 1; else if(mat2 == -2 && (fabs(data->x[i]-xmax) < eps)) boundnodes[i] = 1; } if(data->dim >= 3) { if(mat2 == -5 && fabs(data->z[i]-zmin) < eps) boundnodes[i] = 1; else if(mat2 == -6 && fabs(data->z[i]-zmax) < eps) boundnodes[i] = 1; } } } else { printf("FindBulkBoundary: unknown option %d for finding a side\n",mat2); return(2); } *noboundnodes = 0; for(i=1;i<=data->noknots;i++) if(boundnodes[i]) *noboundnodes += 1; if(info) printf("Located %d nodes at the interval between materials %d and %d\n", *noboundnodes,mat1,mat2); free_Ivector(visited,1,data->noknots); return(0); } int FindBoundaryBoundary(struct FemType *data,struct BoundaryType *bound,int mat1,int mat2, int *boundnodes,int *noboundnodes,int info) { int i,j,k,l; int hits,nonodes,nocorners,maxnodes,minnodes,elemtype,material,bounddim; Real ds,xmin,xmax,ymin,ymax,zmin,zmax,eps,dx1,dx2,dy1,dy2,dz1,dz2,ds1,ds2,dotprod; Real *anglesum; int *visited,sideind[MAXNODESD2],elemind[MAXNODESD2]; eps = 1.0e-4; *noboundnodes = 0; if(mat1 < 1 && mat2 < 1) { printf("FindBoundaryBoundaty: Either of the boundaries must be positive\n"); return(1); } else if(mat1 < 1) { i = mat1; mat1 = mat2; mat2 = i; } if(info) printf("Finding nodes between boundary elements of type %d and %d\n",mat1,mat2); visited = Ivector(1,data->noknots); for(i=1;i<=data->noknots;i++) visited[i] = 0; for(i=1;i<=data->noknots;i++) boundnodes[i] = 0; bounddim = 0; /* Set a tag to all nodes that are part of the other boundary */ for(j=0;j < MAXBOUNDARIES;j++) { if(!bound[j].created) continue; for(i=1; i <= bound[j].nosides; i++) { if(bound[j].types[i] == mat1) { GetElementSide(bound[j].parent[i],bound[j].side[i],bound[j].normal[i], data,sideind,&elemtype); nonodes = elemtype % 100; nocorners = elemtype / 100; for(k=0;knoknots); for(k=1;k<=data->noknots;k++) anglesum[k] = 0.0; bounddim = 2; } nonodes = nocorners; for(k=0;kx[sideind[(k+1)%nonodes]] - data->x[sideind[k]]; dy1 = data->y[sideind[(k+1)%nonodes]] - data->y[sideind[k]]; dz1 = data->z[sideind[(k+1)%nonodes]] - data->z[sideind[k]]; dx2 = data->x[sideind[(k+nonodes-1)%nonodes]] - data->x[sideind[k]]; dy2 = data->y[sideind[(k+nonodes-1)%nonodes]] - data->y[sideind[k]]; dz2 = data->z[sideind[(k+nonodes-1)%nonodes]] - data->z[sideind[k]]; ds1 = sqrt(dx1*dx1+dy1*dy1+dz1*dz1); ds2 = sqrt(dx2*dx2+dy2*dy2+dz2*dz2); dotprod = dx1*dx2 + dy1*dy2 + dz1*dz2; anglesum[sideind[k]] += acos(dotprod / (ds1*ds2)); } } } } } maxnodes = minnodes = abs(visited[1]); for(i=1;i<=data->noknots;i++) { j = abs( visited[i] ); maxnodes = MAX( j, maxnodes ); minnodes = MIN( j, minnodes ); } if(info) printf("There are from %d to %d hits per node\n",minnodes,maxnodes); if(maxnodes < 2) { printf("FindBulkBoundary: Nodes must belong to more than %d elements.\n",maxnodes); return(2); } if(bounddim == 2) { /* For corner nodes eliminate the ones with full angle */ for(i=1;i<=data->noknots;i++) { anglesum[i] /= 2.0 * FM_PI; if(anglesum[i] > 0.99) visited[i] = 0; if(anglesum[i] > 1.01) printf("FindBulkBoundary: surpricingly large angle %.3le in node %d\n",anglesum[i],i); } free_Rvector(anglesum,1,data->noknots); /* For higher order nodes eliminate the ones with more than one hits */ k = 0; for(i=1;i<=data->noknots;i++) { if(visited[i] == -1) visited[i] = 1; else if(visited[i] < -1) { k++; visited[i] = 0; } } if(k && info) printf("Removed %d potential higher order side nodes from the list.\n",k); } if(bounddim == 1) { if(visited[i] == maxnodes || visited[i] < 0) visited[i] = 0; } /* Neighbour to anaything */ if(mat2 == 0) { for(k=1;k<=data->noknots;k++) if(visited[k]) boundnodes[k] = 1; } /* Neighbour to other BCs */ else if(mat2 == -11 || mat2 == -12 || mat2 == -10 || mat2 > 0) { for(j=0;j < MAXBOUNDARIES;j++) { if(!bound[j].created) continue; for(i=1; i <= bound[j].nosides; i++) { material = bound[j].types[i]; if(material == mat1) continue; if(mat2 > 0 && material != mat2) continue; if(mat2 == -11 && material < mat1) continue; if(mat2 == -12 && material > mat1) continue; GetElementSide(bound[j].parent[i],bound[j].side[i],bound[j].normal[i], data,sideind,&elemtype); nonodes = elemtype%100; for(k=0;k= -2*data->dim && mat2 <= -1) { for(j=0;j < MAXBOUNDARIES;j++) { if(!bound[j].created) continue; for(i=1; i <= bound[j].nosides; i++) { material = bound[j].types[i]; if(material != mat1) continue; GetElementSide(bound[j].parent[i],bound[j].side[i],bound[j].normal[i], data,sideind,&elemtype); nonodes = elemtype%100; hits = 0; for(k=0;kx[l]; if(data->dim >= 2) ymax = ymin = data->y[l]; if(data->dim >= 3) zmax = zmin = data->z[l]; for(k=1;kx[l] > xmax) xmax = data->x[l]; if(data->x[l] < xmin) xmin = data->x[l]; if(data->dim >= 2) { if(data->y[l] > ymax) ymax = data->y[l]; if(data->y[l] < ymin) ymin = data->y[l]; } if(data->dim >= 3) { if(data->z[l] > zmax) zmax = data->z[l]; if(data->z[l] < zmin) zmin = data->z[l]; } } ds = (xmax-xmin)*(xmax-xmin); if(data->dim >= 2) ds += (ymax-ymin)*(ymax-ymin); if(data->dim >= 3) ds += (zmax-zmin)*(zmax-zmin); ds = sqrt(ds); eps = 1.0e-3 * ds; for(k=0;kdim == 1) { if(mat2 == -1 && fabs(data->x[l]-xmin) < eps) boundnodes[l] = 1; else if(mat2 == -2 && fabs(data->x[l]-xmax) < eps) boundnodes[l] = 1; } if(data->dim >= 2) { if(mat2 == -1 && (fabs(data->y[l]-ymin) < eps)) elemind[l] = 1; else if(mat2 == -3 && (fabs(data->y[l]-ymax) < eps)) elemind[l] = 1; else if(mat2 == -4 && (fabs(data->x[l]-xmin) < eps)) elemind[l] = 1; else if(mat2 == -2 && (fabs(data->x[l]-xmax) < eps)) elemind[l] = 1; } if(data->dim >= 3) { if(mat2 == -5 && fabs(data->z[l]-zmin) < eps) elemind[l] = 1; else if(mat2 == -6 && fabs(data->z[l]-zmax) < eps) elemind[l] = 1; } } if(data->dim > 1) { hits = 0; for(k=0;k 1) for(k=0;knoknots;i++) if(boundnodes[i]) *noboundnodes += 1; if(info) printf("Located %d nodes at the interval between boundaries %d and %d\n", *noboundnodes,mat1,mat2); free_Ivector(visited,1,data->noknots); return(0); } int IncreaseElementOrder(struct FemType *data,int info) { int i,j,k,l,side,element,maxcon,totcon,con,newknots,ind,ind2; int noelements,noknots,nonodes,maxnodes,maxelemtype,hit,node; int edge,elemtype; int **newnodetable,inds[2],**newtopo; Real *newx,*newy,*newz; if(info) printf("Trying to increase the element order of current elements\n"); CreateDualGraph(data,FALSE,info); noknots = data->noknots; noelements = data->noelements; maxcon = data->dualmaxconnections; newnodetable = Imatrix(0,maxcon-1,1,noknots); for(i=1;i<=noknots;i++) for(j=0;jdualgraph[j][i]; if(con > i) { newknots++; newnodetable[j][i] = noknots + newknots; } } } if(info) printf("There will be %d new nodes in the elements\n",newknots); newx = Rvector(1,noknots+newknots); newy = Rvector(1,noknots+newknots); if(data->dim == 3) newz = Rvector(1,noknots+newknots); for(i=1;i<=noknots;i++) { newx[i] = data->x[i]; newy[i] = data->y[i]; if(data->dim == 3) newz[i] = data->z[i]; } for(i=1;i<=noknots;i++) { for(j=0;jdualgraph[j][i]; ind = newnodetable[j][i]; if(con && ind) { newx[ind] = 0.5*(data->x[i] + data->x[con]); newy[ind] = 0.5*(data->y[i] + data->y[con]); if(data->dim == 3) newz[ind] = 0.5*(data->z[i] + data->z[con]); } } } maxelemtype = GetMaxElementType(data); if(maxelemtype <= 303) maxnodes = 6; else if(maxelemtype == 404) maxnodes = 8; else if(maxelemtype == 504) maxnodes = 10; else if(maxelemtype == 605) maxnodes = 13; else if(maxelemtype == 706) maxnodes = 15; else if(maxelemtype == 808) maxnodes = 20; else { printf("Not implemented for elementtype %d\n",maxelemtype); bigerror("IncreaseElementOrder: Cant continue the subroutine"); } if(info) printf("New leading elementtype is %d\n",100*(maxelemtype/100)+maxnodes); newtopo = Imatrix(1,noelements,0,maxnodes-1); for(element=1;element<=noelements;element++) { elemtype = data->elementtypes[element]; for(i=0;itopology[element][i]; } for(element=1;element<=data->noelements;element++) { elemtype = data->elementtypes[element]; nonodes = data->elementtypes[element] % 100; for(side=0;;side++) { hit = GetElementGraph(element,side,data,inds); if(!hit) break; if(inds[0] > inds[1]) { ind = inds[1]; ind2 = inds[0]; } else { ind = inds[0]; ind2 = inds[1]; } for(j=0;jdualgraph[j][ind]; if(con == ind2) { node = newnodetable[j][ind]; newtopo[element][nonodes+side] = node; } } } elemtype = 100*(elemtype/100)+nonodes+side; data->elementtypes[element] = elemtype; } DestroyDualGraph(data,info); free_Rvector(data->x,1,data->noknots); free_Rvector(data->y,1,data->noknots); free_Rvector(data->z,1,data->noknots); free_Imatrix(data->topology,1,data->noelements,0,data->maxnodes); free_Imatrix(newnodetable,0,maxcon-1,1,noknots); data->x = newx; data->y = newy; if(data->dim == 3) data->z = newz; data->topology = newtopo; data->noknots += newknots; data->maxnodes = maxnodes; if(info) printf("Increased the element order from 1 to 2\n"); return(0); } int IncreaseElementOrderOld(struct FemType *data,int info) { int i,j,k,l,side,element,noedges,elemtype,newnode; int noelements,noknots,nonodes,nosides,newbound,maxnodes; int maxelementtype,maxedgenodes,elemedges,maxelemedges,edge,dosides; int **edgetable,sideind[MAXNODESD1],sideelemtype,allocated; int *indx,*revindx,*identical,**newtopo; Real *arrange,*newx,*newy,*newz; if(info) printf("Trying to increase the element order of current elements\n"); maxelementtype = 0; noelements = data->noelements; noknots = data->noknots; maxelementtype = GetMaxElementType(data); if(maxelementtype/100 > 4) { printf("IncreaseElementOrder: Implemented only for 2D elements!\n"); dosides = 0; return(1); } if(maxelementtype/100 <= 2) maxedgenodes = 1; else if(maxelementtype/100 <= 4) maxedgenodes = 2; maxelemedges = maxelementtype/100; allocated = FALSE; edgeloop: edge = 0; for(element=1;element<=data->noelements;element++) { elemedges = data->elementtypes[element]/100; for(side=0;side sideind[1]) { edgetable[edge][0] = sideind[0]; edgetable[edge][1] = sideind[1]; } else { edgetable[edge][1] = sideind[0]; edgetable[edge][0] = sideind[1]; } } } } if(!allocated) { noedges = edge; edgetable = Imatrix(1,noedges,0,maxedgenodes+1); for(i=1;i<=noedges;i++) for(j=0;j<=maxedgenodes+1;j++) edgetable[i][j] = 0; allocated = TRUE; goto edgeloop; } printf("There are alltogether %d edges in the elements\n",noedges); arrange = Rvector(1,noedges); for(i=1;i<=noedges;i++) arrange[i] = 0.0; for(i=1;i<=noedges;i++) arrange[i] = edgetable[i][0]; indx = Ivector(1,noedges); SortIndex(noedges,arrange,indx); #if 0 printf("noknots = %d\n",noknots); for(i=1;i<=noknots;i++) printf("indx[%d]=%d edge=%d arrange[%d] = %lg arrange[indx[%d]] = %lg\n", i,indx[i],edgetable[i][0],i,arrange[i],i,arrange[indx[i]]); #endif #if 0 revindx = Ivector(1,data->noknots); for(i=1;i<=noknots;i++) revindx[indx[i]] = i; #endif allocated = FALSE; identical = Ivector(1,noedges); for(i=1;i<=noedges;i++) identical[i] = 0; nosides = 0; for(i=1;i<=noedges;i++) { if(identical[i] < 0) continue; if(maxedgenodes == 1) { for(j=i+1;j<=noedges && edgetable[indx[i]][0] == edgetable[indx[j]][0];j++) identical[j] = -i; } else if(maxedgenodes == 2) { for(j=i+1;j<=noedges && edgetable[indx[i]][0] == edgetable[indx[j]][0];j++) if(edgetable[indx[i]][1] == edgetable[indx[j]][1]) identical[j] = -i; } identical[i] = ++nosides; } printf("There will be %d new nodes in the elements\n",nosides); newx = Rvector(1,noknots+nosides); newy = Rvector(1,noknots+nosides); if(data->dim == 3) { newz = Rvector(1,noknots+nosides); } for(i=1;i<=noknots;i++) { newx[i] = data->x[i]; newy[i] = data->y[i]; if(data->dim == 3) { newz[i] = data->z[i]; } } if(maxelementtype <= 303) maxnodes = 6; else if(maxelementtype == 404) maxnodes = 8; newtopo = Imatrix(1,noelements,0,maxnodes-1); for(element=1;element<=noelements;element++) { elemtype = data->elementtypes[element]; elemedges = elemtype/100; for(i=0;itopology[element][i]; } for(j=1;j<=noedges;j++) { newnode = identical[j]; if(newnode < 0) newnode = identical[abs(newnode)]; if(newnode <= 0) printf("Newnode = %d Edge = %d\n",newnode,j); newnode += noknots; edge = indx[j]; element = edgetable[edge][maxedgenodes]; side = edgetable[edge][maxedgenodes+1]; GetElementSide(element,side,1,data,sideind,&sideelemtype); elemtype = data->elementtypes[element]; newtopo[element][elemtype/100+side] = newnode; if(elemtype == 303) data->elementtypes[element] = 306; else if(elemtype == 404) data->elementtypes[element] = 408; newx[newnode] = 0.5*(data->x[sideind[0]] + data->x[sideind[1]]); newy[newnode] = 0.5*(data->y[sideind[0]] + data->y[sideind[1]]); if(data->dim == 3) { newz[newnode] = 0.5*(data->z[sideind[0]] + data->z[sideind[1]]); } } free_Rvector(data->x,1,data->noknots); free_Rvector(data->y,1,data->noknots); free_Rvector(data->z,1,data->noknots); free_Imatrix(data->topology,1,data->noelements,0,data->maxnodes); data->x = newx; data->y = newy; if(data->dim == 3) { data->z = newz; } data->topology = newtopo; data->noknots += nosides; data->maxnodes = maxnodes; free_Ivector(indx,1,noedges); free_Ivector(identical,1,noedges); free_Imatrix(edgetable,1,noedges,0,maxedgenodes+1); printf("Created extra nodes in the middle of the edges\n"); return(0); } int SideToBulkElements(struct FemType *data,struct BoundaryType *bound, int *sidebulk,int sidetobulk,int info) { int i,j,k,l,elem,ind; int noelements,addelements,noknots,nonodes,currenttype; int **newtopo,*newmaterial,*newelementtypes,maxnodes; int sideind[MAXNODESD1],sideelemtype; int material,maxmaterial,maxtype; noknots = data->noknots; noelements = data->noelements; maxnodes = data->maxnodes; if(info) printf("Creating bulk elements of BC side elements.\n"); addelements = 0; maxtype = 0; maxmaterial = 0; currenttype = 0; for(k=0;sidebulk[k];k++) { for(j=0;j < MAXBOUNDARIES;j++) { if(!bound[j].created) continue; for(i=1; i <= bound[j].nosides; i++) { currenttype = bound[j].types[i]; if(currenttype == sidebulk[k]) addelements++; if(currenttype > maxtype) maxtype = currenttype; } } } if(!addelements) return(0); if(info) printf("%d side elements will be made bulk elements.\n",addelements); if(data->dim == 2 && !sidetobulk) { for(k=0;sidebulk[k];k++) FindEdges(data,bound,sidebulk[k],maxtype+k+1,info); return(0); } if(info) printf("Creating %d elements with types starting from %d.\n", addelements,maxtype+1); noelements += addelements; newtopo = Imatrix(1,noelements,0,maxnodes-1); newmaterial = Ivector(1,noelements); newelementtypes = Ivector(1,noelements); for(i=1;i<=data->noelements;i++) { newmaterial[i] = data->material[i]; if(data->material[i] > maxmaterial) maxmaterial = data->material[i]; newelementtypes[i] = data->elementtypes[i]; nonodes = newelementtypes[i]%100; for(j=0;jtopology[i][j]; } printf("Starting new element materials from %d.\n",maxmaterial+1); elem = data->noelements; material = maxmaterial; for(k=0;sidebulk[k];k++) { material++; for(j=0;j < MAXBOUNDARIES;j++) { if(!bound[j].created) continue; for(i=1; i <= bound[j].nosides; i++) { currenttype = bound[j].types[i]; if(currenttype == sidebulk[k]) { elem++; GetElementSide(bound[j].parent[i],bound[j].side[i],bound[j].normal[i], data,sideind,&sideelemtype); newmaterial[elem] = material; newelementtypes[elem] = sideelemtype; for(ind = 0; ind < sideelemtype%100; ind++) newtopo[elem][ind] = sideind[ind]; } } } } free_Imatrix(data->topology,1,data->noelements,0,data->maxnodes-1); free_Ivector(data->elementtypes,1,data->noelements); free_Ivector(data->material,1,data->noelements); data->noelements = noelements; data->topology = newtopo; data->material = newmaterial; data->elementtypes = newelementtypes; if(info) printf("Moved the side elements to bulk.\n"); for(k=0;sidebulk[k];k++) FindEdges(data,bound,maxmaterial+k+1,maxtype+k+1,info); if(info) printf("Side elements were set to bulk elements\n"); return(0); } static void CylindricalCoordinateTransformation(struct FemType *data,Real r1,Real r2, int rectangle) { int i,j,k,i2,j2,ind1,ind2,nonodes1,nonodes2; Real x,y,r,f,z,q,x2,y2,z2,dx,dy,dz,eps,mult; int hits,trials,tests; int candidates,*candidatelist,*indx; if(rectangle) { printf("Rectangular geometry with r1=%.4lg for %d nodes.\n", r1,data->noknots); } else { printf("Cylindrical geometry with r1=%.4lg r2=%.4lg for %d nodes.\n", r1,r2,data->noknots); } for(i=1;i<=data->noknots;i++) { r = data->x[i]; z = data->y[i]; f = data->z[i]; data->z[i] = z; if(r >= r2) { data->x[i] = cos(f)*r; data->y[i] = sin(f)*r; } else if(r <= r2) { mult = r/r1; if(r > r1) { q = (r-r1)/(r2-r1); r = r1; } else { q = -1.0; } if(f <= 0.25*FM_PI) { data->x[i] = r; data->y[i] = r1*4*(f-0.00*FM_PI)/FM_PI; } else if(f <= 0.75*FM_PI) { data->y[i] = r; data->x[i] = -r1*4*(f-0.5*FM_PI)/FM_PI; } else if(f <= 1.25*FM_PI) { data->x[i] = -r; data->y[i] = -r1*4*(f-1.0*FM_PI)/FM_PI; } else if(f <= 1.75*FM_PI){ data->y[i] = -r; data->x[i] = r1*4*(f-1.5*FM_PI)/FM_PI; } else { data->x[i] = r; data->y[i] = r1*4*(f-2.0*FM_PI)/FM_PI; } if(!rectangle && q > 0.0) { data->x[i] = (1-q)*data->x[i] + q*cos(f)*r2; data->y[i] = (1-q)*data->y[i] + q*sin(f)*r2; } else if(rectangle && mult > 1.0) { data->y[i] *= mult; data->x[i] *= mult; } } /* r <= r2 */ } eps = 1.0e-3*data->minsize; candidates = 0; candidatelist = Ivector(1,data->noknots); indx = Ivector(1,data->noknots); for(i=1;i<=data->noknots;i++) indx[i] = 0; for(j=1;j<=data->noelements;j++) { nonodes1 = data->elementtypes[j]%100; for(i=0;itopology[j][i]; indx[ind2] = ind2; } } for(i=1;i<=data->noknots;i++) { if(!indx[i]) continue; x = data->x[i]; y = data->y[i]; if(fabs(y) > r1+eps) continue; if((fabs(x) > r1+eps) && (fabs(y) > eps) ) continue; if((fabs(x) > eps) && (fabs(y) > eps) && (fabs(fabs(x)-r1) > eps) && (fabs(fabs(y)-r1) > eps)) continue; candidates++; candidatelist[candidates] = i; } printf("%d/%d candidates for duplicate nodes.\n",candidates,data->noknots); hits = tests = trials = 0; for(j=1;j<=candidates;j++) { ind1 = indx[candidatelist[j]]; x = data->x[ind1]; y = data->y[ind1]; z = data->z[ind1]; for(j2=j+1;j2<=candidates;j2++) { ind2 = indx[candidatelist[j2]]; x2 = data->x[ind2]; y2 = data->y[ind2]; z2 = data->z[ind2]; dx = x-x2; dy = y-y2; dz = z-z2; tests++; if(dx*dx + dy*dy + dz*dz < eps*eps) { if(ind2 != ind1) { indx[candidatelist[j2]] = ind1; hits++; } } } } printf("Found %d double nodes in %d tests.\n",hits,tests); for(j=1;j<=data->noelements;j++) { nonodes1 = data->elementtypes[j]%100; for(i=0;itopology[j][i]; if(ind2 != indx[ind2]) { trials++; data->topology[j][i] = indx[ind2]; } } } free_Ivector(indx,1,data->noknots); free_Ivector(candidatelist,1,data->noknots); printf("Eliminated %d nodes from topology.\n",trials); } static void CylindricalCoordinateImprove(struct FemType *data,Real factor, Real r1,Real r2) { int i; Real x,y,r,q,q2,c,cmin,cmax,eps; printf("Cylindrical coordinate for r1=%.4lg and r2=%.4lg.\n",r1,r2); eps = 1.0e-10; cmin = 1.0/(3.0-sqrt(3.)); cmax = 1.0; if(factor > 1.0) factor=1.0; else if(factor < 0.0) factor=0.0; c = cmin+(cmax-cmin)*factor; if(fabs(c-1.0) < eps) return; printf("Improving cylindrical mesh quality r1=%.4lg, r2=%.4lg and c=%.4lg\n",r1,r2,c); for(i=1;i<=data->noknots;i++) { x = data->x[i]; y = data->y[i]; r = sqrt(x*x+y*y); if(r >= r2) continue; if(r < eps) continue; if(fabs(x) <= r1+eps && fabs(y) <= r1+eps) { if(fabs(x) < fabs(y)) { q = fabs(x/y); data->x[i] = (c*q+(1.-q))*x; data->y[i] = (c*q+(1.-q))*y; } else { q = fabs(y/x); data->x[i] = (c*q+(1.-q))*x; data->y[i] = (c*q+(1.-q))*y; } } else { if(fabs(x) < fabs(y)) { q = fabs(x/y); q2 = (fabs(y)-r1)/(r2*fabs(y/r)-r1); data->x[i] = (c*q+(1.-q)) *x*(1-q2) + q2*x; data->y[i] = (c*q+(1.-q)) *(1-q2)*y + q2*y; } else { q = fabs(y/x); q2 = (fabs(x)-r1)/(r2*fabs(x/r)-r1); data->x[i] = (c*q+(1.-q))*(1-q2)*x + q2*x; data->y[i] = (c*q+(1.-q))*y*(1-q2) + q2*y; } } } } static void CylindricalCoordinateCurve(struct FemType *data, Real zet,Real rad,Real angle) { int i; Real x,y,z; Real z0,z1,f,f0,f1,x0,x1,z2,x2,r0; printf("Cylindrical coordinate curve, zet=%.3lg rad=%.3lg angle=%.3lg\n", zet,rad,angle); r0 = rad; f0 = FM_PI*(angle/180.); z0 = zet; z1 = z0+r0*f0; for(i=1;i<=data->noknots;i++) { x = data->x[i]; y = data->y[i]; z = data->z[i]; if(z <= z0) continue; if(z >= z1) { z2 = z0 + sin(f0)*(r0+x) + cos(f0)*(z-z1); x2 = (cos(f0)-1.0)*r0 + cos(f0)*x - sin(f0)*(z-z1); data->z[i] = z2; data->x[i] = x2; } else { f = (z-z0)/r0; z2 = z0 + sin(f)*(r0+x); x2 = (cos(f)-1.0)*r0 + cos(f)*x; data->z[i] = z2; data->x[i] = x2; } } } void SeparateCartesianBoundaries(struct FemType *data,struct BoundaryType *bound,int info) { int i,j,k,l,type,maxtype,addtype,elemsides,totsides,used,hit; int sideelemtype,sideind[MAXBOUNDARIES]; int axissum,axismax,done; Real x,y,z,sx,sy,sz,sxx,syy,szz,dx,dy,dz; Real bclim[MAXBOUNDARIES]; int bc[MAXBOUNDARIES],bcdim[MAXBOUNDARIES]; Real eps=1.0e-4; maxtype = 0; totsides = 0; for(j=0;jx[sideind[l]]; y = data->y[sideind[l]]; z = data->z[sideind[l]]; sx += x; sy += y; sz += z; sxx += x*x; syy += y*y; szz += z*z; } sx /= elemsides; sy /= elemsides; sz /= elemsides; sxx /= elemsides; syy /= elemsides; szz /= elemsides; dx = sqrt(sxx-sx*sx); dy = sqrt(syy-sy*sy); dz = sqrt(szz-sz*sz); if(sideelemtype < 300 && dz < eps) { if(dx < eps * dy) { hit = FALSE; for(i=0;inoboundaries;j++) { if(!bound[j].created) continue; if(!bound[j].nosides) continue; for(i=1;i<=bound[j].nosides;i++) { for(k=1;k<=bound[j].nosides;k++) if(maxtype < bound[j].types[k]) maxtype = bound[j].types[k]; } } printf("Maximum boundary type is %d\n",maxtype); #if 0 for(j=0;jnoboundaries;j++) { if(!bound[j].created) continue; if(!bound[j].nosides) continue; if(bound[j].type) { bound[j].types = Ivector(1,bound[j].nosides); for(k=1;k<=bound[j].nosides;k++) bound[j].types[k] = bound[j].type; bound[j].type = 0; } } #endif for(j=0;jnoboundaries;j++) { if(!bound[j].created) continue; if(!bound[j].nosides) continue; for(k=0;k<4;k++) axishit[k] = 0; done = 0; omstart: for(k=1;k<=bound[j].nosides;k++) { GetElementSide(bound[j].parent[k],bound[j].side[k],bound[j].normal[k], data,sideind,&sideelemtype); sx = sy = sz = 0.0; sxx = syy = szz = 0.0; elemsides = sideelemtype%100; /* Compute the variance within each axis */ for(l=0;lx[sideind[l]]; y = data->y[sideind[l]]; z = data->z[sideind[l]]; sx += x; sy += y; sz += z; sxx += x*x; syy += y*y; szz += z*z; } sx /= elemsides; sy /= elemsides; sz /= elemsides; sxx /= elemsides; syy /= elemsides; szz /= elemsides; dx = sqrt(sxx-sx*sx); dy = sqrt(syy-sy*sy); dz = sqrt(szz-sz*sz); if(dx < eps*dy && dx < eps*dz) { if(sx > 0.0) { if(done) { if(axistype[0]) bound[j].types[k] = maxtype + axistype[0]; } else axishit[0] += 1; } if(sx < 0.0) { if(done) { if(axistype[1]) bound[j].types[k] = maxtype + axistype[1]; } else axishit[1] += 1; } } else if(dy < eps*dx && dy < eps*dz) { if(sy > 0.0) { if(done) { if(axistype[2]) bound[j].types[k] = maxtype + axistype[2]; } else axishit[2] += 1; } if(sy < 0.0) { if(done) { if(axistype[3]) bound[j].types[k] = maxtype + axistype[3]; } else axishit[3] += 1; } } } /* All this is done to select the sidetype appropriately */ if(!done) { axissum = 0; axismax = 0; addtype = 0; for(k=0;k<4;k++) { axissum += axishit[k]; if(axishit[k]) addtype++; } if(axissum) { for(k=0;k<4;k++) { axismax = 0; for(l=0;l<4;l++) { if(axishit[l] > axishit[axismax]) axismax = l; } axistype[axismax] = k+1; axishit[axismax] = -(k+1); } if(axissum == bound[j].nosides) { for(k=0;k<4;k++) axistype[k] -= 1; addtype--; } if(addtype) { printf("Separating %d rectangular boundaries from boundary %d.\n",addtype,j); done = 1; goto omstart; } else done = 0; } } maxtype += addtype; } } void CreateKnotsExtruded(struct FemType *dataxy,struct BoundaryType *boundxy, struct GridType *grid, struct FemType *data,struct BoundaryType *bound, int info) { int i,j,k,l,m,n,knot0,knot1,knot2,elem0,size,kmax,noknots,origtype; int nonodes3d,nonodes2d; int cellk,element,level,side,parent,parent2,layers,elemtype; int material,material2,ind0,ind1,ind2,*indx,*topo; int sideelemtype,sideind[MAXNODESD1],sidetype,maxsidetype,newbounds; int refmaterial1[10],refmaterial2[10],refsidetype[10],indxlength; Real z,*newx,*newy,*newz,corder[3]; Real meanx,meany,absx,absy; if(grid->rotate) SetElementDivisionCylinder(grid,info); else if(grid->dimension == 3) SetElementDivisionExtruded(grid,info); else { printf("CreateKnotsExtruded: unknown option!\n"); return; } InitializeKnots(data); data->dim = 3; origtype = dataxy->elementtypes[1]; if(origtype == 303) elemtype = 706; else if(origtype == 404) elemtype = 808; else if(origtype == 408) elemtype = 820; else if(origtype == 409) elemtype = 827; else { printf("CreateKnotsExtruded: not implemented for elementtypes %d!\n",origtype); return; } printf("Elementtype %d extruded to type %d.\n",origtype,elemtype); nonodes2d = origtype%100; data->maxnodes = nonodes3d = elemtype%100; if(nonodes3d <= 8) layers = 1; else layers = 2; data->noknots = noknots = dataxy->noknots*(layers*grid->totzelems+1); data->noelements = dataxy->noelements * grid->totzelems; data->coordsystem = dataxy->coordsystem; data->numbering = dataxy->numbering; data->noboundaries = dataxy->noboundaries; data->maxsize = dataxy->maxsize; data->minsize = dataxy->minsize; data->partitionexist = FALSE; data->periodicexist = FALSE; data->connectexist = FALSE; maxsidetype = 0; AllocateKnots(data); indxlength = MAX(data->noknots,data->noelements); indx = Ivector(0,indxlength); for(i=0;i<=indxlength;i++) indx[i] = 0; newbounds = 0; if(grid->dimension == 3) newbounds = grid->zcells+1; else if(grid->rotate) { if(grid->rotateblocks < 4) newbounds = 4; if(grid->rotatecartesian) newbounds += grid->rotateblocks; } for(j=0;jnoboundaries+newbounds;j++) { if(boundxy[j].created || j>=data->noboundaries) { bound[j] = boundxy[j]; bound[j].created = TRUE; bound[j].areasexist = FALSE; size = bound[j].nosides = boundxy[j].nosides * grid->totzelems; if(j >= data->noboundaries) size = dataxy->noelements; bound[j].coordsystem = COORD_CART3; bound[j].side = Ivector(1,size); bound[j].side2 = Ivector(1,size); bound[j].material = Ivector(1,size); bound[j].parent = Ivector(1,size); bound[j].parent2 = Ivector(1,size); bound[j].types = Ivector(1,size); bound[j].normal = Ivector(1,size); for(i=1;i<=size;i++) { bound[j].types[i] = 0; bound[j].side[i] = 0; bound[j].side2[i] = 0; bound[j].parent[i] = 0; bound[j].parent2[i] = 0; bound[j].material[i] = 0; bound[j].normal[i] = 1; } } } knot0 = 0; knot1 = layers*dataxy->noknots; if(layers == 2) knot2 = dataxy->noknots; elem0 = 0; level = 0; for(cellk=1;cellk <= grid->zcells ;cellk++) { kmax = grid->zelems[cellk]; for(k=1;k<=kmax; k++) { if(0) printf("elem0=%d knot0=%d knot1=%d\n",elem0,knot0,knot1); level++; for(element=1;element <= dataxy->noelements;element++) { if(dataxy->material[element] < grid->zfirstmaterial[cellk]) continue; if(dataxy->material[element] > grid->zlastmaterial[cellk]) continue; if(grid->rotate) { meanx = 0.0; for(i=0;ix[dataxy->topology[element][i]]; meanx = fabs(meanx/nonodes2d); } if(grid->rotate && meanx < 0.0) continue; if(grid->rotate && cellk%2==0 && meanx < grid->rotateradius1) continue; elem0++; /* Vector telling the new element order. */ indx[(level-1)*dataxy->noelements+element] = elem0; data->elementtypes[elem0] = elemtype; if(grid->zmaterial[cellk]) data->material[elem0] = grid->zmaterial[cellk]; else data->material[elem0] = dataxy->material[element]; if(elemtype == 706) { for(i=0;i<3;i++) { data->topology[elem0][i] = dataxy->topology[element][i]+knot0; data->topology[elem0][i+3] = dataxy->topology[element][i]+knot1; } } else if(elemtype == 808) { for(i=0;i<4;i++) { data->topology[elem0][i] = dataxy->topology[element][i]+knot0; data->topology[elem0][i+4] = dataxy->topology[element][i]+knot1; } } if(elemtype == 820 || elemtype == 827) { for(i=0;i<4;i++) { data->topology[elem0][i] = dataxy->topology[element][i]+knot0; data->topology[elem0][i+4] = dataxy->topology[element][i]+knot1; data->topology[elem0][i+8] = dataxy->topology[element][i+4]+knot0; data->topology[elem0][i+12] = dataxy->topology[element][i]+knot2; data->topology[elem0][i+16] = dataxy->topology[element][i+4]+knot1; } } if(elemtype == 827) { for(i=0;i<4;i++) data->topology[elem0][20+i] = dataxy->topology[element][4+i]+knot2; data->topology[elem0][24] = dataxy->topology[element][8]+knot0; data->topology[elem0][25] = dataxy->topology[element][8]+knot1; data->topology[elem0][26] = dataxy->topology[element][8]+knot2; } } knot0 += layers*dataxy->noknots; knot1 += layers*dataxy->noknots; knot2 += layers*dataxy->noknots; } } data->noelements = elem0; printf("Extruded mesh has %d elements in %d levels.\n",elem0,level); /* Set the element coordinates. */ knot0 = 0; for(cellk=1;cellk <= grid->zcells ;cellk++) { if(cellk == 1) k=0; else k=1; for(;k<=grid->zelems[cellk]; k++) { if(grid->zlinear[cellk]) { z = grid->z[cellk-1] + k*grid->dz[cellk]; } else if(grid->zexpand[cellk] > 0.0) { z = grid->z[cellk-1] + grid->dz[cellk] * (1.- pow(grid->zratios[cellk],(Real)(k))) / (1.-grid->zratios[cellk]); } else if(grid->zelems[cellk] <= 2) { z = grid->z[cellk-1] + k*grid->dz[cellk]; } else { if(k<=grid->zelems[cellk]/2) { z = grid->z[cellk-1] + grid->dz[cellk] * (1.- pow(grid->zratios[cellk],(Real)(k))) / (1.-grid->zratios[cellk]); } else { z = grid->z[cellk] - grid->dz[cellk] * (1.- pow(grid->zratios[cellk],(Real)(grid->zelems[cellk]-k))) / (1.-grid->zratios[cellk]); } } for(i=1;i <= dataxy->noknots;i++) { data->x[i+knot0] = dataxy->x[i]; data->y[i+knot0] = dataxy->y[i]; data->z[i+knot0] = z; } knot0 += layers * dataxy->noknots; } } /* Set the coordinates for the middle nodes in case of quadratic elements. */ if(elemtype == 820 || elemtype == 827) { for(element=1;element <= data->noelements;element++) { topo = data->topology[element]; for(i=0;i<4;i++) { data->x[topo[i+12]] = 0.5*(data->x[topo[i]]+data->x[topo[i+4]]); data->y[topo[i+12]] = 0.5*(data->y[topo[i]]+data->y[topo[i+4]]); data->z[topo[i+12]] = 0.5*(data->z[topo[i]]+data->z[topo[i+4]]); } if(elemtype == 827) { for(i=0;i<4;i++) { data->x[topo[i+20]] = 0.5*(data->x[topo[12+i]]+data->x[topo[12+(i+1)%4]]); data->y[topo[i+20]] = 0.5*(data->y[topo[12+i]]+data->y[topo[12+(i+1)%4]]); data->z[topo[i+20]] = 0.5*(data->z[topo[12+i]]+data->z[topo[12+(i+1)%4]]); } data->x[topo[26]] = 0.5*(data->x[topo[0]]+data->x[topo[6]]); data->y[topo[26]] = 0.5*(data->y[topo[0]]+data->y[topo[6]]); data->z[topo[26]] = 0.5*(data->z[topo[0]]+data->z[topo[6]]); } } } if(grid->rotate) CylindricalCoordinateTransformation(data,grid->rotateradius1, grid->rotateradius2,grid->rotatecartesian); maxsidetype = 0; sidetype = 0; /* Extrude the 2D boundary conditions. */ for(j=0;jnoboundaries;j++) { if(!bound[j].created) continue; side = 0; level = 0; for(cellk=1;cellk <= grid->zcells ;cellk++) { for(k=1;k<=grid->zelems[cellk]; k++) { level++; #if 0 printf("level=%d j=%d side=%d\n",level,j,side); #endif for(i=1;i<=boundxy[j].nosides;i++){ ind1 = (level-1)*dataxy->noelements + boundxy[j].parent[i]; parent = indx[ind1]; if(parent) material = data->material[parent]; else material = 0; if(boundxy[j].parent2[i]) { ind2 = (level-1)*dataxy->noelements + boundxy[j].parent2[i]; parent2 = indx[ind2]; } else parent2 = 0; if(parent2) material2 = data->material[parent2]; else material2 = 0; #if 0 printf("ind=[%d %d] parent=[%d %d] material=[%d %d]\n", ind1,ind2,parent,parent2,material,material2); #endif if((parent || parent2) && (material != material2)) { side++; sidetype = boundxy[j].types[i]; bound[j].types[side] = sidetype; if(parent) { bound[j].parent[side] = parent; bound[j].parent2[side] = parent2; bound[j].side[side] = boundxy[j].side[i]; bound[j].side2[side] = boundxy[j].side2[i]; bound[j].material[side] = material; } else { bound[j].parent[side] = parent2; bound[j].parent2[side] = parent; bound[j].side[side] = boundxy[j].side2[i]; bound[j].side2[side] = boundxy[j].side[i]; bound[j].material[side] = material2; } } } } } bound[j].nosides = side; if(sidetype > maxsidetype) maxsidetype = sidetype; printf("Extruded BC %d of type %d was created with %d sides.\n", j,sidetype,side); } if(grid->layeredbc) { /* Find the BCs between layers. */ if(grid->dimension == 3 || grid->rotatecartesian) { side = 0; level = 0; j--; for(cellk=1;cellk <= grid->zcells ;cellk++) { int swap,redo; redo = FALSE; redolayer: for(k=1;k<=grid->zelems[cellk]; k++) { level++; if(!(k == 1) && !(cellk == grid->zcells && k==grid->zelems[cellk])) continue; if(cellk == grid->zcells && k == grid->zelems[cellk]) { if(grid->zelems[cellk] == 1) redo = TRUE; else level++; } if(grid->rotatecartesian && cellk%2 == 1) continue; if(grid->rotatecartesian && k != 1) continue; for(i=0;i<10;i++) { refmaterial1[i] = 0; refmaterial2[i] = 0; refsidetype[i] = 0; } side = 0; j++; for(i=1;i<=dataxy->noelements;i++){ ind1 = (level-2)*dataxy->noelements+i; if(ind1 < 1) parent = 0; else parent = indx[ind1]; ind2 = (level-1)*dataxy->noelements+i; if(ind2 > indxlength) parent2 = 0; else parent2 = indx[ind2]; if(parent == 0 && parent2 != 0) { parent = parent2; parent2 = 0; swap = 1; } else swap = 0; if(!parent) continue; material = data->material[parent]; if(parent2) material2 = data->material[parent2]; else material2 = 0; #if 0 printf("level=%d ind=[%d %d] parent=[%d %d] material=[%d %d] swap=%d\n", level,ind1,ind2,parent,parent2,material,material2,swap); #endif if(grid->rotatecartesian && !material2) { if(origtype == 303) GetElementSide(parent,4-swap,1,data,sideind,&sideelemtype); else GetElementSide(parent,5-swap,1,data,sideind,&sideelemtype); meanx = meany = 0.0; if(cellk%4 == 2) { for(l=0;ly[sideind[l]]; meany += data->z[sideind[l]]; } } else { for(l=0;lx[sideind[l]]; meany += data->z[sideind[l]]; } } meanx = fabs(meanx)/(sideelemtype%100); meany = fabs(meany)/(sideelemtype%100); if(fabs(meanx - grid->rotateradius1) > 1.0e-12) { material2 = material; } else { for(m=0;mxcells && grid->x[m]+1.0e-12 < meanx;m++); for(n=0;nycells && grid->y[n]+1.0e-12 < meany;n++); material2 = grid->structure[n][m+1]; } #if 0 printf("cellk=%d meanx=%.3lg meany=%.3lg material2=%d m=%d n=%d\n", cellk,meanx,meany,material2,m,n); #endif } if(material != material2) { side++; bound[j].nosides = side; bound[j].parent[side] = parent; bound[j].parent2[side] = parent2; bound[j].material[side] = material; if(origtype == 303) { bound[j].side[side] = 4-swap; bound[j].side2[side] = 3+swap; } else { bound[j].side[side] = 5-swap; bound[j].side2[side] = 4+swap; } for(m=0;m<10;m++) { if(refmaterial1[m] == material && refmaterial2[m] == material2) { break; } else if(refmaterial1[m] == 0 && refmaterial2[m] == 0) { refmaterial1[m] = material; refmaterial2[m] = material2; sidetype++; refsidetype[m] = sidetype; break; } else if(m==9) { printf("Layer includes more than 9 new BCs!\n"); } } bound[j].types[side] = refsidetype[m]; } } printf("BC %d on layer %d was created with %d sides.\n",j,level,side); if(sidetype > maxsidetype) maxsidetype = sidetype; if(redo == TRUE) goto redolayer; } } j++; } } /* Create four additional boundaries that may be used to force symmetry constraints. These are only created if the object is only partially rotated. */ if(grid->rotate && grid->rotateblocks < 4) { k = 0; for(element=1;element<=data->noelements;element++) { int blocks,indexlist[12]={0,1,2,3,4,5,6,7,0,3,7,4}; Real eps; eps = 1.0e-8; blocks = grid->rotateblocks; for(l=0;l<=2;l++) { meanx = meany = 0.0; absx = absy = 0.0; for(i=0;i<4;i++) { ind1 = data->topology[element][indexlist[i+4*l]]; meanx += data->x[ind1]; meany += data->y[ind1]; absx += fabs(data->x[ind1]); absy += fabs(data->y[ind1]); } meanx /= 4.0; meany /= 4.0; absx /= 4.0; absy /= 4.0; /* BCs for zero angle */ if(l==0 && absy < eps && meanx > 0.0) { if(meanx < grid->rotateradius2) { bound[j].nosides++; k = bound[j].nosides; bound[j].side[k] = 4; bound[j].parent[k] = element; bound[j].types[k] = sidetype+1; } else { bound[j+2].nosides++; k = bound[j+2].nosides; bound[j+2].side[k] = 4; bound[j+2].parent[k] = element; bound[j+2].types[k] = sidetype+3; } } /* BCs for angles 90, 180 or 270. */ if(l==1 || l==2) { if(blocks==1 && absx < eps && meany > 0.0) { if(meany < grid->rotateradius2) { bound[j+1].nosides++; k = bound[j+1].nosides; if(l==1) bound[j+1].side[k] = 5; if(l==2) bound[j+1].side[k] = 3; bound[j+1].parent[k] = element; bound[j+1].types[k] = sidetype+2; } else { bound[j+3].nosides++; k = bound[j+3].nosides; if(l==1) bound[j+3].side[k] = 5; if(l==2) bound[j+3].side[k] = 3; bound[j+3].parent[k] = element; bound[j+3].types[k] = sidetype+4; } } if(blocks==2 && absy < eps && meanx < 0.0) { if(meanx > -grid->rotateradius2) { bound[j+1].nosides++; k = bound[j+1].nosides; if(l==1) bound[j+1].side[k] = 5; if(l==2) bound[j+1].side[k] = 3; bound[j+1].parent[k] = element; bound[j+1].types[k] = sidetype+2; } else { bound[j+3].nosides++; k = bound[j+3].nosides; if(l==1) bound[j+3].side[k] = 5; if(l==2) bound[j+3].side[k] = 3; bound[j+3].parent[k] = element; bound[j+3].types[k] = sidetype+4; } } if(blocks==3 && absx < eps && meany < 0.0) { if(meany > -grid->rotateradius2) { bound[j+1].nosides++; k = bound[j+1].nosides; if(l==1) bound[j+1].side[k] = 5; if(l==2) bound[j+1].side[k] = 3; bound[j+1].parent[k] = element; bound[j+1].types[k] = sidetype+2; } else { bound[j+3].nosides++; k = bound[j+3].nosides; if(l==1) bound[j+3].side[k] = 5; if(l==2) bound[j+3].side[k] = 3; bound[j+3].parent[k] = element; bound[j+3].types[k] = sidetype+4; } } } } } printf("Symmetry BCs [%d %d %d %d] have [%d %d %d %d] sides.\n", j,j+1,j+2,j+3,bound[j].nosides,bound[j+1].nosides, bound[j+2].nosides,bound[j+3].nosides); for(l=0;l<4;l++) if(bound[j+l].nosides < 1) bound[j+l].created = FALSE; } data->noboundaries = j+4; #if 0 for(i=0;inoknots;i++) indx[i] = 0; for(element=1;element<=data->noelements;element++) { nonodes3d = data->elementtypes[element]%100; for(i=0;itopology[element][i]] = 1; } j = 0; for(i=1;i<=data->noknots;i++) if(indx[i]) indx[i] = ++j; if(j < data->noknots) { printf("%d original nodes moved to %d new ones.\n",data->noknots,j); newx = Rvector(1,j); for(i=1;i<=data->noknots;i++) newx[indx[i]] = data->x[i]; newy = data->x; data->x = newx; for(i=1;i<=data->noknots;i++) newy[indx[i]] = data->y[i]; newz = data->y; data->y = newy; for(i=1;i<=data->noknots;i++) newz[indx[i]] = data->z[i]; free_Rvector(data->z,1,data->noknots); data->z = newz; data->noknots = j; for(element=1;element<=data->noelements;element++) { nonodes3d = data->elementtypes[element]%100; for(i=0;itopology[element][i] = indx[data->topology[element][i]]; } } if(grid->rotate) { ReorderElements(data,bound,FALSE,corder,info); CylindricalCoordinateImprove(data,grid->rotateimprove, grid->rotateradius1,grid->rotateradius2); if(grid->rotatecurve) CylindricalCoordinateCurve(data,grid->curvezet, grid->curverad,grid->curveangle); if(grid->rotatecartesian) SeparateMainaxisBoundaries(data,bound); printf("Created %d elements and %d nodes by rotation of %d degrees.\n", data->noelements,data->noknots,90*grid->rotateblocks); } else if(grid->dimension == 3) if(info) printf("Created %d elements and %d nodes by extruding the 2D geometry\n", data->noelements,data->noknots); free_Ivector(indx,0,indxlength); } void ReduceElementOrder(struct FemType *data,int matmin,int matmax) /* Reduces the element order at material interval [matmin,matmax] */ { int i,j,element,material,elemcode1,elemcode2,maxnode,*indx,reduced; Real *newx,*newy,*newz; indx = Ivector(0,data->noknots); for(i=0;i<=data->noknots;i++) indx[i] = 0; reduced = 0; for(element=1;element<=data->noelements;element++) { elemcode1 = data->elementtypes[element]; material = data->material[element]; elemcode2 = elemcode1; if(material >= matmin && material <= matmax) elemcode2 = 101*(elemcode1/100); if(elemcode2 == 505) elemcode2 = 504; /* tetrahedron */ #if 0 printf("element=%d codes=[%d,%d]\n",element,elemcode1,elemcode2); printf("mat=%d interval=[%d,%d]\n",material,matmin,matmax); #endif if(elemcode2 < elemcode1) reduced++; maxnode = elemcode2%100; for(i=0;itopology[element][i]] = 1; data->elementtypes[element] = elemcode2; } printf("The element order is reduced in %d elements at interval [%d,%d]\n", reduced,matmin,matmax); j = 0; for(i=1;i<=data->noknots;i++) if(indx[i]) indx[i] = ++j; printf("%d original nodes moved to %d new ones.\n",data->noknots,j); newx = Rvector(1,j); newy = Rvector(1,j); if(data->dim == 3) newz = Rvector(1,j); for(i=1;i<=data->noknots;i++) { newx[indx[i]] = data->x[i]; newy[indx[i]] = data->y[i]; if(data->dim == 3) newz[indx[i]] = data->z[i]; } free_Rvector(data->x,1,data->noknots); free_Rvector(data->y,1,data->noknots); if(data->dim == 3) free_Rvector(data->z,1,data->noknots); data->x = newx; data->y = newy; if(data->dim == 3) data->z = newz; data->noknots = j; for(element=1;element<=data->noelements;element++) { maxnode = data->elementtypes[element]%100; for(i=0;itopology[element][i] = indx[data->topology[element][i]]; } } void MergeElements(struct FemType *data,struct BoundaryType *bound, int manual,Real corder[],Real eps,int mergebounds,int info) { int i,j,k,l; int noelements,noknots,newnoknots,nonodes; int *mergeindx,*doubles; Real *newx,*newy,*newz; Real cx,cy,cz,dx,dy,dz,cdist,dist; ReorderElements(data,bound,manual,corder,TRUE); /* The known ordering by vector corder[] is used to reduce the cost of finding the merged nodes. */ cx = corder[0]; cy = corder[1]; cz = corder[2]; /* Normalizing for future use */ cdist = sqrt(cx*cx+cy*cy+cz*cz); cx /= cdist; cy /= cdist; cz /= cdist; noelements = data->noelements; noknots = data->noknots; newnoknots = noknots; mergeindx = Ivector(1,noknots); for(i=1;i<=noknots;i++) mergeindx[i] = 0; doubles = Ivector(1,noknots); for(i=1;i<=noknots;i++) doubles[i] = 0; if(info) printf("Merging nodes close (%.3lg) to one another.\n",eps); dz = 0.0; for(i=1;ix[i] - data->x[j]; dy = data->y[i] - data->y[j]; if(data->dim == 3) dz = data->z[i] - data->z[j]; if(fabs(cx*dx+cy*dy+cz*dz) > eps) break; dist = dx*dx + dy*dy + dz*dz; if(dist < eps*eps) { doubles[i] = doubles[j] = TRUE; mergeindx[j] = -i; newnoknots--; } } } if(mergebounds) MergeBoundaries(data,bound,doubles,info); j = 0; for(i=1;i<=noknots;i++) if(mergeindx[i] == 0) mergeindx[i] = ++j; for(i=1;i<=noknots;i++) { if(mergeindx[i] < 0) mergeindx[i] = mergeindx[-mergeindx[i]]; } printf("%d original nodes merged to %d new nodes.\n", noknots,newnoknots); newx = Rvector(1,newnoknots); newy = Rvector(1,newnoknots); if(data->dim == 3) newz = Rvector(1,newnoknots); for(i=1;i<=noknots;i++) { newx[mergeindx[i]] = data->x[i]; newy[mergeindx[i]] = data->y[i]; if(data->dim == 3) newz[mergeindx[i]] = data->z[i]; } #if 0 for(i=1;i<=noknots;i++) printf("i=%d indx=%d merge=%d\n",i,indx[i],mergeindx[i]); #endif free_Rvector(data->x,1,data->noknots); free_Rvector(data->y,1,data->noknots); if(data->dim == 3) free_Rvector(data->z,1,data->noknots); data->x = newx; data->y = newy; if(data->dim == 3) data->z = newz; #if 0 if(info) printf("Merging the topologies.\n"); #endif l = 0; for(j=1;j<=noelements;j++) { nonodes = data->elementtypes[j] % 100; for(i=0;itopology[j][i]; data->topology[j][i] = mergeindx[k]; } } data->noknots = newnoknots; free_Ivector(mergeindx,1,noknots); if(info) printf("Merging of nodes is complete.\n"); } void MergeBoundaries(struct FemType *data,struct BoundaryType *bound,int *doubles,int info) { int i,i2,j,k,l,totsides,newsides,sidenodes,sideelemtype,side; int parent,parent2,newind,sideind[MAXNODESD1]; totsides = 0; newsides = 0; if(info) printf("Eliminating boundaries at joined nodes\n"); for(j=0;jnoelements; noknots = data->noknots; x = data->x; y = data->y; if(data->dim == 3) z = data->z; bcindx = Ivector(1,noknots); for(i=1;i<=noknots;i++) bcindx[i] = FALSE; for(j=0;j < MAXBOUNDARIES;j++) { if(!bound[j].created) continue; for(i=1; i <= bound[j].nosides; i++) { element = bound[j].parent[i]; side = bound[j].side[i]; GetElementSide(element,side,1,data,sideind,&sideelemtype); sidenodes = sideelemtype%100; for(k=0;kelementtypes[j]; topo = data->topology[j]; if(elemtype == 306) { for(i=0;i<3;i++) { if(!bcindx[topo[i+3]]) { x[topo[i+3]] = 0.5*(x[topo[i]]+x[topo[(i+1)%3]]); y[topo[i+3]] = 0.5*(y[topo[i]]+y[topo[(i+1)%3]]); } } } else if(elemtype == 310) { for(i=0;i<3;i++) { if(!bcindx[topo[2*i+3]]) { x[topo[2*i+3]] = (2.0*x[topo[i]]+1.0*x[topo[(i+1)%3]])/3.0; x[topo[2*i+4]] = (1.0*x[topo[i]]+2.0*x[topo[(i+1)%3]])/3.0; y[topo[2*i+3]] = (2.0*y[topo[i]]+1.0*y[topo[(i+1)%3]])/3.0; y[topo[2*i+4]] = (1.0*y[topo[i]]+2.0*y[topo[(i+1)%3]])/3.0; } } x[topo[9]] = (x[topo[0]]+x[topo[1]]+x[topo[2]])/3.0; y[topo[9]] = (y[topo[0]]+y[topo[1]]+y[topo[2]])/3.0; } else if(elemtype == 408 || elemtype == 409) { for(i=0;i<4;i++) { if(!bcindx[topo[i+4]]) { x[topo[i+4]] = 0.5*(x[topo[i]]+x[topo[(i+1)%4]]); y[topo[i+4]] = 0.5*(y[topo[i]]+y[topo[(i+1)%4]]); } } if(elemtype == 409) { x[topo[8]] = 0.25*(x[topo[0]]+x[topo[1]]+x[topo[2]]+x[topo[3]]); y[topo[8]] = 0.25*(y[topo[0]]+y[topo[1]]+y[topo[2]]+y[topo[3]]); } } else if(elemtype == 412 || elemtype == 416) { for(i=0;i<4;i++) { if(!bcindx[topo[2*i+4]]) { x[topo[2*i+4]] = (2.0*x[topo[i]]+1.0*x[topo[(i+1)%4]])/3.0; x[topo[2*i+5]] = (1.0*x[topo[i]]+2.0*x[topo[(i+1)%4]])/3.0; y[topo[2*i+4]] = (2.0*y[topo[i]]+1.0*y[topo[(i+1)%4]])/3.0; y[topo[2*i+5]] = (1.0*y[topo[i]]+2.0*y[topo[(i+1)%4]])/3.0; } } if(elemtype == 416) { Real xmean,ymean; xmean = (x[topo[0]]+x[topo[1]]+x[topo[2]]+x[topo[3]])/4.0; ymean = (y[topo[0]]+y[topo[1]]+y[topo[2]]+y[topo[3]])/4.0; for(i=0;i<4;i++) { x[topo[11+i]] = (2.*xmean + 1.0*x[i]) / 3.0; y[topo[11+i]] = (2.*ymean + 1.0*y[i]) / 3.0; } } } else { printf("IsoparamametricElements: Not implemented for elementtype %d\n",elemtype); } } if(info) printf("The elements were forced to be isoparametric\n"); } void ElementsToBoundaryConditions(struct FemType *data, struct BoundaryType *bound,int info) { int i,j,k,l,sideelemtype,sideelemtype2,elemind,elemind2,parent,sideelem,sameelem; int sideind[MAXNODESD1],sideind2[MAXNODESD1],elemsides,side,hit,same,minelemtype; int sidenodes,sidenodes2,maxelemtype,elemtype,elemdim,sideelements,material; int *moveelement,*parentorder,*possible,**invtopo; int noelements,maxpossible,noknots,maxelemsides,twiceelem,sideelemdim; int debug,unmoved; if(info) printf("Making elements to boundary conditions\n"); for(j=0;j < MAXBOUNDARIES;j++) bound[j].created = FALSE; for(j=0;j < MAXBOUNDARIES;j++) bound[j].nosides = 0; noelements = data->noelements; noknots = data->noknots; maxelemtype = GetMaxElementType(data); if(info) printf("Leading bulk elementtype is %d\n",maxelemtype); minelemtype = GetMinElementType(data); if(info) printf("Trailing bulk elementtype is %d\n",minelemtype); elemdim = GetElementDimension(maxelemtype); if( elemdim - GetElementDimension(minelemtype) == 0) return; moveelement = Ivector(1,noelements); sideelements = 0; maxelemtype = 0; maxelemsides = 0; unmoved = 0; for(i=1;i<=noelements;i++) { moveelement[i] = FALSE; elemsides = data->elementtypes[i]/100; if(elemsides > 4) sideelemdim = 3; else if(elemsides > 2) sideelemdim = 2; else if(elemsides > 1) sideelemdim = 1; else if(elemsides == 1) sideelemdim = 0; moveelement[i] = elemdim - sideelemdim; if(moveelement[i]) sideelements++; } if(info) printf("There are %d (out of %d) lower dimensional elements.\n", sideelements,noelements); if(sideelements == 0) return; AllocateBoundary(bound,sideelements); possible = Ivector(1,noknots); for(i=1;i<=noknots;i++) possible[i] = 0; for(elemind=1;elemind <= data->noelements;elemind++) { if(moveelement[elemind]) continue; for(i=0;ielementtypes[elemind]%100;i++) { j = data->topology[elemind][i]; possible[j] += 1; } } j = 1; maxpossible = possible[1]; for(i=1;i<=noknots;i++) if(maxpossible < possible[i]) { maxpossible = possible[i]; j = i; } if(info) printf("Node %d belongs to maximum of %d elements\n",j,maxpossible); /* Make a table showing to which elements a node belongs to Include only the potential parents which are not to be moved to BCs. */ invtopo = Imatrix(1,noknots,1,maxpossible); for(i=1;i<=noknots;i++) for(j=1;j<=maxpossible;j++) invtopo[i][j] = 0; for(elemind=1;elemind <= data->noelements;elemind++) { if(moveelement[elemind]) continue; elemtype = data->elementtypes[elemind]; for(i=0;itopology[elemind][i]; for(l=1;invtopo[k][l];l++); invtopo[k][l] = elemind; } } sideelem = 0; sameelem = 0; twiceelem = 0; debug = FALSE; for(elemind=1;elemind <= data->noelements;elemind++) { if(!moveelement[elemind]) continue; same = FALSE; sideelemtype = data->elementtypes[elemind]; sidenodes = sideelemtype % 100; for(i=0;itopology[elemind][i]; for(l=1;l<=maxpossible;l++) { elemind2 = invtopo[sideind[0]][l]; if(!elemind2) continue; for(side=0;;side++) { if(debug) printf("elem1=%d l=%d elem2=%d side=%d\n",elemind,l,elemind2,side); GetElementSide(elemind2,side,1,data,&sideind2[0],&sideelemtype2); if(sideelemtype2 < 300 && sideelemtype > 300) break; if(sideelemtype2 < 200 && sideelemtype > 200) break; if(sideelemtype != sideelemtype2) continue; sidenodes2 = sideelemtype2 % 100; if(sidenodes != sidenodes2) continue; if(sidenodes2 == 1 && sidenodes > 1) break; hit = 0; for(i=0;iparent2[sideelem] = elemind2; bound->side2[sideelem] = side; goto foundtwo; } else { sideelem += 1; same = TRUE; if(debug) printf("sideelem=%d %d %d\n",sideelem,side,elemind2); bound->parent[sideelem] = elemind2; bound->side[sideelem] = side; bound->parent2[sideelem] = 0; bound->side2[sideelem] = 0; material = data->material[elemind]; bound->types[sideelem] = material; if(sidenodes == 2) { if((sideind[0]-sideind[1])*(sideind2[0]-sideind2[1])<0) bound->normal[sideelem] = -1; } if(data->bodynamesexist) { data->boundarynamesexist = TRUE; if(material < MAXBODIES && material < MAXBOUNDARIES) strcpy(data->boundaryname[material],data->bodyname[material]); if(!strncmp(data->boundaryname[material],"body",4)) strncpy(data->boundaryname[material],"bnry",4); } if(moveelement[elemind] > 1) goto foundtwo; } } } if(!same) { moveelement[elemind] = FALSE; unmoved += 1; } foundtwo: if(0); } if(twiceelem) printf("Found %d sides that were multiply given\n",twiceelem); if(sameelem) printf("Found %d side elements that have two parents.\n",sameelem); bound->nosides = sideelem; if(sideelem == sideelements) { printf("Found correctly %d side elements.\n",sideelem); } else { printf("Found %d side elements, could have found %d\n",sideelem,sideelements); printf("Leaving %d lower dimensional elements to be bulk elements\n",unmoved); } /* Reorder remaining master elements */ parentorder = Ivector(1,noelements); j = 0; for(i=1;i<=noelements;i++) { if(!moveelement[i]) { parentorder[i] = ++j; data->material[j] = data->material[i]; k = data->elementtypes[j] = data->elementtypes[i]; for(l=0;ltopology[j][l] = data->topology[i][l]; } else parentorder[i] = 0; } data->noelements = j; if(info) printf("Parent elements were reordered up to indx %d.\n",j); /* Reorder boundary to point at the new arrangement of master elements */ for(i=1;i<=bound->nosides;i++) { if(bound->parent[i]) bound->parent[i] = parentorder[bound->parent[i]]; if(bound->parent2[i]) bound->parent2[i] = parentorder[bound->parent2[i]]; } if(info) printf("Moved %d elements (out of %d) to new positions\n",j,noelements); free_Ivector(parentorder,1,noelements); free_Ivector(moveelement,1,noelements); free_Ivector(possible,1,noknots); return; } int FindPeriodicNodes(struct FemType *data,int periodicdim[],int info) { int i,j,k,l,i2,j2,dim,jmin; int noelements,noknots,hit,tothits,dimvisited; int *topbot,*indxper; int botn,topn,*indbot,*revindtop,*revindbot; Real eps,dist,dx,dy,dz,coordmax,coordmin,ds; Real *coord,*toparr,*botarr,epsmin; if(data->dim < 3) periodicdim[2] = 0; if(!periodicdim[0] && !periodicdim[1] && !periodicdim[2]) return(1); if(data->periodicexist) { printf("FindPeriodicNodes: Subroutine is called for second time¡\n"); return(2); } noknots = data->noknots; tothits = 0; dimvisited = FALSE; data->periodicexist = TRUE; indxper = Ivector(1,noknots); data->periodic = indxper; for(i=1;i<=noknots;i++) indxper[i] = i; for(dim=1;dim<=3;dim++) { if(!periodicdim[dim-1]) continue; if(info) printf("Finding periodic nodes in dim=%d\n",dim); if(dim==1) coord = data->x; else if(dim==2) coord = data->y; else if(dim==3) coord = data->z; coordmax = coordmin = coord[1]; for(i=1;i<=data->noknots;i++) { if(coordmax < coord[i]) coordmax = coord[i]; if(coordmin > coord[i]) coordmin = coord[i]; } if(info) printf("Coordinate in dimension %d is at the interval [%.3lg, %.3lg]\n", dim,coordmin,coordmax); if(coordmax-coordmin < 1.0e-10) continue; if(!dimvisited) { topbot = Ivector(1,noknots); } eps = 1.0e-5 * (coordmax-coordmin); topn = botn = 0; for(i=1;i<=data->noknots;i++) { if(fabs(coord[i]-coordmax) < eps) { topn++; topbot[i] = topn; } else if(fabs(coord[i] - coordmin) < eps) { botn++; topbot[i] = -botn; } else { topbot[i] = 0; } } if(topn != botn) { printf("There should be equal number of top and bottom nodes (%d vs. %d)!\n",topn,botn); return(3); } else { if(info) printf("Looking for %d periodic nodes\n",topn); } toparr = Rvector(1,topn); botarr = Rvector(1,botn); revindtop = Ivector(1,topn); revindbot = Ivector(1,botn); topn = botn = 0; for(i=1;i<=noknots;i++) { j = topbot[i]; if(j > 0) { topn++; revindtop[topn] = i; } else if(j < 0) { j = abs(j); botn++; revindbot[botn] = i; } } if(data->dim == 2) { for(i=1;i<=botn;i++) { j = revindbot[i]; hit = FALSE; for(i2=1;i2<=topn;i2++) { j2 = revindtop[i2]; if(dim == 1) dist = fabs(data->y[j] - data->y[j2]); else if(dim == 2) dist = fabs(data->x[j] - data->x[j2]); if(dist < eps) { hit = TRUE; goto hit2d; } } hit2d: if(hit) { tothits++; if(indxper[j] == j) indxper[j2] = j; else if(indxper[indxper[j]]==indxper[j]) { #if 0 printf("case2: j=[%d %d] i=[%d %d]\n",j,j2,i,i2); #endif indxper[j2] = indxper[j]; } else { printf("unknown 2d case!\n"); } } else { printf("Couldn't find a periodic counterpart for node %d at [%.3lg %.3lg]]\n", j,data->x[j],data->y[j]); } } } dx = dy = dz = 0.0; if(data->dim == 3) { for(i=1;i<=botn;i++) { j = revindbot[i]; hit = FALSE; epsmin = coordmax - coordmin; for(i2=1;i2<=topn;i2++) { j2 = revindtop[i2]; if(dim == 1) { dy = data->y[j] - data->y[j2]; dz = data->z[j] - data->z[j2]; } else if(dim == 2) { dx = data->x[j] - data->x[j2]; dz = data->z[j] - data->z[j2]; } else if(dim == 3) { dx = data->x[j] - data->x[j2]; dy = data->y[j] - data->y[j2]; } if(dx*dx+dy*dy+dz*dz < eps*eps) { hit = TRUE; goto hit3d; } } hit3d: if(hit) { tothits++; if(indxper[j] == j) indxper[j2] = j; else if(indxper[indxper[j]]==indxper[j]) { indxper[j2] = indxper[j]; } else if(indxper[indxper[indxper[j]]]==indxper[indxper[j]]) { indxper[j2] = indxper[indxper[j]]; } else { printf("unknown 3d case!\n"); } } else { printf("The periodic counterpart for node %d was not found!\n",j); } } } dimvisited = TRUE; } if(info) printf("Found all in all %d periodic nodes.\n",tothits); #if 0 if(data->noknots < 200) { for(i=1;i<=data->noknots;i++) if(i!=indxper[i]) printf("i=%d per=%d\n",i,indxper[i]); } #endif return(0); } int CreateBoundaryLayer(struct FemType *data,struct BoundaryType *bound, int nolayers, int *layerbounds, int *layernumber, Real *layerratios, Real *layerthickness, int *layerparents, int maxfilters, Real layereps, int info) /* Create Boundary layers that may be used to solve accurately fluid flow problems and similar equations. */ { int i,j,k,l,m,n,i2,i3,nonodes,maxbc,newbc; int noknots,noelements,elemindx,nodeindx,elemtype; int oldnoknots,oldnoelements,maxelemtype,oldmaxnodes; int nonewnodes,nonewelements,dolayer,dim,order,midpoints; int checkmaterials,parent,parent2,use2,second; Real dx,dy,dz,ds,ratio,q,p,rectfactor; Real *newx,*newy,*newz,corder[3],*oldx,*oldy,*oldz,*elemwidth; Real e1x,e1y,e2x,e2y; int sideelemtype,sideind[MAXNODESD1],ind[MAXNODESD2],sidebc[MAXNODESD1]; int *layernode,*newelementtypes,**newtopo,**oldtopo; int *topomap,*newmaterial,*herit,*inside,*nonlin; int endbcs, *endparents, *endtypes, *endnodes, *endnodes2, *endneighbours; if(maxfilters == 1) maxfilters = 1000; if(layereps > 0.1) layereps = 0.001; if(layereps < 1.0e-8) layereps = 0.001; rectfactor = 1.0e2; dim = data->dim; maxelemtype = GetMaxElementType(data); if(maxelemtype > 409) { printf("Subroutine implemented only up to 2nd degree!\n"); return(2); } if(info) printf("Largest elementtype is %d\n",maxelemtype); second = FALSE; checkmaterials = FALSE; for(k=0;knoelements; oldnoknots = noknots = data->noknots; oldmaxnodes = data->maxnodes; layernode = Ivector(1,oldnoknots); for(i=1;i<=oldnoknots;i++) layernode[i] = 0; /* Go through all the boundaries with boundary layer definitions and compute the number of new nodes and new elements. */ nonewnodes = 0; nonewelements = 0; maxbc = 0; /* Go through the layers and check which ones are active */ for(j=0;j maxbc) maxbc = bound[j].types[i]; for(k=0;kmaterial[parent] == layerparents[k]) dolayer = k + 1; else if(parent = bound[j].parent2[i]) { if(data->material[parent] == layerparents[k]) { use2 = TRUE; dolayer = k + 1; } } } else { dolayer = k + 1; } } } if(!dolayer) continue; /* We have found an boundary element to extrude */ if(use2) GetElementSide(bound[j].parent2[i],bound[j].side2[i],bound[j].normal[i], data,ind,&sideelemtype); else GetElementSide(parent,bound[j].side[i],bound[j].normal[i], data,ind,&sideelemtype); nonewelements += layernumber[dolayer-1]; midpoints = FALSE; if(sideelemtype == 202) { order = 1; } else if(sideelemtype == 203) { order = 2; if(maxelemtype > 408) midpoints = TRUE; } for(l=0;l 4) { for(j=1;j<=noelements;j++) { elemtype = data->elementtypes[j]; for(i=elemtype/100;itopology[j][i]; layernode[k] = abs(layernode[k]); } } } /* Negative indexed means that the node is an end node of the newly created boundary */ endbcs = 0; for(i=1;i<=noknots;i++) if(layernode[i] < 0) endbcs++; if(endbcs) { endparents = Ivector(1,endbcs); endtypes = Ivector(1,endbcs); endnodes = Ivector(1,endbcs); endnodes2 = Ivector(1,endbcs); endneighbours = Ivector(1,2*endbcs); for(i=1;i<=endbcs;i++) endparents[i] = endtypes[i] = endnodes[i] = endnodes2[i] = 0; endbcs = 0; for(i=1;i<=noknots;i++) { if(layernode[i] < 0) { endbcs++; endparents[endbcs] = i; } } } /* Check if the new boundary is already connected to some one, however it must be different from the extruded boundary */ for(i2=1;i2<=endbcs;i2++) { for(j=0;jmaterial[parent] == layerparents[k]) dolayer = TRUE; else if(parent = bound[j].parent2[i]) { if(data->material[parent] == layerparents[k]) { dolayer = TRUE; } } } } if(!dolayer) continue; /* Finally check that this is not one of the extruded boundaries */ dolayer = FALSE; for(k=0;kx[ind[0]] - data->x[ind[1]]); dy = fabs(data->y[ind[0]] - data->y[ind[1]]); if(dx < rectfactor * dy && dy < rectfactor * dx) { endnodes[i2] = ind[i3]; if(sideelemtype%100 > 2) endnodes2[i2] = ind[2]; endneighbours[2*i2] = ind[1-i3]; } if(info) printf("Found an existing boundary %d for the single node %d %d\n", bound[j].types[i],endparents[i2],endnodes[i2]); goto foundbc; } } foundbc: if(!endtypes[i2]) { maxbc++; endtypes[i2] = maxbc; } } /* Find the first unused bc */ for(j=0;j i) i = layernumber[k]; if(endbcs) { if(info) { printf("Allocating for additional %d boundary elements into bc %d.\n", bound[newbc].nosides,newbc); } AllocateBoundary(&bound[newbc],i*endbcs); bound[newbc].created = FALSE; bound[newbc].nosides = 0; } /* The size of new mesh */ noknots = data->noknots + nonewnodes; noelements = data->noelements + nonewelements; oldnoelements = data->noelements; oldnoknots = data->noknots; if(info) { printf("Creating additional %d elements and %d nodes.\n",nonewelements,nonewnodes); printf("Boundary layer mesh has %d elements and %d nodes.\n",noelements,noknots); } /* there will be more nodes if the original mesh consists of triangles */ if(maxelemtype <= 303) data->maxnodes = 4; else if(maxelemtype == 306) data->maxnodes = 8; /* Allocate more space for the enlarged data set */ newtopo = Imatrix(1,noelements,0,data->maxnodes-1); newmaterial = Ivector(1,noelements); newelementtypes = Ivector(1,noelements); newx = Rvector(1,noknots); newy = Rvector(1,noknots); if(dim == 3) { newz = Rvector(1,noknots); for(i=1;i<=noknots;i++) newz[i] = 0.0; } elemwidth = Rvector(1,nonewelements); for(i=1;i<=nonewelements;i++) elemwidth[i] = 0.0; herit = Ivector(1,noknots); for(i=1;i<=oldnoknots;i++) herit[i] = i; for(i=oldnoknots+1;i<=noknots;i++) herit[i] = 0; /* Set the old topology */ for(j=1;j<=data->noelements;j++) { newmaterial[j] = data->material[j]; newelementtypes[j] = data->elementtypes[j]; for(i=0;ielementtypes[j]%100;i++) newtopo[j][i] = data->topology[j][i]; } /* Set the old nodes */ for(i=1;i<=data->noknots;i++) { newx[i] = data->x[i]; newy[i] = data->y[i]; } topomap = Ivector(1,noknots); for(i=1;i<=noknots;i++) topomap[i] = i; inside = Ivector(1,noelements); for(i=1;i<=noelements;i++) inside[i] = FALSE; /* Set the new node topology and nodes */ elemindx = data->noelements; for(j=0;jmaterial[parent] == layerparents[k]) { dolayer = k + 1; } else if(parent2) { l = parent; parent = parent2; parent2 = l; if(data->material[parent] == layerparents[k]) { use2 = TRUE; dolayer = k + 1; } } } else dolayer = k + 1; } } if(!dolayer) continue; if(use2) GetElementSide(bound[j].parent2[i],bound[j].side2[i],bound[j].normal[i], data,ind,&sideelemtype); else GetElementSide(bound[j].parent[i],bound[j].side[i],bound[j].normal[i], data,ind,&sideelemtype); inside[parent] = 1; if(sideelemtype == 202) order = 1; else if(sideelemtype == 203) order = 2; /* Check if some node should result into additional BC */ for(i2=0;i2x[ind[1]] - data->x[ind[0]]); dx = data->y[ind[1]] - data->y[ind[0]]; ds = sqrt(dx*dx+dy*dy); dx /= ds; dy /= ds; n = layernumber[dolayer-1]; ds = -layerthickness[dolayer-1]; for(l=0;l < n;l++) { elemindx++; newmaterial[elemindx] = data->material[parent]; inside[elemindx] = 1; if(n <= 1 || fabs(layerratios[dolayer-1]-1.0) < 0.001) { q = (1.0*(l+1))/n; elemwidth[elemindx-oldnoelements] = ds / n; } else { ratio = pow(layerratios[dolayer-1],-1./(n-1.)); q = (1.- pow(ratio,(Real)(l+1))) / (1.-pow(ratio,(Real)(n))); p = (1.- pow(ratio,(Real)(l))) / (1.-pow(ratio,(Real)(n))); elemwidth[elemindx-oldnoelements] = (q-p) * ds; } for(m=0;mx[ind[m]]; e2y = newy[nodeindx] - data->y[ind[m]]; p = (e1x*e2x + e1y*e2y)/(sqrt(e1x*e1x+e1y*e1y)*sqrt(e2x*e2x+e2y*e2y)); newx[nodeindx] += e1x - p * e2x; newy[nodeindx] += e1y - p * e2y; } else { herit[nodeindx] = ind[m]; newx[nodeindx] = data->x[ind[m]] + e1x; newy[nodeindx] = data->y[ind[m]] + e1y; } } /* Create the bulk elements */ if(l==0) { newtopo[elemindx][3] = ind[0]; newtopo[elemindx][2] = ind[1]; if(order == 2) newtopo[elemindx][6] = ind[2]; } else { newtopo[elemindx][3] = layernode[ind[0]] + order*l; newtopo[elemindx][2] = layernode[ind[1]] + order*l; if(order == 2) newtopo[elemindx][6] = layernode[ind[2]] + (midpoints+1)*l; } newtopo[elemindx][0] = layernode[ind[0]] + order*(l+1); newtopo[elemindx][1] = layernode[ind[1]] + order*(l+1); if(order == 2) { newtopo[elemindx][7] = layernode[ind[0]] + order*l+1; newtopo[elemindx][5] = layernode[ind[1]] + order*l+1; newtopo[elemindx][4] = layernode[ind[2]] + (midpoints+1)*(l+1); if(midpoints) newtopo[elemindx][8] = layernode[ind[2]] + 2*l+1; } if(order == 1) { newelementtypes[elemindx] = 404; } else if(midpoints) { newelementtypes[elemindx] = 409; } else { newelementtypes[elemindx] = 408; } if(l == n-1 && parent2) { elemtype = data->elementtypes[parent2]; inside[parent2] = 2; for(i2=0;i2topology[parent2][i2] == ind[i3]) { if(i3 < 2) { topomap[ind[i3]] = layernode[ind[i3]] + order * n; } else { topomap[ind[i3]] = layernode[ind[i3]] + (midpoints+1) * n; } } } } } } /* Finally set the BC to point to the new boundary */ if(use2) { bound[j].side2[i] = 0; bound[j].parent2[i] = elemindx; } else { bound[j].side[i] = 0; bound[j].parent[i] = elemindx; } } } { int *inside2; inside2 = Ivector(1,noknots); for(i=1;i<=noknots;i++) inside2[i] = 0; /* Put a marker to all nodes that belong to elements that are on the outside */ for(j=1;j<=noelements;j++) { if(inside[j] == 2) { elemtype = data->elementtypes[j]; for(i=0;ielementtypes[j]; k = 0; for(i=0;i 1) inside[j] = 2; } } free_Ivector(inside2,1,noknots); /* Still, go trough all elements and if they are not on the list of active materials assume them outside */ if(checkmaterials) { for(j=1;j<=oldnoelements;j++) { dolayer = FALSE; for(k=0;kmaterial[j] == layerparents[k]) dolayer = TRUE; if(!dolayer) { if(inside[j] == 1) printf("Element %d of material %d should be in the inside\n", j,data->material[j]); inside[j] = 2; } } } /* And finally remap the nodes that are on the outside */ for(j=1;j<=noelements;j++) { if(inside[j] == 2) { elemtype = data->elementtypes[j]; for(i=0;itopology[j][i]]; } } } /* Put the pointers to the enlarged data set and destroy the old data */ oldx = data->x; oldy = data->y; oldtopo = data->topology; data->noelements = noelements; data->noknots = noknots; data->x = newx; data->y = newy; if(dim == 3) data->z = newz; free_Ivector(data->elementtypes,1,oldnoelements); data->elementtypes = newelementtypes; free_Ivector(data->material,1,oldnoelements); data->material = newmaterial; data->topology = newtopo; /* In case one wants to fit the mesh inside the original mesh the mesh nodes may be put to new positions using an appropiate filter. */ /* For higher order elements remove the middlenodes from the list of cornernodes */ if(maxelemtype%100 > 4) { if(info) printf("Marking the higher order nodes\n"); nonlin = Ivector(1,noknots); for(i=1;i<=noknots;i++) nonlin[i] = FALSE; for(j=1;j<=noelements;j++) { elemtype = data->elementtypes[j]; for(i=elemtype/100;itopology[j][i]; nonlin[k] = TRUE; } } } if(maxfilters) { int hit,method,iter; int ind1,ind2,ind3,*fixedx,*fixedy; Real *aidx,*aidy,*weights; Real maxerror,minds,dx2,dy2,ds2,fii; /* There are three methods how to put the weight in the filter, 1) 1/s, 2) fii/s, 3) sin(fii)/s, the second option seems to be best. */ method = 2; if(info) printf("Filtering the mesh to meet the original geometry\n"); fixedx = Ivector(1,noknots); fixedy = Ivector(1,noknots); weights = Rvector(1,noknots); aidx = Rvector(1,noknots); aidy = Rvector(1,noknots); /* Set all the fixed boundaries */ for(i=1;i<=noknots;i++) fixedx[i] = fixedy[i] = 0; /* First, make all other materials except the ones with BL to be fixed */ if(checkmaterials) { for(j=1;j<=noelements;j++) { elemtype = data->elementtypes[j]; dolayer = FALSE; for(k=0;kmaterial[j] == layerparents[k]) dolayer = TRUE; for(i=0;itopology[j][i]; if(dolayer && fixedx[ind1]%2 == 0) fixedx[ind1] += 1; if(!dolayer && fixedx[ind1] < 2) fixedx[ind1] += 2; } } for(i=1;i<=noknots;i++) { if(fixedx[i] == 2) fixedy[i] = 2; else fixedx[i] = 0; } } /* Then set all BC:s fixed except the tangential ones */ for(j=0;j rectfactor * dy) { for(l=0;l rectfactor * dx) { for(l=0;lmaterial[parent] == layerparents[k]) { dolayer = k + 1; } else if(parent2) { l = parent; parent = parent2; parent2 = l; if(data->material[parent] == layerparents[k]) { use2 = TRUE; dolayer = k + 1; } } } else dolayer = k + 1; } } if(dolayer) { for(l=0;lelementtypes[j]; nonodes = elemtype / 100; for(i=0;itopology[j][i]; ind2 = data->topology[j][i2]; ind3 = data->topology[j][i3]; if(j<=oldnoelements) { dx = oldx[oldtopo[j][i2]] - oldx[oldtopo[j][i]]; dy = oldy[oldtopo[j][i2]] - oldy[oldtopo[j][i]]; ds = sqrt(dx*dx+dy*dy); } else { ds = fabs(elemwidth[j-oldnoelements]); } if(ds < minds) minds = ds; if(j<=oldnoelements) { dx2 = oldx[oldtopo[j][i2]] - oldx[oldtopo[j][i3]]; dy2 = oldy[oldtopo[j][i2]] - oldy[oldtopo[j][i3]]; ds2 = sqrt(dx2*dx2+dy2*dy2); } else { ds2 = fabs(elemwidth[j-oldnoelements]); } if(j <= oldnoelements && ds * ds2 < 1.0e-50) { printf("problem elem %d and nodes %d (%d %d)\n",j,i2,i,i3); printf("dist ds=%.3le ds2=%.3le\n",ds,ds2); printf("coord: %.3le %.3le\n",oldx[oldtopo[j][i2]], oldy[oldtopo[j][i2]]); continue; } if(abs(method) == 2 && j<=oldnoelements) { fii = acos((dx*dx2+dy*dy2)/(ds*ds2)) / (FM_PI/2.0); } else if(abs(method) == 3 && j<=oldnoelements) { fii = acos((dx*dx2+dy*dy2)/(ds*ds2)); fii = sin(fii); } else { fii = 1.0; } /* Eliminate the very difficult triplenodes */ dolayer = FALSE; for(k=1;k<=endbcs;k++) if(ind2 == endnodes[k]) dolayer = k; if(dolayer) { for(k=1;k<=2;k++) { if(endneighbours[2*(dolayer-1)+k] == ind1) { weights[ind2] += fii / ds; if(!fixedx[ind2]) newx[ind2] += aidx[ind1] * fii / ds; if(!fixedy[ind2]) newy[ind2] += aidy[ind1] * fii / ds; } } for(k=1;k<=2;k++) { if(endneighbours[2*(dolayer-1)+k] == ind3) { weights[ind2] += fii / ds2; if(!fixedx[ind2]) newx[ind2] += aidx[ind3] * fii / ds2; if(!fixedy[ind2]) newy[ind2] += aidy[ind3] * fii / ds2; } } } else { if(ind2 <= oldnoknots || herit[ind1] == herit[ind2]) { weights[ind2] += fii / ds; if(!fixedx[ind2]) newx[ind2] += aidx[ind1] * fii / ds; if(!fixedy[ind2]) newy[ind2] += aidy[ind1] * fii / ds; } if(ind2 <= oldnoknots || herit[ind3] == herit[ind2]) { weights[ind2] += fii / ds2; if(!fixedx[ind2]) newx[ind2] += aidx[ind3] * fii / ds2; if(!fixedy[ind2]) newy[ind2] += aidy[ind3] * fii / ds2; } } } } if(maxelemtype%100 > 4) { for(j=1;j<=noknots;j++) { if(nonlin[j]) continue; if(weights[j] > 1.0e-50) { if(!fixedx[j]) newx[j] /= weights[j]; if(!fixedy[j]) newy[j] /= weights[j]; } else if(iter==1) { printf("no weight for index %d\n",j); } dx = newx[j] - aidx[j]; dy = newy[j] - aidy[j]; ds = dx*dx + dy*dy; if(ds > maxerror) maxerror = ds; } } else { for(j=1;j<=noknots;j++) { if(!fixedx[j]) newx[j] /= weights[j]; if(!fixedy[j]) newy[j] /= weights[j]; dx = newx[j]-aidx[j]; dy = newy[j]-aidy[j]; ds = dx*dx + dy*dy; if(ds > maxerror) maxerror = ds; } } maxerror = sqrt(maxerror) / minds; if(maxerror < layereps) break; } if(info) { printf("Filtered the new node coordinates %d times with final error %.3le.\n", iter-1,maxerror); } /* In higher order elements map the middle nodes so that they lie in between the corner nodes */ if(maxelemtype%100 > 4) { for(j=1;j<=noelements;j++) { elemtype = data->elementtypes[j]; if(elemtype%100 <= elemtype/100) continue; if(elemtype == 306) { for(k=0;k<3;k++) { if(!fixedx[newtopo[j][k+3]]) { newx[newtopo[j][k+3]] = 0.5 * (newx[newtopo[j][k]] + newx[newtopo[j][(k+1)%3]]); } if(!fixedy[newtopo[j][k+3]]) { newy[newtopo[j][k+3]] = 0.5 * (newy[newtopo[j][k]] + newy[newtopo[j][(k+1)%3]]); } } } else if(elemtype == 408 || elemtype == 409) { if(elemtype == 409) { newx[newtopo[j][8]] = 0.0; newy[newtopo[j][8]] = 0.0; } for(k=0;k<4;k++) { if(!fixedx[newtopo[j][k+4]]) { newx[newtopo[j][k+4]] = 0.5 * (newx[newtopo[j][k]] + newx[newtopo[j][(k+1)%4]]); } if(!fixedy[newtopo[j][k+4]]) { newy[newtopo[j][k+4]] = 0.5 * (newy[newtopo[j][k]] + newy[newtopo[j][(k+1)%4]]); } if(elemtype == 409) { newx[newtopo[j][8]] += 0.25 * newx[newtopo[j][k]]; newy[newtopo[j][8]] += 0.25 * newy[newtopo[j][k]]; } } } else { printf("Unknown elementtype %d\n",elemtype); } } } free_Ivector(fixedx,1,noknots); free_Ivector(fixedy,1,noknots); free_Rvector(aidx,1,noknots); free_Rvector(aidy,1,noknots); free_Rvector(weights,1,noknots); } if(bound[newbc].nosides > 0) bound[newbc].created = TRUE; /* In higher order elements map the middle nodes so that they lie in between the corner nodes. Elemtypes must be 408 or 409 since they are created in this subroutine */ if(!maxfilters && maxelemtype%100 > 4) { if(info) printf("Making the higher order nodes to lie in between\n"); for(j=oldnoelements+1;j<=noelements;j++) { elemtype = data->elementtypes[j]; if(elemtype%100 <= elemtype/100) continue; if(elemtype == 408 || elemtype == 409) { if(elemtype == 409) { newx[newtopo[j][8]] = 0.0; newy[newtopo[j][8]] = 0.0; } for(k=0;k<4;k++) { newx[newtopo[j][k+4]] = 0.5 * (newx[newtopo[j][k]] + newx[newtopo[j][(k+1)%4]]); newy[newtopo[j][k+4]] = 0.5 * (newy[newtopo[j][k]] + newy[newtopo[j][(k+1)%4]]); if(elemtype == 409) { newx[newtopo[j][8]] += 0.25 * newx[newtopo[j][k]]; newy[newtopo[j][8]] += 0.25 * newy[newtopo[j][k]]; } } } } } #if 0 ReorderElements(data,bound,FALSE,corder,info); #endif free_Imatrix(oldtopo,1,oldnoelements,0,oldmaxnodes-1); free_Ivector(layernode,1,oldnoknots); free_Rvector(oldx,1,oldnoknots); free_Rvector(oldy,1,oldnoknots); if(info) printf("Boundary layers created successfully.\n"); if(checkmaterials && !second) { for(k=0;kdim; maxelemtype = GetMaxElementType(data); if(maxelemtype > 409) { printf("Subroutine implemented only up to 2nd degree!\n"); return(2); } if(info) printf("Largest elementtype is %d\n",maxelemtype); oldnoelements = noelements = data->noelements; oldnoknots = noknots = data->noknots; oldmaxnodes = data->maxnodes; layernode = Ivector(1,oldnoknots); for(i=1;i<=oldnoknots;i++) layernode[i] = 0; sharednode = Ivector(1,oldnoknots); for(i=1;i<=oldnoknots;i++) sharednode[i] = 0; /* Go through all the boundaries with boundary layer definitions and compute the numbder of nodes at the surface. */ maxbc = 0; /* Go through the layers and check which ones are active */ for(j=0;j maxbc) maxbc = bound[j].types[i]; for(k=0;k 408) midpoints = TRUE; } for(l=0;l 1.0 || slayer < 1.0e-20) slayer = 1.0; bcnodes = 0; maxcon = 0; for(i=1;i<=data->noknots;i++) { if(layernode[i]) bcnodes++; maxcon = MAX(maxcon, layernode[i]); } if(info) printf("Found %d new nodes in the boundary layers!\n",bcnodes); if(!bcnodes) return(0); if(info) printf("There are %d connections at maximum\n",maxcon); /* there will be more nodes if the original mesh consists of triangles */ if(maxelemtype <= 303) data->maxnodes = 4; else if(maxelemtype == 306) data->maxnodes = 8; /* Compute the number of new elements */ nonewelements = 0; for(j=1;j<=data->noelements;j++) { elemhits = 0; elemtype = data->elementtypes[j]; for(i=0;itopology[j][i]; if( layernode[k]) { sharednode[k] += 1; elemhits++; } } if(elemhits) { nonewelements += nlayer ; if(elemhits != 2) nonewelements += nlayer + 1; } } printf("There will %d new elemenets\n",nonewelements); /* This is a conservative estimate */ nonewnodes = 2*nonewelements; edgepairs = Imatrix(1,nonewnodes,1,3); for(j=1;j<=nonewnodes;j++) edgepairs[j][1] = edgepairs[j][2] = edgepairs[j][3] = 0; /* The size of new mesh */ oldnoelements = data->noelements; oldnoknots = data->noknots; oldtopo = data->topology; oldx = data->x; oldy = data->y; if(dim == 3) oldz = data->z; noknots = oldnoknots + nonewnodes; noelements = oldnoelements + nonewelements; if(info) { printf("Creating additional %d elements and %d nodes.\n",nonewelements,nonewnodes); printf("Boundary layer mesh has at maximum %d elements and %d nodes.\n",noelements,noknots); } /* Allocate more space for the enlarged data set */ newtopo = Imatrix(1,noelements,0,data->maxnodes-1); newmaterial = Ivector(1,noelements); newelementtypes = Ivector(1,noelements); newx = Rvector(1,noknots); newy = Rvector(1,noknots); if(dim == 3) newz = Rvector(1,noknots); /* Set the old topology */ for(j=1;j<=data->noelements;j++) { newmaterial[j] = data->material[j]; newelementtypes[j] = data->elementtypes[j]; for(i=0;ielementtypes[j]%100;i++) newtopo[j][i] = data->topology[j][i]; } /* Set the old nodes */ for(i=1;i<=data->noknots;i++) { newx[i] = data->x[i]; newy[i] = data->y[i]; if(dim == 3) newz[i] = data->z[i]; } noelements = data->noelements; elemind = noelements; noknots = data->noknots; /* Go through elements and make the new elements and nodes */ for(j=1;j<=data->noelements;j++) { elemhits = 0; elemtype = data->elementtypes[j]; elemsides = elemtype % 100; for(i=0;itopology[j][i] ]) elemhits++; if(!elemhits) continue; if(elemtype == 404) { elemdone = FALSE; for(side=0;sidex = newx; data->y = newy; data->topology = newtopo; data->material = newmaterial; data->elementtypes = newelementtypes; data->noknots = noknots; data->noelements = elemind; printf("The created boundary layer mesh has at %d elements and %d nodes.\n",noelements,noknots); return(0); } int RotateTranslateScale(struct FemType *data,struct ElmergridType *eg,int info) { int i,j,k; Real x,y,z,xz,yz,yx,zx,zy,xy,cx,cy,cz; Real xmin, xmax, ymin, ymax, zmin, zmax; if(eg->scale) { if(info) printf("Scaling mesh with vector [%.3lg %.3lg %.3lg]\n", eg->cscale[0],eg->cscale[1],eg->cscale[2]); for(i=1;i<=data->noknots;i++) { data->x[i] *= eg->cscale[0]; data->y[i] *= eg->cscale[1]; if(data->dim == 3) data->z[i] *= eg->cscale[2]; } if(0) printf("Scaling of mesh finished.\n"); } if(eg->rotate) { if(info) printf("Rotating mesh with degrees [%.3lg %.3lg %.3lg]\n", eg->crotate[0],eg->crotate[1],eg->crotate[2]); cx = FM_PI * eg->crotate[0]/180.0; cy = FM_PI * eg->crotate[1]/180.0; cz = FM_PI * eg->crotate[2]/180.0; for(i=1;i<=data->noknots;i++) { x = data->x[i]; if(data->dim >= 2) y = data->y[i]; else y = 0.0; if(data->dim >= 3) z = data->z[i]; else z = 0.0; xz = x*cos(cz) + y*sin(cz); yz = -x*sin(cz) + y*cos(cz); if(data->dim == 3) { yx = yz*cos(cx) + z*sin(cx); zx = -yz*sin(cx) + z*cos(cx); zy = zx*cos(cy) + xz*sin(cy); xy = -zx*sin(cy) + xz*cos(cy); data->x[i] = xy; data->y[i] = yx; data->z[i] = zy; } else { data->x[i] = xz; data->y[i] = yz; } } if(0) printf("Rotation of mesh finished.\n"); } if(eg->translate) { if(info) printf("Translating the mesh with vector [%.3lg %.3lg %.3lg]\n", eg->ctranslate[0],eg->ctranslate[1],eg->ctranslate[2]); for(i=1;i<=data->noknots;i++) { data->x[i] += eg->ctranslate[0]; data->y[i] += eg->ctranslate[1]; if(data->dim == 3) data->z[i] += eg->ctranslate[2]; } if(0) printf("Translation of mesh finished.\n"); } if(eg->center) { xmin = xmax = data->x[1]; ymin = ymax = data->y[1]; if(data->dim == 3) zmin = zmax = data->z[1]; for(i=1;i<=data->noknots;i++) { xmax = MAX( xmax, data->x[i] ); xmin = MIN( xmin, data->x[i] ); ymax = MAX( ymax, data->y[i] ); ymin = MIN( ymin, data->y[i] ); if(data->dim == 3) { zmax = MAX( zmax, data->z[i] ); zmin = MIN( zmin, data->z[i] ); } } cx = 0.5 * (xmin + xmax); cy = 0.5 * (ymin + ymax); if(data->dim == 3) cz = 0.5 * (zmin + zmax); if(info) printf("Setting new center to %.3le %.3le %.3le\n",cx,cy,cz); for(i=1;i<=data->noknots;i++) { data->x[i] -= cx; data->y[i] -= cy; if(data->dim == 3) data->z[i] -= cz; } } return(0); } int CreateDualGraph(struct FemType *data,int full,int info) { int i,j,k,l,m,totcon,noelements, noknots,elemtype,nonodes,hit,ind,ind2; int maxcon,percon,edge; printf("Creating a dual graph of the finite element mesh\n"); if(data->dualexists) { printf("The dual graph already exists! You shoule remove the old graph!\n"); } maxcon = 0; totcon = 0; percon = 0; noelements = data->noelements; noknots = data->noknots; for(i=1;i<=noelements;i++) { elemtype = data->elementtypes[i]; if(!full) { int inds[2]; for(edge=0;;edge++) { if( !GetElementGraph(i,edge,data,&inds[0]) ) break; ind = inds[0]; ind2 = inds[1]; hit = FALSE; for(l=0;ldualgraph[l][ind] == ind2) hit = TRUE; if(data->dualgraph[l][ind] == 0) break; } if(!hit) { if(l >= maxcon) { data->dualgraph[maxcon] = Ivector(1,noknots); for(m=1;m<=noknots;m++) data->dualgraph[maxcon][m] = 0; maxcon++; } data->dualgraph[l][ind] = ind2; totcon++; } /* Make also so symmetric connection */ for(l=0;ldualgraph[l][ind2] == ind) hit = TRUE; if(data->dualgraph[l][ind2] == 0) break; } if(!hit) { if(l >= maxcon) { data->dualgraph[maxcon] = Ivector(1,noknots); for(m=1;m<=noknots;m++) data->dualgraph[maxcon][m] = 0; maxcon++; } data->dualgraph[l][ind2] = ind; totcon++; } } } else { nonodes = data->elementtypes[i] % 100; for(j=0;jtopology[i][j]; for(k=0;ktopology[i][k]; if(ind == ind2) continue; hit = FALSE; for(l=0;ldualgraph[l][ind] == ind2) hit = TRUE; if(data->dualgraph[l][ind] == 0) break; } if(!hit) { if(l >= maxcon) { data->dualgraph[maxcon] = Ivector(1,noknots); for(m=1;m<=noknots;m++) data->dualgraph[maxcon][m] = 0; maxcon++; } data->dualgraph[l][ind] = ind2; totcon++; } } } } } if( data->periodicexist ) { for(ind=1;ind<=noknots;ind++) { ind2 = data->periodic[ind]; if(ind == ind2) continue; hit = FALSE; for(l=0;ldualgraph[l][ind] == ind2) hit = TRUE; if(data->dualgraph[l][ind] == 0) break; } if(!hit) { if(l >= maxcon) { data->dualgraph[maxcon] = Ivector(1,noknots); for(m=1;m<=noknots;m++) data->dualgraph[maxcon][m] = 0; maxcon++; } data->dualgraph[l][ind] = ind2; totcon++; percon++; } } } data->dualmaxconnections = maxcon; data->dualexists = TRUE; if(info) printf("There are at maximum %d connections in dual graph.\n",maxcon); if(info) printf("There are at all in all %d connections in dual graph.\n",totcon); if(info && percon) printf("There are %d periodic connections in dual graph.\n",percon); return(0); } int DestroyDualGraph(struct FemType *data,int info) { int i,maxcon, noknots; if(!data->dualexists) { printf("You tried to destroy a non-existing dual graph\n"); return(1); } maxcon = data->dualmaxconnections; noknots = data->noknots; for(i=0;idualgraph[i],1,noknots); data->dualmaxconnections = 0; data->dualexists = FALSE; if(info) printf("The dual graph was destroyed\n"); return(0); } int CreateInverseTopology(struct FemType *data,int info) { int i,j,k,l,m,totcon,noelements,noknots,elemtype,nonodes,hit,ind,maxcon,percon; int *neededby,minneeded,maxneeded; printf("Creating an inverse topology of the finite element mesh\n"); if(data->invtopoexists) { printf("The inverse topology already exists!\n"); smallerror("The inverse topology not done"); } maxcon = 0; noelements = data->noelements; noknots = data->noknots; neededby = Ivector(1,noknots); for(i=1;i<=noknots;i++) neededby[i] = 0; for(i=1;i<=noelements;i++) { elemtype = data->elementtypes[i]; nonodes = data->elementtypes[i] % 100; for(j=0;jtopology[i][j]; neededby[ind] += 1; l = neededby[ind]; if(l > maxcon) { maxcon++; data->invtopo[maxcon] = Ivector(1,noknots); for(m=1;m<=noknots;m++) data->invtopo[maxcon][m] = 0; } data->invtopo[l][ind] = i; } } minneeded = maxneeded = neededby[1]; for(i=1;i<=noknots;i++) { minneeded = MIN( minneeded, neededby[i]); maxneeded = MAX( maxneeded, neededby[i]); } free_Ivector(neededby,1,noknots); if(info) printf("There are from %d to %d connections in the inverse topology.\n",minneeded,maxneeded); data->invtopoexists = TRUE; data->maxinvtopo = maxcon; return(0); } int MeshTypeStatistics(struct FemType *data,int info) { int i,elemtype,maxelemtype,minelemtype,*elemtypes; maxelemtype = minelemtype = data->elementtypes[1]; for(i=1;i<=data->noelements;i++) { elemtype = data->elementtypes[i]; maxelemtype = MAX( maxelemtype, elemtype ); minelemtype = MIN( minelemtype, elemtype ); } elemtypes = Ivector(minelemtype,maxelemtype); for(i=minelemtype;i<=maxelemtype;i++) elemtypes[i] = 0; for(i=1;i<=data->noelements;i++) { elemtype = data->elementtypes[i]; elemtypes[elemtype] += 1; } if(info) { printf("Number of different elementtypes\n"); for(i=minelemtype;i<=maxelemtype;i++) if(elemtypes[i]) printf("\t%d\t%d\n",i,elemtypes[i]); } free_Ivector(elemtypes,minelemtype,maxelemtype); return(0); }