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