From fa6cef6a704d12f18b9d617a482689e572963735 Mon Sep 17 00:00:00 2001 From: Lucas Sawade Date: Fri, 24 Jan 2025 15:01:01 -0500 Subject: [PATCH] Update all arrays in compute sources to be a function of sources.size() and not nspec. --- CMakeLists.txt | 1 + include/compute/sources/source_medium.hpp | 22 +-- include/compute/sources/sources.hpp | 138 +++++++++++------- src/compute/compute_sources.cpp | 85 ++++++----- tests/unit-tests/assembly/sources/sources.cpp | 17 ++- 5 files changed, 162 insertions(+), 101 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 6ed0a2fc..231d6daa 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -435,6 +435,7 @@ add_library( target_link_libraries( compute + enumerations quadrature mesh # material_class diff --git a/include/compute/sources/source_medium.hpp b/include/compute/sources/source_medium.hpp index e5d0feba..49c5d65d 100644 --- a/include/compute/sources/source_medium.hpp +++ b/include/compute/sources/source_medium.hpp @@ -117,7 +117,7 @@ struct source_medium { template KOKKOS_INLINE_FUNCTION void - store_on_device(const int timestep, const IteratorIndexType &iterator_index, + store_on_device(const int timestep, const IteratorIndexType iterator_index, const PointSourceType &point_source) const { /* For the source it is important to remember that we are using the * mapped index to access the element and source indices @@ -135,24 +135,28 @@ struct source_medium { } } - template - void load_on_host(const int timestep, const IndexType index, + template + void load_on_host(const int timestep, const IteratorIndexType iterator_index, PointSourceType &point_source) const { + const auto index = iterator_index.index; + const auto isource = iterator_index.imap; for (int component = 0; component < components; component++) { point_source.stf(component) = - h_source_time_function(timestep, index.ispec, component); + h_source_time_function(timestep, isource, component); point_source.lagrange_interpolant(component) = - h_source_array(index.ispec, component, index.iz, index.ix); + h_source_array(isource, component, index.iz, index.ix); } } - template - void store_on_host(const int timestep, const IndexType index, + template + void store_on_host(const int timestep, const IteratorIndexType iterator_index, const PointSourceType &point_source) const { + const auto index = iterator_index.index; + const auto isource = iterator_index.imap; for (int component = 0; component < components; component++) { - h_source_time_function(timestep, index.ispec, component) = + h_source_time_function(timestep, isource, component) = point_source.stf(component); - h_source_array(index.ispec, component, index.iz, index.ix) = + h_source_array(isource, component, index.iz, index.ix) = point_source.lagrange_interpolant(component); } } diff --git a/include/compute/sources/sources.hpp b/include/compute/sources/sources.hpp index de9d4417..df422912 100644 --- a/include/compute/sources/sources.hpp +++ b/include/compute/sources/sources.hpp @@ -5,6 +5,7 @@ #include "compute/properties/properties.hpp" #include "enumerations/dimension.hpp" #include "enumerations/material_definitions.hpp" +#include "enumerations/medium.hpp" #include "enumerations/wavefield.hpp" #include "kokkos_abstractions.h" #include "point/sources.hpp" @@ -214,24 +215,25 @@ struct sources { #undef SOURCE_INDICES_VARIABLES_NAME - template + template friend KOKKOS_INLINE_FUNCTION void - load_on_device(const IndexType index, + load_on_device(const IteratorIndexType iterator_index, const specfem::compute::sources &sources, PointSourceType &point_source); - template - friend void load_on_host(const IndexType index, + template + friend void load_on_host(const IteratorIndexType iterator_index, const specfem::compute::sources &sources, PointSourceType &point_source); - template + template friend KOKKOS_INLINE_FUNCTION void - store_on_device(const IndexType index, const PointSourceType &point_source, + store_on_device(const IteratorIndexType iterator_index, + const PointSourceType &point_source, const specfem::compute::sources &sources); - template - friend void store_on_host(const IndexType index, + template + friend void store_on_host(const IteratorIndexType iterator_index, const PointSourceType &point_source, const specfem::compute::sources &sources); }; // namespace compute @@ -278,22 +280,28 @@ load_on_device(const IteratorIndexType iterator_index, "PointSourceType must be an acoustic or elastic point source"); #ifndef NDEBUG + + const int isource = iterator_index.imap; + // Checks if the spectral element index is out of bounds if (index.ispec >= sources.nspec) { Kokkos::abort("Invalid spectral element index detected in source"); } - - // Checks if the spectral element has a source associated with it - if (sources.source_domain_index_mapping(index.ispec) < 0) { - Kokkos::abort("Invalid spectral element index detected in source"); - } - - if (sources.medium_types(index.ispec) != PointSourceType::medium_tag) { + if (sources.medium_types(isource) != PointSourceType::medium_tag) { + std::cout << "compute::sources::load_on_device" << std::endl; + std::cout << " Medium type: " + << specfem::element::to_string(sources.medium_types(isource)) + << std::endl; + std::cout << " PointSourceType medium tag: " + << specfem::element::to_string(PointSourceType::medium_tag) + << std::endl; + std::cout << " File: " << __FILE__ << std::endl; + std::cout << " Line: " << __LINE__ << std::endl; Kokkos::abort("Invalid medium detected in source"); } - if (sources.wavefield_types(index.ispec) != PointSourceType::wavefield_tag) { + if (sources.wavefield_types(isource) != PointSourceType::wavefield_tag) { Kokkos::abort("Invalid wavefield type detected in source"); } #endif @@ -331,12 +339,15 @@ load_on_device(const IteratorIndexType iterator_index, * @param sources Source information for the domain * @param point_source Point source object to load source information into */ -template -void load_on_host(const IndexType index, +template +void load_on_host(const IteratorIndexType iterator_index, const specfem::compute::sources &sources, PointSourceType &point_source) { - static_assert(IndexType::using_simd == false, + // Get the mapping from the iterator index + const auto index = iterator_index.index; + + static_assert(index.using_simd == false, "IndexType must not use SIMD when loading sources"); static_assert( @@ -346,7 +357,7 @@ void load_on_host(const IndexType index, static_assert(PointSourceType::dimension == specfem::dimension::type::dim2, "PointSourceType must be a 2D point source type"); - static_assert(IndexType::dimension == specfem::dimension::type::dim2, + static_assert(index.dimension == specfem::dimension::type::dim2, "IndexType must be a 2D index type"); static_assert( @@ -356,36 +367,38 @@ void load_on_host(const IndexType index, "PointSourceType must be an acoustic or elastic point source"); #ifndef NDEBUG - // Checks if the spectral element index is out of bounds - if (index.ispec >= sources.nspec) { - Kokkos::abort("Invalid spectral element index detected in source"); - } + const int isource = iterator_index.imap; - // Checks if the spectral element has a source associated with it - if (sources.h_source_domain_index_mapping(index.ispec) < 0) { + // Checks if the spectral element index is out of bounds + if ((index.ispec < 0) || (sources.nspec <= index.ispec)) { Kokkos::abort("Invalid spectral element index detected in source"); } - if (sources.h_medium_types(index.ispec) != PointSourceType::medium_tag) { + if (sources.h_medium_types(isource) != PointSourceType::medium_tag) { + std::cout << "compute::sources::load_on_host" << std::endl; + std::cout << " Medium type: " + << specfem::element::to_string(sources.medium_types(isource)) + << std::endl; + std::cout << " PointSourceType medium tag: " + << specfem::element::to_string(PointSourceType::medium_tag) + << std::endl; + std::cout << " File: " << __FILE__ << std::endl; + std::cout << " Line: " << __LINE__ << std::endl; Kokkos::abort("Invalid medium detected in source"); } - if (sources.h_wavefield_types(index.ispec) != - PointSourceType::wavefield_tag) { + if (sources.h_wavefield_types(isource) != PointSourceType::wavefield_tag) { Kokkos::abort("Invalid wavefield type detected in source"); } #endif - IndexType lcoord = index; - lcoord.ispec = sources.h_source_domain_index_mapping(index.ispec); - #define SOURCE_MEDIUM_LOAD_ON_HOST(DIMENSION_TAG, MEDIUM_TAG) \ if constexpr (GET_TAG(DIMENSION_TAG) == specfem::dimension::type::dim2) { \ if constexpr (GET_TAG(MEDIUM_TAG) == PointSourceType::medium_tag) { \ sources \ .CREATE_VARIABLE_NAME(source, GET_NAME(DIMENSION_TAG), \ GET_NAME(MEDIUM_TAG)) \ - .load_on_host(sources.timestep, lcoord, point_source); \ + .load_on_host(sources.timestep, iterator_index, point_source); \ } \ } @@ -440,19 +453,28 @@ store_on_device(const IteratorIndexType iterator_index, "PointSourceType must be an acoustic or elastic point source"); #ifndef NDEBUG - if (index.ispec >= sources.nspec) { - Kokkos::abort("Invalid spectral element index detected in source"); - } + const int isource = iterator_index.imap; - if (sources.source_domain_index_mapping(index.ispec) < 0) { + if ((index.ispec < 0) || (sources.nspec <= index.ispec)) { Kokkos::abort("Invalid spectral element index detected in source"); } - if (sources.medium_types(index.ispec) != PointSourceType::medium_tag) { + if (sources.medium_types(isource) != PointSourceType::medium_tag) { + std::cout << "compute::sources::store_on_device" << std::endl; + std::cout << " isource: " << isource << std::endl; + std::cout << " ispec: " << index.ispec << std::endl; + std::cout << " Medium type: " + << specfem::element::to_string(sources.medium_types(isource)) + << std::endl; + std::cout << " PointSourceType medium tag: " + << specfem::element::to_string(PointSourceType::medium_tag) + << std::endl; + std::cout << " File: " << __FILE__ << std::endl; + std::cout << " Line: " << __LINE__ << std::endl; Kokkos::abort("Invalid medium detected in source"); } - if (sources.wavefield_types(index.ispec) != PointSourceType::wavefield_tag) { + if (sources.wavefield_types(isource) != PointSourceType::wavefield_tag) { Kokkos::abort("Invalid wavefield type detected in source"); } #endif @@ -490,11 +512,14 @@ store_on_device(const IteratorIndexType iterator_index, * @param point_source Point source object to load source information into * @param sources Source information for the domain */ -template -void store_on_host(const IndexType index, const PointSourceType &point_source, +template +void store_on_host(const IteratorIndexType iterator_index, + const PointSourceType &point_source, const specfem::compute::sources &sources) { - static_assert(IndexType::using_simd == false, + const auto index = iterator_index.index; + + static_assert(index.using_simd == false, "IndexType must not use SIMD when storing sources"); static_assert( @@ -504,7 +529,7 @@ void store_on_host(const IndexType index, const PointSourceType &point_source, static_assert(PointSourceType::dimension == specfem::dimension::type::dim2, "PointSourceType must be a 2D point source type"); - static_assert(IndexType::dimension == specfem::dimension::type::dim2, + static_assert(index.dimension == specfem::dimension::type::dim2, "IndexType must be a 2D index type"); static_assert( @@ -514,34 +539,39 @@ void store_on_host(const IndexType index, const PointSourceType &point_source, "PointSourceType must be an acoustic or elastic point source"); #ifndef NDEBUG - if (index.ispec >= sources.nspec) { - Kokkos::abort("Invalid spectral element index detected in source"); - } + const int isource = iterator_index.imap; - if (sources.h_source_domain_index_mapping(index.ispec) < 0) { + if ((index.ispec < 0) || (sources.nspec <= index.ispec)) { Kokkos::abort("Invalid spectral element index detected in source"); } - if (sources.h_medium_types(index.ispec) != PointSourceType::medium_tag) { + if (sources.h_medium_types(isource) != PointSourceType::medium_tag) { + std::cout << "compute::sources::store_on_host" << std::endl; + std::cout << " isource: " << isource << std::endl; + std::cout << " ispec: " << index.ispec << std::endl; + std::cout << " Medium type: " + << specfem::element::to_string(sources.h_medium_types(isource)) + << std::endl; + std::cout << " PointSourceType medium tag: " + << specfem::element::to_string(PointSourceType::medium_tag) + << std::endl; + std::cout << " File: " << __FILE__ << std::endl; + std::cout << " Line: " << __LINE__ << std::endl; Kokkos::abort("Invalid medium detected in source"); } - if (sources.h_wavefield_types(index.ispec) != - PointSourceType::wavefield_tag) { + if (sources.h_wavefield_types(isource) != PointSourceType::wavefield_tag) { Kokkos::abort("Invalid wavefield type detected in source"); } #endif - IndexType lcoord = index; - lcoord.ispec = sources.h_source_domain_index_mapping(index.ispec); - #define SOURCE_MEDIUM_STORE_ON_HOST(DIMENSION_TAG, MEDIUM_TAG) \ if constexpr (GET_TAG(DIMENSION_TAG) == specfem::dimension::type::dim2) { \ if constexpr (GET_TAG(MEDIUM_TAG) == PointSourceType::medium_tag) { \ sources \ .CREATE_VARIABLE_NAME(source, GET_NAME(DIMENSION_TAG), \ GET_NAME(MEDIUM_TAG)) \ - .store_on_host(sources.timestep, lcoord, point_source); \ + .store_on_host(sources.timestep, iterator_index, point_source); \ } \ } diff --git a/src/compute/compute_sources.cpp b/src/compute/compute_sources.cpp index 1f68e562..b6ec7bc7 100644 --- a/src/compute/compute_sources.cpp +++ b/src/compute/compute_sources.cpp @@ -19,17 +19,21 @@ namespace { * FOR ONE SOURCE, THIS WILL NOT HAVE AN IMPACT AT ALL, BUT FOR MANY SOURCES * THIS WILL BECOME A BOTTLENECK. * - * The function runs for every material type and returns vector of sources that - * fall into that material domain + * The function runs for every material type and returns a tuple of two vectors + * - the first vector contains the sources that fall into that material domain + * - the second vector contains the global indices of the sources that */ template -std::vector > sort_sources_per_medium( +std::tuple >, + std::vector > +sort_sources_per_medium( const std::vector > &sources, const specfem::compute::element_types &element_types, const specfem::compute::mesh &mesh) { std::vector > sorted_sources; + std::vector source_indices; // Loop over all sources for (int isource = 0; isource < sources.size(); isource++) { @@ -47,10 +51,11 @@ std::vector > sort_sources_per_medium( // the list of sources and indices if it is. if (element_types.get_medium_tag(lcoord.ispec) == MediumTag) { sorted_sources.push_back(source); + source_indices.push_back(isource); } } - return sorted_sources; + return std::make_tuple(sorted_sources, source_indices); } } // namespace @@ -71,23 +76,15 @@ specfem::compute::sources::sources( h_element_indices(Kokkos::create_mirror_view(element_indices)), source_indices("specfem::sources::indeces", sources.size()), h_source_indices(Kokkos::create_mirror_view(source_indices)), - source_domain_index_mapping( - "specfem::sources::source_domain_index_mapping", nspec), - h_source_domain_index_mapping( - Kokkos::create_mirror_view(source_domain_index_mapping)), - medium_types("specfem::sources::medium_types", nspec), + medium_types("specfem::sources::medium_types", sources.size()), h_medium_types(Kokkos::create_mirror_view(medium_types)), - property_types("specfem::sources::property_types", nspec), + property_types("specfem::sources::property_types", sources.size()), h_property_types(Kokkos::create_mirror_view(property_types)), - boundary_types("specfem::sources::boundary_types", nspec), + boundary_types("specfem::sources::boundary_types", sources.size()), h_boundary_types(Kokkos::create_mirror_view(boundary_types)), - wavefield_types("specfem::sources::wavefield_types", nspec), + wavefield_types("specfem::sources::wavefield_types", sources.size()), h_wavefield_types(Kokkos::create_mirror_view(wavefield_types)) { - for (int ispec = 0; ispec < nspec; ispec++) { - h_source_domain_index_mapping(ispec) = -1; - } - // THERE SHOULD BE LOCATE SOURCES HERE, AND SOURCE SHOULD BE POPULATED // WITH THE LOCAL COORDINATES AND THE GLOBAL ELEMENT INDEX @@ -96,8 +93,10 @@ specfem::compute::sources::sources( // and a vector of indices of the sources in the original sources vector // named source_indices__ #define SORT_SOURCES_PER_MEDIUM(DIMENSION_TAG, MEDIUM_TAG) \ - auto CREATE_VARIABLE_NAME(source, GET_NAME(DIMENSION_TAG), \ - GET_NAME(MEDIUM_TAG)) = \ + auto [CREATE_VARIABLE_NAME(source, GET_NAME(DIMENSION_TAG), \ + GET_NAME(MEDIUM_TAG)), \ + CREATE_VARIABLE_NAME(source_indices, GET_NAME(DIMENSION_TAG), \ + GET_NAME(MEDIUM_TAG))] = \ sort_sources_per_medium( \ sources, element_types, mesh); @@ -108,12 +107,17 @@ specfem::compute::sources::sources( #undef SORT_SOURCES_PER_MEDIUM int nsources = 0; + int nsource_indices = 0; // For a sanity check we count the number of sources and source indices // for each medium and dimension #define COUNT_SOURCES(DIMENSION_TAG, MEDIUM_TAG) \ nsources += CREATE_VARIABLE_NAME(source, GET_NAME(DIMENSION_TAG), \ GET_NAME(MEDIUM_TAG)) \ - .size(); + .size(); \ + nsource_indices += \ + CREATE_VARIABLE_NAME(source_indices, GET_NAME(DIMENSION_TAG), \ + GET_NAME(MEDIUM_TAG)) \ + .size(); CALL_MACRO_FOR_ALL_MEDIUM_TAGS( COUNT_SOURCES, @@ -128,6 +132,12 @@ specfem::compute::sources::sources( throw std::runtime_error( "Not all sources were assigned or sources are assigned multiple times"); } + if (nsources != sources.size()) { + std::cout << "nsources: " << nsources << std::endl; + std::cout << "sources.size(): " << sources.size() << std::endl; + throw std::runtime_error( + "Not all sources were assigned or sources are assigned multiple times"); + } // Reminder we already have // vector current_sources = source__ @@ -138,6 +148,8 @@ specfem::compute::sources::sources( /* Gets the sources and global indices for the current source medium */ \ auto current_sources = CREATE_VARIABLE_NAME( \ source, GET_NAME(DIMENSION_TAG), GET_NAME(MEDIUM_TAG)); \ + auto current_source_indices = CREATE_VARIABLE_NAME( \ + source_indices, GET_NAME(DIMENSION_TAG), GET_NAME(MEDIUM_TAG)); \ /* Loops over the current source*/ \ for (int isource = 0; isource < current_sources.size(); isource++) { \ const auto &source = current_sources[isource]; \ @@ -156,15 +168,16 @@ specfem::compute::sources::sources( * "Multiple sources are detected in the same element"); \ * } \ */ \ - /* source_domain index mapping will be removed */ \ - h_source_domain_index_mapping(ispec) = isource; \ + const int global_isource = current_source_indices[isource]; \ /* setting local source to global element mapping */ \ - h_element_indices(isource) = ispec; \ + h_element_indices(global_isource) = ispec; \ assert(element_types.get_medium_tag(ispec) == GET_TAG(MEDIUM_TAG)); \ - h_medium_types(ispec) = GET_TAG(MEDIUM_TAG); \ - h_property_types(ispec) = element_types.get_property_tag(ispec); \ - h_boundary_types(ispec) = element_types.get_boundary_tag(ispec); \ - h_wavefield_types(ispec) = source->get_wavefield_type(); \ + h_medium_types(global_isource) = GET_TAG(MEDIUM_TAG); \ + h_property_types(global_isource) = \ + element_types.get_property_tag(ispec); \ + h_boundary_types(global_isource) = \ + element_types.get_boundary_tag(ispec); \ + h_wavefield_types(global_isource) = source->get_wavefield_type(); \ } \ this->CREATE_VARIABLE_NAME(source, GET_NAME(DIMENSION_TAG), \ GET_NAME(MEDIUM_TAG)) = \ @@ -194,23 +207,23 @@ specfem::compute::sources::sources( /* Loop over the sources */ \ for (int isource = 0; isource < sources.size(); isource++) { \ int ispec = h_element_indices(isource); \ - if ((h_medium_types(ispec) == GET_TAG(MEDIUM_TAG)) && \ - (h_property_types(ispec) == GET_TAG(PROPERTY_TAG)) && \ - (h_boundary_types(ispec) == GET_TAG(BOUNDARY_TAG))) { \ + if ((h_medium_types(isource) == GET_TAG(MEDIUM_TAG)) && \ + (h_property_types(isource) == GET_TAG(PROPERTY_TAG)) && \ + (h_boundary_types(isource) == GET_TAG(BOUNDARY_TAG))) { \ /* Count the number of sources for each wavefield type */ \ - if (h_wavefield_types(ispec) == \ + if (h_wavefield_types(isource) == \ specfem::wavefield::simulation_field::forward) { \ CREATE_VARIABLE_NAME(count_forward, GET_NAME(DIMENSION_TAG), \ GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG), \ GET_NAME(BOUNDARY_TAG)) \ ++; \ - } else if (h_wavefield_types(ispec) == \ + } else if (h_wavefield_types(isource) == \ specfem::wavefield::simulation_field::backward) { \ CREATE_VARIABLE_NAME(count_backward, GET_NAME(DIMENSION_TAG), \ GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG), \ GET_NAME(BOUNDARY_TAG)) \ ++; \ - } else if (h_wavefield_types(ispec) == \ + } else if (h_wavefield_types(isource) == \ specfem::wavefield::simulation_field::adjoint) { \ CREATE_VARIABLE_NAME(count_adjoint, GET_NAME(DIMENSION_TAG), \ GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG), \ @@ -354,10 +367,10 @@ specfem::compute::sources::sources( /* Loop over all sources */ \ for (int isource = 0; isource < sources.size(); isource++) { \ int ispec = h_element_indices(isource); \ - if ((h_medium_types(ispec) == GET_TAG(MEDIUM_TAG)) && \ - (h_property_types(ispec) == GET_TAG(PROPERTY_TAG)) && \ - (h_boundary_types(ispec) == GET_TAG(BOUNDARY_TAG))) { \ - if (h_wavefield_types(ispec) == \ + if ((h_medium_types(isource) == GET_TAG(MEDIUM_TAG)) && \ + (h_property_types(isource) == GET_TAG(PROPERTY_TAG)) && \ + (h_boundary_types(isource) == GET_TAG(BOUNDARY_TAG))) { \ + if (h_wavefield_types(isource) == \ specfem::wavefield::simulation_field::forward) { \ \ /* Assign global ispec to local forward element index array */ \ diff --git a/tests/unit-tests/assembly/sources/sources.cpp b/tests/unit-tests/assembly/sources/sources.cpp index 82f7b2da..fc2ebb5a 100644 --- a/tests/unit-tests/assembly/sources/sources.cpp +++ b/tests/unit-tests/assembly/sources/sources.cpp @@ -20,10 +20,15 @@ void check_store(specfem::compute::assembly &assembly) { const int ngllz = assembly.mesh.ngllz; const int ngllx = assembly.mesh.ngllx; + std::cout << "Getting the element and source indices" << std::endl; const auto [element_indices, source_indices] = assembly.sources.get_sources_on_device(MediumTag, PropertyTag, BoundaryTag, WavefieldType); + std::cout << "element_indices.size() = " << element_indices.size() + << std::endl; + std::cout << "source_indices.size() = " << source_indices.size() << std::endl; + const int nelements = element_indices.size(); constexpr int num_components = @@ -187,7 +192,9 @@ void check_assembly_source_construction( specfem::point::source; - for (auto &source : sources) { + const int nsources = sources.size(); + for (int isource = 0; isource < nsources; isource++) { + const auto &source = sources[isource]; specfem::point::global_coordinates coord(source->get_x(), source->get_z()); @@ -206,12 +213,18 @@ void check_assembly_source_construction( "stf", 1, components); source->compute_source_time_function(1.0, 0.0, 1, stf); + using mapped_chunk_index_type = + specfem::iterator::impl::mapped_chunk_index_type< + false, specfem::dimension::type::dim2>; for (int iz = 0; iz < ngllz; iz++) { for (int ix = 0; ix < ngllx; ix++) { specfem::point::index index(lcoord.ispec, iz, ix); + const auto mapped_iterator_index = + mapped_chunk_index_type(lcoord.ispec, index, isource); PointSourceType point; - specfem::compute::load_on_host(index, assembly.sources, point); + specfem::compute::load_on_host(mapped_iterator_index, assembly.sources, + point); for (int ic = 0; ic < components; ic++) { const auto lagrange_interpolant = point.lagrange_interpolant(ic);