diff --git a/src/Compadre_Basis.hpp b/src/Compadre_Basis.hpp index e4608dcb1..ac846a2d1 100644 --- a/src/Compadre_Basis.hpp +++ b/src/Compadre_Basis.hpp @@ -895,7 +895,8 @@ void createWeightsAndP(const BasisData& data, const member_type& teamMember, scr } // generate weight vector from distances and window sizes - w(i+my_num_neighbors*d) = GMLS::Wab(r, data._epsilons(target_index), data._weighting_type, data._weighting_p, data._weighting_n); + const double cutoff_p = (data._epsilons_from_sources.extent(0) == 0) ? data._epsilons(target_index) : data._epsilons(data._pc.getNeighborIndex(target_index, i)); + w(i+my_num_neighbors*d) = GMLS::Wab(r, cutoff_p, data._weighting_type, data._weighting_p, data._weighting_n); calcPij(data, teamMember, delta.data(), thread_workspace.data(), target_index, i + d*my_num_neighbors, 0 /*alpha*/, dimension, polynomial_order, false /*bool on only specific order*/, V, reconstruction_space, polynomial_sampling_functional); @@ -964,6 +965,7 @@ void createWeightsAndPForCurvature(const BasisData& data, const member_type& tea // ignores V when calculating weights from a point, i.e. uses actual point values double r; + const double cutoff_p = (data._epsilons_from_sources.extent(0) == 0) ? data._epsilons(target_index) : data._epsilons(data._pc.getNeighborIndex(target_index, i)); // get Euclidean distance of scaled relative coordinate from the origin if (V==NULL) { @@ -974,10 +976,10 @@ void createWeightsAndPForCurvature(const BasisData& data, const member_type& tea // generate weight vector from distances and window sizes if (only_specific_order) { - w(i) = GMLS::Wab(r, data._epsilons(target_index), data._curvature_weighting_type, data._curvature_weighting_p, data._curvature_weighting_n); + w(i) = GMLS::Wab(r, cutoff_p, data._curvature_weighting_type, data._curvature_weighting_p, data._curvature_weighting_n); calcPij(data, teamMember, delta.data(), thread_workspace.data(), target_index, i, 0 /*alpha*/, dimension, 1, true /*bool on only specific order*/); } else { - w(i) = GMLS::Wab(r, data._epsilons(target_index), data._curvature_weighting_type, data._curvature_weighting_p, data._curvature_weighting_n); + w(i) = GMLS::Wab(r, cutoff_p, data._curvature_weighting_type, data._curvature_weighting_p, data._curvature_weighting_n); calcPij(data, teamMember, delta.data(), thread_workspace.data(), target_index, i, 0 /*alpha*/, dimension, data._curvature_poly_order, false /*bool on only specific order*/, V); } diff --git a/src/Compadre_Functors.hpp b/src/Compadre_Functors.hpp index ff8f4978c..789fbae3a 100644 --- a/src/Compadre_Functors.hpp +++ b/src/Compadre_Functors.hpp @@ -27,6 +27,7 @@ struct GMLSBasisData { Kokkos::View _source_extra_data; Kokkos::View _target_extra_data; Kokkos::View _epsilons; + Kokkos::View _epsilons_from_sources; Kokkos::View _prestencil_weights; Kokkos::View _curvature_support_operations; Kokkos::View _operations; diff --git a/src/Compadre_GMLS.cpp b/src/Compadre_GMLS.cpp index eb6e438af..8c0865202 100644 --- a/src/Compadre_GMLS.cpp +++ b/src/Compadre_GMLS.cpp @@ -26,6 +26,11 @@ void GMLS::generatePolynomialCoefficients(const int number_of_batches, const boo _h_ss._neighbor_lists = _pc._nla; _h_ss._max_evaluation_sites_per_target = _additional_pc._nla.getMaxNumNeighbors() + 1; + // check that window sizes match number of either target or source sites + compadre_assert_release(this->_epsilons.extent(0)==this->_pc._target_coordinates.extent(0)); + compadre_assert_release(this->_epsilons_from_sources.extent(0)<=0 + || this->_epsilons_from_sources.extent(0)==this->_pc._source_coordinates.extent(0)); + /* * Generate Quadrature */ @@ -526,10 +531,11 @@ const GMLSBasisData GMLS::extractBasisData() const { data._source_extra_data = gmls._source_extra_data; data._target_extra_data = gmls._target_extra_data; data._pc = gmls._pc; - data._epsilons = gmls._epsilons ; - data._prestencil_weights = gmls._prestencil_weights ; + data._epsilons = gmls._epsilons; + data._epsilons_from_sources = gmls._epsilons_from_sources; + data._prestencil_weights = gmls._prestencil_weights; data._additional_pc = gmls._additional_pc; - data._poly_order = gmls._poly_order ; + data._poly_order = gmls._poly_order; data._curvature_poly_order = gmls._curvature_poly_order; data._NP = gmls._NP; data._global_dimensions = gmls._global_dimensions; diff --git a/src/Compadre_GMLS.hpp b/src/Compadre_GMLS.hpp index 753ebd501..1d96f82b1 100644 --- a/src/Compadre_GMLS.hpp +++ b/src/Compadre_GMLS.hpp @@ -104,6 +104,10 @@ friend class Evaluator; //! h supports determined through neighbor search (device) Kokkos::View _epsilons; + + //! optional, h supports determined through neighbor search (device) + //! which is only used in the kernel + Kokkos::View _epsilons_from_sources; //! generated weights for nontraditional samples required to transform data into expected sampling //! functional form (device). @@ -735,6 +739,12 @@ friend class Evaluator; //! Get a view (device) of all window sizes decltype(_epsilons)* getWindowSizes() { return &_epsilons; } + //! Get a view (device) of all window sizes + decltype(_epsilons)* getWindowSizesFromSourceSites() { + compadre_assert_release(this->_epsilons_from_sources.extent(0)>0); + return &_epsilons_from_sources; + } + //! Get a view (device) of all tangent direction bundles. decltype(_T)* getTangentDirections() { return &_T; } @@ -1044,10 +1054,8 @@ friend class Evaluator; //! Sets window sizes, also called the support of the kernel template void setWindowSizes(view_type epsilons) { - // allocate memory on device _epsilons = decltype(_epsilons)("device epsilons", epsilons.extent(0)); - // copy data from host to device Kokkos::deep_copy(_epsilons, epsilons); this->resetCoefficientData(); @@ -1061,6 +1069,28 @@ friend class Evaluator; this->resetCoefficientData(); } + //! (OPTIONAL) + //! Sets window sizes, also called the support of the kernel + template + void setWindowSizesFromSourceSites(view_type epsilons) { + compadre_assert_release(epsilons.extent(0)==this->_pc._source_coordinates.extent(0)); + // allocate memory on device + _epsilons_from_sources = decltype(_epsilons_from_sources)("device epsilons from source sites", epsilons.extent(0)); + // copy data from host to device + Kokkos::deep_copy(_epsilons, epsilons); + this->resetCoefficientData(); + } + + //! (OPTIONAL) + //! Sets window sizes, also called the support of the kernel (device) + template + void setWindowSizesFromSourceSites(decltype(_epsilons_from_sources) epsilons) { + compadre_assert_release(epsilons.extent(0)==this->_pc._source_coordinates.extent(0)); + // allocate memory on device + _epsilons_from_sources = epsilons; + this->resetCoefficientData(); + } + //! (OPTIONAL) //! Sets orthonormal tangent directions for reconstruction on a manifold. The first rank of this 2D array //! corresponds to the target indices, i.e., rows of the neighbor lists 2D array. The second rank is the