femparser.cpp

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

}

void 
femParser::build (noeudPtr s)
{
  int  err, buildflag= s->l2 ? (int) realpart (eval (s->l2)) : 0;

  if ( __graphic_type == FEM_GRAPHIC )
    {
      __graph->showbdy (nbs, cr, nba, arete, hh, waitm);
    }
  
  err = __mesh.create(nbs, (long) realpart (eval (s->l1)), nba, cr, hh,
                      arete, ngbdy, sd, nbsd, &(flag.t), buildflag);
  if (err)
    switch (err)
      {
      case -1:
        erreur ("Out of memory");
      case 1:
        erreur ("Too few or too many bdy points");
      case 2:
        erreur ("Two or more points are identical");
      case 3:
        erreur ("All points are aligned");
      case 7:
        erreur ("Can't identify bdy:  internal bug");
      case 8:
        erreur ("You gave an edge which is too long");
      case 9:
        erreur ("The bdy is shaped like a 8");
      case 10:
        erreur ("One given point is in a given edge");
      case 11:
        erreur ("One subdomain is not referenced");
      case 20:
        erreur ("3 points are identical (internal bug)");
      case 21:
        erreur ("mshptg stack is too small (internal bug)");
      }
  else
    flag.t = 1;

  if ( __graphic_type == FEM_GRAPHIC )
    {
      __graph->showtriangulation(waitm);
    }

  if (flag.param)
    {
      delete __fem;
    }
  initparam ();
}

int 
femParser::setgeom (int cursloc, int iloc, int precise)
{
  int             iglob;

  if (precise)
     {
       cursom = 3 * cursloc + iloc;
       iglob = __mesh.tr[cursloc][iloc];
       (variables.x)->value = 0.999 * __mesh.rp[iglob][0] + 0.001 * (__mesh.rp[__mesh.tr[cursloc][0]][0]
		 + __mesh.rp[__mesh.tr[cursloc][1]][0] + __mesh.rp[__mesh.tr[cursloc][2]][0]) / 3;
       (variables.y)->value = 0.999 * __mesh.rp[iglob][1] + 0.001 * (__mesh.rp[__mesh.tr[cursloc][0]][1]
		 + __mesh.rp[__mesh.tr[cursloc][1]][1] + __mesh.rp[__mesh.tr[cursloc][2]][1]) / 3;
       (variables.region)->value = __mesh.ngt[cursloc];
     }
  else
     {
       iglob = cursom = cursloc;
       (variables.x)->value = __mesh.rp[iglob][0];
       (variables.y)->value = __mesh.rp[iglob][1];
       (variables.region)->value = __fem->getregion (iglob);
     }
  (variables.cursom)->value = cursom;
  (variables.ng)->value = __mesh.ng[iglob];
  (variables.nx)->value = __fem->normlx[cursom];
  (variables.ny)->value = __fem->normly[cursom];
  return iglob;
}

void 
femParser::maketable (noeudPtr s)
{
  int             i,iloc, nloc = 1 + flag.precise * 2;
  int             nquad = flag.precise ? 3 * __mesh.getNumberOfCells() : __mesh.getNumberOfPoints();
  //  femPoint         *curpoint = NULL;

  if (strcmp((s->name)->name,"x")==0)  //move femMesh
  	for(i=0;i<__mesh.getNumberOfPoints();i++) 
  	{
  		setgeom (i, 0, flag.precise);
		__mesh.rp[i][0] = realpart(eval(s->l1));  
	}
  if (strcmp((s->name)->name,"y")==0)
  	for(i=0;i<__mesh.getNumberOfPoints();i++) 
  	{
  		setgeom (i, 0, flag.precise);
		__mesh.rp[i][1] = realpart(eval(s->l1));  
	}
  
  if (!(s->name)->table)
    (s->name)->table = new creal [nquad];
  if(flag.solv>1) //local table with triangle nb in flag.solv-2
  for (iloc = 0; iloc < 3; iloc++)
		 {
			cursloc = __mesh.tr[flag.solv-2][iloc];
			setgeom (cursloc, iloc, flag.precise);
	        ((s->name)->table)[cursloc] = eval (s->l1);
       }
  else{
	  nquad = flag.precise ? __mesh.getNumberOfCells() : __mesh.getNumberOfPoints();
	  for (cursloc = 0; cursloc < nquad; cursloc++)
	    for (iloc = 0; iloc < nloc; iloc++)
			 {
				setgeom (cursloc, iloc, flag.precise);
		      ((s->name)->table)[cursom] = eval (s->l1);
	       }
       }
}


void
femParser::doconddch(int i, int cursloc,int iloc,int* ib,noeudPtr s)
{
  int             thisbdy = 0, j;
  long            im = s->junk;
  creal           aux, aux2;
  float           aux1;
  int             iglob;

	 iglob = flag.precise ? __mesh.tr[cursloc][iloc] : cursloc;
	 for (j = 0; j < i; j++)
	   thisbdy = thisbdy || (__mesh.ng[iglob] == ib[j]);
	 if (thisbdy)
	    {
	      setgeom (cursloc, iloc, flag.precise);
	      aux = eval (s->l1);
	      aux1 = norm2 (aux);
	      aux2 = (1.F + sqrtofminus1) * penal;
	      aux = (aux1 == 0) ? aux2 : aux;
	      if (flag.complexe)
		 {
		   if (N == 1)
		     (param.p1c)[cursom] = aux;
		 }
	      else
		 {
		   if (N == 1)
		     (param.p1)[cursom] = realpart (aux);
		   else if (N == 2)
		     (param.p2)[cursom][im] = realpart (aux);
		 }
	    }
} 
       
void 
femParser::conddch (noeudPtr s)
{
  long            ibv, ibb = (long) (realpart (s->value));
  int             ib[100];
  int             i=0, iloc, nloc = 1 + flag.precise * 2;
  int             nquad = flag.precise ? __mesh.getNumberOfCells() : __mesh.getNumberOfPoints();

  while (ibb > 0)
     {
       ibv = ibb / 100;
       ib[i++] = ibb - ibv * 100;
       ibb = ibv;		// decodage

     }
  if(flag.solv>1) //local table with triangle nb in flag.solv-2
  for (iloc = 0; iloc < 3; iloc++)
		{
			cursloc = __mesh.tr[flag.solv-2][iloc];
			doconddch(i,cursloc,iloc,ib,s);
		}
  else
  for (cursloc = 0; cursloc < nquad; cursloc++)
    for (iloc = 0; iloc < nloc; iloc++)
       doconddch(i,cursloc,iloc,ib,s);
}


void 
femParser::condfrr (noeudPtr s)
{
  long            ibv, ibb = (long) (realpart (s->value));
  int             thesgn = 1, thisbdy, j, i = 0, ib[100];
  long            im = s->junk;
  int             iglob, iloc, nloc = 1 + flag.precise * 2;
  int             nquad = flag.precise ? __mesh.getNumberOfCells() : __mesh.getNumberOfPoints();

  if (ibb < 0)
     {
       ibb = -ibb;
       thesgn = -1;
     }
  imdnu = im;			// used by l1 in opcondlim

  thesgndnu = thesgn;
  while (ibb > 0)
     {
       ibv = ibb / 100;
       ib[i++] = ibb - ibv * 100;
       ibb = ibv;		// decodage

     }
  if (s->l1)
    eval (s->l1);
  for (cursloc = 0; cursloc < nquad; cursloc++)
    for (iloc = 0; iloc < nloc; iloc++)
       {
	 iglob = flag.precise ? __mesh.tr[cursloc][iloc] : cursloc;
	 thisbdy = 0;
	 for (j = 0; j < i; j++)
	   thisbdy = thisbdy || (__mesh.ng[iglob] == ib[j]);
	 if (thisbdy)
	    {
	      setgeom (cursloc, iloc, flag.precise);
	      if (flag.complexe)
		 {
		   if (N == 1)
		     (param.g1c)[cursom] = (float)thesgn * eval (s->l2);
		 }
	      else
		 {
		   if (N == 1)
		     (param.g1)[cursom] = (float)thesgn * realpart (eval (s->l2));
		   if (N == 2)
		     (param.g2)[cursom][im] = (float)thesgn * realpart (eval (s->l2));
		 }
	    }
       }
}

void
femParser::edp (noeudPtr s)
{
  long            im = s->junk;
  int             iloc, nloc = 1 + flag.precise * 2;
  int             nquad = flag.precise ? __mesh.getNumberOfCells() : __mesh.getNumberOfPoints();

  eval (s->l1);
  for (cursloc = 0; cursloc < nquad; cursloc++)
    for (iloc = 0; iloc < nloc; iloc++)
       {	 setgeom (cursloc, iloc, flag.precise);
	 		if (flag.complexe)
		 	{	if (N == 1)
					(param.f1c)[cursom] = eval (s->l2);
		 	}else{
				if (N == 1)
					(param.f1)[cursom] = realpart (eval (s->l2));
				if (N == 2)
					(param.f2)[cursom][im] = realpart (eval (s->l2));
		 	}
		 }
}

void
femParser::varpde(noeudPtr s)
{
  int   i,j,k, iloc, jloc;
  creal a[36], b[6]; //N<3
 
  __fem->initvarmat ( ihowsyst, flag.complexe, N, &param);
  for(k=0;k<__mesh.getNumberOfPoints();k++)
    for(i=0;i<2*N;i++)
      systable[i]->table[k] = 0;
  for(k=0;k<__mesh.getNumberOfCells();k++)
    {
      flag.solv = 2+k; //used by all global instructions
      for(jloc=0;jloc<3;jloc++)
        for(j=0; j<N; j++)
          { 
            systable[N+j]->table[__mesh.tr[k][jloc]] = 1;
            eval(s->l2);
            b[3*j+jloc] = eval (s->l3);
            if(ihowsyst>0)
              for(iloc=0;iloc<3;iloc++)
                for(i=0;i<N;i++)
                  { 
                    systable[i]->table[__mesh.tr[k][iloc]] = 1;
                    eval(s->l2); 
                    a[18*i+9*j+iloc*3+jloc] = eval (s->l3) - b[3*j+jloc];
                    systable[i]->table[__mesh.tr[k][iloc]] = 0;
                  }
            systable[N+j]->table[__mesh.tr[k][jloc]] = 0;
          }
      __fem->assemble(ihowsyst,flag.complexe,N,k,a,b, &param);
    }
  flag.solv = 1;
  __fem->solvevarpde(N, &param, ihowsyst);
  for(k=0;k<__mesh.getNumberOfPoints();k++)
    if(N==1)
      systable[0]->table[k] = param.sol1[k];
    else 
      for(i=0;i<N;i++)
        systable[i]->table[k] = param.sol2[k][i];
}

void
femParser::solve (noeudPtr s)
{
  int             i, k;
  int             nquad = flag.precise ? 3 * __mesh.getNumberOfCells() : __mesh.getNumberOfPoints();
  float           err;

  eval (s->l1);
  eval (s->l2);
  if (s->symb == symb_solv)
    {
      err = __fem->solvePDE (&param, ihowsyst);
      if (fabs (err + 1) <= 1.0e-20)
        erreur ("Wrong matrix number or singular matrix");
    }
  else if (s->symb == sauvetout)
    {
      if (saveparam (&param, &__mesh, saveallname, N))
        {
          erreur ("Please check; disk is full or locked");
        }
    }
  
  for (cursloc = 0; cursloc < nquad; cursloc++)
    if (flag.complexe)
      {
        if (N == 1)
          {
            (param.nuyx1c)[cursloc] = 0;	/* prepare them for next PDE */
            (param.nuxx1c)[cursloc] = 0;
            (param.nuxy1c)[cursloc] = 0;
            (param.nuyy1c)[cursloc] = 0;
            (param.a11c)[cursloc] = 0;
            (param.a21c)[cursloc] = 0;
            (param.b1c)[cursloc] = 0;
            (param.c1c)[cursloc] = 0;
            (param.f1c)[cursloc] = 0;
            (param.p1c)[cursloc] = 0;
            (param.g1c)[cursloc] = 0;
            if (!flag.precise)
              systable[0]->table[cursloc] = param.sol1c[cursloc];
            else
              {
                k = cursloc / 3;
                i = cursloc - 3 * k;
                systable[0]->table[cursloc] = param.sol1c[__mesh.tr[k][i]];
              }
          }
        else if (N == 2)
          {
            /*      (param.nuyx2)[cursloc] = 0;  
                    (param.nuxx2)[cursloc] = 0;
                    (param.nuxy2)[cursloc] = 0;
                    (param.nuyy2)[cursloc] = 0;
                    (param.a12)[cursloc] = 0;
                    (param.a22)[cursloc] = 0;
                    (param.b2)[cursloc] = 0;
                    (param.c2)[cursloc] = 0;
                    (param.f2)[cursloc] = 0;
                    (param.p2)[cursloc] = 0;
                    (param.g2)[cursloc] = 0;
                    for(i=0;i<N;i++)
                    systable[i]->table[cursloc] = param.sol2[cursloc][i];
            */ 
          }
      }
    else
      {
        if (N == 1)
          {
            (param.nuyx1)[cursloc] = 0.F;	/* prepare them for next PDE */
            (param.nuxx1)[cursloc] = 0.F;
            (param.nuxy1)[cursloc] = 0.F;
            (param.nuyy1)[cursloc] = 0.F;
            (param.a11)[cursloc] = 0.F;
            (param.a21)[cursloc] = 0.F;
            (param.b1)[cursloc] = 0.F;
            (param.c1)[cursloc] = 0.F;
            (param.f1)[cursloc] = 0.F;
            (param.p1)[cursloc] = 0.F;
            (param.g1)[cursloc] = 0.F;
            if (!flag.precise)
              systable[0]->table[cursloc] = param.sol1[cursloc];
            else
              {
                k = cursloc / 3;
                i = cursloc - 3 * k;
                systable[0]->table[cursloc] = param.sol1[__mesh.tr[k][i]];
              }
          }
        else if (N == 2)
          {
            (param.nuyx2)[cursloc] = 0.F;	/* prepare them for next PDE */
            (param.nuxx2)[cursloc] = 0.F;
            (param.nuxy2)[cursloc] = 0.F;
            (param.nuyy2)[cursloc] = 0.F;
            (param.a12)[cursloc] = 0.F;
            (param.a22)[cursloc] = 0.F;
            (param.b2)[cursloc] = 0.F;
            (param.c2)[cursloc] = 0.F;
            (param.f2)[cursloc] = 0.F;
            (param.p2)[cursloc] = 0.F;
            (param.g2)[cursloc] = 0.F;
            for (i = 0; i < N; i++)
              if (!flag.precise)
                systable[i]->table[cursloc] = param.sol2[cursloc][i];
              else
                {
                  k = cursloc / 3;
                  systable[i]->table[cursloc] = param.sol2[__mesh.tr[k][cursloc - 3 * k]][i];
                }
          }
      }
  N = 1;
  N2 = 1;
  flag.syst = 0;
}

void 
femParser::opcondlim (noeudPtr s)
{
  long            jm = s->junk;
  long            im = imdnu;
  long            ibv, ibb = (long) (realpart (s->value));
  int             oper = 1, thisbdy, j, ib[100];
  int             i = 0, iglob, iloc, nloc = 1 + flag.precise * 2;
  int             nquad = flag.precise ? __mesh.getNumberOfCells() : __mesh.getNumberOfPoints();

  if (ibb < 0)
     {
       oper = -1;
       ibb = -ibb;
     }
  if (s->l1 != NULL)
    eval (s->l1);
  while (ibb > 0)
     {
       ibv = ibb / 100;
       ib[i++] = ibb - ibv * 100;
       ibb = ibv;		// decodage

     }
  for (cursloc = 0; cursloc < nquad; cursloc++)
    for (iloc = 0; iloc < nloc; iloc++)
       {
	 iglob = flag.precise ? __mesh.tr[cursloc][iloc] : cursloc;
	 thisbdy = 0;
	 for (j = 0; j < i; j++)
	   thisbdy = thisbdy || (__mesh.ng[iglob] == ib[j]);
	 if (thisbdy)
	    {
	      setgeom (cursloc, iloc, flag.precise);
	      if (flag.complexe)
		 {
		   if (N == 1)
		     if (!s->path)
		       (param.c1c)[cursom] = (float)thesgndnu * oper * eval (s->l2);
		     else
		       (param.c1c)[cursom] = (float)thesgndnu * oper / eval (s->l2);
/*    if(N==2)
   if(!s->path) (param.c2)[cursloc](im,jm)= thesgndnu * oper * eval(s->l2);
   else     (param.c2)[cursloc](im,jm)= thesgndnu * oper / eval(s->l2);
 */ 
		 }
	      else
		 {
		   if (N == 1)
		     if (!s->path)
		       (param.c1)[cursom] = thesgndnu * oper * realpart (eval (s->l2));
		     else
		       (param.c1)[cursom] = thesgndnu * oper / realpart (eval (s->l2));
		   if (N == 2)
		     if (!s->path)
		       (param.c2)[cursom] (im, jm) = thesgndnu * oper * realpart (eval (s->l2));
		     else
		       (param.c2)[cursom] (im, jm) = thesgndnu * oper / realpart (eval (s->l2));
		 }
	    }
       }
}
void 
femParser::oppde (noeudPtr s)
{
  Symbol          thesym = s->symb;
  //  rpoint         *curpoint = __mesh.rp;
  int             im, jm, oper = (long) (realpart (s->value));
  long            ofset = s->junk;
  int             iloc, nloc = 1 + flag.precise * 2;
  int             nquad = flag.precise ? __mesh.getNumberOfCells() : __mesh.getNumberOfPoints();

  im = ofset / 100;
  jm = ofset - 100 * im;
  if (oper > 0)
	 oper = 1;
  else
	 oper = -1;
  if (s->l1 != NULL)
	 eval (s->l1);
  for (cursloc = 0; cursloc < nquad; cursloc++)
	 for (iloc = 0; iloc < nloc; iloc++)
		 {
	 setgeom (cursloc, iloc, flag.precise);
	 if (flag.complexe)
	    {
	      if (N == 1)
		switch (thesym)