detray/builders/material_map_builder.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 // Project include(s).
0012 #include "detray/builders/bin_fillers.hpp"
0013 #include "detray/builders/material_map_factory.hpp"
0014 #include "detray/builders/material_map_generator.hpp"
0015 #include "detray/builders/surface_factory_interface.hpp"
0016 #include "detray/builders/volume_builder_interface.hpp"
0017 #include "detray/geometry/concepts.hpp"
0018 #include "detray/geometry/surface.hpp"
0019 #include "detray/material/material_map.hpp"
0020
0021 // System include(s)
0022 #include <cassert>
0023 #include <map>
0024 #include <memory>
0025 #include <stdexcept>
0026 #include <vector>
0027
0028 namespace detray {
0029
0030 namespace detail {
0031
0032 struct material_coll_size;
0033
0034 template <typename material_t>
0035 struct add_sf_material_map;
0036
0037 }  // namespace detail
0038
0039 /// @brief Build the material maps for a given volume.
0040 ///
0041 /// Decorator class to a volume builder that adds material maps to either
0042 /// surfaces or volumes
0043 template <typename detector_t, std::size_t DIM = 2u,
0044           typename mat_map_factory_t =
0045               material_grid_factory<typename detector_t::algebra_type>>
0046 class material_map_builder final : public volume_decorator<detector_t> {
0047   using material_t = typename detector_t::material;
0048
0049  public:
0050   using scalar_type = dscalar<typename detector_t::algebra_type>;
0051   using detector_type = detector_t;
0052   using value_type = material_slab<scalar_type>;
0053   using bin_data_type =
0054       typename fill_by_bin::template bin_data<DIM, value_type>;
0055
0056   /// @param vol_builder volume builder that should be decorated with material
0057   /// maps
0058   DETRAY_HOST
0059   explicit material_map_builder(
0060       std::unique_ptr<volume_builder_interface<detector_t>> vol_builder)
0061       : volume_decorator<detector_t>(std::move(vol_builder)) {
0062     DETRAY_VERBOSE_HOST("Add material map builder to volume: " << this->name());
0063   }
0064
0065   /// @returns the raw materials and their local bin indices that are
0066   /// currently staged in the builder
0067   DETRAY_HOST
0068   auto data() const -> const std::map<dindex, std::vector<bin_data_type>>& {
0069     return m_bin_data;
0070   }
0071
0072   /// Overwrite, to add material maps in addition to surfaces
0073   /// @{
0074   DETRAY_HOST
0075   void add_surfaces(
0076       std::shared_ptr<surface_factory_interface<detector_t>> sf_factory,
0077       typename detector_t::geometry_context ctx = {}) override {
0078     DETRAY_VERBOSE_HOST("Add [material] surface factory:");
0079
0080     // If the factory also holds surface data, call base volume builder
0081     volume_decorator<detector_t>::add_surfaces(sf_factory, ctx);
0082
0083     // Add material and bin data
0084     using sf_factory_t = material_map_factory<detector_t, axis::multi_bin<DIM>>;
0085     auto mat_factory = std::dynamic_pointer_cast<sf_factory_t>(sf_factory);
0086     if (mat_factory) {
0087       DETRAY_VERBOSE_HOST(
0088           "-> found decoration: " << DETRAY_TYPENAME(sf_factory_t));
0089       (*mat_factory)(this->surfaces(), m_bin_data, m_n_bins, m_axis_spans);
0090     }
0091     auto mat_generator =
0092         std::dynamic_pointer_cast<material_map_generator<detector_t>>(
0093             sf_factory);
0094     if (mat_generator) {
0095       DETRAY_VERBOSE_HOST("-> found decoration: " << DETRAY_TYPENAME(
0096                               material_map_generator<detector_t>));
0097       (*mat_generator)(this->surfaces(), this->masks(), m_bin_data, m_n_bins);
0098       return;
0099     }
0100
0101     if (!mat_factory && !mat_generator) {
0102       DETRAY_VERBOSE_HOST("No material found in this surface factory");
0103       DETRAY_VERBOSE_HOST("-> Built non-material surfaces");
0104     }
0105   }
0106   /// @}
0107
0108   /// Add the volume and the material maps to the detector @param det
0109   DETRAY_HOST
0110   auto build(detector_t& det, typename detector_t::geometry_context ctx = {}) ->
0111       typename detector_t::volume_type* override {
0112     DETRAY_VERBOSE_HOST("Build material maps...");
0113
0114     // Ensure the material links are correct BEFORE the surfaces are built
0115     // and potentially added to an acceleration data structure
0116     update_material_links(det);
0117
0118     DETRAY_DEBUG_HOST(
0119         "-> Let underlying builders construct the volume using correct "
0120         "material links...");
0121     // Construct the surfaces
0122     typename detector_t::volume_type* vol_ptr =
0123         volume_decorator<detector_t>::build(det, ctx);
0124     DETRAY_DEBUG_HOST("-> Back in material map building");
0125
0126     // Build the material grid for every surface that has material
0127     dindex sf_idx = 0u;
0128     auto vol = tracking_volume{det, this->vol_index()};
0129     DETRAY_DEBUG_HOST("Build material maps for surfaces...");
0130     for (const auto& sf_desc : vol.surfaces()) {
0131       if (!surface_has_map(sf_idx)) {
0132         sf_idx++;
0133         continue;
0134       }
0135
0136       DETRAY_DEBUG_HOST("-> surface #" << sf_idx << " sf_desc = " << sf_desc);
0137
0138       // Construct and append the material map for a given surface shape
0139       darray<std::vector<scalar_type>, DIM> axis_spans{};
0140       if (auto axis_spans_itr = m_axis_spans.find(sf_idx);
0141           axis_spans_itr != m_axis_spans.end()) {
0142         axis_spans = axis_spans_itr->second;
0143       }
0144
0145       // Construct and append the material map for a given surface shape
0146       auto sf = geometry::surface{det, sf_desc};
0147       [[maybe_unused]] auto [mat_id, mat_idx] =
0148           sf.template visit_mask<detail::add_sf_material_map<material_t>>(
0149               m_factory, m_bin_data.at(sf_idx), m_n_bins.at(sf_idx), axis_spans,
0150               det._materials);
0151
0152       // Make sure the linking was precomputed correctly
0153       std::stringstream err_msg;
0154       if (mat_id != sf_desc.material().id()) {
0155         err_msg << "-> material id mismatch for surface " << sf_idx
0156                 << ": expected " << sf_desc.material().id() << ", got "
0157                 << mat_id;
0158
0159         DETRAY_FATAL_HOST(err_msg.str());
0160         throw std::runtime_error(err_msg.str());
0161       }
0162       if (mat_idx != sf_desc.material().index()) {
0163         err_msg << "-> material index mismatch for "
0164                    "surface "
0165                 << sf_idx << ": expected " << sf_desc.material().index()
0166                 << ", got " << mat_idx;
0167
0168         DETRAY_FATAL_HOST(err_msg.str());
0169         throw std::runtime_error(err_msg.str());
0170       }
0171       sf_idx++;
0172     }
0173
0174     DETRAY_VERBOSE_HOST(
0175         "Successfully built material maps for volume: " << this->name());
0176
0177     // Give the volume to the next decorator
0178     return vol_ptr;
0179   }
0180
0181  private:
0182   /// Check whether a surface with a given index @param sf_idx should receive
0183   /// material from this builder
0184   bool surface_has_map(const dindex sf_idx) const {
0185     return m_bin_data.contains(sf_idx);
0186   }
0187
0188   /// Set the correct global surface material link for this detector
0189   void update_material_links(const detector_t& det) {
0190     DETRAY_VERBOSE_HOST("Precompute material links...");
0191
0192     // The total number of surfaces that will be built by this builder
0193     const dindex n_surfaces{static_cast<dindex>(this->surfaces().size())};
0194
0195     // Count the number of maps per material type to get the correct indices
0196     std::map<typename material_t::id, dindex> mat_type_count;
0197
0198     for (dindex sf_idx = 0u; sf_idx < n_surfaces; ++sf_idx) {
0199       if (!surface_has_map(sf_idx)) {
0200         continue;
0201       }
0202
0203       auto& sf_desc = this->surfaces().at(sf_idx);
0204       const auto id{sf_desc.material().id()};
0205
0206       DETRAY_DEBUG_HOST("-> surface #" << sf_idx << " (mat. id = " << id
0207                                        << ")");
0208
0209       if (!mat_type_count.contains(id)) {
0210         mat_type_count.emplace(id, 0u);
0211       } else {
0212         mat_type_count.at(id)++;
0213       }
0214
0215       DETRAY_DEBUG_HOST("--> set mat. index to: " << mat_type_count.at(id));
0216
0217       sf_desc.material().set_index(mat_type_count.at(id));
0218     }
0219
0220     // Current sizes of the material stores to get global index offsets
0221     std::map<std::size_t, dindex> size_map;
0222     det._materials.template apply<detail::material_coll_size>(
0223         size_map, std::make_index_sequence<material_t::n_types>{});
0224
0225     DETRAY_DEBUG_HOST("Shift material indices by detector offset...");
0226     for (dindex sf_idx = 0u; sf_idx < n_surfaces; ++sf_idx) {
0227       if (!surface_has_map(sf_idx)) {
0228         continue;
0229       }
0230       auto& sf_desc = this->surfaces().at(sf_idx);
0231
0232       auto coll_idx{static_cast<std::size_t>(sf_desc.material().id())};
0233       sf_desc.material() += size_map.at(coll_idx);
0234
0235       DETRAY_DEBUG_HOST("(Shift = " << size_map.at(coll_idx)
0236                                     << "): material link = "
0237                                     << sf_desc.material());
0238     }
0239   }
0240
0241   /// The surface this material map belongs to (index is volume local)
0242   std::map<dindex, std::vector<bin_data_type>> m_bin_data;
0243   /// Number of bins for the material grid axes
0244   std::map<dindex, darray<std::size_t, DIM>> m_n_bins{};
0245   /// The Axis spans for the material grid axes
0246   std::map<dindex, darray<std::vector<scalar_type>, DIM>> m_axis_spans{};
0247   /// Helper to generate empty grids
0248   mat_map_factory_t m_factory{};
0249 };
0250
0251 namespace detail {
0252
0253 /// A functor to obtain the material collection sizes for every collection
0254 struct material_coll_size {
0255   template <typename... coll_ts, std::size_t... I>
0256   DETRAY_HOST inline void operator()(std::map<std::size_t, dindex>& size_map,
0257                                      std::index_sequence<I...> /*seq*/,
0258                                      const coll_ts&... coll) const {
0259     (size_map.emplace(I, static_cast<dindex>(std::size(coll))), ...);
0260   }
0261 };
0262
0263 /// A functor to add a material map to a surface
0264 template <typename material_t>
0265 struct add_sf_material_map {
0266   template <typename mask_coll_t, typename index_range_t,
0267             typename mat_factory_t, typename bin_data_t, std::size_t DIM,
0268             typename material_store_t, concepts::scalar scalar_t>
0269   DETRAY_HOST inline std::pair<typename material_t::id, dindex> operator()(
0270       [[maybe_unused]] const mask_coll_t& mask_coll,
0271       [[maybe_unused]] const index_range_t& index,
0272       [[maybe_unused]] const mat_factory_t& mat_factory,
0273       [[maybe_unused]] std::vector<bin_data_t>& bin_data,
0274       [[maybe_unused]] const darray<std::size_t, DIM>& n_bins,
0275       [[maybe_unused]] const darray<std::vector<scalar_t>, DIM>& axis_spans,
0276       [[maybe_unused]] material_store_t& mat_store) const {
0277     using mask_t = typename mask_coll_t::value_type;
0278
0279     // No material maps for line surfaces
0280     if constexpr (!concepts::line_object<mask_t> && mask_t::shape::dim == DIM) {
0281       // Map a grid onto the surface mask (the boundaries are taken from
0282       // the @c axis_spans variable, if it is not empty)
0283       mask_t sf_mask = {};
0284       if constexpr (concepts::interval<index_range_t>) {
0285         using index_t = typename index_range_t::index_type;
0286
0287         // Find the true surface extent over all masks
0288         sf_mask = mask_coll.at(index.lower());
0289
0290         if (index.size() > 1u) {
0291           const index_range_t other_masks{
0292               index.lower() + 1u, static_cast<index_t>(index.size() - 1u)};
0293
0294           // Merge sub-masks
0295           for (const auto& sub_mask :
0296                detray::ranges::subrange(mask_coll, other_masks)) {
0297             sf_mask = sf_mask + sub_mask;
0298           }
0299         }
0300       } else {
0301         sf_mask = mask_coll.at(index);
0302       }
0303
0304       auto mat_grid = mat_factory.new_grid(sf_mask, n_bins, {}, {}, axis_spans);
0305
0306       // The detector only knows the non-owning grid types
0307       using non_owning_t = typename decltype(mat_grid)::template type<false>;
0308
0309       // Not every mask shape might be used for material maps
0310       if constexpr (types::contains<material_t, non_owning_t>) {
0311         DETRAY_VERBOSE_HOST("Filling material grid...");
0312
0313         // Add the material slabs to the grid
0314         for (const auto& bin : bin_data) {
0315           mat_grid.template populate<replace<>>(bin.local_bin_idx,
0316                                                 bin.single_element);
0317         }
0318
0319         // Add the material grid to the detector
0320         constexpr auto gid{types::id<material_t, non_owning_t>};
0321         mat_store.template push_back<gid>(mat_grid);
0322         DETRAY_VERBOSE_HOST("Built material grid:" << gid << ":\n"
0323                                                    << mat_grid.axes());
0324
0325         // Return the index of the new material map
0326         return {gid, static_cast<dindex>(mat_store.template size<gid>() - 1u)};
0327       } else {
0328         return {material_t::id::e_none, dindex_invalid};
0329       }
0330     } else {
0331       return {material_t::id::e_none, dindex_invalid};
0332     }
0333   }
0334 };
0335
0336 }  // namespace detail
0337
0338 }  // namespace detray