OpenFPM  5.2.0
Project that contain the implementation of distributed structures
pcp_unit_tests.cpp
1 /* Unit tests for the particle closest point method
2  * Date : 29 June 2023
3  * Author : lschulze
4  */
5 
6 #include<iostream>
7 #include <boost/test/unit_test_log.hpp>
8 #define BOOST_TEST_DYN_LINK
9 #include <boost/test/unit_test.hpp>
10 #include <math.h>
11 //#include <sys/_types/_size_t.h>
12 #include "Vector/vector_dist.hpp"
13 #include "DCPSE/Dcpse.hpp"
14 #include "DCPSE/MonomialBasis.hpp"
15 #include "Draw/DrawParticles.hpp"
16 #include "particle_cp.hpp"
17 #include "Operators/Vector/vector_dist_operators.hpp"
18 #include <chrono>
19 #include "pcp_unit_test_helpfunctions.h"
20 
21 typedef struct EllipseParameters{
22  double origin[3];
23  double radiusA;
24  double radiusB;
25  double radiusC;
26  double eccentricity;
27 } EllipseParams;
28 
29 // Generate an ellipsoid initial levelset signed distance function
30 template<typename particles_type, typename iterator_type, size_t sdf, size_t surf_flag, size_t ref_cp>
31 void initializeLSEllipsoid(particles_type &vd, iterator_type particle_it, const EllipseParams &params, double bandwidth, double perturb_factor, double H)
32 {
33  while(particle_it.isNext())
34  {
35  double posx = particle_it.get().get(0);
36  double posy = particle_it.get().get(1);
37  double posz = particle_it.get().get(2);
38 
39  double ref_cpx = 0.0;
40  double ref_cpy = 0.0;
41  double ref_cpz = 0.0;
42 
43  double dist = DistancePointEllipsoid(params.radiusA, params.radiusB, params.radiusC, abs(posx), abs(posy), abs(posz), ref_cpx, ref_cpy, ref_cpz);
44  if (abs(dist) < bandwidth/2.0)
45  {
46  posx = posx + perturb_factor*randMinusOneToOne()*H;
47  posy = posy + perturb_factor*randMinusOneToOne()*H;
48  posz = posz + perturb_factor*randMinusOneToOne()*H;
49 
50  dist = DistancePointEllipsoid(params.radiusA, params.radiusB, params.radiusC, abs(posx), abs(posy), abs(posz), ref_cpx, ref_cpy, ref_cpz);
51  vd.add();
52  vd.getLastPos()[0] = posx;
53  vd.getLastPos()[1] = posy;
54  vd.getLastPos()[2] = posz;
55 
56  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;
57  vd.template getLastProp<sdf>() = phi_val;
58  vd.template getLastProp<ref_cp>()[0] = return_sign(posx)*ref_cpx;
59  vd.template getLastProp<ref_cp>()[1] = return_sign(posy)*ref_cpy;
60  vd.template getLastProp<ref_cp>()[2] = return_sign(posz)*ref_cpz;
61  vd.template getLastProp<surf_flag>() = 0;
62  }
63 
64  ++particle_it;
65  }
66 }
67 
68 BOOST_AUTO_TEST_SUITE( particle_closest_point_test )
69 
70 BOOST_AUTO_TEST_CASE( ellipsoid )
71 {
72  // simulation params
73  constexpr int poly_order = 4;
74  const double H = 1.0/64.0;
75  const double perturb_factor = 0.3;
76 
77  // pcp params
78  const double bandwidth = 12.0*H;
79  const double regression_rcut_factor = 2.4;
80  const double threshold = 1e-13;
81 
82  // initialize domain and particles
83  const double l = 2.0;
84  Box<3, double> domain({-l/2.0, -l/3.0, -l/3.0}, {l/2.0, l/3.0, l/3.0});
85  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)};
86  size_t bc[3] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
87  Ghost<3, double> g(bandwidth);
88 
89  constexpr int sdf = 0;
90  constexpr int cp = 1;
91  constexpr int normal = 2;
92  constexpr int curvature = 3;
93  constexpr int surf_flag = 4;
94  constexpr int ref_cp = 5;
95  typedef vector_dist<3, double, aggregate<double, Point<3, double>, Point<3, double>, double, int, Point<3, double>>> particles;
96  // | | | | | |
97  // SDF closest point | curvature surface flag |
98  // surface normal reference closest point
99 
100  particles vd(0, domain, bc, g, DEC_GRAN(512));
101  openfpm::vector<std::string> names({"sdf","cp","normal","curvature","ref_cp"});
102  vd.setPropNames(names);
103  EllipseParameters params;
104  for (int k = 0; k < 3; k++) params.origin[k] = 0.0;
105  params.radiusA = 0.75;
106  params.radiusB = 0.5;
107  params.radiusC = 0.5;
108 
109  auto particle_it = DrawParticles::DrawBox(vd, sz, domain, domain);
110 
111  // initialize spurious sdf values and reference solution
112  initializeLSEllipsoid<particles, decltype(particle_it), sdf, surf_flag, ref_cp>(vd, particle_it, params, bandwidth, perturb_factor, H);
113 
114  //vd.write("pcpunittest_init");
115  // initialize and run pcp redistancing
116  Redist_options rdistoptions;
117  rdistoptions.minter_poly_degree = poly_order;
118  rdistoptions.H = H;
119  rdistoptions.r_cutoff_factor = regression_rcut_factor;
120  rdistoptions.sampling_radius = 0.75*bandwidth;
121  rdistoptions.tolerance = threshold;
122  rdistoptions.write_cp = 1;
123  rdistoptions.compute_normals = 1;
124  rdistoptions.compute_curvatures = 1;
125  rdistoptions.only_narrowband = 0;
126 
127  static constexpr unsigned int num_coeffs = minter_lp_degree_one_num_coeffs(3, poly_order);
128 
129  particle_cp_redistancing<particles, sdf, cp, normal, curvature, num_coeffs> pcprdist(vd, rdistoptions);
130  std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
131  pcprdist.run_redistancing();
132  std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
133  std::cout << "Time difference for pcp redistancing = " << std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count() << "[ms]" << std::endl;
134 
135  //vd.write("pcpunittest");
136  // iterate through particles and compute error
137  auto part = vd.getDomainIterator();
138  double sdferr = 0.0;
139  double sdfmaxerr = 0.0;
140  double sdfmeanerr = 0.0;
141  double cperr = 0.0;
142  double cpmaxerr = 0.0;
143  double normalerr = 0.0;
144  double normalmaxerr = 0.0;
145  double curvatureerr = 0.0;
146  double curvaturemaxerr = 0.0;
147  int particlecount = 0;
148 
149  while(part.isNext())
150  {
151  auto a = part.get();
152  Point<3, double> pos = vd.getPos(a);
153  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);
154  Point<3, double> reference_normal;
155  double reference_curvature;
156 
157  // compute reference SDF value w/ reference closest point and check computed SDF value
158  sdferr = abs(vd.getProp<sdf>(a) - reference_sdf);
159  if (sdferr > sdfmaxerr) sdfmaxerr = sdferr;
160  sdfmeanerr += sdferr;
161  ++particlecount;
162 
163  // check computed closest point against reference closest point
164  cperr = norm(vd.getProp<ref_cp>(a) - vd.getProp<cp>(a));
165  if (cperr > cpmaxerr) cpmaxerr = cperr;
166 
167  // compute reference normal and check computed normal
168  reference_normal[0] = 2*vd.getProp<ref_cp>(a)[0]/(params.radiusA*params.radiusA);
169  reference_normal[1] = 2*vd.getProp<ref_cp>(a)[1]/(params.radiusB*params.radiusB);
170  reference_normal[2] = 2*vd.getProp<ref_cp>(a)[2]/(params.radiusC*params.radiusC);
171  reference_normal = return_sign(reference_sdf)*reference_normal/norm(reference_normal);
172  normalerr = norm(reference_normal - vd.getProp<normal>(a));
173  if (normalerr > normalmaxerr) normalmaxerr = normalerr;
174 
175  // compute reference curvature and check computed curvature
176  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));
177  curvatureerr = abs(reference_curvature - vd.getProp<curvature>(a));
178  if (curvatureerr > curvaturemaxerr) curvaturemaxerr = curvatureerr;
179 
180  ++part;
181  }
182  sdfmeanerr = sdfmeanerr/particlecount;
183  std::cout<<"Mean error for sdf is: "<<sdfmeanerr<<std::endl;
184  std::cout<<"Maximum error for sdf is: "<<sdfmaxerr<<std::endl;
185  std::cout<<"Maximum error for the closest point is: "<<cpmaxerr<<std::endl;
186  std::cout<<"Maximum error for surface normal is: "<<normalmaxerr<<std::endl;
187  std::cout<<"Maximum error for curvature is: "<<curvaturemaxerr<<std::endl;
188 
189  double tolerance = 1e-7;
190  bool check;
191  if (std::abs(sdfmaxerr) < tolerance)
192  check = true;
193  else
194  check = false;
195 
196  BOOST_TEST( check );
197 
198 }
199 
200 BOOST_AUTO_TEST_SUITE_END()
Box
This class represent an N-dimensional box.
Definition: Box.hpp:60
DrawParticles::DrawBox
static PointIterator< dim, T, typename vd_type::Decomposition_type > DrawBox(vd_type &vd, size_t(&sz)[dim], Box< dim, T > &domain, Box< dim, T > &sub)
Draw particles in a box.
Definition: DrawParticles.hpp:123
Ghost
Definition: Ghost.hpp:40
Point
This class implement the point shape in an N-dimensional space.
Definition: Point.hpp:28
particle_cp_redistancing
Class for reinitializing a level-set function into a signed distance function using the Closest-Point...
Definition: particle_cp.hpp:64
vector_dist
Distributed vector.
Definition: vector_dist.hpp:176
vector_dist::getProp
auto getProp(vect_dist_key_dx vec_key) -> decltype(vPrp.template get< id >(vec_key.getKey()))
Get the property of an element.
Definition: vector_dist.hpp:655
vector_dist::setPropNames
void setPropNames(const openfpm::vector< std::string > &names)
Set the properties names.
Definition: vector_dist.hpp:2873
vector_dist::getPos
auto getPos(vect_dist_key_dx vec_key) -> decltype(vPos.template get< 0 >(vec_key.getKey()))
Get the position of an element.
Definition: vector_dist.hpp:585
vector_dist::getDomainIterator
vector_dist_iterator getDomainIterator() const
Get an iterator that traverse the particles in the domain.
Definition: vector_dist.hpp:1942
vector_dist::add
void add()
Add local particle.
Definition: vector_dist.hpp:923
vector_dist::getLastPos
auto getLastPos() -> decltype(vPos.template get< 0 >(0))
Get the position of the last element.
Definition: vector_dist.hpp:958
cub::int
KeyT const ValueT ValueT OffsetIteratorT OffsetIteratorT int
[in] The number of segments that comprise the sorting data
Definition: dispatch_radix_sort.cuh:336
particle_cp.hpp
Redist_options
Structure to bundle options for redistancing.
Definition: particle_cp.hpp:33