Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
Distributed_Moving_Laser_Heating.h
Go to the documentation of this file.
1
273 *  
274 *  
275 * @endcode
276 *
277 *
278 * <a name="Includefiles"></a>
279 * <h3>Include files</h3>
280 *
281
282 *
283 * All the include files have already been discussed in previous tutorials.
284 *
285 * @code
286 *   #include <deal.II/dofs/dof_handler.h>
287 *   #include <deal.II/dofs/dof_renumbering.h>
288 *   #include <deal.II/grid/grid_generator.h>
289 *   #include <deal.II/grid/tria_accessor.h>
290 *   #include <deal.II/grid/tria_iterator.h>
291 *   #include <deal.II/dofs/dof_accessor.h>
292 *   #include <deal.II/fe/fe_q.h>
293 *   #include <deal.II/dofs/dof_tools.h>
294 *   #include <deal.II/fe/fe_values.h>
295 *   #include <deal.II/base/quadrature_lib.h>
296 *   #include <deal.II/base/function.h>
297 *   #include <deal.II/numerics/vector_tools.h>
298 *   #include <deal.II/numerics/matrix_tools.h>
299 *   #include <deal.II/lac/vector.h>
300 *   #include <deal.II/lac/full_matrix.h>
301 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
302 *  
303 *   #include <deal.II/grid/grid_in.h>
304 *   #include <deal.II/grid/grid_tools.h>
305 *   #include <deal.II/lac/trilinos_vector.h>
306 *   #include <deal.II/lac/trilinos_sparse_matrix.h>
307 *   #include <deal.II/lac/trilinos_solver.h>
308 *   #include <deal.II/lac/trilinos_precondition.h>
309 *  
310 *   #include <deal.II/lac/sparsity_tools.h>
311 *   #include <deal.II/lac/generic_linear_algebra.h>
312 *  
313 *   #include <deal.II/base/conditional_ostream.h>
314 *   #include <deal.II/base/utilities.h>
315 *   #include <deal.II/base/index_set.h>
316 *   #include <deal.II/distributed/tria.h>
317 *  
318 *   #include <deal.II/numerics/data_out.h>
319 *   #include <fstream>
320 *   #include <iostream>
321 *  
322 *  
323 *   #include <deal.II/base/logstream.h>
324 *   #include <deal.II/lac/affine_constraints.h>
325 *   #include <deal.II/base/timer.h>
326 *  
327 * @endcode
328 *
329 * grid refine
330 *
331 * @code
332 *   #include <deal.II/numerics/error_estimator.h>
333 *   #include <deal.II/distributed/grid_refinement.h>
334 *   #include <deal.II/distributed/solution_transfer.h>
335 *  
336 * @endcode
337 *
338 * remember to use the namespace dealii before
339 * defining a class that is inheritaged from Function<dim>
340 *
341 * @code
342 *   using namespace dealii;
343 *  
344 * @endcode
345 *
346 * In general, it is more clear to separate boundary and initial condictions,
347 * as well as the right hand side function from a file holding all the things.
348 * To do so, in this work, the globalPara.h file defines physical constants,
349 * laser parameters, and heat characteristics of materials involved. Boundary
350 * and initial conditions are defined in boundaryInit.h. The rightHandSide.h
351 * defines the heat source, which in this work is a moving Gaussian beam.
352 *
353
354 *
355 *
356 * @code
357 *   #ifndef GLOBAL_PARA
358 *   #define GLOBAL_PARA
359 *   #include "./globalPara.h"
360 *   #include "./boundaryInit.h"
361 *   #include "./rightHandSide.h"
362 *   #endif
363 *  
364 * @endcode
365 *
366 * Now the main class is defined as following
367 *
368 * @code
369 *   template <int dim>
370 *   class LaserHeating
371 *   {
372 *   public:
373 *   LaserHeating ();
374 *   ~LaserHeating ();
375 *   void run ();
376 *  
377 *   private:
378 *   void make_grid ();
379 *   void setup_system();
380 *  
381 *   void assemble_system_matrix_init (double time_step);
382 *   void dynamic_assemble_rhs_T (double time, double time_step);
383 *  
384 *   void solve_T ();
385 *  
386 *   void refine_mesh();
387 *   void output_results (int output_num) const;
388 *  
389 *   MPI_Comm mpi_communicator;
390 *  
392 *   FE_Q<dim> fe;
393 *   DoFHandler<dim> dof_handler;
394 *  
395 *   AffineConstraints<double> constraints_T;
396 *  
397 *  
398 * @endcode
399 *
400 * system_matrix
401 *
402 * @code
403 *   TrilinosWrappers::SparseMatrix system_matrix_T;
404 *  
405 * @endcode
406 *
407 * for storing left matrix
408 *
409 * @code
410 *   TrilinosWrappers::SparseMatrix left_system_matrix_T;
411 *  
412 * @endcode
413 *
414 * for storing right matrix
415 *
416 * @code
417 *   TrilinosWrappers::SparseMatrix right_system_matrix_T;
418 *  
419 * @endcode
420 *
421 * System_rhs, only locally owned cells
422 *
423 * @code
424 *   TrilinosWrappers::MPI::Vector system_rhs_T;
425 *  
426 * @endcode
427 *
428 * Solutions
429 * Old Solutions with ghost cells, for output
430 *
431 * @code
432 *   TrilinosWrappers::MPI::Vector old_solution_T;
433 *  
434 * @endcode
435 *
436 * Old Solutions only with locally owned cells
437 *
438 * @code
439 *   TrilinosWrappers::MPI::Vector old_solution_T_cal;
440 *  
441 * @endcode
442 *
443 * New Solutions only with locally owned cells
444 *
445 * @code
446 *   TrilinosWrappers::MPI::Vector new_solution_T;
447 *  
448 * @endcode
449 *
450 * Dynamic assembling of the righthandside terms
451 *
452 * @code
453 *   TrilinosWrappers::MPI::Vector dynamic_rhs_T;
454 *  
455 *  
456 *   IndexSet locally_owned_dofs;
457 *   IndexSet locally_relevant_dofs;
458 *  
459 *   ConditionalOStream pcout;
460 *   TimerOutput computing_timer;
461 *  
462 *   double theta;
463 *  
464 *  
465 *  
466 *   };
467 *  
468 * @endcode
469 *
470 * the constructor
471 *
472 * @code
473 *   template <int dim>
474 *   LaserHeating<dim>::LaserHeating ()
475 *   :
476 *   mpi_communicator (MPI_COMM_WORLD),
477 *   triangulation (mpi_communicator),
478 *   fe (1),
479 *   dof_handler (triangulation),
480 *   pcout (std::cout,Utilities::MPI::this_mpi_process(mpi_communicator) == 0),
481 *   computing_timer (mpi_communicator, pcout, TimerOutput::summary, TimerOutput::wall_times),
482 *   theta(0.5)
483 *   {}
484 *  
485 * @endcode
486 *
487 * the destructor
488 *
489 * @code
490 *   template <int dim>
491 *   LaserHeating<dim>::~LaserHeating ()
492 *   {
493 *   dof_handler.clear();
494 *   }
495 *  
496 * @endcode
497 *
498 *
499 * <a name="LaserHeatingmake_grid"></a>
500 * <h4>LaserHeating::make_grid</h4>
501 * make grid by importing msh file, and rescale
502 *
503 * @code
504 *   template <int dim>
505 *   void LaserHeating<dim>::make_grid ()
506 *   {
507 *   TimerOutput::Scope t(computing_timer,"make_grid()");
508 *  
509 *   GridIn<dim> grid_in;
510 *   grid_in.attach_triangulation (triangulation);
511 *   std::ifstream input_file ("geometry.msh");
512 *   grid_in.read_msh (input_file);
514 *  
515 *   pcout << " Number of active cells: "
516 *   << triangulation.n_global_active_cells()
517 *   << std::endl
518 *   << " Total number of cells: "
519 *   << triangulation.n_cells()
520 *   << std::endl;
521 *   }
522 *  
523 * @endcode
524 *
525 *
526 * <a name="LaserHeatingsetup_system"></a>
527 * <h4>LaserHeating::setup_system</h4>
528 * initialization
529 *
530 * @code
531 *   template <int dim>
532 *   void LaserHeating<dim>::setup_system ()
533 *   {
534 *  
535 *   TimerOutput::Scope t(computing_timer,"setup_system()");
536 *  
537 *   dof_handler.distribute_dofs (fe);
538 *  
539 *   pcout << " Number of degrees of freedom: "
540 *   << dof_handler.n_dofs()
541 *   << std::endl;
542 *  
543 *   locally_owned_dofs = dof_handler.locally_owned_dofs();
544 *   DoFTools::extract_locally_relevant_dofs (dof_handler, locally_relevant_dofs);
545 *  
546 * @endcode
547 *
548 * we want to output solution, so here should have ghost cells
549 *
550 * @code
551 *   old_solution_T.reinit (locally_owned_dofs,locally_relevant_dofs,mpi_communicator);
552 *  
553 * @endcode
554 *
555 * locally owned cells
556 *
557 * @code
558 *   old_solution_T_cal.reinit (locally_owned_dofs,mpi_communicator);
559 *   new_solution_T.reinit (locally_owned_dofs,mpi_communicator);
560 *   dynamic_rhs_T.reinit (locally_owned_dofs,mpi_communicator);
561 *   system_rhs_T.reinit (locally_owned_dofs,mpi_communicator);
562 *  
563 *   constraints_T.clear();
564 *   constraints_T.reinit (locally_relevant_dofs);
565 *   DoFTools::make_hanging_node_constraints (dof_handler, constraints_T);
566 *   constraints_T.close();
567 *  
568 *   const unsigned int myid = Utilities::MPI::this_mpi_process (mpi_communicator);
569 *   DynamicSparsityPattern dsp_T(locally_relevant_dofs);
570 *   DoFTools::make_sparsity_pattern (dof_handler, dsp_T, constraints_T, false, myid);
571 *  
573 *   locally_owned_dofs,
574 *   mpi_communicator,
575 *   locally_relevant_dofs);
576 *  
577 *   left_system_matrix_T.reinit (locally_owned_dofs,locally_owned_dofs,dsp_T,mpi_communicator);
578 *   right_system_matrix_T.reinit (locally_owned_dofs,locally_owned_dofs,dsp_T,mpi_communicator);
579 *   system_matrix_T.reinit (locally_owned_dofs,locally_owned_dofs,dsp_T,mpi_communicator);
580 *  
581 *   }
582 *  
583 *  
584 * @endcode
585 *
586 *
587 * <a name="LaserHeatingassemble_system_matrix"></a>
588 * <h4>LaserHeating::assemble_system_matrix</h4>
589 *
590
591 *
592 *
593 * @code
594 *   template <int dim>
595 *   void LaserHeating<dim>::assemble_system_matrix_init (double time_step)
596 *   {
597 *  
598 *   TimerOutput::Scope t(computing_timer,"assemble_system_matrix_init()");
599 *  
600 *   QGauss<dim> quadrature_formula(2);
601 *  
602 *  
603 *   const InitialValues<dim> initial_value_func_T;
604 *  
605 *  
606 *   const RhoC<dim> rho_C_fun_T;
607 *   const K_T<dim> k_fun_T;
608 *  
609 *  
610 *   FEValues<dim> fe_values (fe, quadrature_formula,
613 *  
614 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
615 *   const unsigned int n_q_points = quadrature_formula.size();
616 *  
617 *  
618 *   FullMatrix<double> local_init_matrix (dofs_per_cell, dofs_per_cell);
619 *   Vector<double> local_init_T_rhs (dofs_per_cell);
620 *  
621 *   FullMatrix<double> local_rho_c_T_matrix (dofs_per_cell, dofs_per_cell);
622 *   FullMatrix<double> local_k_T_matrix (dofs_per_cell, dofs_per_cell);
623 *  
624 *  
625 *   std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
626 *  
627 *   system_matrix_T = 0;
628 *   left_system_matrix_T = 0;
629 *   right_system_matrix_T = 0;
630 *  
631 *   system_rhs_T = 0;
632 *  
633 *   for (const auto &cell : dof_handler.active_cell_iterators())
634 *   if(cell->is_locally_owned())
635 *   {
636 *  
637 *   fe_values.reinit (cell);
638 *   local_init_matrix = 0;
639 *  
640 *   local_rho_c_T_matrix = 0;
641 *   local_k_T_matrix = 0;
642 *  
643 *   local_init_T_rhs = 0;
644 *  
645 *   for (unsigned int q=0; q<n_q_points; ++q)
646 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
647 *   {
648 *   const Tensor<1,dim> div_phi_i_u = fe_values.shape_grad (i,q);
649 *   const double phi_i_u = fe_values.shape_value (i,q);
650 *  
651 *   for (unsigned int j=0; j<dofs_per_cell; ++j)
652 *   {
653 *  
654 *   const Tensor<1,dim> div_phi_j_u = fe_values.shape_grad(j,q);
655 *   const double phi_j_u = fe_values.shape_value (j,q);
656 *  
657 *   local_init_matrix(i,j) += (phi_i_u *
658 *   phi_j_u *
659 *   fe_values.JxW (q));
660 *  
661 *   local_rho_c_T_matrix(i,j) += (rho_C_fun_T.value(fe_values.quadrature_point(q)) *
662 *   phi_i_u *
663 *   phi_j_u *
664 *   fe_values.JxW (q))
665 *   +
666 *   time_step * (theta) *
667 *   (k_fun_T.value(fe_values.quadrature_point(q)) *
668 *   div_phi_i_u *
669 *   div_phi_j_u *
670 *   fe_values.JxW (q));
671 *  
672 *   local_k_T_matrix(i,j) += (rho_C_fun_T.value(fe_values.quadrature_point(q)) *
673 *   phi_i_u *
674 *   phi_j_u *
675 *   fe_values.JxW (q))
676 *   -
677 *   time_step * (1.0-theta) *
678 *   (k_fun_T.value(fe_values.quadrature_point(q)) *
679 *   div_phi_i_u *
680 *   div_phi_j_u *
681 *   fe_values.JxW (q));
682 *  
683 *   }
684 *  
685 *   local_init_T_rhs(i) += (phi_i_u *
686 *   initial_value_func_T.value (fe_values.quadrature_point (q)) *
687 *   fe_values.JxW (q));
688 *  
689 *   }
690 *  
691 *   cell->get_dof_indices (local_dof_indices);
692 *  
693 * @endcode
694 *
695 * copy to system_matrix_T and system_rhs_T for projecting initial values
696 *
697 * @code
698 *   constraints_T.distribute_local_to_global(local_init_matrix,
699 *   local_init_T_rhs,
700 *   local_dof_indices,
701 *   system_matrix_T,
702 *   system_rhs_T);
703 *  
704 *  
705 * @endcode
706 *
707 * store M + dt*theta*A as the left_system_matrix
708 *
709
710 *
711 *
712 * @code
713 *   constraints_T.distribute_local_to_global(local_rho_c_T_matrix,
714 *   local_dof_indices,
715 *   left_system_matrix_T);
716 *  
717 * @endcode
718 *
719 * store M - dt*(1-theta)*A as the right_system_matrix
720 *
721 * @code
722 *   constraints_T.distribute_local_to_global(local_k_T_matrix,
723 *   local_dof_indices,
724 *   right_system_matrix_T);
725 *  
726 *  
727 *   }
728 *  
729 *   system_matrix_T.compress(VectorOperation::add);
730 *   left_system_matrix_T.compress(VectorOperation::add);
731 *   right_system_matrix_T.compress(VectorOperation::add);
732 *   system_rhs_T.compress(VectorOperation::add);
733 *  
734 *   }
735 *  
736 *  
737 * @endcode
738 *
739 *
740 * <a name="LaserHeatingdynamic_assemble_rhs_T"></a>
741 * <h4>LaserHeating::dynamic_assemble_rhs_T</h4>
742 * The right hand side is assembled each time during running, which is necessary as
743 * the laser source is moving. To separate the heat source and the right hand
744 * side assembling, the right hand side function is defined as RightHandside<dim>.
745 *
746 * @code
747 *   template <int dim>
748 *   void LaserHeating<dim>::dynamic_assemble_rhs_T (double time, double time_step)
749 *   {
750 *  
751 *   TimerOutput::Scope t(computing_timer,"assemble_rhs_T()");
752 *  
753 *   QGauss<dim> quadrature_formula(2);
754 *  
755 *   RightHandside<dim> rhs_func_T_1;
756 *   rhs_func_T_1.set_time(time);
757 *  
758 *   RightHandside<dim> rhs_func_T_2;
759 *   rhs_func_T_2.set_time(time-time_step);
760 *  
761 *   FEValues<dim> fe_values (fe, quadrature_formula,
762 *   update_values |
764 *  
765 *  
766 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
767 *   const unsigned int n_q_points = quadrature_formula.size();
768 *  
769 *   Vector<double> local_rhs_vector_T (dofs_per_cell);
770 *  
771 *   std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
772 *  
773 *   std::vector<double> Rnp_cal_assemble (n_q_points);
774 *  
775 *  
776 *   dynamic_rhs_T = 0 ;
777 *  
779 *   cell = dof_handler.begin_active(),
780 *   endc = dof_handler.end();
781 *  
782 *   for (; cell!=endc; ++cell)
783 *   if(cell->is_locally_owned())
784 *   {
785 *   fe_values.reinit (cell);
786 *   local_rhs_vector_T = 0;
787 *  
788 *   for (unsigned int q=0; q<n_q_points; ++q)
789 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
790 *   {
791 *   const double phi_i_u = fe_values.shape_value (i,q);
792 *  
793 *  
794 *   local_rhs_vector_T(i) += time_step * theta *
795 *   (phi_i_u *
796 *   rhs_func_T_1.value_v2 (fe_values.quadrature_point (q)) *
797 *   fe_values.JxW (q))
798 *   +
799 *   time_step * (1.0 - theta) *
800 *   (phi_i_u *
801 *   rhs_func_T_2.value_v2 (fe_values.quadrature_point (q)) *
802 *   fe_values.JxW (q));
803 *   }
804 *  
805 *   cell->get_dof_indices (local_dof_indices);
806 *  
807 *  
808 *   constraints_T.distribute_local_to_global(local_rhs_vector_T,
809 *   local_dof_indices,
810 *   dynamic_rhs_T);
811 *  
812 *   }
813 *  
814 *   dynamic_rhs_T.compress(VectorOperation::add);
815 *  
816 *  
817 *   }
818 *  
819 *  
820 * @endcode
821 *
822 *
823 * <a name="LaserHeatingsolve"></a>
824 * <h4>LaserHeating::solve</h4>
825 * Solving the equation is direct. Recall that we have defined several matrices and vectors,
826 * to avoid ambiguous, here, we only use system_matrix_T as the system matrix, system_rhs_T
827 * as the system right hand side. The vector completely_distributed_solution is used to store
828 * the obtained solution.
829 *
830 * @code
831 *   template <int dim>
832 *   void LaserHeating<dim>::solve_T ()
833 *   {
834 *   TimerOutput::Scope t(computing_timer,"solve_T()");
835 *  
836 *   TrilinosWrappers::MPI::Vector completely_distributed_solution (locally_owned_dofs,mpi_communicator);
837 *   SolverControl solver_control (1*system_rhs_T.size(),1e-12*system_rhs_T.l2_norm(),true);
838 *  
839 *   TrilinosWrappers::SolverCG solver (solver_control);
840 *  
841 *   TrilinosWrappers::PreconditionAMG preconditioner;
843 *  
844 *   preconditioner.initialize(system_matrix_T,data);
845 *  
846 *   solver.solve (system_matrix_T,completely_distributed_solution,system_rhs_T,preconditioner);
847 *  
848 *  
849 * @endcode
850 *
851 * Print the number of iterations by hand.
852 *
853
854 *
855 *
856 * @code
857 *   pcout << " " << solver_control.last_step()
858 *   << " CG iterations needed to obtain convergence." << std::endl
859 *   << "\t initial convergence value = " << solver_control.initial_value() << std::endl
860 *   << "\t final convergence value = " << solver_control.last_value() << std::endl
861 *   << std::endl;
862 *  
863 *   constraints_T.distribute (completely_distributed_solution);
864 *   new_solution_T = completely_distributed_solution;
865 *  
866 *   }
867 *  
868 *  
869 * @endcode
870 *
871 *
872 * <a name="LaserHeatingrefine_mesh"></a>
873 * <h4>LaserHeating::refine_mesh</h4>
874 *
875
876 *
877 *
878 * @code
879 *   template <int dim>
880 *   void LaserHeating<dim>::refine_mesh()
881 *   {
882 *  
883 *   TimerOutput::Scope t(computing_timer,"refine_mesh_at_beginning");
884 *  
885 *  
886 *   QGauss<dim> quadrature_formula(2);
887 *  
888 *   FEValues<dim> fe_values (fe, quadrature_formula,update_quadrature_points);
889 *  
890 *  
891 * @endcode
892 *
893 * only refine mesh 5um above the TiO2 and glass interface
894 *
895
896 *
897 *
898 * @code
900 *   cell = triangulation.begin_active();
901 *   cell != triangulation.end(); ++cell)
902 *   if(cell->is_locally_owned())
903 *   {
904 *   fe_values.reinit(cell);
905 *   if(std::abs(fe_values.quadrature_point(0)[1]) <= global_film_thickness+5e-6)
906 *   {
907 *   cell->set_refine_flag();
908 *   }
909 *   else
910 *   {
911 *   cell->clear_refine_flag();
912 *   }
913 *   }
914 *   triangulation.execute_coarsening_and_refinement();
915 *  
916 *   }
917 *  
918 *  
919 * @endcode
920 *
921 *
922 * <a name="LaserHeatinghoutput_results"></a>
923 * <h4>LaserHeatingh::output_results</h4>
924 *
925
926 *
927 *
928 * @code
929 *   template <int dim>
930 *   void LaserHeating<dim>::output_results (int output_num) const
931 *   {
932 *  
933 *   DataOut<dim> data_out;
934 *  
935 *   data_out.attach_dof_handler (dof_handler);
936 *   data_out.add_data_vector (old_solution_T, "T");
937 *  
938 * @endcode
939 *
940 * the length of output numbering
941 *
942
943 *
944 *
945 * @code
946 *   int step_N = 7;
947 *  
948 *   data_out.build_patches ();
949 *  
950 *   const std::string filename = ("solution-" +
951 *   Utilities::int_to_string (output_num, step_N) +
952 *   "." +
953 *   Utilities::int_to_string (triangulation.locally_owned_subdomain(),4) +
954 *   ".vtu");
955 *   std::ofstream output (filename.c_str());
956 *   data_out.write_vtu (output);
957 *  
958 *  
959 * @endcode
960 *
961 * output the overall solution
962 *
963
964 *
965 *
966 * @code
967 *   if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
968 *   {
969 *   std::vector<std::string> filenames;
970 *   for (unsigned int i=0;i<Utilities::MPI::n_mpi_processes(mpi_communicator);++i)
971 *   filenames.push_back ("solution-" +
972 *   Utilities::int_to_string (output_num, step_N) +
973 *   "." +
974 *   Utilities::int_to_string (i,4) +
975 *   ".vtu");
976 *   std::ofstream master_output (("solution-" +
977 *   Utilities::int_to_string (output_num,step_N)+
978 *   ".pvtu").c_str());
979 *   data_out.write_pvtu_record (master_output,filenames);
980 *  
981 *   }
982 *  
983 *   }
984 *  
985 *  
986 *  
987 *  
988 * @endcode
989 *
990 *
991 * <a name="LaserHeatingrun"></a>
992 * <h4>LaserHeating::run</h4>
993 *
994
995 *
996 * This is the function which has the top-level control over everything. Apart
997 * from one line of additional output, it is the same as for the previous
998 * example.
999 *
1000 * @code
1001 *   template <int dim>
1002 *   void LaserHeating<dim>::run ()
1003 *   {
1004 *   pcout << "Solving problem in " << dim << " space dimensions." << std::endl;
1005 *  
1006 *  
1007 *   make_grid();
1008 *   refine_mesh();
1009 *   setup_system ();
1010 *   assemble_system_matrix_init (global_simulation_time_step);
1011 *  
1012 * @endcode
1013 *
1014 * projection of initial conditions by solving.
1015 * solution stored in new_solution_T;
1016 *
1017
1018 *
1019 *
1020 * @code
1021 *   solve_T ();
1022 *  
1023 *   old_solution_T = new_solution_T;
1024 *   old_solution_T_cal = new_solution_T;
1025 *  
1026 * @endcode
1027 *
1028 * reinitialization
1029 *
1030
1031 *
1032 *
1033 * @code
1034 *   system_matrix_T = 0;
1035 *   system_rhs_T = 0;
1036 *  
1037 *  
1038 *   double time_step = global_simulation_time_step;
1039 *   double time = 0;
1040 *   int timestep_number = 0;
1041 *  
1042 * @endcode
1043 *
1044 * output initial values; need ghost cells
1045 *
1046 * @code
1047 *   output_results (0);
1048 *  
1049 *  
1050 *   while(time < global_simulation_end_time)
1051 *   {
1052 *  
1053 *   time += time_step;
1054 *   timestep_number ++;
1055 *  
1056 *   pcout << "Time step " << timestep_number
1057 *   << " at t=" << time
1058 *   << " time_step = " << time_step
1059 *   << std::endl;
1060 *  
1061 * @endcode
1062 *
1063 * the dynamic solving part
1064 *
1065 * @code
1066 *   {
1067 *  
1068 *   right_system_matrix_T.vmult(system_rhs_T,old_solution_T_cal);
1069 *  
1070 *   dynamic_assemble_rhs_T (time,time_step);
1071 *   system_rhs_T.add(1,dynamic_rhs_T);
1072 *  
1073 *   system_matrix_T.copy_from (left_system_matrix_T);
1074 *  
1075 *   {
1076 *   BoundaryValues<dim> boundary_values_function;
1077 *   std::map<types::global_dof_index,double> boundary_values;
1078 *  
1080 *   BOUNDARY_NUM,
1081 *   boundary_values_function,
1082 *   boundary_values);
1083 *  
1084 *   MatrixTools::apply_boundary_values (boundary_values,
1085 *   system_matrix_T,
1086 *   new_solution_T,
1087 *   system_rhs_T,
1088 *   false);
1089 *   }
1090 *  
1091 *  
1092 *   solve_T ();
1093 *  
1094 * @endcode
1095 *
1096 * old_solution_T is used for output, holding ghost cells
1097 * old_solution_T_cal is used for calculation, holding only
1098 * locally owned cells.
1099 *
1100 * @code
1101 *   old_solution_T = new_solution_T;
1102 *   old_solution_T_cal = new_solution_T;
1103 *  
1104 *   if (Utilities::MPI::n_mpi_processes(mpi_communicator) <= 96 && (timestep_number % 50 == 0 ))
1105 *   {
1106 *   TimerOutput::Scope t(computing_timer,"output");
1107 *   output_results (timestep_number);
1108 *   }
1109 *  
1110 *   computing_timer.print_summary ();
1111 *   computing_timer.reset();
1112 *  
1113 *   pcout << std::endl;
1114 *  
1115 *   }
1116 *   }
1117 *  
1118 *   }
1119 *  
1120 *  
1121 * @endcode
1122 *
1123 *
1124 * <a name="Thecodemaincodefunction"></a>
1125 * <h3>The <code>main</code> function</h3>
1126 *
1127
1128 *
1129 *
1130 * @code
1131 *   int main (int argc, char *argv[])
1132 *   {
1133 *   try
1134 *   {
1135 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
1136 *  
1137 *   LaserHeating<2> laserHeating_2d;
1138 *   laserHeating_2d.run ();
1139 *  
1140 *  
1141 *   }
1142 *   catch (std::exception &exc)
1143 *   {
1144 *   std::cerr << std::endl
1145 *   << std::endl
1146 *   << "----------------------------------------------------"
1147 *   << std::endl;
1148 *   std::cerr << "Exception on processing: " << std::endl
1149 *   << exc.what() << std::endl
1150 *   << "Aborting!" << std::endl
1151 *   << "----------------------------------------------------"
1152 *   << std::endl;
1153 *  
1154 *   return 1;
1155 *   }
1156 *   catch (...)
1157 *   {
1158 *   std::cerr << std::endl
1159 *   << std::endl
1160 *   << "----------------------------------------------------"
1161 *   << std::endl;
1162 *   std::cerr << "Unknown exception!" << std::endl
1163 *   << "Aborting!" << std::endl
1164 *   << "----------------------------------------------------"
1165 *   << std::endl;
1166 *   return 1;
1167 *   }
1168 *  
1169 *  
1170 *   return 0;
1171 *   }
1172 * @endcode
1173
1174
1175<a name="ann-boundaryInit.h"></a>
1176<h1>Annotated version of boundaryInit.h</h1>
1177 *
1178 *
1179 *
1180 *
1181
1182 *
1183 * #----------------------------------------------------------
1184 * #
1185 * # This file defines the boundary and initial conditions
1186 * #
1187 * #----------------------------------------------------------
1188 *
1189
1190 *
1191 *
1192 * @code
1193 *   #ifndef GLOBAL_PARA
1194 *   #define GLOBAL_PARA
1195 *   #include "./globalPara.h"
1196 *   #endif
1197 *  
1198 * @endcode
1199 *
1200 * #----------------------------------------------------------
1201 * # Declaration
1202 *
1203
1204 *
1205 *
1206 * @code
1207 *   template <int dim>
1208 *   class BoundaryValues : public Function<dim>
1209 *   {
1210 *   public:
1211 *   BoundaryValues () : Function<dim>() {}
1212 *  
1213 *   virtual double value (const Point<dim> &p,
1214 *   const unsigned int component = 0) const override;
1215 *   };
1216 *  
1217 *   template <int dim>
1218 *   class InitialValues : public Function<dim>
1219 *   {
1220 *   public:
1221 *   InitialValues () : Function<dim>() {}
1222 *  
1223 *   virtual double value (const Point<dim> &p,
1224 *   const unsigned int component = 0) const override;
1225 *   };
1226 *  
1227 *  
1228 * @endcode
1229 *
1230 * #----------------------------------------------------------
1231 * # Implementation
1232 *
1233
1234 *
1235 *
1236 * @code
1237 *   template <int dim>
1238 *   double BoundaryValues<dim>::value (const Point<dim> &/*p*/,
1239 *   const unsigned int /*component*/) const
1240 *   {
1241 *   return 293;
1242 *   }
1243 *  
1244 *  
1245 *   template <int dim>
1246 *   double InitialValues<dim>::value (const Point<dim> &/*p*/,
1247 *   const unsigned int /*component*/) const
1248 *   {
1249 *   return 293;
1250 *   }
1251 *  
1252 *  
1253 *  
1254 *  
1255 *  
1256 * @endcode
1257 *
1258 * #----------------------------------------------------------
1259 * # Declaration and Implementation
1260 * # mass density and heat capacity
1261 *
1262
1263 *
1264 *
1265 * @code
1266 *   template <int dim>
1267 *   class RhoC : public Function<dim>
1268 *   {
1269 *   public:
1270 *   RhoC () : Function<dim>() {}
1271 *  
1272 *   virtual double value (const Point<dim> &p,
1273 *   const unsigned int component = 0) const override;
1274 *   };
1275 *  
1276 *   template <int dim>
1277 *   double RhoC<dim>::value (const Point<dim> &p,
1278 *   const unsigned int /*component*/) const
1279 *   {
1280 * @endcode
1281 *
1282 * # p stores the xyz coordinates at each vertex
1283 * # for 2D problems in xy, we assume the non-uniform is in y-axis.
1284 *
1285 * @code
1286 *   if ( p[1] >= -global_film_thickness )
1287 *   {
1288 *   return global_rho_Tio2 * global_C_Tio2;
1289 *   }
1290 *   else
1291 *   return global_rho_glass * global_C_glass;
1292 *   }
1293 *  
1294 *  
1295 *  
1296 *  
1297 * @endcode
1298 *
1299 * #----------------------------------------------------------
1300 * # Declaration and Implementation
1301 * # thermal conductivity
1302 *
1303
1304 *
1305 *
1306 * @code
1307 *   template <int dim>
1308 *   class K_T : public Function<dim>
1309 *   {
1310 *   public:
1311 *   K_T () : Function<dim>() {}
1312 *  
1313 *   virtual double value (const Point<dim> &p,
1314 *   const unsigned int component = 0) const override;
1315 *   };
1316 *  
1317 *   template <int dim>
1318 *   double K_T<dim>::value (const Point<dim> &p,
1319 *   const unsigned int /*component*/) const
1320 *   {
1321 * @endcode
1322 *
1323 * # p stores the xyz coordinates at each vertex
1324 * # for 2D problems in xy, we assume the non-uniform is in y-axis.
1325 *
1326 * @code
1327 *   if ( p[1] >= -global_film_thickness)
1328 *   {
1329 *   return global_k_Tio2;
1330 *   }
1331 *   else
1332 *   return global_k_glass;
1333 *   }
1334 *  
1335 * @endcode
1336
1337
1338<a name="ann-globalPara.h"></a>
1339<h1>Annotated version of globalPara.h</h1>
1340 *
1341 *
1342 *
1343 *
1344
1345 *
1346 * # physics constants
1347 *
1348 * @code
1349 *   double global_PI = 3.1415927;
1350 *  
1351 * @endcode
1352 *
1353 * # Laser
1354 *
1355 * @code
1356 *   double global_Pow_laser = 0.4; // power [W]
1357 *   double global_spotsize_at_e_2 = 20e-6; // laser spot size at e^(-2) [m]
1358 *   double global_c_laser = global_spotsize_at_e_2 / 4.0; // C parameter in Gaussian func, [m]
1359 *   double global_c_hwhm = global_c_laser * 2.35482 / 2; // HWHM, [m]
1360 *   double global_V_scan_x = 10e-3; // scan speed, [m/s]
1361 *  
1362 *   double global_init_position_x0 = -50e-6; // initial spot center position
1363 *  
1364 * @endcode
1365 *
1366 * # material
1367 * thin film
1368 *
1369 * @code
1370 *   double global_rho_Tio2 = 4200; // mass density, [kg/m^3]
1371 *   double global_C_Tio2 = 690; // heat capacity, [J/kg/K]
1372 *   double global_k_Tio2 = 4.8; // thermal conductivity, [W/m/K]
1373 * @endcode
1374 *
1375 * substrate
1376 *
1377 * @code
1378 *   double global_rho_glass = 2200;
1379 *   double global_C_glass = 700;
1380 *   double global_k_glass = 1.8;
1381 *  
1382 *   double global_film_thickness = 400e-9; // film thickness, [m]
1383 *  
1384 * @endcode
1385 *
1386 * # simulation time
1387 *
1388 * @code
1389 *   double global_simulation_time_step = 1e-5; // 10 [us]
1390 *   double global_simulation_end_time = 100e-6 / global_V_scan_x; // 100 [um] / scan speed
1391 *  
1392 * @endcode
1393 *
1394 * # about the MESH
1395 *
1396 * @code
1397 *   #define BOUNDARY_NUM 11
1398 * @endcode
1399
1400
1401<a name="ann-rightHandSide.h"></a>
1402<h1>Annotated version of rightHandSide.h</h1>
1403 *
1404 *
1405 *
1406 * #----------------------------------------------------------
1407 * #
1408 * # This file defines the boundary and initial conditions
1409 * #
1410 * #----------------------------------------------------------
1411 *
1412
1413 *
1414 *
1415 * @code
1416 *   #ifndef GLOBAL_PARA
1417 *   #define GLOBAL_PARA
1418 *   #include "./globalPara.h"
1419 *   #endif
1420 *  
1421 * @endcode
1422 *
1423 * #----------------------------------------------------------
1424 * # Declaration
1425 *
1426
1427 *
1428 *
1429 * @code
1430 *   template <int dim>
1431 *   class RightHandside : public Function<dim>
1432 *   {
1433 *   public:
1434 *   RightHandside () : Function<dim>() {}
1435 *  
1436 *   double value_v2 (const Point<dim> &p);
1437 *   };
1438 *  
1439 * @endcode
1440 *
1441 * #----------------------------------------------------------
1442 * # Implementation
1443 *
1444
1445 *
1446 *
1447 * @code
1448 *   template <int dim>
1449 *   double RightHandside<dim>::value_v2 (const Point<dim> &p)
1450 *   {
1451 *  
1452 *   double alpha_abs = 1e4;
1453 *  
1454 *  
1455 *   if(p[1] >= -global_film_thickness)
1456 *   {
1457 *  
1458 *   double P00 = global_Pow_laser / global_PI / global_c_laser / global_c_laser / 2.0;
1459 *  
1460 *   double I00 = P00 * std::exp(-
1461 *   (
1462 *   (p[0] - global_V_scan_x * this->get_time()-global_init_position_x0) *
1463 *   (p[0] - global_V_scan_x * this->get_time()-global_init_position_x0)
1464 *   ) /
1465 *   (2.0 * global_c_laser * global_c_laser) );
1466 *  
1467 *  
1468 *   return alpha_abs * I00 * std::exp(-alpha_abs*(0 - p[1]));
1469 *  
1470 *   }
1471 *  
1472 *   else
1473 *   {
1474 * @endcode
1475 *
1476 * # no absorption
1477 *
1478 * @code
1479 *   return 0;
1480 *   }
1481 *  
1482 *   }
1483 *  
1484 *  
1485 * @endcode
1486
1487
1488*/
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
Definition fe_q.h:551
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
void attach_triangulation(Triangulation< dim, spacedim > &tria)
Definition grid_in.cc:156
Definition point.h:112
@ wall_times
Definition timer.h:652
unsigned int level
Definition grid_out.cc:4618
typename ActiveSelector::active_cell_iterator active_cell_iterator
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
const Event initial
Definition event.cc:65
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
@ matrix
Contents is actually a matrix.
@ general
No special properties.
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:150
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:161
std::string get_time()
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:471
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask=ComponentMask())
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void copy(const T *begin, const T *end, U *dest)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
const ::Triangulation< dim, spacedim > & tria