📄 ibgg0.c
字号:
we regularise by refinement and goto nextbreg;
*/
upper=upper0; lower = 0;
if(cndef(nur)){
nuu=nn;
if(rr<cgdim){dr=rr; dy=cnx(nur)[dr]-cnx(nuu)[dr];}
else {dr=rr-cgdim; dy=cnx(nuu)[dr]-cnx(nur)[dr];}
if(dy>ibgtang){
if(cndef(ibgg0.eref(nuu,rr))) goto nextbreg;
}
if(dy*dy<ibggBLObtuse*dx*(dl-dx)){
/* nn shifted along the longer edge and gives an obtuse angle: */
if(cndef(ibgg0.eref(nn,ro))) goto nextbreg;
}
if(cnu(nur)!=nm){
if(cnb(nur)==0){
/* the side nuu-nur intersects a boundary that is not detected.
there can be an obtuse angle on the other side. */
if(dy<dl) goto changeface;
if(cnund(cndr(nuu,ru))){
if(cndef(ibgg0.elot(nuu,ru,rr))) goto nextbreg;
}
if(cnund(cndr(nur,ru))){
if(cndef(ibgg0.elot(nur,ru,rr))) goto nextbreg;
if(cgdim==3){
if(g3qref2(cnq(nuu,croq(ru,o)),rr))
goto nextbreg;
}
ibgmessage(ibgMWGridCoarse);
}
if(ibggBctg*dy<dx){
if(dy>1.5*dl||cndr(nur,rr)==outside){
if(cndef(ibgg0.eref(nuu,rr))) goto nextbreg;
}
bb=g0BNew; if(bb==0) return 0;
rc = ibgiFaceWithEdge(ibGeo,cbd(bb),cnd(nur),cnd(nuu));
ibgassert(rc==ibgFaceFound);
if(rr<cgdim){
ibgassert(cbx(bb)[dr]>=cnx(nuu)[dr]-ibGeo->Delta);
ibgassert(cbx(bb)[dr]<=cnx(nur)[dr]+ibGeo->Delta);
if(cbx(bb)[dr]<cnx(nuu)[dr]+ibGeo->Delta)
cbx(bb)[dr]=cnx(nuu)[dr]+ibGeo->Delta;
if(cbx(bb)[dr]>cnx(nur)[dr]-ibGeo->Delta)
cbx(bb)[dr]=cnx(nur)[dr]-ibGeo->Delta;
}else{
ibgassert(cbx(bb)[dr]<=cnx(nuu)[dr]+ibGeo->Delta);
ibgassert(cbx(bb)[dr]>=cnx(nur)[dr]-ibGeo->Delta);
if(cbx(bb)[dr]>cnx(nuu)[dr]-ibGeo->Delta)
cbx(bb)[dr]=cnx(nuu)[dr]-ibGeo->Delta;
if(cbx(bb)[dr]<cnx(nur)[dr]+ibGeo->Delta)
cbx(bb)[dr]=cnx(nur)[dr]+ibGeo->Delta;
}
for(o1=0;o1<ibgDIM;o1++) if(o1!=dr){
ibgassert(cbx(bb)[o1]<cnx(nuu)[o1]+ibGeo->Delta);
ibgassert(cbx(bb)[o1]>cnx(nuu)[o1]-ibGeo->Delta);
cbx(bb)[o1]=cnx(nuu)[o1];
}
if(cndr(nur,rr)==outside) {g0BDel(bb);}
else{
ibgassert(cnb(nur)==0);
cnb(nur) = bb; cbn(bb) = nur;
/* potential bug - better put nur into g0breg-stack and return */
while(g0breg(bb));
goto nextbreg;
}
} /* of ctg*dy < dx */
}else if((cbt(cnb(nur))==ibgSFace) &&
(ibgiStatusOfEdge(ibGeo,cnd(nur),cbd(bb))==ibgIncorrect)){
br = cnb(nur);
ur = cbu(br);
if(ibgiStatusOfEdge(ibGeo,cnd(nn),cbd(br))==ibgIncorrect){
if(dy<dl) goto changeface;
if(cndef(ibgg0.eref(nuu,rr))) goto nextbreg;
}
if(cnund(cndr(nuu,ru))){
if(cndef(ibgg0.elot(nuu,ru,rr))) goto nextbreg;
}
}
} /* of cnu(nur)!=nm */
if(0){
changeface: rc = ibgiFaceWithEdge(ibGeo,cbd(bb),cnd(nn),cnd(nur));
ibgassert(rc==ibgFaceFound);
for(d=0;d<cgdim;d++){
if(d!=dr){
cbx(bb)[d] = cnx(nn)[d];
}else if(rispos(rr)){
if(cbx(bb)[dr]<cnx(nn)[dr]+ibGeo->Delta)
cbx(bb)[dr]=cnx(nn)[dr]+ibGeo->Delta;
if(cbx(bb)[dr]>cnx(nur)[dr]-ibGeo->Delta)
cbx(bb)[dr]=cnx(nur)[dr]-ibGeo->Delta;
}else{
if(cbx(bb)[dr]>cnx(nn)[dr]-ibGeo->Delta)
cbx(bb)[dr]=cnx(nn)[dr]-ibGeo->Delta;
if(cbx(bb)[dr]<cnx(nur)[dr]+ibGeo->Delta)
cbx(bb)[dr]=cnx(nur)[dr]+ibGeo->Delta;
}
}
goto nextbreg;
}
if(cnund(nor=cndr(no,rr))){
if(cnb(no)==0){
if(cndef(ibgg0.elot(no,rr,ro))) goto nextbreg;
ibgfatal;
}
noo = cndr(no, ro); ibgassert(cndef(noo));
nor = cndr(noo,rr); ibgassert(cndef(nor));
g0boundtest(upper,noo,bm);
}else noo = no;
}else{ nuu = cndr(nn,ru); ibgassert(cndef(nuu));
g0boundtest(lower,nn,bm);
/* the main edge is refined. */
du = cnx(nn)[dd]-cnx(nuu)[dd]; if(du<0) du = -du;
nur= cndr(nuu,rr); ibgassert(cndef(nur));
if(rispos(rr)){dr=rr; dy=cnx(nur)[dr]-cnx(nuu)[dr];}
else {dr=rr-cgdim; dy=cnx(nuu)[dr]-cnx(nur)[dr];}
if(dy>ibgtang){
if(cndef(ibgg0.eref(nuu,rr))) goto nextbreg;
}
if(dy*dy<ibggBLObtuse*(du+dx)*(dl-dx)){
/* an obtuse angle on the other side possible. */
if(cndef(ibgg0.elot(nn,rr,ro))) goto nextbreg;
}
noo = no; nor=cndr(noo,rr); ibgassert(cndef(nor));
g0boundtest(lower,nuu,bm);
}
g0boundtest(upper,nor,bm);
g0boundtest(lower,nur,bm);
/* now nuu, noo and nur are defined. Let's find the rest: */
if(cnund(nr=cndr(nur,ro))){
if(cnb(nur)==0){
if(cndef(ibgg0.elot(nur,ro,rr))) goto nextbreg; ibgfatal;
}
nur = cndr(nuh=nur,rr); ibgassert(cndef(nur));
g0boundtest(lower,nur,bm);
nr = cndr(nur,ro); ibgassert(cndef(nr));
}else nuh = nuu;
if(cnund( cndr(nor,ru))){
if(cnb(nor)==0){
if(cndef(ibgg0.elot(nor,ru,rr))) goto nextbreg; ibgfatal;
}
nor = cndr(noh=nor,rr); ibgassert(cndef(nor));
g0boundtest(upper,nor,bm);
}else noh = noo;
if(nor!=nr){
if(cnb(nr)==0){
/* es kann auch generell elot gemacht werden: */
if(!lower || dx > dl/2){
if(cndef(ibgg0.elot(nr,rother(rr),ro))) goto nextbreg;
ibgfatal;
}
}
g0boundtest(upper,nr,bm);
ibgassert(cndr(nr,ro)==nor);
}
g0boundtest(lower,nr,bm);
/* now all corners nuu, noo, nor, nur are defined. If there are refinement
points on the edges, they are boundary points (nn,nuh,nr,noh).
If a good intersection was found we can continue with another side. */
if(lower==1 || upper==1) continue;
/* potential bug: we really have to count the number of intersections. */
if(!lower){
g0boundtest(lower,nor,bm);
g0boundtest(lower,noh,bm); if(lower==1) continue;
}
if(!upper){
g0boundtest(upper,nur,bm);
g0boundtest(upper,nuh,bm); if(upper==1) continue;
}
/* the correct other intersection is not found. search for Edge-intersection */
b3=g0BNew; if(b3==0) return 0;
rc=ibgiLineWithRectangle(ibGeo,cbd(b3),cbd(bb)
,cnd(noo),cnd(nuu),cnd(nur),cnd(nor));
if(rc!=ibgLineFound) {g0BDel(b3); continue;}
if(ro<cgdim){
if(rr<cgdim) {n1=noo; n2=nuu; n3=nur; n4=nor;}
else {n1=nor; n2=nur; n3=nuu; n4=noo;}
}else{
if(rr<cgdim) {n1=nuu; n2=noo; n3=nor; n4=nur;}
else {n1=nur; n2=nor; n3=noo; n4=nuu;}
}
if(cnx(n1)[dd]-cnx(n2)[dd]>cnx(n3)[dr]-cnx(n2)[dr]) hoch=1;
else hoch=0;
/* find the quarter of the rectangle containing the Edge point: */
dx1=(cbx(b3)[dd]-cnx(n2)[dd])/(cnx(n1)[dd]-cnx(n2)[dd]);
dx2=(cbx(b3)[dr]-cnx(n2)[dr])/(cnx(n3)[dr]-cnx(n2)[dr]);
if(dx1<0.5){r1=dd;
if(dx2<0.5) {nv0=n2; nv1=n1; nv2=n3; nv3=n4; r2=dr;}
else {dx2=1-dx2; nv0=n3; nv1=n4; nv2=n2; nv3=n1; r2=dr+cgdim;}
}else{ r1=dd+cgdim; dx1=1-dx1;
if(dx2<0.5) {nv0=n1; nv1=n2; nv2=n4; nv3=n3; r2=dr;}
else {dx2=1-dx2; nv0=n4; nv1=n3; nv2=n1; nv3=n2; r2=dr+cgdim;}
}
/* find point to shift in the quarter: */
nv=nv0;
if((nvo=cndr(nv0,r1))!=nv1)
if(dx1>0.25) nv=nvo;
if((nvr=cndr(nv0,r2))!=nv2){
if(dx2>0.25){
if(nv==nvo){
if(hoch) nv=nvr;
}else nv=nvr;
}
}
if(cndr(nv,rother(r1))==outside){
nv = nvo;
if(cndr(nv,rother(r2))==outside){
nv = nv3;
}
}else if(cndr(nv,rother(r2))==outside){
nv = nvr;
}
/* point nv found, let's shift him: */
ibgassert(cndef(nv));
bh = cnb(nv);
cnb(nv) = b3; cbn(b3) = nv;
if(bh!=0){
ibgassert(cbt(bh)==ibgSFace);
{g0BDel(bh);}
}
if(rispos(r1)){
ibgassert(cbx(b3)[dd]>=cnx(nv0)[dd]-ibGeo->Delta);
if(cbx(b3)[dd]<cnx(nv0)[dd]+ibGeo->Delta) cbx(b3)[dd]=cnx(nv0)[dd]+ibGeo->Delta;
}else{
ibgassert(cbx(b3)[dd]<=cnx(nv0)[dd]+ibGeo->Delta);
if(cbx(b3)[dd]>cnx(nv0)[dd]-ibGeo->Delta) cbx(b3)[dd]=cnx(nv0)[dd]-ibGeo->Delta;
}
if(rispos(r2)){
ibgassert(cbx(b3)[dr]>=cnx(nv0)[dr]-ibGeo->Delta);
if(cbx(b3)[dr]<cnx(nv0)[dr]+ibGeo->Delta) cbx(b3)[dr]=cnx(nv0)[dr]+ibGeo->Delta;
}else{
ibgassert(cbx(b3)[dr]<=cnx(nv0)[dr]+ibGeo->Delta);
if(cbx(b3)[dr]>cnx(nv0)[dr]-ibGeo->Delta) cbx(b3)[dr]=cnx(nv0)[dr]-ibGeo->Delta;
}
if(cbx(b3)[dd]>cnx(nv)[dd]){
if(cbx(b3)[dd]<cnx(nv)[dd]+ibGeo->Delta) cbx(b3)[dd]=cnx(nv)[dd]+ibGeo->Delta;
}else{
if(cbx(b3)[dd]>cnx(nv)[dd]-ibGeo->Delta) cbx(b3)[dd]=cnx(nv)[dd]-ibGeo->Delta;
}
if(cbx(b3)[dr]>cnx(nv)[dr]){
if(cbx(b3)[dr]<cnx(nv)[dr]+ibGeo->Delta) cbx(b3)[dr]=cnx(nv)[dr]+ibGeo->Delta;
}else{
if(cbx(b3)[dr]>cnx(nv)[dr]-ibGeo->Delta) cbx(b3)[dr]=cnx(nv)[dr]-ibGeo->Delta;
}
for(di=0;di<cgdim;di++){
if(di!=dd && di!=dr){
ibgassert(cbx(b3)[di]<cnx(nv)[di]+ibGeo->Delta);
ibgassert(cbx(b3)[di]>cnx(nv)[di]-ibGeo->Delta);
cbx(b3)[di] = cnx(nv)[di];
}
}
/* potential bug - better put nv into g0breg-stack and return */
while(g0breg(b3));
if(bh==bb) return 0;
goto nextbreg;
}
return 0;
case ibgSLine:
ibgrefline(cbd(bb),&ibgNorm,&ibgtang);
for(d=0;d<cgdim;d++){
if(cbx(bb)[d]!=cnx(nn)[d]){
continue;
}else{
ibgassert(cgdim==3);
for(o=0;o<MAXO3;o++){int qr;
if((qq=cnq(nn,croq(d,o)))==(qr=cnq(nn,croq(3+d,o)))){
if(qq>0){
if(g3qref2(qq,d)) goto nextbreg;
ibgfatal;
}
}
if(qq>0){
if(cnx(cqn(qq,croq(d,o)))[d]-cnx(nn)[d]>ibgtang){
if(g3qref2(qq,d)) goto nextbreg;
}else{
if(g3qref1(qq,d)) goto nextbreg;
}}
if(qr>0){
if(cnx(nn)[d]-cnx(cqn(qr,croq(3+d,o)))[d]>ibgtang){
if(g3qref2(qr,d)) goto nextbreg;
}else{
if(g3qref1(qr,d)) goto nextbreg;
}}
}
}
}
for(rr=0;rr<ibgg0.maxr;rr++){
d = rposit(rr);
if(cnund(no = cndr(nn,rr))) continue;
if(cnb(no)!=0){
bo = cnb(no);
del= 0;
dx = cbx(bo)[0]-cbx(bb)[0]; if(dx<0) del -= dx; else del += dx;
dx = cbx(bo)[1]-cbx(bb)[1]; if(dx<0) del -= dx; else del += dx;
dx = cbx(bo)[2]-cbx(bb)[2]; if(dx<0) del -= dx; else del += dx;
if(del<6*ibGeo->Delta){
/* wenn es hier knallt noch mal "uberlegen ob es so richtig ist: */
if(( rispos(rr)&&(cbx(bo)[rr]<cnx(no)[rr]))
|| (!rispos(rr)&&(cbx(bo)[ d]>cnx(no)[ d]))){
if(cbt(bo)==ibgSNode){
ibgassert(cnb(nn)==bb);
cnb(nn) = 0; g0BDel(bb);
g0bshift(nn,0); bb=bo;
}else{
ibgassert(cnb(no)==bo);
cnb(no) = 0; g0BDel(bo);
g0bshift(no,1);
}
goto nextbreg;
} /* else error??? */
}
}
}
for(d=0;d<cgdim;d++){
if(cnx(nn)[d]<cbx(bb)[d]){
if(cndef(no = cndp(nn,d)))
if(cnx(no)[d]-cnx(nn)[d]>ibgNorm)
{if(cndef(ibgg0.eref(nn,d))) goto nextbreg;}
if(cndef(nu = cndn(nn,d)))
if(cnx(nn)[d]-cnx(nu)[d]>ibgNorm)
{if(cndef(ibgg0.eref(nu,d))) goto nextbreg;}
}else if(cnx(nn)[d]>cbx(bb)[d]){
if(cndef(no = cndn(nn,d)))
if(cnx(nn)[d]-cnx(no)[d]>ibgNorm)
{if(cndef(ibgg0.eref(no,d))) goto nextbreg;}
if(cndef(nu = cndp(nn,d)))
if(cnx(nu)[d]-cnx(nn)[d]>ibgNorm)
{if(cndef(ibgg0.eref(nn,d))) goto nextbreg;}
}else continue;
for(o=0;o<ibgg0.maxo;o++){
ro = ibgg0.ror[d][o];
{register int d1 = rposit(ro);
if(cnx(nn)[d1]==cbx(bb)[d1]) continue;
}
if(cnund(cndr(nn,ro)))
{if(cnund(ibgg0.elot(nn,ro,d)))ibgwarning; else goto nextbreg;}
if(cndef(no)) if(cnund(cndr(no,ro)))
{if(cnund(ibgg0.elot(no,ro,d)))ibgwarning; else goto nextbreg;}
if(cndef(nu)) if(cnund(cndr(nu,ro)))
{if(cnund(ibgg0.elot(nu,ro,d)))ibgwarning; else goto nextbreg;}
}
}
return 0;
case ibgSNode:
return 0;
}
nextbreg:
if(bb<=g0bfirst && ibgg0.bstack[bb]==0){
if(g0bstb<=0) g0bstb = bb;
else{
ibgassert(ibgg0.bstack[g0bstl] == -1);
ibgg0.bstack[g0bstl] = bb;
}
ibgg0.bstack[bb] = -1; g0bstl = bb;
}
return 0;
moverflow: ibgfatal; return 0;
}
static int g1eref(int nu, unsigned ro)
{int b,d,dh,ru;
int no,nh,nur,nor,nr,po,pu;
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);}
/* 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)[d] = cnx(nu)[d];
#ifdef GLEVEL
cnl(nh)[dh] = cnl(nu)[dh];
#endif
}else{ cnx(nh)[dh] = 0.5*(cnx(nu)[dh]+cnx(no)[dh]);
cndp(nh,dh) = cndn(nh,dh) = nothing;
#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));
if(ibgg0.nb!=ibgNULL) {g0bjump(nh);}
return nh;
}/* end of g1eref */
ibgFloat g1cgh;
static ibgFloat g1cgo,g1cgu;
static int ibgg0New(ibGrid0 *g0)
{int n;
ibgg0= *g0; ibgg = *(g0->gg);
g0BNShift();
g0QFree();
g0CAlloc(ibgg.lastPoint);
ibgg.lastCell = ibgg.lastPoint;
for(n=1;n<=ibgg.lastPoint;n++){
ibgg.CPointFirst[n] = n;
ibgg.CSideFirst[n] = 0;
ibgg.CPoint[n] = n;
cct(n) = ibgc0Node;
ccu(n) = 0;
}
*ibgg0.gg = ibgg;
*g0 = ibgg0;
return ibgSuccess;
}
static int ibgg1New(ibGrid0 *g0)
{ibGrid *gg;
ibgPoint nd; ibgSegmentType t; int nn,bb,c,d,na,ne;
ibgFloat xm,dx,big;
int n,c0,*cn0,*cs0,nc;
ibgg0= *g0; ibgg = *(g0->gg);
g0BNShift();
g0QFree();
g0CAlloc(ibgg.lastPoint+10);
ibgg.lastCell=ibgg.lastPoint+1;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -