femparser.cpp

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

		   {
			case symb_lapl:
		     if (!s->path)
			{
			  (param.nuxx1c)[cursom] += -(float)oper * eval (s->l2);
			  (param.nuyy1c)[cursom] += -(float)oper * eval (s->l2);
			}
		     else
			{
			  (param.nuxx1c)[cursom] += -(float)oper / eval (s->l2);
			  (param.nuyy1c)[cursom] += -(float)oper / eval (s->l2);
			}
		     break;
		   case div_x:
		     if (!s->path)
		       (param.a11c)[cursom] += (float)oper * (eval (s->l2));
		     else
		       (param.a11c)[cursom] += (float)oper / (eval (s->l2));
		     break;
		   case div_y:
		     if (!s->path)
		       (param.a21c)[cursom] += (float)oper * (eval (s->l2));
		     else
		       (param.a21c)[cursom] += (float)oper / (eval (s->l2));
		     break;
		   case d_xx:
			  if (!s->path)
		       (param.nuxx1c)[cursom] += -(float)oper * (eval (s->l2));
		     else
		       (param.nuxx1c)[cursom] += -(float)oper / (eval (s->l2));
		     break;
		   case d_xy:
		     if (!s->path)
		       (param.nuxy1c)[cursom] += -(float)oper * (eval (s->l2));
		     else
		       (param.nuxy1c)[cursom] += -(float)oper / (eval (s->l2));
		     break;
		   case d_yx:
		     if (!s->path)
		       (param.nuyx1c)[cursom] += -(float)oper * (eval (s->l2));
		     else
		       (param.nuyx1c)[cursom] += -(float)oper / (eval (s->l2));
		     break;
		   case d_yy:
		     if (!s->path)
		       (param.nuyy1c)[cursom] += -(float)oper * (eval (s->l2));
		     else
		       (param.nuyy1c)[cursom] += -(float)oper / (eval (s->l2));
		     break;
		   case symb_id:
		     if (!s->path)
				 (param.b1c)[cursom] += (float)oper * eval (s->l2);
		     else
		       (param.b1c)[cursom] += (float)oper / eval (s->l2);
		     break;
		   default:
		     break;
		   }
/*  if(N==2)
   switch(thesym)
   {
   case symb_lapl:
   if(!s->path){
   (param.nuxx2)[cursloc](im,jm) += - oper * eval(s->l2);
   (param.nuyy2)[cursloc](im,jm) += - oper * eval(s->l2);
   }
   else    {
   (param.nuxx2)[cursloc](im,jm) += - oper / eval(s->l2);
   (param.nuyy2)[cursloc](im,jm) += - oper / eval(s->l2);
   }
   break;
   case div_x:
   if(!s->path) (param.a12)[cursloc](im,jm)+= oper * (eval(s->l2));
   else     (param.a12)[cursloc](im,jm)+= oper / (eval(s->l2));
   break;
   case div_y:
	if(!s->path) (param.a22)[cursloc](im,jm)+= oper *  (eval(s->l2));
   else     (param.a22)[cursloc](im,jm)+= oper /  (eval(s->l2)); 
   break;
   case d_xx:
   if(!s->path) (param.nuxx2)[cursloc](im,jm) += -oper *  (eval(s->l2)); 
   else     (param.nuxx2)[cursloc](im,jm) += -oper /  (eval(s->l2)); 
   break;
   case d_xy:
   if(!s->path) (param.nuxy2)[cursloc](im,jm)+= -oper *  (eval(s->l2)); 
   else     (param.nuxy2)[cursloc](im,jm)+= -oper /  (eval(s->l2)); 
   break;
   case d_yx:
   if(!s->path) (param.nuyx2)[cursloc](im,jm)+= -oper *  (eval(s->l2)); 
   else     (param.nuyx2)[cursloc](im,jm)+= -oper /  (eval(s->l2)); 
   break;
   case d_yy:
   if(!s->path) (param.nuyy2)[cursloc](im,jm) += -oper *  (eval(s->l2)); 
   else     (param.nuyy2)[cursloc](im,jm) += -oper /  (eval(s->l2)); 
   break;
   case symb_id:
   if(!s->path) (param.b2)[cursloc](im,jm)+= oper * eval(s->l2);
   else     (param.b2)[cursloc](im,jm)+= oper / eval(s->l2);
   break;
   }
 */ 
		 }
	 else
	    {
	      if (N == 1)
		switch (thesym)
		   {
		   case symb_lapl:
		     if (!s->path)
			{
			  (param.nuxx1)[cursom] += -(float)oper * realpart (eval (s->l2));
			  (param.nuyy1)[cursom] += -(float)oper * realpart (eval (s->l2));
			}
		     else
			{
			  (param.nuxx1)[cursom] += -(float)oper / realpart (eval (s->l2));
			  (param.nuyy1)[cursom] += -(float)oper / realpart (eval (s->l2));
			}
		     break;
		   case div_x:
		     if (!s->path)
		       (param.a11)[cursom] += (float)oper * (realpart (eval (s->l2)));
		     else
		       (param.a11)[cursom] += (float)oper / (realpart (eval (s->l2)));
		     break;
		   case div_y:
			  if (!s->path)
		       (param.a21)[cursom] += (float)oper * (realpart (eval (s->l2)));
		     else
		       (param.a21)[cursom] += (float)oper / (realpart (eval (s->l2)));
		     break;					  
		   case d_xx:
		     if (!s->path)
		       (param.nuxx1)[cursom] += -(float)oper * (realpart (eval (s->l2)));
		     else
		       (param.nuxx1)[cursom] += -(float)oper / (realpart (eval (s->l2)));
		     break;
		   case d_xy:
		     if (!s->path)
		       (param.nuxy1)[cursom] += -(float)oper * (realpart (eval (s->l2)));
		     else
		       (param.nuxy1)[cursom] += -(float)oper / (realpart (eval (s->l2)));
		     break;
		   case d_yx:
		     if (!s->path)
		       (param.nuyx1)[cursom] += -(float)oper * (realpart (eval (s->l2)));
		     else
		       (param.nuyx1)[cursom] += -(float)oper / (realpart (eval (s->l2)));
		     break;
		   case d_yy:
		     if (!s->path)
				 (param.nuyy1)[cursom] += -(float)oper * (realpart (eval (s->l2)));
		     else
		       (param.nuyy1)[cursom] += -(float)oper / (realpart (eval (s->l2)));
		     break;
		   case symb_id:
		     if (!s->path)
		       (param.b1)[cursom] += (float)oper * realpart (eval (s->l2));
		     else
		       (param.b1)[cursom] += (float)oper / realpart (eval (s->l2));
		     break;
		   default:
		     break;
		   }
	      if (N == 2)
		switch (thesym)
		   {
		   case symb_lapl:
		     if (!s->path)
			{
			  (param.nuxx2)[cursom] (im, jm) += -(float)oper * realpart (eval (s->l2));
			  (param.nuyy2)[cursom] (im, jm) += -(float)oper * realpart (eval (s->l2));
			}
		     else
			{
			  (param.nuxx2)[cursom] (im, jm) += -(float)oper / realpart (eval (s->l2));
			  (param.nuyy2)[cursom] (im, jm) += -(float)oper / realpart (eval (s->l2));
			}
		     break;
		   case div_x:
		     if (!s->path)
		       (param.a12)[cursom] (im, jm) += (float)oper * (realpart (eval (s->l2)));
		     else
		       (param.a12)[cursom] (im, jm) += (float)oper / (realpart (eval (s->l2)));
		     break;
		   case div_y:
		     if (!s->path)
		       (param.a22)[cursom] (im, jm) += (float)oper * (realpart (eval (s->l2)));
		     else
		       (param.a22)[cursom] (im, jm) += (float)oper / (realpart (eval (s->l2)));
		     break;
		   case d_xx:
		     if (!s->path)
		       (param.nuxx2)[cursom] (im, jm) += -(float)oper * (realpart (eval (s->l2)));
		     else
		       (param.nuxx2)[cursom] (im, jm) += -(float)oper / (realpart (eval (s->l2)));
		     break;
		   case d_xy:
		     if (!s->path)
		       (param.nuxy2)[cursom] (im, jm) += -(float)oper * (realpart (eval (s->l2)));
		     else
				 (param.nuxy2)[cursom] (im, jm) += -(float)oper / (realpart (eval (s->l2)));
			  break;
			case d_yx:
			  if (!s->path)
				 (param.nuyx2)[cursom] (im, jm) += -(float)oper * (realpart (eval (s->l2)));
			  else
				 (param.nuyx2)[cursom] (im, jm) += -(float)oper / (realpart (eval (s->l2)));
			  break;
			case d_yy:
			  if (!s->path)
				 (param.nuyy2)[cursom] (im, jm) += -(float)oper * (realpart (eval (s->l2)));
			  else
				 (param.nuyy2)[cursom] (im, jm) += -(float)oper / (realpart (eval (s->l2)));
			  break;
			case symb_id:
			  if (!s->path)
				 (param.b2)[cursom] (im, jm) += (float)oper * realpart (eval (s->l2));
			  else
				 (param.b2)[cursom] (im, jm) += (float)oper / realpart (eval (s->l2));
			  break;
			default:
		     break;
		   }
	    }
       }
}

void 
femParser::sauvefct (noeudPtr s)
{
  creal          *f = NULL, *g = NULL;
  int             iglob, iloc, nloc = 1 + flag.precise * 2;
  int             nquad = flag.precise ? __mesh.getNumberOfCells() : __mesh.getNumberOfPoints();

  //  rpoint         *curpoint = __mesh.rp;

  char buffer[256];

  if (s->l2 != NULL)
    sprintf(buffer, "%s-%d", s->path, int(realpart(eval(s->l2))));
  else
    sprintf(buffer, "%s", s->path);

  if (realpart (s->value) > 0)
    saveconst (eval (s->l1), buffer);
  else
     {
       f = new creal[nquad*nloc];
       if (flag.precise)
	 g = new creal[__mesh.getNumberOfPoints()];
       for (cursloc = 0; cursloc < nquad; cursloc++)
	 for (iloc = 0; iloc < nloc; iloc++)
	    {
	      iglob = setgeom (cursloc, iloc, flag.precise);
	      f[cursom] = eval (s->l1);
	    }
       if (flag.precise)
	  {
	    for (cursloc = 0; cursloc < __mesh.getNumberOfPoints(); cursloc++)
	      g[cursloc] = __fem->P1ctoP1 (f, cursloc);
	    for (cursloc = 0; cursloc < __mesh.getNumberOfPoints(); cursloc++)
	      f[cursloc] = g[cursloc];
	  }
       if (savefct (f, __mesh.getNumberOfPoints(), buffer))
	  {
	    sprintf (errbuf, "Disk is full\n");
	    erreur (errbuf);
	  }
       delete [] f;f = NULL;
       if (flag.precise)
	 delete [] g;g =  NULL;
     }
}


void 
femParser::chargfct (noeudPtr s)
{
  int             ret;
  char buffer[256];

  if (s->l1 != NULL)
    sprintf(buffer, "%s-%d", s->path, int(realpart(eval(s->l1))));
  else
    sprintf(buffer, "%s", s->path);
  (s->name)->table = new creal[__mesh.getNumberOfPoints()];
  switch (ret = loadfct ((s->name)->table, __mesh.getNumberOfPoints(), buffer), ret)
     {
     case 2:
       if (OPTION)
	  {
	    sprintf (errbuf, "Not enough memory\n");
	    erreur (errbuf);
	  }
       (variables.ne)->value = 1;
       break;

     case 0:
       //if (OPTION)
	  {
	    sprintf (errbuf, "Can't find file: %s\n",buffer);
	    erreur (errbuf);
	  }
       (variables.ne)->value = 1;
       break;

     default:
       (variables.ne)->value = 0;
     }
}

void 
femParser::plot (noeudPtr s)
{
  int             iglob, iloc, nloc = 1 + flag.precise * 2;
  int             nquad = (flag.precise ? __mesh.getNumberOfCells() : __mesh.getNumberOfPoints());

#ifdef NOXGFEM
  /* record the name of the (new) function */
  FunctionRecord (nquad, &s->l1->name->name); 
  Function<float> *f = new Function<float>(nquad, s->l1->name->name);
#endif /* NOXGFEM */

   for (cursloc = 0; cursloc < nquad; cursloc++)
     for (iloc = 0; iloc < nloc; iloc++)
        {
          iglob = setgeom (cursloc, iloc, flag.precise);
          param.fplot[iglob] = realpart (eval (s->l1));
	 
          //(*f)[iglob] = param.fplot[iglob];
	 
#ifdef NOXGFEM
         /* copy the value  of the current function the function array */
         FunctionCopyValue (iglob, param.fplot[iglob]); 
#endif /* NOXGFEM */
        }
   
   //gfem_internals->functions.push_back(f);
  
#if defined(XGFEM)
   SlaveSendFunction(nquad, s->l1->name->name, param.fplot);
#endif /* XGFEM */
   
#ifdef  NOXGFEM
   /* update the function array */
   FunctionUpdate (); 
#endif /* NOXGFEM */

   if ( __graphic_type == FEM_GRAPHIC )
     {
       // graphics for freeFEM
       __graph->equpot(__mesh.ng, param.fplot, 20, waitm);
     }
}


void 
femParser::plot3d (noeudPtr s)
{
//   int             iglob, iloc, nloc = 1 + flag.precise * 2;
//   int             nquad = flag.precise ? __mesh.getNumberOfCells() : __mesh.getNumberOfPoints();

// #ifdef NOXGFEM
//   /* record the name of the (new) function */
//   FunctionRecord (nquad, &s->l1->name->name); 
// #endif /* NOXGFEM */
  
//   for (cursloc = 0; cursloc < nquad; cursloc++)
//     for (iloc = 0; iloc < nloc; iloc++)
//        {
//          iglob = setgeom (cursloc, iloc, flag.precise);
//          param.fplot[iglob] = realpart (eval (s->l1));
// #ifdef NOXGFEM
//          /* copy the value  of the current function the function array */
//          FunctionCopyValue (iglob, realpart (eval (s->l1))); 
// #endif /* NOXGFEM */
//        }
// #if defined(XGFEM)
//   SlaveSendFunction(nquad, s->l1->name->name, param.fplot);
// #endif /* XGFEM */

// #ifdef  NOXGFEM
//    /* update the function array */
//   FunctionUpdate (); 
// #endif /* not NOXGFEM */

//   // graphics for FreeFEM
//   graph3d (param.fplot, waitm);
}

void 
femParser::chartrig (noeudPtr s)
{
  char buffer[256];
  if (s->l1 != NULL)
    sprintf(buffer, "%s-%d", s->path, int(realpart(eval(s->l1))));
  else
    sprintf(buffer, "%s", s->path);
  switch (loadtriangulation (&__mesh, buffer))
     {
     case 1:
       if (OPTION)
	  {
	    sprintf (errbuf, "This file does not exist\n");
	    erreur (errbuf);
	  }
       (variables.ne)->value = 1;
       break;

     case 2:
       sprintf (errbuf, "Not enough memory\n");
       erreur (errbuf);
       break;

     default:
       (variables.ne)->value = 0;

       if ( __graphic_type == FEM_GRAPHIC )
	 {
	   __graph->showtriangulation (waitm);
	 }
     }
  if (flag.param)
    {
      delete __fem;
    }
  
  initparam ();
}


void 
femParser::sauvtrig (noeudPtr s)
{
  char buffer[256];
  if (s->l1 != NULL)
    sprintf(buffer, "%s-%d", s->path, int(realpart(eval(s->l1))));
  else
    sprintf(buffer, "%s", s->path);
  if (savetriangulation (&__mesh, buffer))
     {
       sprintf (errbuf, "Not enough disk space\n");
       erreur (errbuf);
     }
}

/*
 * adapt : adaptation routine
 *  - l1,l2... : functions which are used to adapt the mesh
 *
 * note : these functions must be concatened in one array
 *
 * aniso : anisotropic mesh (0 or 1)
 * c0 : error tolerance
 * lmin,lmax : min-max edge length
 * refwall, hwall : reference on the wall + h on the wall
 * Fluid_NS : Navier-Stokes equations (strong imposition)
 *********** Advanced parameters
 * nit = numbre of iterations during the adaptation process
 * securite : save mesh after each loop (0 or 1)
 * sol_interp : name of the file which contains the interpolated functions
 * angulo :angle
 * reg_ini : initial regularisation = swapping
 * reg_fin : final regularisation = filtre
 * rel_mtr : metric relaxation
 *
 */
void
femParser::adapt (noeudPtr tree)
{
#if defined(ADAPT)
  Scalar *solution, *func, *interp_func;
  Scalar **sol_interp ;

  Mallado_T0* malla;
  FemMesh* t_fin,*t_cad;
  CAD *cad;
  int s,i,j,k,nsol,nfunc;
  double *tiempo,axmin,axmax,aymin,aymax,length;
  int err;
  float h1,hmin = 1.0e10,hmax = 0.0;

  noeudPtr l1 = tree->l1,l2 = tree->l2,l3 = tree->l3,l4 = tree->l4;
  int iglob;

  /* declaration of external functions */
  extern Scalar* read_param(FemMesh*,FemMesh*,parameter&);
  extern int lecmtr (Metrica*,char*,Scalar,int);
  extern void cal_metrica(Mallado_T0*, Scalar*, Metrica*,Boolean,Scalar,Scalar,
                   Scalar,int,int,Scalar,Scalar,int);
  extern void escmtr2D(char*,Metrica*,int,Scalar);
  extern void escsol(int,int,char*,Scalar*);
  FemMesh*
    regenera_malla(Mallado_T0*,Scalar*,char*,int,int,char*,int,int,Boolean,double*,Boolean,Boolean, Boolean,Boolean,Scalar,int,CAD*,int,Scalar**, Scalar *, int, Scalar**);
  
  extern void the_clock(double *);
  
  

#if defined(DEBUG)
  cerr << "adaptation process called" << endl;
#endif

  nsol = 4;
  if (l1 == NULL) // must ha