ryujin 2.1.1 revision ef0fcd4010d109b860652ece4a7b8963fb7d46b1
geometry_geotiff_profile.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
8#include <compile_time_options.h>
9
11#include "geotiff_reader.h"
12
13namespace ryujin
14{
15 namespace Geometries
16 {
24 template <int dim, typename Callable>
25 class ProfileManifold : public dealii::ChartManifold<dim>
26 {
27 public:
28 ProfileManifold(const Callable &callable)
29 : callable_(callable)
30 {
31 }
32
33 dealii::Point<dim>
34 pull_back(const dealii::Point<dim> &space_point) const final
35 {
36 auto chart_point = space_point;
37
38 if constexpr (dim >= 2) {
39 /* transform y-direction (2D) or z-direction (3D): */
40 chart_point[dim - 1] -= callable_(space_point);
41 }
42
43 return chart_point;
44 }
45
46 dealii::Point<dim>
47 push_forward(const dealii::Point<dim> &chart_point) const final
48 {
49 auto space_point = chart_point;
50
51 if constexpr (dim >= 2) {
52 /* transform y-direction (2D) or z-direction (3D): */
53 space_point[dim - 1] += callable_(space_point);
54 }
55
56 return space_point;
57 }
58
59 std::unique_ptr<dealii::Manifold<dim, dim>> clone() const final
60 {
61 return std::make_unique<ProfileManifold<dim, Callable>>(callable_);
62 }
63
64 private:
65 const Callable callable_;
66 };
67
68
69 template <int dim, typename Callable>
70 ProfileManifold<dim, Callable>
71 make_profile_manifold(const Callable &callable)
72 {
73 return {callable};
74 }
75
76
80 template <int dim>
81 class GeoTIFFProfile : public Geometry<dim>
82 {
83 public:
84 GeoTIFFProfile(const std::string subsection)
85 : Geometry<dim>("geotiff profile", subsection)
86 , geotiff_reader_(subsection + "/geotiff profile")
87 {
88 this->add_parameter("position bottom left",
89 point_left_,
90 "Position of bottom left corner");
91
92 for (unsigned int d = 0; d < dim; ++d)
93 point_right_[d] = 20.0;
94 this->add_parameter(
95 "position top right", point_right_, "Position of top right corner");
96
97 subdivisions_x_ = 1;
98 subdivisions_y_ = 1;
99 subdivisions_z_ = 1;
100 boundary_back_ = Boundary::dirichlet;
101 boundary_bottom_ = Boundary::dirichlet;
102 boundary_front_ = Boundary::dirichlet;
103 boundary_left_ = Boundary::dirichlet;
104 boundary_right_ = Boundary::dirichlet;
105 boundary_top_ = Boundary::dirichlet;
106
107 this->add_parameter("subdivisions x",
108 subdivisions_x_,
109 "number of subdivisions in x direction");
110 this->add_parameter(
111 "boundary condition left",
112 boundary_left_,
113 "Type of boundary condition enforced on the left side of the "
114 "domain (faces with normal in negative x direction)");
115 this->add_parameter(
116 "boundary condition right",
117 boundary_right_,
118 "Type of boundary condition enforced on the right side of the "
119 "domain (faces with normal in positive x direction)");
120
121 if constexpr (dim >= 2) {
122 this->add_parameter("subdivisions y",
123 subdivisions_y_,
124 "number of subdivisions in y direction");
125 this->add_parameter(
126 "boundary condition bottom",
127 boundary_bottom_,
128 "Type of boundary condition enforced on the bottom side of the "
129 "domain (faces with normal in negative y direction)");
130 this->add_parameter(
131 "boundary condition top",
132 boundary_top_,
133 "Type of boundary condition enforced on the top side of the "
134 "domain (faces with normal in positive y direction)");
135 }
136
137 if constexpr (dim == 2) {
138 reference_y_coordinate_ = 0.;
139 this->add_parameter(
140 "reference y coordinate",
141 reference_y_coordinate_,
142 "GeoTIFF: select the value for y-coordinate in 2D. That is, the "
143 "1D profile for the lower boundary is queried from the 2D "
144 "geotiff image at coordinates (x, y=constant)");
145 }
146
147 if constexpr (dim == 3) {
148 this->add_parameter("subdivisions z",
149 subdivisions_z_,
150 "number of subdivisions in z direction");
151 this->add_parameter(
152 "boundary condition back",
153 boundary_back_,
154 "Type of boundary condition enforced on the back side of the "
155 "domain (faces with normal in negative z direction)");
156 this->add_parameter(
157 "boundary condition front",
158 boundary_front_,
159 "Type of boundary condition enforced on the front side of the "
160 "domain (faces with normal in positive z direction)");
161 }
162 }
163
164
166 typename dealii::Triangulation<dim> &triangulation) final
167 {
168 /* create mesh: */
169
170 dealii::Triangulation<dim, dim> tria1;
171 tria1.set_mesh_smoothing(triangulation.get_mesh_smoothing());
172
173 if constexpr (dim == 1) {
174 dealii::GridGenerator::subdivided_hyper_rectangle<dim, dim>(
175 tria1, {subdivisions_x_}, point_left_, point_right_);
176 } else if constexpr (dim == 2) {
177 dealii::GridGenerator::subdivided_hyper_rectangle(
178 tria1,
179 {subdivisions_x_, subdivisions_y_},
180 point_left_,
181 point_right_);
182 } else if constexpr (dim == 3) {
183 dealii::GridGenerator::subdivided_hyper_rectangle(
184 tria1,
185 {subdivisions_x_, subdivisions_y_, subdivisions_z_},
186 point_left_,
187 point_right_);
188 }
189
190 triangulation.copy_triangulation(tria1);
191 triangulation.reset_all_manifolds();
192 /* manifold id 0 for transfinite interpolation manifold */
193 triangulation.set_all_manifold_ids(0);
194
195 /* set boundary and manifold ids: */
196
197 for (auto cell : triangulation.active_cell_iterators()) {
198 for (auto f : cell->face_indices()) {
199 auto face = cell->face(f);
200 if (!face->at_boundary())
201 continue;
202 const auto position = face->center();
203
204 if (position[0] < point_left_[0] + 1.e-8) {
205 face->set_boundary_id(boundary_left_);
206 face->set_manifold_id(dealii::numbers::flat_manifold_id);
207 }
208
209 if (position[0] > point_right_[0] - 1.e-8) {
210 face->set_boundary_id(boundary_right_);
211 face->set_manifold_id(dealii::numbers::flat_manifold_id);
212 }
213
214
215 if constexpr (dim == 2) {
216 if (position[1] < point_left_[1] + 1.e-8) {
217 face->set_boundary_id(boundary_bottom_);
218 /* manifold id 1 for ProfileManifold: */
219 face->set_manifold_id(1);
220 }
221 if (position[1] > point_right_[1] - 1.e-8) {
222 face->set_boundary_id(boundary_top_);
223 face->set_manifold_id(dealii::numbers::flat_manifold_id);
224 }
225 }
226
227 if constexpr (dim == 3) {
228 if (position[1] < point_left_[1] + 1.e-8) {
229 face->set_boundary_id(boundary_bottom_);
230 face->set_manifold_id(dealii::numbers::flat_manifold_id);
231 }
232
233 if (position[1] > point_right_[1] - 1.e-8) {
234 face->set_boundary_id(boundary_top_);
235 face->set_manifold_id(dealii::numbers::flat_manifold_id);
236 }
237
238 /*
239 * The lower boundary at z = point_left_[2] is the profile
240 * manifold ant the upper boundary boundary at z =
241 * point_right_[2] is flat:
242 */
243
244 if (position[2] < point_left_[2] + 1.e-8) {
245 face->set_boundary_id(boundary_back_);
246 /* manifold id 1 for ProfileManifold: */
247 face->set_manifold_id(1);
248 }
249
250 if (position[2] > point_right_[2] - 1.e-8) {
251 face->set_boundary_id(boundary_front_);
252 face->set_manifold_id(dealii::numbers::flat_manifold_id);
253 }
254 }
255 } /*for*/
256 } /*for*/
257
258 const auto profile =
259 make_profile_manifold<dim>([&](dealii::Point<dim> point) {
260 /*
261 *
262 */
263 if constexpr (dim == 1) {
264 return 0.;
265 } else if constexpr (dim == 2) {
266 /*
267 * Set the second coordinate to a constant when querying
268 * height information in 2D.
269 */
270 point[1] = reference_y_coordinate_;
271 return geotiff_reader_.compute_height(point);
272 } else if constexpr (dim == 3) {
273 return geotiff_reader_.compute_height(point);
274 }
275 });
276 triangulation.set_manifold(1, profile);
277
278 dealii::TransfiniteInterpolationManifold<dim> transfinite_interpolation;
279 transfinite_interpolation.initialize(triangulation);
280 triangulation.set_manifold(0, transfinite_interpolation);
281 }
282
283 private:
284 GeoTIFFReader<dim> geotiff_reader_;
285
286 dealii::Point<dim> point_left_;
287 dealii::Point<dim> point_right_;
288
289 double reference_y_coordinate_;
290
291 unsigned int subdivisions_x_;
292 unsigned int subdivisions_y_;
293 unsigned int subdivisions_z_;
294
295 Boundary boundary_back_;
296 Boundary boundary_bottom_;
297 Boundary boundary_front_;
298 Boundary boundary_left_;
299 Boundary boundary_right_;
300 Boundary boundary_top_;
301 };
302 } /* namespace Geometries */
303} /* namespace ryujin */
void create_triangulation(typename dealii::Triangulation< dim > &triangulation) final
GeoTIFFProfile(const std::string subsection)
dealii::Point< dim > pull_back(const dealii::Point< dim > &space_point) const final
dealii::Point< dim > push_forward(const dealii::Point< dim > &chart_point) const final
std::unique_ptr< dealii::Manifold< dim, dim > > clone() const final
ProfileManifold< dim, Callable > make_profile_manifold(const Callable &callable)