Skip to content

Commit 3fa4a30

Browse files
michaeldmurphy1claude
michaeldmurphy1
and
claude
committed
Update topologies for Bn radial refinement
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
1 parent 6b151a3 commit 3fa4a30

5 files changed

Lines changed: 342 additions & 34 deletions

File tree

‎src/Domain/CreateInitialElement.cpp‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99
#include "DataStructures/Index.hpp"
1010
#include "DataStructures/IndexIterator.hpp"
1111
#include "Domain/Block.hpp"
12+
#include "Domain/Structure/CreateInitialMesh.hpp"
1213
#include "Domain/Structure/Direction.hpp"
1314
#include "Domain/Structure/Element.hpp"
1415
#include "Domain/Structure/ElementId.hpp"
@@ -228,9 +229,9 @@ Element<VolumeDim> create_initial_element(
228229
compute_element_neighbor_in_same_block(upper_direction));
229230
}
230231
}
232+
const auto topologies = refine_Bn_topology(block.topologies(), element_id);
231233
return Element<VolumeDim>(ElementId<VolumeDim>(element_id),
232-
std::move(neighbors_of_element),
233-
block.topologies());
234+
std::move(neighbors_of_element), topologies);
234235
}
235236
} // namespace domain
236237

‎src/Domain/Structure/CreateInitialMesh.cpp‎

Lines changed: 65 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33

44
#include "Domain/Structure/CreateInitialMesh.hpp"
55

6+
#include <algorithm>
67
#include <array>
78
#include <cstddef>
89
#include <vector>
@@ -106,51 +107,92 @@ std::array<Spectral::Quadrature, Dim> make_quadrature(
106107
}
107108

108109
template <size_t Dim>
109-
bool is_radially_refined_b2(const std::array<domain::Topology, Dim>& topologies,
110-
const ElementId<Dim>& element_id) {
111-
const auto it = alg::find(topologies, domain::Topology::B2Radial);
112-
if (it == topologies.end()) {
113-
return false;
114-
}
115-
return gsl::at(element_id.refinement_levels(),
116-
std::distance(topologies.begin(), it)) != 0;
110+
bool is_angularly_refined(const std::array<domain::Topology, Dim>& topologies,
111+
const ElementId<Dim>& element_id) {
112+
constexpr std::array angular_topologies{
113+
domain::Topology::S1, domain::Topology::B2Angular,
114+
domain::Topology::S2Longitude, domain::Topology::S2Colatitude,
115+
domain::Topology::B3Longitude, domain::Topology::B3Colatitude};
116+
117+
return std::ranges::any_of(
118+
angular_topologies, [&topologies, &element_id](const auto topology) {
119+
const auto it = alg::find(topologies, topology);
120+
121+
return it != topologies.end() and
122+
gsl::at(element_id.refinement_levels(),
123+
std::distance(topologies.begin(), it)) != 0;
124+
});
117125
}
118126
} // namespace
119127

120128
namespace domain {
129+
template <size_t Dim>
130+
std::array<domain::Topology, Dim> refine_Bn_topology(
131+
const std::array<domain::Topology, Dim>& topologies,
132+
const ElementId<Dim>& element_id) {
133+
if constexpr (Dim == 2) {
134+
if (element_id.segment_id(0).index() != 0 and
135+
topologies == domain::topologies::disk) {
136+
return domain::topologies::annulus;
137+
}
138+
} else if constexpr (Dim == 3) {
139+
if (element_id.segment_id(0).index() != 0) {
140+
if (topologies == domain::topologies::cartoon_sphere_inner) {
141+
return domain::topologies::cartoon_sphere;
142+
} else if (topologies == domain::topologies::cartoon_cylinder_inner) {
143+
return domain::topologies::cartoon_cylinder;
144+
} else if (topologies == domain::topologies::full_cylinder) {
145+
return domain::topologies::cylindrical_shell;
146+
} else if (topologies == domain::topologies::full_sphere) {
147+
return domain::topologies::spherical_shell;
148+
}
149+
}
150+
}
151+
return topologies;
152+
}
153+
121154
template <size_t Dim>
122155
Mesh<Dim> create_initial_mesh(
123156
const std::vector<std::array<size_t, Dim>>& initial_extents,
124157
const Element<Dim>& element, const Spectral::Basis i1_basis,
125158
const Spectral::Quadrature i1_quadrature) {
159+
const auto topologies =
160+
refine_Bn_topology(element.topologies(), element.id());
161+
ASSERT(not is_angularly_refined(topologies, element.id()),
162+
"Angular dimensions cannot be angularly refined");
126163
return {initial_extents[element.id().block_id()],
127-
make_basis(element.topologies(), i1_basis),
128-
make_quadrature(element.topologies(), i1_quadrature)};
164+
make_basis(topologies, i1_basis),
165+
make_quadrature(topologies, i1_quadrature)};
129166
}
130167

131168
template <size_t Dim>
132169
Mesh<Dim> create_initial_mesh(
133170
const std::vector<std::array<size_t, Dim>>& initial_extents,
134171
const Block<Dim>& block, [[maybe_unused]] const ElementId<Dim>& element_id,
135172
const Spectral::Basis i1_basis, const Spectral::Quadrature i1_quadrature) {
136-
ASSERT(not is_radially_refined_b2(block.topologies(), element_id),
137-
"Splitting Topology::B2Radial is not yet supported");
138-
return {initial_extents[block.id()], make_basis(block.topologies(), i1_basis),
139-
make_quadrature(block.topologies(), i1_quadrature)};
173+
const auto topologies = refine_Bn_topology(block.topologies(), element_id);
174+
ASSERT(not is_angularly_refined(topologies, element_id),
175+
"Angular dimensions cannot be angularly refined");
176+
return {initial_extents[block.id()], make_basis(topologies, i1_basis),
177+
make_quadrature(topologies, i1_quadrature)};
140178
}
141179
} // namespace domain
142180

143181
#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
144182

145-
#define INSTANTIATE(_, data) \
146-
template Mesh<DIM(data)> domain::create_initial_mesh<DIM(data)>( \
147-
const std::vector<std::array<size_t, DIM(data)>>& initial_extents, \
148-
const Element<DIM(data)>& element, const Spectral::Basis i1_basis, \
149-
const Spectral::Quadrature i1_quadrature); \
150-
template Mesh<DIM(data)> domain::create_initial_mesh<DIM(data)>( \
151-
const std::vector<std::array<size_t, DIM(data)>>& initial_extents, \
152-
const Block<DIM(data)>& block, const ElementId<DIM(data)>& element_id, \
153-
const Spectral::Basis i1_basis, const Spectral::Quadrature _quadrature);
183+
#define INSTANTIATE(_, data) \
184+
template Mesh<DIM(data)> domain::create_initial_mesh<DIM(data)>( \
185+
const std::vector<std::array<size_t, DIM(data)>>& initial_extents, \
186+
const Element<DIM(data)>& element, const Spectral::Basis i1_basis, \
187+
const Spectral::Quadrature i1_quadrature); \
188+
template Mesh<DIM(data)> domain::create_initial_mesh<DIM(data)>( \
189+
const std::vector<std::array<size_t, DIM(data)>>& initial_extents, \
190+
const Block<DIM(data)>& block, const ElementId<DIM(data)>& element_id, \
191+
const Spectral::Basis i1_basis, const Spectral::Quadrature _quadrature); \
192+
template std::array<domain::Topology, DIM(data)> \
193+
domain::refine_Bn_topology<DIM(data)>( \
194+
const std::array<domain::Topology, DIM(data)>& topologies, \
195+
const ElementId<DIM(data)>& element_id);
154196

155197
GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))
156198

‎src/Domain/Structure/CreateInitialMesh.hpp‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,9 @@ template <size_t Dim>
1717
struct ElementId;
1818
template <size_t Dim>
1919
class Mesh;
20+
namespace domain {
21+
enum class Topology : uint8_t;
22+
} // namespace domain
2023
namespace Spectral {
2124
enum class Basis : uint8_t;
2225
enum class Quadrature : uint8_t;
@@ -55,4 +58,16 @@ Mesh<Dim> create_initial_mesh(
5558
const std::vector<std::array<size_t, Dim>>& initial_extents,
5659
const Block<Dim>& block, const ElementId<Dim>& element_id,
5760
Spectral::Basis i1_basis, Spectral::Quadrature i1_quadrature);
61+
62+
/// \ingroup InitializationGroup
63+
/// \brief Create a new array of topologies to account for potential radial
64+
/// refinement of Bn topologies.
65+
///
66+
/// \param topologies initial topologies of the Element
67+
/// \param element_id the ElementId of the Element
68+
template <size_t Dim>
69+
std::array<domain::Topology, Dim> refine_Bn_topology(
70+
const std::array<domain::Topology, Dim>& topologies,
71+
const ElementId<Dim>& element_id);
72+
5873
} // namespace domain

‎tests/Unit/Domain/Structure/Test_CreateInitialMesh.cpp‎

Lines changed: 144 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -118,13 +118,23 @@ SPECTRE_TEST_CASE("Unit.Domain.Structure.CreateInitialMesh", "[Domain][Unit]") {
118118
{Spectral::Basis::Legendre, Spectral::Basis::Chebyshev}) {
119119
for (const auto& i1_quadrature :
120120
{Spectral::Quadrature::GaussLobatto, Spectral::Quadrature::Gauss}) {
121-
CHECK(create_initial_mesh({{{3, 4}}}, disk, element_id_2d, i1_basis,
122-
i1_quadrature) ==
121+
CHECK(create_initial_mesh(
122+
{{{3, 4}}}, disk,
123+
ElementId<2>{0, {{SegmentId{1, 0}, SegmentId{0, 0}}}},
124+
i1_basis, i1_quadrature) ==
123125
Mesh<2>{{{3, 4}},
124126
std::array{Spectral::Basis::ZernikeB2,
125127
Spectral::Basis::ZernikeB2},
126128
std::array{Spectral::Quadrature::GaussRadauUpper,
127129
Spectral::Quadrature::Equiangular}});
130+
CHECK(create_initial_mesh(
131+
{{{3, 4}}}, disk,
132+
ElementId<2>{0, {{SegmentId{1, 1}, SegmentId{0, 0}}}},
133+
i1_basis, i1_quadrature) ==
134+
Mesh<2>{{{3, 4}},
135+
std::array{i1_basis, Spectral::Basis::Fourier},
136+
std::array{i1_quadrature,
137+
Spectral::Quadrature::Equiangular}});
128138
}
129139
}
130140
}
@@ -154,19 +164,32 @@ SPECTRE_TEST_CASE("Unit.Domain.Structure.CreateInitialMesh", "[Domain][Unit]") {
154164
{Spectral::Basis::Legendre, Spectral::Basis::Chebyshev}) {
155165
for (const auto& i1_quadrature :
156166
{Spectral::Quadrature::GaussLobatto, Spectral::Quadrature::Gauss}) {
157-
CHECK(create_initial_mesh({{{3, 2, 4}}}, full_cylinder, element_id_3d,
158-
i1_basis, i1_quadrature) ==
167+
CHECK(create_initial_mesh(
168+
{{{3, 2, 4}}}, full_cylinder,
169+
ElementId<3>{
170+
0, {{SegmentId{1, 0}, SegmentId{0, 0}, SegmentId{0, 0}}}},
171+
i1_basis, i1_quadrature) ==
159172
Mesh<3>{{{3, 2, 4}},
160173
std::array{Spectral::Basis::ZernikeB2,
161174
Spectral::Basis::ZernikeB2, i1_basis},
162175
std::array{Spectral::Quadrature::GaussRadauUpper,
163176
Spectral::Quadrature::Equiangular,
164177
i1_quadrature}});
178+
CHECK(
179+
create_initial_mesh(
180+
{{{3, 2, 4}}}, full_cylinder,
181+
ElementId<3>{
182+
0, {{SegmentId{1, 1}, SegmentId{0, 0}, SegmentId{0, 0}}}},
183+
i1_basis, i1_quadrature) ==
184+
Mesh<3>{{{3, 2, 4}},
185+
std::array{i1_basis, Spectral::Basis::Fourier, i1_basis},
186+
std::array{i1_quadrature, Spectral::Quadrature::Equiangular,
187+
i1_quadrature}});
165188
}
166189
}
167190
}
168191
{
169-
INFO("spheriical_shell");
192+
INFO("spherical_shell");
170193
const Element<3> spherical_shell(element_id_3d, {},
171194
domain::topologies::spherical_shell);
172195
for (const auto& i1_basis :
@@ -191,15 +214,28 @@ SPECTRE_TEST_CASE("Unit.Domain.Structure.CreateInitialMesh", "[Domain][Unit]") {
191214
{Spectral::Basis::Legendre, Spectral::Basis::Chebyshev}) {
192215
for (const auto& i1_quadrature :
193216
{Spectral::Quadrature::GaussLobatto, Spectral::Quadrature::Gauss}) {
194-
CHECK(create_initial_mesh({{{3, 2, 4}}}, full_sphere, element_id_3d,
195-
i1_basis, i1_quadrature) ==
217+
CHECK(create_initial_mesh(
218+
{{{3, 2, 4}}}, full_sphere,
219+
ElementId<3>{
220+
0, {{SegmentId{1, 0}, SegmentId{0, 0}, SegmentId{0, 0}}}},
221+
i1_basis, i1_quadrature) ==
196222
Mesh<3>{{{3, 2, 4}},
197223
std::array{Spectral::Basis::ZernikeB3,
198224
Spectral::Basis::ZernikeB3,
199225
Spectral::Basis::ZernikeB3},
200226
std::array{Spectral::Quadrature::GaussRadauUpper,
201227
Spectral::Quadrature::Gauss,
202228
Spectral::Quadrature::Equiangular}});
229+
CHECK(create_initial_mesh(
230+
{{{3, 2, 4}}}, full_sphere,
231+
ElementId<3>{
232+
0, {{SegmentId{1, 1}, SegmentId{0, 0}, SegmentId{0, 0}}}},
233+
i1_basis, i1_quadrature) ==
234+
Mesh<3>{{{3, 2, 4}},
235+
std::array{i1_basis, Spectral::Basis::SphericalHarmonic,
236+
Spectral::Basis::SphericalHarmonic},
237+
std::array{i1_quadrature, Spectral::Quadrature::Gauss,
238+
Spectral::Quadrature::Equiangular}});
203239
}
204240
}
205241
}
@@ -239,14 +275,113 @@ SPECTRE_TEST_CASE("Unit.Domain.Structure.CreateInitialMesh", "[Domain][Unit]") {
239275
}
240276
}
241277
}
278+
{
279+
const ElementId<3> element_id_on_axis{
280+
0, std::array{SegmentId{1, 0}, SegmentId{0, 0}, SegmentId{0, 0}}};
281+
const ElementId<3> element_id_off_axis{
282+
0, std::array{SegmentId{1, 1}, SegmentId{0, 0}, SegmentId{0, 0}}};
283+
{
284+
INFO("cartoon_sphere_inner");
285+
const Element<3> cartoon_sphere_on_axis(
286+
element_id_on_axis, {}, domain::topologies::cartoon_sphere_inner);
287+
const Element<3> cartoon_sphere_off_axis(
288+
element_id_off_axis, {}, domain::topologies::cartoon_sphere_inner);
289+
for (const auto& i1_basis :
290+
{Spectral::Basis::Legendre, Spectral::Basis::Chebyshev}) {
291+
for (const auto& i1_quadrature : {Spectral::Quadrature::GaussLobatto,
292+
Spectral::Quadrature::Gauss}) {
293+
CHECK(create_initial_mesh({{{3, 1, 1}}}, cartoon_sphere_on_axis,
294+
i1_basis, i1_quadrature) ==
295+
Mesh<3>{{{3, 1, 1}},
296+
std::array{Spectral::Basis::ZernikeB1,
297+
Spectral::Basis::Cartoon,
298+
Spectral::Basis::Cartoon},
299+
std::array{Spectral::Quadrature::GaussRadauUpper,
300+
Spectral::Quadrature::SphericalSymmetry,
301+
Spectral::Quadrature::SphericalSymmetry}});
302+
CHECK(create_initial_mesh({{{3, 1, 1}}}, cartoon_sphere_off_axis,
303+
i1_basis, i1_quadrature) ==
304+
Mesh<3>{{{3, 1, 1}},
305+
std::array{i1_basis, Spectral::Basis::Cartoon,
306+
Spectral::Basis::Cartoon},
307+
std::array{i1_quadrature,
308+
Spectral::Quadrature::SphericalSymmetry,
309+
Spectral::Quadrature::SphericalSymmetry}});
310+
}
311+
}
312+
}
313+
{
314+
INFO("cartoon_cylinder_inner");
315+
const Block<3> cartoon_cylinder_block(
316+
nullptr, 0, {}, "", domain::topologies::cartoon_cylinder_inner);
317+
for (const auto& i1_basis :
318+
{Spectral::Basis::Legendre, Spectral::Basis::Chebyshev}) {
319+
for (const auto& i1_quadrature : {Spectral::Quadrature::GaussLobatto,
320+
Spectral::Quadrature::Gauss}) {
321+
CHECK(create_initial_mesh({{{3, 2, 1}}}, cartoon_cylinder_block,
322+
element_id_on_axis, i1_basis,
323+
i1_quadrature) ==
324+
Mesh<3>{{{3, 2, 1}},
325+
std::array{Spectral::Basis::ZernikeB1, i1_basis,
326+
Spectral::Basis::Cartoon},
327+
std::array{Spectral::Quadrature::GaussRadauUpper,
328+
i1_quadrature,
329+
Spectral::Quadrature::AxialSymmetry}});
330+
CHECK(
331+
create_initial_mesh({{{3, 2, 1}}}, cartoon_cylinder_block,
332+
element_id_off_axis, i1_basis,
333+
i1_quadrature) ==
334+
Mesh<3>{{{3, 2, 1}},
335+
std::array{i1_basis, i1_basis, Spectral::Basis::Cartoon},
336+
std::array{i1_quadrature, i1_quadrature,
337+
Spectral::Quadrature::AxialSymmetry}});
338+
}
339+
}
340+
}
341+
}
342+
242343
#ifdef SPECTRE_DEBUG
243344
CHECK_THROWS_WITH(
244345
create_initial_mesh(
245346
{{{3, 4}}}, Block<2>{nullptr, 0, {}, "", domain::topologies::disk},
246-
ElementId<2>{0, {{SegmentId{1, 0}, SegmentId{0, 0}}}},
347+
ElementId<2>{0, {{SegmentId{1, 0}, SegmentId{1, 0}}}},
348+
Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto),
349+
Catch::Matchers::ContainsSubstring(
350+
"Angular dimensions cannot be angularly refined"));
351+
CHECK_THROWS_WITH(
352+
create_initial_mesh(
353+
{{{3, 4}}}, Block<2>{nullptr, 0, {}, "", domain::topologies::disk},
354+
ElementId<2>{0, {{SegmentId{1, 0}, SegmentId{1, 1}}}},
355+
Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto),
356+
Catch::Matchers::ContainsSubstring(
357+
"Angular dimensions cannot be angularly refined"));
358+
CHECK_THROWS_WITH(
359+
create_initial_mesh(
360+
{{{3, 4, 5}}},
361+
Block<3>{nullptr, 0, {}, "", domain::topologies::full_cylinder},
362+
ElementId<3>{0,
363+
{{SegmentId{0, 0}, SegmentId{1, 1}, SegmentId{0, 0}}}},
364+
Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto),
365+
Catch::Matchers::ContainsSubstring(
366+
"Angular dimensions cannot be angularly refined"));
367+
CHECK_THROWS_WITH(
368+
create_initial_mesh(
369+
{{{3, 4, 5}}},
370+
Block<3>{nullptr, 0, {}, "", domain::topologies::full_sphere},
371+
ElementId<3>{0,
372+
{{SegmentId{1, 0}, SegmentId{1, 1}, SegmentId{0, 0}}}},
373+
Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto),
374+
Catch::Matchers::ContainsSubstring(
375+
"Angular dimensions cannot be angularly refined"));
376+
CHECK_THROWS_WITH(
377+
create_initial_mesh(
378+
{{{3, 4, 5}}},
379+
Block<3>{nullptr, 0, {}, "", domain::topologies::spherical_shell},
380+
ElementId<3>{0,
381+
{{SegmentId{1, 0}, SegmentId{0, 0}, SegmentId{1, 1}}}},
247382
Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto),
248383
Catch::Matchers::ContainsSubstring(
249-
"Splitting Topology::B2Radial is not yet supported"));
384+
"Angular dimensions cannot be angularly refined"));
250385
#endif // SPECTRE_DEBUG
251386
}
252387
} // namespace domain

0 commit comments

Comments
 (0)