Distributed_Moving_Laser_Heating.h

Go to the documentation of this file.

 
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Includefiles"></a> 
 * <h3>Include files</h3>
 * 
 
 * 
 * All the include files have already been discussed in previous tutorials.
 * 
 * @code
 *   #include <deal.II/dofs/dof_handler.h>
 *   #include <deal.II/dofs/dof_renumbering.h>
 *   #include <deal.II/grid/grid_generator.h>
 *   #include <deal.II/grid/tria_accessor.h>
 *   #include <deal.II/grid/tria_iterator.h>
 *   #include <deal.II/dofs/dof_accessor.h>
 *   #include <deal.II/fe/fe_q.h>
 *   #include <deal.II/dofs/dof_tools.h>
 *   #include <deal.II/fe/fe_values.h>
 *   #include <deal.II/base/quadrature_lib.h>
 *   #include <deal.II/base/function.h>
 *   #include <deal.II/numerics/vector_tools.h>
 *   #include <deal.II/numerics/matrix_tools.h>
 *   #include <deal.II/lac/vector.h>
 *   #include <deal.II/lac/full_matrix.h>
 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
 *   
 *   #include <deal.II/grid/grid_in.h>
 *   #include <deal.II/grid/grid_tools.h>
 *   #include <deal.II/lac/trilinos_vector.h>
 *   #include <deal.II/lac/trilinos_sparse_matrix.h>
 *   #include <deal.II/lac/trilinos_solver.h>
 *   #include <deal.II/lac/trilinos_precondition.h>
 *   
 *   #include <deal.II/lac/sparsity_tools.h>
 *   #include <deal.II/lac/generic_linear_algebra.h>
 *   
 *   #include <deal.II/base/conditional_ostream.h>
 *   #include <deal.II/base/utilities.h>
 *   #include <deal.II/base/index_set.h>
 *   #include <deal.II/distributed/tria.h>
 *   
 *   #include <deal.II/numerics/data_out.h>
 *   #include <fstream>
 *   #include <iostream>
 *   
 *   
 *   #include <deal.II/base/logstream.h>
 *   #include <deal.II/lac/affine_constraints.h>
 *   #include <deal.II/base/timer.h>
 *   
 * @endcode
 * 
 * grid refine
 * 
 * @code
 *   #include <deal.II/numerics/error_estimator.h>
 *   #include <deal.II/distributed/grid_refinement.h>
 *   #include <deal.II/distributed/solution_transfer.h>
 *   
 * @endcode
 * 
 * remember to use the namespace dealii before
 * defining a class that is inheritaged from Function<dim>
 * 
 * @code
 *   using namespace dealii;
 *   
 * @endcode
 * 
 * In general, it is more clear to separate boundary and initial condictions,
 * as well as the right hand side function from a file holding all the things.
 * To do so, in this work, the globalPara.h file defines physical constants, 
 * laser parameters, and heat characteristics of materials involved. Boundary
 * and initial conditions are defined in boundaryInit.h. The rightHandSide.h
 * defines the heat source, which in this work is a moving Gaussian beam.
 * 
 
 * 
 * 
 * @code
 *   #ifndef GLOBAL_PARA
 *   #define GLOBAL_PARA
 *   #include "./globalPara.h"
 *   #include "./boundaryInit.h"
 *   #include "./rightHandSide.h"
 *   #endif
 *   
 * @endcode
 * 
 * Now the main class is defined as following
 * 
 * @code
 *   template <int dim>
 *   class LaserHeating
 *   {
 *   public:
 *     LaserHeating ();
 *     ~LaserHeating ();
 *     void run ();
 *   
 *   private:
 *     void make_grid ();
 *     void setup_system();
 *   
 *     void assemble_system_matrix_init (double time_step);
 *     void dynamic_assemble_rhs_T  (double time, double time_step);
 *   
 *     void solve_T ();
 *   
 *     void refine_mesh();
 *     void output_results (int output_num) const;
 *     
 *     MPI_Comm                                    mpi_communicator;
 *   
 *     parallel::distributed::Triangulation<dim>   triangulation;
 *     FE_Q<dim>            fe;
 *     DoFHandler<dim>      dof_handler;
 *   
 *     AffineConstraints<double>     constraints_T;
 *   
 *   
 * @endcode
 * 
 * system_matrix
 * 
 * @code
 *     TrilinosWrappers::SparseMatrix system_matrix_T;
 *   
 * @endcode
 * 
 * for storing left matrix
 * 
 * @code
 *     TrilinosWrappers::SparseMatrix left_system_matrix_T;
 *   
 * @endcode
 * 
 * for storing right matrix
 * 
 * @code
 *     TrilinosWrappers::SparseMatrix right_system_matrix_T;
 *   
 * @endcode
 * 
 * System_rhs, only locally owned cells
 * 
 * @code
 *     TrilinosWrappers::MPI::Vector       system_rhs_T;
 *   
 * @endcode
 * 
 * Solutions
 * Old Solutions with ghost cells, for output
 * 
 * @code
 *     TrilinosWrappers::MPI::Vector       old_solution_T; 
 *   
 * @endcode
 * 
 * Old Solutions only with locally owned cells     
 * 
 * @code
 *     TrilinosWrappers::MPI::Vector       old_solution_T_cal; 
 *   
 * @endcode
 * 
 * New Solutions only with locally owned cells      
 * 
 * @code
 *     TrilinosWrappers::MPI::Vector       new_solution_T; 
 *   
 * @endcode
 * 
 * Dynamic assembling of the righthandside terms
 * 
 * @code
 *     TrilinosWrappers::MPI::Vector       dynamic_rhs_T;
 *   
 *   
 *     IndexSet              locally_owned_dofs;
 *     IndexSet              locally_relevant_dofs;
 *   
 *     ConditionalOStream    pcout;
 *     TimerOutput           computing_timer;
 *   
 *     double theta;
 *   
 *   
 *   
 *   };
 *   
 * @endcode
 * 
 * the constructor
 * 
 * @code
 *   template <int dim>
 *   LaserHeating<dim>::LaserHeating ()
 *     :
 *     mpi_communicator (MPI_COMM_WORLD),
 *     triangulation (mpi_communicator),
 *     fe (1),
 *     dof_handler (triangulation),
 *     pcout (std::cout,Utilities::MPI::this_mpi_process(mpi_communicator) == 0),
 *     computing_timer (mpi_communicator, pcout, TimerOutput::summary, TimerOutput::wall_times),
 *     theta(0.5)
 *   {}
 *   
 * @endcode
 * 
 * the destructor
 * 
 * @code
 *   template <int dim>
 *   LaserHeating<dim>::~LaserHeating ()
 *   {
 *       dof_handler.clear();
 *   }
 *   
 * @endcode
 * 
 * 
 * <a name="LaserHeatingmake_grid"></a> 
 * <h4>LaserHeating::make_grid</h4>
 * make grid by importing msh file, and rescale
 * 
 * @code
 *   template <int dim>
 *   void LaserHeating<dim>::make_grid ()
 *   {
 *     TimerOutput::Scope t(computing_timer,"make_grid()");
 *   
 *     GridIn<dim> grid_in;
 *     grid_in.attach_triangulation (triangulation);
 *     std::ifstream input_file ("geometry.msh");
 *     grid_in.read_msh (input_file);
 *     GridTools::scale (1e-6,triangulation);
 *   
 *     pcout << "   Number of active cells: "
 *               << triangulation.n_global_active_cells()
 *               << std::endl
 *               << "   Total number of cells: "
 *               << triangulation.n_cells()
 *               << std::endl;
 *   }
 *   
 * @endcode
 * 
 * 
 * <a name="LaserHeatingsetup_system"></a> 
 * <h4>LaserHeating::setup_system</h4>
 * initialization
 * 
 * @code
 *   template <int dim>
 *   void LaserHeating<dim>::setup_system ()
 *   {
 *   
 *     TimerOutput::Scope t(computing_timer,"setup_system()");
 *   
 *     dof_handler.distribute_dofs (fe);
 *   
 *     pcout << "   Number of degrees of freedom: "
 *               << dof_handler.n_dofs()
 *               << std::endl;
 *   
 *     locally_owned_dofs = dof_handler.locally_owned_dofs();
 *     DoFTools::extract_locally_relevant_dofs (dof_handler, locally_relevant_dofs);
 *   
 * @endcode
 * 
 * we want to output solution, so here should have ghost cells
 * 
 * @code
 *     old_solution_T.reinit     (locally_owned_dofs,locally_relevant_dofs,mpi_communicator);
 *   
 * @endcode
 * 
 * locally owned cells
 * 
 * @code
 *     old_solution_T_cal.reinit (locally_owned_dofs,mpi_communicator);
 *     new_solution_T.reinit     (locally_owned_dofs,mpi_communicator);
 *     dynamic_rhs_T.reinit      (locally_owned_dofs,mpi_communicator);
 *     system_rhs_T.reinit       (locally_owned_dofs,mpi_communicator);
 *   
 *     constraints_T.clear();
 *     constraints_T.reinit (locally_relevant_dofs);
 *     DoFTools::make_hanging_node_constraints (dof_handler, constraints_T);
 *     constraints_T.close();
 *   
 *     const unsigned int myid = Utilities::MPI::this_mpi_process (mpi_communicator);
 *     DynamicSparsityPattern dsp_T(locally_relevant_dofs);
 *     DoFTools::make_sparsity_pattern (dof_handler, dsp_T, constraints_T, false, myid);
 *   
 *     SparsityTools::distribute_sparsity_pattern (dsp_T,
 *                                                 locally_owned_dofs,
 *                                                 mpi_communicator,
 *                                                 locally_relevant_dofs);
 *   
 *     left_system_matrix_T.reinit  (locally_owned_dofs,locally_owned_dofs,dsp_T,mpi_communicator);
 *     right_system_matrix_T.reinit (locally_owned_dofs,locally_owned_dofs,dsp_T,mpi_communicator);
 *     system_matrix_T.reinit       (locally_owned_dofs,locally_owned_dofs,dsp_T,mpi_communicator);
 *   
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaserHeatingassemble_system_matrix"></a> 
 * <h4>LaserHeating::assemble_system_matrix</h4>
 * 
 
 * 
 * 
 * @code
 *   template <int dim>
 *   void LaserHeating<dim>::assemble_system_matrix_init (double time_step)
 *   {
 *   
 *     TimerOutput::Scope t(computing_timer,"assemble_system_matrix_init()");
 *   
 *     QGauss<dim>  quadrature_formula(2);
 *   
 *   
 *     const InitialValues<dim>   initial_value_func_T;
 *   
 *   
 *     const RhoC<dim>       rho_C_fun_T;
 *     const K_T<dim>        k_fun_T;
 *    
 *   
 *     FEValues<dim> fe_values (fe, quadrature_formula,
 *                              update_values   | update_gradients |
 *                              update_quadrature_points | update_JxW_values);
 *   
 *     const unsigned int   dofs_per_cell = fe.dofs_per_cell;
 *     const unsigned int   n_q_points    = quadrature_formula.size();
 *   
 *   
 *     FullMatrix<double>   local_init_matrix    (dofs_per_cell, dofs_per_cell);
 *     Vector<double>       local_init_T_rhs     (dofs_per_cell);
 *   
 *     FullMatrix<double>   local_rho_c_T_matrix (dofs_per_cell, dofs_per_cell);
 *     FullMatrix<double>   local_k_T_matrix     (dofs_per_cell, dofs_per_cell);
 *   
 *   
 *     std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
 *   
 *     system_matrix_T = 0;
 *     left_system_matrix_T = 0;
 *     right_system_matrix_T = 0;
 *   
 *     system_rhs_T = 0;
 *   
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *         if(cell->is_locally_owned())
 *       {
 *   
 *         fe_values.reinit (cell);
 *         local_init_matrix = 0;
 *   
 *         local_rho_c_T_matrix = 0;
 *         local_k_T_matrix = 0;
 *   
 *         local_init_T_rhs = 0;
 *   
 *         for (unsigned int q=0; q<n_q_points; ++q)
 *           for (unsigned int i=0; i<dofs_per_cell; ++i)
 *             {
 *               const Tensor<1,dim>     div_phi_i_u = fe_values.shape_grad (i,q);
 *               const double            phi_i_u = fe_values.shape_value (i,q);
 *   
 *               for (unsigned int j=0; j<dofs_per_cell; ++j)
 *               {
 *   
 *                 const Tensor<1,dim>       div_phi_j_u = fe_values.shape_grad(j,q);
 *                 const double              phi_j_u = fe_values.shape_value (j,q);
 *   
 *                 local_init_matrix(i,j) +=    (phi_i_u *
 *                                               phi_j_u *
 *                                               fe_values.JxW (q));
 *   
 *                 local_rho_c_T_matrix(i,j) += (rho_C_fun_T.value(fe_values.quadrature_point(q)) *
 *                                               phi_i_u *
 *                                               phi_j_u *
 *                                               fe_values.JxW (q))
 *                                               +
 *                                               time_step * (theta) *  
 *                                               (k_fun_T.value(fe_values.quadrature_point(q)) *
 *                                               div_phi_i_u *
 *                                               div_phi_j_u *
 *                                               fe_values.JxW (q));
 *   
 *                 local_k_T_matrix(i,j) += (rho_C_fun_T.value(fe_values.quadrature_point(q)) *
 *                                               phi_i_u *
 *                                               phi_j_u *
 *                                               fe_values.JxW (q))
 *                                               -
 *                                               time_step * (1.0-theta) * 
 *                                               (k_fun_T.value(fe_values.quadrature_point(q)) *
 *                                               div_phi_i_u *
 *                                               div_phi_j_u *
 *                                               fe_values.JxW (q));
 *                
 *               }
 *   
 *               local_init_T_rhs(i) += (phi_i_u *
 *                                       initial_value_func_T.value (fe_values.quadrature_point (q)) *
 *                                       fe_values.JxW (q));
 *   
 *             }
 *   
 *         cell->get_dof_indices (local_dof_indices);
 *   
 * @endcode
 * 
 * copy to system_matrix_T and system_rhs_T for projecting initial values
 * 
 * @code
 *         constraints_T.distribute_local_to_global(local_init_matrix,
 *                                                local_init_T_rhs,
 *                                                local_dof_indices,
 *                                                system_matrix_T,
 *                                                system_rhs_T);
 *   
 *   
 * @endcode
 * 
 * store M + dt*theta*A as the left_system_matrix
 * 
 
 * 
 * 
 * @code
 *         constraints_T.distribute_local_to_global(local_rho_c_T_matrix,
 *                                                local_dof_indices,
 *                                                left_system_matrix_T);
 *   
 * @endcode
 * 
 * store M - dt*(1-theta)*A as the right_system_matrix
 * 
 * @code
 *         constraints_T.distribute_local_to_global(local_k_T_matrix,
 *                                                local_dof_indices,
 *                                                right_system_matrix_T);
 *   
 *   
 *       }
 *   
 *     system_matrix_T.compress(VectorOperation::add);
 *     left_system_matrix_T.compress(VectorOperation::add);
 *     right_system_matrix_T.compress(VectorOperation::add);
 *     system_rhs_T.compress(VectorOperation::add);
 *   
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaserHeatingdynamic_assemble_rhs_T"></a> 
 * <h4>LaserHeating::dynamic_assemble_rhs_T</h4>
 * The right hand side is assembled each time during running, which is necessary as
 * the laser source is moving. To separate the heat source and the right hand 
 * side assembling, the right hand side function is defined as RightHandside<dim>.
 * 
 * @code
 *   template <int dim>
 *   void LaserHeating<dim>::dynamic_assemble_rhs_T (double time, double time_step)
 *   {
 *   
 *     TimerOutput::Scope t(computing_timer,"assemble_rhs_T()");
 *   
 *     QGauss<dim>  quadrature_formula(2);
 *   
 *     RightHandside<dim> rhs_func_T_1;
 *     rhs_func_T_1.set_time(time);
 *   
 *     RightHandside<dim> rhs_func_T_2;
 *     rhs_func_T_2.set_time(time-time_step);
 *   
 *     FEValues<dim> fe_values (fe, quadrature_formula,
 *                              update_values   |
 *                              update_quadrature_points | update_JxW_values);
 *   
 *   
 *     const unsigned int   dofs_per_cell = fe.dofs_per_cell;
 *     const unsigned int   n_q_points    = quadrature_formula.size();
 *   
 *     Vector<double>       local_rhs_vector_T (dofs_per_cell);
 *   
 *     std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
 *   
 *     std::vector<double>   Rnp_cal_assemble (n_q_points);
 *   
 *   
 *     dynamic_rhs_T = 0 ;
 *   
 *     typename DoFHandler<dim>::active_cell_iterator
 *     cell = dof_handler.begin_active(),
 *     endc = dof_handler.end();
 *   
 *     for (; cell!=endc; ++cell)
 *         if(cell->is_locally_owned())
 *       {
 *         fe_values.reinit (cell);
 *         local_rhs_vector_T = 0;
 *   
 *         for (unsigned int q=0; q<n_q_points; ++q)
 *           for (unsigned int i=0; i<dofs_per_cell; ++i)
 *             {
 *               const double            phi_i_u = fe_values.shape_value (i,q);
 *   
 *   
 *               local_rhs_vector_T(i) +=     time_step * theta *
 *                                            (phi_i_u *
 *                                            rhs_func_T_1.value_v2 (fe_values.quadrature_point (q)) *
 *                                            fe_values.JxW (q))
 *                                            +
 *                                            time_step * (1.0 - theta) *
 *                                            (phi_i_u *
 *                                            rhs_func_T_2.value_v2 (fe_values.quadrature_point (q)) *
 *                                            fe_values.JxW (q));
 *             }
 *   
 *         cell->get_dof_indices (local_dof_indices);
 *   
 *   
 *         constraints_T.distribute_local_to_global(local_rhs_vector_T,
 *                                                local_dof_indices,
 *                                                dynamic_rhs_T);
 *   
 *       }
 *   
 *     dynamic_rhs_T.compress(VectorOperation::add);
 *   
 *   
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaserHeatingsolve"></a> 
 * <h4>LaserHeating::solve</h4>
 * Solving the equation is direct. Recall that we have defined several matrices and vectors,
 * to avoid ambiguous, here, we only use system_matrix_T as the system matrix, system_rhs_T
 * as the system right hand side. The vector completely_distributed_solution is used to store
 * the obtained solution.
 * 
 * @code
 *   template <int dim>
 *   void LaserHeating<dim>::solve_T ()
 *   {
 *     TimerOutput::Scope t(computing_timer,"solve_T()");
 *   
 *     TrilinosWrappers::MPI::Vector   completely_distributed_solution (locally_owned_dofs,mpi_communicator);
 *     SolverControl     solver_control (1*system_rhs_T.size(),1e-12*system_rhs_T.l2_norm(),true);
 *   
 *     TrilinosWrappers::SolverCG      solver (solver_control);
 *   
 *     TrilinosWrappers::PreconditionAMG  preconditioner;
 *     TrilinosWrappers::PreconditionAMG::AdditionalData  data;
 *   
 *     preconditioner.initialize(system_matrix_T,data);
 *   
 *     solver.solve (system_matrix_T,completely_distributed_solution,system_rhs_T,preconditioner);
 *   
 *   
 * @endcode
 * 
 * Print the number of iterations by hand.
 * 
 
 * 
 * 
 * @code
 *     pcout     << "   " << solver_control.last_step()
 *               << " CG iterations needed to obtain convergence." << std::endl
 *               << "\t initial convergence value = " << solver_control.initial_value() << std::endl
 *               << "\t final convergence value = " << solver_control.last_value() << std::endl
 *               << std::endl;
 *   
 *     constraints_T.distribute (completely_distributed_solution);
 *     new_solution_T = completely_distributed_solution;
 *   
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaserHeatingrefine_mesh"></a> 
 * <h4>LaserHeating::refine_mesh</h4>
 * 
 
 * 
 * 
 * @code
 *   template <int dim>
 *   void LaserHeating<dim>::refine_mesh()
 *   {
 *   
 *     TimerOutput::Scope t(computing_timer,"refine_mesh_at_beginning");
 *   
 *   
 *     QGauss<dim>  quadrature_formula(2);
 *   
 *     FEValues<dim> fe_values (fe, quadrature_formula,update_quadrature_points);
 *   
 *   
 * @endcode
 * 
 * only refine mesh 5um above the TiO2 and glass interface 
 * 
 
 * 
 * 
 * @code
 *     for (typename Triangulation<dim>::active_cell_iterator
 *             cell = triangulation.begin_active();
 *             cell != triangulation.end(); ++cell)
 *         if(cell->is_locally_owned())
 *         {
 *             fe_values.reinit(cell);
 *             if(std::abs(fe_values.quadrature_point(0)[1]) <= global_film_thickness+5e-6)
 *             {
 *                 cell->set_refine_flag();
 *             }
 *             else
 *             {
 *                 cell->clear_refine_flag();
 *             }
 *         }
 *     triangulation.execute_coarsening_and_refinement();
 *     
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaserHeatinghoutput_results"></a> 
 * <h4>LaserHeatingh::output_results</h4>
 * 
 
 * 
 * 
 * @code
 *   template <int dim>
 *   void LaserHeating<dim>::output_results (int output_num) const
 *    {
 *   
 *      DataOut<dim> data_out;
 *   
 *      data_out.attach_dof_handler (dof_handler);
 *      data_out.add_data_vector (old_solution_T,         "T");
 *   
 * @endcode
 * 
 * the length of output numbering
 * 
 
 * 
 * 
 * @code
 *      int step_N = 7;
 *   
 *      data_out.build_patches ();
 *   
 *      const std::string filename = ("solution-" +
 *                                   Utilities::int_to_string (output_num, step_N) +
 *                                   "." +
 *                                   Utilities::int_to_string (triangulation.locally_owned_subdomain(),4) +
 *                                   ".vtu");
 *      std::ofstream output (filename.c_str());
 *      data_out.write_vtu (output);
 *   
 *   
 * @endcode
 * 
 * output the overall solution
 * 
 
 * 
 * 
 * @code
 *      if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
 *      {
 *          std::vector<std::string> filenames;
 *          for (unsigned int i=0;i<Utilities::MPI::n_mpi_processes(mpi_communicator);++i)
 *              filenames.push_back ("solution-" +
 *                                   Utilities::int_to_string (output_num, step_N) +
 *                                   "." +
 *                                   Utilities::int_to_string (i,4) +
 *                                   ".vtu");
 *          std::ofstream master_output (("solution-" + 
 *                                       Utilities::int_to_string (output_num,step_N)+
 *                                       ".pvtu").c_str());
 *          data_out.write_pvtu_record (master_output,filenames);
 *   
 *      }
 *   
 *    }
 *   
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaserHeatingrun"></a> 
 * <h4>LaserHeating::run</h4>
 * 
 
 * 
 * This is the function which has the top-level control over everything. Apart
 * from one line of additional output, it is the same as for the previous
 * example.
 * 
 * @code
 *   template <int dim>
 *   void LaserHeating<dim>::run ()
 *   {
 *       pcout << "Solving problem in " << dim << " space dimensions." << std::endl;
 *   
 *       
 *       make_grid();
 *       refine_mesh();
 *       setup_system ();
 *       assemble_system_matrix_init (global_simulation_time_step);
 *   
 * @endcode
 * 
 * projection of initial conditions by solving.
 * solution stored in new_solution_T; 
 * 
 
 * 
 * 
 * @code
 *       solve_T ();
 *   
 *       old_solution_T        = new_solution_T;
 *       old_solution_T_cal    = new_solution_T;
 *   
 * @endcode
 * 
 * reinitialization
 * 
 
 * 
 * 
 * @code
 *       system_matrix_T = 0;
 *       system_rhs_T = 0;
 *   
 *   
 *       double time_step = global_simulation_time_step;
 *       double time = 0;
 *       int timestep_number = 0;
 *   
 * @endcode
 * 
 * output initial values; need ghost cells 
 * 
 * @code
 *       output_results (0);
 *   
 *   
 *         while(time < global_simulation_end_time)
 *         {
 *   
 *           time += time_step;
 *           timestep_number ++;
 *   
 *           pcout << "Time step " << timestep_number
 *                     << " at t=" << time
 *                     << " time_step = " << time_step
 *                     << std::endl;
 *   
 * @endcode
 * 
 * the dynamic solving part
 * 
 * @code
 *           {
 *   
 *               right_system_matrix_T.vmult(system_rhs_T,old_solution_T_cal);
 *   
 *               dynamic_assemble_rhs_T (time,time_step);
 *               system_rhs_T.add(1,dynamic_rhs_T);
 *   
 *               system_matrix_T.copy_from (left_system_matrix_T);
 *   
 *               {
 *                   BoundaryValues<dim> boundary_values_function;
 *                   std::map<types::global_dof_index,double> boundary_values;
 *   
 *                   VectorTools::interpolate_boundary_values (dof_handler,
 *                                                               BOUNDARY_NUM,
 *                                                               boundary_values_function,
 *                                                               boundary_values);
 *   
 *                   MatrixTools::apply_boundary_values (boundary_values,
 *                                                           system_matrix_T,
 *                                                           new_solution_T,
 *                                                           system_rhs_T,
 *                                                           false);
 *               }
 *   
 *   
 *               solve_T ();
 *   
 * @endcode
 * 
 * old_solution_T is used for output, holding ghost cells
 * old_solution_T_cal is used for calculation, holding only
 * locally owned cells.
 * 
 * @code
 *               old_solution_T          = new_solution_T;
 *               old_solution_T_cal      = new_solution_T;
 *   
 *               if (Utilities::MPI::n_mpi_processes(mpi_communicator) <= 96 && (timestep_number % 50  == 0 ))
 *               {
 *                   TimerOutput::Scope t(computing_timer,"output");
 *                   output_results (timestep_number);
 *               }
 *   
 *               computing_timer.print_summary ();
 *               computing_timer.reset();
 *               
 *               pcout << std::endl;
 *   
 *           }
 *         }
 *   
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * 
 * @code
 *   int main (int argc, char *argv[])
 *   {
 *       try
 *       {
 *         Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
 *   
 *         LaserHeating<2> laserHeating_2d;
 *         laserHeating_2d.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;
 *   }
 * @endcode
 
 
<a name="ann-boundaryInit.h"></a>
<h1>Annotated version of boundaryInit.h</h1>
 * 
 * 
 * 
 * 
 
 * 
 * #----------------------------------------------------------
 * #
 * # This file defines the boundary and initial conditions
 * #
 * #----------------------------------------------------------
 * 
 
 * 
 * 
 * @code
 *   #ifndef GLOBAL_PARA
 *   #define GLOBAL_PARA
 *   #include "./globalPara.h"
 *   #endif
 *   
 * @endcode
 * 
 * #----------------------------------------------------------
 * # Declaration
 * 
 
 * 
 * 
 * @code
 *   template <int dim>
 *   class BoundaryValues : public Function<dim>
 *   {
 *   public:
 *     BoundaryValues () : Function<dim>() {}
 *   
 *     virtual double value (const Point<dim>   &p,
 *                           const unsigned int  component = 0) const override;
 *   };
 *   
 *   template <int dim>
 *   class InitialValues : public Function<dim>
 *   {
 *   public:
 *     InitialValues () : Function<dim>() {}
 *   
 *     virtual double value (const Point<dim>   &p,
 *                           const unsigned int  component = 0) const override;
 *   };
 *   
 *   
 * @endcode
 * 
 * #----------------------------------------------------------
 * # Implementation
 * 
 
 * 
 * 
 * @code
 *   template <int dim>
 *   double BoundaryValues<dim>::value (const Point<dim> &/*p*/,
 *                                      const unsigned int /*component*/) const
 *   {
 *     return 293;
 *   }
 *   
 *   
 *   template <int dim>
 *   double InitialValues<dim>::value (const Point<dim> &/*p*/,
 *                                      const unsigned int /*component*/) const
 *   {
 *     return 293;
 *   }
 *   
 *   
 *   
 *   
 *   
 * @endcode
 * 
 * #----------------------------------------------------------
 * # Declaration and Implementation
 * # mass density and heat capacity
 * 
 
 * 
 * 
 * @code
 *   template <int dim>
 *   class RhoC : public Function<dim>
 *   {
 *   public:
 *     RhoC () : Function<dim>() {}
 *   
 *     virtual double value (const Point<dim>   &p,
 *                           const unsigned int  component = 0) const override;
 *   };
 *   
 *   template <int dim>
 *   double RhoC<dim>::value (const Point<dim> &p,
 *                                      const unsigned int /*component*/) const
 *   {
 * @endcode
 * 
 * # p stores the xyz coordinates at each vertex
 * # for 2D problems in xy, we assume the non-uniform is in y-axis.
 * 
 * @code
 *     if ( p[1] >= -global_film_thickness )
 *     {
 *         return global_rho_Tio2 * global_C_Tio2;
 *     }
 *     else
 *         return global_rho_glass * global_C_glass;
 *   }
 *   
 *   
 *   
 *   
 * @endcode
 * 
 * #----------------------------------------------------------
 * # Declaration and Implementation
 * # thermal conductivity
 * 
 
 * 
 * 
 * @code
 *   template <int dim>
 *   class K_T : public Function<dim>
 *   {
 *   public:
 *     K_T () : Function<dim>() {}
 *   
 *     virtual double value (const Point<dim>   &p,
 *                           const unsigned int  component = 0) const override;
 *   };
 *   
 *   template <int dim>
 *   double K_T<dim>::value (const Point<dim> &p,
 *                                      const unsigned int /*component*/) const
 *   {
 * @endcode
 * 
 * # p stores the xyz coordinates at each vertex
 * # for 2D problems in xy, we assume the non-uniform is in y-axis.
 * 
 * @code
 *     if ( p[1] >= -global_film_thickness)
 *     {
 *         return global_k_Tio2;
 *     }
 *     else
 *         return global_k_glass;
 *   }
 *   
 * @endcode
 
 
<a name="ann-globalPara.h"></a>
<h1>Annotated version of globalPara.h</h1>
 * 
 * 
 * 
 * 
 
 * 
 * # physics constants
 * 
 * @code
 *   double global_PI = 3.1415927;
 *   
 * @endcode
 * 
 * # Laser
 * 
 * @code
 *   double global_Pow_laser = 0.4;      // power [W]
 *   double global_spotsize_at_e_2 = 20e-6;  // laser spot size at e^(-2) [m]
 *   double global_c_laser = global_spotsize_at_e_2 / 4.0;   // C parameter in Gaussian func, [m]
 *   double global_c_hwhm = global_c_laser * 2.35482 / 2;    // HWHM, [m]
 *   double global_V_scan_x = 10e-3;   // scan speed, [m/s]
 *   
 *   double global_init_position_x0 = -50e-6;    // initial spot center position
 *   
 * @endcode
 * 
 * # material
 * thin film
 * 
 * @code
 *   double global_rho_Tio2 = 4200;      // mass density, [kg/m^3]
 *   double global_C_Tio2 = 690;         // heat capacity, [J/kg/K]
 *   double global_k_Tio2 = 4.8;         // thermal conductivity, [W/m/K]
 * @endcode
 * 
 * substrate
 * 
 * @code
 *   double global_rho_glass = 2200;
 *   double global_C_glass = 700;
 *   double global_k_glass = 1.8;
 *   
 *   double global_film_thickness = 400e-9;  // film thickness, [m]
 *   
 * @endcode
 * 
 * # simulation time
 * 
 * @code
 *   double global_simulation_time_step = 1e-5;          // 10 [us]
 *   double global_simulation_end_time = 100e-6 / global_V_scan_x; //    100 [um] / scan speed
 *   
 * @endcode
 * 
 * # about the MESH
 * 
 * @code
 *   #define BOUNDARY_NUM 11
 * @endcode
 
 
<a name="ann-rightHandSide.h"></a>
<h1>Annotated version of rightHandSide.h</h1>
 * 
 * 
 * 
 * #----------------------------------------------------------
 * #
 * # This file defines the boundary and initial conditions
 * #
 * #----------------------------------------------------------
 * 
 
 * 
 * 
 * @code
 *   #ifndef GLOBAL_PARA
 *   #define GLOBAL_PARA
 *   #include "./globalPara.h"
 *   #endif
 *   
 * @endcode
 * 
 * #----------------------------------------------------------
 * # Declaration
 * 
 
 * 
 * 
 * @code
 *   template <int dim>
 *   class RightHandside : public Function<dim>
 *   {
 *   public:
 *     RightHandside () : Function<dim>() {}
 *   
 *     double value_v2 (const Point<dim>   &p);
 *   };
 *   
 * @endcode
 * 
 * #----------------------------------------------------------
 * # Implementation
 * 
 
 * 
 * 
 * @code
 *   template <int dim>
 *   double RightHandside<dim>::value_v2 (const Point<dim> &p)
 *   {
 *   
 *       double alpha_abs = 1e4;
 *   
 *   
 *       if(p[1] >= -global_film_thickness)
 *       {
 *   
 *       double P00 = global_Pow_laser / global_PI / global_c_laser / global_c_laser / 2.0;
 *   
 *       double I00 = P00 * std::exp(-
 *                   ( 
 *                       (p[0] - global_V_scan_x * this->get_time()-global_init_position_x0) * 
 *                       (p[0] - global_V_scan_x * this->get_time()-global_init_position_x0)
 *                   ) /
 *                   (2.0 * global_c_laser * global_c_laser) );
 *   
 *   
 *       return alpha_abs * I00 *  std::exp(-alpha_abs*(0 - p[1]));
 *   
 *       }
 *   
 *       else
 *       {
 * @endcode
 * 
 * # no absorption
 * 
 * @code
 *           return 0;
 *       }
 *   
 *   }
 *   
 *   
 * @endcode
 
 
*/