/*
 Solves the wave propagation in a rect hollow waveguide problem
 Governing differential equation is the well known vector wave equation:
 curl( (1/mur) curl(E)) - k_0^2 * epsilon_r * E = -j*k_0*Z_0*J,

 with boundary conditions
 n X E = n X F,                    on all the boundaries
 or n X curl(E) = n X curl(F),  on all the boundaries

the exact solution is given by the FEM book (Chapter8, E^inc)
here we want to use the mesh generated from Gmsh to evaluate the electric field
Current problem is when we use the mesh generated using GridGenerator::hyper_rectangle (triangulation, p1, p2)
it works well, when we change the mesh to be that from Gmsh, it does not work

 Author: Jianan Zhang
 Date:   15/12/2017
*/

#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/convergence_table.h>
#include <deal.II/base/numbers.h>

#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/constraint_matrix.h>

#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_boundary_lib.h>

#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_renumbering.h>

#include <deal.II/fe/fe_nedelec.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_system.h>

#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>

#include <fstream>
#include <iostream>
#include <sstream>

using namespace dealii;

// Exact solution CLASS.
class ExactSolution : public Function<3>
{
public:
    ExactSolution();// (Vector<double> k_in, Vector<double> p_in);

    virtual double value (const Point<3> &p,
                          const unsigned int component) const;
    virtual void vector_value (const Point<3> &p,
                               Vector<double> &result) const;
    virtual void value_list (const std::vector<Point<3> > &points,
                             std::vector<double> &values,
                             const unsigned int component) const;
    virtual void vector_value_list (const std::vector<Point<3> > &points,
                                    std::vector<Vector<double> >   &values) const;
    double curl_value(const Point<3> &p,
                      const unsigned int component);
    void curl_value_list(const std::vector<Point<3> > &points,
                         std::vector<Vector<double> > &value_list);
private:
  double dotprod(const Vector<double> &A, const Vector<double> &B) const;
  double dotprod(const Vector<double> &A, const Point<3> &B) const;

  // Material Parameters:
  double param_omega = 2*numbers::PI*10e9; // angular frequency, rads/sec, @10GHz
  double param_mu0 = 4.0*numbers::PI*1e-7;    //permeability of free space, unit henry/meter
  double param_epsilon0 = 8.854*1e-12;    //permittivity of free space, unit farad/meter

  // Waveguide geometrical parameters
  double param_a = 19.05e-3;          //a, unit meter
  double param_b = 9.524e-3;         //b, unit meter
  double param_L = 5.08e-3;       //ridge length in the waveguide, unit meter
  double param_W = 1.016e-3;       //ridge width in the waveguide, unit meter

  // Define other necessary parameters
  double k0_sq;
  double k_z10_re;
  double k_z10_im;
  double gamma_re;
  double gamma_im;
};

ExactSolution::ExactSolution()
: Function<3> (3+3)  //dim+dim is the number of components, by default is 1, i.e. a scalar function
{
    // Evaluate necessary parameters
    k0_sq = param_omega * param_omega * param_epsilon0 * param_mu0;
    if(k0_sq >= (numbers::PI/param_a)*(numbers::PI/param_a)){
     	k_z10_re = sqrt(k0_sq - (numbers::PI/param_a)*(numbers::PI/param_a));
       	k_z10_im = 0.0;
    }
    else{
      	k_z10_re = 0.0;
      	k_z10_im = -sqrt((numbers::PI/param_a)*(numbers::PI/param_a) - k0_sq);
    }

    gamma_re = -k_z10_im;
    gamma_im = k_z10_re;
}

// Exact solution members:
double ExactSolution::dotprod(const Vector<double> &A, const Vector<double> &B) const
{
    double return_val = 0;
    for(unsigned int k = 0; k < 3; k++) {
        return_val += A(k)*B(k);
    }
    return return_val;
}

double ExactSolution::dotprod(const Vector<double> &A, const Point<3> &B) const
{
    double return_val = 0;
    for(unsigned int k = 0; k < 3; k++) {
        return_val += A(k)*B(k);
    }

    return return_val;
}

double ExactSolution::value(const Point<3> &p,
                            unsigned int component) const
{
    AssertIndexRange(component, 3+3);

    double val;
    /* My exact solution */
    if (component == 0 || component == 2 || component == 3 || component == 5)
      val = 0.0;
    else if (component == 1)
      val = sin(numbers::PI*p[0]/param_a)*cos(k_z10_re*p[2])*exp(k_z10_im*p[2]);
    else
      val = -sin(numbers::PI*p[0]/param_a)*sin(k_z10_re*p[2])*exp(k_z10_im*p[2]);
    /* Ross' exact solution */
//    if (component == 0 || component == 2 || component == 3 || component == 5)
//      val = 0.0;
//    else if (component == 1)
//      val = cos(0.1*p[2]);
//    else
//      val = sin(0.1*p[2]);

    return val;
}

void ExactSolution::vector_value (const Point<3> &p,
                                        Vector<double> &values) const
{
    Assert (values.size() == 3, ExcDimensionMismatch (values.size(), 3));   //why not values.size() == 6?
    /* My exact solution */
    // Real:
    values(0) = 0.0;
    values(1) = sin(numbers::PI*p[0]/param_a)*cos(k_z10_re*p[2])*exp(k_z10_im*p[2]);
    values(2) = 0.0;
    // Imaginary:
    values(3) = 0.0;
    values(4) = -sin(numbers::PI*p[0]/param_a)*sin(k_z10_re*p[2])*exp(k_z10_im*p[2]);
    values(5) = 0.0;

    /* Ross' exact solution */
    // Real:
//    values(0) = 0.0;
//    values(1) = cos(0.1*p[2]);
//    values(2) = 0.0;
//    // Imaginary:
//    values(3) = 0.0;
//    values(4) = sin(0.1*p[2]);
//    values(5) = 0.0;
}

void ExactSolution::value_list (const std::vector<Point<3>> &points,
                                      std::vector<double> &values,
                                      const unsigned int component) const
{
    Assert (values.size() == points.size(), ExcDimensionMismatch(values.size(), points.size()));
    AssertIndexRange(component, 3);

    /* My exact solution */
    for (unsigned int i=0; i<points.size(); ++i) {
      if (component == 0 || component == 2 || component == 3 || component == 5)
         values[i] = 0.0;
      else if (component == 1)
         values[i] = sin(numbers::PI*points[i][0]/param_a)*cos(k_z10_re*points[i][2])*exp(k_z10_im*points[i][2]);
      else
         values[i] = -sin(numbers::PI*points[i][0]/param_a)*sin(k_z10_re*points[i][2])*exp(k_z10_im*points[i][2]);
    }
    /* Ross' exact solution */
//    for (unsigned int i=0; i<points.size(); ++i) {
//      if (component == 0 || component == 2 || component == 3 || component == 5)
//         values[i] = 0.0;
//      else if (component == 1)
//         values[i] = cos(0.1*points[i][2]);
//      else
//         values[i] = sin(0.1*points[i][2]);
//    }
}

void ExactSolution::vector_value_list (const std::vector< Point<3> > &points,
                                             std::vector< Vector<double> > &value_list) const
{
    Assert(value_list.size() == points.size(), ExcDimensionMismatch(value_list.size(), points.size()));
    const unsigned int n_points = points.size();
    /* My exact solution */
    for (unsigned int i=0; i<n_points; ++i)
    {
        // Real:
        value_list[i](0) = 0.0;
        value_list[i](1) = sin(numbers::PI*points[i][0]/param_a)*cos(k_z10_re*points[i][2])*exp(k_z10_im*points[i][2]);
        value_list[i](2) = 0.0;
        // Imaginary:
        value_list[i](3) = 0.0;
        value_list[i](4) = -sin(numbers::PI*points[i][0]/param_a)*sin(k_z10_re*points[i][2])*exp(k_z10_im*points[i][2]);
        value_list[i](5) = 0.0;
    }

    /* Ross' exact solution */
//    for (unsigned int i=0; i<n_points; ++i)
//    {
//        // Real:
//        value_list[i](0) = 0.0;
//        value_list[i](1) = cos(0.1*points[i][2]);
//        value_list[i](2) = 0.0;
//        // Imaginary:
//        value_list[i](3) = 0.0;
//        value_list[i](4) = sin(0.1*points[i][2]);
//        value_list[i](5) = 0.0;
//    }
}
// Additional functions to create Neumann conditions.
// evaluate curl(E), in 2D, curl(E) is a complex-valued scalar which has 2 components
double ExactSolution::curl_value(const Point<3> &p,
                                 unsigned int component)
{
    AssertIndexRange(component, 3);

    Vector<double> values(3 + 3);

    /* My exact solution */
    // Real:
    values(0) = -sin(numbers::PI*p[0]/param_a) * exp(k_z10_im*p[2]) * (k_z10_im*cos(k_z10_re*p[2]) - k_z10_re*sin(k_z10_re*p[2]) );
    values(1) = 0.0;
    values(2) = cos(k_z10_re*p[2])*exp(k_z10_im*p[2])*cos(numbers::PI*p[0]/param_a)*numbers::PI/param_a;

    // Imaginary:
    values(3) = sin(numbers::PI*p[0]/param_a) * (sin(k_z10_re*p[2])*exp(k_z10_im*p[2])*k_z10_im + exp(k_z10_im*p[2])*cos(k_z10_re*p[2])*k_z10_re);
    values(4) = 0.0;
    values(5) = -sin(k_z10_re*p[2])*exp(k_z10_im*p[2])*cos(numbers::PI*p[0]/param_a)*numbers::PI/param_a;

    /* Ross' exact solution */
    // Real:
//    values(0) = 0.1*sin(0.1*p[2]);
//    values(1) = 0.0;
//    values(2) = 0.0;
//
//    // Imaginary:
//    values(3) = -0.1*cos(0.1*p[2]);
//    values(4) = 0.0;
//    values(5) = 0.0;

    double val = values(component);

    return val;
}

void ExactSolution::curl_value_list(const std::vector<Point<3>> &points,
                                          std::vector<Vector<double> > &value_list)
{
  Assert(value_list.size() == points.size(), ExcDimensionMismatch(value_list.size(), points.size()));
  const unsigned int n_points = points.size();

  /* My exact solution */
  for (unsigned int i=0; i<n_points; ++i)
  {
      // Real:
      value_list[i](0) = -sin(numbers::PI*points[i][0]/param_a) * exp(k_z10_im*points[i][2]) * (k_z10_im*cos(k_z10_re*points[i][2]) - k_z10_re*sin(k_z10_re*points[i][2]) );
      value_list[i](1) = 0.0;
      value_list[i](2) = cos(k_z10_re*points[i][2])*exp(k_z10_im*points[i][2])*cos(numbers::PI*points[i][0]/param_a)*numbers::PI/param_a;

      // Imaginary:
      value_list[i](3) = sin(numbers::PI*points[i][0]/param_a) * (sin(k_z10_re*points[i][2])*exp(k_z10_im*points[i][2])*k_z10_im + exp(k_z10_im*points[i][2])*cos(k_z10_re*points[i][2])*k_z10_re);
      value_list[i](4) = 0.0;
      value_list[i](5) = -sin(k_z10_re*points[i][2])*exp(k_z10_im*points[i][2])*cos(numbers::PI*points[i][0]/param_a)*numbers::PI/param_a;
  }

  /* Ross' exact solution */
//  for (unsigned int i=0; i<n_points; ++i)
//  {
//      const Point<3> &p = points[i];
//
//      // Real:
//      value_list[i](0) = 0.1*sin(0.1*points[i][2]);
//      value_list[i](1) = 0.0;
//      value_list[i](2) = 0.0;
//      // Imaginary:
//      value_list[i](3) = -0.1*cos(0.1*points[i][2]);
//      value_list[i](4) = 0.0;
//      value_list[i](5) = 0.0;
//  }
}
// END EXACT SOLUTION

// COMPUTE E_re CLASS
class ComputeEr : public DataPostprocessorVector<3>
{
  public:
    ComputeEr ();
    virtual void compute_derived_quantities_vector (
        const std::vector< Vector< double > > &uh,
        const std::vector< std::vector< Tensor< 1, 3 > > > &duh,
        const std::vector< std::vector< Tensor< 2, 3 > > > &dduh,
        const std::vector< Point< 3 > > &normals,
        const std::vector<Point<3> > &evaluation_points,
        std::vector< Vector< double > > &computed_quantities) const;
};

ComputeEr::ComputeEr ()
  : DataPostprocessorVector<3> ("Er", update_values)
{}

void ComputeEr::compute_derived_quantities_vector (
    const std::vector< Vector< double > >                  &uh,
    const std::vector< std::vector< Tensor< 1, 3 > > >  & /*duh*/,
    const std::vector< std::vector< Tensor< 2, 3 > > >  & /*dduh*/,
    const std::vector< Point< 3 > >                     & /*normals*/,
    const std::vector<Point<3> >                        & /*evaluation_points*/,
    std::vector< Vector< double > >                        &computed_quantities
    ) const
{
  Assert(computed_quantities.size() == uh.size(),
      ExcDimensionMismatch (computed_quantities.size(), uh.size()));

  for (unsigned int i = 0; i < computed_quantities.size(); i++) {
    Assert(computed_quantities[i].size() == 3,
        ExcDimensionMismatch (computed_quantities[i].size(), 3));
    Assert(uh[i].size() == 3 * 2, ExcDimensionMismatch (uh[i].size(), 3 * 2));
    for (unsigned int j = 0; j < 3; ++j) {
      computed_quantities[i][j] = uh[i][j];
    }
  }
}

// COMPUTE E_im CLASS
class ComputeEi : public DataPostprocessorVector<3>
{
  public:
    ComputeEi ();
    virtual void compute_derived_quantities_vector (
        const std::vector< Vector< double > > &uh,
        const std::vector< std::vector< Tensor< 1, 3 > > > &duh,
        const std::vector< std::vector< Tensor< 2, 3 > > > &dduh,
        const std::vector< Point< 3 > > &normals,
        const std::vector<Point<3> > &evaluation_points,
        std::vector< Vector< double > > &computed_quantities) const;
};

ComputeEi::ComputeEi ()
  : DataPostprocessorVector<3> ("Ei", update_values)
{}

void ComputeEi::compute_derived_quantities_vector (
    const std::vector< Vector< double > >                  &uh,
    const std::vector< std::vector< Tensor< 1, 3 > > >  & /*duh*/,
    const std::vector< std::vector< Tensor< 2, 3 > > >  & /*dduh*/,
    const std::vector< Point< 3 > >                     & /*normals*/,
    const std::vector<Point<3> >                        & /*evaluation_points*/,
    std::vector< Vector< double > >                        &computed_quantities
    ) const
{
  Assert(computed_quantities.size() == uh.size(),
      ExcDimensionMismatch (computed_quantities.size(), uh.size()));

  for (unsigned int i = 0; i < computed_quantities.size(); i++) {
    Assert(computed_quantities[i].size() == 3,
        ExcDimensionMismatch (computed_quantities[i].size(), 3));
    Assert(uh[i].size() == 3 * 2, ExcDimensionMismatch (uh[i].size(), 3 * 2));
    for (unsigned int j = 0; j < 3; ++j) {
      computed_quantities[i][j] = uh[i][j + 3];
    }
  }
}

// COMPUTE E_mag CLASS
class ComputeEmag : public DataPostprocessorVector<3>
{
  public:
    ComputeEmag ();
    virtual void compute_derived_quantities_vector (
        const std::vector< Vector< double > > &uh,
        const std::vector< std::vector< Tensor< 1, 3 > > > &duh,
        const std::vector< std::vector< Tensor< 2, 3 > > > &dduh,
        const std::vector< Point< 3 > > &normals,
        const std::vector<Point<3> > &evaluation_points,
        std::vector< Vector< double > > &computed_quantities) const;
};

ComputeEmag::ComputeEmag ()
  : DataPostprocessorVector<3> ("Emag", update_values)
{}

void ComputeEmag::compute_derived_quantities_vector (
    const std::vector< Vector< double > >                  &uh,
    const std::vector< std::vector< Tensor< 1, 3 > > >  & /*duh*/,
    const std::vector< std::vector< Tensor< 2, 3 > > >  & /*dduh*/,
    const std::vector< Point< 3 > >                     & /*normals*/,
    const std::vector<Point<3> >                        & /*evaluation_points*/,
    std::vector< Vector< double > >                        &computed_quantities
    ) const
{
  Assert(computed_quantities.size() == uh.size(),
      ExcDimensionMismatch (computed_quantities.size(), uh.size()));

  for (unsigned int i = 0; i < computed_quantities.size(); i++) {
    Assert(computed_quantities[i].size() == 3,
        ExcDimensionMismatch (computed_quantities[i].size(), 3));
    Assert(uh[i].size() == 3 * 2, ExcDimensionMismatch (uh[i].size(), 3 * 2));
    for (unsigned int j = 0; j < 3; ++j) {
      computed_quantities[i][j] = sqrt(uh[i][j]*uh[i][j] + uh[i][j+3]*uh[i][j+3]);
    }
  }
}

// MAIN Waveguide CLASS
class WaveguideProblem
{
public:
    WaveguideProblem (const unsigned int order);
    ~WaveguideProblem ();
    void run ();

private:
    double dotprod(const Tensor<1,3> &A, const Tensor<1,3> &B) const;
    double dotprod(const Tensor<1,3> &A, const Vector<double> &B) const;
    //void read_in_mesh (std::string mesh_file);
    void setup_system ();
    void assemble_system ();
    void solve ();
    void process_solution(const unsigned int cycle);
    void output_results_vtk (const unsigned int cycle) const;
    double calcErrorHcurlNorm();

    Triangulation<3>   triangulation;
    DoFHandler<3>      dof_handler;
    FESystem<3>        fe;
    ConstraintMatrix     constraints;
    SparsityPattern      sparsity_pattern;
    SparseMatrix<double> system_matrix;
    Vector<double>       solution;
    Vector<double>       system_rhs;
    
    ExactSolution        exact_solution;

    ConvergenceTable     convergence_table;
    
    // Input parameters
    unsigned int p_order;
    unsigned int quad_order;
    
    // Material Parameters:
    double param_omega = 2*numbers::PI*10e9; // angular frequency, rads/sec, @10GHz
    double param_mu0 = 4.0*numbers::PI*1e-7;    //permeability of free space, unit henry/meter
    double param_epsilon0 = 8.854*1e-12;    //permittivity of free space, unit farad/meter

    // Waveguide geometrical parameters
    double param_a = 19.05e-3;          //a, unit meter
    double param_b = 9.524e-3;         //b, unit meter
    double param_L = 5.08e-3;       //ridge length in the waveguide, unit meter
    double param_W = 1.016e-3;       //ridge width in the waveguide, unit meter

    // Define other necessary parameters
    double k0_sq;
    double k_z10_re;
    double k_z10_im;
    double gamma_re;
    double gamma_im;
};

WaveguideProblem::WaveguideProblem (const unsigned int order)
:
dof_handler (triangulation),
// Defined as FESystem, and we need 2 FE_Nedelec - first (0) is real part, second (1) is imaginary part.
fe (FE_Nedelec<3>(order), 1,
	FE_Nedelec<3>(order), 1)
{
    p_order = order;
    quad_order = p_order+2;

    // Evaluate necessary parameters
    k0_sq = param_omega * param_omega * param_epsilon0 * param_mu0;

    if(k0_sq >= (numbers::PI/param_a)*(numbers::PI/param_a)){
     	k_z10_re = sqrt(k0_sq - (numbers::PI/param_a)*(numbers::PI/param_a));
     	std::cout << "k_z10_re is: " << k_z10_re << std::endl;
       	k_z10_im = 0.0;
    }
    else{
      	k_z10_re = 0.0;
      	k_z10_im = -sqrt((numbers::PI/param_a)*(numbers::PI/param_a) - k0_sq);
      	std::cout << "k_z10_im is: " << k_z10_im << std::endl;
    }

    gamma_re = -k_z10_im;
    gamma_im = k_z10_re;
}

WaveguideProblem::~WaveguideProblem ()
{
    dof_handler.clear ();
}

double WaveguideProblem::dotprod(const Tensor<1,3> &A, const Tensor<1,3> &B) const
{
    double return_val = 0;
    for(unsigned int k = 0; k < 3; k++) {
        return_val += A[k]*B[k];
    }
    return return_val;
}

double WaveguideProblem::dotprod(const Tensor<1,3> &A, const Vector<double> &B) const
{
    double return_val = 0;
    for(unsigned int k = 0; k < 3; k++) {
        return_val += A[k]*B(k);
    }
    return return_val;
}

double WaveguideProblem::calcErrorHcurlNorm()
{
     QGauss<3>  quadrature_formula(quad_order);
     const unsigned int n_q_points = quadrature_formula.size();

     FEValues<3> fe_values (fe, quadrature_formula,
                              update_values    |  update_gradients |
                              update_quadrature_points  |  update_JxW_values);

     // Extractors to real and imaginary parts
     const FEValuesExtractors::Vector E_re(0);
     const FEValuesExtractors::Vector E_im(3);

     const unsigned int dofs_per_cell = fe.dofs_per_cell;

     std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);

     // storage for exact sol:
     std::vector<Vector<double> > exactsol(n_q_points, Vector<double>(fe.n_components()));
     std::vector<Vector<double> > exactcurlsol(n_q_points, Vector<double>(fe.n_components()));
     Tensor<1,3>  exactsol_re;
     Tensor<1,3>  exactsol_im;
     Tensor<1,3>  exactcurlsol_re;
     Tensor<1,3>  exactcurlsol_im;

     // storage for computed sol:
     std::vector<Tensor<1,3>> sol_re(n_q_points);
     std::vector<Tensor<1,3>> sol_im(n_q_points);
     Tensor<1,3> curlsol_re;
     Tensor<1,3> curlsol_im;

     double h_curl_norm=0.0;

     unsigned int block_index_i;

     typename DoFHandler<3>::active_cell_iterator

     cell = dof_handler.begin_active(),
     endc = dof_handler.end();

     for (; cell!=endc; ++cell)
     {
         fe_values.reinit (cell);

         // Store exact values of E and curlE:
         exact_solution.vector_value_list(fe_values.get_quadrature_points(), exactsol);
         exact_solution.curl_value_list(fe_values.get_quadrature_points(), exactcurlsol);

         // Store computed values at quad points:
         fe_values[E_re].get_function_values(solution, sol_re);
         fe_values[E_im].get_function_values(solution, sol_im);

         // Calc values of curlE from fe solution:
         cell->get_dof_indices (local_dof_indices);
         // Loop over quad points to calculate solution:
         for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
         {
             // Split exact solution into real/imaginary parts:
             for (unsigned int component=0;component<3;component++)
             {
                 exactsol_re[component] = exactsol[q_point][component];
                 exactsol_im[component] = exactsol[q_point][component+3];
                 exactcurlsol_re[component] = exactcurlsol[q_point][component];
                 exactcurlsol_im[component] = exactcurlsol[q_point][component+3];
             }
             // Loop over DoFs to calculate curl of solution @ quad point
             curlsol_re=0.0;
             curlsol_im=0.0;
             for (unsigned int i=0; i<dofs_per_cell; ++i)
             {
                 block_index_i = fe.system_to_block_index(i).first;
                 // Construct local curl value @ quad point
                 if (block_index_i==0)
                 {
                     curlsol_re += solution(local_dof_indices[i])*fe_values[E_re].curl(i,q_point);
                 }
                 else if (block_index_i==1)
                 {
                     curlsol_im += solution(local_dof_indices[i])*fe_values[E_im].curl(i,q_point);
                 }
             }
             // Integrate difference at each point:
             h_curl_norm += ( (exactsol_re-sol_re[q_point])*(exactsol_re-sol_re[q_point])
                             + (exactsol_im-sol_im[q_point])*(exactsol_im-sol_im[q_point])
                             + (exactcurlsol_re-curlsol_re)*(exactcurlsol_re-curlsol_re)
                             + (exactcurlsol_im-curlsol_im)*(exactcurlsol_im-curlsol_im)
                             )*fe_values.JxW(q_point);
         }
     }

    return sqrt(h_curl_norm);
}

void WaveguideProblem::setup_system ()
{
    dof_handler.distribute_dofs (fe);
    std::cout << "Number of degrees of freedom: "
                << dof_handler.n_dofs()
                << std::endl;
    DoFRenumbering::block_wise (dof_handler);   //Sort the degrees of freedom by vector block
    solution.reinit (dof_handler.n_dofs());
    system_rhs.reinit (dof_handler.n_dofs());
    
    std::fstream sys_rhs_output_file1;
    sys_rhs_output_file1.open("sys_rhs_before_project.txt", std::fstream::out);
    system_rhs.print(sys_rhs_output_file1);
    sys_rhs_output_file1.close();

    constraints.clear ();   
    DoFTools::make_hanging_node_constraints (dof_handler, constraints);

    // Implement Dirichlet boundary condition on the waveguide walls
    // Real part (begins at 0)
    VectorTools::project_boundary_values_curl_conforming_l2(dof_handler,
							 	 	 	 	 	 	 	    0,
							 	 	 	 	 	 	 	    exact_solution,
														    0,
														    constraints);
    // Imaginary part (begins at 3)
    VectorTools::project_boundary_values_curl_conforming_l2(dof_handler,
														    3,
														    exact_solution,
														    0,
														    constraints);
    constraints.close ();
    
    std::fstream sys_rhs_output_file2;
    sys_rhs_output_file2.open("sys_rhs_after_project.txt", std::fstream::out);
    system_rhs.print(sys_rhs_output_file2);
    sys_rhs_output_file2.close();

    DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
    DoFTools::make_sparsity_pattern(dof_handler,
                                    dsp,
                                    constraints,
                                    false);
    
    sparsity_pattern.copy_from(dsp);
    system_matrix.reinit (sparsity_pattern);
}

void WaveguideProblem::assemble_system ()
{
    QGauss<3>  quadrature_formula(quad_order);   //quad_order = p_order + 2, = 2
    QGauss<2>  face_quadrature_formula(quad_order);
    
    FEValues<3> fe_values (fe, quadrature_formula,
                             update_values | update_gradients |
                             update_quadrature_points | update_JxW_values);
    FEFaceValues<3> fe_face_values(fe, face_quadrature_formula,
                                     update_values | update_quadrature_points |
                                     update_normal_vectors | update_JxW_values);
    
    const unsigned int dofs_per_cell = fe.dofs_per_cell;
    const unsigned int n_q_points    = quadrature_formula.size();
    const unsigned int n_face_q_points = face_quadrature_formula.size();

    FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
    Vector<double> cell_rhs (dofs_per_cell);
    
    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
   
    // Extractors to real and imaginary parts
    const FEValuesExtractors::Vector E_re(0);
    const FEValuesExtractors::Vector E_im(3);

    // Neumann storage
    std::vector<Vector<double>> neumann_value_list(n_face_q_points, Vector<double>(fe.n_components()));
    Tensor<1,3> neumann_value_list_re;
    Tensor<1,3> neumann_value_list_im;
    Tensor<1,3> neumann_value_re;
    Tensor<1,3> neumann_value_im;
    Tensor<1,3> normal_vector;

    unsigned int block_index_i;
    unsigned int block_index_j;
    
    typename DoFHandler<3>::active_cell_iterator
    cell = dof_handler.begin_active(),
    endc = dof_handler.end();

    for (; cell!=endc; ++cell)
    {
    	fe_values.reinit (cell);
        cell_matrix = 0.;
        cell_rhs = 0.;

        // Loop over quad points:
        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
        {
            for (unsigned int i=0; i<dofs_per_cell; ++i)
            {
                block_index_i = fe.system_to_block_index(i).first;
                // Construct local matrix:
                for (unsigned int j=0; j<dofs_per_cell; ++j)
                {
                    block_index_j = fe.system_to_block_index(j).first;
                    // block 0 = real, block 1 = imag
                    // free space, mur = 1, epsilonr = 1
                    if (block_index_i == block_index_j)
                    {
                    	/* Ross assembly*/
//                        cell_matrix(i,j) += ( (fe_values[E_re].curl(i,q_point)*fe_values[E_re].curl(j,q_point)
//                                                                 + fe_values[E_im].curl(i,q_point)*fe_values[E_im].curl(j,q_point))
//                                             - 0.01*( fe_values[E_re].value(i,q_point)*fe_values[E_re].value(j,q_point)
//                                                                 + fe_values[E_im].value(i,q_point)*fe_values[E_im].value(j,q_point))
//                                             )*fe_values.JxW(q_point);
                        /* My assembly */
                        cell_matrix(i,j) += ( (fe_values[E_re].curl(i,q_point)*fe_values[E_re].curl(j,q_point)
                                                                 + fe_values[E_im].curl(i,q_point)*fe_values[E_im].curl(j,q_point))
                                             - k0_sq*( fe_values[E_re].value(i,q_point)*fe_values[E_re].value(j,q_point)
                                                                 + fe_values[E_im].value(i,q_point)*fe_values[E_im].value(j,q_point))
                                             )*fe_values.JxW(q_point);
                    }
                    else
                    {
                        if (block_index_i == 0)
                        {
                        	// previously, Kappa = Kappa_re + j*Kappa_im = -omega^2*epsilon + j*omega*sigma,
                            cell_matrix(i,j) += -0.0*(  fe_values[E_re].value(i,q_point)*fe_values[E_im].value(j,q_point)
                                                                   + fe_values[E_im].value(i,q_point)*fe_values[E_re].value(j,q_point)
                                                                   )*fe_values.JxW(q_point);
                        }
                        else if (block_index_i == 1)
                        {
                            cell_matrix(i,j) += 0.0*(  fe_values[E_re].value(i,q_point)*fe_values[E_im].value(j,q_point)
                                                                  + fe_values[E_im].value(i,q_point)*fe_values[E_re].value(j,q_point)
                                                                  )*fe_values.JxW(q_point);
                        }
                    }
                }
            }
        }

        // Loop over faces for neumann condition:
        for (unsigned int face_number = 0; face_number < GeometryInfo<3>::faces_per_cell; ++face_number)
        {
            fe_face_values.reinit (cell, face_number);
            if (cell->face(face_number)->at_boundary() && cell->face(face_number)->boundary_id() == 1 )
            {
                exact_solution.curl_value_list(fe_face_values.get_quadrature_points(), neumann_value_list);
                for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                {
                    for (unsigned int component=0; component<3; component++)
                    {
                        neumann_value_list_re[component] = 1.0*neumann_value_list[q_point](component);
                        neumann_value_list_im[component] = 1.0*neumann_value_list[q_point](component+3);
                        normal_vector[component] = fe_face_values.normal_vector(q_point)[component];
                    }
                    neumann_value_re = cross_product_3d(normal_vector, neumann_value_list_re);
                    neumann_value_im = cross_product_3d(normal_vector, neumann_value_list_im);

                    for (unsigned int i=0; i<dofs_per_cell; ++i)
                    {
                        cell_rhs(i) += -(neumann_value_re*fe_face_values[E_re].value(i,q_point)*fe_face_values.JxW(q_point));
                        cell_rhs(i) += -(neumann_value_im*fe_face_values[E_im].value(i,q_point)*fe_face_values.JxW(q_point));
                    }
                }
            }
        }

        cell->get_dof_indices (local_dof_indices);
        constraints.distribute_local_to_global(cell_matrix,
                                               cell_rhs,
                                               local_dof_indices,
                                               system_matrix, system_rhs);

    }
}

void WaveguideProblem::solve ()
{
    /* Direct */
    SparseDirectUMFPACK A_direct;
    A_direct.initialize(system_matrix);

    A_direct.vmult (solution, system_rhs);
    constraints.distribute (solution);

	/* CG */
//	SolverControl           solver_control (20000, 1e-6);
//	SolverCG<>              solver (solver_control);
//	PreconditionSSOR<> preconditioner;
//	preconditioner.initialize(system_matrix, 1.2);
////	std::cout << "system_rhs:" <<std::endl;
////	system_rhs.print(std::cout);
//	solver.solve (system_matrix,
//			      solution,
//			      system_rhs,
//	              preconditioner);
//	std::cout << "   " << solver_control.last_step()
//	          << " CG iterations needed to obtain convergence."
//	          << std::endl;
//	constraints.distribute (solution);
}

void WaveguideProblem::process_solution(const unsigned int cycle)
{
    // Masks for real & imaginary parts
    const ComponentSelectFunction<3> E_re_mask (std::make_pair(0,3), 3+3);
    const ComponentSelectFunction<3> E_im_mask (std::make_pair(3, 3+3), 3+3);
    
    Vector<double> diff_per_cell_re(triangulation.n_active_cells());
    Vector<double> diff_per_cell_im(triangulation.n_active_cells());
    
    VectorTools::integrate_difference(dof_handler, solution, exact_solution,
                                      diff_per_cell_re, QGauss<3>(quad_order+1),
                                      VectorTools::L2_norm,
                                      &E_re_mask);
    
    
    VectorTools::integrate_difference(dof_handler, solution, exact_solution,
                                      diff_per_cell_im, QGauss<3>(quad_order+1),
                                      VectorTools::L2_norm,
                                      &E_im_mask);
    
    const double L2_error_re = diff_per_cell_re.l2_norm();
    const double L2_error_im = diff_per_cell_im.l2_norm();
    const double L2_error = sqrt(L2_error_re*L2_error_re + L2_error_im*L2_error_im);
    
    double Hcurl_error = calcErrorHcurlNorm();

    convergence_table.add_value("cycle", cycle);
    convergence_table.add_value("cells", triangulation.n_active_cells());
    convergence_table.add_value("dofs", dof_handler.n_dofs());
    convergence_table.add_value("L2 Error", L2_error);
    convergence_table.add_value("H(curl) Error", Hcurl_error);
}


void WaveguideProblem::output_results_vtk (const unsigned int cycle) const
{
	if(cycle == 0)
	{
		ComputeEr Er;
		ComputeEi Ei;
		ComputeEmag Emag;

		DataOut<3, DoFHandler<3>> data_out;

		data_out.attach_dof_handler(dof_handler);
		data_out.add_data_vector(solution, Er);
		data_out.add_data_vector(solution, Ei);
		data_out.add_data_vector(solution, Emag);
		data_out.build_patches();

		const std::string filename = "waveguide_solution_10G.vtk";
		std::ofstream output(filename.c_str());
		data_out.write_vtk(output);
	}
	else
	{ }
}

void WaveguideProblem::run ()
{
    for (unsigned int cycle=0; cycle<1; ++cycle)
    {
        std::cout << "Cycle " << cycle << ':';
        if (cycle == 0)
        {
			/* Read the mesh in*/
//			GridIn<3> grid_in;
//			grid_in.attach_triangulation (triangulation);
//			std::ifstream input_file("hollow_v2.msh");
//			grid_in.read_msh(input_file);
//			std::cout << "number of cells:" <<triangulation.n_active_cells() << std::endl;

        	/*Create a rectangle 3D mesh*/
            const Point<3> p1 (0, 0, 0);
            const Point<3> p2 (param_a, param_b, 3*param_L);
            GridGenerator::hyper_rectangle (triangulation, p1, p2);

			//Set boundaries to neumann (boundary_id = 1 for port 1, boundary_id = 2 for port 2)
//			typename Triangulation<3>::cell_iterator
//			cell = triangulation.begin (),
//			endc = triangulation.end();
//			for (; cell!=endc; ++cell)
//			{
//				for (unsigned int f_number=0; f_number < GeometryInfo<3>::faces_per_cell; ++f_number)
//				{
//					if (cell->face(f_number)->at_boundary())
//					{
//						cell->face(f_number)->set_boundary_id (1);
//					}
//				}
//			}

            triangulation.refine_global (2);
//            std::cout << "number of cells:" <<triangulation.n_active_cells() << std::endl;
//            std::ofstream out2 ("refined_grid.eps");
//            GridOut grid_out2;
//            grid_out2.write_eps (triangulation, out2);
//            std::cout << "Grid written to refined_grid.eps" << std::endl;
        }
        else
        {
           triangulation.refine_global (1);
        }

        setup_system ();
        assemble_system ();
        solve ();
        process_solution (cycle);
        output_results_vtk (cycle);
        std::cout << " Done" << std::endl;
    }

    convergence_table.set_precision("L2 Error",8);
    convergence_table.set_scientific("L2 Error",true);

    convergence_table.set_precision("H(curl) Error",8);
    convergence_table.set_scientific("H(curl) Error",true);

    convergence_table.write_text(std::cout);
}
// END Waveguide CLASS

int main ()
{
    try
    {
    	unsigned int p_order = 1;
    	WaveguideProblem maxwell(p_order);
    	maxwell.run ();
    }
    catch (std::exception &exc)
      {
        std::cerr << std::endl << std::endl
                  << "----------------------------------------------------"
                  << std::endl;
        std::cerr << "Exception on processing: " << std::endl
                  << exc.what() << std::endl
                  << "Aborting!" << std::endl
                  << "----------------------------------------------------"
                  << std::endl;
        return 1;
      }
    catch (...)
      {
        std::cerr << std::endl << std::endl
                  << "----------------------------------------------------"
                  << std::endl;
        std::cerr << "Unknown exception!" << std::endl
                  << "Aborting!" << std::endl
                  << "----------------------------------------------------"
                  << std::endl;
        return 1;
      }

    return 0;
}