448 *
template<
int spacedim>
479 *
template <
int dim,
int spacedim, MeshType MESH, InitialConditionType ICTYPE>
506 *
,
double end_time = 0.5);
553 *
const unsigned int degree;
663 *
const double end_time;
672 *
template <
int spacedim>
681 *
const unsigned int component = 0)
const override;
693 *
template <
int spacedim>
695 *
const unsigned int component)
const
716 *
template<
int spacedim, MeshType MESH, InitialConditionType ICTYPE>
760 *
for(
int i = 0; i < 10; ++i){
772 *
const double radius)
777 *
for(
int i = 0; i < 10; ++i){
814 *
const unsigned int component)
const
816 *
if(component == 0){
817 *
if(p.square() <= radius){
837 *
const unsigned int component)
const
839 *
if(component == 0){
841 *
const Point<3> compare(p - center);
842 *
if(compare.square() <= radius){
862 *
const unsigned int component)
const
864 *
if(component == 0){
865 *
const Point<3> center(18.41988074, 0, 0);
866 *
const Point<3> compare(p - center);
867 *
if(compare.square() <= radius){
887 *
const unsigned int component)
const
889 *
if(component == 0){
891 *
const Point<3> compare(p - center);
892 *
if(compare.square() <= radius){
912 *
const unsigned int component)
const
914 *
if(component == 0){
916 *
const Point<3> compare(p - center);
917 *
if(compare.square() <= radius){
937 *
const unsigned int component)
const
939 *
if(component == 0){
942 *
for(
int i=0; i < 10; ++i){
962 *
const unsigned int component)
const
964 *
if(component == 0){
967 *
double width = ((std::atan2(p(1),p(2)) - 3.1415926)/3.1415926)*18.84955592;
968 *
for(
int i=0; i < 10; ++i){
973 *
return x_val*w_val;
988 *
const unsigned int component)
const
990 *
if(component == 0){
993 *
for(
int i=0; i < 10; ++i){
1013 *
const unsigned int component)
const
1015 *
if(component == 0){
1018 *
for(
int i=0; i < 10; ++i){
1038 *
const unsigned int component)
const
1040 *
if(component == 0){
1043 *
double width = ((std::atan2(p(1),p(2)) - 3.1415926)/3.1415926)*18.84955592;
1044 *
for(
int i=0; i < 10; ++i){
1049 *
return x_val*w_val;
1064 *
const unsigned int component)
const
1066 *
if(component == 0){
1082 *
const unsigned int component)
const
1084 *
if(component == 0){
1100 *
const unsigned int component)
const
1102 *
if(component == 0){
1118 *
const unsigned int component)
const
1120 *
if(component == 0){
1136 *
const unsigned int component)
const
1138 *
if(component == 0){
1146 *
template <
int dim,
int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1161 *
template <
int dim,
int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1179 *
, end_time(end_time)
1189 *
template <
int dim,
int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1192 *
dof_handler.distribute_dofs(fe);
1199 * const std::vector<types::global_dof_index> dofs_per_component =
1200 * DoFTools::count_dofs_per_fe_component(dof_handler);
1201 * const unsigned int n_u = dofs_per_component[0],
1202 * n_v = dofs_per_component[1];
1204 * std::cout << "Number of active cells: " << triangulation.n_active_cells()
1206 * << "Total number of cells: " << triangulation.n_cells()
1208 * << "Number of degrees of freedom: " << dof_handler.n_dofs()
1209 * << " (" << n_u << '+
' << n_v << ')
' << std::endl;
1211 * DynamicSparsityPattern dsp(dof_handler.n_dofs());
1213 * DoFTools::make_sparsity_pattern(dof_handler,
1215 * sparsity_pattern.copy_from(dsp);
1217 * system_matrix.reinit(sparsity_pattern);
1219 * solution.reinit(dof_handler.n_dofs());
1220 * old_solution.reinit(dof_handler.n_dofs());
1221 * system_rhs.reinit(dof_handler.n_dofs());
1225 * /** @brief Uses a direct solver to invert the system matrix, then multiplies the RHS vector by the inverted matrix to get the solution.
1226 * * Also includes a timer feature, which is currently commented out, but can be helpful to compute how long a run will take
1227 * * @tparam dim The dimension of the manifold
1228 * * @tparam spacedim The dimension of the ambient space
1229 * * @tparam MESH The type of mesh being used, doesn't
change how this function
works
1232 * template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1233 * void SHEquation<dim, spacedim, MESH, ICTYPE>::solve_time_step()
1237 * std::cout << "Solving linear system" << std::endl;
1244 * SparseDirectUMFPACK direct_solver;
1246 * direct_solver.initialize(system_matrix);
1248 * direct_solver.vmult(solution, system_rhs);
1253 * std::cout << "done (" << timer.cpu_time() << " s)" << std::endl;
1260 * /** @brief Converts the solution vector into a .vtu file and labels the outputs as u and v
1261 * * @tparam dim The dimension of the manifold
1262 * * @tparam spacedim The dimension of the ambient space
1263 * * @tparam MESH The type of mesh being used, doesn't
change how this function
works
1266 * template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1267 * void SHEquation<dim, spacedim, MESH, ICTYPE>::output_results() const
1269 * std::vector<std::string> solution_names(1, "u");
1270 * solution_names.emplace_back("v");
1271 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
1273 * DataComponentInterpretation::component_is_scalar);
1274 * interpretation.push_back(DataComponentInterpretation::component_is_scalar);
1276 * DataOut<dim, spacedim> data_out;
1277 * data_out.add_data_vector(dof_handler,
1280 * interpretation /*,
1281 * DataOut<dim, spacedim>::type_dof_data*/);
1283 * data_out.build_patches(degree + 1);
1287 * Takes the output_file_name string and appends timestep_number with up to three leading 0's
1294 *
data_out.write_vtu(output);
1400 *
template <
int dim,
int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1457 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1467 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1482 *
for(
const auto &cell : dof_handler.active_cell_iterators()){
1486 *
fe_values.reinit(cell);
1488 *
cell->get_dof_indices(local_dof_indices);
1490 *
for(
const unsigned int q_index : fe_values.quadrature_point_indices()){
1492 *
for(
const unsigned int i : fe_values.dof_indices()){
1498 *
const double phi_i_u = fe_values[
u].value(i, q_index);
1500 *
const double phi_i_v = fe_values[v].value(i, q_index);
1503 *
for(
const unsigned int j : fe_values.dof_indices())
1512 *
const double phi_j_v = fe_values[v].value(
j, q_index);
1533 *
for(
unsigned int i : fe_values.dof_indices()){
1534 *
for(
unsigned int j : fe_values.dof_indices()){
1535 *
system_matrix.add(local_dof_indices[i],
1536 *
local_dof_indices[
j],
1547 *
while (time <= end_time)
1579 *
for(
const auto &cell : dof_handler.active_cell_iterators()){
1592 *
fe_values.reinit(cell);
1594 *
cell->get_dof_indices(local_dof_indices);
1601 *
for(
const unsigned int q_index : fe_values.quadrature_point_indices()){
1614 *
for(
const unsigned int i : fe_values.dof_indices()){
1615 *
Un1 +=
old_solution(local_dof_indices[i])*fe_values[
u].value(i, q_index);
1623 *
for(
const unsigned int i : fe_values.dof_indices()){
1625 *
*fe_values[
u].
value(i, q_index)*fe_values.JxW(q_index);
1634 *
for(
unsigned int i : fe_values.dof_indices()){
1676 *
template <
int dim,
int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1706 *
template <
int dim,
int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1728 * mesh is fine enough to resolve this
1731 * GridTools::transform(transform_function<spacedim>, triangulation);
1734 * template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1735 * void SHEquation<dim, spacedim, MESH, ICTYPE>::make_sphere()
1737 * GridGenerator::hyper_sphere(triangulation, Point<3>(0, 0, 0), 18.41988074);
1738 * triangulation.refine_global(refinement_number);
1741 * template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1742 * void SHEquation<dim, spacedim, MESH, ICTYPE>::make_torus()
1744 * GridGenerator::torus(triangulation, 9., 4.);
1745 * triangulation.refine_global(refinement_number);
1747 * template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
1748 * void SHEquation<dim, spacedim, MESH, ICTYPE>::make_hypercube()
1750 * GridGenerator::hyper_cube(triangulation, -18.84955592, 18.84955592);
1751 * triangulation.refine_global(refinement_number);
1753 * } // namespace SwiftHohenbergSolver
1759 * using namespace SwiftHohenbergSolver;
1763 * An array of mesh types. We iterate over this to allow for longer runs without having to stop the code
1766 * MeshType mesh_types[5] = {HYPERCUBE, CYLINDER, SPHERE, TORUS, SINUSOID};
1769 * An array of initial condition types. We iterate this as well, for the same reason
1772 * InitialConditionType ic_types[3] = {HOTSPOT, PSUEDORANDOM, RANDOM};
1776 * Controls how long the code runs
1779 * const double end_time = 100.;
1783 * The number of times we refine the hypercube mesh
1786 * const unsigned int ref_num = 6;
1790 * The timestep will be 1/timestep_denominator
1793 * const unsigned int timestep_denominator = 25;
1797 * Loops over mesh types, then initial condition types, then loops over values of g_1
1800 * for(const auto MESH : mesh_types){
1801 * for(const auto ICTYPE: ic_types){
1802 * for(int i = 0; i < 8; ++i){
1805 * The value of g_1 passed to the solver object
1808 * const double g_constant = 0.2*i;
1812 * Used to distinguish the start of each run
1815 * std::cout<< std::endl << std::endl;
1820 * Switch statement that determines what template parameters are used by the solver object. Template parameters must be known at compile time, so we cannot
1821 * pass this as a variable unfortunately. In each case, we create a filename string (named appropriately for the particular case), output to the console what
1822 * we are running, create the solver object, and call run(). Note that for the cylinder, sphere, and sinusoid we decrease the refinement number by 1. This keeps
1823 * the number of dofs used in these cases comparable to the number of dofs on the 2D hypercube (otherwise the number of dofs is much larger). For the torus, we
1824 * decrease the refinement number by 2.
1833 * std::string filename = "HYPERCUBE-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1834 * std::cout << "Running: " << filename << std::endl << std::endl;
1836 * SHEquation<2, 2, HYPERCUBE, HOTSPOT> heat_equation_solver(1, timestep_denominator,
1837 * ref_num, 0.3, g_constant,
1838 * filename, end_time);
1839 * heat_equation_solver.run();
1843 * case PSUEDORANDOM:
1845 * std::string filename = "HYPERCUBE-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1846 * std::cout << "Running: " << filename << std::endl << std::endl;
1848 * SHEquation<2, 2, HYPERCUBE, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
1849 * ref_num, 0.3, g_constant,
1850 * filename, end_time);
1851 * heat_equation_solver.run();
1857 * std::string filename = "HYPERCUBE-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1858 * std::cout << "Running: " << filename << std::endl << std::endl;
1860 * SHEquation<2, 2, HYPERCUBE, RANDOM> heat_equation_solver(1, timestep_denominator,
1861 * ref_num, 0.3, g_constant,
1862 * filename, end_time);
1863 * heat_equation_solver.run();
1872 * std::string filename = "CYLINDER-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1873 * std::cout << "Running: " << filename << std::endl << std::endl;
1875 * SHEquation<2, 3, CYLINDER, HOTSPOT> heat_equation_solver(1, timestep_denominator,
1876 * ref_num-1, 0.3, g_constant,
1877 * filename, end_time);
1878 * heat_equation_solver.run();
1882 * case PSUEDORANDOM:
1884 * std::string filename = "CYLINDER-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1885 * std::cout << "Running: " << filename << std::endl << std::endl;
1887 * SHEquation<2, 3, CYLINDER, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
1888 * ref_num-1, 0.3, g_constant,
1889 * filename, end_time);
1890 * heat_equation_solver.run();
1896 * std::string filename = "CYLINDER-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1897 * std::cout << "Running: " << filename << std::endl << std::endl;
1899 * SHEquation<2, 3, CYLINDER, RANDOM> heat_equation_solver(1, timestep_denominator,
1900 * ref_num-1, 0.3, g_constant,
1901 * filename, end_time);
1902 * heat_equation_solver.run();
1911 * std::string filename = "SPHERE-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1912 * std::cout << "Running: " << filename << std::endl << std::endl;
1914 * SHEquation<2, 3, SPHERE, HOTSPOT> heat_equation_solver(1, timestep_denominator,
1915 * ref_num-1, 0.3, g_constant,
1916 * filename, end_time);
1917 * heat_equation_solver.run();
1921 * case PSUEDORANDOM:
1923 * std::string filename = "SPHERE-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1924 * std::cout << "Running: " << filename << std::endl << std::endl;
1926 * SHEquation<2, 3, SPHERE, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
1927 * ref_num-1, 0.3, g_constant,
1928 * filename, end_time);
1929 * heat_equation_solver.run();
1935 * std::string filename = "SPHERE-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1936 * std::cout << "Running: " << filename << std::endl << std::endl;
1938 * SHEquation<2, 3, SPHERE, RANDOM> heat_equation_solver(1, timestep_denominator,
1939 * ref_num-1, 0.3, g_constant,
1940 * filename, end_time);
1941 * heat_equation_solver.run();
1950 * std::string filename = "TORUS-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1951 * std::cout << "Running: " << filename << std::endl << std::endl;
1953 * SHEquation<2, 3, TORUS, HOTSPOT> heat_equation_solver(1, timestep_denominator,
1954 * ref_num-2, 0.3, g_constant,
1955 * filename, end_time);
1956 * heat_equation_solver.run();
1960 * case PSUEDORANDOM:
1962 * std::string filename = "TORUS-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1963 * std::cout << "Running: " << filename << std::endl << std::endl;
1965 * SHEquation<2, 3, TORUS, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
1966 * ref_num-2, 0.3, g_constant,
1967 * filename, end_time);
1968 * heat_equation_solver.run();
1974 * std::string filename = "TORUS-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1975 * std::cout << "Running: " << filename << std::endl << std::endl;
1977 * SHEquation<2, 3, TORUS, RANDOM> heat_equation_solver(1, timestep_denominator,
1978 * ref_num-2, 0.3, g_constant,
1979 * filename, end_time);
1980 * heat_equation_solver.run();
1989 * std::string filename = "SINUSOID-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
1990 * std::cout << "Running: " << filename << std::endl << std::endl;
1992 * SHEquation<2, 3, SINUSOID, HOTSPOT> heat_equation_solver(1, timestep_denominator,
1993 * ref_num-1, 0.3, g_constant,
1994 * filename, end_time);
1995 * heat_equation_solver.run();
1999 * case PSUEDORANDOM:
2001 * std::string filename = "SINUSOID-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
2002 * std::cout << "Running: " << filename << std::endl << std::endl;
2004 * SHEquation<2, 3, SINUSOID, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
2005 * ref_num-1, 0.3, g_constant,
2006 * filename, end_time);
2007 * heat_equation_solver.run();
2013 * std::string filename = "SINUSOID-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
2014 * std::cout << "Running: " << filename << std::endl << std::endl;
2016 * SHEquation<2, 3, SINUSOID, RANDOM> heat_equation_solver(1, timestep_denominator,
2017 * ref_num-1, 0.3, g_constant,
2018 * filename, end_time);
2019 * heat_equation_solver.run();
2028 * catch (std::exception &exc)
2030 * std::cout << "An error occurred" << std::endl;
2031 * std::cerr << std::endl
2033 * << "----------------------------------------------------"
2035 * std::cerr << "Exception on processing: " << std::endl
2036 * << exc.what() << std::endl
2037 * << "Aborting!" << std::endl
2038 * << "----------------------------------------------------"
2045 * std::cout << "Error occurred, made it past first catch" << std::endl;
2046 * std::cerr << std::endl
2048 * << "----------------------------------------------------"
2050 * std::cerr << "Unknown exception!" << std::endl
2051 * << "Aborting!" << std::endl
2052 * << "----------------------------------------------------"
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
@ 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 random(DoFHandler< dim, spacedim > &dof_handler)
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
return_type extract_boundary_mesh(const MeshType< dim, spacedim > &volume_mesh, MeshType< dim - 1, spacedim > &surface_mesh, const std::set< types::boundary_id > &boundary_ids=std::set< types::boundary_id >())
void torus(Triangulation< dim, spacedim > &tria, const double centerline_radius, const double inner_radius, const unsigned int n_cells_toroidal=6, const double phi=2.0 *numbers::PI)
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.
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
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)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation