8#include <compile_time_options.h>
24 template <
int dim,
typename Callable>
34 pull_back(
const dealii::Point<dim> &space_point)
const final
36 auto chart_point = space_point;
38 if constexpr (dim >= 2) {
40 chart_point[dim - 1] -= callable_(space_point);
49 auto space_point = chart_point;
51 if constexpr (dim >= 2) {
53 space_point[dim - 1] += callable_(space_point);
59 std::unique_ptr<dealii::Manifold<dim, dim>>
clone() const final
61 return std::make_unique<ProfileManifold<dim, Callable>>(callable_);
65 const Callable callable_;
69 template <
int dim,
typename Callable>
70 ProfileManifold<dim, Callable>
85 :
Geometry<dim>(
"geotiff profile", subsection)
86 , geotiff_reader_(subsection +
"/geotiff profile")
88 this->add_parameter(
"position bottom left",
90 "Position of bottom left corner");
92 for (
unsigned int d = 0; d < dim; ++d)
93 point_right_[d] = 20.0;
95 "position top right", point_right_,
"Position of top right corner");
107 this->add_parameter(
"subdivisions x",
109 "number of subdivisions in x direction");
111 "boundary condition left",
113 "Type of boundary condition enforced on the left side of the "
114 "domain (faces with normal in negative x direction)");
116 "boundary condition right",
118 "Type of boundary condition enforced on the right side of the "
119 "domain (faces with normal in positive x direction)");
121 if constexpr (dim >= 2) {
122 this->add_parameter(
"subdivisions y",
124 "number of subdivisions in y direction");
126 "boundary condition bottom",
128 "Type of boundary condition enforced on the bottom side of the "
129 "domain (faces with normal in negative y direction)");
131 "boundary condition top",
133 "Type of boundary condition enforced on the top side of the "
134 "domain (faces with normal in positive y direction)");
137 if constexpr (dim == 2) {
138 reference_y_coordinate_ = 0.;
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)");
147 if constexpr (dim == 3) {
148 this->add_parameter(
"subdivisions z",
150 "number of subdivisions in z direction");
152 "boundary condition back",
154 "Type of boundary condition enforced on the back side of the "
155 "domain (faces with normal in negative z direction)");
157 "boundary condition front",
159 "Type of boundary condition enforced on the front side of the "
160 "domain (faces with normal in positive z direction)");
166 typename dealii::Triangulation<dim> &triangulation)
final
170 dealii::Triangulation<dim, dim> tria1;
171 tria1.set_mesh_smoothing(triangulation.get_mesh_smoothing());
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(
179 {subdivisions_x_, subdivisions_y_},
182 }
else if constexpr (dim == 3) {
183 dealii::GridGenerator::subdivided_hyper_rectangle(
185 {subdivisions_x_, subdivisions_y_, subdivisions_z_},
190 triangulation.copy_triangulation(tria1);
191 triangulation.reset_all_manifolds();
193 triangulation.set_all_manifold_ids(0);
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())
202 const auto position = face->center();
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);
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);
215 if constexpr (dim == 2) {
216 if (position[1] < point_left_[1] + 1.e-8) {
217 face->set_boundary_id(boundary_bottom_);
219 face->set_manifold_id(1);
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);
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);
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);
244 if (position[2] < point_left_[2] + 1.e-8) {
245 face->set_boundary_id(boundary_back_);
247 face->set_manifold_id(1);
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);
259 make_profile_manifold<dim>([&](dealii::Point<dim> point) {
263 if constexpr (dim == 1) {
265 }
else if constexpr (dim == 2) {
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);
276 triangulation.set_manifold(1, profile);
278 dealii::TransfiniteInterpolationManifold<dim> transfinite_interpolation;
279 transfinite_interpolation.initialize(triangulation);
280 triangulation.set_manifold(0, transfinite_interpolation);
286 dealii::Point<dim> point_left_;
287 dealii::Point<dim> point_right_;
289 double reference_y_coordinate_;
291 unsigned int subdivisions_x_;
292 unsigned int subdivisions_y_;
293 unsigned int subdivisions_z_;
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
ProfileManifold(const Callable &callable)
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)