femgraphicdeviceindependent.cpp

FreeFEM is an implementation of the GFEM language dedicated to the finite element method. It provid

// Emacs will be in -*- Mode: c++ -*-
//
// ********** DO NOT REMOVE THIS BANNER **********
//
// SUMMARY: Language for a Finite Element Method
//
// AUTHORS:  C. Prud'homme
// ORG    :          
// E-MAIL :  [e-mail]
//
// ORIG-DATE:     June-94
// LAST-MOD:  6-Apr-02 at 13:10:02 by Christophe Prud'homme
//
// DESCRIPTION:
//   This program is free software; you can redistribute it and/or modify
//   it under the terms of the GNU General Public License as published by
//   the Free Software Foundation; either version 2 of the License, or
//   (at your option) any later version.
  
//   This program is distributed in the hope that it will be useful,
//   but WITHOUT ANY WARRANTY; without even the implied warranty of
//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//   GNU General Public License for more details.
  
//   You should have received a copy of the GNU General Public License
//   along with this program; if not, write to the Free Software
//   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
// DESCRIP-END.
//
#ifdef __GNUG__
#pragma implementation
#endif

/*
 * FreeFem includes
 */
#include <femGraphic.hpp>
#include <femMisc.hpp>

#include <femGraphicDeviceIndependent.hpp>

namespace fem
{
DECLARE_CLASS(femMesh);

femGraphicDeviceIndependent::femGraphicDeviceIndependent( femMeshPtr t )
  :
  __t( t )
{
}
femGraphicDeviceIndependent::~femGraphicDeviceIndependent()
{
}

void
femGraphicDeviceIndependent::setTriangulation( femMeshPtr t )
{
  __t = t;
}


void 
femGraphicDeviceIndependent::Init (femPoint* rp, int ns, char *s)
{
  int             i;
  float           xmn = rp[0][0];
  float           xmx = rp[0][0];
  float           ymn = rp[0][1];
  float           ymx = rp[0][1];
  float           xctr, yctr, ray;
  float           x, y;

  for (i = 1; i < ns; i++)
     {
       x = (rp[i])[0];
       y = (rp[i])[1];
       xmx = x > xmx ? x : xmx;
       xmn = xmn > x ? x : xmn;
       ymx = y > ymx ? y : ymx;
       ymn = ymn > y ? y : ymn;
     }
  xctr = (xmx + xmn) / 2;
  yctr = (ymx + ymn) / 2;
  ray = (xmx - xctr) > (ymx - yctr) ? xmx - xctr : ymx - yctr;
  reffecran ();
  if (*s == 'o')
    cadreortho (xctr, yctr, ray);
  else
    cadre (xmn, xmx, ymn, ymx);

}


void 
femGraphicDeviceIndependent::quicksort (float *tb, int *o, int n)
{
  int             i, j;
  int             b;
  float           x, y;

  while (n > 1)
     {
       x = tb[n / 2];
       for (i = 0, j = n - 1; i <= j; i++, j--)
	  {
	    while (tb[i] > x)
	      i++;
	    while (tb[j] < x)
	      j--;
	    if (i > j)
	      break;
	    y = tb[i];
	    tb[i] = tb[j];
	    tb[j] = y;
	    b = o[i];
	    o[i] = o[j];
	    o[j] = b;
	  }
       n -= i;
       if (n > j + 1)
	  {
	    quicksort (tb, o, j + 1);
	    tb += i;
	    o += i;
	  }
       else
	  {
	    quicksort (tb + i, o + i, n);
	    n = j + 1;
	  }
     }
}


void 
femGraphicDeviceIndependent::initt (void)
{
  femPoint         *r = __t->rp;
  femTriangle       *tt = __t->tr;
  int             n = __t->getNumberOfCells();
  int             i;

  table = new float[n];

//    closegraphique();
  ordre = new int[n];

  for (i = 0; i < n; i++)
    ordre[i] = i;
  for (i = 0; i < n; i++)
    table[i] = (r[tt[i][0]][1] + r[tt[i][1]][1] + r[tt[i][2]][1]) / 3;
}

void 
femGraphicDeviceIndependent::projection (float *f)
{
  int             i;
  femPoint         *rp = __t->rp;

  proj = new femPoint[__t->getNumberOfPoints()];
  for (i = 0; i < __t->getNumberOfPoints(); i++)
     {
       proj[i][0] = rp[i][0] / 5;
       proj[i][1] = f[i] + rp[i][1] / 3;
     }
}

void 
femGraphicDeviceIndependent::graph3d (float *fct, int waitm)
{

  int             i, j, coul;
  int             coul1 = 11, coul2 = 3;
  int             p1, p2, p3;
  float           test;
  femTriangle       *tt = __t->tr;
  femPoint         *pp = __t->rp;

  initt ();
  quicksort (table, ordre, __t->getNumberOfCells());
  delete[]table;
  table = NULL;
  projection (fct);
  Init (proj, __t->getNumberOfPoints(), "n");
  for (i = 0; i < __t->getNumberOfCells(); i++)
     {
       float           poly[6];
       int             ii = ordre[i];

       for (j = 0; j < 3; j++)
	  {
	    poly[2 * j] = proj[tt[ii][j]][0];
	    poly[2 * j + 1] = proj[tt[ii][j]][1];
	  }
       if (pp[tt[ii][0]][1] < pp[tt[ii][1]][1])
	  {
	    p1 = tt[ii][0];
	    p2 = tt[ii][1];
	  }
       else
	  {
	    p2 = tt[ii][0];
	    p1 = tt[ii][1];
	  }
       if (pp[tt[ii][2]][1] < pp[p2][1])
	 if (pp[tt[ii][2]][1] < pp[p1][1])
	    {
	      p3 = p2;
	      p2 = p1;
	      p1 = tt[ii][2];
	    }
	 else
	    {
	      p3 = p2;
	      p2 = tt[ii][2];
	    }
       else
	 p3 = tt[ii][2];
       test = proj[p3][1] * (proj[p1][0] - proj[p2][0]) +
	 proj[p3][0] * (proj[p2][1] - proj[p1][1]) +
	 proj[p2][0] * proj[p1][1] - proj[p1][0] * proj[p2][1];
       if (test <= 0)
	 if (proj[p2][0] < proj[p1][0])
	   coul = coul2;
	 else
	   coul = coul1;
       else if (proj[p2][0] > proj[p1][0])
	 coul = coul2;
       else
	 coul = coul1;
       if (i == 0 && coul == coul2)
	  {
	    int             exchg;

	    exchg = coul2;
	    coul = coul2 = coul1;
	    coul1 = exchg;
	  }
/*     couleurint(coul); */
       /*    raffpoly(3,poly); */
/*     couleur(4); */
       rmoveto (proj[tt[ii][2]][0], proj[tt[ii][2]][1]);
       for (j = 0; j < 3; j++)
	  {
	    rlineto (proj[tt[ii][j]][0], proj[tt[ii][j]][1]);
	  }
     }
  rattente (waitm);
  delete[]proj;
  proj = NULL;
  delete[]ordre;
  ordre = NULL;
}


void 
femGraphicDeviceIndependent::showtriangulation (int waitm)
{
  int             i, j;
  femTriangle       *tt = __t->tr;
  femPoint         *pp = __t->rp;

  Init (pp, __t->getNumberOfPoints(), "o");
  couleur (6);
  for (i = 0; i < __t->getNumberOfCells(); i++)
     {
       rmoveto (pp[tt[i][2]][0], pp[tt[i][2]][1]);
       for (j = 0; j < 3; j++)
	 rlineto (pp[tt[i][j]][0], pp[tt[i][j]][1]);
     }
  rattente (waitm);
}


void 
femGraphicDeviceIndependent::contour (int ng[], int coul)
{
  int             i = 0, j, jp;
  int             tj, tjp;
  femTriangle       *tt = __t->tr;
  femPoint         *pp = __t->rp;

  couleur (coul);
  for (; i < __t->getNumberOfCells(); i++)
    for (j = 0; j < 3; j++)
       {
	 jp = j + 1 == 3 ? 0 : j + 1;
	 tj = tt[i][j];
	 tjp = tt[i][jp];
	 if (ng[tj] != 0 && ng[tjp] != 0)
	    {
	      rmoveto (pp[tj][0], pp[tj][1]);
	      rlineto (pp[tjp][0], pp[tjp][1]);
	    }
       }
}


void 
femGraphicDeviceIndependent::equpot (int *ng, float *f, int nl, int waitm)
{
  float           qp[2][5];
  int             ns = __t->getNumberOfPoints(), nt = __t->getNumberOfCells();
  femTriangle       *tr = __t->tr;
  femPoint         *rp = __t->rp;
  float           xln, fm, xf, xfm, fi, fj, xlam;
  int             im, i, j, k, l, ik, jk;

  Init (rp, ns, "o");
  contour (ng, 11);
  couleur (2);

  fm = f[0];			/*   search fmin and fmax */
  xfm = fm;
  for (i = 1; i <= ns; i++)
     {
       if (f[i - 1] > fm)
	 fm = f[i - 1];
       if (f[i - 1] < xfm)
	 xfm = f[i - 1];
     }
  if ((float) fabs (fm - xfm) < 1.0e-15)
    nl = 1;
  for (l = 1; l <= nl; l++)	/*    loop on the level curves */
     {
       if (nl == 1)
	 xln = 0.5F;
       else
	 xln = (l - 1.0F) / (nl - 1.0F);
       xf = xfm + (fm - xfm) * xln;
       for (k = 1; k <= nt; k++)	/*   loop on each femTriangle */
	  {
	    im = 0;
	    for (i = 0; i <= 2; i++)
	       {
		 j = ((i + 1) == 3 ? 0 : i + 1);
		 ik = tr[k - 1][i];
		 jk = tr[k - 1][j];
		 fi = realpart (f[ik]);
		 fj = realpart (f[jk]);
		 if (((fi <= xf) && (fj >= xf)) || ((fi >= xf) && (fj <= xf)))
		    {
		      if ((float) fabs (fi - fj) <= 0.1e-10)	/* one side must be drawn */
			 {
			   rmoveto (rp[ik][0], rp[ik][1]);
			   rlineto (rp[jk][0], rp[jk][1]);
			 }
		      else
			 {
			   xlam = (fi - xf) / (fi - fj);
			   im = im + 1;
			   qp[0][im] = rp[ik][0] * (1.0F - xlam) + rp[jk][0] * xlam;
			   qp[1][im] = rp[ik][1] * (1.0F - xlam) + rp[jk][1] * xlam;
			 }
		    }
	       }
	    if (im >= 2)	/*    draw one segment */
	       {
		 rmoveto (qp[0][1], qp[1][1]);
		 rlineto (qp[0][2], qp[1][2]);
	       }
	  }
     }
  xfmin = xfm; /* solution min */
  xfmax = fm; /* solution max */
  //std::cout << xfmin << " " << xfmax << "\n";
  contour (ng, 11);
  rattente (waitm, xfmin, xfmax);
}

void 
femGraphicDeviceIndependent::equpotd (int *ng, float *f, int nl, int waitm)
{
  float           qp[2][5];
  int             ns = __t->getNumberOfPoints();
  int             nt = __t->getNumberOfCells();
  femTriangle*       tr = __t->tr;
  femPoint*         rp = __t->rp;
  float           xln, fm, xf, xfm, fi, fj, xlam;
  int             im, i, j, k, l, ik, jk;

  Init (rp, ns, "o");
  contour (ng, 11);
  couleur (2);

  fm = f[0];			/*   search fmin and fmax */
  xfm = fm;
  for (i = 1; i <= 3 * nt; i++)
     {
       if (f[i - 1] > fm)
	 fm = f[i - 1];
       if (f[i - 1] < xfm)
	 xfm = f[i - 1];
     }
  if ((float) fabs (fm - xfm) < 1.0e-15)
    nl = 1;
  for (k = 1; k <= nt; k++)	/*   loop on each femTriangle */
     {
       for (l = 1; l <= nl; l++)	/*    loop on the level curves */
	  {
	    if (nl == 1)
	      xln = 0.5F;
	    else
	      xln = (l - 1.0F) / (nl - 1.0F);
	    xf = xfm + (fm - xfm) * xln;
	    im = 0;
	    for (i = 0; i <= 2; i++)
	       {
		 j = ((i + 1) == 3 ? 0 : i + 1);
		 ik = tr[k - 1][i];
		 jk = tr[k - 1][j];
		 fi = realpart (f[3 * (k - 1) + i]);
		 fj = realpart (f[3 * (k - 1) + j]);
		 if (((fi <= xf) && (fj >= xf)) || ((fi >= xf) && (fj <= xf)))
		    {
		      if ((float) fabs (fi - fj) <= 0.1e-10)	/* one side must be drawn */
			 {
			   rmoveto (rp[ik][0], rp[ik][1]);
			   rlineto (rp[jk][0], rp[jk][1]);
			 }
		      else
			 {
			   xlam = (fi - xf) / (fi - fj);
			   im = im + 1;
			   qp[0][im] = rp[ik][0] * (1.0F - xlam) + rp[jk][0] * xlam;
			   qp[1][im] = rp[ik][1] * (1.0F - xlam) + rp[jk][1] * xlam;
			 }
		    }
	       }
	    if (im >= 2)	/*    draw one segment */
	       {
		 rmoveto (qp[0][1], qp[1][1]);
		 rlineto (qp[0][2], qp[1][2]);
	       }
	  }
     }
  contour (ng, 11);
  rattente (waitm);
}

void 
femGraphicDeviceIndependent::showbdy (long nbs, float *cr, long nba, long *arete, float *hh, int waitm)
{
  long            i;
  float           d, x, y, xmx = -(float)1e10, xmn = (float)1e10, ymx = -(float)1e10, ymn = (float)1e10;
  float           xctr, yctr, ray;

  for (i = 0; i < nbs; i++)
     {
       x = cr[2 * i];
       y = cr[2 * i + 1];
       xmx = x > xmx ? x : xmx;
       xmn = xmn > x ? x : xmn;
       ymx = y > ymx ? y : ymx;
       ymn = ymn > y ? y : ymn;
     }
  xctr = (xmx + xmn) / 2;
  yctr = (ymx + ymn) / 2;
  ray = (xmx - xctr) > (ymx - yctr) ? xmx - xctr : ymx - yctr;
  reffecran ();
  cadreortho (xctr, yctr, ray);
  for (i = 0; i < nbs; i++)
     {
       d = 0.1F * hh[i];
       rmoveto (cr[2 * i], cr[2 * i + 1]);
       rlineto (d + cr[2 * i], cr[2 * i + 1]);
       rlineto (d + cr[2 * i], d + cr[2 * i + 1]);
       rlineto (cr[2 * i], d + cr[2 * i + 1]);
       rlineto (cr[2 * i], cr[2 * i + 1]);
     }
  for (i = 0; i < nba; i++)
     {
       rmoveto (cr[2 * arete[2 * i]], cr[2 * arete[2 * i] + 1]);
       rlineto (cr[2 * arete[2 * i + 1]], cr[2 * arete[2 * i + 1] + 1]);
     }
  rattente(waitm);
}
}