diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 0101c12e5c91f607f6b1e934c97ac818bbaf66b9..82881e8d9cc69bdfdaebc618ce8a2932706a29ab 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -49,6 +49,7 @@ add_executable(numerics ${OPENFPM_INIT_FILE} ${CUDA_SOURCES}
# 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
)
add_dependencies(numerics ofpmmemory)
@@ -92,12 +93,14 @@ target_include_directories (numerics PUBLIC ${Boost_INCLUDE_DIRS})
target_include_directories (numerics PUBLIC ${Vc_INCLUDE_DIR})
target_include_directories (numerics PUBLIC ${ALPAKA_ROOT}/include)
target_include_directories (numerics PUBLIC ${MPI_C_INCLUDE_DIRS})
+target_include_directories (numerics PUBLIC /datapot/lschulze/Applications/minter/src)
if(EIGEN3_FOUND)
target_include_directories (numerics PUBLIC /usr/local/include)
target_include_directories (numerics PUBLIC ${EIGEN3_INCLUDE_DIR})
endif()
link_directories(${PARMETIS_ROOT} ${METIS_ROOT})
+link_directories(/datapot/lschulze/Applications/minter/build)
target_link_libraries(numerics ${Boost_LIBRARIES})
target_link_libraries(numerics ${HDF5_LIBRARIES})
target_link_libraries(numerics ${PARMETIS_LIBRARIES})
@@ -105,6 +108,8 @@ 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})
+target_link_libraries(numerics minter)
+
if(PETSC_FOUND)
target_include_directories (numerics PUBLIC ${PETSC_INCLUDES})
target_link_libraries(numerics ${PETSC_LIBRARIES})
@@ -129,8 +134,8 @@ if (TEST_COVERAGE)
endif()
## Regression
-target_include_directories(numerics PUBLIC /Users/lschulze/minter/src/)
-target_link_libraries(numerics -L/Users/lschulze/minter/build minter)
+target_include_directories(numerics PUBLIC /datapot/lschulze/Applications/minter/src/)
+target_link_libraries(numerics -L/datapot/lschulze/Applications/minter/build minter)
##
# Request that particles be built with -std=c++11
diff --git a/src/DCPSE/SupportBuilder.hpp b/src/DCPSE/SupportBuilder.hpp
index 7702d4345c21e0363a70dec8919634c89176a761..9cd3a8e8f4ac9190bb79034cba69615cc677df48 100644
--- a/src/DCPSE/SupportBuilder.hpp
+++ b/src/DCPSE/SupportBuilder.hpp
@@ -17,7 +17,8 @@
enum support_options
{
N_PARTICLES,
- RADIUS
+ RADIUS,
+ AT_LEAST_N_PARTICLES
};
template<typename vector_type>
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..30aa41b2b819428586d71f4b291beeb26d7d70f1
--- /dev/null
+++ b/src/level_set/particle_cp/pcp_unit_tests.cpp
@@ -0,0 +1,196 @@
+/* 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 surf_flag = 4;
+ constexpr int ref_cp = 5;
+ typedef vector_dist<3, double, aggregate<double, Point<3, double>, Point<3, double>, double, int, Point<3, double>>> particles;
+ // | | | | | |
+ // SDF closest point | curvature surf flag |
+ // surface normal reference closest point
+
+ particles vd(0, domain, bc, g, DEC_GRAN(512));
+ openfpm::vector<std::string> names({"sdf","cp","normal","curvature","C flag", "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, taylor4, sdf, cp, normal, curvature, num_coeffs> pcprdist(vd, rdistoptions);
+ pcprdist.run_redistancing();
+
+ //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))/(2*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/regression.hpp b/src/regression/regression.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..30aefa14642c33e582f9634a647ac8111ee11ecf
--- /dev/null
+++ b/src/regression/regression.hpp
@@ -0,0 +1,461 @@
+/*
+ * 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 "/datapot/lschulze/Applications/minter/include/minter.h"
+
+
+template<typename vector_type_support, typename NN_type>
+class RegressionSupport
+{
+ //std::vector<size_t> keys;
+ 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_ */