netgen/libsrc/meshing/netrule2.cpp

#include <mystdlib.h>
#include "meshing.hpp"

namespace netgen
{

netrule :: netrule ()
{
  name = new char[1];
  name[0] = char(0);
  quality = 0;
}

netrule ::  ~netrule()
{
  delete [] name;
  for(int i = 0; i < oldutofreearea_i.Size(); i++)
    delete oldutofreearea_i[i];
  for(int i = 0; i < freezone_i.Size(); i++)
    delete freezone_i[i];
}



void netrule :: SetFreeZoneTransformation (const Vector & devp, int tolclass)
{
  double lam1 = 1.0/tolclass;
  double lam2 = 1.-lam1;

  double mem1[100], mem2[100], mem3[100];

  int vs = oldutofreearea.Height();
  FlatVector devfree(vs, mem1);

  int fzs = freezone.Size();
  transfreezone.SetSize (fzs);

  if (tolclass <= oldutofreearea_i.Size())
    {
      oldutofreearea_i[tolclass-1] -> Mult (devp, devfree);

      auto& fzi = *freezone_i[tolclass-1];
      for (int i = 0; i < fzs; i++)
	{
	  transfreezone[i][0] = fzi[i][0] + devfree[2*i];
	  transfreezone[i][1] = fzi[i][1] + devfree[2*i+1];
	}
    }
  else
    {
      FlatVector devfree1(vs, mem2);
      FlatVector devfree2(vs, mem3);

      oldutofreearea.Mult (devp, devfree1);
      oldutofreearealimit.Mult (devp, devfree2);
      devfree.Set2 (lam1, devfree1, lam2, devfree2);

      for (int i = 0; i < fzs; i++)
	{
	  transfreezone[i][0] = lam1 * freezone[i][0] + lam2 * freezonelimit[i][0] + devfree[2*i];
	  transfreezone[i][1] = lam1 * freezone[i][1] + lam2 * freezonelimit[i][1] + devfree[2*i+1];
	}
    }


  if (fzs > 0)
    {
      fzmaxx = fzminx = transfreezone[0][0];
      fzmaxy = fzminy = transfreezone[0][1];
    }

  for (int i = 1; i < fzs; i++)
    {
      if (transfreezone[i][0] > fzmaxx) fzmaxx = transfreezone[i][0];
      if (transfreezone[i][0] < fzminx) fzminx = transfreezone[i][0];
      if (transfreezone[i][1] > fzmaxy) fzmaxy = transfreezone[i][1];
      if (transfreezone[i][1] < fzminy) fzminy = transfreezone[i][1];
    }

  for (int i = 0; i < fzs; i++)
    {
      Point2d p1 = transfreezone[i];
      Point2d p2 = transfreezone[(i+1) % fzs];

      Vec2d vn (p2.Y() - p1.Y(), p1.X() - p2.X());

      double len2 = vn.Length2();

      if (len2 < 1e-10)
	{
	  freesetinequ(i, 0) = 0;
	  freesetinequ(i, 1) = 0;
	  freesetinequ(i, 2) = -1;
	}
      else
	{
	  vn /= sqrt (len2);    // scaling necessary ?

	  freesetinequ(i,0) = vn.X(); 
	  freesetinequ(i,1) = vn.Y(); 
	  freesetinequ(i,2) = -(p1.X() * vn.X() + p1.Y() * vn.Y());
	}
    }
}


/*
int netrule :: IsInFreeZone2 (const Point2d & p) const
{
  for (int i = 0; i < transfreezone.Size(); i++)
    {
      if (freesetinequ(i, 0) * p.X() + 
	  freesetinequ(i, 1) * p.Y() +
	  freesetinequ(i, 2) > 0) return 0;
    }
  return 1;
}
*/

int netrule :: IsLineInFreeZone2 (const Point2d & p1, const Point2d & p2) const
{
  if ( (p1.X() > fzmaxx && p2.X() > fzmaxx) ||
       (p1.X() < fzminx && p2.X() < fzminx) ||
       (p1.Y() > fzmaxy && p2.Y() > fzmaxy) ||
       (p1.Y() < fzminy && p2.Y() < fzminy) ) return 0;

  for (int i = 1; i <= transfreezone.Size(); i++)
    {
      if (freesetinequ.Get(i, 1) * p1.X() + freesetinequ.Get(i, 2) * p1.Y() +
	  freesetinequ.Get(i, 3) > -1e-8 &&    // -1e-6
	  freesetinequ.Get(i, 1) * p2.X() + freesetinequ.Get(i, 2) * p2.Y() +
	  freesetinequ.Get(i, 3) > -1e-8       // -1e-6
	  ) return 0;
    }

  double nx =  (p2.Y() - p1.Y());
  double ny = -(p2.X() - p1.X());
  double nl = sqrt (nx * nx + ny * ny);
  if (nl > 1e-8)
    {
      nx /= nl;
      ny /= nl;
      double c = - (p1.X() * nx + p1.Y() * ny);

      bool allleft = true;
      bool allright = true;

      for (int i = 1; i <= transfreezone.Size(); i++)
	{
	  bool left  = transfreezone.Get(i)[0] * nx + transfreezone.Get(i)[1] * ny + c <  1e-7;
          bool right = transfreezone.Get(i)[0] * nx + transfreezone.Get(i)[1] * ny + c > -1e-7;
	  if (!left) allleft = false;
	  if (!right) allright = false;
	}
      if (allleft || allright) return false;
    }

  return true;
}

int netrule :: ConvexFreeZone () const
{
  int n = transfreezone.Size();
  for (int i = 1; i <= n; i++)
    {
      const bool counterclockwise = CCW (transfreezone.Get(i), 
					 transfreezone.Get(i % n + 1),
					 transfreezone.Get( (i+1) % n + 1 ),
					 1e-7);
      //(*testout) << "ccw " << counterclockwise << endl << " p1 " << transfreezone.Get(i) << " p2 " << transfreezone.Get(i % n + 1)
      //		 << " p3 " << transfreezone.Get( (i+1) % n + 1 ) << endl;
      if (!counterclockwise )
	return 0;
    }
  return 1;
}


/*
float netrule :: CalcPointDist (int pi, const Point2d & p) const
{
  float dx = p.X() - points.Get(pi).X();
  float dy = p.Y() - points.Get(pi).Y();
  const threefloat * tf = &tolerances.Get(pi);

  return tf->f1 * dx * dx + tf->f2 * dx * dy + tf->f3 * dy * dy;
}
*/

float netrule :: CalcLineError (int li, const Vec<2> & v) const
{
  float dx = v[0] - linevecs.Get(li)[0];
  float dy = v[1] - linevecs.Get(li)[1];

  const threefloat * ltf = &linetolerances.Get(li);
  return ltf->f1 * dx * dx + ltf->f2 * dx * dy + ltf->f3 * dy * dy;
}
} // namespace netgen