📄 ibgdefault.c
字号:
/* last edit: Ilja Schmelzer -------------- 18-OCT-1994 09:08:19.82 */
/************************************************************************/
/* */
/* <<< I B G >>> - Intersection - Based Grid generation package */
/* */
/* Version 1.1 by Ilja Schmelzer schmelzer@iaas-berlin.d400.de */
/* */
/* to be distributed under IBG license conditions (see "readme.ibg") */
/* */
/************************************************************************/
/* <<< IBGDefault >>> - Intersection-Based Geometry - Default functions */
#include "ibgdefault.h"
#include "ibglib.h"
#include "ibgi.h"
#include "ibgd0.h"
ibgFloat ibgDMaxFace = 0.1; /* distance for boundary test */
ibgFloat ibgDMaxLine = 0.1;
ibgFloat ibgDMaxNode = 0.1;
int ibgVBorderref = 0;
int ibgDefaultFree(ibGeometry g)
{
if(g != (ibGeometry)ibgNULL) free(g);
return ibgSuccess;
}
/* bisection algorithm based on the Region-function: */
int ibgDefaultFaceWithEdge(ibGeometry g, ibgPoint *nint,
ibgPoint *n1, ibgPoint *n2)
{ibgPoint c1,c2,cm,cf;
ibgFloat delta=0;
int m1=ibgpSegment(*n1), m2=ibgpSegment(*n2), d, rc;
c1= *n1; c2= *n2;
for(d=0;d<ibgDIM;d++){
register ibgFloat dx = ibgpX(c1)[d] - ibgpX(c2)[d];
delta += dx>0 ? dx : -dx;
}
while(delta>g->Delta){ delta /= 2;
for(d=0;d<ibgDIM;d++){
ibgpX(cm)[d] = 0.5*(ibgpX(c1)[d] + ibgpX(c2)[d]);
}
ibgiRegionOfPoint(g,&cm,&c1);
if (ibgpSegment(cm)==m1){
if(delta>ibgDMaxFace){
rc = ibgDefaultFaceWithEdge(g,nint,&c1,&cm);
if(rc==ibgFaceFound){return ibgFaceFound;}
rc = ibgDefaultFaceWithEdge(g,nint,&cm,&c2);
return rc;
}else if(m1==m2){
ibgpType(*nint) = ibgSNothing;
return ibgNotFound;
}
c1 = cm;
}else if(ibgpSegment(cm)==m2){
c2 = cm;
}else {
c2 = cm;
m2 = ibgpSegment(cm);
}
}
for(d=0;d<ibgDIM;d++){
ibgpX(*nint)[d] = 0.5*(ibgpX(c1)[d] + ibgpX(c2)[d]);
}
ibgpSegment(*nint)=ibgDefaultFaceNr(m1,m2);
ibgpType(*nint)=ibgSFace;
return ibgFaceFound;
}
int ibgDefaultLineWithTriangle(ibGeometry g,
ibgPoint *nint, ibgPoint *nface,
ibgPoint *n1, ibgPoint *n2, ibgPoint *n3)
{ibgPoint n12,n23,n31,nold;
int mf,m1,m2,m3,m4,b1,i,d,s,sn;
ibgFloat dd,xm,xp;
dd = 0;
#define max3(x1,x2,x3) x2>x3 ? x1>x2? x1:x2 : x1>x3? x1:x3
#define min3(x1,x2,x3) x2<x3 ? x1<x2? x1:x2 : x1<x3? x1:x3
for(d=0;d<ibgDIM;d++) {
xp = max3(ibgpX(*n1)[d],ibgpX(*n2)[d],ibgpX(*n3)[d]);
xm = min3(ibgpX(*n1)[d],ibgpX(*n2)[d],ibgpX(*n3)[d]);
dd += xp - xm;
}
if(dd < g->Delta) {
for(d=0;d<ibgDIM;d++) {
ibgpX(*nint)[d] = (ibgpX(*n1)[d]+ibgpX(*n2)[d]+ibgpX(*n3)[d])/3;
}
ibgiRegionOfPoint(g,nint,n1);
ibgpSegment(*nint) = ibgDefaultLineNr;
ibgpType(*nint) = ibgSLine;
return ibgLineFound;
}
if(dd > ibgDMaxLine) goto linehalf;
mf=ibgpSegment(*nface); m1=ibgpSegment(*n1); m2=ibgpSegment(*n2);
if(ibgIncorrect!=ibgiStatusOfEdge(g,n1,nface)){
if((m3=ibgpSegment(*n3))==m1){
if(m2==m1) goto linehalf;
ibgiFaceWithEdge(g,nint,n3,n2);
if(ibgpSegment(*nint)==mf) return ibgiOnSide23;
goto linehalf;
}
ibgiFaceWithEdge(g,nint,n1,n3);
if(ibgpSegment(*nint)==mf) return ibgiOnSide31;
goto linehalf;
}else if(ibgIncorrect!=ibgiStatusOfEdge(g,n2,nface)){
if((m3=ibgpSegment(*n3))==m2){
if(m1==m2) goto linehalf;
ibgiFaceWithEdge(g,nint,n3,n1);
if(ibgpSegment(*nint)==mf) return ibgiOnSide31;
goto linehalf;
}
b1=ibgiFaceWithEdge(g,nint,n2,n3);
if(ibgpSegment(*nint)==mf) return ibgiOnSide23;
goto linehalf;
}
linehalf:
for(d=0;d<ibgDIM;d++) {
ibgpX(n12)[d] = (ibgpX(*n1)[d]+ibgpX(*n2)[d])/2;
ibgpX(n23)[d] = (ibgpX(*n2)[d]+ibgpX(*n3)[d])/2;
ibgpX(n31)[d] = (ibgpX(*n3)[d]+ibgpX(*n1)[d])/2;
}
ibgiRegionOfPoint(g,&n12,n1);
ibgiRegionOfPoint(g,&n23,n2);
ibgiRegionOfPoint(g,&n31,n3);
dd=0;
for(d=0;d<ibgDIM;d++) {
dd += (ibgpX(*nface)[d] - ibgpX(n12)[d])
* (ibgpX( n12)[d] - ibgpX(*n1)[d]);
}
if(dd<0) s= -1; else s= -2;
for(i=0;i<ibgdLSteps;i++){
switch(s){
case -1:
sn=ibgDefaultLineWithTriangle(g,nint,nface,n1,&n12,&n31);
switch(sn){
case ibgiOnSide12: s = ibgiOnSide12; break;
case ibgiOnSide23: s = -10; break;
case ibgiOnSide31: s = ibgiOnSide31; break;
}break;
case -2:
sn=ibgDefaultLineWithTriangle(g,nint,nface,&n12,n2,&n23);
switch(sn){
case ibgiOnSide12: s = ibgiOnSide12; break;
case ibgiOnSide23: s = ibgiOnSide23; break;
case ibgiOnSide31: s = -20; break;
}break;
case -23:
sn=ibgDefaultLineWithTriangle(g,nint,&nold,&n12,&n31,n1);
switch(sn){
case ibgiOnSide12: s = -10; break;
case ibgiOnSide23: s = ibgiOnSide31; break;
case ibgiOnSide31: s = ibgiOnSide12; break;
}break;
case -31:
sn=ibgDefaultLineWithTriangle(g,nint,&nold,&n23,&n12,n2);
switch(sn){
case ibgiOnSide12: s = -20; break;
case ibgiOnSide23: s = ibgiOnSide12; break;
case ibgiOnSide31: s = ibgiOnSide23; break;
}break;
case -12:
sn=ibgDefaultLineWithTriangle(g,nint,&nold,&n31,&n23,n3);
switch(sn){
case ibgiOnSide12: s = -30; break;
case ibgiOnSide23: s = ibgiOnSide23; break;
case ibgiOnSide31: s = ibgiOnSide31; break;
}break;
case -10:
sn=ibgDefaultLineWithTriangle(g,nint,&nold,&n31,&n12,&n23);
switch(sn){
case ibgiOnSide12: s = -23; break;
case ibgiOnSide23: s = -31; break;
case ibgiOnSide31: s = -12; break;
}break;
case -20:
sn=ibgDefaultLineWithTriangle(g,nint,&nold,&n12,&n23,&n31);
switch(sn){
case ibgiOnSide12: s = -31; break;
case ibgiOnSide23: s = -12; break;
case ibgiOnSide31: s = -23; break;
}break;
case -30:
sn=ibgDefaultLineWithTriangle(g,nint,&nold,&n23,&n31,&n12);
switch(sn){
case ibgiOnSide12: s = -12; break;
case ibgiOnSide23: s = -23; break;
case ibgiOnSide31: s = -31; break;
}break;
}
if(sn == ibgLineFound){
return ibgLineFound;
}
if(s>0) return s;
nold = *nint;
}
ibgpSegment(*nint) = ibgDefaultLineNr;
ibgpType(*nint) = ibgSNothing;
return ibgNotFound;
}
static int ibgpeq(ibGeometry g, ibgPoint *n1,ibgPoint *n2)
{ibgFloat *x1,*x2;int d;
if(ibgpSegment(*n1)!=ibgpSegment(*n2)) return 1;
x1=ibgpX(*n1); x2=ibgpX(*n2);
for(d=0;d<ibgDIM;d++){register x;
if((x = x1[d]-x2[d]) > g->Delta) return 1;
if((x ) < -g->Delta) return 1;
}
return 0;
}
static int snnr[4][3] = {{1,2,3},{4,3,2},{1,3,4},{4,2,1}};
static int ssnr[4][3] = {{1,2,3},{0,3,2},{1,3,0},{0,2,1}};
static int snr[4] = {ibgiOnSide123,ibgiOnSide234,ibgiOnSide134,ibgiOnSide124};
int ibgDefaultNodeInTetrahedron(ibGeometry g,
ibgPoint *nint, ibgPoint *nline,
ibgPoint *n1, ibgPoint *n2, ibgPoint *n3, ibgPoint *n4)
{
#define MAXPOINTS 6
#define MAXTRIANGLES 16
ibgPoint point[MAXPOINTS],*nn1,*nn2,*n;
ibgPoint *(sn[MAXTRIANGLES][3]);
int (ss[MAXTRIANGLES][3]);
int (sr[MAXTRIANGLES]);
ibgPoint *(nn[5]);
ibgPoint *n12,*n13,*n14,*n23,*n24,*n34,nold,nface;
int d,i,sb[3],s0,t0,t1,t2,i0,i1,i2,st;
ibgSegment m00,m01,m10,m11,ul;
int ibgNodeBorderRef0=ibgVBorderref;
ibgFloat dd,d1,d2,d3,xm,xp;
ibgFloat x0[ibgDIM],x1[ibgDIM],x2[ibgDIM],x3[ibgDIM];
int first=1,slast,nlast,sh,s,rc;
/* Initialization */
/* diameter computation */
dd = 0;
#define max4(x1,x2,x3,x4) x3>x4? x2>x3 ? x1>x2? x1:x2 : x1>x3? x1:x3:\
x2>x4 ? x1>x2? x1:x2 : x1>x4? x1:x4
#define min4(x1,x2,x3,x4) x3<x4? x2<x3 ? x1<x2? x1:x2 : x1<x3? x1:x3:\
x2<x4 ? x1<x2? x1:x2 : x1<x4? x1:x4
for(d=0;d<ibgDIM;d++) {
xp = max4(ibgpX(*n1)[d],ibgpX(*n2)[d],ibgpX(*n3)[d],ibgpX(*n4)[d]);
xm = min4(ibgpX(*n1)[d],ibgpX(*n2)[d],ibgpX(*n3)[d],ibgpX(*n4)[d]);
dd += xp - xm;
}
if(dd < g->Delta) {
for(d=0;d<ibgDIM;d++) {
ibgpX(*nint)[d] = (ibgpX(*n1)[d]+ibgpX(*n2)[d]
+ ibgpX(*n3)[d]+ibgpX(*n4)[d])/4;
}
ibgiRegionOfPoint(g,nint,n1);
ibgpSegment(*nint) = ibgDefaultLineNr;
ibgpType(*nint) = ibgSNode;
return ibgNodeFound;
}
nn[1] = n1; nn[2] = n2; nn[3] = n3; nn[4] = n4;
for(s0=0;s0<4;s0++){
sr[s0] = snr[s0];
sn[s0][0] = nn[snnr[s0][0]];
sn[s0][1] = nn[snnr[s0][1]];
sn[s0][2] = nn[snnr[s0][2]];
ss[s0][0] = ssnr[s0][0];
ss[s0][1] = ssnr[s0][1];
ss[s0][2] = ssnr[s0][2];
}
slast = 4; nlast = 0;
nodebordertest:
/* weitere Verfeinerung des Randes */
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -