From a25a9fce14b5b7384f9e077aea26ce771f018a3a Mon Sep 17 00:00:00 2001 From: Krishna Kumar Date: Wed, 1 Apr 2020 15:36:26 -0500 Subject: [PATCH] :art: Format THM code --- clang-tools/format.sh | 1 - include/boundaryvalues.h | 24 ++-- include/geothermal.h | 247 ++++++++++++++++---------------------- include/globalvariables.h | 5 +- include/initialvalues.h | 58 ++++----- include/sourceterm.h | 86 +++++++------ src/main.cc | 18 +-- 7 files changed, 189 insertions(+), 250 deletions(-) diff --git a/clang-tools/format.sh b/clang-tools/format.sh index 8826a05..a5324db 100644 --- a/clang-tools/format.sh +++ b/clang-tools/format.sh @@ -1,4 +1,3 @@ #!/bin/bash find ../include/ -iname *.h -o -iname *.tcc | xargs clang-format -i find ../src/ -iname *.h -o -iname *.cc | xargs clang-format -i -find ../tests/ -iname *.h -o -iname *.cc | xargs clang-format -i diff --git a/include/boundaryvalues.h b/include/boundaryvalues.h index aaf8cba..623ee61 100644 --- a/include/boundaryvalues.h +++ b/include/boundaryvalues.h @@ -13,10 +13,8 @@ namespace EquationData { template class PressureDirichletBoundaryValues : public Function { public: - PressureDirichletBoundaryValues() - : Function() - , period(0.2) - {} // 之前的没有 + PressureDirichletBoundaryValues() + : Function(), period(0.2) {} // 之前的没有 virtual double value(const Point& p, const unsigned int component = 0) const; // boundary @@ -41,10 +39,8 @@ double PressureDirichletBoundaryValues::value( template class TemperatureDirichletBoundaryValues : public Function { public: - TemperatureDirichletBoundaryValues() - : Function() - , period(0.2) - {} // 之前的没有 + TemperatureDirichletBoundaryValues() + : Function(), period(0.2) {} // 之前的没有 virtual double value(const Point& p, const unsigned int component = 0) const; // boundary @@ -79,9 +75,7 @@ template class PressureNeumanBoundaryValues : public Function { public: PressureNeumanBoundaryValues() - : Function() - , period(0.2) - {} // 之前的没有 + : Function(), period(0.2) {} // 之前的没有 virtual double value(const Point& p, const unsigned int component = 0) const; // boundary @@ -106,10 +100,8 @@ double PressureNeumanBoundaryValues::value( template class TemperatureNeumanBoundaryValues : public Function { public: - TemperatureNeumanBoundaryValues() - : Function() - , period(0.2) - {} // 之前的没有 + TemperatureNeumanBoundaryValues() + : Function(), period(0.2) {} // 之前的没有 virtual double value(const Point& p, const unsigned int component = 0) const; // boundary // virtual void vector_value(const Point& p, //放在这里没啥用 @@ -129,7 +121,7 @@ double TemperatureNeumanBoundaryValues::value( const double time = this->get_time(); return 0.; // boundary value is set to zero in - // this case + // this case } // template diff --git a/include/geothermal.h b/include/geothermal.h index a269e49..526685b 100644 --- a/include/geothermal.h +++ b/include/geothermal.h @@ -46,13 +46,12 @@ using namespace dealii; template -class CoupledTH -{ -public: +class CoupledTH { + public: CoupledTH(const unsigned int degree); void run(); -private: + private: void make_grid_and_dofs(); void assemble_T_system(); void assemble_P_system(); @@ -60,41 +59,40 @@ class CoupledTH void linear_solve_T(); void output_results() const; - Triangulation triangulation; // grid - const unsigned int T_degree; //element degree + Triangulation triangulation; // grid + const unsigned int T_degree; // element degree const unsigned int P_degree; - FE_Q P_fe; // element type - FE_Q T_fe; // element type - DoFHandler dof_handler; // grid<->eleemnt - + FE_Q P_fe; // element type + FE_Q T_fe; // element type + DoFHandler dof_handler; // grid<->eleemnt // ConstraintMatrix constraints; // hanging node - SparsityPattern sparsity_pattern; // sparsity - SparseMatrix P_mass_matrix; // M_P - SparseMatrix T_mass_matrix; // M_T - SparseMatrix P_stiffness_matrix; // K_P - SparseMatrix T_stiffness_matrix; // K_T - SparseMatrix P_system_matrix; // M_P + k*theta*K_P - SparseMatrix T_system_matrix; // M_T + k*theta*K_T - - Vector P_solution; // P solution at n - Vector T_solution; // T solution at n - Vector old_P_solution; // P solution at n-1 - Vector old_T_solution; // T solution at n-1 - Vector P_system_rhs; // right hand side of P system - Vector T_system_rhs; // right hand side of T system - - double time; //t - double time_step; //dt - unsigned int timestep_number; //n_t + SparsityPattern sparsity_pattern; // sparsity + SparseMatrix P_mass_matrix; // M_P + SparseMatrix T_mass_matrix; // M_T + SparseMatrix P_stiffness_matrix; // K_P + SparseMatrix T_stiffness_matrix; // K_T + SparseMatrix P_system_matrix; // M_P + k*theta*K_P + SparseMatrix T_system_matrix; // M_T + k*theta*K_T + + Vector P_solution; // P solution at n + Vector T_solution; // T solution at n + Vector old_P_solution; // P solution at n-1 + Vector old_T_solution; // T solution at n-1 + Vector P_system_rhs; // right hand side of P system + Vector T_system_rhs; // right hand side of T system + + double time; // t + double time_step; // dt + unsigned int timestep_number; // n_t const double theta; }; template -CoupledTH::CoupledTH(const unsigned int degree) // initialization +CoupledTH::CoupledTH(const unsigned int degree) // initialization : T_degree(degree), P_degree(degree), P_fe(P_degree), @@ -102,20 +100,17 @@ CoupledTH::CoupledTH(const unsigned int degree) // initialization dof_handler(triangulation), time(0.0), - time_step(0.5 / 20), // a time step constant at 1/500 (remember that one - // period of the source on the right hand side was - // set to 0.2 above, so we resolve each period with - // 100 time steps) + time_step(0.5 / 20), // a time step constant at 1/500 (remember that one + // period of the source on the right hand side was + // set to 0.2 above, so we resolve each period with + // 100 time steps) timestep_number(0), - theta(0.5) -{ -} + theta(0.5) {} // CHECK THE GRID INPUT template -void print_mesh_info(const Triangulation &triangulation, - const std::string &filename) -{ +void print_mesh_info(const Triangulation& triangulation, + const std::string& filename) { std::cout << "Mesh info:" << std::endl << " dimension: " << dim << std::endl << " no. of cells: " << triangulation.n_active_cells() << std::endl; @@ -124,11 +119,9 @@ void print_mesh_info(const Triangulation &triangulation, typename Triangulation::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end(); - for (; cell != endc; ++cell) - { + for (; cell != endc; ++cell) { for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; - ++face) - { + ++face) { if (cell->face(face)->at_boundary()) boundary_count[cell->face(face)->boundary_id()]++; } @@ -137,8 +130,7 @@ void print_mesh_info(const Triangulation &triangulation, std::cout << " boundary indicators: "; for (std::map::iterator it = boundary_count.begin(); - it != boundary_count.end(); ++it) - { + it != boundary_count.end(); ++it) { std::cout << it->first << "(" << it->second << " times) "; } std::cout << std::endl; @@ -147,33 +139,30 @@ void print_mesh_info(const Triangulation &triangulation, std::ofstream out(filename.c_str()); GridOut grid_out; grid_out.write_eps(triangulation, out); - std::cout << " written to " << filename << std::endl - << std::endl; + std::cout << " written to " << filename << std::endl << std::endl; } template -void CoupledTH::make_grid_and_dofs() -{ - GridIn gridin; //instantiate a gridinput +void CoupledTH::make_grid_and_dofs() { + GridIn gridin; // instantiate a gridinput gridin.attach_triangulation(triangulation); std::ifstream f("mesh.msh"); gridin.read_msh(f); print_mesh_info(triangulation, "grid-1.eps"); - dof_handler.distribute_dofs(P_fe); // distribute dofs to grid - dof_handler.distribute_dofs(T_fe); // distribute dofs to grid + dof_handler.distribute_dofs(P_fe); // distribute dofs to grid + dof_handler.distribute_dofs(T_fe); // distribute dofs to grid std::cout << "Number of active cells: " << triangulation.n_active_cells() << " (on " << triangulation.n_levels() << " levels)" << std::endl - << "Number of degrees of freedom: " - << 2*dof_handler.n_dofs() << " (" - << dof_handler.n_dofs() - << '+' << dof_handler.n_dofs() << ')' << std::endl + << "Number of degrees of freedom: " << 2 * dof_handler.n_dofs() + << " (" << dof_handler.n_dofs() << '+' << dof_handler.n_dofs() + << ')' << std::endl << std::endl; // sparsity pattern for T - DynamicSparsityPattern T_dsp(dof_handler.n_dofs()); // sparsity + DynamicSparsityPattern T_dsp(dof_handler.n_dofs()); // sparsity DoFTools::make_sparsity_pattern( dof_handler, T_dsp /*constraints, keep_constrained_dofs = ,true*/); sparsity_pattern.copy_from(T_dsp); @@ -187,7 +176,7 @@ void CoupledTH::make_grid_and_dofs() T_system_rhs.reinit(dof_handler.n_dofs()); // sparsity pattern for P - DynamicSparsityPattern P_dsp(dof_handler.n_dofs()); // sparsity + DynamicSparsityPattern P_dsp(dof_handler.n_dofs()); // sparsity DoFTools::make_sparsity_pattern( dof_handler, P_dsp /*constraints, keep_constrained_dofs = ,true*/); sparsity_pattern.copy_from(P_dsp); @@ -208,8 +197,7 @@ void CoupledTH::make_grid_and_dofs() } template -void CoupledTH::assemble_P_system() -{ +void CoupledTH::assemble_P_system() { // reset matreix to zero // p_mass_matrix = 0; @@ -220,9 +208,9 @@ void CoupledTH::assemble_P_system() QGauss T_face_quadrature_formula(T_degree + 1); - FEValues T_fe_values( T_fe, T_quadrature_formula, - update_values | update_gradients | - update_gradients | update_JxW_values); + FEValues T_fe_values( + T_fe, T_quadrature_formula, + update_values | update_gradients | update_gradients | update_JxW_values); FEFaceValues T_fe_face_values(T_fe, T_face_quadrature_formula, update_values | update_normal_vectors | @@ -236,8 +224,7 @@ void CoupledTH::assemble_P_system() FEValues P_fe_values( P_fe, P_quadrature_formula, - update_values | update_gradients | - update_gradients | update_JxW_values); + update_values | update_gradients | update_gradients | update_JxW_values); FEFaceValues P_fe_face_values(P_fe, P_face_quadrature_formula, update_values | update_normal_vectors | @@ -278,11 +265,10 @@ void CoupledTH::assemble_P_system() EquationData::PressureDirichletBoundaryValues P_boundary; // loop for cell - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell != endc; ++cell) - { + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell != endc; ++cell) { // initialization P_local_mass_matrix = 0; P_local_stiffness_matrix = 0; @@ -290,7 +276,8 @@ void CoupledTH::assemble_P_system() T_fe_values.reinit(cell); P_fe_values.reinit(cell); - // get the values at gauss point old_T_sol_values from the system old_T_solution + // get the values at gauss point old_T_sol_values from the system + // old_T_solution T_fe_values.get_function_values(old_T_solution, old_T_sol_values); T_fe_values.get_function_gradients(old_T_solution, old_T_sol_grads); P_fe_values.get_function_values(old_P_solution, old_P_sol_values); @@ -298,17 +285,15 @@ void CoupledTH::assemble_P_system() // get source term value at the gauss point P_source_term.set_time(time); - P_source_term.value_list(P_fe_values.get_quadrature_points(), P_source_values); // 一列q个 + P_source_term.value_list(P_fe_values.get_quadrature_points(), + P_source_values); // 一列q个 // loop for q_point ASSMBLING CELL METRIX - for (unsigned int q = 0; q < P_n_q_points; ++q) - { - for (unsigned int i = 0; i < P_dofs_per_cell; ++i) - { + for (unsigned int q = 0; q < P_n_q_points; ++q) { + for (unsigned int i = 0; i < P_dofs_per_cell; ++i) { const Tensor<1, dim> grad_phi_i_P = P_fe_values.shape_grad(i, q); const double phi_i_P = P_fe_values.shape_value(i, q); - for (unsigned int j = 0; j < P_dofs_per_cell; ++j) - { + for (unsigned int j = 0; j < P_dofs_per_cell; ++j) { const Tensor<1, dim> grad_phi_j_P = P_fe_values.shape_grad(j, q); const double phi_j_P = P_fe_values.shape_value(j, q); P_local_mass_matrix(i, j) += (phi_i_P * phi_j_P * P_fe_values.JxW(q)); @@ -325,22 +310,19 @@ void CoupledTH::assemble_P_system() // APPLIED BOUNDARY CONDITION for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; - ++face_no) - { - if (cell->at_boundary(face_no)) - { + ++face_no) { + if (cell->at_boundary(face_no)) { P_fe_face_values.reinit(cell, face_no); // get boundary condition QP_boundary.set_time(time); - QP_boundary.value_list(P_fe_face_values.get_quadrature_points(), QP_bd_values); - - for (unsigned int q = 0; q < P_n_face_q_points; ++q) - { - for (unsigned int i = 0; i < P_dofs_per_cell; ++i) - { - P_local_rhs(i) += -(P_fe_face_values.value(i, q) * - QP_bd_values[q] * P_fe_face_values.JxW(q)); + QP_boundary.value_list(P_fe_face_values.get_quadrature_points(), + QP_bd_values); + + for (unsigned int q = 0; q < P_n_face_q_points; ++q) { + for (unsigned int i = 0; i < P_dofs_per_cell; ++i) { + P_local_rhs(i) += -(P_fe_face_values.value(i, q) * QP_bd_values[q] * + P_fe_face_values.JxW(q)); } } } @@ -349,10 +331,8 @@ void CoupledTH::assemble_P_system() // local ->globe cell->get_dof_indices(P_local_dof_indices); - for (unsigned int i = 0; i < P_dofs_per_cell; ++i) - { - for (unsigned int j = 0; j < P_dofs_per_cell; ++j) - { + for (unsigned int i = 0; i < P_dofs_per_cell; ++i) { + for (unsigned int j = 0; j < P_dofs_per_cell; ++j) { P_mass_matrix.add(P_local_dof_indices[i], P_local_dof_indices[j], P_local_mass_matrix(i, j)); P_stiffness_matrix.add(P_local_dof_indices[i], P_local_dof_indices[j], @@ -360,8 +340,8 @@ void CoupledTH::assemble_P_system() P_system_matrix.copy_from(P_mass_matrix); P_system_matrix.add( theta * time_step, - P_stiffness_matrix); // P_mass_matrix + - // theta*time_step*P_stiffness_matrix + P_stiffness_matrix); // P_mass_matrix + + // theta*time_step*P_stiffness_matrix } P_system_rhs(P_local_dof_indices[i]) += P_local_rhs(i); } @@ -379,8 +359,7 @@ void CoupledTH::assemble_P_system() } template -void CoupledTH::assemble_T_system() -{ +void CoupledTH::assemble_T_system() { // reset matreix to zero // T_mass_matrix = 0; @@ -446,17 +425,16 @@ void CoupledTH::assemble_T_system() EquationData::TemperatureDirichletBoundaryValues T_boundary; // loop for cell - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell != endc; ++cell) - { + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell != endc; ++cell) { // initialization T_local_mass_matrix = 0; T_local_stiffness_matrix = 0; T_local_rhs = 0; T_fe_values.reinit(cell); - P_fe_values.reinit(cell); //// may combine in the same cell + P_fe_values.reinit(cell); //// may combine in the same cell // get teh values at gauss point T_fe_values.get_function_values(old_T_solution, old_T_sol_values); T_fe_values.get_function_gradients(old_T_solution, old_T_sol_grads); @@ -464,17 +442,15 @@ void CoupledTH::assemble_T_system() P_fe_values.get_function_gradients(old_P_solution, old_P_sol_grads); // get right hand side T_source_term.set_time(time); - T_source_term.value_list(T_fe_values.get_quadrature_points(), T_source_values); + T_source_term.value_list(T_fe_values.get_quadrature_points(), + T_source_values); // loop for q_point ASSMBLING CELL METRIX - for (unsigned int q = 0; q < T_n_q_points; ++q) - { - for (unsigned int i = 0; i < T_dofs_per_cell; ++i) - { + for (unsigned int q = 0; q < T_n_q_points; ++q) { + for (unsigned int i = 0; i < T_dofs_per_cell; ++i) { const Tensor<1, dim> grad_phi_i_T = T_fe_values.shape_grad(i, q); const double phi_i_T = T_fe_values.shape_value(i, q); - for (unsigned int j = 0; j < T_dofs_per_cell; ++j) - { + for (unsigned int j = 0; j < T_dofs_per_cell; ++j) { const Tensor<1, dim> grad_phi_j_T = T_fe_values.shape_grad(j, q); const double phi_j_T = T_fe_values.shape_value(j, q); T_local_mass_matrix(i, j) += (phi_i_T * phi_j_T * T_fe_values.JxW(q)); @@ -491,19 +467,15 @@ void CoupledTH::assemble_T_system() // APPLIED NEUMAN BOUNDARY CONDITION for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; - ++face_no) - { - if (cell->at_boundary(face_no)) - { + ++face_no) { + if (cell->at_boundary(face_no)) { T_fe_face_values.reinit(cell, face_no); // set boudnary condition QT_boundary.set_time(time); QT_boundary.value_list(T_fe_face_values.get_quadrature_points(), QT_bd_values); - for (unsigned int q = 0; q < T_n_face_q_points; ++q) - { - for (unsigned int i = 0; i < T_dofs_per_cell; ++i) - { + for (unsigned int q = 0; q < T_n_face_q_points; ++q) { + for (unsigned int i = 0; i < T_dofs_per_cell; ++i) { T_local_rhs(i) += -(T_fe_face_values.shape_value(i, q) * QT_bd_values[q] * T_fe_face_values.JxW(q)); } @@ -512,10 +484,8 @@ void CoupledTH::assemble_T_system() } // local ->globe cell->get_dof_indices(T_local_dof_indices); - for (unsigned int i = 0; i < T_dofs_per_cell; ++i) - { - for (unsigned int j = 0; j < T_dofs_per_cell; ++j) - { + for (unsigned int i = 0; i < T_dofs_per_cell; ++i) { + for (unsigned int j = 0; j < T_dofs_per_cell; ++j) { T_mass_matrix.add(T_local_dof_indices[i], T_local_dof_indices[j], T_local_mass_matrix(i, j)); T_stiffness_matrix.add(T_local_dof_indices[i], T_local_dof_indices[j], @@ -523,8 +493,8 @@ void CoupledTH::assemble_T_system() T_system_matrix.copy_from(T_mass_matrix); T_system_matrix.add( theta * time_step, - T_stiffness_matrix); // T_mass_matrix + - // theta*time_step*T_stiffness_matrix + T_stiffness_matrix); // T_mass_matrix + + // theta*time_step*T_stiffness_matrix } T_system_rhs(T_local_dof_indices[i]) += T_local_rhs(i); } @@ -541,16 +511,15 @@ void CoupledTH::assemble_T_system() } template -void CoupledTH::linear_solve_T() -{ +void CoupledTH::linear_solve_T() { SolverControl solver_control( 1000, - 1e-8 * T_system_rhs.l2_norm()); // setting for cg - SolverCG<> cg(solver_control); // config cg - PreconditionSSOR<> preconditioner; // precond - preconditioner.initialize(T_system_matrix, 1.0); // initialize precond + 1e-8 * T_system_rhs.l2_norm()); // setting for cg + SolverCG<> cg(solver_control); // config cg + PreconditionSSOR<> preconditioner; // precond + preconditioner.initialize(T_system_matrix, 1.0); // initialize precond cg.solve(T_system_matrix, T_solution, T_system_rhs, - preconditioner); // solve eq + preconditioner); // solve eq // constraints.distribute(solution); // make sure if the value is // consistent at // the constraint point @@ -563,8 +532,7 @@ void CoupledTH::linear_solve_T() // // Neither is there anything new in generating graphical output: template -void CoupledTH::output_results() const -{ +void CoupledTH::output_results() const { DataOut T_data_out; T_data_out.attach_dof_handler(dof_handler); @@ -579,20 +547,17 @@ void CoupledTH::output_results() const } template -void CoupledTH::run() -{ +void CoupledTH::run() { make_grid_and_dofs(); // assemble_T_system(); VectorTools::interpolate(dof_handler, EquationData::TemperatureInitialValues(), old_T_solution); - VectorTools::interpolate(dof_handler, - EquationData::PressureInitialValues(), - old_P_solution); + VectorTools::interpolate( + dof_handler, EquationData::PressureInitialValues(), old_P_solution); - do - { + do { std::cout << "Timestep " << timestep_number << std::endl; assemble_T_system(); diff --git a/include/globalvariables.h b/include/globalvariables.h index 7a68b02..5d6d591 100644 --- a/include/globalvariables.h +++ b/include/globalvariables.h @@ -1,7 +1,6 @@ #pragma once -namespace EquationData -{ +namespace EquationData { const double perm = 1e-9; const double eta = 1; const double kappa = 1000; @@ -10,4 +9,4 @@ const double density = 1; const double T0 = 0; const double P0 = 100000.; -} \ No newline at end of file +} // namespace EquationData \ No newline at end of file diff --git a/include/initialvalues.h b/include/initialvalues.h index abad470..a8d9c7a 100644 --- a/include/initialvalues.h +++ b/include/initialvalues.h @@ -1,35 +1,31 @@ #pragma once // #include +#include "globalvariables.h" #include #include -#include #include -#include "globalvariables.h" +#include using namespace dealii; -namespace EquationData -{ +namespace EquationData { template -class PressureInitialValues : public Function -{ -public: - PressureInitialValues() : Function(1) {} +class PressureInitialValues : public Function { + public: + PressureInitialValues() : Function(1) {} - virtual double value(const Point &p, - const unsigned int component = 0) const; - // virtual void vector_value(const Point &p, - // Vector &value) const; + virtual double value(const Point& p, + const unsigned int component = 0) const; + // virtual void vector_value(const Point &p, + // Vector &value) const; }; template -double -PressureInitialValues::value(const Point &, - const unsigned int) const -{ - return P0; +double PressureInitialValues::value(const Point&, + const unsigned int) const { + return P0; } // template @@ -41,30 +37,28 @@ PressureInitialValues::value(const Point &, // } template -class TemperatureInitialValues : public Function -{ -public: - TemperatureInitialValues() : Function(1) {} +class TemperatureInitialValues : public Function { + public: + TemperatureInitialValues() : Function(1) {} - virtual double value(const Point &p, - const unsigned int component = 0) const; - // virtual void vector_value(const Point &p, - // Vector &value) const; + virtual double value(const Point& p, + const unsigned int component = 0) const; + // virtual void vector_value(const Point &p, + // Vector &value) const; }; template -double -TemperatureInitialValues::value(const Point &, - const unsigned int) const -{ - return T0; +double TemperatureInitialValues::value(const Point&, + const unsigned int) const { + return T0; } // template // void TemperatureInitialValues::vector_value(const Point &p, -// Vector &values) const +// Vector &values) +// const // { // for (unsigned int c = 0; c < this->n_components; ++c) // values(c) = TemperatureInitialValues::value(p, c); // } -} // namespace EquationData \ No newline at end of file +} // namespace EquationData \ No newline at end of file diff --git a/include/sourceterm.h b/include/sourceterm.h index aa7e40c..f7883ef 100644 --- a/include/sourceterm.h +++ b/include/sourceterm.h @@ -1,44 +1,41 @@ #pragma once // #include +#include "globalvariables.h" #include #include -#include #include -#include "globalvariables.h" +#include using namespace dealii; -namespace EquationData -{ +namespace EquationData { template -class PressureSourceTerm : public Function -{ -public: - PressureSourceTerm() : Function(), period(0.2) {} +class PressureSourceTerm : public Function { + public: + PressureSourceTerm() : Function(), period(0.2) {} - virtual double value(const Point &p, - const unsigned int component = 0) const; // assign rhs - // equal to zero - // at each point + virtual double value(const Point& p, + const unsigned int component = 0) const; // assign rhs + // equal to + // zero at each + // point - // virtual void vector_value(const Point &p, - // Vector &value) const; + // virtual void vector_value(const Point &p, + // Vector &value) const; -private: - const double period; + private: + const double period; }; template -double -PressureSourceTerm::value(const Point &p, - const unsigned int component) const -{ - // (void)component; - // Assert(component == 0, ExcIndexRange(component, 0, 1)); // for debug - // Assert(dim == 3, ExcNotImplemented()); - const double time = this->get_time(); // get time - return 0.; +double PressureSourceTerm::value(const Point& p, + const unsigned int component) const { + // (void)component; + // Assert(component == 0, ExcIndexRange(component, 0, 1)); // for debug + // Assert(dim == 3, ExcNotImplemented()); + const double time = this->get_time(); // get time + return 0.; } // template @@ -50,38 +47,37 @@ PressureSourceTerm::value(const Point &p, // } template -class TemperatureSourceTerm : public Function -{ -public: - TemperatureSourceTerm() : Function(), period(0.2) {} +class TemperatureSourceTerm : public Function { + public: + TemperatureSourceTerm() : Function(), period(0.2) {} - virtual double value(const Point& p, - const unsigned int component = 0) const; + virtual double value(const Point& p, + const unsigned int component = 0) const; - // virtual void vector_value(const Point &p, - // Vector &value) const; + // virtual void vector_value(const Point &p, + // Vector &value) const; -private: - const double period; + private: + const double period; }; template -double TemperatureSourceTerm::value(const Point &p, - const unsigned int ) const -{ - // (void)component; - // Assert(component == 0, ExcIndexRange(component, 0, 1)); // for debug - // Assert(dim == 3, ExcNotImplemented()); - // const double time = this->get_time(); // get time - return 0.; +double TemperatureSourceTerm::value(const Point& p, + const unsigned int) const { + // (void)component; + // Assert(component == 0, ExcIndexRange(component, 0, 1)); // for debug + // Assert(dim == 3, ExcNotImplemented()); + // const double time = this->get_time(); // get time + return 0.; } // template // void TemperatureSourceTerm::vector_value(const Point &p, -// Vector &values) const +// Vector &values) +// const // { // for (unsigned int c = 0; c < this->n_components; ++c) // values(c) = TemperatureSourceTerm::value(p, c); // } -} // namespace EquationData \ No newline at end of file +} // namespace EquationData \ No newline at end of file diff --git a/src/main.cc b/src/main.cc index b7d75c8..074e1c6 100644 --- a/src/main.cc +++ b/src/main.cc @@ -1,17 +1,13 @@ -#include "initialvalues.h" #include "boundaryvalues.h" -#include "sourceterm.h" #include "geothermal.h" -int main() -{ - try - { +#include "initialvalues.h" +#include "sourceterm.h" +int main() { + try { using namespace dealii; CoupledTH<3> coupled_TH_solver(1); coupled_TH_solver.run(); - } - catch (std::exception &exc) - { + } catch (std::exception& exc) { std::cerr << std::endl << std::endl << "----------------------------------------------------" @@ -22,9 +18,7 @@ int main() << "----------------------------------------------------" << std::endl; return 1; - } - catch (...) - { + } catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------"