cogenchart.cxx

Delaunay三角形的网格剖分程序

#include "cogenchart.hxx"

cogChartStackEntry*	cogChartStackEntry::base = 0;

cogenChart	CogenChart::cogenchart() const {return (CogenChart *) this;}

CogenChart::CogenChart(cogenerator gen, wzMap* c)
  :CogenDecorator(gen),
   coordinates(c),
   originalGeometry(original->geometry())
{
  original->endInitialization();
}

wzFloat CogenChart::chi(const wzPoint& p0) const
{
  wzPoint p;
  coordinates->inverse(p.X,p0.X);
  return original->chi(p);
}

static wzFloat ppmin = 1/(wzInfty*wzInfty);

void CogenChart::getMetric(wzMetricData& data, const wzPoint& p0) const
{
  wzIndex i,j,d=coordinates->dimension();
  wzFloat Delta[3][3],pp,c;
  coordinates->dx(p0.X,Delta[0],Delta[1],Delta[2]);
  for(i=0;i<d;i++){
    for(j=0,pp=0;j<d;j++){
      c = Delta[j][i]; pp += c*c;
    }
    if(pp<ppmin) pp=ppmin;
    ((wzDiagonalMetricData*)(&data))->g_ii[i] = pp;
  }
}

void CogenChart::getMetric(wzMetricData& data, const cogFlag1& f) const
{
  wzIndex i,j,d=coordinates->dimension();
  wzFloat Delta[3][3],pp,c;
  coordinates->dx(f.p1.X,Delta[0],Delta[1],Delta[2]);
  for(i=0;i<d;i++){
    for(j=0,pp=0;j<d;j++){
      c = Delta[j][i]; pp += c*c;
    }
    if(pp<ppmin) pp=ppmin;
    ((wzDiagonalMetricData*)(&data))->g_ii[i] = pp;
  }
}

wzFloat CogenChart::g_ii(const wzMetricData& data,
			 const wzPoint& q, wzIndex i) const
{
  return ((wzDiagonalMetricData*)(&data))->g_ii[i];
}
wzFloat CogenChart::g_ii(const wzMetricData& data,
			 const cogFlag1& q, wzIndex i) const
{
  return ((wzDiagonalMetricData*)(&data))->g_ii[i];
}

cogIndex CogenChart::Point(wzPoint& p0) const
{
  wzPoint p;
  wzIndex rc;
  cogChartStackEntry* base = base->base;
  if(base && base->marker == (void*) this){ // *(void**)& coordinates){
    base->base = base->next;
    rc = originalGeometry->Point(*base->point);
    p0.segment() = base->point->segment();
    base->base = base;
  }else{
    coordinates->inverse(p.X,p0.X);
    rc = originalGeometry->Point(p);
    p0.segment() = p.segment();
  }
  return rc;
}

cogIndex CogenChart::Line (cogFlag1& f0, const cogLine& s0) const
{
  wzPoint fp1,fp0,fpo,sp0,sp1;
  cogLine s(sp0,sp1);
  cogFlag1 f(fp1,fp0,fpo);
  wzIndex rc = Cogeometry::Line(f0,s0);
  if(rc == cogRCFaceFound){
    sp0 = f0.p0; coordinates->inverse(sp0.X,f0.p0.X);
    sp1 = f0.po; coordinates->inverse(sp1.X,f0.po.X);
    //  sp0 = s0.c0; coordinates->inverse(sp0.X,s0.c0.X);
    //  sp1 = s0.c1; coordinates->inverse(sp1.X,s0.c1.X);
    rc = original->geometry()->Line(f,s);
    f0.p1 = fp1; coordinates->call(f0.p1.X,fp1.X);
    f0.p0 = fp0; coordinates->call(f0.p0.X,fp0.X);
    f0.po = fpo; coordinates->call(f0.po.X,fpo.X);
    //    f0.p1 = fp1; coordinates->call(f0.p1.X,fp1.X);
    //    f0.p0 = fp0; coordinates->call(f0.p0.X,fp0.X);
    //    f0.po = fpo; coordinates->call(f0.po.X,fpo.X);
  }
  return rc;
}

cogIndex CogenChart::BoundaryCondition(cogFlag1& f) const
{
  wzPoint p1,p0,po; cogFlag1 pp(p1,p0,po);
  coordinates->inverse(p1.X,f.p1.X);
  coordinates->inverse(p0.X,f.p0.X);
  coordinates->inverse(po.X,f.po.X);
  original->geometry()->BoundaryCondition(pp);
  f.p1.segment() = p1.segment();
  return cogRCConditionFound;
}

CogeometryChart::CogeometryChart(cogeometry old, wzMap* c, void* mark)
  :original(old),
   coordinates(c),
   marker(mark)
{}

cogIndex CogeometryChart::Point(wzPoint& p0) const
{
  wzIndex rc; wzPoint p;
  coordinates->call(p.X,p0.X);
  cogChartStackEntry e;
  e.next = e.base;
  e.point = &p0;
  e.marker = marker; //*(void**)& coordinates;
  e.base = &e;
  try{
    rc = original->Point(p);
  }catch(...){e.base = e.next; throw;}
  e.base = e.next;
  p0.segment() = p.segment();
  return rc;
}

cogIndex CogeometryChart::BoundaryCondition(cogFlag1& f) const
{
  wzPoint p1,p0,po; cogFlag1 pp(p1,p0,po);
  coordinates->call(p1.X,f.p1.X);
  coordinates->call(p0.X,f.p0.X);
  coordinates->call(po.X,f.po.X);
  original->BoundaryCondition(pp);
  f.p1.segment() = p1.segment();
  return cogRCConditionFound;
}

cogIndex CogeometryChart::Line (cogFlag1& f0, const cogLine& s0) const
{
  wzPoint fp1,fp0,fpo,sp0,sp1;
  cogLine s(sp0,sp1);
  cogFlag1 f(fp1,fp0,fpo);
  wzIndex rc = Cogeometry::Line(f0,s0);
  /*
    Dieser Teil ist leider nicht mit dem Stackmechanismus kompatibel.
    Muss nochmal durchdacht werden.

  if(rc == cogRCFaceFound){
    sp0 = f0.p0; coordinates->call(sp0.X,f0.p0.X);
    sp1 = f0.po; coordinates->call(sp1.X,f0.po.X);
    rc = original->Line(f,s);
    f0.p1 = fp1; coordinates->inverse(f0.p1.X,fp1.X);
    f0.p0 = fp0; coordinates->inverse(f0.p0.X,fp0.X);
    f0.po = fpo; coordinates->inverse(f0.po.X,fpo.X);
  } 
  */ 
  return rc;
}