femparser.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: 24-Oct-01 at 18:49:36 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.
//
#include <list>

#include <iostream>
#include <fstream>


#include <femCommon.hpp> 
#include <femMisc.hpp> 

#include <femIdentifier.hpp>
#include <femMesh.hpp>
#include <femLexical.hpp> 
#include <femGraphic.hpp> 
#include <femGraphicDeviceIndependent.hpp> 
#include <femDisk.hpp> 
#include <femParser.hpp> 

#include <femSolver.hpp> 


#if defined(ADAPT)
# include <header.hxx>
# include <params.hxx>
# include <triangl.hxx>
# include <m_t0.hxx>
# include <cad.hxx>

parameter adapt_param;	
#endif /* ADAPT */

#include <femFunction.hpp>

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


#define sqr(x) ((x)*(x))
#define Mini(x, y) (x) < (y) ? (x) : (y)
#define Maxi(x, y)  (x) > (y) ? (x) : (y)
#define MAXSD 200
#define MAXTREENODES 2000	/* of analyser */
#define MAXBDYVERTEX 3000	/* of contour */
#define MAXSYS 100		/* nb of equations in systems */

namespace fem
{
int             N= 1;
int             N2 = 1;		// nb of variables of system 

ident          *systable[MAXSYS];	// table of pointers to variables of the system
char           errbuf[1024]; // 1024 seems to be big enough to hold the messages



#define penal (float)1.e-14
#define xi1 1/sqrt(3)
#define xi2 1/sqrt(3)

/*
  for adapt
 */
  int            tabref[refbdy][refbdy];

  

extern int      bug;
extern creal    sqrtofminus1;

char           *mesg[102] =
{"'('", "')'", "'{'", "'}'", "'constant'",
 "'new variable'", "'old variable'", "'+'",
 "'-'", "'*'", "'/'", "'%'","'<'", "'<='", "'>'", "'>='",
 "'=='", "'!='", "','", "';'", "'sin'", "'cos'",
 "'log'", "'exp'", "'sqrt'", "'abs'", "'^'", "'acos'",
 "'asin'", "'tan'", "'and'", "'cosh'", "sinh'", "'tanh'",
 "'or'", "'min'", "'max'", "'dx'", "'dy'", "'if'",
 "'then'", "'else'", "'iter'", "'error'", "'final'", "':='",
 "'fonction'", "'bdy'", "'buildmesh'", "'atan'", "'='",
 "'solve'", "'id'", "'dnu'", "'id'", "'laplace'", "'div'",
 "'plot'", "'changewait'", "'plot3d'", "'''", "'save'", "'load'",
 "'savemesh'", "'loadmesh'", "'halt'", "'include'", "'dx'",
 "'dy'", "'convect'", "'evalfct'", "'exec'", "'saveall'",
 "'user'", "'Re'", "'Im'", "'system'", "'pde'", "'id_bdy'",
 "'dnu_bdy'", "'dxx'", "'dyy'", "'dxy'", "'dyx'", "'complex'", "'precise'"
 ,"'scal'", "'nx'", "'ny'", "'one'", "'wait'", "'nowait'", "'rhsconvect'",
 "'adaptmesh'","'polygon'","'intt'","'int'","']'","'['", "'varsolve'","'penal'","':'"};


femParser::femParser(  )
  :
  __tree(),
  __function_list(),
  __mesh(),
  __graph( new femGraphicDeviceIndependent( &__mesh ) ),
  pt( 0 ),
  nbsd( 0 ),
  nbs( 0 ),
  nba( 0 ),
  Iter( 1 ),
  waitm( 0 ),
  __text( 0 ),
  __graphic_type( FEM_GRAPHIC )
{
  numnoeuds = 0;
  waitm = 1;
  pt = NULL;
  noeuds = new noeudPtr[MAXTREENODES];
  sd = new long[2*MAXSD];
  arete = new long[2 * MAXBDYVERTEX];
  ngbdy = new int[MAXBDYVERTEX];
  cr = new float[2*MAXBDYVERTEX+1];
  hh = new float[MAXBDYVERTEX];
  memset (hh, 0, MAXBDYVERTEX*sizeof(float));
  flag.si = 0;
  flag.syst = 0;
  flag.param = 0;
  flag.complexe = 0;
  flag.precise = 0;
}

femParser::~femParser()
{
  bucheron(__tree);
  libere();

  if ( __graphic_type == FEM_GRAPHIC )
    {
      closegraphique();
    }
  delete __graph;
}


/*
 * parse the buffer
 */
void
femParser::parse ()
{
  
  initlex( __text );
  __tree = instruction();
  chvar();

  if ( __graphic_type == FEM_GRAPHIC )
    {
      initgraphique();
    }
  
  eval(__tree);
}

void 
femParser::chvar ()
{
  int             i;
  char           *c;
  Complex I(0,1);

  for (i = 0; i < numidents; i++)
     {
       c = idents[i].name;
       if (!strcmp (c, "x"))
	 variables.x = idents + i;
       if (!strcmp (c, "y"))
	 variables.y = idents + i;
       if (!strcmp (c, "t"))
	 variables.t = idents + i;
       if (!strcmp (c, "ib"))
	 variables.ng = idents + i;
       if (!strcmp (c, "region"))
	 variables.region = idents + i;
       if (!strcmp (c, "iv"))
	 variables.cursom = idents + i;
       if (!strcmp (c, "nx"))
	 variables.nx = idents + i;
       if (!strcmp (c, "ny"))
	 variables.ny = idents + i;
       if (!strcmp (c, "nexist"))
	 variables.ne = idents + i;
       if (!strcmp (c, "I"))
	 variables.I = idents + i;
       if (!strcmp (c, "pi"))
	 variables.pi = idents + i;
     }
  (variables.I)->value = I;
  (variables.pi)->value = 4.0F * atan (1.0F);
}

void 
femParser::plante (noeudPtr * res, Symbol s, creal v, long j, ident * n, char *pt,
	noeud * L1, noeud * L2, noeud * L3, noeud * L4)
{
  noeudPtr p = new noeud;
  
  if (numnoeuds == MAXTREENODES)
    erreur ("Tree is too big...");
  p->symb = s;
  p->value = v;
  p->junk = j;
  p->name = n;
  p->l1 = L1;
  p->l2 = L2;
  p->l3 = L3;
  p->l4 = L4;
  if (pt)
     {
       p->path = new char[(strlen (pt) + 1)];
       strcpy (p->path,pt);
     }
  else
    p->path = NULL;
  *res = noeuds[numnoeuds++] = p;
}




void 
femParser::match (Symbol s)
{
  if (s == cursym)
    nextsym ();
  else
     {
       sprintf (errbuf, "line %d: Unexpected symbol: %s Instead of %s",
		numligne, mesg[cursym], mesg[s]);
       erreur (errbuf);
     }
}

noeudPtr           expr ();

noeudPtr 
femParser::facteur ()
{
  Symbol          thesym = cursym;
  ident          *theid = curident;
  noeudPtr           l1, l2, l3, l4, res = NULL,l[5]={NULL,NULL,NULL,NULL,NULL};
  int count;

  switch (thesym)
     {
     case lpar:
       nextsym ();
       res = expr ();
       match (rpar);
       break;

     case symb_user:
       nextsym ();
       match (lpar);
       l1 = expr ();
       if (cursym == comma)
	  {
	    nextsym ();
	    l2 = expr ();
	  }
       else
	 l2 = NULL;
       plante (&res, thesym, 0.F, 0, NULL, pt, l1, l2, NULL, NULL);
       match (rpar);
       break;

     case one:
     case sine:
     case cosine:
     case atane:
     case exponential:
     case logarithm:
     case absolute:
     case root:
     case acose:
     case asine:
     case tane:
     case coshe:
     case sinhe:
     case tanhe:
     case partreal:
     case partimag:
     case partial_x:
     case partial_y:
     case penall:
       nextsym ();
       match (lpar);
       if ((cursym != fdecl) && ((thesym == partial_x) || (thesym == partial_y)))
	  {
	    sprintf (errbuf, "Line %d : Array-function name expected", numligne);
	    erreur (errbuf);
	  }
       l1 = expr ();
       plante (&res, thesym, 0.F, 0, NULL, pt, l1, NULL, NULL, NULL);
       match (rpar);
       break;
     case gint:
     case bint:
       if (cursym == bint)
         {
           /*
            * integrale de bord: il faut specifier au moins un bord
            */
           nextsym();
           match (lpar);
         }
       else
         {
           /*
            * integrale globale: si on ne specifie pas de sous domaines on
            *                    fait l'integrale sur tout le domaine
            */
           nextsym();
           if (cursym == lpar) nextsym();
         }
       if (cursym != bracketl)
         {
           count = 0;
           while (cursym != rpar)
             {
               if (count == 3)
                 {
                   sprintf (errbuf, "Line %d : More than 3 boundaries.", numligne);
                   erreur (errbuf);
                 }
               l[count++] = expr();
               if (cursym == comma) nextsym ();
             }
           if (count == 0 && cursym == bint)
             {
               sprintf (errbuf, "Line %d : at least one bdy ref is expected.", numligne);
               erreur (errbuf);
             }
           nextsym();
         }
       match (bracketl);
       if (cursym != fdecl)
         	{
           		sprintf (errbuf, "Line %d : Array-function name expected", numligne);
           		erreur (errbuf);
         	}
        l[3] = expr();
       	theid = NULL;
       	if (cursym == comma) 
       		{
       			nextsym ();
       		    if (cursym != fdecl)
         		{
           			sprintf (errbuf, "Line %d : Array-function name expected", numligne);
           			erreur (errbuf);
         		}
         		theid=curident;
         		nextsym();
         	}	
       plante (&res, thesym, 0.F, 0, theid, pt, l[0], l[1], l[2], l[3]);
       match (bracketr);
       break;
     case mini:
     case maxi:
     case prodscal:
       nextsym ();
       match (lpar);
       if ((thesym == prodscal) && (cursym != fdecl))
	  {
	    sprintf (errbuf, "Line %d : Array-function name expected", numligne);
	    erreur (errbuf);
	  }
       l1 = expr ();
       match (comma);
       if ((thesym == prodscal) && (cursym != fdecl))
	  {
	    sprintf (errbuf, "Line %d : Array-function name expected", numligne);
	    erreur (errbuf);
	  }
       l2 = expr ();
       plante (&res, thesym, 0.F, 0, NULL, pt, l1, l2, NULL, NULL);
       match (rpar);
       break;

     case symb_convect:
     case rhsconvect:
       nextsym ();
       match (lpar);
       if (cursym != fdecl)
	  {
	    sprintf (errbuf, "Line %d : Array-function name expected", numligne);
	    erreur (errbuf);
	  }
       l1 = expr ();
       match (comma);
       if (cursym != fdecl)
	  {
	    sprintf (errbuf, "Line %d : Array-function name expected", numligne);
	    erreur (errbuf);
	  }
       l2 = expr ();
       match (comma);
       if (cursym != fdecl)
	  {
	    sprintf (errbuf, "Line %d : Array-function name expected", numligne);
	    erreur (errbuf);
	  }
       l3 = expr ();
       match (comma);
       l4 = expr ();
       plante (&res, thesym, 0.F, 0, NULL, pt, l1, l2, l3, l4);
       match (rpar);
       break;

     case cste:
       plante (&res, thesym, curcst, 0, NULL, pt, NULL, NULL, NULL, NULL);
       nextsym ();
       break;

     case newvar:
       match (oldvar);
     case oldvar:
       plante (&res, thesym, 0.F, 0, curident, pt, NULL, NULL, NULL, NULL);
       nextsym ();
       break;

     case fdecl:
       nextsym ();
       if (cursym == lpar)	/* for f(x,y) */
	  {
	    match (lpar);
	    l1 = expr ();
	    match (comma);
	    l2 = expr ();
	    plante (&res, evalfct, 0.F, 0, theid, pt, l1, l2, NULL, NULL);
	    match (rpar);
	  }
       else if (!flag.si)
	 plante (&res, thesym, 0.F, 0, theid, pt, NULL, NULL, NULL, NULL);
       else
	  {
	    sprintf (errbuf, "line %d: Array-functions are not allowed in the logical expression of an if statement(use max,min...)", numligne);
	    erreur (errbuf);
	  }
       break;
     default :
       break;
     }
  if (cursym == expo)
     {
       nextsym ();
       l2 = facteur ();
       plante (&res, expo, 0.F, 0, NULL, pt, res, l2, NULL, NULL);
     }
  return res;
}

noeudPtr 
femParser::terme ()
{
  noeudPtr           l2, res = facteur ();

  while ((cursym == star) || (cursym == slash) || (cursym == modulo))
     {
       Symbol          thesym = cursym;

       nextsym ();
       l2 = facteur ();
       plante (&res, thesym, 0.F, 0, NULL, pt, res, l2, NULL, NULL);
     }
  return res;
}


noeudPtr 
femParser::exprarith ()
{
  noeudPtr           l2, res;

  switch (cursym)
     {
     case op_plus:
       nextsym ();