femparser.cpp

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

}

noeudPtr 
femParser::instruction ()
{
  noeudPtr           l1 = NULL, l2 = NULL, l3 = NULL, res = NULL;
  Symbol          thesym;
  ident          *theid = NULL;
  float           thecst;
  char           *thechaine;

  switch (cursym)
     {
     case si:
       nextsym ();
       flag.si = 1;		/* this forbids vectors in logical expression in if then else statements */
       l1 = expr ();
       flag.si = 0;
       match (alors);
       l2 = instruction ();
       if (cursym == autrement)
	  {
	    nextsym ();
	    l3 = instruction ();
	  }
       else
	 l3 = NULL;
       plante (&res, si, 0.F, 0, NULL, pt, l1, l2, l3, NULL);
       break;

     case loop:
       nextsym ();
       match (lpar);
       thecst = curcst;
       match (cste);
       match (rpar);
       l1 = instruction ();
       plante (&res, loop, thecst, 0, NULL, pt, l1, NULL, NULL, NULL);
       break;

     case newvar:
     case oldvar:
     case fdecl:
       theid = curident;
       nextsym ();
       if (cursym == fctdef)
	  {
  		if (strcmp(theid->name,"x")&&strcmp(theid->name,"y")) // moving the mesh
	    	theid->symb = fdecl;
	    if (!flag.fct)
	      match (becomes);
	    match (fctdef);
	    l1 = expr ();
	    plante (&res, fctdef, 0.F, 0, theid, pt, l1, NULL, NULL, NULL);
	  }
       else
	  {
	    if (theid->symb == fdecl)
	      theid->symb = oldvar;
	    match (becomes);
	    l1 = expr ();
	    plante (&res, becomes, 0.F, 0, theid, pt, l1, NULL, NULL, NULL);
	  }
       break;

     case fctfile:		/* when part of the program is stored (not documented) */
       res = fctfileproc ();
       break;

     case symb_bdy:
     case polygon:
       res = symb_bdyproc ();
       break;

     case symb_build:
       if (!flag.build)
	  {
	    sprintf (errbuf, "line %d: Illegal use of symbol 'buildmesh'", numligne);
	    erreur (errbuf);
	  }
       nba = 0;
       nbsd = 0;
       nbs = 0;
       nextsym ();
       match (lpar);
       l1 = expr ();
       if(cursym==comma) {nextsym(); l2 = expr();} else l2=NULL;
       plante (&res, symb_build, 0.F, 0, NULL, pt, l1, l2, NULL, NULL);
       match (rpar);
       flag.build = 0;
       flag.bdy = 0;
       flag.fct = 1;
       flag.onbdy = 1;
       flag.solv = 1;
       break;

     case chargmsh:
     case sauvmsh:
     case symb_exec:
       res = diskmshproc ();
       break;

     case charge:
       if (!flag.fct)
	  {
	    sprintf (errbuf, "Illegal use of symbol %s\n", mesg[cursym]);
	    erreur (errbuf);
	  }
       thesym = cursym;
       nextsym ();
       match (lpar);
       thechaine = curchaine;
       nextsym ();
       match (comma);
       theid = curident;
       switch (theid->symb)
	  {
	  case newvar:
	  case oldvar:
	    theid->symb = fdecl;
	    nextsym ();
	    break;

	  default:
	    match (fdecl);
	  }
       if (cursym == comma)
	 {
	   nextsym();
	   l1 = expr();
	 }
       match (rpar);
       plante (&res, thesym, 0.F, 0, theid, thechaine, l1, NULL, NULL, NULL);
       break;

     case sauve:
       if (!flag.fct)
	  {
	    sprintf (errbuf, "line %d: Unexpected symbol %s\n", numligne, mesg[cursym]);
	    erreur (errbuf);
	  }
       thesym = cursym;
       nextsym ();
       match (lpar);
       thechaine = curchaine;
       nextsym ();
       match (comma);
       if (cursym == oldvar)
	  {
	    thecst = 1.F;
	    std::ofstream        file (thechaine);

	    file << ' ' << std::endl;
	    file.close ();
	  }			// clean the file
       else
	 thecst = 0.F;
       l1 = expr ();
       if (cursym == comma)
	 {
	   nextsym();
	   l2 = expr();
	 }
       match (rpar);
       plante (&res, thesym, thecst, 0, NULL, thechaine, l1, l2, NULL, NULL);
       break;

     case trace:
     case trace3d:
       if (!flag.fct)
	  {
	    sprintf (errbuf, "line %d: Unexpected symbol %s\n", numligne, mesg[cursym]);
	    erreur (errbuf);
	  }
       thesym = cursym;
       nextsym ();
       match (lpar);
       l1 = expr ();
       plante (&res, thesym, 0.F, 0, NULL, pt, l1, NULL, NULL, NULL);
       match (rpar);
       break;

     case adaptmesh:
       {
         int i = 0;
       noeudPtr li[4] = {NULL,NULL,NULL,NULL};
       nextsym ();
       match (lpar);
       li[i++] = expr ();
       while (cursym == comma)
         {
           nextsym();
           li[i] = expr();
           i++;
         }
       plante (&res, adaptmesh, 0.F, 0, NULL, pt, li[0], li[1], li[2], li[3]);
       match (rpar);
       }
       break;
       
     case symb_dch:
       if (flag.syst)
	 res = symb_dchproc ();
       else
	  {
	    sprintf (errbuf, "line %d: onbdy must be nested within solve(){...}\n", numligne);
	    erreur (errbuf);
	  }
       break;


     case symb_pde:
       res = symb_pdeproc ();
       break;

     case symb_solv:
     case sauvetout:
       thesym = cursym;
       if (flag.syst)
	 		erreur ("Embedded solvers are not allowed");
       flag.syst = 1;
       res = preparesolve ();
       l2 = instruction ();
       plante (&res, thesym, 0.F, N, NULL, NULL, res, l2, NULL, NULL);
       flag.syst = 0;
       N = 1;
       break;

     case varsolve:
     //case sauvetout:
       thesym = cursym;
       if (flag.syst)
	 		erreur ("Embedded solvers are not allowed");
       flag.syst = 1;
       res = prepvarsolve ();
       l2 = instruction ();
       match(colon);
       l3 = expr();
       plante (&res, colon, 0.F, N, NULL, NULL, res, l2, l3, NULL);
       flag.syst = 0;
       N = 1;
       break;

     case arret:
     case changewait:
     case wait_state:
     case nowait:
     case symb_complex:
     case symb_precise:
       plante (&res, cursym, 0.F, 0, NULL, pt, NULL, NULL, NULL, NULL);
       nextsym ();
       break;

     case lbrace:
       do
	  {
	    nextsym ();
	    l1 = instruction ();
	    if (res == NULL)
	      res = l1;
	    else if (l1)
	      plante (&res, lbrace, 0.F, 0, NULL, pt, res, l1, NULL, NULL);
	  }
       while (cursym == semicolon);
       match (rbrace);
       break;

     case semicolon:
       nextsym ();
     case rbrace:
     case _end:
       break;

     default:
       sprintf (errbuf, "line %d: Cannot use this symbol to begin an expression: %s",
		numligne, mesg[cursym]);
       erreur (errbuf);
     }
  return res;
}

void 
femParser::defbdy (noeudPtr s)
{
  int             i, j, j0=0, j1=0; 
  float           gch, dte, dcr;
  static int      refl0 = 0,refl1,refs;
  static int      first=0;
  
  if (first==0) {
    first=1;
    for (i = 0;i < refbdy;i++)
      for (j = 0;j < refbdy;j++)
        tabref[i][j]  = 0;
  }
  
  gch = realpart (eval (s->l1));
  dte = realpart (eval (s->l2));
  dcr = realpart (eval (s->l3)) - 1;
  (variables.ng)->value = (int) (realpart (s->value));
  refl1 =  int(realpart ((variables.ng)->value));
  refs  = int(s->junk);
  if (refl0 != 0 && refs != 0)
    tabref[refl0][refs] = refl0;
  if (refl1 !=0 && refs != 0)
    tabref[refs][refl1] = refl1;
  for (i = 0; i <= dcr; i++)
     {
       (variables.t)->value = gch + i * (dte - gch) / dcr;
       eval (s->l4);
       cr[2 * nbs] = realpart ((variables.x)->value);
       cr[2 * nbs + 1] = realpart ((variables.y)->value);
       if (i == 0 && nbs == 0 && int(s->junk) > 0)
	 ngbdy[nbs] = int(s->junk);
       else
	 ngbdy[nbs] = (int) realpart ((variables.ng)->value);
       if (j1 = __mesh.check(cr, nbs), j1 == -1)
	 j1 = nbs++;
       else 
	 if (i!=0)
	   tabref[refl1][ngbdy[j1]]=refl1;
       if (i == 0 && j1 != -1 && nbs != 1 && int(s->junk) > 0)
	 ngbdy[j1] = int(s->junk);
       if (i != 0)
	  {
	    arete[2 * nba] = j0;
	    arete[2 * nba + 1] = j1;
	    nba++;
	    hh[j0] = 0.66F * sqrt (::pow (cr[2 * j1] - cr[2 * j0], 2) + ::pow (cr[2 * j1 + 1] - cr[2 * j0 + 1], 2));
	  }
       j0 = j1;
     }
  hh[j1] = hh[j0];
  sd[2 * nbsd] = nba - 1;
  sd[2 * nbsd + 1] = nbsd + 1;
  nbsd++;
  refl0 = refl1;
}

void 
femParser::defbdybypoint (noeudPtr s)
{
  int  dte, i = -1, j0=0, j1=0;
  int   ngb = (int) (realpart (s->value)); 
  float xl;
  float* cr1 = new float[2*MAXBDYVERTEX+1];
  float* cr2 = new float[2*MAXBDYVERTEX+1];
  int dcr=readpoints(s->path, cr1, MAXBDYVERTEX);
  if( dcr <=0 )       
  {
    if(dcr<0) 
      sprintf (errbuf, "Too many points in file %s ", s->path); 
    else 
      sprintf (errbuf, "Could not read points from file %s ", s->path); 
    erreur (errbuf);
  } 
  if(!(s->l1))
    dte = 0; 
  else
    dte = (int)(realpart (eval (s->l1)));
  dcr--;
  for(j1=0; j1<dcr;j1++)
  	for(j0=1;j0<=dte+1;j0++)
	{ 
		i=j0+ (dte+1)*j1;
		xl = 1.0F - (j0-1.0F)/(dte+1.0F);
		cr2[2*i-2] = cr1[2*j1]*xl + (1.0F-xl)*cr1[2*j1+2];
		cr2[2*i-1] = cr1[2*j1+1]*xl + (1.0F-xl)*cr1[2*j1+3];
	}
  cr2[2*i] = cr1[2*dcr]; cr2[2*i+1] = cr1[2*dcr+1];
  for (i = 0; i <= (1+dte)*dcr; i++)
     {
       (variables.ng)->value = ngb;
       cr[2 * nbs] = cr2[2*i];
       cr[2 * nbs + 1] = cr2[2*i + 1];
       ngbdy[nbs] = ngb;
       if (j1 = __mesh.check(cr, nbs), j1 == -1)  j1 = nbs++;
       if (i != 0)
	  {
	    arete[2 * nba] = j0;
	    arete[2 * nba + 1] = j1;
	    nba++;
	    hh[j0] = 0.66F * (float)sqrt ( (cr[2 * j1] - cr[2 * j0]) * (cr[2 * j1] - cr[2 * j0]) 
	    				+  (cr[2 * j1 + 1] - cr[2 * j0 + 1]) *  (cr[2 * j1 + 1] - cr[2 * j0 + 1]));
	  }
       j0 = j1;
     }
  hh[j1] = hh[j0];
  sd[2 * nbsd] = nba - 1;
  sd[2 * nbsd + 1] = nbsd + 1;
  nbsd++;
  delete [] cr1;
  delete [] cr2;
}

void 
femParser::initparam ()
{
  long tnp = flag.precise ? 3 * (long) __mesh.getNumberOfCells() : (long) __mesh.getNumberOfPoints();

  if (flag.param == 0)
    {
      __fem = new FEM( &__mesh, flag.precise );
    }
  if (flag.complexe)
     {
       if ((N == 1) && ((flag.param == 0) || (flag.param == 2)))
	  {
	    param.p1c = new creal[tnp];
	    param.c1c = new creal[tnp];
	    param.g1c = new creal[tnp];
	    param.f1c = new creal[tnp];
	    param.b1c = new creal[tnp];
	    param.nuyx1c = new creal[tnp];
	    param.nuxx1c = new creal[tnp];
	    param.nuxy1c = new creal[tnp];
	    param.nuyy1c = new creal[tnp];
	    param.a11c = new creal[tnp];
	    param.a21c = new creal[tnp];
	    param.sol1c = new creal[tnp];
            if (param.fplot)
              {delete [] param.fplot;param.fplot = NULL;}
            param.fplot = new float[tnp];
            
	    for (int i = 0; i < tnp; i++)
	       {
		 param.p1c[i] = 0;
		 param.c1c[i] = 0;
		 param.g1c[i] = 0;
		 param.f1c[i] = 0;
		 param.b1c[i] = 0;
		 param.nuyx1c[i] = 0;
		 param.a11c[i] = 0;
		 param.a21c[i] = 0;
		 param.nuxx1c[i] = 0;
		 param.nuxy1c[i] = 0;
		 param.nuyy1c[i] = 0;
		 param.sol1c[i] = 0;
	       }
	    flag.param += N;
	  }
     }
  else if ((N == 1) && ((flag.param == 0) || (flag.param == 2)))
     {
       param.p1 = new float[tnp];
       param.c1 = new float[tnp];
       param.g1 = new float[tnp];
       param.f1 = new float[tnp];
       param.b1 = new float[tnp];
       param.nuyx1 = new float[tnp];
       param.nuxx1 = new float[tnp];
       param.nuxy1 = new float[tnp];
       param.nuyy1 = new float[tnp];
       param.a11 = new float[tnp];
       param.a21 = new float[tnp];
       param.sol1 = new float[tnp];
       if (param.fplot)
         delete [] param.fplot;
       param.fplot = new float[tnp];
       for (int i = 0; i < tnp; i++)
	  {
	    param.p1[i] = 0.F;
	    param.c1[i] = 0.F;
	    param.g1[i] = 0.F;
	    param.f1[i] = 0.F;
	    param.b1[i] = 0.F;
	    param.nuyx1[i] = 0.F;
	    param.a11[i] = 0.F;
	    param.a21[i] = 0.F;
	    param.nuxx1[i] = 0.F;
	    param.nuxy1[i] = 0.F;
	    param.nuyy1[i] = 0.F;
	    param.sol1[i] = 0.F;
	  }
       flag.param += N;
     }
  else if ((N == 2) && (flag.param < 2))
     {
       param.p2.init (tnp);	//= new cvect[tnp];//(cvect*)safecalloc(tnp,sizecvect);
       param.c2.init (tnp);	//= new cmat[tnp];//(cmat*)safecalloc(tnp,sizecmat);
       param.g2.init (tnp);	//=new cvect[tnp];//(cvect*)safecalloc(tnp,sizecvect);
       param.f2.init (tnp);	//=new cvect[tnp];//(cvect*)safecalloc(tnp,sizecvect);
       param.b2.init (tnp);	//=new cmat[tnp];//(cmat*)safecalloc(tnp,sizecmat);
       param.nuyx2.init (tnp);	//=new cmat[tnp];//(cmat*)safecalloc(tnp,sizecmat);
       param.nuxx2.init (tnp);	//=new cmat[tnp];//(cmat*)safecalloc(tnp,sizecmat);
       param.nuxy2.init (tnp);	//=new cmat[tnp];//(cmat*)safecalloc(tnp,sizecmat);
       param.nuyy2.init (tnp);	//=new cmat[tnp];//(cmat*)safecalloc(tnp,sizecmat);
       param.a12.init (tnp);	//=new cmat[tnp];//(cmat*)safecalloc(tnp,sizecmat);
       param.a22.init (tnp);	//=new cmat[tnp];//(cmat*)safecalloc(tnp,sizecmat);
       param.sol2.init (tnp);	//=new cvect[tnp];//(cvect*)safecalloc(tnp,sizecvect);
       if (param.fplot)
         delete [] param.fplot;
       param.fplot = new float[tnp];
       for (int i = 0; i < tnp; i++)
	  {
	    param.p2[i] = 0.F;
	    param.c2[i] = 0.F;
	    param.g2[i] = 0.F;
	    param.f2[i] = 0.F;
	    param.b2[i] = 0.F;
	    param.nuyx2[i] = 0.F;
	    param.a12[i] = 0.F;
	    param.a22[i] = 0.F;
	    param.nuxx2[i] = 0.F;
	    param.nuxy2[i] = 0.F;
	    param.nuyy2[i] = 0.F;
	    param.sol2[i] = 0.F;
          }
       flag.param += N;
     }