femmesh.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: 12-Jul-01 at 10:01:04 by 
//
// 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.
//

#include <cmath>
#include <cstdio>
#include <cstdlib>

#include <femMesh.hpp>
#include <femGraphic.hpp> 

#define i_nint(x) ((long) x >= 0 ? x + 0.5 : x - 0.5)
#define amin(a, b) ((a) < (b) ? (a) : (b))
#define amax(a, b) ((a) < (b) ? (b) : (a))
#define aabs(a) ((a) < 0 ? -(a) : (a))

namespace fem
{
static long     nothing = -1073741824L;

femMesh::femMesh()
    :
    np( 0 ),
    nt( 0 ),
    rp( 0 ),
    tr( 0 ),
    ngt( 0 ),
    ng( 0 )
{
}

femMesh::femMesh( femMesh const& __t )
    :
    np( __t.getNumberOfPoints() ),
    nt( __t.getNumberOfCells() ),
    rp( new femPoint[np] ),
    tr( new femTriangle[nt] ),
    ngt( new int[nt] ),
    ng( new int[np] )
{
    for( int __i = 0; __i < np; __i++)
        {
            rp[__i][0] = __t.rp[__i][0];
            rp[__i][1] = __t.rp[__i][1];
            ng[__i] = __t.ng[__i];
        }
    for( int __i = 0; __i < nt; __i++)
        {
            tr[__i][0] = __t.tr[__i][0];
            tr[__i][1] = __t.tr[__i][1];
            tr[__i][2] = __t.tr[__i][2];
            ngt[__i] = __t.ngt[__i];
        }
    
    
}


void
femMesh::setDimensions( int p, int t )
{
  np = p;
  nt = t;
  
  rp = new femPoint[p];
  tr = new femTriangle[t];
  ng = new int[p];
  ngt = new int[t];
  

}

void 
femMesh::swapWithNeibhor ( int k )
{
  int             k1, j1, j;
  int             next[3];	//={1,2,0};  

  next[0] = 1;
  next[1] = 2;
  next[2] = 0;
  for (k1 = 0; k1 < nt; k1++)
    for (j1 = 0; j1 <= 2; j1++)
      for (j = 0; j <= 2; j++)
	if ((tr[k1][j1] == tr[k][next[j]]) && (tr[k1][next[j1]] == tr[k][j]))
	   {
	     tr[k1][next[j1]] = tr[k][3 - j - next[j]];
	     tr[k][next[j]] = tr[k1][3 - j1 - next[j1]];
	     return;
	   }
}

void 
femMesh::removeBdyT()
{				/* swap diagonals when a femTriangle has 3 bdy sides and reverse order if area<0 */
  int             k, j;
  float           area2;
  femPoint         *q = rp;

  for (k = 0; k < nt; k++)
     {
       area2 = (q[tr[k][1]][0] - q[tr[k][0]][0]) * (q[tr[k][2]][1] - q[tr[k][0]][1])
	 - (q[tr[k][1]][1] - q[tr[k][0]][1]) * (q[tr[k][2]][0] - q[tr[k][0]][0]);
       if (area2 <= 0)
         {
	    j = tr[k][1];
	    tr[k][1] = tr[k][2];
	    tr[k][2] = j;
	  }
       for (k = 0; k < nt; k++)
	 if ((ng[tr[k][0]] != 0) && (ng[tr[k][1]] != 0) && (ng[tr[k][2]] != 0))
	   swapWithNeibhor( k );
     }
}


/*****************************************************************************
I. Fabrique une triangulation (de nbsmax sommets maxi.) a partir de:
--------------------------------------------------------------------
    -cr: tableau de taille 2*nbs (abscisse & ordonnees des points en entree)
    -arete: tableau de taille 2*nba (depart & arrivee des aretes)
      (decrire le contour dans le sens direct: domaine a
            gauche des aretes exterieures)
    -h: tableau de taille nbs (precision relative aux points en entree)


II. Alloue et affecte une structure t de type triangulation:
------------------------------------------------------------
     -t.getNumberOfPoints(): nombre de sommets en sortie
     -t.getNumberOfCells(): nombre de femTriangles en sortie
     -t.rp[i][0]:} abscisse
      t.rp[i][1]:} & ordonnee du ieme point du maillage
     -t.tr[i][j]: numero du jeme point du ieme femTriangle du maillage


III. Renvoie le type d'erreur:
------------------------------
     -  0 : pas d'erreur
     -(-1): pas assez de memoire
     - >0 : erreur mshptg_
  mshptg erreurs
1: mshptg := 'le nb de point est < 3 or > que le nb max de point';
2: mshptg := 'il y des points confondus ';
3: mshptg := 'tout les points sont align閟 ';
7: mshptg := 'bug dans l''algorithme pour retrouver la fronti