cogrefine.cxx

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

#include "cog.hxx"

cogFloat xon(const wzPoint& x,const wzPoint& y,const wzPoint& z);

// This file defines the standard binary refinement algorithm 
/*
cogIndex Cogeometry::refine(cogFlag1 &f,
			    const cogLine& s) const;
cogIndex Cogeometry::refine(cogFlag2& f, const cogFlag1& o,
			    const cogTriangle& s) const;
cogIndex Cogeometry::refine(cogFlag3& f, const cogFlag2& o,
			    const cogTetrahedron& s) const;
			    */

cogIndex Cogeometry::refine(cogFlag1 &f,
			    const cogLine& s) const
{
  cogPoint m; cogLine(s.c0,s.c1).barycentre(m); Point(m);
  cogIndex rc = Line(f,cogLine(s.c0,m));
  if(rc==cogRCFaceFound) return rc;
  return Line(f,cogLine(m,s.c1));
}

cogIndex Cogeometry::refine(cogFlag2& f, const cogFlag1& o,
			    const cogTriangle& s) const
{
  cogPoint m2,m0,m1;
  cogIndex rc,rt;
  cogLine(s.c1,s.c2).barycentre(m0); Point(m0);
  cogLine(s.c2,s.c0).barycentre(m1); Point(m1);
  cogLine(s.c0,s.c1).barycentre(m2); Point(m2);
  cogFloat sc = m2.scalar(s.c0,o.p1);
  cogIndex firsttry=1;
nexttry:
  try{
    if(sc > 0){
      rc=Triangle(f,o,cogTriangle(s.c0,m2,m1));
      switch(rc){
      case cogRCEdgeFound: 	return rc;
      case cogRCFaceFoundOn01: 	return rc;
      case cogRCFaceFoundOn02:	return rc;
      case cogRCFaceFoundOn12: 	rt = 12; break;
      }
    }else{    
      rc=Triangle(f,o,cogTriangle(m2,s.c1,m0));
      switch(rc){
      case cogRCEdgeFound: 	return rc;
      case cogRCFaceFoundOn01: 	return rc;
      case cogRCFaceFoundOn12:	return rc;
      case cogRCFaceFoundOn02: 	rt = 20; break;
      }
    }
  }catch(cogTriangleRoundTripFailed){
    if(s.c0.near(o.p1,Delta2D)) throw;
    if(s.c1.near(o.p1,Delta2D)) throw;
    if(m2.near(o.p1,Delta2D)){
      if(!firsttry){
	wzAssert(sc==0);
	m2.affine(0.4,s.c0,s.c1);
	sc = m2.scalar(s.c0,o.p1);
	goto nexttry;
      }
      firsttry = 0; sc = -sc;
      goto nexttry;
    }
    throw; //wzAssert(0);
  }
  cogPoint h1,h0,ho; cogFlag1 t(h1,h0,ho); int i;
  i=0;
  for(;i<30;i++){
    t = f;
    switch(rt){
    case 01:
      rc=Triangle(f,t,cogTriangle(m0,m1,m2));
      switch(rc){
      case cogRCEdgeFound: 	return rc;
      case cogRCFaceFoundOn01: 	rt = 10; continue;
      case cogRCFaceFoundOn12: 	rt = 21; continue;
      case cogRCFaceFoundOn02: 	rt = 02; continue;
      }
      break;
    case 10:
      rc=Triangle(f,t,cogTriangle(m1,m0,s.c2));
      switch(rc){
      case cogRCEdgeFound: 	return rc;
      case cogRCFaceFoundOn01: 	rt = 01; continue;
      case cogRCFaceFoundOn12: 	return cogRCFaceFoundOn12;
      case cogRCFaceFoundOn02: 	return cogRCFaceFoundOn02;
      }break;
    case 20:
      rc=Triangle(f,t,cogTriangle(m2,m0,m1));
      switch(rc){
      case cogRCEdgeFound: 	return rc;
      case cogRCFaceFoundOn01: 	rt = 02; continue;
      case cogRCFaceFoundOn12: 	rt = 10; continue;
      case cogRCFaceFoundOn02: 	rt = 21; continue;
      }break;
    case 02:
      rc=Triangle(f,t,cogTriangle(m0,m2,s.c1));
      switch(rc){
      case cogRCEdgeFound: 	return rc;
      case cogRCFaceFoundOn01: 	rt = 20; continue;
      case cogRCFaceFoundOn12: 	return cogRCFaceFoundOn01;
      case cogRCFaceFoundOn02: 	return cogRCFaceFoundOn12;
      }break;
    case 12:
      rc=Triangle(f,t,cogTriangle(m1,m2,m0));
      switch(rc){
      case cogRCEdgeFound: 	return rc;
      case cogRCFaceFoundOn01: 	rt = 21; continue;
      case cogRCFaceFoundOn12: 	rt = 02; continue;
      case cogRCFaceFoundOn02: 	rt = 10; continue;
      }
      break;
    case 21:
      rc=Triangle(f,t,cogTriangle(m2,m1,s.c0));
      switch(rc){
      case cogRCEdgeFound: 	return rc;
      case cogRCFaceFoundOn01: 	rt = 12; continue;
      case cogRCFaceFoundOn12: 	return cogRCFaceFoundOn02;
      case cogRCFaceFoundOn02: 	return cogRCFaceFoundOn01;
      }break;
    default: wzAssert(0);
    }
  }
  throw cogErrorFaceLost();
}

cogIndex Cogeometry::refine(cogFlag3& f, const cogFlag2& o,
			    const cogTetrahedron& s) const
{
  cogPoint m4,m5,m6,m7,m8,m9;
  cogIndex rc,rt,incoming;
  cogFloat vol=s.volume()/2;
  cogLine(s.c1,s.c2).barycentre(m4); Point(m4);
  cogLine(s.c2,s.c0).barycentre(m5); Point(m5);
  cogLine(s.c0,s.c1).barycentre(m6); Point(m6);
  cogLine(s.c0,s.c3).barycentre(m7); Point(m7);
  cogLine(s.c1,s.c3).barycentre(m8); Point(m8);
  cogLine(s.c2,s.c3).barycentre(m9); Point(m9);
  cogTetrahedron tb0(o.p2,s.c1,s.c2,s.c3); if(tb0.volume()>vol){
    incoming = 0;
  }else{
    cogTetrahedron tb1(s.c0,o.p2,s.c2,s.c3); if(tb1.volume()>vol){
      incoming = 1;
    }else{
      cogTetrahedron tb2(s.c0,s.c1,o.p2,s.c3); if(tb2.volume()>vol){
	incoming = 2;
      }else{
	incoming = 3;
      }
    }
  }
  cogTetrahedron tb3(s.c0,s.c1,s.c2,o.p2);
  wzAssert(tb3.volume()<0.01*vol);
  wzAssert(tb3.volume()>-0.01*vol);
  switch(incoming){
  case 0:
    rc=Tetrahedron(f,o,cogTetrahedron(s.c0,m6,m5,m7));
    switch(rc){
    case cogRCVertexFound: 	return rc;
    case cogRCEdgeFoundOn012: 	return rc;
    case cogRCEdgeFoundOn013:	return rc;
    case cogRCEdgeFoundOn023:	return rc;
    case cogRCEdgeFoundOn123: 	rt = 576; break;
    }break;
  case 1:
    rc=Tetrahedron(f,o,cogTetrahedron(m6,s.c1,m4,m8));
    switch(rc){
    case cogRCVertexFound: 	return rc;
    case cogRCEdgeFoundOn012: 	return rc;
    case cogRCEdgeFoundOn013:	return rc;
    case cogRCEdgeFoundOn123:	return rc;
    case cogRCEdgeFoundOn023: 	rt = 468; break;
    }break;
  case 2:
    rc=Tetrahedron(f,o,cogTetrahedron(m5,m4,s.c2,m9));
    switch(rc){
    case cogRCVertexFound: 	return rc;
    case cogRCEdgeFoundOn012: 	return rc;
    case cogRCEdgeFoundOn023:	return rc;
    case cogRCEdgeFoundOn123:	return rc;
    case cogRCEdgeFoundOn013: 	rt = 495; break;
    }break;
  case 3:
    rc=Tetrahedron(f,o,cogTetrahedron(m6,m4,m5,m8));
    switch(rc){
    case cogRCVertexFound: 	return rc;
    case cogRCEdgeFoundOn012: 	return rc;
    case cogRCEdgeFoundOn013:	rt = 486; break;
    case cogRCEdgeFoundOn023:	rt = 568; break;
    case cogRCEdgeFoundOn123: 	rt = 458; break;
    }break;
  }
  cogPoint h2,h1,h0,ho; cogFlag2 t(h2,h1,h0,ho);
  for(int i=0;i<50;i++){
    t = f;
    switch(rt){
    case 458:
      rc=Tetrahedron(f,t,cogTetrahedron(m4,m5,m8,m9));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 485; continue;
      case cogRCEdgeFoundOn013:	rt = 459; continue;
      case cogRCEdgeFoundOn023:	return cogRCEdgeFoundOn123;
      case cogRCEdgeFoundOn123:	rt = 589; continue;
      }break;
    case 459:
      rc=Tetrahedron(f,t,cogTetrahedron(m4,m5,m9,s.c2));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 495; continue;
      case cogRCEdgeFoundOn013:	return cogRCEdgeFoundOn012;
      case cogRCEdgeFoundOn023:	return cogRCEdgeFoundOn123;
      case cogRCEdgeFoundOn123:	return cogRCEdgeFoundOn023;
      }break;
    case 468:
      rc=Tetrahedron(f,t,cogTetrahedron(m4,m6,m8,m5));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 486; continue;
      case cogRCEdgeFoundOn013:	return cogRCEdgeFoundOn012;
      case cogRCEdgeFoundOn023:	rt = 458; continue;
      case cogRCEdgeFoundOn123:	rt = 568; continue;
      }break;
    case 485:
      rc=Tetrahedron(f,t,cogTetrahedron(m4,m8,m5,m6));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 458; continue;
      case cogRCEdgeFoundOn013:	rt = 486; continue;
      case cogRCEdgeFoundOn023:	return cogRCEdgeFoundOn012;
      case cogRCEdgeFoundOn123:	rt = 568; continue;
      }break;
    case 486:
      rc=Tetrahedron(f,t,cogTetrahedron(m4,m8,m6,s.c1));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 468; continue;
      case cogRCEdgeFoundOn013:	return cogRCEdgeFoundOn123;
      case cogRCEdgeFoundOn023:	return cogRCEdgeFoundOn012;
      case cogRCEdgeFoundOn123:	return cogRCEdgeFoundOn013;
      }break;
    case 495:
      rc=Tetrahedron(f,t,cogTetrahedron(m4,m9,m5,m8));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 459; continue;
      case cogRCEdgeFoundOn013:	return cogRCEdgeFoundOn123;
      case cogRCEdgeFoundOn023:	rt = 485; continue;
      case cogRCEdgeFoundOn123:	rt = 589; continue;
      }break;
    case 567:
      rc=Tetrahedron(f,t,cogTetrahedron(m5,m6,m7,s.c0));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 576; continue;
      case cogRCEdgeFoundOn013:	return cogRCEdgeFoundOn012;
      case cogRCEdgeFoundOn023:	return cogRCEdgeFoundOn023;
      case cogRCEdgeFoundOn123:	return cogRCEdgeFoundOn013;
      }break;
    case 568:
      rc=Tetrahedron(f,t,cogTetrahedron(m5,m6,m8,m7));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 586; continue;
      case cogRCEdgeFoundOn013:	rt = 567; continue;
      case cogRCEdgeFoundOn023:	rt = 578; continue;
      case cogRCEdgeFoundOn123:	return cogRCEdgeFoundOn013;
      }break;
    case 576:
      rc=Tetrahedron(f,t,cogTetrahedron(m5,m7,m6,m8));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 567; continue;
      case cogRCEdgeFoundOn013:	rt = 578; continue;
      case cogRCEdgeFoundOn023:	rt = 586; continue;
      case cogRCEdgeFoundOn123:	return cogRCEdgeFoundOn013;
      }break;
    case 578:
      rc=Tetrahedron(f,t,cogTetrahedron(m5,m7,m8,m9));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 587; continue;
      case cogRCEdgeFoundOn013:	return cogRCEdgeFoundOn023;
      case cogRCEdgeFoundOn023:	rt = 598; continue;
      case cogRCEdgeFoundOn123:	rt = 789; continue;
      }break;
    case 586:
      rc=Tetrahedron(f,t,cogTetrahedron(m5,m8,m6,m4));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 568; continue;
      case cogRCEdgeFoundOn013:	rt = 458; continue;
      case cogRCEdgeFoundOn023:	return cogRCEdgeFoundOn012;
      case cogRCEdgeFoundOn123:	rt = 486; continue;
      }break;
    case 587:
      rc=Tetrahedron(f,t,cogTetrahedron(m5,m8,m7,m6));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 578; continue;
      case cogRCEdgeFoundOn013:	rt = 586; continue;
      case cogRCEdgeFoundOn023:	rt = 567; continue;
      case cogRCEdgeFoundOn123:	return cogRCEdgeFoundOn013;
      }break;
    case 589:
      rc=Tetrahedron(f,t,cogTetrahedron(m5,m8,m9,m7));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 598; continue;
      case cogRCEdgeFoundOn013:	rt = 587; continue;
      case cogRCEdgeFoundOn023:	return cogRCEdgeFoundOn023;
      case cogRCEdgeFoundOn123:	rt = 789; continue;
      }break;
    case 598:
      rc=Tetrahedron(f,t,cogTetrahedron(m5,m9,m8,m4));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 589; continue;
      case cogRCEdgeFoundOn013:	rt = 459; continue;
      case cogRCEdgeFoundOn023:	rt = 485; continue;
      case cogRCEdgeFoundOn123:	return cogRCEdgeFoundOn123;
      }break;
    case 789:
      rc=Tetrahedron(f,t,cogTetrahedron(m7,m8,m9,s.c3));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 879; continue;
      case cogRCEdgeFoundOn013:	return cogRCEdgeFoundOn013;
      case cogRCEdgeFoundOn023:	return cogRCEdgeFoundOn023;
      case cogRCEdgeFoundOn123:	return cogRCEdgeFoundOn123;
      }break;
    case 879:
      rc=Tetrahedron(f,t,cogTetrahedron(m8,m7,m9,m5));
      switch(rc){
      case cogRCVertexFound: 	return rc;
      case cogRCEdgeFoundOn012:	rt = 789; continue;
      case cogRCEdgeFoundOn013:	rt = 587; continue;
      case cogRCEdgeFoundOn023:	rt = 598; continue;
      case cogRCEdgeFoundOn123:	return cogRCEdgeFoundOn023;
      }break;
    default:
      wzAssert(0);
    }
  }
  throw cogErrorEdgeLost();
}