⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 ibgg0.c

📁 有限元学习研究用源代码(老外的),供科研人员参考
💻 C
📖 第 1 页 / 共 5 页
字号:
 if(noru==cndr(nurr,ro)){
 	cndr(nn,rr) = noru;
	cndr(noru,rother(rr)) = nn;
			   	return noru;
 }
 else	 			return nothing;
}

static int g3elo2(int nn, unsigned rr, unsigned ro)
{int rc,q0;
 if(cndef(rc=g3elot(nn,rr,ro)))	return rc;
 q0=cnq(nn,ibgg0.rrq[rr][ro]);
 ibgassert(q0==cnq(nn,ibgg0.rrq[ro][rr]));
 g3qref2(q0,ibgg0.rrr[ro][rr]);
 if(cndef(rc=cndr(nn,rr)))	return rc;
 if(cndef(rc=g3elot(nn,rr,ro)))	return rc;
 if(cndef(rc=g3elot(nn,rr,ibgg0.rrr[ro][rr])))	return rc;
 ibgfatal;
}

/* parameter-controlled refinement of a cuboid */
static int g3qref(int qq)
{int nn,np,rc=0;
 ibgFloat dx,dy,dz;
beg:
 nn = cqn(qq,ibgg0.firstq); np = cqn(qq,ibgg0.lastq);
 dx = cnx(np)[0]-cnx(nn)[0];
 dy = cnx(np)[1]-cnx(nn)[1];
 dz = cnx(np)[2]-cnx(nn)[2];
/* Omax - maximal allowed anisotropy of a cuboid */
 if(dx>ibggOmax*dy)	{if((rc=g3qref2(qq,0))>0) goto beg;}
 if(dx>ibggOmax*dz)	{if((rc=g3qref2(qq,0))>0) goto beg;}
 if(dy>ibggOmax*dz)	{if((rc=g3qref2(qq,1))>0) goto beg;}
 if(dy>ibggOmax*dx)	{if((rc=g3qref2(qq,1))>0) goto beg;}
 if(dz>ibggOmax*dx)	{if((rc=g3qref2(qq,2))>0) goto beg;}
 if(dz>ibggOmax*dy)	{if((rc=g3qref2(qq,2))>0) goto beg;}
 return rc;
}

static void g3nref(int nn)
{int ref=1,i,q,d,u,o,rv,dv,rq,dq,e,no,noo,nu,nuu,nr,nur,nor,qq,qo;
 ibgFloat dx,dy,dz;
 ibgFloat xx[ibgDIM],lmax;
 ibgrefreg(cnd(nn),&lmax);
/* isotropic refinement (maxlength) in positive directions: */
 xx[0]=cnx(nn)[0]+lmax;
 xx[1]=cnx(nn)[1]+lmax;
 xx[2]=cnx(nn)[2]+lmax;
 while(ref){ref=0;
	for(d=0;d<cgdim;d++){
		if(cndef(no=cndr(nn,d))){
		       	if(cnx(no)[d]<=xx[d])		continue;
			if(cndef(g3eref(nn,d))) ref=1;
		}else{
			e =croq(d,0);
			if(cnund(qq=cnq(nn,e)))		continue;
		       	if(cnx(cqn(qq,e))[d]<=xx[d])	continue;
			if(g3qref2(qq,d)) ref=1;
		}
	}
 }
/* anisotropic refinement in positive directions: */
 for(d=0;d<cgdim;d++){
 	if(cnund(no=cndr(nn,d)))	continue;
	while(ibgrefedge(cnd(nn),cnd(no))) if(cnund(no=g3eref(nn,d))) break;
 }
/* isotropic refinement (maxlength) in negative directions: */
 xx[0]=cnx(nn)[0]-lmax;
 xx[1]=cnx(nn)[1]-lmax;
 xx[2]=cnx(nn)[2]-lmax;
 while(ref){ref=0;
	for(d=0;d<cgdim;d++){
		if(cndef(no=cndr(nn,u=cgdim+d))){
		       	if(cnx(no)[d]>=xx[d])		continue;
			if(cndef(g3eref(no,d))) ref=1;
		}else{
			e=croq(u,0);
			if(cnund(qq=cnq(nn,e)))		continue;
		       	if(cnx(no=cqn(qq,e))[d]>=xx[d])	continue;
			if(g3qref2(qq,d)) ref=1;
		}
  	}
 }
/* geometric refinement tests for refined edges: */
 for(d=0;d<cgdim;d++){u=cgdim+d;
	if(cnund(no=cndr(nn,d))) continue;
	if(cnund(nu=cndr(nn,u))) continue;
	dx=cnx(no)[d]-cnx(nu)[d];
	qq=cnq(nn,croq(d,0)); /* 3 == MAXO-1 */
	for(o=0;o<MAXO3;o++){qo=qq;
		if(cnund(qq = cnq(nn,e=croq(d,o+1)))){
			if(cnund(qo))				continue;
			if(qo!=cnq(nn,croq(u,o)))		continue;
       		}else  	if(qq!=cnq(nn,croq(u,o+1)))		continue;
/* nu-no is a refined edge of qq --- start tests for qq: */
		if(qq!=qo){
		      	if(!cnnothing(cndr(nn,rv=cror(d,o))))	continue;
			nr=cndr(nu,rv);
			if(rv<cgdim){
				dv = rv;
				dy = (cnx(nr)[dv]-cnx(nu)[dv])/dx;
			}else{	dv = rv-cgdim;
				dy = (cnx(nu)[dv]-cnx(nr)[dv])/dx;
			}
			if(dy<ibggOref) g3elot(nn,rv,d);
			else if(dy>ibggOinv){
				nor=cndr(no,rv); ibgassert(cndef(nor));
				if(cnund(cndr(nor,u))) g3elot(nor,u,rv);
				nur=cndr(nu,rv); ibgassert(cndef(nur));
				if(cnund(cndr(nur,d))) g3elot(nur,d,rv);
			}
		}else{
			nr=cqn(qq,e);
		      	rv=cror(d,o);
			if(rv<cgdim){
				dv = rv;
				dy = (cnx(nr)[dv]-cnx(nu)[dv])/dx;
			}else{	dv = rv-cgdim;
				dy = (cnx(nu)[dv]-cnx(nr)[dv])/dx;
			}
		      	rq=cror(d,o+1);
			if(rq<cgdim){
				dq = rq;
				dz = (cnx(nr)[dq]-cnx(nu)[dq])/dx;
			}else{	dq = rq-cgdim;
				dz = (cnx(nu)[dq]-cnx(nr)[dq])/dx;
			}
			if(dy<ibggOref){
				if(dz<ibggOref)
			      		g3qref2(qq,d);
				else{
					g3qref2(qq,dq);
					g3elot(nn,rv,d);
				}
			}else if(dy*dy<ibggOref2*(1.0+4*dz*dz)){
				if(dz<ibggOref)
			      		g3qref2(qq,d);
				else{
					g3qref2(qq,dq);
					g3elot(nn,rv,d);
				}
			}else if(dy>ibggOinv){
				g3qref1(qq,dq);
			}
		}
	}
 }
/* refinement of parallel neighbour edges in the directions d and u:
	ibggPrefL - number of neighbours to control
	ibggPref  - maximal ratio of length of the edge to dx
*/
 for(d=0;d<cgdim;d++){
	if(cnund(no = cndr(nn,d))) continue;
	dx = cnx(no)[d] - cnx(nn)[d];
	u  = cgdim + d;
	e  = croq(d,0);
	for(i=0;i<ibggPrefL;i++){
		if(cnund(noo=cndr(no,d))){
		   	if(cnund(qq=cnq(no,e)))	break;
			noo=cqn(qq,e);
			while(cnx(noo)[d]-cnx(no)[d]>ibggPref*dx){
				if(!g3qref2(qq,d)) break;
				noo=cqn(qq=cnq(no,e),e);
				ibgassert(cndef(noo));
			}
  			if(cnund(noo=cndr(no,d))) noo=cqn(qq=cnq(no,e),e);
		}else{
			while(cnx(noo)[d]-cnx(no)[d]>ibggPref*dx){
				if(cnund(noo=g3eref(no,d))) break;
			}
		}
		no = noo;
	}
	nu = nn;
	e  = croq(u,0);
	for(i=0;i<ibggPrefL;i++){
		if(cnund(nuu=cndr(nu,u))){
		   	if(cnund(qq=cnq(nu,e)))	break;
			nuu=cqn(qq,e);
			while(cnx(nu)[d]-cnx(nuu)[d]>ibggPref*dx){
				if(!g3qref2(qq,d)) break;
				nuu=cqn(qq=cnq(nu,e),e);
				ibgassert(cndef(nuu));
			}
			if(cnund(nuu=cndr(nu,u))) nuu=cqn(qq=cnq(nu,e),e);
		}else{
			while(cnx(nu)[d]-cnx(nuu)[d]>ibggPref*dx){
				if(cnund(nuu=g3eref(nu,u))) break;
			}
		}
		nu = nuu;
	}
 }
 for(d=0;d<cgdim;d++){
	u  = cgdim + d;
	if(cnund(nu = cndr(nn,u))) continue;
	dx = cnx(nn)[d] - cnx(nu)[d];
	e  = croq(u,0);
	for(i=0;i<ibggPrefL;i++){
		if(cnund(nuu=cndr(nu,u))){
		   	if(cnund(qq=cnq(nu,e)))	break;
			nuu=cqn(qq,e);
			while(cnx(nu)[d]-cnx(nuu)[d]>ibggPref*dx){
				if(!g3qref2(qq,d)) break;
				nuu=cqn(qq=cnq(nu,e),e);
				ibgassert(cndef(nuu));
			}
  			if(cnund(nuu=cndr(nu,u))) nuu=cqn(qq=cnq(nu,e),e);
		}else{
			while(cnx(nu)[d]-cnx(nuu)[d]>ibggPref*dx){
				if(cnund(nuu=g3eref(nu,u))) break;
			}
		}
		nu = nuu;
	}
	no = nn;
	e  = croq(d,0);
	for(i=0;i<ibggPrefL;i++){
		if(cnund(noo=cndr(no,d))){
		   	if(cnund(qq=cnq(no,e)))	break;
			noo=cqn(qq,e);
			while(cnx(noo)[d]-cnx(no)[d]>ibggPref*dx){
				if(!g3qref2(qq,d)) break;
				noo=cqn(qq=cnq(no,e),e);
				ibgassert(cndef(noo));
			}
			if(cnund(noo=cndr(no,d))) noo=cqn(qq=cnq(no,e),e);
		}else{
			while(cnx(noo)[d]-cnx(no)[d]>ibggPref*dx){
				if(cnund(noo=g3eref(no,d))) break;
			}
		}
		no = noo;
	}
 }
}

static void g2nref(int nn)
{int ref=1,i,d,u,o,rv,dv,rq,dq,e,no,noo,nu,nuu,nr,nur,nor;
 ibgFloat dx,dy,dz;
 ibgFloat xx[ibgDIM],lmax;
 ibgrefreg(cnd(nn),&lmax);
/* isotropic refinement (maxlength) in positive directions: */
 xx[0]=cnx(nn)[0]+lmax;
 xx[1]=cnx(nn)[1]+lmax;
 while(ref){ref=0;
	for(d=0;d<cgdim;d++){
		if(cndef(no=cndr(nn,d))){
		       	if(cnx(no)[d]<=xx[d])		continue;
			if(cndef(g2eref(nn,d))) ref=1;
		}
	}
 }
/* anisotropic refinement in positive directions: */
 for(d=0;d<cgdim;d++){
 	if(cnund(no=cndr(nn,d)))	continue;
	while(ibgrefedge(cnd(nn),cnd(no))) if(cnund(no=g2eref(nn,d))) break;
 }
/* isotropic refinement (maxlength) in negative directions: */
 xx[0]=cnx(nn)[0]-lmax;
 xx[1]=cnx(nn)[1]-lmax;
 while(ref){ref=0;
	for(d=0;d<cgdim;d++){
		if(cndef(no=cndn(nn,d))){
		       	if(cnx(no)[d]>=xx[d])		continue;
			if(cndef(g2eref(no,d))) ref=1;
		}
  	}
 }
/* geometric refinement tests for refined edges: */
 for(d=0;d<cgdim;d++){u=cgdim+d;
	if(cnund(no=cndr(nn,d))) continue;
	if(cnund(nu=cndr(nn,u))) continue;
	dx=cnx(no)[d]-cnx(nu)[d];
	for(o=0;o<MAXO2;o++){
	      	if(!cnnothing(cndr(nn,rv=cror(d,o))))	continue;
		nr=cndr(nu,rv);
		if(rv<cgdim){
			dv = rv;
			dy = (cnx(nr)[dv]-cnx(nu)[dv])/dx;
		}else{	dv = rv-cgdim;
			dy = (cnx(nu)[dv]-cnx(nr)[dv])/dx;
		}
		if(dy<ibggOref) g2elot(nn,rv,d);
		else if(dy>ibggOinv){
			nor=cndr(no,rv); ibgassert(cndef(nor));
			if(cnund(cndr(nor,u))) g2elot(nor,u,rv);
			nur=cndr(nu,rv); ibgassert(cndef(nur));
			if(cnund(cndr(nur,d))) g2elot(nur,d,rv);
		}
	}
 }
/* refinement of parallel neighbour edges in the directions d and u:
	ibggPrefL - number of neighbours to control
	ibggPref  - maximal ratio of length of the edge to dx
*/
 for(d=0;d<cgdim;d++){
	if(cnund(no = cndr(nn,d))) continue;
	dx = cnx(no)[d] - cnx(nn)[d];
	u  = cgdim + d;
	for(i=0;i<ibggPrefL;i++){
		if(cnund(noo=cndr(no,d))){
		}else{
			while(cnx(noo)[d]-cnx(no)[d]>ibggPref*dx){
				if(cnund(noo=g2eref(no,d))) break;
			}
		}
		no = noo;
	}
	nu = nn;
	for(i=0;i<ibggPrefL;i++){
		if(cnund(nuu=cndr(nu,u))){
		}else{
			while(cnx(nu)[d]-cnx(nuu)[d]>ibggPref*dx){
				if(cnund(nuu=g2eref(nu,u))) break;
			}
		}
		nu = nuu;
	}
 }
 for(d=0;d<cgdim;d++){
	u  = rnegat(d);
	if(cnund(nu = cndr(nn,u))) continue;
	dx = cnx(nn)[d] - cnx(nu)[d];
	for(i=0;i<ibggPrefL;i++){
		if(cnund(nuu=cndr(nu,u))){
		}else{
			while(cnx(nu)[d]-cnx(nuu)[d]>ibggPref*dx){
				if(cnund(nuu=g2eref(nu,u))) break;
			}
		}
		nu = nuu;
	}
	no = nn;
	for(i=0;i<ibggPrefL;i++){
		if(cnund(noo=cndr(no,d))){
		}else{
			while(cnx(noo)[d]-cnx(no)[d]>ibggPref*dx){
				if(cnund(noo=g2eref(no,d))) break;
			}
		}
		no = noo;
	}
 }
}

static void g1nref(int nn)
{int ref=1,i,q,d=0,u,o,rv,dv,rq,dq,e,no,noo,nu,nuu,nr,nur,nor,qq,qo;
 ibgFloat dx,dy,dz;
 ibgFloat xx[ibgDIM],lmax;
/* isotropic refinement (maxlength) in positive directions: */
 ibgrefreg(cnd(nn),&lmax);
 xx[0]=cnx(nn)[0]+lmax;
 while(ref){ref=0;
	if(cndef(no=cndp(nn,0))){
	       	if(cnx(no)[d]<=xx[0])		continue;
		if(g1eref(nn,0)) ref=1;
	}
 }
/* anisotropic refinement in positive directions: */
 if(cndef(no=cndp(nn,0))){
 	while(ibgrefedge(cnd(nn),cnd(no))) if(cnund(no=g1eref(nn,0))) break;
 }
/* isotropic refinement (maxlength) in negative directions: */
 xx[0]=cnx(nn)[0]-lmax;
 while(ref){ref=0;
	if(cndef(no=cndn(nn,0))){
	       	if(cnx(no)[0]>=xx[0])		continue;
		if(g1eref(no,0)) ref=1;
  	}
 }
/* refinement of parallel neighbour edges in the directions d and u:
	ibggPrefL - number of neighbours to control
	ibggPref  - maximal ratio of length of the edge to dx
*/
 if(cndef(no = cndr(nn,0))){
	dx = cnx(no)[0] - cnx(nn)[0];
	u  = rnegat(d);
	for(i=0;i<ibggPrefL;i++){
		if(cndef(noo=cndr(no,d))){
			while(cnx(noo)[d]-cnx(no)[d]>ibggPref*dx){
				if(cnund(noo=g1eref(no,d))) break;
			}
		}
		no = noo;
	}
	nu = nn;
	for(i=0;i<ibggPrefL;i++){
		if(cndef(nuu=cndr(nu,u))){
			while(cnx(nu)[d]-cnx(nuu)[d]>ibggPref*dx){
				if(cnund(nuu=g1eref(nu,u))) break;
			}
		}
		nu = nuu;
	}
 }
 u  = rnegat(d);
 if(cndef(nu = cndr(nn,u))){
	dx = cnx(nn)[d] - cnx(nu)[d];
	for(i=0;i<ibggPrefL;i++){
		if(cnund(nuu=cndr(nu,u))){
		}else{
			while(cnx(nu)[d]-cnx(nuu)[d]>ibggPref*dx){
				if(cnund(nuu=g1eref(nu,u))) break;
			}
		}
		nu = nuu;
	}
	no = nn;
	for(i=0;i<ibggPrefL;i++){
		if(cnund(noo=cndr(no,d))){
		}else{
			while(cnx(noo)[d]-cnx(no)[d]>ibggPref*dx){
				if(cnund(noo=g1eref(no,d))) break;
			}
		}
		no = noo;
	}
 }
}

/* Delaunay-Triangulierung des Gitters */
static ibgFloat big;

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -