ibg2simplex.c

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

 }else if(ibgridGDIM(*gg)==2){
	cpoints = 3*rcells + 2*fcells + 1;
	csides = 3*rcells + 4*fcells + 1;
	rtyp = ibgc2Triangle;	rpoints = 3;
	ftyp = ibgc2Edge;	fpoints = 2;
 }else if(ibgridGDIM(*gg)==3){
	cpoints = 4*rcells + 3*fcells + 1;
	csides = 4*rcells + 5*fcells + 1;
	rtyp = ibgc3Tetrahedron;rpoints = 4;
     	ftyp = ibgc3Triangle;	fpoints = 3;
 }		     
 gg->freeCPoint	= cpoints-1;
 gg->freeCSide	= csides-1;
 gg->maxCPoint	= cpoints;
 gg->maxCSide	= csides;
 gg->CellType	= malloc((cells+1)*sizeof(ibgCellType));
 gg->CellSegment        = malloc((cells+1)*sizeof(ibgSegment));
 gg->CPointFirst = malloc((cells+1)*sizeof(int));
 gg->CSideFirst = malloc((cells+1)*sizeof(int));
 gg->CPoint	= malloc((cpoints+1)*sizeof(int));
 gg->CSide	= malloc((csides+1)*sizeof(int));
 gg->lastCOut	= gg->maxCOut = 0;
/* reading point information: */
 for(n=1;n<=points;n++){
	fgets(buffer,buflen,file);
	if(ibgridGDIM(*gg)==1){
	    	rc = sscanf(buffer,"%f",&x); y=0; z=0;
		if(rc != 1) goto readerror;
	}else if(ibgridGDIM(*gg)==2){
	    	rc = sscanf(buffer,"%f %f",&x,&y); z=0;
		if(rc != 2) goto readerror;
	}else if(ibgridGDIM(*gg)==3){
	    	rc = sscanf(buffer,"%f %f %f",&x,&y,&z);
		if(rc != 3) goto readerror;
	}
	point = &ibgridPoint(*gg,n);
	if(x<xmin)	xmin = x;
	if(y<ymin)	ymin = y;
	if(z<zmin)	zmin = z;
	if(x>xmax)	xmax = x;
	if(y>ymax)	ymax = y;
	if(z>zmax)	zmax = z;
	ibgpX(*point)[0] = x;
	ibgpX(*point)[1] = y;
	ibgpX(*point)[2] = z;
 }
 ibgridXMin(*gg,0) = xmin;
 ibgridXMin(*gg,1) = ymin;
 ibgridXMin(*gg,2) = ymin;
 ibgridXMax(*gg,0) = xmax;
 ibgridXMax(*gg,1) = ymax;
 ibgridXMax(*gg,2) = zmax;
 ccn=0; ccs=0;
 for(c=1;c<=rcells;c++){	/* at first the rcells region cells: */
	fgets(buffer,buflen,file);
	rc = sscanf(buffer,"%d %d %d %d %d %d %d %d %d ",
		&num[0],&num[1],&num[2],&num[3],
		&num[4],&num[5],&num[6],&num[7],&num[8]);
	ibgridCellType(*gg,c) = rtyp;
	gg->CPointFirst[c] = ccn;
	gg->CSideFirst[c] = ccs;
	for(i=0;i<rpoints;i++){
		gg->CPoint[ccn++] = num[i];	/* at first the points */
	}
	ibgridCellSegment(*gg,c) = num[i++];	/* now the segment */
	for(j=0;j<rpoints;i++,j++){		/* now the neighbour cells: */
			/* positive -> neighbour is a region cell: */
		if(num[i]>0)		gg->CSide[ccs++] = num[i];
			/* negative -> neighbour is a boundary face cell */
		else if(num[i]<0)	gg->CSide[ccs++] = rcells - num[i];
		else 			gg->CSide[ccs++] = 0;
	}
 }
/* now the boundary face cells from rcells + 1 to rcells + fcells */
 for(cb=1;cb<=fcells;c++,cb++){
	for(i=0;i<9;i++) num[i]=0;
	fgets(buffer,buflen,file);
	rc = sscanf(buffer,"%d %d %d %d %d %d %d %d %d ",
		&num[0],&num[1],&num[2],&num[3],
		&num[4],&num[5],&num[6],&num[7],&num[8]);
	gg->CellType[c] = ftyp;
	ibgassert(gg->CellType[c] == ftyp);
	gg->CPointFirst[c] = ccn;	/* begin of point list */
	gg->CSideFirst[c] = ccs;	/* begin of neighbour list */
	for(i=0;i<fpoints;i++){
	    	gg->CPoint[ccn++] = num[i];
	}
	u  = num[i++];		/* segment  number */
	cl = num[i++];		/* left  neighbour (cell or segment number) */
	cr = num[i++];		/* right neighbour (cell or segment number) */
		/* neighbour cells in the boundary face grid: */
	for(j=0;j<fpoints;j++){
		gg->CSide[ccs++] = rcells+num[i++];
	}
	gg->CSide[ccs++] = cl;
	gg->CSide[ccs++] = cr;
	if(u==0){ /* automatical generation of boundary face segment number */
	    if(cl>0) ul = ibgridCellSegment(*gg,cl);
	    else     ul = -cl;
	    if(cr>0) ur = ibgridCellSegment(*gg,cr);
	    else     ur = -cr;
	    u = ibgDefaultFaceNr(ul,ur);
	}
	ibgridCellSegment(*gg,c) = u;
 }
 return gg;
readerror:
 ibGridFree(gg);
 return ibgNULL;
}

int ibgToSimplexU(ibGrid *gg, ibgSegmentType typ, ibgSegment u0,
				ibgBoolean bound, FILE *file)
{int gdim,d,ic,cc,u,in,nn,n,*cs,*cn,*nnum,*cnum,i,nmax,creg,cfac,rc;
 ibgCellType t,ts;
 ibgPoint *point;
 rc = ibgSuccess;
 fprintf(file,ibgNformatSimplex);
 fprintf(file," ");
 fprintf(file,ibgNversionSimplex);
 fprintf(file,"\n");
 fprintf(file,"\n");
 fprintf(file,"%d dimensional grid \n",gdim=gg->gridDimension);
 cnum = calloc((gg->lastCell+1),sizeof(int));
 nnum = calloc((gg->lastPoint+1),sizeof(int));
 creg = 0;
 if(typ == ibgSRegion){
	foribgRegionCells(*gg,ic,t,u){
		if(u != u0) continue;
 		cnum[ic] = ++creg;
		cn = ibgridCPointList(*gg,ic);
		for(n=0;n<ibgcPoints(t);n++){
			nnum[cn[n]] = 1;
		}
	}endibgRegionCells(*gg,ic,t,u);
 }
 cfac = 0;
 if(typ == ibgSFace){
  	foribgFaceCells(*gg,ic,t,u){
	      	if(u != u0) continue;
		cnum[ic] = ++cfac;
		cn = ibgridCPointList(*gg,ic);
		for(n=0;n<ibgcPoints(t);n++){
			nnum[cn[n]] = 1;
		}
	}endibgFaceCells(*gg,ic,t,u);
 }else if(bound){
  	foribgFaceCells(*gg,ic,t,u){
		cc = ibgridCellLeft(*gg,ic,t);
	 	if(u0 != ibgridCellSegment(*gg,cc)){
			cc = ibgridCellRight(*gg,ic,t);
		 	if(u0 != ibgridCellSegment(*gg,cc)) continue;
		}
		cnum[ic] = ++cfac;
		cn = ibgridCPointList(*gg,ic);
		for(n=0;n<ibgcPoints(t);n++){
			nnum[cn[n]] = 1;
		}
	}endibgFaceCells(*gg,ic,t,u);
 }
 nmax = 0;
 foribgPoints(*gg,in,t,u,point){
 	if(nnum[in]) nnum[in] = ++nmax;
 }endibgPoints(*gg,in,t,u,point);
 fprintf(file,"%d points \n",nmax);
 fprintf(file,"%d cells \n",creg);
 fprintf(file,"%d boundary cells \n",cfac);
 foribgPoints(*gg,in,t,u,point){
       	if(nnum[in]==0) continue;
	for(d=0;d<gdim;d++) fprintf(file,"%f ",(float)(ibgpX(*point)[d]));
	fprintf(file,"\n");
 }endibgPoints(*gg,in,t,u,point);
 if(creg){
	foribgRegionCells(*gg,ic,t,u){
		if(cnum[ic]==0) continue;
		if(gg->gridDimension==2)
			{if(t!=ibgc2Triangle)	{rc=ibgError;continue;}}
		else if(gg->gridDimension==3)
			{if(t!=ibgc3Tetrahedron){rc=ibgError;continue;}}
		cn = ibgridCPointList(*gg,ic);
		for(i=0;i<ibgcPoints(t);i++){
			nn = cn[i];
			ibgassert(nn>0);
			nn = nnum[nn];
			ibgassert(nn>0);
		 	fprintf(file,"%d ",nn);
		}
	 	fprintf(file,"%d ",u);
		cs = ibgridCSideList(*gg,ic);
		for(i=0;i<ibgcSides(t);i++){
			cc = cs[i]; ts = ibgridCellType(*gg,cc);
			if(ibgcCODIM(ts)==ibgSRegion){
			 	fprintf(file,"%d ",cnum[cc]);
			}else if(ibgcCODIM(ts)==ibgSFace){
			 	fprintf(file,"-%d ",cnum[cc]);
			}else	ibgfatal;
		}
		fprintf(file,"\n");
	}endibgRegionCells(*gg,ic,t,u);
 }
 if(cfac){
     	foribgFaceCells(*gg,ic,t,u){
		if(cnum[ic]==0) continue;
		if(gg->gridDimension==2)
			{if(t!=ibgc2Edge)	{rc=ibgError;continue;}}
		else if(gg->gridDimension==3)
			{if(t!=ibgc3Triangle)	{rc=ibgError;continue;}}
		cn = ibgridCPointList(*gg,ic);
		for(i=0;i<ibgcPoints(t);i++){
			nn = cn[i];		/* nn is number in grid */
			ibgassert(nn>0);
			nn = nnum[nn];		/* nn is now number in file */
			ibgassert(nn>0);
		 	fprintf(file,"%d ",nn);
		}
/* user-defined face segment number: */
		if(u < ibgSMAX)	fprintf(file,"%d ",u);
/* automatically created face segment number (of type left*ibgSMAX + right segment) */
		else	       	fprintf(file,"%d ",0);
/* inner (region) neighbours (left and right) of the face cell:
		cc = ibgridCellLeft(*gg,ic,t);
	 	fprintf(file,"%d ",ibgridCellSegment(*gg,cc));
		cc = ibgridCellRight(*gg,ic,t);
	 	fprintf(file,"%d ",ibgridCellSegment(*gg,cc));
/* neighbours of the face cell in the boundary face grid */
		cs = ibgridCSideList(*gg,ic);
		for(i=0;i<ibgcSides(t);i++){
			cc = cs[i]; ts = ibgridCellType(*gg,cc);
			if(ibgcCODIM(ts)==ibgSFace){
			 	fprintf(file,"%d ",cnum[cc]);
			}else if(ibgcCODIM(ts)==ibgSLine){
	/* we have a boundary line cell - the boundary of a boundary face */
			 	fprintf(file,"0 ");
			}else	ibgfatal;
		}
		fprintf(file,"\n");
	}endibgFaceCells(*gg,ic,t,u);
 }
 free(nnum);
 free(cnum);
 return rc;
}