ibggmake.c

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

 }
 cxx = cnx(ns); xs = cxx[0]-x;   ys = cxx[1]-y;   zs = cxx[2]-z;
 cxx = cnx(n0); x0 = cxx[0]-x;   y0 = cxx[1]-y;   z0 = cxx[2]-z;
 cxx = cnx(n1); x1 = cxx[0]-x;   y1 = cxx[1]-y;   z1 = cxx[2]-z;
 cxx = cnx(n2); x2 = cxx[0]-x;   y2 = cxx[1]-y;   z2 = cxx[2]-z;
 det0= x0*(y1*z2-z1*y2) + y0*(z1*x2-x1*z2) + z0*(x1*y2-y1*x2);
 vol0= ibggVRel*det0 - ibggVAbs;
 det = x0*(y1*zs-z1*ys) + y0*(z1*xs-x1*zs) + z0*(x1*ys-y1*xs);
 if(det > vol0){
	if(g3segment1edge++ >= 1) return 0;
	if(g3EDelete(cb,ibgcSEdge(tb,sb,0)))	return 2;
	return 0;
 }
 det = x1*(y2*zs-z2*ys) + y1*(z2*xs-x2*zs) + z1*(x2*ys-y2*xs);
 if(det > vol0){
	if(g3segment1edge++ >= 1) return 0;
	if(g3EDelete(cb,ibgcSEdge(tb,sb,1)))	return 2;
	return 0;
 }
 det = x2*(y0*zs-z0*ys) + y2*(z0*xs-x0*zs) + z2*(x0*ys-y0*xs);
 if(det > vol0){
	if(g3segment1edge++ >= 1) return 0;
	if(g3EDelete(cb,ibgcSEdge(tb,sb,2)))	return 2;
	return 0;
 }
 for(i=0;i<ibgcSSize(ts,ss);i++){
	if(n0==(nh=ccn(cs)[sh=ibgcSPoint1(ts,ss,i)])) ss0=sh;
	else if(n1==nh) ss1=sh;
	else {ibgassert(n2==nh); ss2=sh;}
 }
 g0CGet(&cn,ibgc3Tetrahedron); if(cn<0) return 0;
 ccn(cn)[0]   = nb;
 ccn(cn)[1]   = n1;
 ccn(cn)[2]   = n2;
 ccn(cn)[3]   = ns;
 ccs(cn)[1]   = cs;
 ccs(cn)[2]   = cb;
 ccn(cb)[sb2] = ns;
 ccn(cs)[ss1] = nb;
 co=ccs(cb)[sb1]; ccs(cb)[sb1]=cs; ccs(cs)[ss]= co;
 for(i=0;i<ibgcSides(cct(co));i++) if(ccs(co)[i]==cb) {ccs(co)[i]=cs;break;}
 co=ccs(cb)[sb0]; ccs(cb)[sb0]=cn; ccs(cn)[3] = co;
 for(i=0;i<ibgcSides(cct(co));i++) if(ccs(co)[i]==cb) {ccs(co)[i]=cn;break;}
 co=ccs(cs)[ss2]; ccs(cs)[ss2]=cb; ccs(cb)[sb]= co;
 for(i=0;i<ibgcSides(cct(co));i++) if(ccs(co)[i]==cs) {ccs(co)[i]=cb;break;}
 co=ccs(cs)[ss0]; ccs(cs)[ss0]=cn; ccs(cn)[0] = co;
 for(i=0;i<ibgcSides(cct(co));i++) if(ccs(co)[i]==cs) {ccs(co)[i]=cn;break;}
 if(ibgcRegion(cb)==ibgNotFound)if(st[cb]==0){st[cb] = nosegment0; nosegment0 = cb;}
 if(ibgcRegion(cs)==ibgNotFound)if(st[cs]==0){st[cs] = nosegment0; nosegment0 = cs;}
 if(ibgcRegion(cn)==ibgNotFound){
	if(cn>=mst){mst=(3*mst)/2; st1=(int*)realloc(st,mst*sizeof(int));
 		if(st1!=ibgNULL) st=st1; else ibgfatal;
	}
	if(st[cn]==0){st[cn] = nosegment0; nosegment0 = cn;}
 }
 g3cvol(cb);	g3cvol(cs);	g3cvol(cn);
 return 1;
}

int g34Swap(int c1, int c2, int c3, int c4)
{int i,co,s12,s13,s1u,s1o,s21,s24,s2u,s2o,s31,s34,s3u,s3o,s42,s43,s4u,s4o;
 int u1,u2,u3,u4,n1,n2,n3,n4,nu,no;
 ibgFloat *cxx;
 ibgFloat x,y,z,x2,y2,z2,x3,y3,z3,x4,y4,z4,xu,yu,zu,xo,yo,zo,det,det0,vol0;
 ibgCellType t1=cct(c1);
 for(i=0;i<4;i++){
	if(ccs(c1)[i]==c2) {n1=ccn(c1)[i];s12=i;}
	if(ccs(c2)[i]==c1) {n2=ccn(c2)[i];}
	if(ccs(c1)[i]==c3) {n3=ccn(c1)[i];s13=i;}
	if(ccs(c3)[i]==c1) {n4=ccn(c3)[i];}
 }
 if(ibgIncorrect == ibgiStatusOfEdge(ibGeo,cnd(n1),cnd(n2))) return 0;
 for(i=0;i<ibgcSSize(t1,s12);i++){
	if(ibgcSSide(t1,s12,i)==s13){
		nu = ccn(c1)[ibgcSPoint1(t1,s12,i)];
		no = ccn(c1)[ibgcSPoint1(t1,s12,i+1)];
		break;
	}
 }
 cxx = cnx(n1); x  = cxx[0];	 y  = cxx[1];	  z  = cxx[2];
 cxx = cnx(n2); x2 = cxx[0]-x;   y2 = cxx[1]-y;   z2 = cxx[2]-z;
 cxx = cnx(nu); xu = cxx[0]-x;   yu = cxx[1]-y;   zu = cxx[2]-z;
 cxx = cnx(n4); x4 = cxx[0]-x;   y4 = cxx[1]-y;   z4 = cxx[2]-z;
 cxx = cnx(no); xo = cxx[0]-x;   yo = cxx[1]-y;   zo = cxx[2]-z;
 cxx = cnx(n3); x3 = cxx[0]-x;   y3 = cxx[1]-y;   z3 = cxx[2]-z;
 vol0= -ibggVAbs;
 det = x2*(yu*z4-zu*y4) + y2*(zu*x4-xu*z4) + z2*(xu*y4-yu*x4);
 if(det > vol0) return 0;
 det = x2*(y4*zo-z4*yo) + y2*(z4*xo-x4*zo) + z2*(x4*yo-y4*xo);
 if(det > vol0) return 0;
 det = x2*(yo*z3-zo*y3) + y2*(zo*x3-xo*z3) + z2*(xo*y3-yo*x3);
 if(det > vol0) return 0;
 det = x2*(y3*zu-z3*yu) + y2*(z3*xu-x3*zu) + z2*(x3*yu-y3*xu);
 if(det > vol0) return 0;
 for(i=0;i<4;i++){
	if(ccn(c1)[i]==nu) s1u=i;
	if(ccn(c1)[i]==no) s1o=i;
	if(ccn(c2)[i]==n2) s21=i;
	if(ccn(c2)[i]==n3) s24=i;
	if(ccn(c2)[i]==nu) s2u=i;
	if(ccn(c2)[i]==no) s2o=i;
	if(ccn(c3)[i]==n1) s34=i;
	if(ccn(c3)[i]==n4) s31=i;
	if(ccn(c3)[i]==nu) s3u=i;
	if(ccn(c3)[i]==no) s3o=i;
	if(ccn(c4)[i]==n2) s43=i;
	if(ccn(c4)[i]==n4) s42=i;
	if(ccn(c4)[i]==nu) s4u=i;
	if(ccn(c4)[i]==no) s4o=i;
 }
 ccn(c1)[s1o] = n2;
 ccn(c2)[s2u] = n1;
 ccn(c3)[s3o] = n2;
 ccn(c4)[s4u] = n1;
 co=ccs(c1)[s1u]; ccs(c1)[s1u] = c2; ccs(c2)[s21] = co;
 for(i=0;i<ibgcSides(cct(co));i++) if(ccs(co)[i]==c1) {ccs(co)[i]=c2;break;}
 co=ccs(c2)[s2o]; ccs(c2)[s2o] = c1; ccs(c1)[s12] = co;
 for(i=0;i<ibgcSides(cct(co));i++) if(ccs(co)[i]==c2) {ccs(co)[i]=c1;break;}
 co=ccs(c3)[s3u]; ccs(c3)[s3u] = c4; ccs(c4)[s43] = co;
 for(i=0;i<ibgcSides(cct(co));i++) if(ccs(co)[i]==c3) {ccs(co)[i]=c4;break;}
 co=ccs(c4)[s4o]; ccs(c4)[s4o] = c3; ccs(c3)[s34] = co;
 for(i=0;i<ibgcSides(cct(co));i++) if(ccs(co)[i]==c4) {ccs(co)[i]=c3;break;}
 if(ibgcRegion(c1)==ibgNotFound)if(st[c1]==0){st[c1] = nosegment0; nosegment0 = c1;}
 if(ibgcRegion(c2)==ibgNotFound)if(st[c2]==0){st[c2] = nosegment0; nosegment0 = c2;}
 if(ibgcRegion(c3)==ibgNotFound)if(st[c3]==0){st[c3] = nosegment0; nosegment0 = c3;}
 if(ibgcRegion(c4)==ibgNotFound)if(st[c4]==0){st[c4] = nosegment0; nosegment0 = c4;}
 g3cvol(c1);	g3cvol(c2);	g3cvol(c3);	g3cvol(c4);
 return 1;
}

static int g3EDelete(int cb, int eb)
{ibgCellType tb=cct(cb),tl,tr;
 int ub,ul,cl,cr,co,sbl,sbr,sbo,sbu,slb,slr,slo,slu,sl4,sr4,srb,srl,sro,sru;
 int nu,no,nl,nr,ns,i;
 ibgFloat *cxx;
 ibgFloat x,y,z,xl,yl,zl,xr,yr,zr,xo,yo,zo,xu,yu,zu,det1,det2,vol0;
 int rc=0;
 int c4=0;
beg:
 cl=ccs(cb)[sbl=ibgcESideL(tb,eb)]; if(ccdef(cl)) return 0; tl=cct(cl);
 cr=ccs(cb)[sbr=ibgcESideR(tb,eb)]; if(ccdef(cr)) return 0; tr=cct(cr);
 nu=ccn(cb)[sbu=ibgcEPoint1(tb,eb)];
 no=ccn(cb)[sbo=ibgcEPoint2(tb,eb)];
 nl=ccn(cb)[sbl];
 nr=ccn(cb)[sbr];
 for(i=0;i<ibgcSides(tl);i++){
	if(ccs(cl)[i]==cb)	{slb=i; ns=ccn(cl)[slb];}
	else if(ccn(cl)[i]==nr)	{sl4=i;
		if(ccs(cl)[i]!=cr) {
/* four or more tetrahedrons aroung the given edge: */
		      	for(srb=0;srb<4;srb++) if(ccs(cr)[srb]==cb) break;
			if(g3SDelete(cb,sbl)){
				cb=ccs(cr)[srb];
				for(eb=0;eb<6;eb++){
				    if(ccn(cb)[ibgcEPoint1(tb,eb)]==nu)
					if(ccn(cb)[ibgcEPoint2(tb,eb)]==no)
						break;
				    if(ccn(cb)[ibgcEPoint1(tb,eb)]==no)
					if(ccn(cb)[ibgcEPoint2(tb,eb)]==nu)
						break;
				}
				ibgassert(eb<6);
				rc=1;
    				goto beg;
			}
		      	for(slb=0;slb<4;slb++) if(ccs(cl)[slb]==cb) break;
			if(g3SDelete(cb,sbr)){
				cb=ccs(cl)[slb];
				for(eb=0;eb<6;eb++){
				    if(ccn(cb)[ibgcEPoint1(tb,eb)]==nu)
					if(ccn(cb)[ibgcEPoint2(tb,eb)]==no)
						break;
				    if(ccn(cb)[ibgcEPoint1(tb,eb)]==no)
					if(ccn(cb)[ibgcEPoint2(tb,eb)]==nu)
						break;
				}
				ibgassert(eb<6);
				rc=1;
    				goto beg;
			}
/* special case: edge with 4 Tetrahedrons */
			c4=ccs(cl)[sl4];
			for(sr4=0;sr4<ibgcSides(tr);sr4++){
				if(ccn(cr)[sr4]==nl)	{
					if(ccs(cr)[sr4]!=c4){
/* five or more tetrahedrons around the edge: */
				      		if(g3SDelete(cr,sr4)) return 2;
				      		if(g3SDelete(cl,sl4)) return 2;
    						return rc;
					}
					break;
				}
			}
	 		if(ccout(cl)){
				if(ccdef(c4)) return rc;
				if(ccout(cr)) return rc;
				if(g34Swap(cb,cr,cl,c4)) return 1;
				return rc;
			}
	 		if(ccout(cr)){
				if(ccdef(c4)) return rc;
				if(ccout(cl)) return rc;
				if(g34Swap(cb,cl,cr,c4)) return 1;
				return rc;
			}
			if(ccout(c4)) return rc;
			if(g34Swap(cb,cl,cr,c4)) return 1;
			if(g34Swap(cr,cb,c4,cl)) return 1;
			return rc;
		}
		slr=i;
	}
	else if(ccn(cl)[i]==nu)	{slu=i;}
	else if(ccn(cl)[i]==no)	{slo=i;}
	else ibgfatal;
 }
 if(ccout(cl)) return 0;
 if(ccout(cr)) return 0;
 cxx = cnx(ns); x  = cxx[0];	 y  = cxx[1];	  z  = cxx[2];
 cxx = cnx(nu); xu = cxx[0]-x;   yu = cxx[1]-y;   zu = cxx[2]-z;
 cxx = cnx(nr); xr = cxx[0]-x;   yr = cxx[1]-y;   zr = cxx[2]-z;
 cxx = cnx(no); xo = cxx[0]-x;   yo = cxx[1]-y;   zo = cxx[2]-z;
 cxx = cnx(nl); xl = cxx[0]-x;   yl = cxx[1]-y;   zl = cxx[2]-z;
 vol0= ibggVAbs;
 det1 = xl*(yo*zr-zo*yr) + yl*(zo*xr-xo*zr) + zl*(xo*yr-yo*xr);
 det2 = xr*(yu*zl-zu*yl) + yr*(zu*xl-xu*zl) + zr*(xu*yl-yu*xl);
 if(det1 < ibggVRel*det2 + vol0) return rc;
 if(det2 < ibggVRel*det1 + vol0) return rc;
 for(i=0;i<ibgcSides(tr);i++){
	if(ccs(cr)[i]==cb)	{srb=i; ibgassert(ns==ccn(cr)[srb]);}
	else if(ccn(cr)[i]==nl)	{ibgassert(ccs(cr)[i]==cl); srl=i;}
	else if(ccn(cr)[i]==nu)	{sru=i;}
	else if(ccn(cr)[i]==no)	{sro=i;}
	else ibgfatal;
 }
 ccn(cb)[sbo] = ns;
 ccn(cl)[slu] = nl;
 co=ccs(cb)[sbu]; ccs(cb)[sbu]=cl; ccs(cl)[slb] = co;
 for(i=0;i<ibgcSides(cct(co));i++) if(ccs(co)[i]==cb) {ccs(co)[i]=cl;break;}
 co=ccs(cl)[slo]; ccs(cl)[slo]=cb; ccs(cb)[sbl] = co;
 for(i=0;i<ibgcSides(cct(co));i++) if(ccs(co)[i]==cl) {ccs(co)[i]=cb;break;}
 co=ccs(cr)[sru]; ccs(cl)[slr] = co;
 for(i=0;i<ibgcSides(cct(co));i++) if(ccs(co)[i]==cr) {ccs(co)[i]=cl;break;}
 co=ccs(cr)[sro]; ccs(cb)[sbr] = co;
 for(i=0;i<ibgcSides(cct(co));i++) if(ccs(co)[i]==cr) {ccs(co)[i]=cb;break;}
 if(ibgcRegion(cb)==ibgNotFound)if(st[cb]==0){st[cb] = nosegment0; nosegment0 = cb;}
 if(ibgcRegion(cl)==ibgNotFound)if(st[cl]==0){st[cl] = nosegment0; nosegment0 = cl;}
 g3cvol(cb);	g3cvol(cl);
 if(ccu(cr)<0) ibgridNFree(ibgg,-ccu(cr));
 g0CDel(cr);
 return 1;
}

int g3segment1(int c)
{ibgCellType t=cct(c);
 int rc,n,n0,e, *cn0=ccn(c);
beg:
 if(t!=ibgc3Tetrahedron) return 0;
 if(ibgcRegion(c)==ibgFound) return 1;
 g3segment1edge = 1;
 for(e=0;e<ibgcEdges(t);e++){
	if(ibgIncorrect != ibgiStatusOfEdge(ibGeo,cnd(cn0[ibgcEPoint1(t,e)])
				       ,cnd(cn0[ibgcEPoint2(t,e)]))) continue;
	if(g3EDelete(c,e)) return 1;
 }
 g3segment1edge = 0;
 n0 = -1;
 for(n=0;n<ibgcPoints(t);n++){
  	if(ibgpType(*cnd(cn0[n]))==ibgSLine){
		if(n0==-1) n0=n;
     		else goto next;
	}
  	if(ibgpType(*cnd(cn0[n]))==ibgSNode){
		n0 = n; goto vert;
	}
 }
 if(n0>=0){
oneline:
	rc = g3SDelete(c,n0);
	if(rc==1) return 1;
    	if(rc==2) return 2;
 }
vert:
next:
 return 0;
}
int g3segment2(int c)
{ibgCellType t=cct(c);
 if(t!=ibgc3Tetrahedron) return 0;
/* if(ibgcRegion(c)==ibgFound) return 1; */
 if(g3segment1(c)) return 1;
 if(g3CHalf(c)) return 1;
 return 0;
}


void ibgridMakeRegions(ibGrid *gg, ibGeometry g)
{int c,u;
 ibgg = *gg;
 ibGeo = g;
 ibgtCopy(&(ibgg.Topology),&(g->top));
 unc=0; mec=0; mbc=0; mrc=0;
 nosegmentb = -1;
 nosegment0 = -1;
 nosegment1 = -1;
 nosegment2 = -1;
 ibgcurrc   = 0;
 if((mst = ibgg.maxCell)<=0) return;
 st  = (int*) calloc(mst,sizeof(int));
 if(st==ibgNULL) ibgfatal;
/* at this stage undefined cells have a negative segment and a cell type
	ibgcNothing */
 for(c=1;c<=ibgg.lastCell;c++){
	if(ccund(c)) continue;
	if(ibgcRegion(c)==ibgNotFound){
	    st[c] = nosegmentb; nosegmentb = c; ibgcurrc++;
	}
 }
 if(ibgg.gridDimension==2){
        nosegment0 = nosegmentb;
	do{
		while(nosegment0>0){
		      	c = nosegment0; nosegment0 = st[c]; st[c] = 0;
			if(cnosegment(ccu(c))){
			        g2segment(c);
			        if(ccu(c)<0){ibgassert(st[c]==0);
			        	st[c] = nosegment2; nosegment2 = c;
			        }
			}
		}
		if(nosegment1>0){
			continue;
		}
		break;
	}while(1);
 }else if(ibgg.gridDimension==3){
	while(nosegmentb>0){
                c = nosegmentb; nosegmentb = st[c]; st[c] = 0;
	        if(ccund(c)) continue;
	        while(cnosegment(ccu(c))){
	            u = g3segment1(c);
	            if(st[c]!=0) break;
	            if(u==0){
	        	st[c] = nosegment0; nosegment0 = c;
	                break;
	            }
	        }
	}
	do{
		while(nosegment0>0){
		    if(ibgcurrc<0) break;
		    c = nosegment0; nosegment0 = st[c]; st[c] = 0;
		    if(ccund(c)) continue;
		    while(cnosegment(ccu(c))){
			ibgcurrc--;
			u = g3segment1(c);
			if(st[c]!=0) break;
			if(u==0){
			    st[c] = nosegment1; nosegment1 = c;
			    break;
			}
		    }
		}
		if(nosegment1>0){
	  	        c = nosegment1; nosegment1 = st[c]; st[c] = 0;
	   	        if(ccund(c)) continue;
		        while(cnosegment(ccu(c))){
		        	u = g3segment2(c);
		        	if(st[c]!=0) break;
		             	if(u==0){
		        		st[c] = nosegment2; nosegment2 = c;
		        		break;
		        	}
		        }
		        continue;
		}
		break;
	}while(1);
 }
 while(nosegment2>0){
      	c = nosegment2; nosegment2 = st[c]; st[c] = 0;
	if(ccund(c)) continue;
   	if(cnosegment(ccu(c))){
		ibgcSetRegion(c);
	}
 }
 free(st);
 ibgg.firstPoint   = 1;
 ibgg.firstCell   = 1;
 ibgg.firstRegion = 1;
 ibgg.lastRegion  = ibgg.lastCell;
 free(ibgg.CellData);
 ibgg.CellData = ibgNULL;
 *gg = ibgg;
}