ibgg0.c

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

 c0=0; cn0=ibgg.CPoint; cs0=ibgg.CSide; nc=0; na=1;
 for(n=1;n<=ibgg.lastPoint;n++){
	c0++;
	*(cn0++) = n;
	*(cs0++) = cndr(n,1);
 	if((*(cn0++) = *(cs0++) = cndr(n,0))==outside) ne=n;
	ibgg.CPointFirst[c0] = ibgg.CSideFirst[c0] = nc;
	nc += 2;
 	cct(c0) = ibgc1Edge;
 	ccu(c0) = 0;
 }
 ibgassert(cndn(na,0)==outside);
 ibgassert(cndp(ne,0)==outside);
 ibgg.lastPoint += 2; if(g0NRealloc()==0) return ibgError;
 xm=0.5*(cgmin(0)+cgmax(0)); dx=(cgmax(0)-cgmin(0));
 big=0.01*dx*dx/ibGeo->Delta;
 cnx(ibgg.lastPoint-1)[0] = xm-big;
 cnx(ibgg.lastPoint  )[0] = xm+big;
 cnu(ibgg.lastPoint-1)    = nothing;
 cnu(ibgg.lastPoint  )    = nothing;
 cnt(ibgg.lastPoint-1)    = ibgSRegion;
 cnt(ibgg.lastPoint  )    = ibgSRegion;
 c0++;
 ibgg.CPointFirst[c0] = ibgg.CSideFirst[c0] = nc;
 nc += 2;
 cct(c0) = ibgc1Edge;
 ccn(c0)[1] = na;
 ccs(c0)[1] = na;		ccs(na)[0] = c0;
 ccn(c0)[0] = ibgg.lastPoint-1;	ccs(c0)[0] = 0;
 ccn(ne)[1] = ibgg.lastPoint;	ccs(ne)[1] = 0;
 ibgassert(nc<ibgg.maxCPoint);
 ibgassert(nc<ibgg.maxCSide);
 ibgassert(c0<ibgg.maxCell);
 ibgg.freeCPoint = ibgg.freeCSide = nc;
 ibgg.lastCell  = c0;
 ibgassert(c0==ibgg.lastPoint-1);
 *ibgg0.gg = ibgg;
 *g0 = ibgg0;
 return ibgSuccess;
}

/* Lot auf Kante in Richtung rr in <rr,ro>-Fl"ache. */
static int g2elot(int nn, unsigned rr, unsigned ro)
{
 int no,nu,nor,nur,noru,nuro,norr,nurr,dd;
 no =cndr(nn,ro);		if(cnund(no))	return nothing;
 nu =cndr(nn,rother(ro));	if(cnund(nu))	return nothing;
 nor=cndr(no,rr);		if(cnund(nor))	return nothing;
 nur=cndr(nu,rr);		if(cnund(nur))	return nothing;
 noru=cndr(nor,rother(ro));
 if(cnnothing(noru)){
#ifdef GLEVEL
 	dd=rposit(rr); if(cnl(no)[dd]>=cnl(nor)[dd]) return nothing;
#endif
	norr=cndr(nor,rr);
	if(cnund(norr))		return nothing;
	noru=cndr(norr,rother(ro));
	if(cnnothing(noru))	return nothing;
	if(noru==nur)		return g2eref(noru,ro);
	if(noru==cndr(nur,rr)){
		cndr(nor,rother(ro)) = nur;
		cndr(nur,ro)	  = nor;
				return g2eref(nur,ro);
	}
	nurr=cndr(nur,rr);
	if(cnund(nurr))		return nothing;
	if(noru==cndr(nurr,ro)){
		cndr(nor,rother(ro)) = nur;
		cndr(nur,ro)	  = nor;
				return g2eref(nur,ro);
	}
	else 			return nothing;
 }
 if(noru==nur)			return g2eref(noru,ro);
 if(noru==cndr(nur,rr))		return g2eref(noru,ro);
 if(noru==cndr(nur,ro)){
 	cndr(nn,rr) = noru;
 	cndr(noru,rother(rr)) = nn;
				return noru;
 }
#ifdef GLEVEL
 dd=rposit(rr); if(cnl(nu)[dd]>=cnl(nur)[dd]) return nothing;
#endif
 nurr=cndr(nur,rr);
 if(cnund(nurr))		        return nothing;
 if(noru==cndr(nurr,ro)){
 	cndr(nn,rr) = noru;
	cndr(noru,rother(rr)) = nn;
				return noru;
 }
 else	 			return nothing;
}

static int g2eref(int nu, unsigned ro)
{int o,u,b,bu,d,dh,dr,eo,e1,rr,r1,ru;
 int no,nh,nur,nor,nr,noo,nuu;
 ibgFloat dx,dy;
 if(rispos(ro))	{d=ro;		ru=ro+cgdim;	no=cndr(nu,ro);}
 else		{d=ro-cgdim;	ru=ro;ro=d;	no=nu;nu=cndr(no,ru);}
 if((cnx(no)[d]-cnx(nu)[d])<=ibggDmin[d]) return nothing;

/* Schaffung der Voraussetzungen f"ur eine Unterteilung: */
 o=cror(d,0);
 if(cnund(cndr(nu,o))) g2elot(nu,o,d);
 if(cnund(cndr(no,o))) g2elot(no,o,d);
 u=rother(o);
 if(cnund(cndr(nu,u))) g2elot(nu,u,d);
 if(cnund(cndr(no,u))) g2elot(no,u,d);

/* Eigentliche Halbierung der Kante: */
 nh=g0NNew;
 cndr(nh,0) = cndr(nh,1) = cndr(nh,2) = cndr(nh,3) = nothing;
 cndr(nu,ro) = nh; cndr(nh,ru) = nu; cndr(no,ru) = nh; cndr(nh,ro) = no;
 if(ibgg0.nb!=ibgNULL) {cnb(nh) = 0;}

/* Definieren der Koordinaten und Level: */
 cnx(nh)[o] = cnx(nu)[o];
 cnx(nh)[2] = cnx(nu)[2];
 cnx(nh)[d] = 0.5*(cnx(nu)[d]+cnx(no)[d]);
#ifdef GLEVEL
 cnl(nh)[o] = cnl(nu)[o];
 cnl(nh)[2] = 0;
 cnl(nh)[d] = cnl(nu)[d]+1;
 if(cnl(nh)[d]<=cnl(no)[d]) cnl(nh)[d] = cnl(no)[d]+1;
#endif

/* Definieren der Daten auf dem neuen Knoten: */
 ibgiRegionOfPoint(ibGeo,cnd(nh),cnd(nu));

/* Hilfsstack der Reihenfolge der Knotenerzeugung: */
 ibgg0.no[nh]=nu;

/* Herstellen der Verbindungen zu den direkt sichtbaren Nachbarn: */
 for(o=0;o<MAXO2;o++){rr=g2ror[ro][o];
        nur=cndr(nu,rr);
        if(nur==outside)	cndr(nh,rr) = outside;
        else{
        	ibgassert(cndef(nur));
        	nr=cndr(nur,ro);
        	if(cnnothing(nr)) {
        		nur=cndr(nur,rr);	ibgassert(cndef(nur));
        		nr =cndr(nur,ro);	ibgassert(cndef(nr));
        	}
        	nor=cndr(no,rr); ibgassert(cndef(nor));
        	if(nr==nor)		 cndr(nh,rr) = -nr;
        	else if(nr==cndr(nor,rr))	 cndr(nh,rr) = -nr;
        	else {
        		ibgassert((cndr(nr,ro)==nor)||(cndr(nr,ro)==cndr(nor,rr)));
        		cndr(nr,rother(rr)) = nh; cndr(nh,rr) = nr;
        	}
        }
 }
 dx=cnx(no)[d]-cnx(nu)[d];
 for(o=0;o<MAXO2;o++){rr=g2ror[ro][o]; dr=rposit(rr);
        if(cnnothing(nr=cndr(nh,rr))){
        	nr=cndr(nu,rr);
        	if(rr>=cgdim)	dy=(cnx(nu)[dr]-cnx(nr)[dr])/dx;
        	else		dy=(cnx(nr)[dr]-cnx(nu)[dr])/dx;
        	if(dy<ibggOref) g2elot(nh,rr,ro);
        }else{
        	if(rr>=cgdim)	dy=(cnx(nu)[dr]-cnx(nr)[dr])/dx;
        	else		dy=(cnx(nr)[dr]-cnx(nu)[dr])/dx;
        	if(dy>ibggOinv){
        		if(cnund(cndr(nr=cndr(nu,rr),ro)))
        			g2elot(nr,ro,rr);
        		else if(cnund(cndr(nr=cndr(no,rr),ru)))
        			g2elot(nr,ru,rr);
        	}
        }
 }
 if(cndef(noo=cndr(no,ro))){
        if(cnx(noo)[d]-cnx(no)[d]>ibggPref*dx)
        			g2eref(no,ro);
 }
 if(cndef(nuu=cndr(nu,ru))){
        if(cnx(nu)[d]-cnx(nuu)[d]>ibggPref*dx)
        			g2eref(nuu,ro);
 }
 if(ibgg0.nb!=ibgNULL) {g0bjump(nh);}
 return nh;
}/* end of g2eref */


ibgFloat g2cgh=1.0;
static ibgFloat g2cgo,g2cgu,g2ctg,g2tan,g2hlf;

/* Die eigentliche Halbierungroutine eines W黵fels.  Halbiert nur wenn
alle Halbierungspunkte in dieser Richtung bereits vorhanden sind. */

int g3qref0(int qq, unsigned ro)
{
 int nu,no,nh,n0,n1,n2,n3,nr,nl,nur,nhr,nor,nurr,ou,oo,oh;
 int qn,o,w,eo,eu,fo,fu,go,gu;
 int rr,ru=rother(ro);
/* Test der Ausf"uhrbarkeit */
 if(cnund(qq)) return 0;
 nu = cqn(qq,(croq(ru,3)));  no = cqn(qq,(int)croq(ro,3));
 for(o=0;o<MAXO3;o++){
        ou = nu; nu = cqn(qq,croq(ru,o));
        oo = no; no = cqn(qq,croq(ro,o));
        if(no==cndr(nu,ro))	return 0;
        if((nor=cndr(oo,cror(ro,o)))==no) continue;
        if((nur=cndr(ou,cror(ru,o)))==nu) continue;
/* why the following is possible? Unknown, but has happened.
        if(nor==nothing) return 0;
        if(nur==nothing) return 0;
	*/
        if(nor==cndr(nur,ro))	return 0;
 }
/* Eigentliche Unterteilung des Quaders: */
 qn = g0QNew;
 eu = croq(ru,3); eo = croq(ro,3);
 nu = cqn(qq,eu); nh = cndr(nu,ro); ibgAssert(cndef(nh));
 no = cqn(qq,eo); ibgAssert(cndr(nh,ro) == no);
 fu = cqother(eo); fo = cqother(eu);
 for(o=0;o<MAXO3;o++){
        ou = nu; gu = fu; oh = nh; go = fo;
        eu = croq(ru,o);  eo = croq(ro,o);
        fu = cqother(eo); fo = cqother(eu);
        nu = cqn(qq,eu);  nh = cndr(nu,ro); ibgAssert(cndef(nh));
	no = cqn(qq,eo);  ibgAssert(cndr(nh,ro) == no);
	cqn(qn,eu) = nu;  cqn(qn,eo) = nh; cqn(qq,eu) = nh;
        cnq(nu,fo) = qn;  cnq(nh,fu) = qn; ibgAssert(cnq(nh,fo)==qq);
	rr = cror(ru,o);
	if((nhr=cndr(oh,rr))!=nh) {
	    ibgAssert(cndef(nhr));
            ibgAssert(cndr(nhr,rr)==nh);
	    cnq(nhr,fu)=qn; cnq(nhr,gu)=qn;
	}
	if((nur=cndr(ou,rr))!=nu) {cnq(nur,fo)=qn; cnq(nur,go)=qn;
	        if(o>1)continue;
	        if(cnund(nurr=cndr(nur,cror(ru,o+1)))) continue;
	        for(w=0;w<MAXO3;w++){
	        	if(cnq(nurr,croq(ro,w))==qq)
	        	   cnq(nurr,croq(ro,w))= qn;
	        }
	}
 }
/* evtl. Eintragung einer Verbindung auf der Halbierungsfl"ache: */
 n0 = cqn(qq,croq(ru,0));
 n1 = cqn(qq,croq(ru,1));
 n2 = cqn(qq,croq(ru,2));
 n3 = cqn(qq,croq(ru,3));
 if(n1==(nr=cndr(n0,cror(ru,1)))){
	if(n2==(nr=cndr(n1,cror(ru,2)))) return 1;
	if(n0==(nl=cndr(n3,cror(ru,0)))) return 1;
	cndr(nr,cror(ru,3)) = nl;
	cndr(nl,cror(ru,1)) = nr;	 return 1;}
 else	if(n3==(nl=cndr(n2,cror(ru,3)))){
	if(n2==(nr=cndr(n1,cror(ru,2)))) return 1;
	if(n0==(nl=cndr(n3,cror(ru,0)))) return 1;
	cndr(nr,cror(ru,3)) = nl;
	cndr(nl,cror(ru,1)) = nr;	 return 1;}
 else  {cndr(nr,cror(ru,2)) = nl;
	cndr(nl,cror(ru,0)) = nr;	 return 1;}
}

static int g3qref2(int qq, unsigned ro)
{
 int nu,no,nn,np,nur,nor,ou,oo;
 int o,ru=rother(ro);
 if(cnund(qq)) return 0;
 nn = cqn(qq,ibgg0.firstq);		np = cqn(qq,ibgg0.lastq);
 nu = cqn(qq,(croq(ru,3)));	no = cqn(qq,(int)croq(ro,3));
 for(o=0;o<MAXO3;o++){
 	ou = nu; nu = cqn(qq,croq(ru,o));
 	oo = no; no = cqn(qq,croq(ro,o));
 	if(no==cndr(nu,ro)){
 		g3eref(nu,ro);
 	 	if(nn!=cqn(qq,ibgg0.firstq)) return 1;
 		if(np!=cqn(qq,ibgg0.lastq)) return 1;
 	}
 	if((nor=cndr(oo,cror(ro,o)))==no) continue;
 	if((nur=cndr(ou,cror(ru,o)))==nu) continue;
        if(nor==nothing) return 0;
	if(nur==nothing) return 0;
  	if(nor==cndr(nur,ro)){
 		g3eref(nur,ro);
 	 	if(nn!=cqn(qq,ibgg0.firstq)) return 1;
 		if(np!=cqn(qq,ibgg0.lastq)) return 1;
 	}
 }
 return g3qref0(qq,ro);
}

static int g3qref1(int qq, unsigned ro)
{
 int nu,no,nur,nor,ou,oo;
 int o,ru=rother(ro);
 if(cnund(qq)) return 0;
 nu = cqn(qq,(croq(ru,3)));  no = cqn(qq,(int)croq(ro,3));
 for(o=0;o<MAXO3;o++){
 	ou = nu; nu = cqn(qq,croq(ru,o));
 	oo = no; no = cqn(qq,croq(ro,o));
        if(no!=cndr(nu,ro))	return g3qref2(qq,ro);
        if((nor=cndr(oo,cror(ro,o)))==no) continue;
        if((nur=cndr(ou,cror(ru,o)))==nu) continue;
        if(nor!=cndr(nur,ro))	return g3qref2(qq,ro);
 }
 return 0;
}

/* Die zentrale Kantenverfeinerungsroutine */

static int g3eref(int nu, unsigned ro)
{int qq,qo,o,b,d,dh,dr,eo,e1,rr,r1,ru;
 int no,nh,nur,nurr,nor,noo,nuu,nr,lo;
 ibgFloat dx,dy;
 if(rispos(ro))	{d=ro;		ru=ro+cgdim;	no=cndr(nu,ro);}
 else		{d=ro-cgdim;	ru=ro;ro=d;	no=nu;nu=cndr(no,ru);}
 if(cnund(no)) return nothing;
/* Schaffung der Voraussetzungen f"ur eine Unterteilung: */
 for(o=0;o<MAXO3;o++){
        if(cnund(qq= cnq(nu,lo=croq(ro,o))))continue;
        if     (qq==cnq(no,lo))
        	{if(!g3qref2(qq,ro)) return nothing;}
        else if(qq==cnq(nu,croq(ru,o)))
        	{if(!g3qref2(qq,ro)) return nothing;}
 }
/* Hartes externes Verbot einer weiteren Unterteilung: */
 if((cnx(no)[d]-cnx(nu)[d])<=ibggDmin[d]) return nothing;
/* Eigentliche Halbierung der Kante: */
 nh=g0NNew; if(nh==0) {return nothing;}
 cndr(nu,ro) = nh; cndr(nh,ru) = nu; cndr(no,ru) = nh; cndr(nh,ro) = no;
 if(ibgg0.nb!=ibgNULL) {cnb(nh) = 0;}
/* Definieren der Koordinaten und Level: */
 for(dh=0;dh<ibgDIM;dh++){
        if(d!=dh){
        	cnx(nh)[dh] = cnx(nu)[dh];
                cndp(nh,dh) = cndn(nh,dh) = nothing;
#ifdef GLEVEL
         	cnl(nh)[dh] = cnl(nu)[dh];
#endif
        }else{	cnx(nh)[dh] = 0.5*(cnx(nu)[dh]+cnx(no)[dh]);
#ifdef GLEVEL
        	cnl(nh)[dh] = cnl(nu)[dh]+1;
        	if(cnl(nh)[dh]<=cnl(no)[dh]) cnl(nh)[dh] = cnl(no)[dh]+1;
#endif
        }
 }
/* Definieren der Daten auf dem neuen Knoten: */
 ibgiRegionOfPoint(ibGeo,cnd(nh),cnd(nu));

/* Hilfsstack der Reihenfolge der Knotenerzeugung: */
 ibgg0.no[nh]=nu;
/* Definieren der nq-Zeiger: */
 for(o=0;o<MAXO3;o++){
        eo = croq(ro,o);
        cnq(nh,eo) = (cnq(nh,croq(ru,o)) = cnq(nu,eo));
 }
/* Herstellen der Verbindungen zu den direkt sichtbaren Nachbarn: */
 for(o=0;o<MAXO3;o++){
        rr = cror(ro,o);
        if(cnund(nur=cndr(nu,rr)))		{cndr(nh,rr)=  nur;continue;}
        if(cnund(nr=cndr(nur,ro))){
        	if(cnq(nh ,croq(ro,o)) !=
        	   cnq(nur,croq(ro,o)))		{cndr(nh,rr)= -nur;continue;}
        	if(cnund(nurr=cndr(nur,rr)))	{cndr(nh,rr)= -nur;continue;}
        	if(cnund(nr  =cndr(nurr,ro))) 	{cndr(nh,rr)= -nur;continue;}
        }
        if(cnund(nor=cndr(no,rr)))		{cndr(nh,rr)=  nor;continue;}
        if(nr==nor)		 		{cndr(nh,rr)= -nor;continue;}
        if(nr==cndr(nor,rr))			{cndr(nh,rr)= -nor;continue;}
        if((cndr(nr,ro)!=nor)&&(cndr(nr,ro)!=cndr(nor,rr)))
        					{cndr(nh,rr)= -nor;continue;}
        eo=croq(ro,o);	e1=croq(ro,o+1);
        if(cnq(nu,eo)==cnq(nu,e1))		{cndr(nh,rr)= -nor;continue;}
        cndr(nr,rother(rr)) = nh; cndr(nh,rr) = nr;
 }
/* Ende der eigentlichen Verfeinerung. Die direkten Nachbarn sind definiert */
/* Kontrolle m"oglicher Quaderhalbierungen: */
 qq=cnq(nh,croq(ro,3)); /* 3 == MAXO-1 */
 for(o=0;o<MAXO3;o++){qo=qq;
        if(qo==(qq=cnq(nh,croq(ro,o)))) continue;
        if(g3qref0(qq,ro))	qq=cnq(nh,croq(ro,o));
 }
 if(ibgg0.nb!=ibgNULL) {g0bjump(nh);}
 return nh;
}/* end of g3eref */

static int g3elot(int nn, unsigned rr, unsigned ro)
{
 int no,nu,nor,nur,noru,nuro,norr,nurr,dd;
 no =cndr(nn,ro);		if(cnund(no))	return nothing;
 nu =cndr(nn,rother(ro));	if(cnund(nu))	return nothing;
 nor=cndr(no,rr);		if(cnund(nor))	return nothing;
 nur=cndr(nu,rr);		if(cnund(nur))	return nothing;
 noru=cndr(nor,rother(ro));
 if(cnnothing(noru)){
#ifdef GLEVEL
 	dd=posit(rr); if(cnl(no)[dd]>=cnl(nor)[dd]) return nothing;
#endif
	norr=cndr(nor,rr);
	if(cnund(norr))		return nothing;
	noru=cndr(norr,rother(ro));
	if(cnnothing(noru))	return nothing;
	if(noru==nur)		return g3eref(noru,ro);
	if(noru==cndr(nur,rr)){
/*		cndr(nor,rother(ro)) = nur;
		cndr(nur,ro)	  = nor;
				return g3eref(nur,ro);
*/				return nothing;
	}
	nurr=cndr(nur,rr);
	if(cnund(nurr))		return nothing;
	if(noru==cndr(nurr,ro)){
/*		cndr(nor,rother(ro)) = nur;
		cndr(nur,ro)	  = nor;
				return g3eref(nur,ro);
*/				return nothing;
	}
	else 			return nothing;
 }
 if(noru==nur)			return g3eref(noru,ro);
 if(noru==cndr(nur,rr))		return g3eref(noru,ro);
 if(noru==cndr(nur,ro)){
 	cndr(nn,rr) = noru;
 	cndr(noru,rother(rr)) = nn;
				return noru;
 }
#ifdef GLEVEL
 dd=posit(rr); if(cnl(nu)[dd]>=cnl(nur)[dd]) return nothing;
#endif
 nurr=cndr(nur,rr);
 if(cnund(nurr))			return nothing;