ibgdwasy.c

有限元学习研究用源代码(老外的),供科研人员参考

/* last edit: Ilja Schmelzer -------------- 20-OCT-1994 11:01:23.12	*/
/************************************************************************/
/*                                                                      */
/*  <<< I B G >>> - Intersection - Based Grid generation package 	*/
/*                                                                      */
/*  Version 1.1 by Ilja Schmelzer   [e-mail]       */
/*                                                                      */
/*  to be distributed under IBG license conditions (see "readme.ibg")	*/
/*                                                                      */
/************************************************************************/
/* <<< IBGDWASY >>> - Intersection-Based Geometry Description

	        Using the input of the WASY GmbH

We use here the spline package "ibgdSpline.h". So, this file can be used as an
example of a geometry definition using a set of splines to describe the
boundary.

The spline package was developed to allow the usage of existent geometry
descriptions with boundary grids or splines in IBGG. In general, we don't
recommend to use boundary splines for geometry description. Even a small
mismatch in the boundary definition or rounding errors can lead to global
errors in the  geometry description so that the grid generation is not very
stable. We recommend to use small random shifts if such problems occur. */

#include "ibg.h"
#include "ibgi.h"
#include "ibglib.h"
#include "ibgg.h"
#include "ibgd.h"
#include "ibgd0.h"
#include "ibgdspline.h"

#define big1 1.e6

ibGeometry ibgdPolygons(char* filename, ibgSegment outside,
        	        ibgFloat *xmin, ibgFloat *xmax,
        	        ibgFloat *ymin, ibgFloat *ymax)
{ibgdSplines *s;
 ibgdSplineLine *s2;
 FILE *file;
#define buflen 200
 char buffer[buflen];
 int i, last, rc, regs=0, segm, reg, ireg=0, num=0;
 ibgFloat *xx,*yy,x0,y0,xb,yb,x1,y1;
 ibgFloat area[ibgSMAX], a, epsx,epsy;
 ibgSegment *sl,*sr;
 ibgInt *cc,oldp,firstp,nextp,search;
 file = fopen(filename,"r");
 while(1){
        fgets(buffer,buflen,file);
        rc = sscanf(buffer,"%d",&reg);
        if(rc != 1) goto ende1;
        regs++; segm=0; a=0;
        while(1){
            fgets(buffer,buflen,file);
            rc = sscanf(buffer,"%f %f",&x1, &y1);
            if(rc != 2){
                if(segm){
                    num++;
                }
                break;
            }
            if(segm){
                num++;
                a += x0*y1 - x1*y0;
            }else{xb = x1; yb = y1; segm = 1;}
            x0 = x1; y0 = y1;
        }
        a += x1*yb - xb*y1;
        area[ireg] = a;
        ibgprintf("region %d: area %f\n",reg,a);
        ireg++;
 }
ende1:
 ibgprintf("geometry with %d regions and %d polygon segments\n",regs,num);
 fclose(file);
 file = fopen(filename,"r");
 *xmin=*xmax=xb; *ymin=*ymax=yb;
 num += 5;
 xx = (ibgFloat *) malloc(num*sizeof(ibgFloat));
 yy = (ibgFloat *) malloc(num*sizeof(ibgFloat));
 sl = (ibgSegment *)  malloc(num*sizeof(ibgSegment));
 sr = (ibgSegment *)  malloc(num*sizeof(ibgSegment));
 cc = (ibgInt *)  malloc(num*sizeof(ibgInt));
 last=0; sl[last] = 0;
 while(1){
        fgets(buffer,buflen,file);
        rc = sscanf(buffer,"%d",&reg);
        if(rc != 1)	    goto ende2;
        regs++;
        oldp=firstp=nextp=0;
        while(1){
            fgets(buffer,buflen,file);
            rc = sscanf(buffer,"%f %f",&x1, &y1);
            if(rc != 2){
                if(firstp!=oldp){
                    nextp = firstp;
		    if((xx[firstp]==xx[oldp])&&(yy[firstp]==yy[oldp])) break;
	            x1 = xx[firstp];
                    y1 = yy[firstp];
                }else{
                    nextp = 0;
		    break;
                }
            }
            search = 1;
          testold:
            if (oldp){
	        oldp--;
                if((x1==xx[oldp]) && (y1==yy[oldp]) && (sl[oldp])){
                    sr[oldp] = reg; search = 0;
                }else{
	            oldp++;
                    if(cc[oldp] > oldp) {
                        oldp = cc[oldp]; goto testold;
                    }else{
                        if(last > oldp){
                            last++;
                            xx[last] = xx[oldp]; yy[last] = yy[oldp];
                            cc[last] = cc[oldp]; cc[oldp] = last; oldp = last;
                        }
	                last++;
	                xx[last] = x1; yy[last] = y1; cc[last] = last;
	                sl[oldp] = reg; sr[oldp] = sl[last] = sr[last] = 0;
	                oldp = last;
	            }
	        }
	    }
	    if(search){
	        if(nextp){
	            cc[last] = oldp = nextp;
	            while(cc[nextp]>nextp) nextp=cc[nextp];
	            cc[nextp] = last;
	            break;
	        }
	        for(i=1;i<last;i++){
                    if((x1==xx[i]) && (y1==yy[i])){ /* found!! */
	                if(oldp){
	                    cc[last] = oldp = i; while(cc[i]>i) i=cc[i];
	                    cc[i] = last;
	        	}else{
	        	    firstp = oldp = i;
	        	}
	        	break;
	            }
	        }
                if(i>=last){
	            if(oldp){
	            }else{
	                last++;
	                xx[last] = x1; yy[last] = y1; cc[last] = last;
                        sl[last] = sr[last] = 0;
	                firstp = oldp = last;
	            }
	        }

            }
            if(*xmax<x1) *xmax=x1;
            if(*ymax<y1) *ymax=y1;
            if(*xmin>x1) *xmin=x1;
            if(*ymin>y1) *ymin=y1;
            if(rc!=2){
/*	        for(i=1;i<=last;i++){
	            printf("(%f,%f) %d %d %d\n",xx[i],yy[i],sl[i],sr[i],cc[i]);
	        }
		*/
	        break;
	    }
        }
 }
ende2:
 fclose(file);
/* creation of the spline list, space allocation for the splines: */
 s  = ibgdSplineListNew(2,last+10,sizeof(ibgdSplineLine),ibgdSplineLineTest);
/* this is the array of splines we have to initialize now: */
 s2 = (ibgdSplineLine*) s->spline;
 num = 0;
 epsx = 0.273529e-6 * (*xmax-*xmin);
 epsy = 0.192164e-6 * (*ymax-*ymin);
 for(i=1;i<=last;i++){
     if(sr[i]==0) sr[i] = outside;
     printf("(%f,%f) %d %d %d\n",xx[i],yy[i],sl[i],sr[i],cc[i]);
     if(sl[i]==0) continue;
     if(sl[i]==sr[i]) continue;
     ibgdSplineLineNew(&s2[num],(xx[i]+epsx)/1e3,(yy[i]+epsy)/1e3,
        	                (xx[i+1]+epsx)/1e3,(yy[i+1]+epsy)/1e3,
		                 sl[i],sr[i],0);
     num++;
 }
 *xmax /= 1e3;
 *xmin /= 1e3;
 *ymax /= 1e3;
 *ymin /= 1e3;
 s->num = num;
 s->outside = outside;
/* Now we create the geometry description defined by the spline list.
Remark that this routine assumes that the splines in the list are already
initialized. For example, a search tree for fast search may be created in this
function. */

 {ibGeometry geo = ibgdSplineList(s);
  geo->Delta *= (*xmax - *xmin);
  return geo;
 }
}