ibgdsplitr.c

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

 }else{	p23 = f3/(f3-f2);
	q23 = f2/(f2-f3);
 }
 if(p23<0.1)	{p23=0.1; q23=0.9;}
 if(q23<0.1)	{q23=0.1; p23=0.9;}
 if(f1*f3>=0){
	p31 = q31 = 0.5;
 }else{	p31 = f1/(f1-f3);
	q31 = f3/(f3-f1);
 }
 if(p31<0.1)	{p31=0.1; q31=0.9;}
 if(q31<0.1)	{q31=0.1; p31=0.9;}
 for(d=0;d<ibgDIM;d++){
	ibgpX(n12)[d]=(p12*ibgpX(*n1)[d]+q12*ibgpX(*n2)[d]);
	ibgpX(n23)[d]=(p23*ibgpX(*n2)[d]+q23*ibgpX(*n3)[d]);
	ibgpX(n31)[d]=(p31*ibgpX(*n3)[d]+q31*ibgpX(*n1)[d]);
 }
 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=SplitRegionLine0(splt,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=SplitRegionLine0(splt,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=SplitRegionLine0(splt,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=SplitRegionLine0(splt,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=SplitRegionLine0(splt,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=SplitRegionLine0(splt,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=SplitRegionLine0(splt,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=SplitRegionLine0(splt,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;
	}else if(sn== ibgSplitFound){
		return ibgSplitFound;
  	}else if(s>0){
	      	return s;
	}
	nold = *nint;
 }
 ibgpSegment(*nint) = ibgDefaultLineNr;
 ibgpType(*nint) = ibgSNothing;
 ibgfatal;	/* besser man weiss es */
 return ibgNotFound;
}

/* nbase is here a region point over nface (in rold or in old rnew-part) */
static int SplitRegionLine1(ibgSplitRegion splt,
	   	      	ibgPoint *nint, ibgPoint *nface, ibgPoint *nbase,
			ibgPoint *n1, ibgPoint *n2, ibgPoint *n3)
{int rc,s12,s23,s31;
 ibgPoint nh;
 rc = ibgiFaceWithEdge(splt->gold,&nh,nbase,n2);
 if(rc==ibgNotFound){ rc = ibgiFaceWithEdge(splt->gold,&nh,n2,n3);
  if(rc==ibgNotFound){ rc = ibgiFaceWithEdge(splt->gold,&nh,n3,n1);
   if(rc==ibgNotFound){ rc = ibgiFaceWithEdge(splt->gold,&nh,n1,nbase);
    if(rc==ibgNotFound){;
    }else{
	rc = SplitRegionLine0(splt,nint,&nh,n1,n2,n3);
   	s12 = ibgiOnSide12; s23 = ibgiOnSide23; s31 = ibgiOnSide31;
    }
   }else{
	rc = SplitRegionLine0(splt,nint,&nh,n3,n1,n2);
   	s12 = ibgiOnSide31; s23 = ibgiOnSide12; s31 = ibgiOnSide23;
   }
  }else{
	rc = SplitRegionLine0(splt,nint,&nh,n2,n3,n1);
   	s12 = ibgiOnSide23; s23 = ibgiOnSide31; s31 = ibgiOnSide12;
  }
 }else{
	rc = SplitRegionLine0(splt,nint,&nh,n1,n2,n3);
   	s12 = ibgiOnSide12; s23 = ibgiOnSide23; s31 = ibgiOnSide31;
 }
 switch(rc){
 case ibgiOnSide12:	return s12;
 case ibgiOnSide23:	return s23;
 case ibgiOnSide31:	return s31;
 default:		return rc;
 }
}

static int LineWithTriangle(ibGeometry This, ibgPoint *nint, ibgPoint *nface,
		const ibgPoint *n1, const ibgPoint *n2, const ibgPoint *n3)
{ibgPoint nh1,nh2,nh3,nhf,nhb,nh;
 ibgSplitRegion splt = (ibgSplitRegion) This;
 int rc,d;
 ibgSegment fu;
 ibgFloat dd,dr,dx;
 nh1 = *n1; nh2 = *n2; nh3 = *n3; nhf = *nface;
 if(ibgpSegment(*n1) == splt->rnew){
	rc=ibgiRegionOfPoint(splt->gold,&nh1,n1);
	if(rc==ibgError) ibgpSegment(nh1) = splt->rold;
 }
 if(ibgpSegment(*n2) == splt->rnew){
	rc=ibgiRegionOfPoint(splt->gold,&nh2,n2);
	if(rc==ibgError) ibgpSegment(nh2) = splt->rold;
 }
 if(ibgpSegment(*n3) == splt->rnew){
	rc=ibgiRegionOfPoint(splt->gold,&nh3,n3);
	if(rc==ibgError) ibgpSegment(nh3) = splt->rold;
 }
 fu = ibgpSegment(nhf);
 if(ibgCorrect==ibgtStatusRF(&splt->g.top,splt->rnew,fu)){
	if(fu==splt->face){
		rc=ibgiRegionOfPoint(splt->gold,&nhf,&nh1);
	 	/* search along the new boundary */
		rc = SplitRegionLine1(splt,nint,nface,&nhf,&nh1,&nh2,&nh3);
		if(rc==ibgNotFound){
			rc = ibgDefaultLineWithTriangle(This,
		     				nint,nface,n1,n2,n3);
		 	if(rc != ibgLineFound) return rc;
			ibgwarning;
		/* very strange - Default has found a line, Line1 not. */
		}
		if(rc==ibgSplitFound){
			ibgpSegment(*nint) =
				splitCreateLineNr(splt,ibgpSegment(*nint));
			ibgpType(*nint) = ibgSLine;
			return ibgLineFound;
		}
		return rc;
	}
	dd = 0;
	nhb = nh1;
facerecover:
	rc = ibgiFaceWithEdge(splt->gold,&nh,&nhb,&nh2);
	if(rc!=ibgFaceFound){
		if((dd>0) && (dd<20*splt->g.Delta)){
      	/* uncertain old geometry description!!! use nh instead of nhf */
			ibgmessage(ibgMWUncGeom);
			goto facefound;
		}else{
			rc = ibgDefaultLineWithTriangle(This,
						nint,nface,n1,n2,n3);
		 	return rc;
		}
	}
	dd = dr = 0;
	for(d=0;d<ibgDIM;d++){
	   	dx = ibgpX(nh)[d]-ibgpX(nhf)[d];
		dd += dx>0 ? dx : -dx;
		dr += dx*(ibgpX(nh)[d]-ibgpX(nhb)[d]);
	}
	if(dd>2*splt->g.Delta){
		if(dr>=0){
			for(d=0;d<ibgDIM;d++){
			   	ibgpX(nhb)[d] = (ibgpX(nh)[d]+ibgpX(nhf)[d])/2;
			}
			ibgiRegionOfPoint(splt->gold,&nhb,&nh);
			goto facerecover;
		}else{
			;
		}
	}
facefound:
	nhf = nh;
	fu = ibgpSegment(nhf);
 }
 if(ibgCorrect==ibgtStatusRF(splt->told,splt->rold,fu)){
	rc = SplitRegionLine0(splt,nint,&nhf,&nh1,&nh2,&nh3);
 }else{
	rc = ibgiLineWithTriangle(splt->gold,nint,&nhf,&nh1,&nh2,&nh3);
 }

 if(rc==ibgLineFound){
    if(ibgCorrect==ibgtStatusRL(splt->told,splt->rold,ibgpSegment(*nint))){
    	if((splt->func)(splt->fThis,nint)>0){
  		ibgpSegment(*nint) = splitChangeLineNr(splt,ibgpSegment(*nint));
	}
    }
 }else if(rc==ibgSplitFound){
	ibgpSegment(*nint) = splitCreateLineNr(splt,ibgpSegment(*nint));
	ibgpType(*nint) = ibgSLine;
	return ibgLineFound;
 }else{
    if(ibgCorrect==ibgtStatusRF(splt->told,splt->rold,ibgpSegment(*nint))){
    	if((splt->func)(splt->fThis,nint)>0){
  		ibgpSegment(*nint) = splitChangeFaceNr(splt,ibgpSegment(*nint));
	}
    }
 }
 return  rc;
}

static int Free(ibGeometry This)
{ibgSplitRegion splt=((ibgSplitRegion)This);
 ibgdFree(splt->gold);
 return ibgSuccess;
}
ibGeometryClassRec	ibgSplitRegionClass;

ibGeometry ibgdSplitRegion(ibGeometry 	gold,
			ibgSegment		rnew,
			ibgSegment		rold,
			ibgFloat	(*func)(ibgPtObject fThis, ibgPoint *n),
			ibgPtObject	fThis,
			ibgFloat	Dfunc)
{ibgSplitRegion splt=(ibgSplitRegion) malloc(sizeof(ibgSplitRegionRec));
 ibgdInitialize((ibGeometry)splt,&ibgSplitRegionClass);
 ibgtCopy(&(splt->g.top),&(gold->top));
 splt->gold = gold;
 splt->told = &(gold->top);
 splt->rold = rold; splt->rnew=rnew;
 splt->face = ibgDefaultFaceNr(rold,rnew);
 splt->func = func; splt->fThis=fThis;
 splt->epsilon = 0.5*splt->g.Delta*Dfunc;
 ibgSplitRegionClass.RegionOfPoint	= RegionOfPoint;
 ibgSplitRegionClass.FaceWithEdge	= FaceWithEdge;
 ibgSplitRegionClass.LineWithTriangle	= LineWithTriangle;
 ibgSplitRegionClass.Free		= Free;
 return (ibGeometry) splt;
}