diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 049f2ed8f13525a47bb4070cb1338c6a820c8858..612d508f3bc9948726daf556bd630c65de579295 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -77,6 +77,8 @@ if ( CUDA_ON_BACKEND STREQUAL "HIP" AND HIP_FOUND )
else()
add_executable(numerics ${OPENFPM_INIT_FILE} ${CUDA_SOURCES}
+ regression/regression_test.cpp
+ regression/poly_levelset_test.cpp
OdeIntegrators/tests/OdeIntegratores_base_tests.cpp
OdeIntegrators/tests/OdeIntegrator_grid_tests.cpp
DCPSE/DCPSE_op/tests/DCPSE_op_subset_test.cpp
@@ -105,14 +107,16 @@ else()
Operators/Vector/vector_dist_operators_unit_tests.cpp
Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.cpp
../../src/lib/pdata.cpp
- DCPSE/DCPSE_op/tests/DCPSE_op_Surface_tests.cpp
-# BoundaryConditions/tests/method_of_images_cylinder_unit_test.cpp
-# level_set/closest_point/closest_point_unit_tests.cpp
- level_set/redistancing_Sussman/tests/redistancingSussman_fast_unit_test.cpp
-# level_set/redistancing_Sussman/tests/help_functions_unit_test.cpp
- level_set/redistancing_Sussman/tests/narrowBand_unit_test.cpp
-# level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp
-# level_set/redistancing_Sussman/tests/convergence_test.cpp
+ BoundaryConditions/tests/method_of_images_cylinder_unit_test.cpp
+# level_set/closest_point/closest_point_unit_tests.cpp
+# level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp
+# level_set/redistancing_Sussman/tests/convergence_test.cpp
+ level_set/particle_cp/pcp_unit_tests.cpp
+ DCPSE/DCPSE_op/tests/DCPSE_op_Surface_tests.cpp
+# BoundaryConditions/tests/method_of_images_cylinder_unit_test.cpp
+ level_set/redistancing_Sussman/tests/redistancingSussman_fast_unit_test.cpp
+# level_set/redistancing_Sussman/tests/help_functions_unit_test.cpp
+ level_set/redistancing_Sussman/tests/narrowBand_unit_test.cpp
)
set_property(TARGET numerics PROPERTY CUDA_ARCHITECTURES OFF)
@@ -168,6 +172,7 @@ target_include_directories (numerics PUBLIC ${BLITZ_ROOT}/include)
target_include_directories (numerics PUBLIC ${ALGOIM_ROOT}/include)
target_include_directories (numerics PUBLIC ${ALPAKA_ROOT}/include)
target_include_directories (numerics PUBLIC ${MPI_C_INCLUDE_DIRS})
+target_include_directories (numerics PUBLIC ${MINTER_ROOT}/include)
if(EIGEN3_FOUND)
target_include_directories (numerics PUBLIC /usr/local/include)
target_include_directories (numerics PUBLIC ${EIGEN3_INCLUDE_DIR})
@@ -181,6 +186,7 @@ target_link_libraries(numerics -L${LIBHILBERT_LIBRARY_DIRS} ${LIBHILBERT_LIBRARI
target_link_libraries(numerics ${Vc_LIBRARIES})
target_link_libraries(numerics ${MPI_C_LIBRARIES})
target_link_libraries(numerics ${MPI_CXX_LIBRARIES})
+
if(PETSC_FOUND)
target_include_directories (numerics PUBLIC ${PETSC_INCLUDES})
target_link_libraries(numerics ${PETSC_LIBRARIES})
@@ -209,6 +215,10 @@ if (TEST_COVERAGE)
target_link_libraries(numerics -lgcov)
endif()
+## Regression
+target_include_directories(numerics PUBLIC ${MINTER_ROOT}/include)
+##
+
# Request that particles be built with -std=c++11
# As this is a public compile feature anything that links to particles
# will also build with -std=c++11
@@ -353,9 +363,13 @@ install(FILES DMatrix/EMatrix.hpp
DESTINATION openfpm_numerics/include/DMatrix
COMPONENT OpenFPM)
+install(FILES level_set/particle_cp/particle_cp.hpp
+ DESTINATION openfpm_numerics/include/level_set/particle_cp
+ COMPONENT OpenFPM)
+
install(FILES level_set/closest_point/closest_point.hpp
- DESTINATION openfpm_numerics/include/level_set/closest_point
- COMPONENT OpenFPM)
+ DESTINATION openfpm_numerics/include/level_set/closest_point
+ COMPONENT OpenFPM)
#if(BUILD_TESTING)
diff --git a/src/DCPSE/MonomialBasis.hpp b/src/DCPSE/MonomialBasis.hpp
index 381b975aabf878431222aa0a10aa7aa44c7a411e..dad5866b30da7f21974a9f26f700c6a8ed672a8a 100644
--- a/src/DCPSE/MonomialBasis.hpp
+++ b/src/DCPSE/MonomialBasis.hpp
@@ -23,6 +23,8 @@ public:
MonomialBasis(const vector_type<unsigned int, Args...> °rees, unsigned int convergenceOrder);
+ MonomialBasis(unsigned int orderLimit);
+
MonomialBasis(unsigned int degrees[dim], unsigned int convergenceOrder);
// explicit MonomialBasis(Point<dim, unsigned int> degrees, unsigned int convergenceOrder);
@@ -60,6 +62,8 @@ public:
return lhs;
}
+
+
private:
void generateBasis(vector_type<unsigned int, Args...> m, unsigned int r);
};
@@ -138,7 +142,6 @@ void MonomialBasis<dim, T, vector_type, Args...>::generateBasis(vector_type<unsi
alphaMin = 0;
}
//std::cout<<"AlphaMin: "<<alphaMin<<std::endl;
- //unsigned char alphaMin = 0; // we want to always have 1 in the basis
while (it.isNext())
{
diff --git a/src/DCPSE/SupportBuilder.hpp b/src/DCPSE/SupportBuilder.hpp
index 8d83f7f8a6f590963b4e87d9152379ab6985c376..50c69e2684eccb11329385607399e8fd3837efa1 100644
--- a/src/DCPSE/SupportBuilder.hpp
+++ b/src/DCPSE/SupportBuilder.hpp
@@ -19,7 +19,8 @@ enum support_options
N_PARTICLES,
RADIUS,
LOAD,
- ADAPTIVE
+ ADAPTIVE,
+ AT_LEAST_N_PARTICLES
};
diff --git a/src/level_set/particle_cp/particle_cp.hpp b/src/level_set/particle_cp/particle_cp.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c6f344f41e0c9f1b38dde1cca61aa03b4f5d904d
--- /dev/null
+++ b/src/level_set/particle_cp/particle_cp.hpp
@@ -0,0 +1,598 @@
+//
+// Created by lschulze on 2021-10-12
+//
+/**
+ * @file particle_cp.hpp
+ * @class particle_cp_redistancing
+ *
+ * @brief Class for reinitializing a level-set function into a signed distance function using the Closest-Point
+ * method on particles.
+ *
+ * @details The redistancing scheme here is based on Saye, Robert. "High-order methods for computing distances
+ * to implicitly defined surfaces." Communications in Applied Mathematics and Computational Science 9.1 (2014):
+ * 107-141.
+ * Within this file, it has been extended so that it also works on arbitrarily distributed particles. It mainly
+ * consists of these steps: Firstly, the particle distribution is assessed and a subset is determined, which will
+ * provide a set of sample points. For these, an interpolation of the old signed-distance function values in their
+ * support yields a polynomial whose zero level-set reconstructs the surface. The interpolation polynomials are
+ * then used to create a collection of sample points on the interface. Finally, a Newton-optimisation is performed
+ * in order to determine the closest point on the surface to a given query point, allowing for an update of the
+ * signed distance function.
+ *
+ *
+ * @author Lennart J. Schulze
+ * @date October 2021
+ */
+
+#include <math.h>
+#include "Vector/vector_dist.hpp"
+#include "regression/regression.hpp"
+
+template<unsigned int dim, unsigned int n_c> using particles_surface = vector_dist<dim, double, aggregate<int, int, double, double[dim], double[n_c]>>;
+struct Redist_options
+{
+ size_t max_iter = 1000; // params for the Newton algorithm for the solution of the constrained
+ double tolerance = 1e-11; // optimization problem, and for the projection of the sample points on the interface
+
+ double H; // interparticle spacing
+ double r_cutoff_factor; // radius of neighbors to consider during interpolation, as factor of H
+ double sampling_radius;
+
+ int compute_normals = 0;
+ int compute_curvatures = 0;
+
+ int write_sdf = 1;
+ int write_cp = 0;
+
+ unsigned int minter_poly_degree = 4;
+ float minter_lp_degree = 1.0;
+
+ int verbose = 0;
+ int min_num_particles = 0; // if 0, cutoff radius is used, else the int is the minimum number of particles for the regression.
+ // If min_num_particles is used, the cutoff radius is used for the cell list. It is sensible to use
+ // a smaller rcut for the celllist, to promote symmetric supports (otherwise it is possible that the
+ // particle is at the corner of a cell and hence only has neighbors in certain directions)
+ float r_cutoff_factor_min_num_particles; // this is the rcut for the celllist
+ int only_narrowband = 1; // only redistance particles with phi < sampling_radius, or all particles if only_narrowband = 0
+};
+
+template <typename particles_in_type, size_t phi_field, size_t closest_point_field, size_t normal_field, size_t curvature_field, unsigned int num_minter_coeffs>
+class particle_cp_redistancing
+{
+public:
+ particle_cp_redistancing(particles_in_type & vd, Redist_options &redistOptions) : redistOptions(redistOptions),
+ vd_in(vd),
+ vd_s(vd.getDecomposition(), 0),
+ r_cutoff2(redistOptions.r_cutoff_factor*redistOptions.r_cutoff_factor*redistOptions.H*redistOptions.H),
+ minterModelpcp(RegressionModel<dim, vd_s_sdf>(
+ redistOptions.minter_poly_degree, redistOptions.minter_lp_degree))
+ {
+ // Constructor
+ }
+
+ void run_redistancing()
+ {
+ if (redistOptions.verbose)
+ {
+ std::cout<<"Verbose mode. Make sure the vd.getProp<4>(a) is an integer that pcp can write surface flags onto."<<std::endl;
+ }
+
+ detect_surface_particles();
+
+ interpolate_sdf_field();
+
+ find_closest_point();
+ }
+
+private:
+ static constexpr size_t num_neibs = 0;
+ static constexpr size_t vd_s_close_part = 1;
+ static constexpr size_t vd_s_sdf = 2;
+ static constexpr size_t vd_s_sample = 3;
+ static constexpr size_t minter_coeff = 4;
+ // static constexpr size_t vd_s_minter_model = 5;
+ static constexpr size_t vd_in_sdf = phi_field; // this is really required in the vd_in vector, so users need to know about it.
+ static constexpr size_t vd_in_close_part = 4; // this is not needed by the method, but more for debugging purposes, as it shows all particles for which
+ // interpolation and sampling is performed.
+ static constexpr size_t vd_in_normal = normal_field;
+ static constexpr size_t vd_in_curvature = curvature_field;
+ static constexpr size_t vd_in_cp = closest_point_field;
+
+ Redist_options redistOptions;
+ particles_in_type &vd_in;
+ static constexpr unsigned int dim = particles_in_type::dims;
+ static constexpr unsigned int n_c = num_minter_coeffs;
+
+ int dim_r = dim;
+ int n_c_r = n_c;
+
+ particles_surface<dim, n_c> vd_s;
+ double r_cutoff2;
+ RegressionModel<dim, vd_s_sdf> minterModelpcp;
+
+ int return_sign(double phi)
+ {
+ if (phi > 0) return 1;
+ if (phi < 0) return -1;
+ return 0;
+ }
+
+ // this function lays the groundwork for the interpolation step. It classifies particles according to two possible classes:
+ // (a) close particles: These particles are close to the surface and will be the central particle in the subsequent interpolation step.
+ // After the interpolation step, they will also be projected on the surface. The resulting position on the surface will be referred to as a sample point,
+ // which will provide the initial guesses for the final Newton optimization, which finds the exact closest-point towards any given query point.
+ // Further, the number of neighbors in the support of close particles are stored, to efficiently initialize the Vandermonde matrix in the interpolation.
+ // (b) surface particles: These particles have particles in their support radius, which are close particles (a). In order to compute the interpolation
+ // polynomials for (a), these particles need to be stored. An alternative would be to find the relevant particles in the original particle vector.
+
+ void detect_surface_particles()
+ {
+ vd_in.template ghost_get<vd_in_sdf>();
+
+ auto NN = vd_in.getCellList(sqrt(r_cutoff2) + redistOptions.H);
+ auto part = vd_in.getDomainIterator();
+
+ while (part.isNext())
+ {
+ vect_dist_key_dx akey = part.get();
+ // depending on the application this can spare computational effort
+ if ((redistOptions.only_narrowband) && (std::abs(vd_in.template getProp<vd_in_sdf>(akey)) > redistOptions.sampling_radius))
+ {
+ ++part;
+ continue;
+ }
+ int surfaceflag = 0;
+ int sgn_a = return_sign(vd_in.template getProp<vd_in_sdf>(akey));
+ Point<dim,double> xa = vd_in.getPos(akey);
+ int num_neibs_a = 0;
+ double min_sdf = abs(vd_in.template getProp<vd_in_sdf>(akey));
+ vect_dist_key_dx min_sdf_key = akey;
+ if (redistOptions.verbose) vd_in.template getProp<vd_in_close_part>(akey) = 0;
+ int isclose = 0;
+
+ auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(xa));
+ while (Np.isNext())
+ {
+ vect_dist_key_dx bkey = Np.get();
+ int sgn_b = return_sign(vd_in.template getProp<vd_in_sdf>(bkey));
+ Point<dim, double> xb = vd_in.getPos(bkey);
+ Point<dim,double> dr = xa - xb;
+ double r2 = norm2(dr);
+
+ // check if the particle will provide a polynomial interpolation and a sample point on the interface
+ if ((sqrt(r2) < (1.5*redistOptions.H)) && (sgn_a != sgn_b)) isclose = 1;
+
+ // count how many particles are in the neighborhood and will be part of the interpolation
+ if (r2 < r_cutoff2) ++num_neibs_a;
+
+ // decide whether this particle will be required for any interpolation. This is the case if it has a different sign in its slightly extended neighborhood
+ // its up for debate if this is stable or not. Possibly, this could be avoided by simply using vd_in in the interpolation step.
+ if ((sqrt(r2) < ((redistOptions.r_cutoff_factor + 1.5)*redistOptions.H)) && (sgn_a != sgn_b)) surfaceflag = 1;
+ ++Np;
+
+ }
+
+ if (surfaceflag) // these particles will play a role in the subsequent interpolation: Either they carry interpolation polynomials,
+ { // or they will be taken into account in interpolation
+ vd_s.add();
+ for(int k = 0; k < dim; k++) vd_s.getLastPos()[k] = vd_in.getPos(akey)[k];
+ vd_s.template getLastProp<vd_s_sdf>() = vd_in.template getProp<vd_in_sdf>(akey);
+ vd_s.template getLastProp<num_neibs>() = num_neibs_a;
+
+ if (isclose) // these guys will carry an interpolation polynomial and a resulting sample point
+ {
+ vd_s.template getLastProp<vd_s_close_part>() = 1;
+ if (redistOptions.verbose) vd_in.template getProp<vd_in_close_part>(akey) = 1;
+ }
+ else // these particles will not carry an interpolation polynomial
+ {
+ vd_s.template getLastProp<vd_s_close_part>() = 0;
+ if (redistOptions.verbose) vd_in.template getProp<vd_in_close_part>(akey) = 0;
+ }
+ }
+ ++part;
+ }
+ }
+
+ void interpolate_sdf_field()
+ {
+ int message_insufficient_support = 0;
+ int message_projection_fail = 0;
+
+ vd_s.template ghost_get<vd_s_sdf>();
+ double r_cutoff_celllist = sqrt(r_cutoff2);
+ if (redistOptions.min_num_particles != 0) r_cutoff_celllist = redistOptions.r_cutoff_factor_min_num_particles*redistOptions.H;
+ auto NN_s = vd_s.getCellList(r_cutoff_celllist);
+ auto part = vd_s.getDomainIterator();
+
+ // iterate over particles that will get an interpolation polynomial and generate a sample point
+ while (part.isNext())
+ {
+ vect_dist_key_dx a = part.get();
+
+ // only the close particles (a) will get the full treatment (interpolation + projection)
+ if (vd_s.template getProp<vd_s_close_part>(a) != 1)
+ {
+ ++part;
+ continue;
+ }
+
+ const int num_neibs_a = vd_s.template getProp<num_neibs>(a);
+ Point<dim, double> xa = vd_s.getPos(a);
+ int neib = 0;
+ int k_project = 0;
+
+ if(redistOptions.min_num_particles == 0)
+ {
+ auto regSupport = RegressionSupport<decltype(vd_s), decltype(NN_s)>(vd_s, part, sqrt(r_cutoff2), RADIUS, NN_s);
+ if (regSupport.getNumParticles() < n_c) message_insufficient_support = 1;
+ minterModelpcp.computeCoeffs(vd_s, regSupport);
+ }
+ else
+ {
+ auto regSupport = RegressionSupport<decltype(vd_s), decltype(NN_s)>(vd_s, part, n_c + 3, AT_LEAST_N_PARTICLES, NN_s);
+ if (regSupport.getNumParticles() < n_c) message_insufficient_support = 1;
+ minterModelpcp.computeCoeffs(vd_s, regSupport);
+ }
+
+ auto& minterModel = minterModelpcp.model;
+ for(int k = 0; k < n_c; k++) vd_s.template getProp<minter_coeff>(a)[k] = minterModel->getCoeffs()[k];
+
+ double grad_p_minter_mag2;
+
+ EMatrix<double, Eigen::Dynamic, 1> grad_p_minter(dim_r, 1);
+ EMatrix<double, Eigen::Dynamic, 1> x_minter(dim_r, 1);
+ for(int k = 0; k < dim_r; k++) x_minter[k] = xa[k];
+
+ double p_minter = get_p_minter(x_minter, minterModel);
+ k_project = 0;
+ while ((abs(p_minter) > redistOptions.tolerance) && (k_project < redistOptions.max_iter))
+ {
+ grad_p_minter = get_grad_p_minter(x_minter, minterModel);
+ grad_p_minter_mag2 = grad_p_minter.dot(grad_p_minter);
+ x_minter = x_minter - p_minter*grad_p_minter/grad_p_minter_mag2;
+ p_minter = get_p_minter(x_minter, minterModel);
+ //incr = sqrt(p*p*dpdx*dpdx/(grad_p_mag2*grad_p_mag2) + p*p*dpdy*dpdy/(grad_p_mag2*grad_p_mag2));
+ ++k_project;
+ }
+ // store the resulting sample point on the surface at the central particle.
+ for(int k = 0; k < dim; k++) vd_s.template getProp<vd_s_sample>(a)[k] = x_minter[k];
+ if (k_project == redistOptions.max_iter)
+ {
+ if (redistOptions.verbose) std::cout<<"didnt work for "<<a.getKey()<<std::endl;
+ message_projection_fail = 1;
+ }
+
+ ++part;
+ }
+
+ if (message_insufficient_support) std::cout<<"Warning: less number of neighbours than required for interpolation"
+ <<" for some particles. Consider using at least N particles function."<<std::endl;
+ if (message_projection_fail) std::cout<<"Warning: Newton-style projections towards the interface do not satisfy"
+ <<" given tolerance for some particles"<<std::endl;
+
+ }
+
+ // This function now finds the exact closest point on the surface to any given particle.
+ // For this, it iterates through all particles in the input vector (since all need to be redistanced),
+ // and finds the closest sample point on the surface. Then, it uses this sample point as an initial guess
+ // for a subsequent optimization problem, which finds the exact closest point on the surface.
+ // The optimization problem is formulated as follows: Find the point in space with the minimal distance to
+ // the given query particle, with the constraint that it needs to lie on the surface. It is realized with a
+ // Lagrange multiplier that ensures that the solution lies on the surface, i.e. p(x)=0.
+ // The polynomial that is used for checking if the constraint is fulfilled is the interpolation polynomial
+ // carried by the sample point.
+
+ void find_closest_point()
+ {
+ // iterate over all particles, i.e. do closest point optimisation for all particles, and initialize
+ // all relevant variables.
+
+ vd_s.template ghost_get<vd_s_close_part,vd_s_sample,minter_coeff>();
+
+ auto NN_s = vd_s.getCellList(redistOptions.sampling_radius);
+ auto part = vd_in.getDomainIterator();
+
+ int message_step_limitation = 0;
+ int message_convergence_problem = 0;
+
+ // do iteration over all query particles
+ while (part.isNext())
+ {
+ vect_dist_key_dx a = part.get();
+
+ if ((redistOptions.only_narrowband) && (std::abs(vd_in.template getProp<vd_in_sdf>(a)) > redistOptions.sampling_radius))
+ {
+ ++part;
+ continue;
+ }
+
+ // initialise all variables specific to the query particle
+ EMatrix<double, Eigen::Dynamic, 1> xa(dim_r, 1);
+ for(int k = 0; k < dim; k++) xa[k] = vd_in.getPos(a)[k];
+
+ double val;
+ // spatial variable that is optimized
+ EMatrix<double, Eigen::Dynamic, 1> x(dim_r, 1);
+ // Increment variable containing the increment of spatial variables + 1 Lagrange multiplier
+ EMatrix<double, Eigen::Dynamic, 1> dx(dim_r + 1, 1);
+ for(int k = 0; k < (dim + 1); k++) dx[k] = 1.0;
+
+ // Lagrange multiplier lambda
+ double lambda = 0.0;
+ // values regarding the gradient of the Lagrangian and the Hessian of the Lagrangian, which contain
+ // derivatives of the constraint function also.
+ double p = 0;
+ EMatrix<double, Eigen::Dynamic, 1> grad_p(dim_r, 1);
+ for(int k = 0; k < dim; k++) grad_p[k] = 0.0;
+
+ EMatrix<double, Eigen::Dynamic, Eigen::Dynamic> H_p(dim_r, dim_r);
+ for(int k = 0; k < dim; k++)
+ {
+ for(int l = 0; l < dim; l++) H_p(k, l) = 0.0;
+ }
+
+ EMatrix<double, Eigen::Dynamic, 1> c(n_c_r, 1);
+ for (int k = 0; k < n_c; k++) c[k] = 0.0;
+
+ // f(x, lambda) is the Lagrangian, initialize its gradient and Hessian
+ EMatrix<double, Eigen::Dynamic, 1> nabla_f(dim_r + 1, 1);
+ EMatrix<double, Eigen::Dynamic, Eigen::Dynamic> H_f(dim_r + 1, dim_r + 1);
+ for(int k = 0; k < (dim + 1); k++)
+ {
+ nabla_f[k] = 0.0;
+ for(int l = 0; l < (dim + 1); l++) H_f(k, l) = 0.0;
+ }
+
+ // Initialise variables for searching the closest sample point. Note that this is not a very elegant
+ // solution to initialise another x_a vector, but it is done because of the two different types EMatrix and
+ // point vector, and can probably be made nicer.
+ Point<dim, double> xaa = vd_in.getPos(a);
+
+ // Now we will iterate over the sample points, which means iterating over vd_s.
+ auto Np = NN_s.template getNNIterator<NO_CHECK>(NN_s.getCell(vd_in.getPos(a)));
+
+ vect_dist_key_dx b_min = get_closest_neighbor<decltype(NN_s)>(xaa, vd_s, NN_s);
+
+ // set x0 to the sample point which was closest to the query particle
+ for(int k = 0; k < dim; k++) x[k] = vd_s.template getProp<vd_s_sample>(b_min)[k];
+
+ // these two variables are mainly required for optional subroutines in the optimization procedure and for
+ // warnings if the solution strays far from the neighborhood.
+ EMatrix<double, Eigen::Dynamic, 1> x00(dim_r,1);
+ for(int k = 0; k < dim; k++) x00[k] = x[k];
+
+ EMatrix<double, Eigen::Dynamic, 1> x00x(dim_r,1);
+ for(int k = 0; k < dim; k++) x00x[k] = 0.0;
+
+ auto& model = minterModelpcp.model;
+ EVectorXd temp(n_c, 1);
+ for(int k = 0; k < n_c; k++) temp[k] = vd_s.template getProp<minter_coeff>(b_min)[k];
+ Point<dim, int> derivOrder;
+ Point<dim, double> grad_p_minter;
+ model->setCoeffs(temp);
+
+ if(redistOptions.verbose)
+ {
+ std::cout<<std::setprecision(16)<<"VERBOSE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for particle "<<a.getKey()<<std::endl;
+ std::cout<<"x_poly: "<<vd_s.getPos(b_min)[0]<<", "<<vd_s.getPos(b_min)[1]<<"\nxa: "<<xa[0]<<", "<<xa[1]<<"\nx_0: "<<x[0]<<", "<<x[1]<<"\nc: "<<c<<std::endl;
+
+ std::cout<<"interpol_i(x_0) = "<<get_p_minter(x, model)<<std::endl;
+ }
+
+ // variable containing updated distance at iteration k
+ EMatrix<double, Eigen::Dynamic, 1> xax = x - xa;
+ // iteration counter k
+ int k_newton = 0;
+
+ // calculations needed specifically for k_newton == 0
+ p = get_p_minter(x, model);
+ grad_p = get_grad_p_minter(x, model);
+
+ // this guess for the Lagrange multiplier is taken from the original paper by Saye and can be done since
+ // p is approximately zero at the sample point.
+ lambda = -xax.dot(grad_p)/grad_p.dot(grad_p);
+
+ for(int k = 0; k < dim; k++) nabla_f[k] = xax[k] + lambda*grad_p[k];
+ nabla_f[dim] = p;
+
+ nabla_f(Eigen::seq(0, dim - 1)) = xax + lambda*grad_p;
+ nabla_f[dim] = p;
+ // The exit criterion for the Newton optimization is on the magnitude of the gradient of the Lagrangian,
+ // which is zero at a solution.
+ double nabla_f_norm = nabla_f.norm();
+
+ // Do Newton iterations.
+ while((nabla_f_norm > redistOptions.tolerance) && (k_newton<redistOptions.max_iter))
+ {
+ // gather required derivative values
+ H_p = get_H_p_minter(x, model);
+
+ // Assemble Hessian matrix, grad f has been computed at the end of the last iteration.
+ H_f(Eigen::seq(0, dim_r - 1), Eigen::seq(0, dim_r - 1)) = lambda*H_p;
+ for(int k = 0; k < dim; k++) H_f(k, k) = H_f(k, k) + 1.0;
+ H_f(Eigen::seq(0, dim_r - 1), dim_r) = grad_p;
+ H_f(dim_r, Eigen::seq(0, dim_r - 1)) = grad_p.transpose();
+ H_f(dim_r, dim_r) = 0.0;
+
+ // compute Newton increment
+ dx = - H_f.inverse()*nabla_f;
+
+ // prevent Newton algorithm from leaving the support radius by scaling step size. It is invoked in case
+ // the increment exceeds 50% of the cutoff radius and simply scales the step length down to 10% until it
+ // does not exceed 50% of the cutoff radius.
+
+ while((dx.dot(dx)) > 0.25*r_cutoff2)
+ {
+ message_step_limitation = 1;
+ dx = 0.1*dx;
+ }
+
+ // apply increment and compute new iteration values
+ x = x + dx(Eigen::seq(0, dim - 1));
+ lambda = lambda + dx[dim];
+
+ // prepare values for next iteration and update the exit criterion
+ xax = x - xa;
+ p = get_p_minter(x, model);
+ grad_p = get_grad_p_minter(x, model);
+
+ nabla_f(Eigen::seq(0, dim - 1)) = xax + lambda*grad_p;
+ nabla_f[dim] = p;
+
+ //norm_dx = dx.norm(); // alternative criterion: incremental tolerance
+ nabla_f_norm = nabla_f.norm();
+
+ ++k_newton;
+
+ if(redistOptions.verbose)
+ {
+ std::cout<<"dx: "<<dx[0]<<", "<<dx[1]<<std::endl;
+ std::cout<<"H_f:\n"<<H_f<<"\nH_f_inv:\n"<<H_f.inverse()<<std::endl;
+ std::cout<<"x:"<<x[0]<<", "<<x[1]<<"\nc:"<<std::endl;
+ std::cout<<c<<std::endl;
+ std::cout<<"dpdx: "<<grad_p<<std::endl;
+ std::cout<<"k = "<<k_newton<<std::endl;
+ std::cout<<"x_k = "<<x[0]<<", "<<x[1]<<std::endl;
+ std::cout<<x[0]<<", "<<x[1]<<", "<<nabla_f_norm<<std::endl;
+ }
+ }
+
+ // Check if the Newton algorithm achieved the required accuracy within the allowed number of iterations.
+ // If not, throw a warning, but proceed as usual (even if the accuracy dictated by the user cannot be reached,
+ // the result can still be very good.
+ if (k_newton == redistOptions.max_iter) message_convergence_problem = 1;
+
+ // And finally, compute the new sdf value as the distance of the query particle x_a and the result of the
+ // optimization scheme x. We conserve the initial sign of the sdf, since the particle moves with the phase
+ // and does not cross the interface.
+ if (redistOptions.write_sdf) vd_in.template getProp<vd_in_sdf>(a) = return_sign(vd_in.template getProp<vd_in_sdf>(a))*xax.norm();
+ if (redistOptions.write_cp) for(int k = 0; k <dim; k++) vd_in.template getProp<vd_in_cp>(a)[k] = x[k];
+ // This is to avoid loss of mass conservation - in some simulations, a particle can get assigned a 0-sdf, so
+ // that it lies on the surface and cannot be distinguished from the one phase or the other. The reason for
+ // this probably lies in the numerical tolerance. As a quick fix, we introduce an error of the tolerance,
+ // but keep the sign.
+ if ((k_newton == 0) && (xax.norm() < redistOptions.tolerance))
+ {
+ vd_in.template getProp<vd_in_sdf>(a) = return_sign(vd_in.template getProp<vd_in_sdf>(a))*redistOptions.tolerance;
+ }
+
+ // debug optimisation
+ if(redistOptions.verbose)
+ {
+ //std::cout<<"new sdf: "<<vd.getProp<sdf>(a)<<std::endl;
+ std::cout<<"x_final: "<<x[0]<<", "<<x[1]<<std::endl;
+ std::cout<<"p(x_final) :"<<get_p_minter(x, model)<<std::endl;
+ std::cout<<"nabla p(x_final)"<<get_grad_p_minter(x, model)<<std::endl;
+
+ std::cout<<"lambda: "<<lambda<<std::endl;
+ }
+
+ if (redistOptions.compute_normals)
+ {
+ for(int k = 0; k<dim; k++) vd_in.template getProp<vd_in_normal>(a)[k] = return_sign(vd_in.template getProp<vd_in_sdf>(a))*grad_p(k)*1/grad_p.norm();
+ }
+
+ if (redistOptions.compute_curvatures)
+ {
+ H_p = get_H_p_minter(x, model);
+ // divergence of normalized gradient field:
+ if (dim == 2)
+ {
+ vd_in.template getProp<vd_in_curvature>(a) = (H_p(0,0)*grad_p(1)*grad_p(1) - 2*grad_p(1)*grad_p(0)*H_p(0,1) + H_p(1,1)*grad_p(0)*grad_p(0))/std::pow(sqrt(grad_p(0)*grad_p(0) + grad_p(1)*grad_p(1)),3);
+ }
+ else if (dim == 3)
+ { // Mean curvature is 0.5*fluid mechanical curvature (see https://link.springer.com/article/10.1007/s00466-021-02128-9), but here we get fluidmechanical curvature, which is the divergence of the normalized gradient field
+ vd_in.template getProp<vd_in_curvature>(a) =
+ ((H_p(1,1) + H_p(2,2))*std::pow(grad_p(0), 2) + (H_p(0,0) + H_p(2,2))*std::pow(grad_p(1), 2) + (H_p(0,0) + H_p(1,1))*std::pow(grad_p(2), 2)
+ - 2*grad_p(0)*grad_p(1)*H_p(0,1) - 2*grad_p(0)*grad_p(2)*H_p(0,2) - 2*grad_p(1)*grad_p(2)*H_p(1,2))*std::pow(std::pow(grad_p(0), 2) + std::pow(grad_p(1), 2) + std::pow(grad_p(2), 2), -1.5);
+ }
+
+ }
+ ++part;
+ }
+ if (redistOptions.verbose and message_step_limitation)
+ {
+ std::cout<<"Step size limitation invoked"<<std::endl;
+ }
+ if (message_convergence_problem)
+ {
+ std::cout<<"Warning: Newton algorithm has reached maximum number of iterations, does not converge for some particles"<<std::endl;
+ }
+
+ }
+
+ template<typename NNlist_type> vect_dist_key_dx get_closest_neighbor(Point<dim, double> & xa, particles_surface<dim, n_c> & vd_surface, NNlist_type & NN_s)
+ {
+ auto Np = NN_s.template getNNIterator<NO_CHECK>(NN_s.getCell(xa));
+
+ // distance is the the minimum distance to be beaten
+ double distance = 1000000.0;
+ // dist_calc is the variable for the distance that is computed per sample point
+ double dist_calc = 1000000000.0;
+ // b_min is the handle referring to the particle that carries the sample point with the minimum distance so far
+ vect_dist_key_dx b_min = Np.get();
+
+ while(Np.isNext())
+ {
+ vect_dist_key_dx b = Np.get();
+
+ if (!vd_s.template getProp<vd_s_close_part>(b))
+ {
+ ++Np;
+ continue;
+ }
+ Point<dim, double> xb = vd_s.template getProp<vd_s_sample>(b);
+ dist_calc = norm(xa - xb);
+ if (dist_calc < distance)
+ {
+ distance = dist_calc;
+ b_min = b;
+ }
+
+ ++Np;
+ }
+ return(b_min);
+ }
+
+ // minterface
+ template<typename PolyType>
+ inline double get_p_minter(EMatrix<double, Eigen::Dynamic, 1> xvector, PolyType model)
+ {
+ return(model->eval(xvector.transpose())(0));
+ }
+
+
+ template<typename PolyType>
+ inline EMatrix<double, Eigen::Dynamic, 1> get_grad_p_minter(EMatrix<double, Eigen::Dynamic, 1> xvector, PolyType model)
+ {
+ EMatrix<double, Eigen::Dynamic, 1> grad_p(dim, 1);
+ std::vector<int> derivOrder(dim, 0);
+ for(int k = 0; k < dim; k++){
+ std::fill(derivOrder.begin(), derivOrder.end(), 0);
+ derivOrder[k] = 1;
+ grad_p[k] = model->deriv_eval(xvector.transpose(), derivOrder)(0);
+ }
+ return(grad_p);
+ }
+
+ template<typename PolyType>
+ inline EMatrix<double, Eigen::Dynamic, Eigen::Dynamic> get_H_p_minter(EMatrix<double, Eigen::Dynamic, 1> xvector, PolyType model)
+ {
+ EMatrix<double, Eigen::Dynamic, Eigen::Dynamic> H_p(dim, dim);
+ std::vector<int> derivOrder(dim, 0);
+
+ for(int k = 0; k < dim; k++){
+ for(int l = 0; l < dim; l++)
+ {
+ std::fill(derivOrder.begin(), derivOrder.end(), 0);
+ derivOrder[k]++;
+ derivOrder[l]++;
+ H_p(k, l) = model->deriv_eval(xvector.transpose(), derivOrder)(0);
+ }
+ }
+ return(H_p);
+ }
+
+};
+
diff --git a/src/level_set/particle_cp/pcp_unit_test_helpfunctions.h b/src/level_set/particle_cp/pcp_unit_test_helpfunctions.h
new file mode 100644
index 0000000000000000000000000000000000000000..6f13bc7b694640a739393ccaeccb4eebd133c50f
--- /dev/null
+++ b/src/level_set/particle_cp/pcp_unit_test_helpfunctions.h
@@ -0,0 +1,196 @@
+#include <math.h>
+// This is a collection of helpfunctions that were used to run convergence tests and benchmarks for the ellipse/ellipsoid
+// in the particle closest point draft. Computing the theoretical closest points and distances from a given query point
+// is done using the first four functions which were adopted from David Eberly "Distance from a Point to an Ellipse, an
+// Ellipsoid, or a Hyperellipsoid", 2013.
+//
+// Created by lschulze
+double GetRoot (double r0, double z0, double z1, double g)
+ {
+ const int maxIter = 100;
+ double n0 = r0*z0;
+ double s0 = z1 - 1;
+ double s1 = ( g < 0 ? 0 : sqrt(n0*n0+z1*z1) - 1 ) ;
+ double s = 0;
+ for ( int i = 0; i < maxIter; ++i ){
+ if (i == (maxIter - 1)) std::cout<<"distance point ellipse algorithm did not converge."<<std::endl;
+ s = ( s0 + s1 ) / 2 ;
+ if ( s == s0 || s == s1 ) {break; }
+ double ratio0 = n0 /( s + r0 );
+ double ratio1 = z1 /( s + 1 );
+ g = ratio0*ratio0 + ratio1*ratio1 - 1 ;
+ if (g > 0) {s0 = s;} else if (g < 0) {s1 = s ;} else {break ;}
+ }
+ return s;
+ }
+
+double GetRoot(double r0, double r1, double z0, double z1, double z2, double g)
+{
+ const int maxIter = 100;
+ double n0 = r0*z0;
+ double n1 = r1*z1;
+ double s0 = z2 - 1;
+ double s1 = (g < 0 ? 0 : sqrt(n0*n0 + n1*n1 + z2*z2) - 1) ;
+ double s = s;
+ for(int i = 0 ; i < maxIter ; ++i )
+ {
+ if (i == (maxIter - 1)) std::cout<<"distance point ellipse algorithm did not converge."<<std::endl;
+ s = ( s0 + s1 ) / 2 ;
+ if ( s == s0 || s == s1 ) {break; }
+ double ratio0 = n0 / ( s + r0 );
+ double ratio1 = n1 / ( s + r1 );
+ double ratio2 = z2 / ( s + 1 );
+ g = ratio0*ratio0 + ratio1*ratio1 +ratio2*ratio2 - 1;
+ if ( g > 0 ) { s0 = s ;}
+ else if ( g < 0 ) { s1 = s ; }
+ else {break;}
+ }
+ return (s);
+}
+
+double DistancePointEllipse(double e0, double e1, double y0, double y1, double& x0, double& x1)
+ {
+ double distance;
+ if ( y1 > 0){
+ if ( y0 > 0){
+ double z0 = y0 / e0;
+ double z1 = y1 / e1;
+ double g = z0*z0+z1*z1 - 1;
+ if ( g != 0){
+ double r0 = (e0/e1)*(e0/e1);
+ double sbar = GetRoot(r0 , z0 , z1 , g);
+ x0 = r0 * y0 /( sbar + r0 );
+ x1 = y1 /( sbar + 1 );
+ distance = sqrt( (x0-y0)*(x0-y0) + (x1-y1)*(x1-y1) );
+ }else{
+ x0 = y0;
+ x1 = y1;
+ distance = 0;
+ }
+ }
+ else // y0 == 0
+ {x0 = 0 ; x1 = e1 ; distance = abs( y1 - e1 );}
+ }else{ // y1 == 0
+ double numer0 = e0*y0 , denom0 = e0*e0 - e1*e1;
+ if ( numer0 < denom0 ){
+ double xde0 = numer0/denom0;
+ x0 = e0*xde0 ; x1 = e1*sqrt(1 - xde0*xde0 );
+ distance = sqrt( (x0-y0)*(x0-y0) + x1*x1 );
+ }else{
+ x0 = e0;
+ x1 = 0;
+ distance = abs( y0 - e0 );
+ }
+ }
+ return distance;
+ }
+
+double DistancePointEllipsoid(double e0, double e1, double e2, double y0, double y1, double y2, double& x0, double& x1, double& x2)
+{
+ double distance;
+ if( y2 > 0 )
+ {
+ if( y1 > 0 )
+ {
+ if( y0 > 0 )
+ {
+ double z0 = y0 / e0;
+ double z1 = y1 / e1;
+ double z2 = y2 / e2;
+ double g = z0*z0 + z1*z1 + z2*z2 - 1 ;
+ if( g != 0 )
+ {
+ double r0 = (e0/e2)*(e0/e2);
+ double r1 = (e1/e2)*(e1/e2);
+ double sbar = GetRoot ( r0 , r1 , z0 , z1 , z2 , g );
+ x0 = r0 *y0 / ( sbar + r0 );
+ x1 = r1 *y1 / ( sbar + r1 );
+ x2 = y2 / ( sbar + 1 );
+ distance = sqrt( (x0 - y0)*(x0 - y0) + (x1 - y1)*(x1 - y1) + (x2 - y2)*(x2 - y2));
+ }
+ else
+ {
+ x0 = y0;
+ x1 = y1;
+ x2 = y2;
+ distance = 0;
+ }
+ }
+ else // y0 == 0
+ {
+ x0 = 0;
+ distance = DistancePointEllipse( e1 , e2 , y1 , y2, x1, x2);
+ }
+ }
+ else // y1 == 0
+ {
+ if( y0 > 0 )
+ {
+ x1 = 0;
+ distance = DistancePointEllipse( e0 , e2 , y0 , y2, x0, x2);
+ }
+ else // y0 == 0
+ {
+ x0 = 0;
+ x1 = 0;
+ x2 = e2;
+ distance = abs(y2 - e2);
+ }
+ }
+ }
+ else // y2 == 0
+ {
+ double denom0 = e0*e0 - e2*e2;
+ double denom1 = e1*e1 - e2*e2;
+ double numer0 = e0*y0;
+ double numer1 = e1*y1;
+ bool computed = false;
+ if((numer0 < denom0) && (numer1 < denom1))
+ {
+ double xde0 = numer0/denom0;
+ double xde1 = numer1/denom1 ;
+ double xde0sqr = xde0 *xde0;
+ double xde1sqr = xde1 * xde1 ;
+ double discr = 1 - xde0sqr - xde1sqr;
+ if( discr > 0 )
+ {
+ x0 = e0*xde0;
+ x1 = e1*xde1;
+ x2 = e2*sqrt(discr);
+ distance = sqrt((x0 - y0)*(x0 - y0) + (x1 - y1)*(x1 - y1) + x2*x2);
+ computed = true;
+ }
+ }
+ if( !computed )
+ {
+ x2 = 0;
+ distance = DistancePointEllipse(e0 , e1 , y0 , y1, x0, x1);
+ }
+ }
+ return distance;
+}
+
+constexpr unsigned int factorial(unsigned int x)
+{
+ unsigned int fact = 1;
+ for(int i = 1; i < (x + 1); i++) fact = fact*i;
+ return(fact);
+}
+constexpr unsigned int minter_lp_degree_one_num_coeffs(unsigned int dims, unsigned int poly_degree)
+{
+ return(factorial(dims + poly_degree)/(factorial(dims)*factorial(poly_degree)));
+}
+
+int return_sign(double phi)
+{
+ if (phi > 0) return 1;
+ if (phi < 0) return -1;
+ return 0;
+}
+
+double randMinusOneToOne()
+{
+ double temp = rand() / (RAND_MAX + 1.0);
+ //std::cout<<(2.0*temp - 1.0)<<std::endl;
+ return(2.0*temp - 1.0);
+}
diff --git a/src/level_set/particle_cp/pcp_unit_tests.cpp b/src/level_set/particle_cp/pcp_unit_tests.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..28ece4c5b2ea4b953fcb322e8446e4b4840dfc95
--- /dev/null
+++ b/src/level_set/particle_cp/pcp_unit_tests.cpp
@@ -0,0 +1,198 @@
+/* Unit tests for the particle closest point method
+ * Date : 29 June 2023
+ * Author : lschulze
+ */
+
+#include<iostream>
+#include <boost/test/unit_test_log.hpp>
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+#include <math.h>
+//#include <sys/_types/_size_t.h>
+#include "Vector/vector_dist.hpp"
+#include "DCPSE/Dcpse.hpp"
+#include "DCPSE/MonomialBasis.hpp"
+#include "Draw/DrawParticles.hpp"
+#include "particle_cp.hpp"
+#include "Operators/Vector/vector_dist_operators.hpp"
+#include <chrono>
+#include "pcp_unit_test_helpfunctions.h"
+
+typedef struct EllipseParameters{
+ double origin[3];
+ double radiusA;
+ double radiusB;
+ double radiusC;
+ double eccentricity;
+} EllipseParams;
+
+// Generate an ellipsoid initial levelset signed distance function
+template<typename particles_type, typename iterator_type, size_t sdf, size_t ref_cp>
+void initializeLSEllipsoid(particles_type &vd, iterator_type particle_it, const EllipseParams ¶ms, double bandwidth, double perturb_factor, double H)
+{
+ while(particle_it.isNext())
+ {
+ double posx = particle_it.get().get(0);
+ double posy = particle_it.get().get(1);
+ double posz = particle_it.get().get(2);
+
+ double ref_cpx = 0.0;
+ double ref_cpy = 0.0;
+ double ref_cpz = 0.0;
+
+ double dist = DistancePointEllipsoid(params.radiusA, params.radiusB, params.radiusC, abs(posx), abs(posy), abs(posz), ref_cpx, ref_cpy, ref_cpz);
+ if (abs(dist) < bandwidth/2.0)
+ {
+ posx = posx + perturb_factor*randMinusOneToOne()*H;
+ posy = posy + perturb_factor*randMinusOneToOne()*H;
+ posz = posz + perturb_factor*randMinusOneToOne()*H;
+
+ dist = DistancePointEllipsoid(params.radiusA, params.radiusB, params.radiusC, abs(posx), abs(posy), abs(posz), ref_cpx, ref_cpy, ref_cpz);
+ vd.add();
+ vd.getLastPos()[0] = posx;
+ vd.getLastPos()[1] = posy;
+ vd.getLastPos()[2] = posz;
+
+ double phi_val = sqrt(((posx - params.origin[0])/params.radiusA)*((posx - params.origin[0])/params.radiusA) + ((posy - params.origin[1])/params.radiusB)*((posy - params.origin[1])/params.radiusB) + ((posz - params.origin[2])/params.radiusC)*((posz - params.origin[2])/params.radiusC)) - 1.0;
+ vd.template getLastProp<sdf>() = phi_val;
+ vd.template getLastProp<ref_cp>()[0] = return_sign(posx)*ref_cpx;
+ vd.template getLastProp<ref_cp>()[1] = return_sign(posy)*ref_cpy;
+ vd.template getLastProp<ref_cp>()[2] = return_sign(posz)*ref_cpz;
+ }
+
+ ++particle_it;
+ }
+}
+
+BOOST_AUTO_TEST_SUITE( particle_closest_point_test )
+
+BOOST_AUTO_TEST_CASE( ellipsoid )
+{
+ // simulation params
+ constexpr int poly_order = 4;
+ const double H = 1.0/64.0;
+ const double perturb_factor = 0.3;
+
+ // pcp params
+ const double bandwidth = 12.0*H;
+ const double regression_rcut_factor = 2.4;
+ const double threshold = 1e-13;
+
+ // initialize domain and particles
+ const double l = 2.0;
+ Box<3, double> domain({-l/2.0, -l/3.0, -l/3.0}, {l/2.0, l/3.0, l/3.0});
+ size_t sz[3] = {(size_t)(l/H + 0.5), (size_t)((2.0/3.0)*l/H + 0.5), (size_t)((2.0/3.0)*l/H + 0.5)};
+ size_t bc[3] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
+ Ghost<3, double> g(bandwidth);
+
+ constexpr int sdf = 0;
+ constexpr int cp = 1;
+ constexpr int normal = 2;
+ constexpr int curvature = 3;
+ constexpr int ref_cp = 4;
+ typedef vector_dist<3, double, aggregate<double, Point<3, double>, Point<3, double>, double, Point<3, double>>> particles;
+ // | | | | |
+ // SDF closest point | curvature |
+ // surface normal reference closest point
+
+ particles vd(0, domain, bc, g, DEC_GRAN(512));
+ openfpm::vector<std::string> names({"sdf","cp","normal","curvature","ref_cp"});
+ vd.setPropNames(names);
+ EllipseParameters params;
+ for (int k = 0; k < 3; k++) params.origin[k] = 0.0;
+ params.radiusA = 0.75;
+ params.radiusB = 0.5;
+ params.radiusC = 0.5;
+
+ auto particle_it = DrawParticles::DrawBox(vd, sz, domain, domain);
+
+ // initialize spurious sdf values and reference solution
+ initializeLSEllipsoid<particles, decltype(particle_it), sdf, ref_cp>(vd, particle_it, params, bandwidth, perturb_factor, H);
+
+ //vd.write("pcpunittest_init");
+ // initialize and run pcp redistancing
+ Redist_options rdistoptions;
+ rdistoptions.minter_poly_degree = poly_order;
+ rdistoptions.H = H;
+ rdistoptions.r_cutoff_factor = regression_rcut_factor;
+ rdistoptions.sampling_radius = 0.75*bandwidth;
+ rdistoptions.tolerance = threshold;
+ rdistoptions.write_cp = 1;
+ rdistoptions.compute_normals = 1;
+ rdistoptions.compute_curvatures = 1;
+ rdistoptions.only_narrowband = 0;
+
+ static constexpr unsigned int num_coeffs = minter_lp_degree_one_num_coeffs(3, poly_order);
+
+ particle_cp_redistancing<particles, sdf, cp, normal, curvature, num_coeffs> pcprdist(vd, rdistoptions);
+ std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
+ pcprdist.run_redistancing();
+ std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
+ std::cout << "Time difference for pcp redistancing = " << std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count() << "[ms]" << std::endl;
+
+ //vd.write("pcpunittest");
+ // iterate through particles and compute error
+ auto part = vd.getDomainIterator();
+ double sdferr = 0.0;
+ double sdfmaxerr = 0.0;
+ double sdfmeanerr = 0.0;
+ double cperr = 0.0;
+ double cpmaxerr = 0.0;
+ double normalerr = 0.0;
+ double normalmaxerr = 0.0;
+ double curvatureerr = 0.0;
+ double curvaturemaxerr = 0.0;
+ int particlecount = 0;
+
+ while(part.isNext())
+ {
+ auto a = part.get();
+ Point<3, double> pos = vd.getPos(a);
+ double reference_sdf = return_sign(sqrt(((pos[0] - params.origin[0])/params.radiusA)*((pos[0] - params.origin[0])/params.radiusA) + ((pos[1] - params.origin[1])/params.radiusB)*((pos[1] - params.origin[1])/params.radiusB) + ((pos[2] - params.origin[2])/params.radiusC)*((pos[2] - params.origin[2])/params.radiusC)) - 1.0)*norm(vd.getProp<ref_cp>(a) - pos);
+ Point<3, double> reference_normal;
+ double reference_curvature;
+
+ // compute reference SDF value w/ reference closest point and check computed SDF value
+ sdferr = abs(vd.getProp<sdf>(a) - reference_sdf);
+ if (sdferr > sdfmaxerr) sdfmaxerr = sdferr;
+ sdfmeanerr += sdferr;
+ ++particlecount;
+
+ // check computed closest point against reference closest point
+ cperr = norm(vd.getProp<ref_cp>(a) - vd.getProp<cp>(a));
+ if (cperr > cpmaxerr) cpmaxerr = cperr;
+
+ // compute reference normal and check computed normal
+ reference_normal[0] = 2*vd.getProp<ref_cp>(a)[0]/(params.radiusA*params.radiusA);
+ reference_normal[1] = 2*vd.getProp<ref_cp>(a)[1]/(params.radiusB*params.radiusB);
+ reference_normal[2] = 2*vd.getProp<ref_cp>(a)[2]/(params.radiusC*params.radiusC);
+ reference_normal = return_sign(reference_sdf)*reference_normal/norm(reference_normal);
+ normalerr = norm(reference_normal - vd.getProp<normal>(a));
+ if (normalerr > normalmaxerr) normalmaxerr = normalerr;
+
+ // compute reference curvature and check computed curvature
+ reference_curvature = (std::abs(vd.getProp<ref_cp>(a)[0]*vd.getProp<ref_cp>(a)[0] + vd.getProp<ref_cp>(a)[1]*vd.getProp<ref_cp>(a)[1] + vd.getProp<ref_cp>(a)[2]*vd.getProp<ref_cp>(a)[2] - params.radiusA*params.radiusA - params.radiusB*params.radiusB - params.radiusC*params.radiusC))/(params.radiusA*params.radiusA*params.radiusB*params.radiusB*params.radiusC*params.radiusC*std::pow(vd.getProp<ref_cp>(a)[0]*vd.getProp<ref_cp>(a)[0]/std::pow(params.radiusA, 4) + vd.getProp<ref_cp>(a)[1]*vd.getProp<ref_cp>(a)[1]/std::pow(params.radiusB, 4) + vd.getProp<ref_cp>(a)[2]*vd.getProp<ref_cp>(a)[2]/std::pow(params.radiusC, 4), 1.5));
+ curvatureerr = abs(reference_curvature - vd.getProp<curvature>(a));
+ if (curvatureerr > curvaturemaxerr) curvaturemaxerr = curvatureerr;
+
+ ++part;
+ }
+ sdfmeanerr = sdfmeanerr/particlecount;
+ std::cout<<"Mean error for sdf is: "<<sdfmeanerr<<std::endl;
+ std::cout<<"Maximum error for sdf is: "<<sdfmaxerr<<std::endl;
+ std::cout<<"Maximum error for the closest point is: "<<cpmaxerr<<std::endl;
+ std::cout<<"Maximum error for surface normal is: "<<normalmaxerr<<std::endl;
+ std::cout<<"Maximum error for curvature is: "<<curvaturemaxerr<<std::endl;
+
+ double tolerance = 1e-7;
+ bool check;
+ if (std::abs(sdfmaxerr) < tolerance)
+ check = true;
+ else
+ check = false;
+
+ BOOST_TEST( check );
+
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/regression/poly_levelset.hpp b/src/regression/poly_levelset.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..30581e8004f9b6faf38d72fd40760b0f4330e523
--- /dev/null
+++ b/src/regression/poly_levelset.hpp
@@ -0,0 +1,118 @@
+/*
+ * Regression module
+ * Obtains polynomial models for data from vector_dist
+ * author : sachin (sthekke@mpi-cbg.de)
+ * date : 18.01.2023
+ *
+ */
+#ifndef OPENFPM_NUMERICS_POLYLEVELSET_HPP
+#define OPENFPM_NUMERICS_POLYLEVELSET_HPP
+
+#include "Vector/map_vector.hpp"
+#include "Space/Shape/Point.hpp"
+#include "DMatrix/EMatrix.hpp"
+
+#include "minter.h"
+
+
+template<int spatial_dim, typename MatType = EMatrixXd, typename VecType = EVectorXd>
+class PolyLevelset
+{
+ minter::LevelsetPoly<spatial_dim, MatType, VecType> *model;
+
+public:
+ template<typename vector_type>
+ PolyLevelset(vector_type &vd, double tol)
+ {
+ constexpr int dim = vector_type::dims;
+
+ MatType points(vd.size_local(), dim);
+
+ auto it = vd.getDomainIterator();
+ int i = 0;
+ while(it.isNext())
+ {
+ auto key = it.get();
+ for(int j = 0;j < dim;++j)
+ points(i,j) = vd.getPos(key)[j];
+
+ ++i;
+ ++it;
+ }
+
+ // construct polynomial model
+ model = new minter::LevelsetPoly<spatial_dim, MatType, VecType>(points, tol);
+ }
+
+
+ ~PolyLevelset()
+ {
+ if(model)
+ delete model;
+ }
+
+ // T : Point<vector_type::dims, typename vector_type::stype>
+ template<typename T>
+ double eval(T pos)
+ {
+ int dim = pos.dims;
+ MatType point(1,dim);
+ for(int j = 0;j < dim;++j)
+ point(0,j) = pos.get(j);
+
+ return model->eval(point)(0);
+ }
+
+ // T1 : Point<vector_type::dims, typename vector_type::stype>
+ // T2 : Point<vector_type::dims, int>
+ template<typename T1, typename T2>
+ double deriv(T1 pos, T2 deriv_order)
+ {
+ int dim = pos.dims;
+ MatType point(1,dim);
+ for(int j = 0;j < dim;++j)
+ point(0,j) = pos.get(j);
+
+ std::vector<int> order;
+ for(int j = 0;j < dim;++j)
+ order.push_back(deriv_order.get(j));
+
+ return model->deriv_eval(point, order)(0);
+ }
+
+ // T : Point<vector_type::dims, typename vector_type::stype>
+ template<typename T>
+ T estimate_normals_at(T pos)
+ {
+ int dim = pos.dims;
+ MatType point(1,dim);
+ for(int j = 0;j < dim;++j)
+ point(0,j) = pos.get(j);
+
+ T normal;
+ auto normal_minter = model->estimate_normals_at(point);
+
+ for(int j = 0;j < dim;++j)
+ normal.get(j) = normal_minter(0,j);
+
+ return normal;
+ }
+
+ // T : Point<vector_type::dims, typename vector_type::stype>
+ template<typename T>
+ double estimate_mean_curvature_at(T pos)
+ {
+ int dim = pos.dims;
+ MatType point(1,dim);
+ for(int j = 0;j < dim;++j)
+ point(0,j) = pos.get(j);
+
+ auto mc = model->estimate_mean_curvature_at(point);
+
+ return mc(0);
+ }
+};
+
+
+
+#endif /* POLYLEVELSET_HPP_ */
\ No newline at end of file
diff --git a/src/regression/poly_levelset_test.cpp b/src/regression/poly_levelset_test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e2cc89a9bced067424f575595d6776f6adb983f2
--- /dev/null
+++ b/src/regression/poly_levelset_test.cpp
@@ -0,0 +1,92 @@
+/*
+ * Unit tests for the Regression module PolyLevelset submodule
+ * author : Sachin (sthekke@mpi-cbg.de)
+ * date : 18.01.2023
+ *
+ */
+
+
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+#include "poly_levelset.hpp"
+#include "Vector/vector_dist.hpp"
+#include "DMatrix/EMatrix.hpp"
+
+BOOST_AUTO_TEST_SUITE( Regression_test )
+
+
+
+BOOST_AUTO_TEST_CASE ( PolyLevelset_Sphere )
+{
+ Box<3,double> domain({-2.0,-2.0,-2.0},{2.0,2.0,2.0});
+
+ // Here we define the boundary conditions of our problem
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+ // extended boundary around the domain, and the processor domain
+ Ghost<3,double> g(0.01);
+
+ using vectorType = vector_dist<3,double, aggregate<double, double> >;
+
+ vectorType vd(1024,domain,bc,g);
+
+ constexpr int mean_curvature = 0;
+ constexpr int gauss_curvature = 1;
+
+ // Initialize points on sphere
+ auto it = vd.getDomainIterator();
+ while (it.isNext())
+ {
+ auto key = it.get();
+ double theta = ((double)rand() / RAND_MAX) * M_PI;
+ double phi = ((double)rand() / RAND_MAX) * 2.0 * M_PI;
+
+ vd.getPos(key)[0] = cos(theta) * sin(phi);
+ vd.getPos(key)[1] = cos(theta) * cos(phi);
+ vd.getPos(key)[2] = sin(theta);
+
+ vd.template getProp<mean_curvature>(key) = 0.0;
+ vd.template getProp<gauss_curvature>(key) = 0.0;
+
+ ++it;
+ }
+ vd.map();
+
+ // auto model = PolyLevelset<3>(vd, 1e-4);
+/*
+ double max_err = -1.0;
+ auto it2 = vd.getDomainIterator();
+ while (it2.isNext())
+ {
+ auto key = it2.get();
+ Point<3, double> pos = {vd.getPos(key)[0], vd.getPos(key)[1], vd.getPos(key)[2]};
+ // vd.template getProp<mean_curvature>(key) = model.estimate_mean_curvature_at(pos);
+ // vd.template getProp<gauss_curvature>(key) = model->estimate_gauss_curvature_at(vd.getPos(key));
+
+ double val = vd.getProp<mean_curvature>(key);
+ double actual = 1.0;
+ double err = std::abs(actual - val);
+ if (err > max_err) max_err = err;
+
+ ++it2;
+ }
+
+
+ double tolerance = 1e-4;
+ bool check;
+ if (std::abs(max_err) < tolerance)
+ check = true;
+ else
+ check = false;
+ std::cout<<"Max err (poly level) = "<<max_err<<"\n";
+
+ BOOST_TEST( check );
+ */
+ // if(model)
+ // delete model;
+
+}
+
+
+BOOST_AUTO_TEST_SUITE_END()
+
diff --git a/src/regression/regression.hpp b/src/regression/regression.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..bb52ba209512e58188e77467ffef67f7f5760ec3
--- /dev/null
+++ b/src/regression/regression.hpp
@@ -0,0 +1,460 @@
+/*
+ * Regression module
+ * Obtains polynomial models for data from vector_dist
+ * author : sachin (sthekke@mpi-cbg.de)
+ * date : 28.04.2022
+ *
+ */
+#ifndef OPENFPM_NUMERICS_REGRESSION_HPP
+#define OPENFPM_NUMERICS_REGRESSION_HPP
+
+#include "Vector/map_vector.hpp"
+#include "Space/Shape/Point.hpp"
+#include "DMatrix/EMatrix.hpp"
+#include "DCPSE/SupportBuilder.hpp"
+#include "minter.h"
+
+
+template<typename vector_type_support, typename NN_type>
+class RegressionSupport
+{
+ openfpm::vector<size_t> keys;
+
+public:
+
+ template<typename iterator_type>
+ RegressionSupport(vector_type_support &vd, iterator_type itPoint, unsigned int requiredSize, support_options opt, NN_type &cellList_in) : domain(vd),
+ cellList(cellList_in)
+ {
+ rCut = cellList_in.getCellBox().getHigh(0);
+ // Get spatial position from point iterator
+ vect_dist_key_dx p = itPoint.get();
+ vect_dist_key_dx pOrig = itPoint.getOrig();
+ Point<vector_type_support::dims, typename vector_type_support::stype> pos = domain.getPos(p.getKey());
+
+ // Get cell containing current point and add it to the set of cell keys
+ grid_key_dx<vector_type_support::dims> curCellKey = cellList.getCellGrid(pos); // Here get the key of the cell where the current point is
+ std::set<grid_key_dx<vector_type_support::dims>> supportCells;
+ supportCells.insert(curCellKey);
+
+ // Make sure to consider a set of cells providing enough points for the support
+ enlargeSetOfCellsUntilSize(supportCells, requiredSize + 1,
+ opt); // NOTE: this +1 is because we then remove the point itself
+
+ // Now return all the points from the support into a vector
+ keys = getPointsInSetOfCells(supportCells, p, pOrig, requiredSize, opt);
+ }
+
+ auto getKeys()
+ {
+ return keys;
+ }
+
+ auto getNumParticles()
+ {
+ return keys.size();
+ }
+
+private:
+
+ vector_type_support &domain;
+ NN_type &cellList;
+ typename vector_type_support::stype rCut;
+
+
+ void enlargeSetOfCellsUntilSize(std::set<grid_key_dx<vector_type_support::dims>> &set, unsigned int requiredSize,
+ support_options opt) {
+ if (opt == support_options::RADIUS) {
+ auto cell = *set.begin();
+ grid_key_dx<vector_type_support::dims> middle;
+ int n = std::ceil(rCut / cellList.getCellBox().getHigh(0));
+ size_t sz[vector_type_support::dims];
+ for (int i = 0; i < vector_type_support::dims; i++) {
+ sz[i] = 2 * n + 1;
+ middle.set_d(i, n);
+ }
+ grid_sm<vector_type_support::dims, void> g(sz);
+ grid_key_dx_iterator<vector_type_support::dims> g_k(g);
+ while (g_k.isNext()) {
+ auto key = g_k.get();
+ key = cell + key - middle;
+ if (isCellKeyInBounds(key)) {
+ set.insert(key);
+ }
+ ++g_k;
+ }
+ }
+ else if (opt == support_options::AT_LEAST_N_PARTICLES) {
+
+ int done = 0;
+ int n = 1;
+ auto cell = *set.begin();
+ grid_key_dx<vector_type_support::dims> middle;
+ size_t sz[vector_type_support::dims];
+
+ while(true) // loop for the number of cells enlarged per dimension
+ {
+ std::set<grid_key_dx<vector_type_support::dims>> temp_set;
+ for (int i = 0; i < vector_type_support::dims; i++) {
+ sz[i] = 2 * n + 1;
+ middle.set_d(i, n);
+ }
+
+ grid_sm<vector_type_support::dims, void> g(sz);
+ grid_key_dx_iterator<vector_type_support::dims> g_k(g);
+ while (g_k.isNext()) {
+ auto key = g_k.get();
+ key = cell + key - middle;
+ if (isCellKeyInBounds(key)) {
+ temp_set.insert(key);
+ }
+ ++g_k;
+ }
+ if (getNumElementsInSetOfCells(temp_set) < requiredSize) n++;
+ else
+ {
+ set = temp_set;
+ break;
+ }
+ }
+ }
+ else {
+ while (getNumElementsInSetOfCells(set) <
+ 5.0 * requiredSize) //Why 5*requiredSize? Becasue it can help with adaptive resolutions.
+ {
+ auto tmpSet = set;
+ for (const auto el: tmpSet) {
+ for (unsigned int i = 0; i < vector_type_support::dims; ++i) {
+ const auto pOneEl = el.move(i, +1);
+ const auto mOneEl = el.move(i, -1);
+ if (isCellKeyInBounds(pOneEl)) {
+ set.insert(pOneEl);
+ }
+ if (isCellKeyInBounds(mOneEl)) {
+ set.insert(mOneEl);
+ }
+ }
+ }
+
+ }
+ }
+ }
+
+ openfpm::vector<size_t> getPointsInSetOfCells(std::set<grid_key_dx<vector_type_support::dims>> set,
+ vect_dist_key_dx &p,
+ vect_dist_key_dx &pOrig,
+ size_t requiredSupportSize,
+ support_options opt) {
+ struct reord {
+ typename vector_type_support::stype dist;
+ size_t offset;
+
+ bool operator<(const reord &p) const { return this->dist < p.dist; }
+ };
+
+ openfpm::vector<reord> rp;
+ openfpm::vector<size_t> points;
+ Point<vector_type_support::dims, typename vector_type_support::stype> xp = domain.getPos(p);
+ for (const auto cellKey: set) {
+ const size_t cellLinId = getCellLinId(cellKey);
+ const size_t elemsInCell = getNumElementsInCell(cellKey);
+ for (size_t k = 0; k < elemsInCell; ++k) {
+ size_t el = cellList.get(cellLinId, k);
+
+ Point<vector_type_support::dims, typename vector_type_support::stype> xq = domain.getPosOrig(el);
+ //points.push_back(el);
+
+ reord pr;
+
+ pr.dist = xp.distance(xq);
+ pr.offset = el;
+ rp.add(pr);
+ }
+ }
+
+ if (opt == support_options::RADIUS) {
+ for (int i = 0; i < rp.size(); i++) {
+ if (rp.get(i).dist < rCut) {
+ points.add(rp.get(i).offset);
+ }
+ }
+ /* #ifdef SE_CLASS1
+ if (points.size()<requiredSupportSize)
+ {
+ std::cerr<<__FILE__<<":"<<__LINE__<<"Note that the DCPSE neighbourhood doesn't have asked no. particles (Increase the rCut or reduce the over_sampling factor)";
+ std::cout<<"Particels asked (minimum*oversampling_factor): "<<requiredSupportSize<<". Particles Possible with given options:"<<points.size()<<"."<<std::endl;
+ }
+ #endif*/
+ }
+ else if (opt == support_options::AT_LEAST_N_PARTICLES) {
+ for (int i = 0; i < rp.size(); i++) points.add(rp.get(i).offset);
+ }
+ else {
+ rp.sort();
+ for (int i = 0; i < requiredSupportSize; i++) {
+ points.add(rp.get(i).offset);
+ }
+ }
+ //MinSpacing=MinSpacing/requiredSupportSize
+ return points;
+ }
+
+ size_t getCellLinId(const grid_key_dx<vector_type_support::dims> &cellKey) {
+ mem_id id = cellList.getGrid().LinId(cellKey);
+ return static_cast<size_t>(id);
+ }
+
+ size_t getNumElementsInCell(const grid_key_dx<vector_type_support::dims> &cellKey) {
+ const size_t curCellId = getCellLinId(cellKey);
+ size_t numElements = cellList.getNelements(curCellId);
+ return numElements;
+ }
+ size_t getNumElementsInSetOfCells(const std::set<grid_key_dx<vector_type_support::dims>> &set)
+ {
+ size_t tot = 0;
+ for (const auto cell : set)
+ {
+ tot += getNumElementsInCell(cell);
+ }
+ return tot;
+}
+
+ bool isCellKeyInBounds(grid_key_dx<vector_type_support::dims> key)
+ {
+ const size_t *cellGridSize = cellList.getGrid().getSize();
+ for (size_t i = 0; i < vector_type_support::dims; ++i)
+ {
+ if (key.value(i) < 0 || key.value(i) >= cellGridSize[i])
+ {
+ return false;
+ }
+ }
+ return true;
+ }
+};
+
+
+template<int spatial_dim, unsigned int prp_id, typename MatType = EMatrixXd, typename VecType = EVectorXd>
+class RegressionModel
+{
+
+public:
+ minter::PolyModel<spatial_dim, MatType, VecType> *model = nullptr;
+ minter::PolyModel<spatial_dim, MatType, VecType> *deriv_model[spatial_dim];
+
+ RegressionModel(unsigned int poly_degree, float lp_degree) {
+ // construct polynomial model
+ model = new minter::PolyModel<spatial_dim, MatType, VecType>();
+ model->setDegree(poly_degree, lp_degree);
+
+ // placeholder for derivatives
+ for(int i = 0;i < spatial_dim;++i)
+ deriv_model[i] = nullptr;
+ }
+
+ template<typename vector_type, typename reg_support_type>
+ RegressionModel(vector_type &vd, reg_support_type &support, unsigned int poly_degree, float lp_degree = 2.0)
+ {
+ int num_particles = support->getNumParticles();
+ int dim = vector_type::dims;
+
+ MatType points(num_particles, dim);
+ VecType values(num_particles);
+
+ auto keys = support->getKeys();
+ for(int i = 0;i < num_particles;++i)
+ {
+ for(int j = 0;j < dim;++j)
+ points(i,j) = vd.getPos(keys.get(i))[j];
+ values(i) = vd.template getProp<prp_id>(keys.get(i));
+ }
+
+ // construct polynomial model
+ model = new minter::PolyModel<spatial_dim, MatType, VecType>(points, values, poly_degree, lp_degree);
+
+ // placeholder for derivatives
+ for(int i = 0;i < dim;++i)
+ deriv_model[i] = nullptr;
+ }
+
+
+ // Constructor for all points in a proc (domain + ghost) and a specified poly_degree
+ template<typename vector_type>
+ RegressionModel(vector_type &vd, unsigned int poly_degree, float lp_degree = 2.0)
+ {
+ int num_particles = vd.size_local_with_ghost();
+ int dim = vector_type::dims;
+
+ MatType points(num_particles, dim);
+ VecType values(num_particles);
+
+ auto it = vd.getDomainAndGhostIterator();
+ int i = 0;
+ while (it.isNext())
+ {
+ auto key = it.get();
+ for(int j = 0;j < dim;++j)
+ points(i,j) = vd.getPos(key)[j];
+
+ values(i) = vd.template getProp<prp_id>(key);
+
+ ++it;
+ ++i;
+ }
+
+ // construct polynomial model
+ model = new minter::PolyModel<spatial_dim, MatType, VecType>(points, values, poly_degree, lp_degree);
+
+ for(i = 0;i < dim;++i)
+ deriv_model[i] = nullptr;
+ }
+
+ // Constructor for all points in a proc (domain + ghost) within a tolerance
+ template<typename vector_type>
+ RegressionModel(vector_type &vd, double tolerance)
+ {
+ int num_particles = vd.size_local_with_ghost();
+ int dim = vector_type::dims;
+
+ MatType points(num_particles, dim);
+ VecType values(num_particles);
+
+ auto it = vd.getDomainAndGhostIterator();
+ int i = 0;
+ while (it.isNext())
+ {
+ auto key = it.get();
+ for(int j = 0;j < dim;++j)
+ points(i,j) = vd.getPos(key)[j];
+
+ values(i) = vd.template getProp<prp_id>(key);
+
+ ++it;
+ ++i;
+ }
+
+ int poly_degree = 1;
+ double error = -1.0;
+ minter::PolyModel<spatial_dim, MatType, VecType> *mdl = nullptr;
+
+ do
+ {
+ ++poly_degree;
+ if(mdl)
+ delete mdl;
+
+ // construct polynomial model
+ mdl = new minter::PolyModel<spatial_dim, MatType, VecType>(points, values, poly_degree, 2.0);
+
+ // check if linf_error is within the tolerance
+ error = -1.0;
+ for(i = 0;i < num_particles;++i)
+ {
+ double pred = mdl->eval(points.block(i,0,1,dim))(0); // evaluated for one point
+ double err = std::abs(pred - values(i));
+ if (err > error)
+ error = err;
+ }
+ // std::cout<<"Fit of degree "<<poly_degree<<" with error = "<<error<<std::endl;
+
+ }while(error > tolerance);
+
+ model = mdl;
+
+ for(i = 0;i < dim;++i)
+ deriv_model[i] = nullptr;
+
+ }
+
+ template<typename vector_type, typename reg_support_type>
+ void computeCoeffs(vector_type& vd, reg_support_type& support)
+ {
+
+ unsigned int dim = vector_type::dims;
+ auto num_particles = support.getNumParticles();
+
+ Eigen::MatrixXd points(num_particles, dim);
+ Eigen::VectorXd values(num_particles);
+
+ auto keys = support.getKeys();
+ for(int i = 0;i < num_particles;++i)
+ {
+ for(int j = 0;j < dim;++j)
+ points(i,j) = vd.getPos(keys.get(i))[j];
+ values(i) = vd.template getProp<prp_id>(keys.get(i));
+ }
+
+ model->computeCoeffs(points, values);
+ }
+
+ ~RegressionModel()
+ {
+
+ if(model)
+ delete model;
+
+ for(int i = 0;i < spatial_dim;++i)
+ {
+ if(deriv_model[i])
+ delete deriv_model[i];
+ }
+ }
+
+ template<typename T> // Typical: Point<vector_type::dims, typename vector_type::stype>
+ double eval(T pos)
+ {
+ int dim = pos.dims;
+ MatType point(1,dim);
+ for(int j = 0;j < dim;++j)
+ point(0,j) = pos.get(j);
+
+ return model->eval(point)(0);
+ }
+
+ // T1 : Point<vector_type::dims, typename vector_type::stype>
+ // T2 : Point<vector_type::dims, int>
+ template<typename T1, typename T2>
+ double deriv(T1 pos, T2 deriv_order)
+ {
+
+ int dim = pos.dims;
+ MatType point(1,dim);
+ for(int j = 0;j < dim;++j)
+ point(0,j) = pos.get(j);
+
+ std::vector<int> order;
+ for(int j = 0;j < dim;++j)
+ order.push_back(deriv_order.get(j));
+
+ return model->deriv_eval(point, order)(0);
+ }
+
+ void compute_grad()
+ {
+ for(int i = 0;i < spatial_dim;++i)
+ {
+ std::vector<int> ord(spatial_dim, 0);
+ ord[i] = 1;
+ deriv_model[i] = model->derivative(ord);
+ }
+ }
+
+ // T: Point<vector_type::dims, typename vector_type::stype>
+ template<typename T>
+ T eval_grad(T pos)
+ {
+ T res;
+
+ if(!deriv_model[0])
+ compute_grad();
+
+ for(int i = 0;i < spatial_dim;++i)
+ res.get(i) = deriv_model[i]->eval(pos);
+
+ return res;
+ }
+
+};
+
+
+#endif /* REGRESSION_HPP_ */
diff --git a/src/regression/regression_test.cpp b/src/regression/regression_test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0c96fd81923d25f1ba4fa73ec364b296148046fc
--- /dev/null
+++ b/src/regression/regression_test.cpp
@@ -0,0 +1,183 @@
+/*
+ * Unit tests for the regression module
+ * author : Sachin (sthekke@mpi-cbg.de)
+ * date : 28.04.2022
+ *
+ */
+
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+#include "regression.hpp"
+#include "Vector/vector_dist.hpp"
+#include "DMatrix/EMatrix.hpp"
+
+BOOST_AUTO_TEST_SUITE( Regression_test )
+
+
+
+BOOST_AUTO_TEST_CASE ( Regression_local )
+{
+ Box<2,float> domain({0.0,0.0},{1.0,1.0});
+ // Here we define the boundary conditions of our problem
+ size_t bc[2]={PERIODIC,PERIODIC};
+ // extended boundary around the domain, and the processor domain
+ Ghost<2,float> g(0.01);
+
+ using vectorType = vector_dist<2,float, aggregate<double> >;
+ vectorType vd(2048,domain,bc,g);
+ // the scalar is the element at position 0 in the aggregate
+ const int scalar = 0;
+
+ auto it = vd.getDomainIterator();
+ while (it.isNext())
+ {
+ auto key = it.get();
+ double posx = (double)rand() / RAND_MAX;
+ double posy = (double)rand() / RAND_MAX;
+
+ vd.getPos(key)[0] = posx;
+ vd.getPos(key)[1] = posy;
+
+ // Use synthetic data x*y
+ vd.template getProp<scalar>(key) = sin(posx*posy);
+ ++it;
+ }
+ vd.map();
+
+ auto it2 = vd.getDomainIterator();
+ auto NN = vd.getCellList(0.1);
+ /*
+ auto support = RegressionSupport<vectorType, decltype(NN)>(vd, it2, 10, N_PARTICLES, NN);
+
+ std::cout<<"Initialized domain with "<<support.getNumParticles()<<" particles.\n";
+
+ auto model = RegressionModel<2, 0, vectorType>(vd, dom, 6, 2.0);
+
+ double max_err = -1.0;
+ for(double x = 0.75; x < 0.85;x+=0.01)
+ {
+ for(double y = 0.75; y < 0.85; y+=0.01)
+ {
+ Point<2, double> pos{x,y};
+ double val = model->eval(pos);
+ double actual = sin(x*y);
+ double err = std::abs(actual - val);
+ if (err > max_err) max_err = err;
+ }
+ }
+
+ double tolerance = 1e-5;
+ bool check;
+ if (std::abs(max_err) < tolerance)
+ check = true;
+ else
+ check = false;
+ std::cout<<"Max err = "<<max_err<<"\n";
+ BOOST_TEST( check );
+ */
+}
+
+
+
+BOOST_AUTO_TEST_CASE ( Regression_without_domain_initialization)
+{
+ Box<2,float> domain({0.0,0.0},{1.0,1.0});
+ // Here we define the boundary conditions of our problem
+ size_t bc[2]={PERIODIC,PERIODIC};
+ // extended boundary around the domain, and the processor domain
+ Ghost<2,float> g(0.01);
+
+ using vectorType = vector_dist<2,float, aggregate<double> >;
+
+ vectorType vd_orig(1024,domain,bc,g);
+
+ vectorType vd_test(vd_orig.getDecomposition(), 200);
+
+ Vcluster<> & v_cl = create_vcluster();
+ long int N_prc = v_cl.getProcessingUnits();
+
+ if (v_cl.getProcessUnitID() == 0)
+ std::cout<<"Running regression test with "<<N_prc<<" procs.\n";
+
+ // the scalar is the element at position 0 in the aggregate
+ const int scalar = 0;
+
+ // Creating a grid with synthetic data
+ {
+ auto it = vd_orig.getDomainIterator();
+ while (it.isNext())
+ {
+ auto key = it.get();
+ double posx = (double)rand() / RAND_MAX;
+ double posy = (double)rand() / RAND_MAX;
+
+ vd_orig.getPos(key)[0] = posx;
+ vd_orig.getPos(key)[1] = posy;
+
+ // Use synthetic data x*y
+ vd_orig.template getProp<scalar>(key) = sin(posx*posy);
+ ++it;
+ }
+ vd_orig.map();
+ }
+
+ // Creating test points
+ {
+ auto it = vd_test.getDomainIterator();
+ while (it.isNext())
+ {
+ auto key = it.get();
+ double posx = (double)rand() / RAND_MAX;
+ double posy = (double)rand() / RAND_MAX;
+
+ vd_test.getPos(key)[0] = posx;
+ vd_test.getPos(key)[1] = posy;
+
+ vd_test.template getProp<scalar>(key) = 0.0;
+
+ ++it;
+ }
+ vd_test.map();
+ }
+
+ auto model = RegressionModel<2, 0>(vd_orig, 1e-6);
+
+ double max_err = -1.0;
+ // Checking the error
+ {
+ auto it = vd_test.getDomainIterator();
+ while (it.isNext())
+ {
+ auto key = it.get();
+ Point<2, double> pos = {vd_test.getPos(key)[0], vd_test.getPos(key)[1]};
+ double val = model.eval(pos);
+ double actual = sin(pos[0]*pos[1]);
+ double err = std::abs(actual - val);
+ if (err > max_err) max_err = err;
+
+ vd_test.template getProp<scalar>(key) = val;
+
+ ++it;
+ }
+ vd_test.ghost_get<scalar>();
+ }
+
+ v_cl.max(max_err);
+ v_cl.execute();
+ if (v_cl.getProcessUnitID() == 0)
+ std::cout << "Maximum error: " << max_err << "\n";
+
+ double tolerance = 1e-5;
+ bool check;
+ if (std::abs(max_err) < tolerance)
+ check = true;
+ else
+ check = false;
+
+ BOOST_TEST( check );
+
+
+}
+
+
+BOOST_AUTO_TEST_SUITE_END()