#include <solve.hpp>

namespace ngsolve
{
  using namespace ngsolve;


  // parser for pde file


  enum TOKEN_TYPE 
    { 
      UNDEF = 0,
      NUMBER = 1, 
      STRING = 2,
      END = 3,  
  
      PLUS = '+', 
      MINUS = '-', 
      MUL = '*', DIV = '/',
      LP = '(', RP = ')', EQUAL = '=', COMMA = ',',
      LSB = '[', RSB = ']',
      COMMENT = '#',
      KW_DEFINE = 100, KW_GEOMETRY, KW_MESH,
      KW_CONSTANT, KW_VARIABLE, KW_COEFFICIENT, KW_FESPACE, KW_GRIDFUNCTION, 
      KW_BILINEARFORM, KW_LINEARFORM, KW_PRECONDITIONER, KW_BEMELEMENT,
      KW_INTEGRATOR = 200, KW_NUMPROC_ID,
      KW_NUMPROC = 300, 
    };

  struct kwstruct
  {
    TOKEN_TYPE kw; 
    char * name;
  };

  static kwstruct defkw[] =
    {
      { KW_DEFINE,      "define" },
      { KW_GEOMETRY,    "geometry" },
      { KW_MESH,        "mesh" },
      { KW_CONSTANT,    "constant" },
      { KW_VARIABLE,    "variable" },
      { KW_COEFFICIENT, "coefficient" },
      { KW_FESPACE,     "fespace" },
      { KW_GRIDFUNCTION, "gridfunction" },
      { KW_BILINEARFORM, "bilinearform" },
      { KW_BEMELEMENT,  "bemelement" },
      { KW_LINEARFORM,  "linearform" },
      { KW_PRECONDITIONER, "preconditioner" },
      { KW_NUMPROC,     "numproc" },
      { TOKEN_TYPE(0) }
    };



  class PDEScanner
  {
    TOKEN_TYPE token;
    double num_value;
    string string_value;
    SymbolTable<TOKEN_TYPE> keywords;
    SymbolTable<int> integrators;
    SymbolTable<int> numprocs;

    int linenum;

  public:
    istream * scanin;

    PDEScanner (istream * ascanin);

    TOKEN_TYPE GetToken() const
    { return token; }

    double GetNumValue() const
    { return num_value; }

    const char * GetStringValueC() const
    { return string_value.c_str(); }
    string GetStringValue() const
    { return string_value; }

    void ReadNext();

    void Error (const string & err);
  };


  PDEScanner :: PDEScanner (istream * ascanin)
  {
    scanin = ascanin;
    token = END;
    num_value = 0;
    linenum = 1;

    // initialize keywords:

    kwstruct * kwp;
    kwp = defkw;
    while (kwp->kw)
      {
	keywords.Set (kwp->name, kwp->kw);
	kwp++;
      }


    // initialize integrator and numproc keywords:

    Integrators & itgs = GetIntegrators();
    for (int i = 0; i < itgs.GetBFIs().Size(); i++)
      integrators.Set (itgs.GetBFIs()[i]->name, 1);

    for (int i = 0; i < itgs.GetLFIs().Size(); i++)
      integrators.Set (itgs.GetLFIs()[i]->name, 1);

    NumProcs & nps = GetNumProcs();
    for (int i = 0; i < nps.GetNumProcs().Size(); i++)
      numprocs.Set (nps.GetNumProcs()[i]->name, 1);
  }


  void PDEScanner :: ReadNext ()
  {
    char ch;
  

    // skip whitespaces
    do
      { 
	scanin->get(ch);

	if (ch == '\n') 
	  linenum++;

	// end of file reached
	if (scanin->eof())
	  {
	    token = END;
	    return;
	  }

	// skip comment line
	if (ch == '#')
	  {
	    while (ch != '\n')
	      {
		scanin->get(ch);
		if (scanin->eof())
		  {
		    token = END;
		    return;
		  }
	      }
	    linenum++;
	  }	
      }
    while (isspace(ch));
  
  
    switch (ch)
      {
	// case '*': case '/': 
	// case '+': 
      case '-': 
      case '(': case ')':
      case '[': case ']':
      case '=': case ',':
	{
	  token = TOKEN_TYPE (ch);
	  break;
	}
      
      default:
	{
	  if (isdigit (ch) || ch == '.')
	    {
	      scanin->putback (ch);
	      (*scanin) >> num_value;
	      token = NUMBER;
	      return;
	    }

	  (*scanin).putback (ch);
	  (*scanin) >> string_value;

	  if (integrators.Used (string_value))
	    {
	      token = KW_INTEGRATOR;
	      return;
	    }

	  if (numprocs.Used (string_value))
	    {
	      token = KW_NUMPROC_ID;
	      return;
	    }
	
	  if (keywords.Used (string_value))
	    {
	      token = keywords[string_value];
	      return;
	    }
	
	  token = STRING;
	}
      }
  }

  void PDEScanner :: Error (const string & err)
  {
    stringstream errstr;
    errstr << "Parsing error in line " << linenum << ": " << endl << err << endl;
    errstr << "input continues with <<<";
    for (int i = 0; i < 50; i++)
      {
	char ch;
	scanin->get(ch);
	errstr << ch;
	if (scanin->eof())
	  {
	    errstr << "(end of file)";
	    break;
	  }
      }
    errstr << endl << ">>> stop parsing" << endl;
    throw Exception (errstr.str());
  }
  


  void CommandList ();
  void DefineCommand ();
  void NumProcCommand ();
  void ParseFlags (Flags & flags);

  static PDEScanner * scan;
  static PDE * pde;
  


  void CommandList ()
  {
    while (scan->GetToken() != END)
      {
	switch (scan->GetToken())
	  {

	  case KW_GEOMETRY:
	    {
	      scan->ReadNext();
	      if (scan->GetToken() != '=')
		scan->Error ("Expected =");
	      scan->ReadNext();
	      cout << "Load Geometry from File " << scan->GetStringValue() << endl;
	      // Ng_LoadGeometry ((char*)scan->GetStringValue());
	      Ng_LoadGeometry (const_cast<char*> (scan->GetStringValueC()));
	      scan->ReadNext();
	      break;
	    }

	  case KW_MESH:
	    {
	      scan->ReadNext();
	      if (scan->GetToken() != '=')
		scan->Error ("Expected '='");
	      scan->ReadNext();
	      cout << "Load Mesh from File " << scan->GetStringValue() << endl;
	      Ng_LoadMesh ((char*)scan->GetStringValueC());
	      scan->ReadNext();

	      if (!pde->GetMeshAccess().GetNP())
		throw Exception ("No mesh or empty mesh file\n");
	      break;
	    }

	  case KW_DEFINE:
	    {
	      DefineCommand ();
	      break;
	    }

	  case KW_NUMPROC:
	    {
	      scan->ReadNext();
	      string npid = scan->GetStringValue();

	      scan->ReadNext();
	      string name = scan->GetStringValue();

	      scan->ReadNext();
	      Flags flags;
	      ParseFlags (flags);
	      
	      if (GetNumProcs().GetNumProc(npid))
		pde -> AddNumProc (name, GetNumProcs().GetNumProc(npid)->creator(*pde, flags));
	      else
		throw Exception (string("Undefined numproc ") + npid);
	      break;
	    }

	  default:
	    {
	      stringstream errstr;
	      errstr << "Unknown command";
	      if (scan->GetToken() == STRING)
		errstr << ": " << scan->GetStringValue() << endl;
	      else
		errstr << ", token = " << scan->GetToken() << endl;
	      scan -> Error (errstr.str());

	      scan->ReadNext();
	    }
	  }
      }
    cout << "End of file reached" << endl << endl << endl;
  }






  void DefineCommand ()
  {
    scan->ReadNext();
    switch (scan->GetToken())
      {
      case KW_CONSTANT:
	{
	  scan->ReadNext();

	  string name = scan->GetStringValue ();

	  scan->ReadNext();
	  if (scan->GetToken() != '=')
	    scan->Error ("Expected =");

	  scan->ReadNext();

	  double val;

	  if (scan->GetToken() == LP)
	    {
	      EvalFunction * fun = new EvalFunction ();
	      for (int i = 0; i < pde->GetConstantTable().Size(); i++)
		fun->DefineConstant (pde->GetConstantTable().GetName(i),
				     pde->GetConstantTable()[i]);
	      for (int i = 0; i < pde->GetVariableTable().Size(); i++)
		fun->DefineGlobalVariable (pde->GetVariableTable().GetName(i),
					   &pde->GetVariableTable()[i]);
	      
	      fun->Parse (*scan->scanin);
	      if (fun->IsConstant())
		val = fun->Eval (NULL);
	      else
		scan->Error  ("Expression not constant");
	    }
	  else if (scan->GetToken() == NUMBER)
	    val = scan->GetNumValue();
	  else if (scan->GetToken() == MINUS)
	    {
	      scan->ReadNext();
	      if (scan->GetToken() != NUMBER)
		scan -> Error ("expected a number");
	      val = -scan->GetNumValue();
	    }
	  else
	    scan->Error ("Expected a number");

	  pde->AddConstant (name, val);
	  scan->ReadNext();
	  break;
	}

      case KW_VARIABLE:
	{
	  scan->ReadNext();

	  string name = scan->GetStringValue ();

	  scan->ReadNext();
	  if (scan->GetToken() != '=')
	    scan->Error ("Expected =");

	  scan->ReadNext();

	  double val;
	  
	  if (scan->GetToken() == LP)
	    {
	      EvalFunction * fun = new EvalFunction ();
	      for (int i = 0; i < pde->GetConstantTable().Size(); i++)
		fun->DefineConstant (pde->GetConstantTable().GetName(i),
				     pde->GetConstantTable()[i]);
	      for (int i = 0; i < pde->GetVariableTable().Size(); i++)
		fun->DefineGlobalVariable (pde->GetVariableTable().GetName(i),
					   &pde->GetVariableTable()[i]);
	      
	      fun->Parse (*scan->scanin);
	      if (fun->IsConstant())
		val = fun->Eval (NULL);
	      else
		scan->Error  ("Expression not constant");
	    }
	  else if (scan->GetToken() == NUMBER)
	    val = scan->GetNumValue();
	  else if (scan->GetToken() == MINUS)
	    {
	      scan->ReadNext();
	      if (scan->GetToken() != NUMBER)
		scan -> Error ("expected a number");
	      val = -scan->GetNumValue();
	    }
	  else
	    scan->Error ("Expected a number");

	  pde->AddVariable (name, val);
	  scan->ReadNext();
	  break;
	}

      case KW_COEFFICIENT:
	{
	  ARRAY<double> dcoeffs;
	  ARRAY<EvalFunction*> coeffs;
	  ARRAY < ARRAY < ARRAY<double>* >* > polycoeffs;
	  ARRAY < ARRAY<double>* > polybounds;
	  scan->ReadNext();

	  string name = scan->GetStringValue ();
	  scan->ReadNext();

	  if (strcmp (scan->GetStringValueC(), "material") == 0)
	    {
	      // define values by material

	      scan->ReadNext();
	      SymbolTable<double> values;
	      while (scan->GetToken() == STRING)
		{
		  string mat = scan->GetStringValue();
		  scan->ReadNext();

		  double val = 1;

		  if (scan->GetToken() == LP)
		    {
		      EvalFunction * fun = new EvalFunction ();
		      for (int i = 0; i < pde->GetConstantTable().Size(); i++)
			fun->DefineConstant (pde->GetConstantTable().GetName(i),
					     pde->GetConstantTable()[i]);
		      for (int i = 0; i < pde->GetVariableTable().Size(); i++)
			fun->DefineGlobalVariable (pde->GetVariableTable().GetName(i),
						   &pde->GetVariableTable()[i]);
		      
		      fun->Parse (*scan->scanin);
		      if (fun->IsConstant())
			val = fun->Eval (NULL);
		      else
			scan->Error  ("Expression not constant");
		    }
		  else 
		    {
		      if (scan->GetToken() == MINUS)
			{
			  val = -1;
			  scan->ReadNext();
			}
		      if (scan->GetToken() != NUMBER)
			scan->Error ("number expected");
		      val *= scan->GetNumValue();
		    }
		  scan->ReadNext();

		  values.Set (mat, val);
		}

	      cout << "values = " << endl << values << endl;

	      int maxdom = -1;
	      int ne = pde->GetMeshAccess().GetNE();
	      for (int i = 0; i < ne; i++)
		maxdom = max2 (maxdom, pde->GetMeshAccess().GetElIndex(i));
	      maxdom++;

	      dcoeffs.SetSize(maxdom);
	      dcoeffs = 0;
	      for (int i = 0; i < ne; i++)
		{
		  string mat = pde->GetMeshAccess().GetElMaterial(i);
		  int index = pde->GetMeshAccess().GetElIndex(i);
		  // cout << "mat = " << mat << ", ind = " << index << endl;
		  double val;
		  if (values.Used (mat))
		    val = values[mat];
		  else if (values.Used ("default"))
		    val = values["default"];
		  else
		    throw Exception (string ("No value defined for material ")+mat);
		  
		  dcoeffs[index] = val;
		}
	      (*testout) << "material coefficients = " << endl << dcoeffs << endl;
	      pde->AddCoefficientFunction
		(name, new DomainConstantCoefficientFunction(dcoeffs));
	    }

	  else
	    
	    {
	      while (scan->GetToken() == NUMBER || 
		     scan->GetToken() == MINUS ||
		     scan->GetToken() == LP ||
		     scan->GetToken() == LSB)
		{
		  if (scan->GetToken() == LP)
		    {
		      EvalFunction * fun = new EvalFunction ();
		      for (int i = 0; i < pde->GetConstantTable().Size(); i++)
			fun->DefineConstant (pde->GetConstantTable().GetName(i),
					     pde->GetConstantTable()[i]);
		      for (int i = 0; i < pde->GetVariableTable().Size(); i++)
			fun->DefineGlobalVariable (pde->GetVariableTable().GetName(i),
						   &pde->GetVariableTable()[i]);
		      
		      fun->Parse (*scan->scanin);
		      coeffs.Append (fun);
		      if (fun->IsConstant())
			dcoeffs.Append (fun->Eval (NULL));
		      scan->ReadNext();
		    }
		  else if (scan->GetToken() == LSB) // polynomial [MW]
		    {
		      cout << "polynomial: ";
		      scan->ReadNext();
		      double val;
		      ARRAY< ARRAY<double>* > *polyco = new ARRAY< ARRAY<double>* >;
		      ARRAY<double> *polyb = new ARRAY<double>;
		      while(scan->GetToken() != RSB)
			{
			  ARRAY<double> *polyc = new ARRAY<double>;
			  while(scan->GetToken() != RSB && scan->GetToken() != LP)
			    {
			      while (scan->GetToken() == COMMA) 
				scan->ReadNext();
			      val = scan->GetNumValue();
			      //cout << "val " << val << endl;
			      if (scan->GetToken() == MINUS)
				{
				  scan->ReadNext();
				  if (scan->GetToken() != NUMBER)
				    scan -> Error ("expected a number");
				  
				  val = - scan->GetNumValue();
				}
			      polyc->Append(val);
			      scan->ReadNext();
			    }

			  cout << (*polyc)[0];
			  for(int i = 1; i<polyc->Size(); i++)
			    {
			      cout << " + " << (*polyc)[i] << "*t^"<<i;
			    }
			  

			  polyco->Append(polyc);
			  if(scan->GetToken() == LP)
			    {
			      scan->ReadNext();
			      val = scan->GetNumValue();
			      if (scan->GetToken() == MINUS)
				{
				  scan->ReadNext();
				  if (scan->GetToken() != NUMBER)
				    scan -> Error ("expected a number");
				  
				  val = - scan->GetNumValue();
				}
			      polyb->Append(val);
			      cout << " until " << val <<", then ";
			      scan->ReadNext(); // RP
			      scan->ReadNext();
			    }
			}
		      polybounds.Append(polyb);
		      polycoeffs.Append(polyco);

		      cout << endl;

		      scan->ReadNext();
		    }
		  else
		    {
		      double val = scan->GetNumValue();
		      if (scan->GetToken() == MINUS)
			{
			  scan->ReadNext();
			  if (scan->GetToken() != NUMBER)
			    scan -> Error ("expected a number");
			  
			  val = - scan->GetNumValue();
			}
		      
		      EvalFunction * fun = new EvalFunction();
		      fun->AddConstant (val);
		      coeffs.Append (fun);
		      dcoeffs.Append (val);
		      scan->ReadNext();
		    }
		  
		  if (scan->GetToken() == COMMA) 
		    scan->ReadNext();
		}
	      
	      bool allconst = 1;
	      for (int i = 0; i < coeffs.Size(); i++)
		if (!coeffs[i]->IsConstant())
		  allconst = 0;
	      
	      if (polycoeffs.Size() > 0)
		{
		  pde->AddCoefficientFunction
		    (name, new PolynomialCoefficientFunction(polycoeffs,polybounds));
		}
	      else if (allconst)
		{
		  pde->AddCoefficientFunction
		    (name, new DomainConstantCoefficientFunction(dcoeffs));
		  for (int hi = 0; hi < coeffs.Size(); hi++)
		    delete coeffs[hi];
		}
	      else
		{
		  if (pde->GetMeshAccess().GetDimension() == 2)
		    {
		      pde->AddCoefficientFunction
			(name, new DomainVariableCoefficientFunction<2>(coeffs));
		    }
		  else
		    pde->AddCoefficientFunction
		      (name, new DomainVariableCoefficientFunction<3>(coeffs));

		  for (int hi = 0; hi < coeffs.Size(); hi++)
		    delete coeffs[hi];
		}
	    }

	  break;
	}

      case KW_FESPACE:
	{
	  scan->ReadNext();

	  string name = scan->GetStringValue ();
	  scan->ReadNext();

	  Flags flags;
	  ParseFlags (flags);
	  pde->AddFESpace (name, flags);

	  break;
	}

      case KW_GRIDFUNCTION:
	{
	  scan->ReadNext();

	  string name = scan->GetStringValue ();
	  scan->ReadNext();

	  Flags flags;
	  ParseFlags (flags);
	  pde->AddGridFunction (name, flags);

	  break;
	}

      case KW_BILINEARFORM:
	{
	  scan->ReadNext();

	  string name = scan->GetStringValue ();
	  scan->ReadNext();

	  Flags flags;
	  ParseFlags (flags);
	  pde->AddBilinearForm (name, flags);
	
	  // read bilinear-form components
	  int ncoeffs;
	  while (1)
	    {
	      ncoeffs = -1;

	      TOKEN_TYPE integrator_token = scan->GetToken();

	      if (integrator_token == KW_INTEGRATOR)
		{
		  const ngfem::Integrators::IntegratorInfo * info =
		    ngfem::GetIntegrators() 
		    . GetBFI(scan->GetStringValue(), 
			     pde->GetMeshAccess().GetDimension());

		  if (info)
		    {
		      ARRAY<CoefficientFunction*> coeffs(info->numcoeffs);
		      scan->ReadNext();
		      for (int i = 0; i < info->numcoeffs; i++)
			{
			  coeffs[i] = pde->GetCoefficientFunction (scan->GetStringValue(), 1);
			  if (!coeffs[i])
			    {
			      T_GridFunction<double> * gf = 
				dynamic_cast<T_GridFunction<double>*> (pde->GetGridFunction (scan->GetStringValue(), 1));
			      if (gf)
				coeffs[i] = new GridFunctionCoefficientFunction (*gf);
			      else
				{
				  throw Exception (string("undefined coefficient ") + scan->GetStringValue());
				}
			    }

			  scan->ReadNext();
			}
		    
		      Flags partflags;
		      ParseFlags (partflags);

		      ngfem::BilinearFormIntegrator * integrator = 
			dynamic_cast<ngfem::BilinearFormIntegrator*> (info->creator(coeffs));
		      
		      if (partflags.NumFlagDefined ("order"))
			{
			  integrator -> 
			    SetIntegrationOrder
			    (int (partflags.GetNumFlag ("order",0)));
			}

		      if (partflags.NumFlagDefined ("comp"))
			{
			  if (dynamic_cast<const CompoundFESpace*> 
			      (&pde->GetBilinearForm (name)->GetFESpace()))
			    {
			      integrator = new CompoundBilinearFormIntegrator
				(*integrator, 
				 int(partflags.GetNumFlag ("comp", 1))-1);
			    }
			  else
			    {
			      integrator = new BlockBilinearFormIntegrator
				(*integrator, 
				 pde->GetBilinearForm (name)->GetFESpace().GetDimension(),
				 int(partflags.GetNumFlag ("comp", 1))-1);
			    }
			}

		      if (partflags.GetDefineFlag ("normal"))
			{
			  integrator = new NormalBilinearFormIntegrator
			    (*integrator);
			}

		      if (partflags.GetDefineFlag ("imag"))
			{
			  integrator = new ComplexBilinearFormIntegrator
			    (*integrator, Complex(0,1));
			}


		      if (partflags.NumFlagDefined ("definedon"))
			{
			  int domain = int(partflags.GetNumFlag("definedon", 0));
			  int size = max2 (pde->GetMeshAccess().GetNDomains(), 
					   pde->GetMeshAccess().GetNBoundaries());
			  BitArray definedon(size);
			  definedon.Clear();
			  definedon.Set(domain-1);
			  (*testout) << "definedon = " << definedon << endl;
			  integrator->SetDefinedOn (definedon);
			}

		      pde->AddBilinearFormIntegrator (name, integrator);
		      continue;
		    }
		  else
		    {
		      throw Exception ("Should have integrator");
		    }
		

		}
	      break;
	 
	    }
		   
	  break;
	}

      case KW_LINEARFORM:
	{
	  scan->ReadNext();
	  string name = scan->GetStringValue ();

	  scan->ReadNext();
	  Flags flags;
	  ParseFlags (flags);

	  pde->AddLinearForm (name, flags);
	
	  // read linear-form components
	  int ncoeffs;
	  while (1)
	    {
	      ncoeffs = -1;

	      TOKEN_TYPE integrator_token = scan->GetToken();

	      if (integrator_token == KW_INTEGRATOR)
		{

		  const ngfem::Integrators::IntegratorInfo * info =
		    ngfem::GetIntegrators() 
		    . GetLFI(scan->GetStringValue(), 
			     pde->GetMeshAccess().GetDimension());
		
		  if (info)
		    {
		      ARRAY<CoefficientFunction*> coeffs(info->numcoeffs);
		      scan->ReadNext();
		      for (int i = 0; i < info->numcoeffs; i++)
			{
			  coeffs[i] = pde->GetCoefficientFunction (scan->GetStringValue());
			  scan->ReadNext();
			}
		    
		      Flags partflags;
		      ParseFlags (partflags);
		    
		      ngfem::LinearFormIntegrator * integrator = 
			dynamic_cast<ngfem::LinearFormIntegrator*> (info->creator(coeffs));

		      (*testout) << "partflags = " << endl << partflags << endl;
		    
		      if (partflags.NumFlagDefined ("order"))
			{
			  integrator -> 
			    SetIntegrationOrder
			    (int (partflags.GetNumFlag ("order",0)));
			}

		      if (partflags.NumFlagDefined ("comp"))
			{
			  if (dynamic_cast<const CompoundFESpace*> 
			      (&pde->GetLinearForm (name)->GetFESpace()))
			    {
			      integrator = new CompoundLinearFormIntegrator
				(*integrator, 
				 int(partflags.GetNumFlag ("comp", 1))-1);
			    }
			  else
			    {
			      integrator = new BlockLinearFormIntegrator
				(*integrator, 
				 pde->GetLinearForm (name)->GetFESpace().GetDimension(),
				 int(partflags.GetNumFlag ("comp", 1))-1);
			    }
			}

		      if (partflags.NumFlagDefined ("normal"))
			{
			  throw Exception ("-normal flag currently not available");
			}
		    
		      if (partflags.GetDefineFlag ("imag"))
			{
			  integrator = new ComplexLinearFormIntegrator
			    (*integrator, Complex(0,1));
			} 

		      if (partflags.NumFlagDefined ("definedon"))
			{
			  int domain = int(partflags.GetNumFlag("definedon", 0));
			  int size = max2 (pde->GetMeshAccess().GetNDomains(), 
					   pde->GetMeshAccess().GetNBoundaries());
			  BitArray definedon(size);
			  definedon.Clear();
			  definedon.Set(domain-1);
			  (*testout) << "definedon = " << domain << endl;
			  integrator->SetDefinedOn (definedon);
			}
		    
		      pde->AddLinearFormIntegrator (name, integrator);
		      continue;
		    }
		  else
		    {
		      throw Exception ("Should have integrator");
		    }

		}

	      break;

	    }
		   
	  break;
	}


      case KW_PRECONDITIONER:
	{
	  scan->ReadNext();
	  string name = scan->GetStringValue ();
	
	  scan->ReadNext();
	  Flags flags;

	  ParseFlags (flags);
	  pde->AddPreconditioner (name, flags);

	  break;	
	}


      case KW_BEMELEMENT:
	{
	  scan->ReadNext();
	  string name = scan->GetStringValue ();

	  scan->ReadNext();
	  Flags flags;

	  ParseFlags (flags);
	  pde->AddBEMElement (name, flags);

	  break;
	}

      case STRING:
	{
	  scan -> Error ("Illegal define command " + string(scan->GetStringValue()));
	  break;
	}
      default:
	{
	  stringstream errstr;
	  errstr << "Illegal define token " << scan->GetStringValue();
	  scan -> Error (errstr.str());
	}
      }
  }





  void ParseFlags (Flags & flags)
  {
    while (scan->GetToken() == '-')
      {
	scan->ReadNext();
	string flag = string("-") + scan->GetStringValue();
	flags.SetCommandLineFlag (flag.c_str());
	scan->ReadNext();
      }
  }




  void PDE :: LoadPDE (const string & filename)
  {
    cout << "Load PDE from file " << filename << endl;
    pde = this;
    
    ifstream infile (filename.c_str());
    if (!infile.good())
      {
	throw Exception (string ("PDE file " + filename + " not found"));
      }


    scan = new PDEScanner (&infile);
    scan->ReadNext();
    CommandList();
    delete scan;
  }


}