ryujin 2.1.1 revision 350e54cc11f3d21282bddcf3e6517944dbc508bf
geometry_rectangular_domain.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2022 - 2023 by the ryujin authors
4//
5
6#pragma once
7
9
10namespace ryujin
11{
12 namespace Geometries
13 {
34 template <int dim>
35 class RectangularDomain : public Geometry<dim>
36 {
37 public:
38 RectangularDomain(const std::string subsection)
39 : Geometry<dim>("rectangular domain", subsection)
40 {
41 this->add_parameter("position bottom left",
42 point_left_,
43 "Position of bottom left corner");
44
45 for (unsigned int d = 0; d < dim; ++d)
46 point_right_[d] = 20.0;
47 this->add_parameter(
48 "position top right", point_right_, "Position of top right corner");
49
50 grading_push_forward_ = dim == 1 ? "x" : (dim == 2 ? "x;y" : "x;y;z;");
51 this->add_parameter("grading push forward",
52 grading_push_forward_,
53 "push forward of grading manifold");
54
55 grading_pull_back_ = grading_push_forward_;
56 this->add_parameter("grading pull back",
57 grading_pull_back_,
58 "pull back of grading manifold");
59
60 subdivisions_x_ = 1;
61 subdivisions_y_ = 1;
62 subdivisions_z_ = 1;
63 boundary_back_ = Boundary::dirichlet;
64 boundary_bottom_ = Boundary::dirichlet;
65 boundary_front_ = Boundary::dirichlet;
66 boundary_left_ = Boundary::dirichlet;
67 boundary_right_ = Boundary::dirichlet;
68 boundary_top_ = Boundary::dirichlet;
69
70 this->add_parameter("subdivisions x",
71 subdivisions_x_,
72 "number of subdivisions in x direction");
73 this->add_parameter(
74 "boundary condition left",
75 boundary_left_,
76 "Type of boundary condition enforced on the left side of the "
77 "domain (faces with normal in negative x direction)");
78 this->add_parameter(
79 "boundary condition right",
80 boundary_right_,
81 "Type of boundary condition enforced on the right side of the "
82 "domain (faces with normal in positive x direction)");
83
84 if constexpr (dim >= 2) {
85 this->add_parameter("subdivisions y",
86 subdivisions_y_,
87 "number of subdivisions in y direction");
88 this->add_parameter(
89 "boundary condition bottom",
90 boundary_bottom_,
91 "Type of boundary condition enforced on the bottom side of the "
92 "domain (faces with normal in negative y direction)");
93 this->add_parameter(
94 "boundary condition top",
95 boundary_top_,
96 "Type of boundary condition enforced on the top side of the "
97 "domain (faces with normal in positive y direction)");
98 }
99
100 if constexpr (dim == 3) {
101 this->add_parameter("subdivisions z",
102 subdivisions_z_,
103 "number of subdivisions in z direction");
104 this->add_parameter(
105 "boundary condition back",
106 boundary_back_,
107 "Type of boundary condition enforced on the back side of the "
108 "domain (faces with normal in negative z direction)");
109 this->add_parameter(
110 "boundary condition front",
111 boundary_front_,
112 "Type of boundary condition enforced on the front side of the "
113 "domain (faces with normal in positive z direction)");
114 }
115 }
116
117
118 void create_triangulation(dealii::Triangulation<dim> &triangulation) final
119 {
120 /* create mesh: */
121
122 dealii::Triangulation<dim, dim> tria1;
123 tria1.set_mesh_smoothing(triangulation.get_mesh_smoothing());
124
125 if constexpr (dim == 1) {
126 dealii::GridGenerator::subdivided_hyper_rectangle<dim, dim>(
127 tria1, {subdivisions_x_}, point_left_, point_right_);
128 } else if constexpr (dim == 2) {
129 dealii::GridGenerator::subdivided_hyper_rectangle(
130 tria1,
131 {subdivisions_x_, subdivisions_y_},
132 point_left_,
133 point_right_);
134 } else if constexpr (dim == 3) {
135 dealii::GridGenerator::subdivided_hyper_rectangle(
136 tria1,
137 {subdivisions_x_, subdivisions_y_, subdivisions_z_},
138 point_left_,
139 point_right_);
140 }
141
142 triangulation.copy_triangulation(tria1);
143
144 /* create grading: */
145
146 // FIXME: We should ideally check the push forward and pull back
147 // for compatiblity.
148 if (grading_push_forward_ != grading_pull_back_) {
149 dealii::FunctionManifold<dim> grading(grading_push_forward_,
150 grading_pull_back_);
151 triangulation.set_all_manifold_ids(1);
152 triangulation.set_manifold(1, grading);
153 }
154
155 /* set boundary ids: */
156
157 for (auto cell : triangulation.active_cell_iterators()) {
158 for (auto f : cell->face_indices()) {
159 auto face = cell->face(f);
160 if (!face->at_boundary())
161 continue;
162 const auto position = face->center();
163
164 if (position[0] < point_left_[0] + 1.e-8)
165 face->set_boundary_id(boundary_left_);
166 if (position[0] > point_right_[0] - 1.e-8)
167 face->set_boundary_id(boundary_right_);
168
169 if constexpr (dim >= 2) {
170 if (position[1] < point_left_[1] + 1.e-8)
171 face->set_boundary_id(boundary_bottom_);
172 if (position[1] > point_right_[1] - 1.e-8)
173 face->set_boundary_id(boundary_top_);
174 }
175
176 if constexpr (dim == 3) {
177 if (position[2] < point_left_[2] + 1.e-8)
178 face->set_boundary_id(boundary_back_);
179 if (position[2] > point_right_[2] - 1.e-8)
180 face->set_boundary_id(boundary_front_);
181 }
182 } /*for*/
183 } /*for*/
184
185 std::vector<int> directions;
186
187 if (boundary_left_ == Boundary::periodic ||
188 boundary_right_ == Boundary::periodic) {
189 AssertThrow(boundary_left_ == boundary_right_,
190 dealii::ExcMessage(
191 "For prescribing periodic boundaries in x-direction, "
192 "both, the left and right boundary conditions must "
193 "be set to periodic"));
194 directions.push_back(0);
195 }
196
197 if (dim >= 2 && (boundary_bottom_ == Boundary::periodic ||
198 boundary_top_ == Boundary::periodic)) {
199 AssertThrow(boundary_bottom_ == boundary_top_,
200 dealii::ExcMessage(
201 "For prescribing periodic boundaries in y-direction, "
202 "both, the bottom and top boundary conditions must "
203 "be set to periodic"));
204 directions.push_back(1);
205 }
206
207 if (dim == 3 && (boundary_back_ == Boundary::periodic ||
208 boundary_front_ == Boundary::periodic)) {
209 AssertThrow(boundary_back_ == boundary_front_,
210 dealii::ExcMessage(
211 "For prescribing periodic boundaries in z-direction, "
212 "both, the back and front boundary conditions must "
213 "be set to periodic"));
214 directions.push_back(2);
215 }
216
217#ifndef BUG_COLLECT_PERIODIC_FACES_INSTANTIATION
218 if (!directions.empty()) {
219 std::vector<dealii::GridTools::PeriodicFacePair<
220 typename dealii::Triangulation<dim>::cell_iterator>>
221 periodic_faces;
222
223 for (const auto direction : directions)
224 dealii::GridTools::collect_periodic_faces(
225 triangulation,
226 /*b_id */ Boundary::periodic,
227 /*direction*/ direction,
228 periodic_faces);
229
230 triangulation.add_periodicity(periodic_faces);
231 }
232#endif
233 }
234
235 private:
236 dealii::Point<dim> point_left_;
237 dealii::Point<dim> point_right_;
238
239 unsigned int subdivisions_x_;
240 unsigned int subdivisions_y_;
241 unsigned int subdivisions_z_;
242
243 std::string grading_push_forward_;
244 std::string grading_pull_back_;
245
246 Boundary boundary_back_;
247 Boundary boundary_bottom_;
248 Boundary boundary_front_;
249 Boundary boundary_left_;
250 Boundary boundary_right_;
251 Boundary boundary_top_;
252 };
253 } /* namespace Geometries */
254} /* namespace ryujin */
void create_triangulation(dealii::Triangulation< dim > &triangulation) final
RectangularDomain(const std::string subsection)