#include "../include/solve.hpp"


namespace ngsolve
{
  using namespace ngsolve;


  NumProc :: NumProc (PDE & apde)
    : NGS_Object (apde.GetMeshAccess(), "numproc"), pde(apde)
  {
    ;
  }
  
  NumProc :: ~NumProc()
  {
    ;
  }

  void NumProc :: Do()
  {
    ;
  }


  void NumProc :: Do (LocalHeap & lh)
  {
    Do();
  }

  void NumProc :: PrintReport (ostream & ost)
  {
    ost << "Base-class NumProc" << endl;
  }

  void NumProc :: PrintDoc (ostream & ost)
  {
    ost << "No documentation available" << endl;
  }








  /* ***************************** Numproc CalcFlux ************************** */



  ///
  class NumProcCalcFlux : public NumProc
  {
  protected:
    ///
    BilinearForm * bfa;
    ///
    GridFunction * gfu;
    ///
    GridFunction * gfflux;
    /// compute flux, not gradient
    bool applyd;
    ///
    bool useall;
    ///
    int domain;
  public:
    ///
    NumProcCalcFlux (PDE & apde, const Flags & flags);
    ///
    virtual ~NumProcCalcFlux();

    static NumProc * Create (PDE & pde, const Flags & flags)
    {
      return new NumProcCalcFlux (pde, flags);
    }

    static void PrintDoc (ostream & ost);

    ///
    virtual void Do(LocalHeap & lh);
    ///
    virtual string GetClassName () const
    {
      return "Calc Flux";
    }


    virtual void PrintReport (ostream & ost)
    {
      ost << GetClassName() << endl
	  << "Bilinear-form    = " << bfa->GetName() << endl
	  << "Differential-Op  = " << bfa->GetIntegrator(0)->Name() << endl
	  << "Gridfunction-In  = " << gfu->GetName() << endl
	  << "Gridfunction-Out = " << gfflux->GetName() << endl
	  << "apply coeffs     = " << applyd << endl
	  << "use all int'rs   = " << useall << endl;
    }
  };


  NumProcCalcFlux :: NumProcCalcFlux (PDE & apde, const Flags & flags)
    : NumProc (apde)
  {
    bfa = pde.GetBilinearForm (flags.GetStringFlag ("bilinearform", NULL));
    if (bfa->NumIntegrators()==0)
      throw Exception ("bilinearform used for CalcFlux needs at least one integrator");

    gfu = pde.GetGridFunction (flags.GetStringFlag ("solution", NULL));
    gfflux = pde.GetGridFunction (flags.GetStringFlag ("flux", NULL));
    applyd = flags.GetDefineFlag ("applyd");
    useall = flags.GetDefineFlag ("useall");
    domain = static_cast<int>(flags.GetNumFlag("domain",0))-1;
  }


  NumProcCalcFlux :: ~NumProcCalcFlux()
  {
    ;
  }

  void NumProcCalcFlux :: PrintDoc (ostream & ost)
  {
    ost << 
      "\n\nNumproc CalcFlux:\n" \
      "-----------------\n" \
      "Computes the natural flux of the bvp:\n\n"\
      "- Heat flux for thermic problems\n"\
      "- Stresses for mechanical problems\n"\
      "- Induction for magnetostatic problems\n\n"\
      "Required flags:\n" 
      "-bilinearform=<bfname>\n" 
      "    the first integrator for the bf computes the flux\n" \
      "-gridfunction=<gfname>\n" \
      "    grid-function providing the primal solution field\n" \
      "-flux=<gfname>\n" \
      "    grid-function used for storing the flux (e.g., vector-valued L2)\n\n" \
      "\nOptional flags:\n" \
      "-applyd   apply coefficient matrix (compute either strains or stresses, B-field or H-field,..\n"\
      "-useall   use all integrators for computing the flux, and add up result\n"\
	<< endl;
  }


  void NumProcCalcFlux :: Do(LocalHeap & lh)
  {
    cout << "Num-proc calc flux" << endl;

    if (!useall)
      {
	if (!bfa->GetFESpace().IsComplex())
	  CalcFluxProject (pde.GetMeshAccess(),
			   dynamic_cast<const S_GridFunction<double>&> (*gfu), 
			   dynamic_cast<S_GridFunction<double>&> (*gfflux), 
			   *bfa->GetIntegrator(0),
			   applyd, domain, lh);
	else
	  CalcFluxProject (pde.GetMeshAccess(),
			   dynamic_cast<const S_GridFunction<Complex>&> (*gfu), 
			   dynamic_cast<S_GridFunction<Complex>&> (*gfflux), 
			   *bfa->GetIntegrator(0),
			   applyd, domain, lh);
      }
    else
      {
	gfflux->GetVector() = 0;
	for (int k = 0; k < bfa->NumIntegrators(); k++)
	  {
	    if (!bfa->GetFESpace().IsComplex())
	      CalcFlux (pde.GetMeshAccess(),
			dynamic_cast<const S_GridFunction<double>&> (*gfu), 
			dynamic_cast<S_GridFunction<double>&> (*gfflux), 
			*bfa->GetIntegrator(k),
			applyd, 1, domain);
	    else
	      CalcFlux (pde.GetMeshAccess(),
			dynamic_cast<const S_GridFunction<Complex>&> (*gfu), 
			dynamic_cast<S_GridFunction<Complex>&> (*gfflux), 
			*bfa->GetIntegrator(k),
			applyd, 1, domain);
	  }
      }
  }



  /* **************************** Numproc Draw Flux ********************************* */


  ///
  class NumProcDrawFlux : public NumProc
  {
  protected:
    netgen::SolutionData * vis;
    ///
    BilinearForm * bfa;
    ///
    GridFunction * gfu;
    /// compute flux, not gradient
    bool applyd;

    BilinearFormIntegrator * bfi2d;
    BilinearFormIntegrator * bfi3d;
    string label;

  public:
    ///
    NumProcDrawFlux (PDE & apde, const Flags & flags);
    ///
    virtual ~NumProcDrawFlux();

    static NumProc * Create (PDE & pde, const Flags & flags)
    {
      return new NumProcDrawFlux (pde, flags);
    }
    ///
    virtual void Do();

    static void PrintDoc (ostream & ost);
    ///
    virtual string GetClassName () const
    {
      return "Draw Flux";
    }


    virtual void PrintReport (ostream & ost)
    {
      ost << GetClassName() << endl;
      if (bfa) ost << "Bilinear-form    = " << bfa->GetName() << endl;
      if (bfa) ost << "Differential-Op  = " << bfa->GetIntegrator(0)->Name() << endl;
      if (gfu) ost << "Gridfunction-In  = " << gfu->GetName() << endl;
      ost << "apply coeffs     = " << applyd << endl;
    }
  };


  NumProcDrawFlux :: NumProcDrawFlux (PDE & apde, const Flags & flags)
    : NumProc (apde)
  {
    bfa = 0;
    gfu = 0;
    
    if (flags.GetDefineFlag ("order"))
      {
	Ng_SolutionData soldata;
	Ng_InitSolutionData (&soldata);
	soldata.name = "order";
	soldata.soltype = NG_SOLUTION_ELEMENT_ORDER;
	Ng_SetSolutionData (&soldata);
	return;
      }
    if (flags.GetDefineFlag ("marked"))
      {
	Ng_SolutionData soldata;
	Ng_InitSolutionData (&soldata);
	soldata.name = "marked";
	soldata.soltype = NG_SOLUTION_MARKED_ELEMENTS;
	Ng_SetSolutionData (&soldata);
	return;
      }


    bfa = pde.GetBilinearForm (flags.GetStringFlag ("bilinearform", NULL));
    gfu = pde.GetGridFunction (flags.GetStringFlag ("solution", NULL));
    applyd = flags.GetDefineFlag ("applyd");
    label = flags.GetStringFlag ("label", "");

    bfi3d = 0;
    bfi2d = 0;
    for (int i = 0; i < bfa->NumIntegrators(); i++)
      {
	if (!bfi3d && bfa->GetIntegrator(i)->DimElement() == 3)
	  bfi3d = bfa->GetIntegrator(i);
	if (!bfi2d && bfa->GetIntegrator(i)->DimElement() == 2)
	  bfi2d = bfa->GetIntegrator(i);
      }

    /*
      if (bfi2d) cout << "bfi2d = " 
      << bfi2d->Name()
      << ", dim = " << bfi2d->DimFlux()
      << endl;
      if (bfi3d) cout << "bfi3d = " << bfi3d->Name()
      << ", dim = " << bfi3d->DimFlux()
      << endl;
    */

    if (!bfa->GetFESpace().IsComplex())
      vis = new VisualizeGridFunction<double> (ma, gfu, bfi2d, bfi3d, applyd);
    else
      vis = new VisualizeGridFunction<Complex> (ma, gfu, bfi2d, bfi3d, applyd);


    Ng_SolutionData soldata;
    Ng_InitSolutionData (&soldata);
  
    // soldata.name = const_cast<char*> (gfu->GetName().c_str());
    soldata.name = (char*)label.c_str();
    soldata.data = 0;
    soldata.components = vis->GetComponents();
    soldata.iscomplex = vis->IsComplex();
    soldata.draw_surface = bfi2d != 0;
    soldata.draw_volume  = bfi3d != 0;
    soldata.dist = 1;
    soldata.soltype = NG_SOLUTION_VIRTUAL_FUNCTION;
    soldata.solclass = vis;
    Ng_SetSolutionData (&soldata);
  }

  NumProcDrawFlux :: ~NumProcDrawFlux()
  {
    ;
  }

  void NumProcDrawFlux :: Do()
  {
    cout << "Num-proc draw flux" << endl;
  }










  void NumProcDrawFlux :: PrintDoc (ostream & ost)

  {
    ost <<
      "\n\nNumproc DrawFlux:\n" \
      "-----------------\n" \
      "Adds the natural flux to the visualization dialogbox:\n"\
      "It takes the first integrator of the bilinear-form\n" \
      "- Heat flux for thermic problems\n"\
      "- Stresses for mechanical problems\n"\
      "- Induction for magnetostatic problems\n\n"\
      "Required flags:\n"
      "-bilinearform=<bfname>\n" \
      "    the first integrator for the bf computes the flux\n" \
      "-solution=<gfname>\n" \
      "    grid-function providing the primal solution field\n" \
      "\nOptional flags:\n" \
      "-applyd\n" \
      "    apply coefficient matrix (compute either strains or stresses, B-field or H-field,..\n"\
      "-label=<name>\n" \
      "    label printed in the visualization dialogbox\n\n";

  }






  ///
  class NumProcEvaluate : public NumProc
  {
  protected:
    ///
    BilinearForm * bfa;
    ///
    LinearForm * lff;
    ///
    GridFunction * gfu;
    ///
    GridFunction * gfv;
    ///
    Vector<double> point;
    ///
    Vector<double> point2;
    ///
    string filename, text;
    ///
    string variablename;
    ///
    bool applyd;
    ///
    bool hermitsch;
  public:
    ///
    NumProcEvaluate (PDE & apde, const Flags & flags);
    ///
    virtual ~NumProcEvaluate();

    static NumProc * Create (PDE & pde, const Flags & flags)
    {
      return new NumProcEvaluate (pde, flags);
    }
    static void PrintDoc (ostream & ost);

    ///
    virtual void Do(LocalHeap & lh);
    ///
    virtual void PrintReport (ostream & ost);

    virtual string GetClassName () const
    {
      return "Evaluate";
    }

  };





  NumProcEvaluate :: NumProcEvaluate (PDE & apde, const Flags & flags)
    : NumProc (apde), point(1), point2(1)
  {
    bfa = pde.GetBilinearForm (flags.GetStringFlag ("bilinearform", ""), 1); 
    lff = pde.GetLinearForm (flags.GetStringFlag ("linearform", ""), 1);
    gfu = pde.GetGridFunction (flags.GetStringFlag ("gridfunction", ""), 0); 
    gfv = pde.GetGridFunction (flags.GetStringFlag ("gridfunction2", ""), 1); 

    variablename = flags.GetStringFlag ("resultvariable", "");

    if (flags.NumListFlagDefined ("point"))
      {
	const ARRAY<double> & p = flags.GetNumListFlag ("point");
	point.SetSize(p.Size());
	for (int i = 0; i < p.Size(); i++)
	  point(i) = p[i];
	//      cout << "point = " << point << endl;
      }
    if (flags.NumListFlagDefined ("point2"))
      {
	const ARRAY<double> & p = flags.GetNumListFlag ("point2");
	point2.SetSize(p.Size());
	for (int i = 0; i < p.Size(); i++)
	  point2(i) = p[i];
	//      cout << "point2 = " << point2 << endl;
      }
    text = flags.GetStringFlag ("text","value");
    filename = flags.GetStringFlag ("filename","");
    applyd = flags.GetDefineFlag ("applyd");
    hermitsch = flags.GetDefineFlag ("hermitsch");
  }

  NumProcEvaluate :: ~NumProcEvaluate()
  {
    ;
  }



  void NumProcEvaluate :: PrintDoc (ostream & ost)

  {
    ost <<
      "\n\nNumproc Evaluate:\n" \
      "-----------------\n" \
      "Evaluates linearform or bilinearform, or pointvalues:\n"\
      "Required flags:\n" \
      "-gridfunction=<gfname>\n" \
      "    gridfunction to evaluate\n" \
      "\nOptional flags:\n" \
      "-bilinearform=<bfname>\n" \
      "    evaluates bilinear-form bfname(gfname,gfname2)\n"\
      "    needs second gridfunction <gfname2>\n" \
      "-gridfunction2=<gfname2>\n" \
      "    second gridfunciton for bilinear-form evaluation\n" \
      "-linearform=<lfname>\n" \
      "    evaluates linearform lfname(gfname)\n" \
      "-point=[x,y,z]\n" \
      "    evaluates diffop applied to gridfunction in point p\n" \
      "    diffop taken from first term in bilinear-form\n" \
      "-point2=[x2,y2,z2]\n" \
      "    evaluates diffop applied to gridfunction in line p1-p2, and writes to file\n" \
      "-applyd\n" \
      "    evaluates flux instead of derivatives\n" \
      "-text=blabla \n" \
      "    prints text \n" \
      " -resultvariabe=<varname> \n" \
      "    stores scalar, real results in variable <varname>\n" \
      "-filename=<fn> \n" \
      "    writes values in file \n";
  }




  void NumProcEvaluate :: Do(LocalHeap & lh)
  {
    int i, j;
    double result(0);

    if (lff)
      {
	if (!gfu)
	  throw Exception ("evaluate linear-form needs an argument -gridfunction=u");

	cout << "<" << lff->GetName() << ", " << gfu->GetName() << "> = " << flush;
	if (!lff->GetFESpace().IsComplex())
	  {
	    result = S_InnerProduct<double>(lff->GetVector(), gfu->GetVector());
	    cout << result << endl;
	  }
	else
	  cout << S_InnerProduct<Complex>(lff->GetVector(), gfu->GetVector()) << endl;
      }
    else if (point.Size() >= 2)
      {
	const BilinearFormIntegrator & bfi = *bfa->GetIntegrator(0);

	if (point2.Size() >= 2)
	  {
	    // plot values along line p1-p2

	    ofstream ofile (filename.c_str());

	    for (i = 0; i <= 1000; i++)
	      {
		lh.CleanUp();
		FlatVector<double> p(point.Size(), lh);
		p = point + (double(i)/1000) * (point2-point);
		if (!gfu->GetFESpace().IsComplex())
		  {
		    FlatVector<double> pflux;
		    CalcPointFlux (ma, *gfu, p,
				   pflux, bfi, applyd, lh);
		  
		    for (j = 0; j < p.Size(); j++)
		      ofile << p(j) << " ";
		    for (j = 0; j < pflux.Size(); j++)
		      ofile << pflux(j) << " ";
		    ofile << endl;
		  }
		else
		  {
		    FlatVector<Complex> pflux;
		    CalcPointFlux (ma, *gfu, p,
				   pflux, bfi, applyd, lh);
		    for (j = 0; j < p.Size(); j++)
		      ofile << p(j) << " ";
		    for (j = 0; j < pflux.Size(); j++)
		      ofile << pflux(j).real() << " " << pflux(j).imag() << " ";
		    ofile << endl;
		  }
	      }
	  }
	else
	  {
	    cout << text << " (" << point(0);
	    for (int i = 1; i < point.Size(); i++)
	      cout << "," << point(i);
	    cout << ") = " << flush;

	    if (!gfu->GetFESpace().IsComplex())
	      {
		FlatVector<double> pflux;
		CalcPointFlux (ma, *gfu, point,
			       pflux, bfi, applyd, lh);
		for (int i = 0; i < pflux.Size(); i++)
		  cout << pflux(i) << "   ";
		cout << endl;
	      }
	    else
	      {
		FlatVector<Complex> pflux;
		CalcPointFlux (ma, *gfu, point,
			       pflux, bfi, applyd, lh);
		for (int i = 0; i < pflux.Size(); i++)
		  cout << pflux(i) << "   ";
		cout << endl;
	      }
	  }
      }
    else if (bfa)
      {
	if (!gfu)
	  throw Exception ("evaluate bilinear-form needs an argument -gridfunction=u");
	if (!gfv)
	  throw Exception ("evaluate bilinear-form needs an argument -gridfunction2=v");

	cout << bfa->GetName() << "(" << gfu->GetName() << ", " << gfv->GetName() << ") = " << flush;
	BaseVector & vecu = gfu->GetVector();
	BaseVector & vecv = gfv->GetVector();
	BaseVector & hv = *vecu.CreateVector();
	bfa->GetMatrix().Mult (vecu, hv);
	if (!bfa->GetFESpace().IsComplex())
	  {
	    result = S_InnerProduct<double>(vecv, hv);
	    cout << setprecision(16) << result << endl;
	    cout << "err = " << setprecision(16) << sqrt(fabs((result-0.2140758026709826))) << endl;
	  }
	else
	  {
	    if (!hermitsch)
	      cout << S_InnerProduct<Complex>(vecv, hv) << endl;
	    else
	      {
		FlatVector<Complex> fvecv = vecv.FVComplex();
		FlatVector<Complex> fhv = hv.FVComplex();
		Complex sum = 0.0;
		for (int i = 0; i < fvecv.Size(); i++)
		  sum += fvecv(i) * conj (fhv(i));
		cout << sum << endl;
	      }
	  }
	delete &hv;
      }
    
    pde.GetVariable(variablename,true) = result;
  }

  void NumProcEvaluate :: PrintReport (ostream & ost)
  {
    ost << "NumProcEvaluate:" << endl;
  }





  /////////////////////////////////////////////////////////////////////////////
  ///
  class NumProcAnalyze : public NumProc
  {
  protected:
    GridFunction * gfu;
    ///
    string variablename;
    ///
    bool nodistinction;
    ///
    bool volanalyze;
    ///
    bool surfanalyze;
    ///
    int component;
    ///
    ARRAY<int> surfdomains;
    ///
    ARRAY<int> voldomains;
  public:
    ///
    NumProcAnalyze (PDE & apde, const Flags & flags);
    ///
    virtual ~NumProcAnalyze();

    static NumProc * Create (PDE & pde, const Flags & flags)
    {
      return new NumProcAnalyze (pde, flags);
    }
    static void PrintDoc (ostream & ost);

    ///
    virtual void Do(LocalHeap & lh);
    ///
    virtual void PrintReport (ostream & ost);

    virtual string GetClassName () const
    {
      return "Analyze";
    }

  };





  NumProcAnalyze :: NumProcAnalyze (PDE & apde, const Flags & flags)
    : NumProc (apde)
  {
    gfu = pde.GetGridFunction (flags.GetStringFlag ("gridfunction", ""), 0); 

    variablename = flags.GetStringFlag ("resultvariable", "");


    volanalyze = flags.GetDefineFlag("volume");
    surfanalyze = flags.GetDefineFlag("surface");
    nodistinction = flags.GetDefineFlag("nodistinction");

    if(!volanalyze && !surfanalyze) volanalyze = true;

    component = static_cast<int>(flags.GetNumFlag("comp",0)) - 1;

    if(flags.NumListFlagDefined("voldomains"))
      {
	voldomains.SetSize(flags.GetNumListFlag("voldomains").Size());
	for(int i=0; i<voldomains.Size(); i++)
	  {
	    voldomains[i] = static_cast<int>(flags.GetNumListFlag("voldomains")[i]);
	  }
      }
    
    if(flags.NumListFlagDefined("surfdomains"))
      {
	surfdomains.SetSize(flags.GetNumListFlag("surfdomains").Size());
	for(int i=0; i<surfdomains.Size(); i++)
	  {
	    surfdomains[i] = static_cast<int>(flags.GetNumListFlag("surfdomains")[i]);
	  }
      }
	    
  }

  NumProcAnalyze :: ~NumProcAnalyze()
  {
    ;
  }



  void NumProcAnalyze :: PrintDoc (ostream & ost)

  {
    ost <<
      "\n\nNumproc Analyze:\n" \
      "--------------------------\n" \
      "\nAnalyzes a gridfunction, i.e. calculates maximum, minimum, and average value\n"\
      "Required flags:\n" \
      "-gridfunction=<gfname>\n" \
      "    gridfunction to analyze\n" \
      "\nOptional flags:\n" \
      "-volume\n"\
      "    analyze inside volumes\n" \
      "-surface\n" \
      "    analyze on surfaces\n" \
      "-comp=<component>\n" \
      "    analyze only one given component\n" \
      "-voldomains=<domainslist>\n" \
      "    analyze only inside given domains\n" \
      "-surfdomains=<domainslist>\n" \
      "    analyze only on given surfaces\n" \
      "-resultvariable=<basename>\n" \
      "    the results are saved to the variables basename.<vol|surf><number>.comp<number>.<min|max|av>\n" \
      "-nodistinction\n" \
      "    if set then there is no distinction of the different components and domains. Results are saved to basename.<min|max|av>\n";
  }




  void NumProcAnalyze :: Do(LocalHeap & lh)
  {
    int i, j, dom;

    bool writevar = (strcmp(variablename.c_str(), "") != 0);
    string actvarname;

    const BilinearFormIntegrator * Evaluator_ptr;
    const BilinearFormIntegrator * BoundaryEvaluator_ptr;
    

    const FESpace & fes = gfu->GetFESpace();

    const int components = fes.GetEvaluator()->DimFlux();
    int ndomains;
    
    string typestring;

    for(int count=0; count < 2; count++)
      {
	if(count == 0)
	  {
	    if(!volanalyze) continue;

	    typestring = ".vol";

	    Evaluator_ptr = fes.GetEvaluator();
	    BoundaryEvaluator_ptr = NULL;
	    ndomains = pde.GetMeshAccess().GetNDomains();

	  }
	if(count == 1)
	  {
	    if(volanalyze) continue;

	    typestring = ".surf";

	    Evaluator_ptr = NULL;
	    BoundaryEvaluator_ptr = fes.GetBoundaryEvaluator();
	    ndomains = pde.GetMeshAccess().GetNBoundaries();
	  }

	ARRAY<int> & domains = ((count == 0) ? voldomains : surfdomains);


	if(domains.Size() == 0)
	  {
	    domains.SetSize(ndomains);
	    for(int i=0; i<domains.Size(); i++)
	      {
		domains[i] = i;
	      }
	  }

	const FESpace & fes = gfu->GetFESpace();

	VisualizeGridFunction<double> vgfu(pde.GetMeshAccess(),gfu,
					   BoundaryEvaluator_ptr,
					   Evaluator_ptr,false);


	ARRAY<double> mini, maxi, average, vol;

	if(component == -1)
	  {
	    mini.SetSize(components*ndomains);
	    maxi.SetSize(components*ndomains);
	    average.SetSize(components*ndomains);
	  }
	else
	  {
	    mini.SetSize(ndomains);
	    maxi.SetSize(ndomains);
	    average.SetSize(ndomains);
	  }
	  

	if(nodistinction)
	  {
	    vol.SetSize(ndomains);
	    vgfu.Analyze(mini,maxi,average,vol,component);
	  }
	else
	  {
	    vgfu.Analyze(mini,maxi,average,component);
	  }	

	cout << endl << gfu->GetName();

	if(nodistinction)
	  {
	    if (count == 0) cout << " (volume)" << endl;
	    else cout << " (surface)" << endl;

	    double gmin(1e100), gmax(-1e100), gav(0), gvol(0);

	    for(int i=0; i<domains.Size(); i++)
	      {
		dom = domains[i];
		gvol += vol[dom];
		
		if(component == -1)
		  {
		    for(int j=0; j<components; j++)
		      {
			if(mini[dom*components*j] < gmin) gmin = mini[dom*components*j];
			if(maxi[dom*components*j] > gmax) gmax = maxi[dom*components*j];
			gav += average[dom*components*j];
		      }
		  }
		else
		  {
		    if(mini[dom] < gmin) gmin = mini[dom];
		    if(maxi[dom] > gmax) gmax = maxi[dom];
		    gav += average[dom];
		  }
	      }
	    gav /= gvol;


	    cout << "min: " << gmin << endl
		 << "max: " << gmax << endl
		 << "av:  " << gav << endl;
	    if(writevar)
	      {
		ostringstream avnmin,avnmax,avnav;

		avnmin << variablename << ".min";
		actvarname = avnmin.str();
		if(!pde.VariableUsed(actvarname)) pde.AddVariable(actvarname,gmin);
		else pde.GetVariable(actvarname) = gmin;

		avnmax << variablename << ".max";
		actvarname = avnmax.str();
		if(!pde.VariableUsed(actvarname)) pde.AddVariable(actvarname,gmax);
		else pde.GetVariable(actvarname) = gmax;

		avnav << variablename << ".av";
		actvarname = avnav.str();
		if(!pde.VariableUsed(actvarname)) pde.AddVariable(actvarname,gav);
		else pde.GetVariable(actvarname) = gav;
	      }
	  }
	else
	  {
	    if(component != -1) cout << ", component " << component+1;
	    cout << endl;
	    for(int i=0; i<domains.Size(); i++)
	      {
		dom = domains[i];
		
		if(count == 0) cout << "on domain " << dom<<": " << endl;
		else cout << "on surfacedomain " << dom<<": " << endl;
		
		if(component == -1)
		  {
		    cout << "min:" << endl;
		    for(int j=0; j<components; j++)
		      {
			cout << mini[dom*components+j] << " ";
			if(writevar)
			  {
			    ostringstream avn;
			    avn << variablename << typestring << dom << ".comp" << j+1 << ".min";
			    actvarname = avn.str();
			    if(!pde.VariableUsed(actvarname)) pde.AddVariable(actvarname,mini[dom*components+j]);
			    else pde.GetVariable(actvarname) = mini[dom*components+j];
			  }
		      }
		    cout << endl << "max:" << endl;
		    for(int j=0; j<components; j++)
		      {
			cout << maxi[dom*components+j] << " ";
			if(writevar)
			  {
			    ostringstream avn;
			    avn << variablename << typestring << dom << ".comp" << j+1 << ".max";
			    actvarname = avn.str();
			    if(!pde.VariableUsed(actvarname)) pde.AddVariable(actvarname,maxi[dom*components+j]);
			    else pde.GetVariable(actvarname) = maxi[dom*components+j];
			  }
		      }
		    cout << endl << "av:" << endl;
		    for(int j=0; j<components; j++)
		      {
			cout << average[dom*components+j] << " ";
			if(writevar)
			  {
			    ostringstream avn;
			    avn << variablename << typestring << dom << ".comp" << j+1 << ".av";
			    actvarname = avn.str();
			    if(!pde.VariableUsed(actvarname)) pde.AddVariable(actvarname,average[dom*components+j]);
			    else pde.GetVariable(actvarname) = average[dom*components+j];
			  }
		      }
		    cout << endl;
		  }
		else
		  {
		    cout << "min: " << mini[dom] << " max: " << maxi[dom] << " av: " << average[dom] << endl;
		    if(writevar)
		      {
			ostringstream avn;
			avn << variablename << typestring << dom << ".comp" << component+1;
			actvarname = avn.str()+".min";
			if(!pde.VariableUsed(actvarname)) pde.AddVariable(actvarname,mini[dom]);
			else pde.GetVariable(actvarname) = mini[dom];
			actvarname = avn.str()+".max";
			if(!pde.VariableUsed(actvarname)) pde.AddVariable(actvarname,maxi[dom]);
			else pde.GetVariable(actvarname) = maxi[dom];
			actvarname = avn.str()+".av";
			if(!pde.VariableUsed(actvarname)) pde.AddVariable(actvarname,average[dom]);
			else pde.GetVariable(actvarname) = average[dom];
		      }
		  }
	      }
	  }
      }
  }

  void NumProcAnalyze :: PrintReport (ostream & ost)
  {
    ost << "NumProcAnalyze:" << endl;
  }









  //////////////////////////////////////////////////////////////////////////////


  class NumProcWarn : public NumProc
  {
  protected:
    ///
    string variablename1,variablename2;
    ///
    double val1,val2;
    ///
    bool less,lessorequal,greater,greaterorequal;
    ///
    string text;
    
  public:
    ///
    NumProcWarn (PDE & apde, const Flags & flags);
    ///
    virtual ~NumProcWarn();

    static NumProc * Create (PDE & pde, const Flags & flags)
    {
      return new NumProcWarn (pde, flags);
    }

    static void PrintDoc (ostream & ost);

    ///
    virtual void Do(LocalHeap & lh);
    ///
    virtual void PrintReport (ostream & ost);

    virtual string GetClassName () const
    {
      return "Warn";
    }
  };


  
  NumProcWarn :: NumProcWarn (PDE & apde, const Flags & flags)
    : NumProc(apde)
  {
    text = flags.GetStringFlag("text","");

    variablename1 = flags.GetStringFlag("var1","");
    variablename2 = flags.GetStringFlag("var2","");

    val1 = flags.GetNumFlag("val1",0);
    val2 = flags.GetNumFlag("val2",0);

    less = flags.GetDefineFlag("less");
    lessorequal = flags.GetDefineFlag("lessorequal");
    greater = flags.GetDefineFlag("greater");
    greaterorequal = flags.GetDefineFlag("greaterorequal");

    if(!less && !lessorequal && !greater && !greaterorequal) less = true;
  }

  NumProcWarn :: ~NumProcWarn(){;}

  void NumProcWarn :: PrintDoc (ostream & ost)
  {
    ost <<
      "\n\nNumproc Warn:\n" \
      "--------------------------\n" \
      "Checks if variables fulfill a given inequality\n"\
      "\nFlags:\n"\
      " -var1=<variablename> or -val1=<value>\n"\
      " -var2=<variablename> or -val2=<value>\n"\
      " -less or -lessorequal or -greater or -greaterorequal\n"\
      "\nOptional flags:\n"\
      " -text=<text>\n"\
      "     prints the text \n"\
      "\nIf var1 resp. val1 is less, lessorequal, greater, greaterorequal than var2 resp. val2, "\
      "then a warning is given.\n";
  }

  void NumProcWarn :: Do(LocalHeap & lh)
  {
    double value1,value2;
    ostringstream warnleft, warnright;
    string warnmiddle;

    if(strcmp(variablename1.c_str(),"") != 0)
      { 
	value1 = pde.GetVariable(variablename1);
	warnleft << variablename1 << " ("<< value1 <<")";
      }
    else 
      {
	value1 = val1;
	warnleft << value1;
      }
    if(strcmp(variablename2.c_str(),"") != 0)
      { 
	value2 = pde.GetVariable(variablename2);
	warnright << variablename2 << " ("<< value2 <<")";
      }
    else 
      {
	value2 = val2;
	warnright << value2;
      }

    bool warn;

    if(less)
      {
	warn = (value1 < value2);
	warnmiddle = " < ";
      }
    else if(lessorequal)
      {
	warn = (value1 <= value2);
	warnmiddle = " <= ";
      }
    else if(greater)
      {      
	warn = (value1 > value2);
	warnmiddle = " > ";
      }
    else if(greaterorequal)
      {
	warn = (value1 >= value2);
	warnmiddle = " >= ";
      }

    if(warn)
      {
	cout << "Warning: " << text << endl << warnleft.str() << warnmiddle << warnright.str() << endl;
	ostringstream tclstring;
	
	tclstring << "set w .warning;"
		  << "toplevel $w;" 
		  << "message $w.mes -aspect 2000 -text \""
		  << "WARNING:\n" << text << "\n" << warnleft.str() << warnmiddle << warnright.str()<< "\n"
		  << "\"" << endl
		  << "button $w.done -text \"Done\" -command \"destroy $w\"" << endl
		  << "pack $w.mes" << endl
		  << "pack $w.done" << endl
		  << "wm withdraw $w\n"
	  //<< "wm geom $w +200+100;"
		  << "wm deiconify $w;"
		  << "wm title $w \"Warning\"\n";

	char *dummy; dummy = new char[tclstring.str().size()+1];
	strcpy(dummy,tclstring.str().c_str());

	Tcl_Eval(pde.tcl_interpreter,dummy);

	delete [] dummy;
		 
      }
  }


  void NumProcWarn :: PrintReport (ostream & ost)
  {
    ost << "NumProcWarn:" << endl;
  }


  //////////////////////////////////////////////////////////////////////////////


  class NumProcTclMenu : public NumProc
  {
  protected:


  public:
    ///
    NumProcTclMenu (PDE & apde, const Flags & flags);
    ///
    virtual ~NumProcTclMenu();

    static NumProc * Create (PDE & pde, const Flags & flags)
    {
      return new NumProcTclMenu (pde, flags);
    }

    static void PrintDoc (ostream & ost);

    ///
    virtual void Do(LocalHeap & lh);
    ///
    virtual void PrintReport (ostream & ost);

    virtual string GetClassName () const
    {
      return "TclMenu";
    }
  };


  
  NumProcTclMenu :: NumProcTclMenu (PDE & apde, const Flags & flags)
    : NumProc(apde)
  {

    bool newmenu = flags.GetDefineFlag("newmenu");

    string menuname (flags.GetStringFlag("menuname",""));

    string text (flags.GetStringFlag("text",""));


    ARRAY<double> centerpoint;
    bool center = flags.NumListFlagDefined("centerpoint");
    if(center) centerpoint = flags.GetNumListFlag("centerpoint");

    
    ARRAY<double> rotation_pars; // (alpha, v1, v2, v3)
    bool rotation = flags.NumListFlagDefined("rotation");
    if(rotation) rotation_pars = flags.GetNumListFlag("rotation");
    

    ARRAY<double> clipvec;
    bool clip = flags.NumListFlagDefined("clipvec");
    if(clip) clipvec = flags.GetNumListFlag("clipvec");

    string fieldname (flags.GetStringFlag("fieldname",""));
    int component = static_cast<int>(flags.GetNumFlag("comp",1));
    string evaluate (flags.GetStringFlag("evaluate",""));
    if(evaluate != "") component = 0;

    double deformationscale = flags.GetNumFlag("deformationscale",0);
    bool deformationoff = ( flags.NumFlagDefined("deformationscale") && fabs(deformationscale) < 1.e-6 );
    bool deformationon = ( !deformationoff && flags.NumFlagDefined("deformationscale") );

    double light = flags.GetNumFlag("light",-1);
    if(light > 1) light = 1;

    double minval, maxval;
    bool autoscale = flags.GetDefineFlag("autoscale");
    bool noautoscale = ( flags.NumFlagDefined("minval") && flags.NumFlagDefined("maxval") );
    if(noautoscale)
      {
	minval = flags.GetNumFlag("minval",0);
	maxval = flags.GetNumFlag("maxval",0);
      }

    

    // Build menus

    
    ostringstream tclstring;

    bool ng_setvispar(false), ng_vissetpar(false);

    if(newmenu)
      {
	tclstring << ".ngmenu add cascade -label \"" << text 
		  <<"\" -menu .ngmenu." << menuname <<" -underline 0\n"
		  << "menu .ngmenu." << menuname << endl;
      }
    else
      {
	tclstring << ".ngmenu." << menuname << " add command -label \"" << text << "\" \\" << endl
		  << "-command {" << endl;
	if(center)
	  {
	    for(int i = centerpoint.Size()-1; i<3; i++) centerpoint.Append(0);
	    tclstring << "set viewoptions.usecentercoords 1" << endl
		      << "set viewoptions.centerx " << centerpoint[0] << endl
		      << "set viewoptions.centery " << centerpoint[1] << endl
		      << "set viewoptions.centerz " << centerpoint[2] << endl
		      << "set dummy $selectvisual" << endl
		      << "set selectvisual \"mesh\"" << endl
		      << "Ng_SetVisParameters; Ng_Center" << endl
		      << "if {$dummy != \"mesh\"} {set selectvisual $dummy; Ng_SetVisParameters}" << endl;
	  }
	if(clip)
	  {
	    for(int i = centerpoint.Size()-1; i<3; i++) clipvec.Append(0);
	    tclstring << "set viewoptions.clipping.enable 1" << endl
		      << "set viewoptions.clipping.nx " << clipvec[0] << endl
		      << "set viewoptions.clipping.ny " << clipvec[1] << endl
		      << "set viewoptions.clipping.nz " << clipvec[2] << endl
		      << "set viewoptions.clipping.dist 0" << endl
		      << "clippingdialog"<<endl;
	    ng_setvispar = true;
	  }
	if(rotation)
	  {
	    for(int i = rotation_pars.Size(); i<4; i++) rotation_pars.Append(0);
	    tclstring << "Ng_ArbitraryRotation";
	    for(int i=0; i<rotation_pars.Size(); i++) tclstring << " " << rotation_pars[i];
	    tclstring << ";" << endl;
	  }

	if(fieldname != "" && !deformationon)
	  {
	    if(evaluate != "") 
	      tclstring << "set visoptions.evaluate " << evaluate<< endl;
	    tclstring << "set visoptions.scalfunction " << fieldname<< "." << component << endl;
	    if(clip)
	      tclstring << "set visoptions.clipsolution scal" << endl;
	    ng_vissetpar = true;
	  }

	if(deformationoff)
	  {
	    tclstring << "set visoptions.deformation 0" << endl;
	    ng_vissetpar = true;
	  }
	if(deformationon && fieldname != "")
	  {
	    tclstring << "set visoptions.deformation 1" << endl
		      << "set visoptions.scaledeform1 " << deformationscale << endl
		      << "set visoptions.scaledeform2 1" << endl
		      << "set visoptions.vecfunction " << fieldname << endl;
	    ng_vissetpar = true;
	  }
	if(light >= 0)
	  {
	    tclstring << "set viewoptions.light.amb " << light << endl;
	    ng_setvispar = true;
	  }

	if(autoscale)
	  {
	    tclstring << "set visoptions.autoscale 1" << endl;
	    ng_vissetpar = true;
	  }

	if(noautoscale)
	  {
	    tclstring << "set visoptions.autoscale 0" << endl
		      << "set visoptions.mminval " << minval << endl
		      << "set visoptions.mmaxval " << maxval << endl;
	    ng_vissetpar = true;
	  }




	if(ng_setvispar)
	  tclstring << "Ng_SetVisParameters" << endl;
	if(ng_vissetpar)
	  tclstring << "Ng_Vis_Set parameters" << endl;


	tclstring << "redraw" << endl
		  << "}" << endl;
	//cout << tclstring.str() << endl;
	
      }


    // the following is a lengthy workaround for a visual-c++ problem
    char *dummy; dummy = new char[tclstring.str().size()+1];
    strcpy(dummy,tclstring.str().c_str());
    
    Tcl_Eval(pde.tcl_interpreter,dummy);
    
    delete [] dummy;
  }

  NumProcTclMenu :: ~NumProcTclMenu(){;}

  void NumProcTclMenu :: PrintDoc (ostream & ost)
  {
    ost <<
      "\n\nNumproc TclMenu:\n" \
      "--------------------------\n" \
      "Adds menu entries for specific visualization\n"\
      "\nFlags:\n"\
      " -menuname=<name>\n"\
      "      menu identifier\n"\
      "\nOptional Flags:\n"\
      " -newmenu\n"\
      "      adds a new menu\n"\
      " -text=<text>\n"\
      "      name of the menu or the menu-entry\n"\
      " -centerpoint=[x,y,z]\n"\
      "      sets the center point\n"\
      " -clipvec=[x,y,z]\n"\
      "      turns on a clipping plane with the given normal-vector and opens the clipping dialog-box\n"\
      " -rotation=[alpha1,vx1,vy1,vz1,alpha2,vx2,vy2,vz2,...]\n"\
      "      rotates the standard view by the angle alpha1 round the vector (vx1,vy1,vz1), then by alpha2 round (vx2,vy2,vz2) etc.\n"\
      " -fieldname=<name>\n"\
      "      name of the field to be displayed\n"\
      " -comp=<component>\n"\
      "      component\n"\
      " -evaluate = < abs | abstens | mises | main >\n"
      " -deformationscale = <value>\n"\
      "      display deformation scaled by given value\n"\
      " -light = <value>\n"\
      "      sets the ambient light to a value between 0 and 1\n"\
      " -autoscale\n"\
      "      turns on autoscale\n"\
      " -minval=<value> -maxval=<value>\n"\
      "      turns off autoscale to use the given minimal and maximal values\n";
  }

  void NumProcTclMenu :: Do(LocalHeap & lh)
  {
    
    
  }


  void NumProcTclMenu :: PrintReport (ostream & ost)
  {
    ost << "NumProcTclMenu:" << endl;
  }



  //////////////////////////////////////////////////////////////////////////////













  NumProcs::NumProcInfo::
  NumProcInfo (const string & aname,
	       NumProc* (*acreator)(PDE & pde, const Flags & flags),
	       void (*aprintdoc) (ostream & ost) )
    : name(aname), creator(acreator), printdoc(aprintdoc)
  {
    ;
  }
  
  NumProcs :: NumProcs ()
  {
    ;
  }

  NumProcs :: ~NumProcs()
  {
    for (int i = 0; i < npa.Size(); i++)
      delete npa[i];
  }
  
  void NumProcs :: 
  AddNumProc (const string & aname,
	      NumProc* (*acreator)(PDE & pde, const Flags & flags),
	      void (*printdoc) (ostream & ost) )
  {
    npa.Append (new NumProcInfo(aname, acreator, printdoc));
  }

  const NumProcs::NumProcInfo * 
  NumProcs::GetNumProc(const string & name)
  {
    for (int i = 0; i < npa.Size(); i++)
      {
	if (name == npa[i]->name)
	  return npa[i];
      }
    return 0;
  }

  void NumProcs :: Print (ostream & ost) const
  {
    ost << endl << "NumProcs:" << endl;
    ost <<         "---------" << endl;
    ost << setw(20) << "Name" << endl;
    for (int i = 0; i < npa.Size(); i++)
      ost << setw(20) << npa[i]->name << endl;
  }



 
  NumProcs & GetNumProcs ()
  {
    static NumProcs nps;
    return nps;
  }


  // standard numprocs:




  namespace {


    class Init
    { 
    public: 
      Init ();
    };
    
    Init::Init()
    {
      GetNumProcs().AddNumProc ("calcflux", NumProcCalcFlux::Create, NumProcCalcFlux::PrintDoc);
      GetNumProcs().AddNumProc ("drawflux", NumProcDrawFlux::Create, NumProcDrawFlux::PrintDoc);
      GetNumProcs().AddNumProc ("evaluate", NumProcEvaluate::Create, NumProcEvaluate::PrintDoc);
      GetNumProcs().AddNumProc ("analyze", NumProcAnalyze::Create, NumProcAnalyze::PrintDoc);
      GetNumProcs().AddNumProc ("warn", NumProcWarn::Create, NumProcWarn::PrintDoc);
      GetNumProcs().AddNumProc ("tclmenu", NumProcTclMenu::Create, NumProcTclMenu::PrintDoc);

      //      GetNumProcs().AddNumProc ("setvisual", NumProcSetVisual::Create);
    }

    
    Init init;
  }
}