270 * <a name=
"Distributed_Moving_Laser_Heating.cc-Includefiles"></a>
271 * <h3>Include files</h3>
275 * All the include files have already been discussed in previous tutorials.
278 * #include <deal.II/dofs/dof_handler.h>
279 * #include <deal.II/dofs/dof_renumbering.h>
280 * #include <deal.II/grid/grid_generator.h>
281 * #include <deal.II/grid/tria_accessor.h>
282 * #include <deal.II/grid/tria_iterator.h>
283 * #include <deal.II/dofs/dof_accessor.h>
284 * #include <deal.II/fe/fe_q.h>
285 * #include <deal.II/dofs/dof_tools.h>
286 * #include <deal.II/fe/fe_values.h>
287 * #include <deal.II/base/quadrature_lib.h>
288 * #include <deal.II/base/function.h>
289 * #include <deal.II/numerics/vector_tools.h>
290 * #include <deal.II/numerics/matrix_tools.h>
291 * #include <deal.II/lac/vector.h>
292 * #include <deal.II/lac/full_matrix.h>
293 * #include <deal.II/lac/dynamic_sparsity_pattern.h>
295 * #include <deal.II/grid/grid_in.h>
296 * #include <deal.II/grid/grid_tools.h>
297 * #include <deal.II/lac/trilinos_vector.h>
298 * #include <deal.II/lac/trilinos_sparse_matrix.h>
299 * #include <deal.II/lac/trilinos_solver.h>
300 * #include <deal.II/lac/trilinos_precondition.h>
302 * #include <deal.II/lac/sparsity_tools.h>
303 * #include <deal.II/lac/generic_linear_algebra.h>
305 * #include <deal.II/base/conditional_ostream.h>
306 * #include <deal.II/base/utilities.h>
307 * #include <deal.II/base/index_set.h>
308 * #include <deal.II/distributed/tria.h>
310 * #include <deal.II/numerics/data_out.h>
312 * #include <iostream>
315 * #include <deal.II/base/logstream.h>
316 * #include <deal.II/lac/affine_constraints.h>
317 * #include <deal.II/base/timer.h>
324 * #include <deal.II/numerics/error_estimator.h>
325 * #include <deal.II/distributed/grid_refinement.h>
326 * #include <deal.II/distributed/solution_transfer.h>
330 * remember to use the
namespace dealii before
331 * defining a
class that is inheritaged from
Function<dim>
338 * In
general, it is more clear to separate boundary and
initial condictions,
339 * as well as the right hand side function from a file holding all the things.
340 * To
do so, in
this work, the globalPara.h file defines physical constants,
341 * laser parameters, and heat characteristics of materials involved. Boundary
342 * and
initial conditions are defined in boundaryInit.h. The rightHandSide.h
343 * defines the heat source, which in
this work is a moving Gaussian beam.
349 * #ifndef GLOBAL_PARA
350 * #define GLOBAL_PARA
351 * #include
"./globalPara.h"
352 * #include
"./boundaryInit.h"
353 * #include
"./rightHandSide.h"
358 * Now the main
class is defined as following
371 *
void setup_system();
373 *
void assemble_system_matrix_init (
double time_step);
374 *
void dynamic_assemble_rhs_T (
double time,
double time_step);
378 *
void refine_mesh();
379 *
void output_results (
int output_num)
const;
406 *
for storing right
matrix
413 * System_rhs, only locally owned cells
421 * Old Solutions with ghost cells,
for output
428 * Old Solutions only with locally owned cells
435 * New Solutions only with locally owned cells
442 * Dynamic assembling of the righthandside terms
466 * LaserHeating<dim>::LaserHeating ()
468 * mpi_communicator (MPI_COMM_WORLD),
483 * LaserHeating<dim>::~LaserHeating ()
485 * dof_handler.clear();
491 * <a name=
"Distributed_Moving_Laser_Heating.cc-LaserHeatingmake_grid"></a>
492 * <h4>LaserHeating::make_grid</h4>
493 * make grid by importing msh file, and rescale
497 *
void LaserHeating<dim>::make_grid ()
503 * std::ifstream input_file (
"geometry.msh");
504 * grid_in.read_msh (input_file);
507 * pcout <<
" Number of active cells: "
510 * <<
" Total number of cells: "
518 * <a name=
"Distributed_Moving_Laser_Heating.cc-LaserHeatingsetup_system"></a>
519 * <h4>LaserHeating::setup_system</h4>
524 *
void LaserHeating<dim>::setup_system ()
529 * dof_handler.distribute_dofs (fe);
531 * pcout <<
" Number of degrees of freedom: "
532 * << dof_handler.n_dofs()
535 * locally_owned_dofs = dof_handler.locally_owned_dofs();
540 * we want to output solution, so here should have ghost cells
543 * old_solution_T.reinit (locally_owned_dofs,locally_relevant_dofs,mpi_communicator);
547 * locally owned cells
550 * old_solution_T_cal.reinit (locally_owned_dofs,mpi_communicator);
551 * new_solution_T.reinit (locally_owned_dofs,mpi_communicator);
552 * dynamic_rhs_T.reinit (locally_owned_dofs,mpi_communicator);
553 * system_rhs_T.reinit (locally_owned_dofs,mpi_communicator);
555 * constraints_T.clear();
556 * constraints_T.reinit (locally_relevant_dofs);
558 * constraints_T.close();
565 * locally_owned_dofs,
567 * locally_relevant_dofs);
569 * left_system_matrix_T.reinit (locally_owned_dofs,locally_owned_dofs,dsp_T,mpi_communicator);
570 * right_system_matrix_T.reinit (locally_owned_dofs,locally_owned_dofs,dsp_T,mpi_communicator);
571 * system_matrix_T.reinit (locally_owned_dofs,locally_owned_dofs,dsp_T,mpi_communicator);
579 * <a name=
"Distributed_Moving_Laser_Heating.cc-LaserHeatingassemble_system_matrix"></a>
580 * <h4>LaserHeating::assemble_system_matrix</h4>
587 *
void LaserHeating<dim>::assemble_system_matrix_init (
double time_step)
595 *
const InitialValues<dim> initial_value_func_T;
598 *
const RhoC<dim> rho_C_fun_T;
599 *
const K_T<dim> k_fun_T;
606 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
607 *
const unsigned int n_q_points = quadrature_formula.size();
617 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
619 * system_matrix_T = 0;
620 * left_system_matrix_T = 0;
621 * right_system_matrix_T = 0;
625 *
for (
const auto &cell : dof_handler.active_cell_iterators())
626 * if(cell->is_locally_owned())
629 * fe_values.
reinit (cell);
630 * local_init_matrix = 0;
632 * local_rho_c_T_matrix = 0;
633 * local_k_T_matrix = 0;
635 * local_init_T_rhs = 0;
637 *
for (
unsigned int q=0; q<n_q_points; ++q)
638 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
640 *
const Tensor<1,dim> div_phi_i_u = fe_values.shape_grad (i,q);
641 *
const double phi_i_u = fe_values.shape_value (i,q);
643 *
for (
unsigned int j=0; j<dofs_per_cell; ++j)
646 *
const Tensor<1,dim> div_phi_j_u = fe_values.shape_grad(j,q);
647 *
const double phi_j_u = fe_values.shape_value (j,q);
649 * local_init_matrix(i,j) += (phi_i_u *
651 * fe_values.JxW (q));
653 * local_rho_c_T_matrix(i,j) += (rho_C_fun_T.value(fe_values.quadrature_point(q)) *
658 * time_step * (theta) *
659 * (k_fun_T.value(fe_values.quadrature_point(q)) *
662 * fe_values.JxW (q));
664 * local_k_T_matrix(i,j) += (rho_C_fun_T.value(fe_values.quadrature_point(q)) *
669 * time_step * (1.0-theta) *
670 * (k_fun_T.value(fe_values.quadrature_point(q)) *
673 * fe_values.JxW (q));
677 * local_init_T_rhs(i) += (phi_i_u *
678 * initial_value_func_T.value (fe_values.quadrature_point (q)) *
679 * fe_values.JxW (q));
683 * cell->get_dof_indices (local_dof_indices);
690 * constraints_T.distribute_local_to_global(local_init_matrix,
699 * store M + dt*theta*
A as the left_system_matrix
705 * constraints_T.distribute_local_to_global(local_rho_c_T_matrix,
707 * left_system_matrix_T);
711 * store M - dt*(1-theta)*A as the right_system_matrix
714 * constraints_T.distribute_local_to_global(local_k_T_matrix,
716 * right_system_matrix_T);
732 * <a name=
"Distributed_Moving_Laser_Heating.cc-LaserHeatingdynamic_assemble_rhs_T"></a>
733 * <h4>LaserHeating::dynamic_assemble_rhs_T</h4>
734 * The right hand side is assembled each time during running, which is necessary as
735 * the laser source is moving. To separate the heat source and the right hand
736 * side assembling, the right hand side function is defined as RightHandside<dim>.
740 *
void LaserHeating<dim>::dynamic_assemble_rhs_T (
double time,
double time_step)
747 * RightHandside<dim> rhs_func_T_1;
748 * rhs_func_T_1.set_time(time);
750 * RightHandside<dim> rhs_func_T_2;
751 * rhs_func_T_2.set_time(time-time_step);
758 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
759 *
const unsigned int n_q_points = quadrature_formula.size();
763 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
765 * std::vector<double> Rnp_cal_assemble (n_q_points);
768 * dynamic_rhs_T = 0 ;
771 * cell = dof_handler.begin_active(),
772 * endc = dof_handler.end();
774 *
for (; cell!=endc; ++cell)
775 *
if(cell->is_locally_owned())
777 * fe_values.reinit (cell);
778 * local_rhs_vector_T = 0;
780 *
for (
unsigned int q=0; q<n_q_points; ++q)
781 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
783 *
const double phi_i_u = fe_values.shape_value (i,q);
786 * local_rhs_vector_T(i) += time_step * theta *
788 * rhs_func_T_1.value_v2 (fe_values.quadrature_point (q)) *
791 * time_step * (1.0 - theta) *
793 * rhs_func_T_2.value_v2 (fe_values.quadrature_point (q)) *
794 * fe_values.JxW (q));
797 * cell->get_dof_indices (local_dof_indices);
800 * constraints_T.distribute_local_to_global(local_rhs_vector_T,
815 * <a name=
"Distributed_Moving_Laser_Heating.cc-LaserHeatingsolve"></a>
816 * <h4>LaserHeating::solve</h4>
817 * Solving the equation is direct. Recall that we have defined several matrices and vectors,
818 * to avoid ambiguous, here, we only use system_matrix_T as the system
matrix, system_rhs_T
819 * as the system right hand side. The vector completely_distributed_solution is used to store
820 * the obtained solution.
824 *
void LaserHeating<dim>::solve_T ()
829 *
SolverControl solver_control (1*system_rhs_T.size(),1e-12*system_rhs_T.l2_norm(),
true);
836 * preconditioner.initialize(system_matrix_T,
data);
838 * solver.solve (system_matrix_T,completely_distributed_solution,system_rhs_T,preconditioner);
843 * Print the number of iterations by hand.
849 * pcout <<
" " << solver_control.last_step()
850 * <<
" CG iterations needed to obtain convergence." << std::endl
851 * <<
"\t initial convergence value = " << solver_control.initial_value() << std::endl
852 * <<
"\t final convergence value = " << solver_control.last_value() << std::endl
855 * constraints_T.distribute (completely_distributed_solution);
856 * new_solution_T = completely_distributed_solution;
864 * <a name=
"Distributed_Moving_Laser_Heating.cc-LaserHeatingrefine_mesh"></a>
865 * <h4>LaserHeating::refine_mesh</h4>
872 *
void LaserHeating<dim>::refine_mesh()
885 * only
refine mesh 5um above the TiO2 and glass
interface
894 *
if(cell->is_locally_owned())
896 * fe_values.reinit(cell);
897 *
if(
std::abs(fe_values.quadrature_point(0)[1]) <= global_film_thickness+5
e-6)
899 * cell->set_refine_flag();
903 * cell->clear_refine_flag();
914 * <a name=
"Distributed_Moving_Laser_Heating.cc-LaserHeatinghoutput_results"></a>
915 * <h4>LaserHeatingh::output_results</h4>
922 *
void LaserHeating<dim>::output_results (
int output_num)
const
928 * data_out.add_data_vector (old_solution_T,
"T");
932 * the length of output numbering
940 * data_out.build_patches ();
942 *
const std::string filename = (
"solution-" +
947 * std::ofstream output (filename.c_str());
948 * data_out.write_vtu (output);
953 * output the overall solution
961 * std::vector<std::string> filenames;
963 * filenames.push_back (
"solution-" +
968 * std::ofstream master_output ((
"solution-" +
971 * data_out.write_pvtu_record (master_output,filenames);
983 * <a name=
"Distributed_Moving_Laser_Heating.cc-LaserHeatingrun"></a>
984 * <h4>LaserHeating::run</h4>
988 * This is the function which has the top-
level control over everything. Apart
989 * from
one line of additional output, it is the same as
for the previous
994 *
void LaserHeating<dim>::run ()
996 * pcout <<
"Solving problem in " << dim <<
" space dimensions." << std::endl;
1002 * assemble_system_matrix_init (global_simulation_time_step);
1006 * projection of
initial conditions by solving.
1007 * solution stored in new_solution_T;
1015 * old_solution_T = new_solution_T;
1016 * old_solution_T_cal = new_solution_T;
1026 * system_matrix_T = 0;
1030 *
double time_step = global_simulation_time_step;
1032 *
int timestep_number = 0;
1039 * output_results (0);
1042 *
while(time < global_simulation_end_time)
1045 * time += time_step;
1046 * timestep_number ++;
1048 * pcout <<
"Time step " << timestep_number
1049 * <<
" at t=" << time
1050 * <<
" time_step = " << time_step
1060 * right_system_matrix_T.vmult(system_rhs_T,old_solution_T_cal);
1062 * dynamic_assemble_rhs_T (time,time_step);
1063 * system_rhs_T.add(1,dynamic_rhs_T);
1065 * system_matrix_T.copy_from (left_system_matrix_T);
1068 * BoundaryValues<dim> boundary_values_function;
1069 * std::map<types::global_dof_index,double> boundary_values;
1073 * boundary_values_function,
1088 * old_solution_T is used
for output, holding ghost cells
1089 * old_solution_T_cal is used
for calculation, holding only
1090 * locally owned cells.
1093 * old_solution_T = new_solution_T;
1094 * old_solution_T_cal = new_solution_T;
1099 * output_results (timestep_number);
1102 * computing_timer.print_summary ();
1103 * computing_timer.reset();
1105 * pcout << std::endl;
1116 * <a name=
"Distributed_Moving_Laser_Heating.cc-Thecodemaincodefunction"></a>
1117 * <h3>The <code>main</code> function</h3>
1123 *
int main (
int argc,
char *argv[])
1129 * LaserHeating<2> laserHeating_2d;
1130 * laserHeating_2d.run ();
1134 *
catch (std::exception &exc)
1136 * std::cerr << std::endl
1138 * <<
"----------------------------------------------------"
1140 * std::cerr <<
"Exception on processing: " << std::endl
1141 * << exc.what() << std::endl
1142 * <<
"Aborting!" << std::endl
1143 * <<
"----------------------------------------------------"
1150 * std::cerr << std::endl
1152 * <<
"----------------------------------------------------"
1154 * std::cerr <<
"Unknown exception!" << std::endl
1155 * <<
"Aborting!" << std::endl
1156 * <<
"----------------------------------------------------"
1167<a name=
"ann-boundaryInit.h"></a>
1168<h1>Annotated version of boundaryInit.h</h1>
1187 * #----------------------------------------------------------
1189 * # This file defines the boundary and
initial conditions
1191 * #----------------------------------------------------------
1197 * #ifndef GLOBAL_PARA
1198 * #define GLOBAL_PARA
1199 * #include
"./globalPara.h"
1204 * #----------------------------------------------------------
1211 *
template <
int dim>
1212 *
class BoundaryValues :
public Function<dim>
1215 * BoundaryValues () :
Function<dim>() {}
1218 *
const unsigned int component = 0)
const override;
1221 *
template <
int dim>
1222 *
class InitialValues :
public Function<dim>
1225 * InitialValues () :
Function<dim>() {}
1228 *
const unsigned int component = 0)
const override;
1234 * #----------------------------------------------------------
1241 *
template <
int dim>
1242 *
double BoundaryValues<dim>::value (
const Point<dim> &,
1243 *
const unsigned int )
const
1249 *
template <
int dim>
1250 *
double InitialValues<dim>::value (
const Point<dim> &,
1251 *
const unsigned int )
const
1262 * #----------------------------------------------------------
1263 * # Declaration and Implementation
1264 * # mass density and heat capacity
1270 *
template <
int dim>
1271 *
class RhoC :
public Function<dim>
1277 *
const unsigned int component = 0)
const override;
1280 *
template <
int dim>
1281 *
double RhoC<dim>::value (
const Point<dim> &p,
1282 *
const unsigned int )
const
1286 * # p stores the xyz coordinates at each vertex
1287 * #
for 2D problems in xy, we assume the non-uniform is in y-axis.
1290 *
if ( p[1] >= -global_film_thickness )
1292 *
return global_rho_Tio2 * global_C_Tio2;
1295 *
return global_rho_glass * global_C_glass;
1303 * #----------------------------------------------------------
1304 * # Declaration and Implementation
1305 * # thermal conductivity
1311 *
template <
int dim>
1318 *
const unsigned int component = 0)
const override;
1321 *
template <
int dim>
1322 *
double K_T<dim>::value (
const Point<dim> &p,
1323 *
const unsigned int )
const
1327 * # p stores the xyz coordinates at each vertex
1328 * #
for 2D problems in xy, we assume the non-uniform is in y-axis.
1331 *
if ( p[1] >= -global_film_thickness)
1333 *
return global_k_Tio2;
1336 *
return global_k_glass;
1342<a name=
"ann-globalPara.h"></a>
1343<h1>Annotated version of globalPara.h</h1>
1362 * # physics constants
1365 *
double global_PI = 3.1415927;
1372 *
double global_Pow_laser = 0.4;
1373 *
double global_spotsize_at_e_2 = 20
e-6;
1374 *
double global_c_laser = global_spotsize_at_e_2 / 4.0;
1375 *
double global_c_hwhm = global_c_laser * 2.35482 / 2;
1376 *
double global_V_scan_x = 10
e-3;
1378 *
double global_init_position_x0 = -50
e-6;
1386 *
double global_rho_Tio2 = 4200;
1387 *
double global_C_Tio2 = 690;
1388 *
double global_k_Tio2 = 4.8;
1394 *
double global_rho_glass = 2200;
1395 *
double global_C_glass = 700;
1396 *
double global_k_glass = 1.8;
1398 *
double global_film_thickness = 400
e-9;
1405 *
double global_simulation_time_step = 1
e-5;
1406 *
double global_simulation_end_time = 100
e-6 / global_V_scan_x;
1413 * #define BOUNDARY_NUM 11
1417<a name=
"ann-rightHandSide.h"></a>
1418<h1>Annotated version of rightHandSide.h</h1>
1437 * #----------------------------------------------------------
1439 * # This file defines the boundary and
initial conditions
1441 * #----------------------------------------------------------
1447 * #ifndef GLOBAL_PARA
1448 * #define GLOBAL_PARA
1449 * #include
"./globalPara.h"
1454 * #----------------------------------------------------------
1461 *
template <
int dim>
1462 *
class RightHandside :
public Function<dim>
1465 * RightHandside () :
Function<dim>() {}
1472 * #----------------------------------------------------------
1479 *
template <
int dim>
1480 *
double RightHandside<dim>::value_v2 (
const Point<dim> &p)
1483 *
double alpha_abs = 1e4;
1486 *
if(p[1] >= -global_film_thickness)
1489 *
double P00 = global_Pow_laser / global_PI / global_c_laser / global_c_laser / 2.0;
1493 * (p[0] - global_V_scan_x * this->
get_time()-global_init_position_x0) *
1494 * (p[0] - global_V_scan_x * this->
get_time()-global_init_position_x0)
1496 * (2.0 * global_c_laser * global_c_laser) );
1499 *
return alpha_abs * I00 *
std::exp(-alpha_abs*(0 - p[1]));
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
void attach_triangulation(Triangulation< dim, spacedim > &tria)
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.
std::vector< index_type > data
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
@ matrix
Contents is actually a matrix.
@ general
No special properties.
constexpr types::blas_int one
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
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