ibgdefault.c

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

 sb[0] = slast;
 sb[1] = slast + slast;
 sb[2] = slast + slast + slast;
 for(s0=0;s0<slast;s0++){
	sr[sb[0]+s0] = sr[sb[1]+s0] = sr[sb[2]+s0] = sr[s0];
	sn[sb[0]+s0][0] = sn[s0][0];
	sn[sb[1]+s0][1] = sn[s0][1];
	sn[sb[2]+s0][2] = sn[s0][2];
  	for(t0=0;t0<3;t0++){
		t1=(t0+1)%3; t2=(t0+2)%3;
	    	if((st=ss[s0][t0])>s0){n = &point[nlast++];
			for(d=0;d<ibgDIM;d++)
				ibgpX(*n)[d] = (ibgpX(*(sn[s0][t1]))[d]
					        +  ibgpX(*(sn[s0][t2]))[d])/2;
			ibgiRegionOfPoint(g,n,sn[s0][t1]);
		}else{
			for(i0=0;i0<3;i0++) if(ss[st][i0]==s0) break;
	       		n = sn[st][i0];
			i1=(i0+1)%3; i2=(i0+2)%3;
		      	ss[sb[t1]+s0][t0] = sb[i2]+st;
			ss[sb[t2]+s0][t0] = sb[i1]+st;
		      	ss[sb[i1]+st][t0] = sb[t2]+s0;
		      	ss[sb[i2]+st][t0] = sb[t1]+s0;
		}
	      	sn[s0][t0] = sn[sb[t1]+s0][t2] = sn[sb[t2]+s0][t1] = n;
	}
  	ss[sb[0]+s0][0] = s0;
	ss[sb[1]+s0][1] = s0;
	ss[sb[2]+s0][2] = s0;
 }
 for(s0=0;s0<slast;s0++){
  	ss[s0][0] = sb[0]+s0;
  	ss[s0][1] = sb[1]+s0;
  	ss[s0][2] = sb[2]+s0;
 }
 slast *= 4;
 if(first){first=0;
/* Unterteilung des Wuerfels */
	n23 = sn[0][0];
	n13 = sn[0][1];
	n12 = sn[0][2];
	n24 = sn[1][1];
	n34 = sn[1][2];
	n14 = sn[2][1];
 }
 if(dd>ibgDMaxNode) goto nodehalf;

 /* Durchlauf durch die Kanten zum Finden von Schnitten */
 ul = ibgpSegment(*nline);
 for(s0=0;s0<slast;s0++){
  	for(t0=0;t0<3;t0++){
	    	if((st=ss[s0][t0])>s0) continue;
		t1=(t0+1)%3; t2=(t0+2)%3;
		nn1 = sn[s0][t1];
		nn2 = sn[s0][t2];
		if(ibgpSegment(*n1)==ibgpSegment(*n2)) continue;
		rc = ibgiFaceWithEdge(g,&nface,nn1,nn2);
		rc = ibgiLineWithTriangle(g,nint,&nface,nn1,nn2,sn[s0][t0]);
		if(ul==ibgpSegment(*nint)){
			if(ibgpeq(g,nint,nline)) return sr[s0];
		}
		for(i0=0;i0<3;i0++) if(ss[st][i0]==s0) break;
		rc = ibgiLineWithTriangle(g,nint,&nface,nn2,nn1,sn[st][i0]);
		if(ul==ibgpSegment(*nint)){
			if(ibgpeq(g,nint,nline)) return sr[s0];
		}
	}
 }
/* refine the tetrahedron: */
nodehalf:
 if(ibgNodeBorderRef0){
	ibgNodeBorderRef0--; goto nodebordertest;
 }
 for(d=0;d<ibgDIM;d++){
	x1[d] = ibgpX(*n23)[d]   - ibgpX(*n4)[d];
	x2[d] = ibgpX(*n13)[d]   - ibgpX(*n4)[d];
	x3[d] = ibgpX(*n12)[d]   - ibgpX(*n4)[d];
	x0[d] = ibgpX(*nline)[d] - ibgpX(*n4)[d];
 }
 d1 = x0[0]*(x2[1]*x3[2]-x3[1]*x2[2])
    + x0[1]*(x2[2]*x3[0]-x3[2]*x2[0])
    + x0[2]*(x2[0]*x3[1]-x3[0]*x2[1]);
 d2 = x0[0]*(x3[1]*x1[2]-x1[1]*x3[2])
    + x0[1]*(x3[2]*x1[0]-x1[2]*x3[0])
    + x0[2]*(x3[0]*x1[1]-x1[0]*x3[1]);
 d3 = x0[0]*(x1[1]*x2[2]-x1[1]*x2[2])
    + x0[1]*(x1[2]*x2[0]-x1[2]*x2[0])
    + x0[2]*(x1[0]*x2[1]-x1[0]*x2[1]);
 if(d1>0)	s = -100;
 else if(d2>0)	s = -200;
 else if(d3>0)	s = -300;
 else 		s = -400;
 for(i=0;i<ibgdVSteps;i++){
	switch(s){
	case -100:
		sh=ibgDefaultNodeInTetrahedron(g,nint,nline,n1,n12,n13,n14);
		switch(sh){
			case ibgiOnSide123:	s = ibgiOnSide123;  break;
			case ibgiOnSide124:	s = ibgiOnSide124;  break;
			case ibgiOnSide134:	s = ibgiOnSide134;  break;
			case ibgiOnSide234:	s = -91;  break;
		}break;
	case -200:
		sh=ibgDefaultNodeInTetrahedron(g,nint,nline,n12,n2,n23,n24);
		switch(sh){
			case ibgiOnSide123:	s = ibgiOnSide123;  break;
			case ibgiOnSide124:	s = ibgiOnSide124;  break;
			case ibgiOnSide134:	s = -92;  break;
			case ibgiOnSide234:	s = ibgiOnSide234;  break;
		}break;
	case -300:
		sh=ibgDefaultNodeInTetrahedron(g,nint,nline,n13,n23,n3,n34);
		switch(sh){
			case ibgiOnSide123:	s = ibgiOnSide123;  break;
			case ibgiOnSide124:	s = -93;  break;
			case ibgiOnSide134:	s = ibgiOnSide134;  break;
			case ibgiOnSide234:	s = ibgiOnSide234;  break;
		}break;
	case -400:
		sh=ibgDefaultNodeInTetrahedron(g,nint,nline,n23,n13,n12,n24);
		switch(sh){
			case ibgiOnSide123:	s = ibgiOnSide123;  break;
			case ibgiOnSide124:	s = -23;  break;
			case ibgiOnSide134:	s = -20;  break;
			case ibgiOnSide234:	s = -21;  break;
		}break;
	case -19:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n14,n13,n12,n1);
		switch(sh){
			case ibgiOnSide123:	s = -91;  break;
			case ibgiOnSide124:	s = ibgiOnSide134;  break;
			case ibgiOnSide134:	s = ibgiOnSide124;  break;
			case ibgiOnSide234:	s = ibgiOnSide123;  break;
		}break;
	case -29:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n12,n23,n24,n2);
		switch(sh){
			case ibgiOnSide123:	s = -92;  break;
			case ibgiOnSide124:	s = ibgiOnSide123;  break;
			case ibgiOnSide134:	s = ibgiOnSide124;  break;
			case ibgiOnSide234:	s = ibgiOnSide234;  break;
		}break;
	case -39:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n34,n23,n13,n3);
		switch(sh){
			case ibgiOnSide123:	s = -93;  break;
			case ibgiOnSide124:	s = ibgiOnSide234;  break;
			case ibgiOnSide134:	s = ibgiOnSide134;  break;
			case ibgiOnSide234:	s = ibgiOnSide123;  break;
		}break;
	case -49:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n14,n24,n34,n4);
		switch(sh){
			case ibgiOnSide123:	s = -94;  break;
			case ibgiOnSide124:	s = ibgiOnSide124;  break;
			case ibgiOnSide134:	s = ibgiOnSide134;  break;
			case ibgiOnSide234:	s = ibgiOnSide234;  break;
		}break;
	case -91:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n12,n13,n14,n24);
		switch(sh){
			case ibgiOnSide123:	s = -19;  break;
			case ibgiOnSide124:	s = -12;  break;
			case ibgiOnSide134:	s = ibgiOnSide124;  break;
			case ibgiOnSide234:	s = -14;  break;
		}break;
	case -21:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n24,n13,n12,n14);
		switch(sh){
			case ibgiOnSide123:	s = -12;  break;
			case ibgiOnSide124:	s = -14;  break;
			case ibgiOnSide134:	s = ibgiOnSide124;  break;
			case ibgiOnSide234:	s = -19;  break;
		}break;
	case -41:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n24,n14,n13,n12);
		switch(sh){
			case ibgiOnSide123:	s = -14;  break;
			case ibgiOnSide124:	s = ibgiOnSide124;  break;
			case ibgiOnSide134:	s = -12;  break;
			case ibgiOnSide234:	s = -19;  break;
		}break;
	case -92:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n24,n23,n12,n13);
		switch(sh){
			case ibgiOnSide123:	s = -29;  break;
			case ibgiOnSide124:	s = -23;  break;
			case ibgiOnSide134:	s = -21;  break;
			case ibgiOnSide234:	s = ibgiOnSide123;  break;
		}break;
	case -12:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n24,n13,n12,n23);
		switch(sh){
			case ibgiOnSide123:	s = -21;  break;
			case ibgiOnSide124:	s = -23;  break;
			case ibgiOnSide134:	s = -29;  break;
			case ibgiOnSide234:	s = ibgiOnSide123;  break;
		}break;
	case -32:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n13,n23,n24,n12);
		switch(sh){
			case ibgiOnSide123:	s = -23;  break;
			case ibgiOnSide124:	s = ibgiOnSide123;  break;
			case ibgiOnSide134:	s = -21;  break;
			case ibgiOnSide234:	s = -29;  break;
		}break;
	case -93:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n13,n23,n34,n24);
		switch(sh){
			case ibgiOnSide123:	s = -39;  break;
			case ibgiOnSide124:	s = -32;  break;
			case ibgiOnSide134:	s = -34;  break;
			case ibgiOnSide234:	s = ibgiOnSide234;  break;
		}break;
	case -23:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n24,n23,n13,n34);
		switch(sh){
			case ibgiOnSide123:	s = -32;  break;
			case ibgiOnSide124:	s = ibgiOnSide234;  break;
			case ibgiOnSide134:	s = -34;  break;
			case ibgiOnSide234:	s = -39;  break;
		}break;
	case -43:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n34,n24,n13,n23);
		switch(sh){
			case ibgiOnSide123:	s = -34;  break;
			case ibgiOnSide124:	s = ibgiOnSide234;  break;
			case ibgiOnSide134:	s = -39;  break;
			case ibgiOnSide234:	s = -32;  break;
		}break;
	case -94:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n34,n24,n14,n13);
		switch(sh){
			case ibgiOnSide123:	s = -49;  break;
			case ibgiOnSide124:	s = -43;  break;
			case ibgiOnSide134:	s = ibgiOnSide134;  break;
			case ibgiOnSide234:	s = -41;  break;
		}break;
	case -14:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n13,n14,n24,n34);
		switch(sh){
			case ibgiOnSide123:	s = -41;  break;
			case ibgiOnSide124:	s = ibgiOnSide134;  break;
			case ibgiOnSide134:	s = -43;  break;
			case ibgiOnSide234:	s = -49;  break;
		}break;
	case -34:
		sh=ibgDefaultNodeInTetrahedron(g,nint,&nold,n13,n24,n34,n14);
		switch(sh){
			case ibgiOnSide123:	s = -43;  break;
			case ibgiOnSide124:	s = -41;  break;
			case ibgiOnSide134:	s = ibgiOnSide134;  break;
			case ibgiOnSide234:	s = -49;  break;
		}break;
	}
	if(sh == ibgNodeFound){
		return ibgNodeFound;
  	}else if(s > 0){
	      	return s;
	}
	nold = *nint;
 }
 ibgpSegment(*nint) = ibgDefaultNodeNr; 
 ibgpType(*nint) = ibgSNothing;
 return ibgNotFound;
}