cogdefault.cxx

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

#include "cog.hxx"

// This file contains the default implementations: 
// - search for boundary continuation;
// - if below Delta the barycentre;
// - else refinement;
/*
cogIndex Cogeometry::Line(cogFlag1 &f,
			  const cogLine& s) const;
cogIndex Cogeometry::Triangle(cogFlag2& f, const cogFlag1& o,
			      const cogTriangle& s) const;
cogIndex Cogeometry::Tetrahedron(cogFlag3& f, const cogFlag2& o,
                                 const cogTetrahedron& s) const;
				 */

cogIndex Cogeometry::Point(cogPoint& p) const
{
  p.segment() = DefaultRegion;
  return cogRCRegionFound;
}

cogIndex Cogeometry::BoundaryCondition(cogFlag1& f) const
{
  f.p1.segment() = DefaultFace;
  return cogRCConditionFound;
}

cogIndex Cogeometry::Line(cogFlag1& f,
			  const cogLine& s) const
{
  if(s.c1.segment() == s.c0.segment()){
    f.p0 = s.c1; return cogRCFaceNotFound;
  }
  if(s.diameter()<Delta){
    f.p0 = s.c0; f.po = s.c1;
    s.barycentre(f.p1);
    BoundaryCondition(f);
    return cogRCFaceFound;
  }
  return refine(f,s);
}

cogIndex Cogeometry::Triangle(cogFlag2& f, const cogFlag1& o,
			      const cogTriangle& s) const
{
  cogIndex rc,rt;
  cogFlag1 r(o.p1,o.po,o.p0);
  if(   Line(f,r,cogLine(s.c1,s.c0))==cogRCFaceFound)	{rc=cogRCFaceFoundOn01;}
  else if(Line(f,cogLine(s.c0,s.c2))==cogRCFaceFound)	{rc=cogRCFaceFoundOn02;}
  else if(Line(f,cogLine(s.c2,s.c1))==cogRCFaceFound)	{rc=cogRCFaceFoundOn12;}
  else if(Line(f,cogLine(s.c1,o.po))==cogRCFaceFound)	{rc=cogRCFaceFoundOn01;}
  else wzAssert(0);
  cogPoint h1,h0,ho; cogFlag1 t(h1,h0,ho); 
  if(   Line(t,o,cogLine(s.c0,s.c1))==cogRCFaceFound)	{rt=cogRCFaceFoundOn01;}
  else if(Line(t,cogLine(s.c1,s.c2))==cogRCFaceFound)	{rt=cogRCFaceFoundOn12;}
  else if(Line(t,cogLine(s.c2,s.c0))==cogRCFaceFound)	{rt=cogRCFaceFoundOn02;}
  else if(Line(t,cogLine(s.c0,o.p0))==cogRCFaceFound)	{rt=cogRCFaceFoundOn01;}
  else wzAssert(0);
  if(rc == cogRCFaceFoundOn01){
    if(rt != cogRCFaceFoundOn01){otherSide(f,t); return rt;}
    if(f.p1.near(o.p1,Delta)){
      if(! t.p1.near(o.p1,Delta)){otherSide(f,t); return rt;}
      if(f.p1.segment() == o.p1.segment()){
	if(o.p1.near(s.c0,Delta) ||o.p1.near(s.c1,Delta))
	  throw cogTriangleRoundTripFailed();
      }
    }
  }
  if(rt == cogRCFaceFoundOn01){
    if(rc != cogRCFaceFoundOn01){return rc;}
    if(t.p1.near(o.p1,Delta)){
      if(! f.p1.near(o.p1,Delta)){return rc;}
      if(t.p1.segment() == o.p1.segment()){
	if(o.p1.near(s.c0,Delta) ||o.p1.near(s.c1,Delta))
	  throw cogTriangleRoundTripFailed();
      }
    }
  }
  if(f.p0.segment() != o.p0.segment()){throw cogTriangleRoundTripFailed();}
  if(t.p0.segment() != o.po.segment()){throw cogTriangleRoundTripFailed();}
  if((rc == rt) &&
     (f.po.segment() == o.po.segment()) &&
     (f.p1.segment() == o.p1.segment()) &&
     (f.p1.near(t.p1,Delta))){
    wzAssert(f.po.segment() == t.p0.segment());
    return rc;
  }
  cogFloat diameter=s.diameter();
  if(diameter<singleIntersectionLimit2D){
    if((f.po.segment() == o.po.segment()) &&
       (f.p1.segment() == o.p1.segment())){
      return rc;
    }
    if((t.po.segment() == o.p0.segment()) &&
       (t.p1.segment() == o.p1.segment())){
      otherSide(f,t);
      return rt;
    }
  }
  if(diameter<Delta2D){
    if(!specialRefinement){
      cogFloat oldDelta = Delta;
      ((Cogeometry*)this)->specialRefinement = 1;
      ((Cogeometry*)this)->setDelta(0.1*Delta*Delta/Delta2D);
      cogPoint h1,h0,ho; cogFlag1 h(h1,h0,ho); Line(h,cogLine(o.p0,o.po));
      if(o.p1.segment()!=h1.segment()){
	Line(f,cogLine(o.po,o.p0)); 
	otherSide(h,f);
	if(o.p1.segment()!=h1.segment()){	// dirty trick
	  h1.segment() = o.p1.segment();
	}
      }
      try{
	rc = refine(f,h,s);
	if(rc==cogRCEdgeFound){
	  f = h;
	}
      }catch(...){
	((Cogeometry*)this)->setDelta(oldDelta);
	((Cogeometry*)this)->specialRefinement = 0;
	throw;
      }
      ((Cogeometry*)this)->setDelta(oldDelta);
      ((Cogeometry*)this)->specialRefinement = 0;
      return rc;
    }
    wzAssert(specialRefinement);
    {
      int onLine;
      wzPoint mm; s.barycentre(mm); Point(mm);
      if(     Line(f,  cogLine(s.c0,  mm))==cogRCFaceFound) {onLine=0;}
      else if(Line(f,  cogLine(  mm,s.c1))==cogRCFaceFound) {onLine=0;}
      else if(Line(f,o,cogLine(s.c0,s.c1))==cogRCFaceFound) {onLine=1;}
      else if(Line(f,r,cogLine(s.c1,s.c0))==cogRCFaceFound) {onLine=2;}
      else throw cogTriangleRoundTripFailed();
      if(onLine){
	wzPoint ml; cogLine(f.p1,o.p1).barycentre(ml); Point(ml);
	if(Line(f,cogLine(ml,mm))!=cogRCFaceFound){
	  if(onLine==1) {
	    if(Line(f,cogLine(ml,s.c1))!=cogRCFaceFound){
	      //	      wzDebugger();
	      throw cogTriangleRoundTripFailed();
	    }
	  }else{
	    if(Line(f,cogLine(ml,s.c1))!=cogRCFaceFound){
	      //	      wzDebugger();
	      throw cogTriangleRoundTripFailed();
	    }
	  }
	}
      }
    }
    f.p2 = f.p1;
    f.p2.segment() = cogEdge();
    f = o;
    wzAssert(f.p0.segment()==o.p0.segment());
    wzAssert(f.po.segment()==o.po.segment());
    return cogRCEdgeFound;
  }
  rc = refine(f,o,s);
  if(rc == cogRCFaceFoundOn01){
    wzAssert(0);
    rc = refine(f,o,s);
    if(f.p1.near(o.p1,Delta)){
      ;
    }
  }
  return rc;
}

cogIndex   Cogeometry::Tetrahedron(cogFlag3& f, const cogFlag2& o,
                                 const cogTetrahedron& s) const
{
  cogIndex rc,rt,i;
  cogPoint h1,ho,h0; cogFlag1 t(h1,h0,ho);
  // test of correct input state
  cogFloat vol=1.01*s.volume(); wzAssert(vol>0);
  cogFloat v,v0 = -0.01*vol;
  cogTetrahedron tb0(o.p2,s.c1,s.c2,s.c3); wzAssert((v=tb0.volume())<vol && (v>v0));
  cogTetrahedron tb1(s.c0,o.p2,s.c2,s.c3); wzAssert((v=tb1.volume())<vol && (v>v0));
  cogTetrahedron tb2(s.c0,s.c1,o.p2,s.c3); wzAssert((v=tb2.volume())<vol && (v>v0));
  cogTetrahedron tb3(s.c0,s.c1,s.c2,o.p2); wzAssert((v=tb3.volume())<-v0 && (v>v0));
  // real start of the program;

  try{
  rc = Triangle(f,o,cogTriangle(s.c0,s.c1,s.c2));
  switch(rc){
  case cogRCEdgeFound:		{rc=cogRCEdgeFoundOn012; goto edgefound;}
  case cogRCFaceFoundOn01:	rt = 10; break;
  case cogRCFaceFoundOn12:	rt = 21; break;
  case cogRCFaceFoundOn02:	rt = 02; break;
  }
  otherSide(t,f);
  //t=f;
  for(i=0;i<15;i++){
    switch(rt){
    case 01:
      rc = Triangle(f,t,cogTriangle(s.c0,s.c1,s.c2));
      switch(rc){
      case cogRCEdgeFound:	{rc=cogRCEdgeFoundOn012; goto edgefound;}
      case cogRCFaceFoundOn01:	rt = 10; break;
      case cogRCFaceFoundOn12:	rt = 21; break;
      case cogRCFaceFoundOn02:	rt = 02; break;
      } break;
    case 12:
      rc = Triangle(f,t,cogTriangle(s.c1,s.c2,s.c0));
      switch(rc){
      case cogRCEdgeFound:	{rc=cogRCEdgeFoundOn012; goto edgefound;}
      case cogRCFaceFoundOn01:	rt = 21; break;
      case cogRCFaceFoundOn12:	rt = 02; break;
      case cogRCFaceFoundOn02:	rt = 10; break;
      } break;
    case 20:
      rc = Triangle(f,t,cogTriangle(s.c2,s.c0,s.c1));
      switch(rc){
      case cogRCEdgeFound:	{rc=cogRCEdgeFoundOn012; goto edgefound;}
      case cogRCFaceFoundOn01:	rt = 02; break;
      case cogRCFaceFoundOn12:	rt = 10; break;
      case cogRCFaceFoundOn02:	rt = 21; break;
      } break;
    case 10:
      rc = Triangle(f,t,cogTriangle(s.c1,s.c0,s.c3));
      switch(rc){
      case cogRCEdgeFound:	{rc=cogRCEdgeFoundOn013; goto edgefound;}
      case cogRCFaceFoundOn01:	rt = 01; break;
      case cogRCFaceFoundOn12:	rt = 30; break;
      case cogRCFaceFoundOn02:	rt = 13; break;
      }break;
    case 03:
      rc = Triangle(f,t,cogTriangle(s.c0,s.c3,s.c1));
      switch(rc){
      case cogRCEdgeFound:	{rc=cogRCEdgeFoundOn013; goto edgefound;}
      case cogRCFaceFoundOn01:	rt = 30; break;
      case cogRCFaceFoundOn12:	rt = 13; break;
      case cogRCFaceFoundOn02:	rt = 01; break;
      }break;
    case 31:
      rc = Triangle(f,t,cogTriangle(s.c3,s.c1,s.c0));
      switch(rc){
      case cogRCEdgeFound:	{rc=cogRCEdgeFoundOn013; goto edgefound;}
      case cogRCFaceFoundOn01:	rt = 13; break;
      case cogRCFaceFoundOn12:	rt = 01; break;
      case cogRCFaceFoundOn02:	rt = 30; break;
      }break;
    case 21:
      rc = Triangle(f,t,cogTriangle(s.c2,s.c1,s.c3));
      switch(rc){
      case cogRCEdgeFound:	{rc=cogRCEdgeFoundOn123; goto edgefound;}
      case cogRCFaceFoundOn01:	rt = 12; break;
      case cogRCFaceFoundOn12:	rt = 31; break;
      case cogRCFaceFoundOn02:	rt = 23; break;
      }break;
    case 13:
      rc = Triangle(f,t,cogTriangle(s.c1,s.c3,s.c2));
      switch(rc){
      case cogRCEdgeFound:	{rc=cogRCEdgeFoundOn123; goto edgefound;}
      case cogRCFaceFoundOn01:	rt = 31; break;
      case cogRCFaceFoundOn12:	rt = 23; break;
      case cogRCFaceFoundOn02:	rt = 12; break;
      }break;
    case 32:
      rc = Triangle(f,t,cogTriangle(s.c3,s.c2,s.c1));
      switch(rc){
      case cogRCEdgeFound:	{rc=cogRCEdgeFoundOn123; goto edgefound;}
      case cogRCFaceFoundOn01:	rt = 23; break;
      case cogRCFaceFoundOn12:	rt = 12; break;
      case cogRCFaceFoundOn02:	rt = 31; break;
      }break;
    case 02:
      rc = Triangle(f,t,cogTriangle(s.c0,s.c2,s.c3));
      switch(rc){
      case cogRCEdgeFound:	{rc=cogRCEdgeFoundOn023; goto edgefound;}
      case cogRCFaceFoundOn01:	rt = 20; break;
      case cogRCFaceFoundOn12:	rt = 32; break;
      case cogRCFaceFoundOn02:	rt = 03; break;
      }break;
    case 23:
      rc = Triangle(f,t,cogTriangle(s.c2,s.c3,s.c0));
      switch(rc){
      case cogRCEdgeFound:	{rc=cogRCEdgeFoundOn023; goto edgefound;}
      case cogRCFaceFoundOn01:	rt = 32; break;
      case cogRCFaceFoundOn12:	rt = 03; break;
      case cogRCFaceFoundOn02:	rt = 20; break;
      }break;
    case 30:
      rc = Triangle(f,t,cogTriangle(s.c3,s.c0,s.c2));
      switch(rc){
      case cogRCEdgeFound:	{rc=cogRCEdgeFoundOn023; goto edgefound;}
      case cogRCFaceFoundOn01:	rt = 03; break;
      case cogRCFaceFoundOn12:	rt = 20; break;
      case cogRCFaceFoundOn02:	rt = 32; break;
      }break;
    }
    t = f;
  }
  wzAssert(0);
  }catch(cogTriangleRoundTripFailed){wzDebugger();}
edgefound:
  if(f.p2.segment() != o.p2.segment()) goto notfound;
  {
    wzPoint ho2,ho1,ho0,hoo; cogFlag2 to(ho2,ho1,ho0,hoo);
    wzPoint hf2,hf1,hf0,hfo; cogFlag2 tf(hf2,hf1,hf0,hfo);
    wzPoint go2,go1,go0,goo; cogFlag2 oo(go2,go1,go0,goo);
    wzPoint gf2,gf1,gf0,gfo; cogFlag2 of(gf2,gf1,gf0,gfo);
    /*	this would be the better code if nextFace is stable enough for iterative calls:
    cogIndex k;
    tf = f; otherSide(to,o);
    for(k=0;k<10;k++){
      if(tf.p0.segment() != to.p0.segment()) goto notfound;
      if(tf.p1.segment() != to.p1.segment()) goto notfound;
      of = tf; oo = to; nextFace(tf,of); nextFace(to,oo);
      if((   tf.p0.segment() == f.p0.segment())
	 && (tf.p1.segment() == f.p1.segment())
	 && (tf.po.segment() == f.po.segment())){
	return rc;
      }
    }	The following code is a replacement: test only the two next neighbour faces:
    */
    
    otherSide(to,f); nextFace(of,to); nextFace(tf,f); 
    otherSide(to,o); nextFace(oo,to); nextFace(to,o);
    if(f.p0.segment() == o.p0.segment()){
      if( f.po.segment() !=  o.po.segment()) goto notfound;
      if(tf.p1.segment() != to.p1.segment()) goto notfound;
      if(of.p1.segment() != oo.p1.segment()) goto notfound;
      if(tf.p0.segment() != to.p0.segment()) goto notfound;
      if(of.p0.segment() != oo.p0.segment()) goto notfound;
      if(tf.po.segment() != to.po.segment()) goto notfound;
      if(of.po.segment() != oo.po.segment()) goto notfound;
    }else if (f.p0.segment() == o.po.segment()){	 // should be checked out later.
      if( f.po.segment() !=  o.p0.segment()) goto notfound;
      if(tf.p1.segment() != oo.p1.segment()) goto notfound;
      if(of.p1.segment() != to.p1.segment()) goto notfound;
      if(tf.p0.segment() != oo.p0.segment()) goto notfound;
      if(of.p0.segment() != to.p0.segment()) goto notfound;
      if(tf.po.segment() != oo.po.segment()) goto notfound;
      if(of.po.segment() != to.po.segment()) goto notfound;
    }else{
      goto notfound;
    }
    return rc;
  }    
notfound:
  if(s.diameter()<Delta3D){
    f = o;
    s.barycentre(f.p3);
    f.p3.segment() = cogVertex();
    return cogRCVertexFound;
  }
  return refine(f,o,s);
}