ibgg0.c

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

 free(ibgg0.bstack);
 return ibgSuccess;
}

static int g0nfirst,g0nstb,g0nstl;
static int g0bfirst,g0bstb,g0bstl;

static int g0bshift(int nn, int firstcall)
{int m,bm,nu,no,nh,nr,lv,lv0,r0,d,u,r,o,bb,bo,bu,ne,be,rc;
 ibgFloat dx,dx0,dxmin;
 int lb[3],rb[3],lb0;
 if(cnb(nn)!=0) return 0;
beg:
/* Feststellung der Vorzugsrichtung r0 zum Verschieben des Knotens */
 m=cnu(nn); lv0=0; lb0=0; dxmin=big1;
 for(d=0;d<cgdim;d++){nu=cndn(nn,d);
notest:	if(cndef(no=cndp(nn,d))){
	        dx=cnx(no)[d]-cnx(nn)[d]; if(dx<dxmin) dxmin=dx;
	        if(cnb(no)==0){if(cnu(no)==m) goto nutest; lv=3;}
	        else if(cbt(bo=cnb(no))==ibgSFace) {
	        	lb0=1;
	        	if(ibgIncorrect!=ibgiStatusOfEdge(ibGeo,cnd(nn),cbd(bo))){
	        		if(cbx(bo)[d]<cnx(no)[d]) goto nutest;
	        	}
	        	if(m==cnu(no)) goto nutest;
	        	lv=2;
	        }else{	lb0=1; goto nutest;
	        }
	        if(nu==outside) lv=1;
	        if(lv>lv0)	{r0=d; dx0=dx; lv0=lv; ne=no;}
	        else if(lv==lv0){if(dx<dx0) {r0=d; dx0=dx; ne=no;}}
	}
nutest:	if(cndef(nu)){
	        dx=cnx(nn)[d]-cnx(nu)[d]; if(dx<dxmin) dxmin=dx;
	        if(cnb(nu)==0){if(cnu(nu)==m) continue; lv=3;}
	        else if(cbt(bu=cnb(nu))==ibgSFace) {
	        	lb0=1;
	        	if(ibgIncorrect !=
	        		ibgiStatusOfEdge(ibGeo,cnd(nn),cbd(bu))){
	        			if(cbx(cnb(nu))[d]>cnx(nu)[d]) continue;
	        	}
	        	if(m==cnu(nu)) continue;
	        	lv=2;
	        }else{	lb0=1; continue;
	        }
	        if(no==outside) lv=1;
	        if(lv>lv0)	{r0=cgdim+d; dx0=dx; lv0=lv; ne=nu;}
	        else if(lv==lv0){if(dx<dx0) {r0=cgdim+d; dx0=dx; ne=nu;}}
	}
 }
/* lv0=1 -> Verschiebung vom Au"senrand weg; lv0=2 -> Verschiebung zu face; */
 if(lv0>0){
	if(r0>=cgdim)	{d=r0-cgdim; r=r0; no=nn; nu=cndn(nn,d);}
	else    	{r=r0+cgdim; d=r0; nu=nn; no=cndp(nn,d);}
	if(lv0==1){
	        if(firstcall)	goto nextbsh;
	        if(cnb(ne)==0){
	        	if(cndr(no,d)==outside){
	        	    if(cndr(nu,r)==outside){
	        		if(cndef(ibgg0.eref(nu,d))) goto nextbsh;
	         		 return 0;
	        	    }
	        	}
	        }else{
	        	if(cndef(ibgg0.eref(nu,d))) goto nextbsh;
	        	 return 0;
	        }
	}
	if(dx0>ibggBLWait*dxmin && firstcall)	goto nextbsh;
	if(dx0>ibggBLHalf*dxmin){
	        if(firstcall)	goto nextbsh;
    	        if(cndef(ibgg0.eref(nu,d))) goto nextbsh;
	}
	bb = g0BNew; if(bb==0) return 0;
   	rc = ibgiFaceWithEdge(ibGeo,cbd(bb),cnd(nn),cnd(ne));
	ibgassert(rc==ibgFaceFound);
	bm = cbu(bb);
	if(cnb(ne)==0){
	        if(cndr(no,d)==outside){
/* Au"senrand no mu"s Au"senrand bleiben !!! */
	            ibgassert(cndr(nu,r)!=outside);
	            if(cbx(bb)[d]>=cnx(no)[d]-2*ibGeo->Delta){
	        	ibgassert(cnb(no)==0);
	        	cnb(no)=bb; cbn(bb)=no;
	             	cbx(bb)[d]=cnx(no)[d];
	        	if(no==nn) rc=0; else rc=nn;
	            }else{
	        	ibgassert(cnb(nu)==0);
	        	cnb(nu)=bb; cbn(bb)=nu;
	        	if(nu==nn) rc=0; else rc=nn;
	            }
	        }else if(cndr(nu,r)==outside){
/* Au"senrand nu mu"s Au"senrand bleiben !!! */
	            if(cbx(bb)[d]<=cnx(nu)[d]+2*ibGeo->Delta){
	        	ibgassert(cnb(nu)==0);
	        	cnb(nu)=bb; cbn(bb)=nu;
	        	cbx(bb)[d]=cnx(nu)[d];
	        	if(nu==nn) rc=0; else rc=nn;
	            }else{
	        	ibgassert(cnb(no)==0);
	        	cnb(no)=bb; cbn(bb)=no;
	        	if(no==nn) rc=0; else rc=nn;
	            }
	        }else if(cbx(bb)[d]<0.5*(cnx(nu)[d]+cnx(no)[d])){
	        	ibgassert(cnb(nu)==0);
	         	cnb(nu)=bb; cbn(bb)=nu;
	        	if(nu==nn) rc=0; else rc=nn;
	        }else{
	        	ibgassert(cnb(no)==0);
	        	cnb(no)=bb; cbn(bb)=no;
	        	if(no==nn) rc=0; else rc=nn;
	        }
 	}else{	/* the other point ne is already boundary point */
	        be=cnb(ne);
	        ibgassert(lv0>1);
	        rc=0;
	        if((cbx(bb)[d]-cnx(nn)[d])/(cbx(be)[d]-cnx(nn)[d])
	        			>=ibggBV){
	         	g0BDel(bb);
	        }else if(cbu(bb)==cbu(be)){
	        	if((cbx(bb)[d]-cnx(nn)[d])*(cbx(be)[d]-cnx(ne)[d])
	        			< 0.0){
/* that means, be is not compatible with cnu(nn)!!! */
	        		ibgmessage(ibgMWIncFFace);
	        		ibgmsg(" %d not face of %d ",cbu(bb),cnu(nn));
	        	 	g0BDel(bb);
	        	}else{
	        		ibgassert(cnb(nn)==0);
	        	 	cnb(nn)=bb; cbn(bb)=nn;
	        	}
	        }else{
	        	ibgassert(cnb(nn)==0);
	         	cnb(nn)=bb; cbn(bb)=nn;
	        }
	        if(r0<cgdim){
#ifdef GDEBUG
	              	if(cbx(bb)[d]>cbx(be)[d]){
	        		ibgassert(cbx(be)[d]>=cnx(ne)[d]-ibGeo->Delta);
	        		ibgassert(cbx(bb)[d]<=cnx(ne)[d]+ibGeo->Delta);
	        	}
#endif
	        	if(cbx(bb)[d]>cbx(be)[d]-ibGeo->Delta)
				cbx(bb)[d]=cbx(be)[d]-ibGeo->Delta;
		}else{
#ifdef GDEBUG
		      	if(cbx(bb)[d]<cbx(be)[d]){
				ibgassert(cbx(be)[d]<=cnx(ne)[d]+ibGeo->Delta);
				ibgassert(cbx(bb)[d]>=cnx(ne)[d]-ibGeo->Delta);
			}
#endif
			if(cbx(bb)[d]<cbx(be)[d]+ibGeo->Delta)
				cbx(bb)[d]=cbx(be)[d]+ibGeo->Delta;
		}
	}
	ibgassert(cbx(bb)[d]>=cnx(nu)[d]-ibGeo->Delta);
	ibgassert(cbx(bb)[d]<=cnx(no)[d]+ibGeo->Delta);
	if(cbx(bb)[d]<cnx(nu)[d]+ibGeo->Delta) cbx(bb)[d]=cnx(nu)[d]+ibGeo->Delta;
	if(cbx(bb)[d]>cnx(no)[d]-ibGeo->Delta) cbx(bb)[d]=cnx(no)[d]-ibGeo->Delta;
	for(o=0;o<ibgDIM;o++) if(o!=d){
		ibgassert(cbx(bb)[o]<cnx(nu)[o]+ibGeo->Delta);
		ibgassert(cbx(bb)[o]>cnx(nu)[o]-ibGeo->Delta);
		cbx(bb)[o]=cnx(nu)[o]; /* Rundungsfehler-Korrektur */
	}
	for(o=0;o<ibgg0.maxo;o++){
		r=cror(d,o);
		if(cnund(nr=cndr(nn,r)))	continue;
		if(cnb(nr)!=0)		continue;
	 	if(cnu(nr)==m)		continue;
/* including nr into the stack */
		if(nr>g0nfirst)		continue;
		if(ibgg0.nstack[nr]!=0)	continue;
		if(firstcall){
			if(nr>nn)	continue;
		}
		if(g0nstb<=0)		   g0nstb  = nr;
		else		ibgg0.nstack[g0nstl] = nr;
		ibgg0.nstack[nr] = -1;       g0nstl  = nr;
	}
	if(rc==nn) goto nextbsh;
 }else if(lb0){
	if(firstcall)	goto nextbsh;
	for(d=0;d<cgdim;d++) lb[d] = 0;
	for(r=0;r<ibgg0.maxr;r++)	if(cndef(no=cndr(nn,r))){
		if(cnb(no)==0)		continue;
		bm=cnb(no);
	   	for(d=0;d<cgdim;d++) if(d!=rposit(r)){
	 		if(cnx(no)[d]!=cbx(bm)[d]) {lb[d] = 1; rb[d] = r;}
		}
	}
	for(r=0;r<ibgg0.maxr;r++)	if(cnnothing(cndr(nn,r))){
		if(lb[rposit(r)]==0)	continue;
	  	if(cnund(ibgg0.elo2(nn,r,rb[rposit(r)]))){
			ibgfatal;
		}
	}
 }
 return 0;
nextbsh:
 if(nn <= g0nfirst && ibgg0.nstack[nn]==0){
	if(g0nstb<=0)	g0nstb = nn;
	else		ibgg0.nstack[g0nstl] = nn;
	ibgg0.nstack[nn] = -1;	   g0nstl  = nn;
 }
 return 0;
error:
 ibgfatal; return 0;
}



static void g0bjump(int n0)
{int r,r0,d,d0,nr,b,nr0,no,b0,i,rc;
 ibgFloat dx,dd,dx0,dd0=big1;
 int m=cnu(n0);
 for(r=0;r<ibgg0.maxr;r++){
	if(cnund(nr=cndr(n0,r)))	continue;
	if(cnb(nr)==0){
	        if(m==cnu(nr))		continue;
	        if(nr<=g0nfirst && ibgg0.nstack[nr]==0){
	        	if(g0nstb<=0)	g0nstb = nr;
	        	else		ibgg0.nstack[g0nstl] = nr;
	        	ibgg0.nstack[nr] = -1;	   g0nstl  = nr;
	        }
	        continue;
	}
	b=cnb(nr); d=rposit(r);
	if(r==d){
	        if(cbx(b)[d]>0.5*(cnx(nr)[d]+cnx(n0)[d])-ibGeo->Delta) continue;
	        dx = cbx(b)[d] - cnx(n0)[d];
	        dd = cnx(nr)[d] - cnx(n0)[d];
	}else{	if(cbx(b)[d]<0.5*(cnx(nr)[d]+cnx(n0)[d])+ibGeo->Delta) continue;
	        dx = cnx(n0)[d] - cbx(b)[d];
	        dd = cnx(n0)[d] - cnx(nr)[d];
	}
	if(dx<0)	{nr0=nr; b0=b; dd0=dd; r0=r; dx0 = dx; break;}
	if(dd<dd0)	{nr0=nr; b0=b; dd0=dd; r0=r; dx0 = dx; }
 }
 if(dd0==big1){
	if(n0<=g0nfirst) g0bshift(n0,1);
moverflow:      return;
 }
 d0 = rposit(r0);
#ifdef GDEBUG
 if(cbt(cnb(nr0))==ibgSFace){
 	for(d=0;d<cgdim;d++) if(d!=d0) ibgassert(cnx(n0)[d]==cbx(b0)[d]);
 }
#endif
/* special regularizations if dx0 is near 0: */
 if(dx0<0){
	if(cbt(cnb(nr0))==ibgSFace)	if(cnu(n0)!=cnu(nr0)){
		 				/* error in b0-computation!!!*/
		rc = ibgiFaceWithEdge(ibGeo,cbd(b0),cnd(n0),cnd(nr0));
		ibgassert(rc==ibgFaceFound);
		for(d=0;d<cgdim;d++)if(d!=d0)(cbx(b0)[d]=cnx(n0)[d]);
		if(rispos(r0)){
	      		if(cbx(b0)[d0]<cnx(n0)[d0]+ibGeo->Delta)
			   cbx(b0)[d0]=cnx(n0)[d0]+ibGeo->Delta;
			if(cbx(b0)[d0]>cnx(nr0)[d0]-ibGeo->Delta)
			   cbx(b0)[d0]=cnx(nr0)[d0]-ibGeo->Delta;
		}else{
			if(cbx(b0)[d0]>cnx(n0)[d0]-ibGeo->Delta)
			   cbx(b0)[d0]=cnx(n0)[d0]-ibGeo->Delta;
			if(cbx(b0)[d0]<cnx(nr0)[d0]+ibGeo->Delta)
			   cbx(b0)[d0]=cnx(nr0)[d0]+ibGeo->Delta;
		}
		goto shift;
      	}
 	if(rispos(r0)){
		if(cbx(b0)[d0]>cnx(n0)[d0]-ibGeo->Delta)
		   cbx(b0)[d0]=cnx(n0)[d0]-ibGeo->Delta;
	}else{
      		if(cbx(b0)[d0]<cnx(n0)[d0]+ibGeo->Delta)
		   cbx(b0)[d0]=cnx(n0)[d0]+ibGeo->Delta;
	}
 }else{
	if(cbt(cnb(nr0))==ibgSFace)	if(cnu(n0)==cnu(nr0)){
				 		/* error in b0-computation!!!*/
		no = cndr(n0,rother(r0));
		ibgassert(cndef(no));
		ibgassert(cnu(no)!=cnu(n0));
		rc = ibgiFaceWithEdge(ibGeo,cbd(b0),cnd(n0),cnd(no));
		ibgassert(rc==ibgFaceFound);
		for(d=0;d<cgdim;d++)if(d!=d0)(cbx(b0)[d]=cnx(n0)[d]);
		if(rispos(r0)){
			if(cbx(b0)[d0]>cnx(n0)[d0]-ibGeo->Delta)
			   cbx(b0)[d0]=cnx(n0)[d0]-ibGeo->Delta;
			if(cbx(b0)[d0]<cnx(no)[d0]+ibGeo->Delta)
			   cbx(b0)[d0]=cnx(no)[d0]+ibGeo->Delta;
		}else{
	      		if(cbx(b0)[d0]<cnx(n0)[d0]+ibGeo->Delta)
			   cbx(b0)[d0]=cnx(n0)[d0]+ibGeo->Delta;
			if(cbx(b0)[d0]>cnx(no)[d0]-ibGeo->Delta)
			   cbx(b0)[d0]=cnx(no)[d0]-ibGeo->Delta;
		}
		goto shift;
	}
	if(rispos(r0)){
      		if(cbx(b0)[d0]<cnx(n0)[d0]+ibGeo->Delta)
		   cbx(b0)[d0]=cnx(n0)[d0]+ibGeo->Delta;
	}else{
		if(cbx(b0)[d0]>cnx(n0)[d0]-ibGeo->Delta)
		   cbx(b0)[d0]=cnx(n0)[d0]-ibGeo->Delta;
	}
 }
shift:
 ibgassert(cnb(nr0)==b0);
 cnb(nr0)= 0;
 ibgassert(cnb(n0)==0);
 cnb(n0) = b0;		cbn(b0)  = n0;
 if(b0<=g0bfirst && ibgg0.bstack[b0]==0){
	if(g0bstb<=0)	g0bstb = b0;
	else{
		ibgassert(ibgg0.bstack[g0bstl] == -1);
		ibgg0.bstack[g0bstl] = b0;
	}
	ibgg0.bstack[b0] = -1;       g0bstl  = b0;
 }
 g0bjump(nr0);
}

/* regularising around the Face-point nn: */

/* potential bug: really we have to test the position of b1 on the edge. */
#define g0boundtest(d,n,bm)\
 if(!d){if(cnb(n)==0);\
	else if(cbt(cnb(n))==ibgSFace){if(cbu(cnb(n))==bm) d=1; else d= -1;}\
	else d=1;}\
 else if(cbt(cnb(n))>=ibgSLine){d=1;}

static int g0breg(int bb)
{int nn,bo,bu,b3,bh,br,dd,d,dr,di,ro,ru,rr,r1,r2,o,o1,qq,qo,nm,bm,om,ur;
 ibgSegmentType tnn;
 int upper0,upper,lower,rc,hoch;
 int no,nu,noo,nuu,noh,nuh,nor,nur,nr,n1,n2,n3,n4,nv,nv0,nv1,nv2,nv3,nvo,nvr;
 ibgFloat dx,dy,dl,del,du,dx1,dx2;
 nn = cbn(bb);
 switch(cbt(bb)){
 	case ibgSFace:
/* find shift direction dd, point no (oben), shift length dx, edge length dl */
 ibgrefface(cbd(bb),&ibgNorm,&ibgtang);
 for(dd=0;dd<cgdim;dd++){
	if((dx=cbx(bb)[dd]-cnx(nn)[dd])==0) continue;
	if(dx>0){ro=dd; ru=cgdim+dd;
		no=cndr(nn,ro); ibgAssert(cndef(no));
		dl=cnx(no)[dd]-cnx(nn)[dd];
		break;
	}else	{ro=cgdim+dd; ru=dd; dx = -dx;
                no=cndr(nn,ro); ibgAssert(cndef(no));
		dl=cnx(nn)[dd]-cnx(no)[dd];
		break;
	}
 }
 ibgassert(cndef(no));
 ibgassert(dx>0);
 if(dl>ibgNorm) {if(cndef(ibgg0.eref(nn,ro))) goto nextbreg;}
/* test correctness of segments */
 bm = cbu(bb); nm = cnu(nn);
 if(ibgiStatusOfEdge(ibGeo,cnd(nn),cbd(bb))==ibgIncorrect)
  	{if(cndef(ibgg0.eref(nn,ro))) goto nextbreg;}
 if(cndef(nu=cndr(nn,ru))){
   if(cnb(nu)!=0){
	if(cnu(nu)!=nm){
		if(ibgiStatusOfEdge(ibGeo,cnd(nn),cbd(cnb(nu)))==ibgIncorrect){
			if(cndef(ibgg0.eref(nn,ru))) goto nextbreg;
		}
	}
   }
 }
/* here we start to detect if we have found another intersection with
   the boundary bm going around the rectangle starting in lower/upper
   direction: lower/upper == 0 - unknown; 1 - found; -1 - not found; */
 if(cnb(no)==0){
	if(ibgiStatusOfEdge(ibGeo,cnd(no),cbd(bb))==ibgIncorrect){
		if(cndef(ibgg0.eref(nn,ro))) goto nextbreg;
		ibgmessage(ibgMWGridCoarse);
	 	upper0 = -1;
	}else{
	 	upper0 = 0;
	}
 }else if(cbt(cnb(no))==ibgSFace){
	if(rispos(ro)){
	  if(cbx(cnb(no))[dd]>cnx(no)[dd]){
/* potential error: refinement too often */
		if(ibgiStatusOfEdge(ibGeo,cnd(no),cbd(bb))==ibgIncorrect){
			if(cndef(ibgg0.eref(nn,ro))) goto nextbreg;
			ibgmessage(ibgMWGridCoarse);
		 	upper0 = -1;
		}else{
		 	upper0 = 0;
		}
	  }else upper0 = 1;
	}else if(cbx(cnb(no))[dd]<cnx(no)[dd]){
		if(ibgiStatusOfEdge(ibGeo,cnd(no),cbd(bb))==ibgIncorrect){
			if(cndef(ibgg0.eref(nn,ro))) goto nextbreg;
		     	ibgmessage(ibgMWGridCoarse);
			upper0 = -1;
		}else{
			upper0 = 0;
		}
	}else	upper0 = 1;
	if(upper0==1){
/* Anderer Schnitt des Rechtecks mit Rand gefunden. Ist es der Richtige? */
		if(bm==cbu(cnb(no))) upper0 = 1;
	 	else if(ibgIncorrect!=
			ibgiStatusOfEdge(ibGeo,cbd(bb),cbd(cnb(no)))){
/* Verschiedene Raender, aber gleiches Nachbargebiet vorhanden */
		     	upper0 = -1;
		}else{
			if(cndef(ibgg0.eref(nn,ro))) goto nextbreg;
		     	ibgmessage(ibgMWGridCoarse);
		     	upper0 = -1;
		}
	}
 }else{	upper0 = 1;
 }
/* look into the sides of the edge: */
 if(cgdim==3)		qq=cnq(nn,croq(ro,0));
 else if(cgdim<2)	return 0;
 for(o=1;o<=ibgg0.maxo;o++){
	if(cgdim==3){	qo=qq; if(qo==(qq=cnq(nn,croq(ro,o)))) continue;
			nur=cndr(nn,rr=cror(ro,o-1));
	}else if((nur=cndr(nn,rr=cror(ro,o-1)))==outside) continue;
/* we have found a side of a cuboid. Let's go around it. We find the
   corners nuu,noo,nur,nor. If there are refined region points on the edges
   or other dangerous situations (obtuse angles)