Skip to content

Commit

Permalink
fixup comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@AlexCarpenter46
AlexCarpenter46 committed Sep 5, 2024
1 parent 3736877 commit 8bff09c
Show file tree
Hide file tree
Showing 12 changed files with 245 additions and 182 deletions.
148 changes: 92 additions & 56 deletions src/Domain/Creators/BinaryCompactObject.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,24 +114,29 @@ BinaryCompactObject<UseWorldtube>::BinaryCompactObject(
// Determination of parameters for domain construction:
const double tan_half_opening_angle = tan(0.5 * opening_angle_);
translation_ = 0.5 * (x_coord_a_ + x_coord_b_);
if (cube_length.has_value()) {
length_inner_cube_ = cube_length.value();
} else {
length_inner_cube_ = x_coord_a_ - x_coord_b_;
length_inner_cube_ = cube_length.value_or(x_coord_a_ - x_coord_b_);
if (length_inner_cube_ < (x_coord_a_ - x_coord_b_)) {
PARSE_ERROR(
context,
"The cube length should be greater than or equal to the initial "
"separation between the two objects.");
}
length_outer_cube_ =
2.0 * envelope_radius_ / sqrt(2.0 + square(tan_half_opening_angle));

// We chose to handle potential roundoff differences here by using equal
// within roundoff instead of exact equality because the wedge map expects
// exact equality when checking for a zero offset.
if (equal_within_roundoff(length_inner_cube_, (x_coord_a_ - x_coord_b_),
std::numeric_limits<double>::epsilon() * 100.0,
length_inner_cube_)) {
offset_x_coord_a_ = 0.0;
offset_x_coord_b_ = 0.0;
} else {
offset_x_coord_a_ =
x_coord_a_ - (x_coord_a_ + x_coord_b_ + length_inner_cube_) / 2.0;
x_coord_a_ - (x_coord_a_ + x_coord_b_ + length_inner_cube_) * 0.5;
offset_x_coord_b_ =
x_coord_b_ - (x_coord_a_ + x_coord_b_ - length_inner_cube_) / 2.0;
x_coord_b_ - (x_coord_a_ + x_coord_b_ - length_inner_cube_) * 0.5;
}

// Calculate number of blocks
Expand Down Expand Up @@ -172,12 +177,6 @@ BinaryCompactObject<UseWorldtube>::BinaryCompactObject(
"malformed. A recommended radius is:\n"
<< suggested_value);
}
if (length_inner_cube_ < (x_coord_a_ - x_coord_b_)) {
PARSE_ERROR(
context,
"The cube length should be greater than or equal to the initial "
"separation between the two objects.");
}
// The following options are irrelevant if the inner regions are covered
// with simple blocks, so we only check them if object_A_ uses the first
// variant type.
Expand Down Expand Up @@ -225,10 +224,8 @@ BinaryCompactObject<UseWorldtube>::BinaryCompactObject(
"or neither.");
}
}
const bool filled_excision_a =
(not use_single_block_a_ and not is_excised_a_);
const bool filled_excision_b =
(not use_single_block_b_ and not is_excised_b_);
const bool filled_excision_a = not(use_single_block_a_ or is_excised_a_);
const bool filled_excision_b = not(use_single_block_b_ or is_excised_b_);
if ((filled_excision_a or filled_excision_b) and
not equal_within_roundoff(offset_x_coord_a_, 0.0)) {
PARSE_ERROR(context,
Expand Down Expand Up @@ -445,24 +442,42 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
// Each object is surrounded by 6 inner wedges that make a sphere, and
// another 6 outer wedges that transition to a cube.
const auto& object_a = std::get<Object>(object_A_);

Maps maps_center_A =
domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(
object_a.inner_radius, object_a.outer_radius,
inner_sphericity_A, 1.0, length_inner_cube_ / 2.0,
{{offset_x_coord_a_, 0.0, 0.0}}, use_equiangular_map_, false,
{}, object_A_radial_distribution),
translation_A);
Maps maps_cube_A =
domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(
object_a.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_,
1.0, 0.0, length_inner_cube_ / 2.0,
{{offset_x_coord_a_, 0.0, 0.0}}, use_equiangular_map_),
translation_A);
Maps maps_center_A{};
Maps maps_cube_A{};
if (offset_x_coord_a_ == 0.0) {
maps_center_A =
domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(
object_a.inner_radius, object_a.outer_radius,
inner_sphericity_A, 1.0, use_equiangular_map_, std::nullopt,
false, {}, object_A_radial_distribution),
translation_A);
maps_cube_A = domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(object_a.outer_radius,
sqrt(3.0) * 0.5 * length_inner_cube_, 1.0,
0.0, use_equiangular_map_),
translation_A);
} else {
maps_center_A = domain::make_vector_coordinate_map_base<
Frame::BlockLogical, Frame::Inertial, 3>(
sph_wedge_coordinate_maps(
object_a.inner_radius, object_a.outer_radius, inner_sphericity_A,
1.0, use_equiangular_map_,
std::pair<double, std::array<double, 3>>{
{length_inner_cube_ * 0.5}, {{offset_x_coord_a_, 0.0, 0.0}}},
false, {}, object_A_radial_distribution),
translation_A);
maps_cube_A = domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(
object_a.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_, 1.0,
0.0, use_equiangular_map_,
std::pair<double, std::array<double, 3>>{
{length_inner_cube_ * 0.5}, {{offset_x_coord_a_, 0.0, 0.0}}}),
translation_A);
}
std::move(maps_center_A.begin(), maps_center_A.end(),
std::back_inserter(maps));
std::move(maps_cube_A.begin(), maps_cube_A.end(), std::back_inserter(maps));
Expand All @@ -485,23 +500,42 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
// Each object is surrounded by 6 inner wedges that make a sphere, and
// another 6 outer wedges that transition to a cube.
const auto& object_b = std::get<Object>(object_B_);
Maps maps_center_B =
domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(
object_b.inner_radius, object_b.outer_radius,
inner_sphericity_B, 1.0, length_inner_cube_ / 2.0,
{{offset_x_coord_b_, 0.0, 0.0}}, use_equiangular_map_, false,
{}, object_B_radial_distribution),
translation_B);
Maps maps_cube_B =
domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(
object_b.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_,
1.0, 0.0, length_inner_cube_ / 2.0,
{{offset_x_coord_b_, 0.0, 0.0}}, use_equiangular_map_),
translation_B);
Maps maps_center_B{};
Maps maps_cube_B{};
if (offset_x_coord_b_ == 0.0) {
maps_center_B =
domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(
object_b.inner_radius, object_b.outer_radius,
inner_sphericity_B, 1.0, use_equiangular_map_, std::nullopt,
false, {}, object_B_radial_distribution),
translation_B);
maps_cube_B = domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(object_b.outer_radius,
sqrt(3.0) * 0.5 * length_inner_cube_, 1.0,
0.0, use_equiangular_map_),
translation_B);
} else {
maps_center_B = domain::make_vector_coordinate_map_base<
Frame::BlockLogical, Frame::Inertial, 3>(
sph_wedge_coordinate_maps(
object_b.inner_radius, object_b.outer_radius, inner_sphericity_B,
1.0, use_equiangular_map_,
std::pair<double, std::array<double, 3>>{
{length_inner_cube_ * 0.5}, {{offset_x_coord_b_, 0.0, 0.0}}},
false, {}, object_B_radial_distribution),
translation_B);
maps_cube_B = domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(
object_b.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_, 1.0,
0.0, use_equiangular_map_,
std::pair<double, std::array<double, 3>>{
{length_inner_cube_ * 0.5}, {{offset_x_coord_b_, 0.0, 0.0}}}),
translation_B);
}
std::move(maps_center_B.begin(), maps_center_B.end(),
std::back_inserter(maps));
std::move(maps_cube_B.begin(), maps_cube_B.end(), std::back_inserter(maps));
Expand All @@ -528,11 +562,13 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
std::back_inserter(maps));

// --- Outer spherical shell (10 blocks) ---
Maps maps_outer_shell = domain::make_vector_coordinate_map_base<
Frame::BlockLogical, Frame::Inertial, 3>(sph_wedge_coordinate_maps(
envelope_radius_, outer_radius_, 1.0, 1.0, 1.0, {{0.0, 0.0, 0.0}},
use_equiangular_map_, true, {}, {radial_distribution_outer_shell_},
ShellWedges::All, opening_angle_));
Maps maps_outer_shell =
domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
sph_wedge_coordinate_maps(envelope_radius_, outer_radius_, 1.0, 1.0,
use_equiangular_map_, std::nullopt, true,
{}, {radial_distribution_outer_shell_},
ShellWedges::All, opening_angle_));
std::move(maps_outer_shell.begin(), maps_outer_shell.end(),
std::back_inserter(maps));

Expand Down
4 changes: 2 additions & 2 deletions src/Domain/Creators/Sphere.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -308,8 +308,8 @@ Domain<3> Sphere::create_domain() const {
sph_wedge_coordinate_maps(
inner_radius_, outer_radius_,
fill_interior_ ? std::get<InnerCube>(interior_).sphericity : 1.0, 1.0,
1.0, {{0.0, 0.0, 0.0}}, use_equiangular_map_, false,
radial_partitioning_, radial_distribution_, which_wedges_),
use_equiangular_map_, std::nullopt, false, radial_partitioning_,
radial_distribution_, which_wedges_),
compression);

std::unordered_map<std::string, ExcisionSphere<3>> excision_spheres{};
Expand Down
104 changes: 54 additions & 50 deletions src/Domain/DomainHelpers.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -595,9 +595,10 @@ size_t which_wedge_index(const ShellWedges& which_wedges) {
std::vector<domain::CoordinateMaps::Wedge<3>> sph_wedge_coordinate_maps(
const double inner_radius, const double outer_radius,
const double inner_sphericity, const double outer_sphericity,
const double cube_half_length, const std::array<double, 3> focal_offset,
const bool use_equiangular_map, const bool use_half_wedges,
const std::vector<double>& radial_partitioning,
const bool use_equiangular_map,
const std::optional<std::pair<double, std::array<double, 3>>>
offset_options,
const bool use_half_wedges, const std::vector<double>& radial_partitioning,
const std::vector<domain::CoordinateMaps::Distribution>&
radial_distribution,
const ShellWedges which_wedges, const double opening_angle) {
Expand All @@ -618,8 +619,56 @@ std::vector<domain::CoordinateMaps::Wedge<3>> sph_wedge_coordinate_maps(
const size_t number_of_layers = 1 + radial_partitioning.size();
double temp_inner_radius = inner_radius;
double temp_inner_sphericity = inner_sphericity;
const bool zero_offset = (magnitude(focal_offset) == 0);
if (zero_offset) {
if (offset_options.has_value()) {
for (size_t layer_i = 0; layer_i < number_of_layers; layer_i++) {
const auto& radial_distribution_this_layer =
radial_distribution.at(layer_i);
std::optional<double> optional_outer_radius{};
if (outer_sphericity != 0.0) {
optional_outer_radius = outer_radius;
} else {
optional_outer_radius = std::nullopt;
}
// Generate wedges/half-wedges a layer at a time.
std::vector<Wedge3DMap> wedges_for_this_layer{};
if (not use_half_wedges) {
for (size_t face_j = which_wedge_index(which_wedges); face_j < 6;
face_j++) {
wedges_for_this_layer.emplace_back(
temp_inner_radius, optional_outer_radius,
offset_options.value().first, offset_options.value().second,
gsl::at(wedge_orientations, face_j), use_equiangular_map,
Halves::Both, radial_distribution_this_layer);
}
} else {
for (size_t i = 0; i < 4; i++) {
wedges_for_this_layer.emplace_back(
temp_inner_radius, optional_outer_radius,
offset_options.value().first, offset_options.value().second,
gsl::at(wedge_orientations, i), use_equiangular_map,
Halves::LowerOnly, radial_distribution_this_layer);
wedges_for_this_layer.emplace_back(
temp_inner_radius, optional_outer_radius,
offset_options.value().first, offset_options.value().second,
gsl::at(wedge_orientations, i), use_equiangular_map,
Halves::UpperOnly, radial_distribution_this_layer);
}
wedges_for_this_layer.emplace_back(
temp_inner_radius, optional_outer_radius,
offset_options.value().first, offset_options.value().second,
gsl::at(wedge_orientations, 4), use_equiangular_map, Halves::Both,
radial_distribution_this_layer);
wedges_for_this_layer.emplace_back(
temp_inner_radius, optional_outer_radius,
offset_options.value().first, offset_options.value().second,
gsl::at(wedge_orientations, 5), use_equiangular_map, Halves::Both,
radial_distribution_this_layer);
}
for (const auto& wedge : wedges_for_this_layer) {
wedges_for_all_layers.push_back(wedge);
}
}
} else {
double temp_outer_radius{};
for (size_t layer_i = 0; layer_i < number_of_layers; layer_i++) {
const auto& radial_distribution_this_layer =
Expand Down Expand Up @@ -678,51 +727,6 @@ std::vector<domain::CoordinateMaps::Wedge<3>> sph_wedge_coordinate_maps(
temp_inner_sphericity = outer_sphericity;
}
}
} else {
for (size_t layer_i = 0; layer_i < number_of_layers; layer_i++) {
const auto& radial_distribution_this_layer =
radial_distribution.at(layer_i);
std::optional<double> optional_outer_radius{};
if (outer_sphericity != 0.0) {
optional_outer_radius = outer_radius;
} else {
optional_outer_radius = std::nullopt;
}
// Generate wedges/half-wedges a layer at a time.
std::vector<Wedge3DMap> wedges_for_this_layer{};
if (not use_half_wedges) {
for (size_t face_j = which_wedge_index(which_wedges); face_j < 6;
face_j++) {
wedges_for_this_layer.emplace_back(
temp_inner_radius, optional_outer_radius, cube_half_length,
focal_offset, gsl::at(wedge_orientations, face_j),
use_equiangular_map, Halves::Both,
radial_distribution_this_layer);
}
} else {
for (size_t i = 0; i < 4; i++) {
wedges_for_this_layer.emplace_back(
temp_inner_radius, optional_outer_radius, cube_half_length,
focal_offset, gsl::at(wedge_orientations, i), use_equiangular_map,
Halves::LowerOnly, radial_distribution_this_layer);
wedges_for_this_layer.emplace_back(
temp_inner_radius, optional_outer_radius, cube_half_length,
focal_offset, gsl::at(wedge_orientations, i), use_equiangular_map,
Halves::UpperOnly, radial_distribution_this_layer);
}
wedges_for_this_layer.emplace_back(
temp_inner_radius, optional_outer_radius, cube_half_length,
focal_offset, gsl::at(wedge_orientations, 4), use_equiangular_map,
Halves::Both, radial_distribution_this_layer);
wedges_for_this_layer.emplace_back(
temp_inner_radius, optional_outer_radius, cube_half_length,
focal_offset, gsl::at(wedge_orientations, 5), use_equiangular_map,
Halves::Both, radial_distribution_this_layer);
}
for (const auto& wedge : wedges_for_this_layer) {
wedges_for_all_layers.push_back(wedge);
}
}
}
return wedges_for_all_layers;
}
Expand Down
5 changes: 3 additions & 2 deletions src/Domain/DomainHelpers.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,8 +181,9 @@ size_t which_wedge_index(const ShellWedges& which_wedges);
*/
std::vector<domain::CoordinateMaps::Wedge<3>> sph_wedge_coordinate_maps(
double inner_radius, double outer_radius, double inner_sphericity,
double outer_sphericity, double cube_half_length,
std::array<double, 3> focal_offset, bool use_equiangular_map,
double outer_sphericity, bool use_equiangular_map,
std::optional<std::pair<double, std::array<double, 3>>> offset_options =
std::nullopt,
bool use_half_wedges = false,
const std::vector<double>& radial_partitioning = {},
const std::vector<domain::CoordinateMaps::Distribution>&
Expand Down
3 changes: 2 additions & 1 deletion tests/Unit/Domain/Creators/Test_Tags.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,8 @@ void test_center_tags() {
const std::unique_ptr<DomainCreator<3>> domain_creator =
std::make_unique<domain::creators::BinaryCompactObject<false>>(
Object{0.2, 5.0, 8.0, true, true}, Object{0.6, 4.0, -5.5, true, true},
std::array<double, 2>{{0.1, 0.2}}, 100.0, 500.0, 13.5, 1_st, 5_st);
std::array<double, 2>{{0.1, 0.2}}, 100.0, 500.0, std::nullopt, 1_st,
5_st);

const auto grid_center_A =
Tags::ObjectCenter<ObjectLabel::A>::create_from_options(domain_creator);
Expand Down
18 changes: 18 additions & 0 deletions tests/Unit/Domain/Creators/TimeDependentOptions/Test_BinaryCompactObject.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -612,6 +612,24 @@ void test_errors() {
Catch::Matchers::ContainsSubstring(
"Trying to build the shape map for object " +
get_output(domain::ObjectLabel::B)));
if (IsCylindrical) {
CHECK_THROWS_WITH(
([&radii]() {
TimeDependentMapOptions<IsCylindrical> time_dep_opts{
1.0,
ExpMapOptions<IsCylindrical>{{1.0, 0.0}, 0.0, 0.01},
RotMapOptions<IsCylindrical>{{0.0, 0.0, 0.0}},
std::nullopt,
ShapeMapAOptions<IsCylindrical>{8, {}},
std::nullopt};
time_dep_opts.build_maps(
std::array{std::array{5.0, 0.0, 0.0}, std::array{-5.0, 0.0, 0.0}},
{{7.5, 0.0, 0.0}}, {{-7.5, 0.0, 0.0}}, radii, radii, 25.0, 100.0);
}()),
Catch::Matchers::ContainsSubstring(
"When using the CylindricalBinaryCompactObject domain creator, "
"the excision centers cannot be offset."));
}
CHECK_THROWS_WITH(
TimeDependentMapOptions<IsCylindrical>{}
.template grid_to_distorted_map<domain::ObjectLabel::A>(true),
Expand Down
Loading

0 comments on commit 8bff09c

Please sign in to comment.