diff --git a/src/level_set/particle_cp/particle_cp.hpp b/src/level_set/particle_cp/particle_cp.hpp
index ca6672071a084692f33fe059e1501886bb3002da..ba1e7f86d2b75e06ccdc4a28b8f00eeb8885ccf0 100644
--- a/src/level_set/particle_cp/particle_cp.hpp
+++ b/src/level_set/particle_cp/particle_cp.hpp
@@ -33,7 +33,7 @@
#include <chrono>
#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], EVectorXd>>;
+template<unsigned int dim, unsigned int n_c> using particles_surface = vector_dist<dim, double, aggregate<int, int, double, double[dim], EVectorXd>>;
struct Redist_options
{
size_t max_iter = 1000; // params for the Newton algorithm for the solution of the constrained
@@ -61,15 +61,7 @@ struct Redist_options
int only_narrowband = 1; // only redistance particles with phi < sampling_radius, or all particles if only_narrowband = 0
};
-// Add new polynomial degrees here in case new interpolation schemes are implemented.
-// Template typenames for initialization of vd_s
-struct quadratic {};
-struct bicubic {};
-struct taylor4 {};
-struct minter_polynomial {};
-
-
-template <typename particles_in_type, typename polynomial_degree, size_t phi_field, size_t closest_point_field, size_t normal_field, size_t curvature_field, unsigned int num_minter_coeffs>
+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:
@@ -80,12 +72,7 @@ public:
minterModelpcp(RegressionModel<dim, vd_s_sdf>(
redistOptions.minter_poly_degree, redistOptions.minter_lp_degree))
{
- // initialize polynomial degrees
- if (std::is_same<polynomial_degree,quadratic>::value) polynomialDegree = "quadratic";
- else if (std::is_same<polynomial_degree,bicubic>::value) polynomialDegree = "bicubic";
- else if (std::is_same<polynomial_degree,taylor4>::value) polynomialDegree = "taylor4";
- else if (std::is_same<polynomial_degree,minter_polynomial>::value) polynomialDegree = "minterpolation";
- else throw std::invalid_argument("Invalid polynomial degree given. Valid choices currently are quadratic, bicubic, taylor4, minter_polynomial.");
+ // Constructor
}
void run_redistancing()
@@ -93,7 +80,6 @@ public:
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;
- std::cout<<"Minterpol variable is "<<minterpol<<std::endl;
}
detect_surface_particles();
@@ -108,8 +94,7 @@ private:
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 interpol_coeff = 4;
- static constexpr size_t minter_coeff = 5;
+ 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
@@ -119,21 +104,15 @@ private:
static constexpr size_t vd_in_cp = closest_point_field;
Redist_options redistOptions;
- static constexpr unsigned int minterpol = (num_minter_coeffs > 0) ? 1 : 0;
particles_in_type &vd_in;
static constexpr unsigned int dim = particles_in_type::dims;
- static constexpr unsigned int n_c = (minterpol == 1) ? num_minter_coeffs :
- (dim == 2 && std::is_same<polynomial_degree,quadratic>::value) ? 6 :
- (dim == 2 && std::is_same<polynomial_degree,bicubic>::value) ? 16 :
- (dim == 2 && std::is_same<polynomial_degree,taylor4>::value) ? 15 :
- (dim == 3 && std::is_same<polynomial_degree,taylor4>::value) ? 35 : 1;
+ 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;
- std::string polynomialDegree;
RegressionModel<dim, vd_s_sdf> minterModelpcp;
int return_sign(double phi)
@@ -227,7 +206,7 @@ private:
vd_s.template ghost_get<vd_s_sdf>();
double r_cutoff_celllist = sqrt(r_cutoff2);
- if (minterpol && (redistOptions.min_num_particles != 0)) r_cutoff_celllist = redistOptions.r_cutoff_factor_min_num_particles*redistOptions.H;
+ 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();
@@ -247,135 +226,42 @@ private:
Point<dim, double> xa = vd_s.getPos(a);
int neib = 0;
int k_project = 0;
- // initialize monomial basis functions m, Vandermonde row builder vrb, Vandermonde matrix V, the right-hand
- // side vector phi and the solution vector that contains the coefficients of the interpolation polynomial, c.
-
- // The functions called here are originally from the DCPSE work. Minter regression is better suited for the application.
- if (!minterpol){
- MonomialBasis <dim> m(4);
-
- if (polynomialDegree == "quadratic") {
- unsigned int degrees[dim];
- for (int k = 0; k < dim; k++) degrees[k] = 1;
- m = MonomialBasis<dim>(degrees, 1);
- } else if (polynomialDegree == "taylor4") {
- unsigned int degrees[dim];
- for (int k = 0; k < dim; k++) degrees[k] = 1;
- unsigned int degrees_2 = 3;
- if (dim == 3) degrees_2 = 2;
- m = MonomialBasis<dim>(degrees, degrees_2);
- } else if (polynomialDegree == "bicubic") {
- // do nothing as it has already been initialized like this
- }
-
- if (m.size() > num_neibs_a) message_insufficient_support = 1;
-
- VandermondeRowBuilder<dim, double> vrb(m);
-
- EMatrix<double, Eigen::Dynamic, Eigen::Dynamic> V(num_neibs_a, m.size());
- EMatrix<double, Eigen::Dynamic, 1> phi(num_neibs_a, 1);
- EMatrix<double, Eigen::Dynamic, 1> c(m.size(), 1);
-
- // iterate over the support of the central particle, in order to build the Vandermonde matrix and the
- // right-hand side vector
- auto Np = NN_s.template getNNIterator<NO_CHECK>(NN_s.getCell(vd_s.getPos(a)));
-
- while (Np.isNext()) {
- vect_dist_key_dx b = Np.get();
- Point<dim, double> xb = vd_s.getPos(b);
- Point<dim, double> dr = xa - xb;
- double r2 = norm2(dr);
-
- if (r2 < r_cutoff2) {
- // Fill phi-vector from the right hand side
- phi[neib] = vd_s.template getProp<vd_s_sdf>(b);
- vrb.buildRow(V, neib, xb, 1.0);
-
- ++neib;
- }
- ++Np;
- }
-
- // solve the system of equations using an orthogonal decomposition (This is a solid choice for an
- // overdetermined system of equations with a bad conditioning, which is usually the case for the
- // Vandermonde matrices.
- c = V.completeOrthogonalDecomposition().solve(phi);
-
- // store the coefficients at the particle
- for (int k = 0; k < m.size(); k++) {
- vd_s.template getProp<interpol_coeff>(a)[k] = c[k];
- }
-
- // after interpolation coefficients have been computed and stored, perform Newton-style projections
- // towards the interface in order to obtain a sample.
- // initialise derivatives of the interpolation polynomial and an interim position vector x. The iterations
- // start at the central particle x_a.
-
- // double incr; // this is an optional variable in case the iterations are aborted depending on the magnitude
- // of the increment instead of the magnitude of the polynomial value and the number of iterations.
- EMatrix<double, Eigen::Dynamic, 1> grad_p(dim_r, 1);
- double grad_p_mag2;
-
- EMatrix<double, Eigen::Dynamic, 1> x(dim_r, 1);
- for(int k = 0; k < dim; k++) x[k] = xa[k];
- double p = get_p(x, c, polynomialDegree);
-
- EMatrix<double, Eigen::Dynamic, 1> x_debug(dim_r, 1);
- for(int k = 0; k < dim_r; k++) x_debug[k] = 0.25;
-
- // do the projections.
- while ((abs(p) > redistOptions.tolerance) && (k_project < redistOptions.max_iter))
- {
- grad_p = get_grad_p(x, c, polynomialDegree);
- grad_p_mag2 = grad_p.dot(grad_p);
-
- x = x - p * grad_p / grad_p_mag2;
- p = get_p(x, c, polynomialDegree);
- //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[k];
- }
-
- if (minterpol)
+
+ if(redistOptions.min_num_particles == 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 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;
- vd_s.template getProp<minter_coeff>(a) = minterModel->getCoeffs();
+ auto& minterModel = minterModelpcp.model;
+ vd_s.template getProp<minter_coeff>(a) = minterModel->getCoeffs();
- double grad_p_minter_mag2;
+ 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];
+ 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];
- }
+ 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;
@@ -407,7 +293,7 @@ private:
// 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,interpol_coeff>();
+ vd_s.template ghost_get<vd_s_close_part,vd_s_sample>();
auto NN_s = vd_s.getCellList(redistOptions.sampling_radius);
auto part = vd_in.getDomainIterator();
@@ -484,29 +370,17 @@ private:
EMatrix<double, Eigen::Dynamic, 1> x00x(dim_r,1);
for(int k = 0; k < dim; k++) x00x[k] = 0.0;
- // take the interpolation polynomial of the sample particle closest to the query particle
- for (int k = 0 ; k < n_c ; k++) c[k] = vd_s.template getProp<interpol_coeff>(b_min)[k];
-
auto& model = minterModelpcp.model;
Point<dim, int> derivOrder;
Point<dim, double> grad_p_minter;
- if (minterpol)
- {
- model->setCoeffs(vd_s.template getProp<minter_coeff>(b_min));
- }
+ model->setCoeffs(vd_s.template getProp<minter_coeff>(b_min));
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;
- if (!minterpol)
- {
- std::cout<<"interpol_i(x_0) = "<<get_p(x, c, polynomialDegree)<<std::endl;
- }
- else
- {
- std::cout<<"interpol_i(x_0) = "<<get_p_minter(x, model)<<std::endl;
- }
+
+ std::cout<<"interpol_i(x_0) = "<<get_p_minter(x, model)<<std::endl;
}
// variable containing updated distance at iteration k
@@ -515,17 +389,8 @@ private:
int k_newton = 0;
// calculations needed specifically for k_newton == 0
- if (!minterpol)
- {
- p = get_p(x, c, polynomialDegree);
- grad_p = get_grad_p(x, c, polynomialDegree);
- }
-
- if (minterpol)
- {
- p = get_p_minter(x, model);
- grad_p = get_grad_p_minter(x, model);
- }
+ 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.
@@ -544,8 +409,7 @@ private:
while((nabla_f_norm > redistOptions.tolerance) && (k_newton<redistOptions.max_iter))
{
// gather required derivative values
- if (!minterpol) H_p = get_H_p(x, c, polynomialDegree);
- if (minterpol) H_p = get_H_p_minter(x, model);
+ 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;
@@ -573,16 +437,9 @@ private:
// prepare values for next iteration and update the exit criterion
xax = x - xa;
- if (!minterpol)
- {
- p = get_p(x, c, polynomialDegree);
- grad_p = get_grad_p(x, c, polynomialDegree);
- }
- if (minterpol)
- {
- p = get_p_minter(x, model);
- grad_p = get_grad_p_minter(x, model);
- }
+ 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;
@@ -627,19 +484,11 @@ private:
if(redistOptions.verbose)
{
//std::cout<<"new sdf: "<<vd.getProp<sdf>(a)<<std::endl;
- //std::cout<<"c:\n"<<vd.getProp<interpol_coeff>(a)<<"\nmin_sdf: "<<vd.getProp<min_sdf>(a)<<" , min_sdf_x: "<<vd.getProp<min_sdf_x>(a)<<std::endl;
std::cout<<"x_final: "<<x[0]<<", "<<x[1]<<std::endl;
- if (!minterpol)
- {
- std::cout<<"p(x_final) :"<<get_p(x, c, polynomialDegree)<<std::endl;
- std::cout<<"nabla p(x_final)"<<get_grad_p(x, c, polynomialDegree)<<std::endl;
- }
- else
- {
- 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;
+ 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)
@@ -649,8 +498,7 @@ private:
if (redistOptions.compute_curvatures)
{
- if (!minterpol) H_p = get_H_p(x, c, polynomialDegree);
- if (minterpol) H_p = get_H_p_minter(x, model);
+ H_p = get_H_p_minter(x, model);
// divergence of normalized gradient field:
if (dim == 2)
{
@@ -710,146 +558,6 @@ private:
return(b_min);
}
- // monomial regression polynomial evaluation
- inline double get_p(EMatrix<double, Eigen::Dynamic, 1> xvector, EMatrix<double, Eigen::Dynamic, 1> c, std::string polynomialdegree)
- {
- const double x = xvector[0];
- const double y = xvector[1];
-
- if (dim == 2)
- {
- if (polynomialdegree == "bicubic")
- {
- return(c[0] + c[1]*x + c[2]*x*x + c[3]*x*x*x + c[4]*y + c[5]*x*y + c[6]*x*x*y + c[7]*x*x*x*y + c[8]*y*y + c[9]*x*y*y + c[10]*x*x*y*y + c[11]*x*x*x*y*y +c[12]*y*y*y + c[13]*x*y*y*y + c[14]*x*x*y*y*y + c[15]*x*x*x*y*y*y);
- }
- if (polynomialdegree == "quadratic")
- {
- return(c[0] + c[1]*x + c[2]*x*x + c[3]*y + c[4]*x*y + c[5]*y*y);
- }
- if (polynomialdegree == "taylor4")
- {
- return(c[0] + c[1]*x + c[2]*x*x + c[3]*x*x*x + c[4]*x*x*x*x + c[5]*y + c[6]*x*y + c[7]*x*x*y + c[8]*x*x*x*y + c[9]*y*y + c[10]*y*y*x + c[11]*x*x*y*y + c[12]*y*y*y + c[13]*y*y*y*x + c[14]*y*y*y*y);
- }
- else throw std::invalid_argument("received unknown polynomial degree");
- }
- if (dim == 3)
- {
- const double z = xvector[2];
- if (polynomialdegree == "taylor4")
- {
- return(c[0] + c[1]*x + c[2]*x*x + c[3]*x*x*x + c[4]*x*x*x*x + c[5]*y + c[6]*x*y + c[7]*x*x*y + c[8]*x*x*x*y + c[9]*y*y + c[10]*y*y*x + c[11]*x*x*y*y + c[12]*y*y*y + c[13]*y*y*y*x + c[14]*y*y*y*y
- + c[15]*z + c[16]*x*z + c[17]*x*x*z + c[18]*x*x*x*z + c[19]*y*z + c[20]*x*y*z + c[21]*x*x*y*z + c[22]*y*y*z + c[23]*x*y*y*z + c[24]*y*y*y*z + c[25]*z*z + c[26]*x*z*z + c[27]*x*x*z*z
- + c[28]*y*z*z + c[29]*x*y*z*z + c[30]*y*y*z*z + c[31]*z*z*z + c[32]*x*z*z*z + c[33]*y*z*z*z + c[34]*z*z*z*z);
- }
- else throw std::invalid_argument("received unknown polynomial degree");
- }
- else throw std::invalid_argument("received unknown input dimension. Spatial variable needs to be either 2D or 3D.");
- }
-
- inline EMatrix<double, Eigen::Dynamic, 1> get_grad_p(EMatrix<double, Eigen::Dynamic, 1> xvector, EMatrix<double, Eigen::Dynamic, 1> c, std::string polynomialdegree)
- {
- const double x = xvector[0];
- const double y = xvector[1];
- EMatrix<double, Eigen::Dynamic, 1> grad_p(dim_r, 1);
-
- if (dim == 2)
- {
- if (polynomialdegree == "quadratic")
- {
- grad_p[0] = c[1] + 2*c[2]*x + c[4]*y;
- grad_p[1] = c[3] + c[4]*x + 2*c[5]*y;
- }
- else if (polynomialdegree == "bicubic")
- {
- grad_p[0] = c[1] + 2*c[2]*x + 3*c[3]*x*x + c[5]*y + 2*c[6]*x*y + 3*c[7]*x*x*y + c[9]*y*y + 2*c[10]*x*y*y + 3*c[11]*x*x*y*y + c[13]*y*y*y + 2*c[14]*x*y*y*y + 3*c[15]*x*x*y*y*y;
- grad_p[1] = c[4] + c[5]*x + c[6]*x*x + c[7]*x*x*x + 2*c[8]*y + 2*c[9]*x*y + 2*c[10]*x*x*y + 2*c[11]*x*x*x*y + 3*c[12]*y*y + 3*c[13]*x*y*y + 3*c[14]*x*x*y*y + 3*c[15]*x*x*x*y*y;
- }
- else if (polynomialdegree == "taylor4")
- {
- grad_p[0] = c[1] + 2*c[2]*x + 3*c[3]*x*x + 4*c[4]*x*x*x + c[6]*y + 2*c[7]*x*y + 3*c[8]*x*x*y + c[10]*y*y + 2*c[11]*x*y*y + c[13]*y*y*y;
- grad_p[1] = c[5] + c[6]*x + c[7]*x*x + c[8]*x*x*x + 2*c[9]*y + 2*c[10]*x*y + 2*c[11]*x*x*y + 3*c[12]*y*y + 3*c[13]*y*y*x + 4*c[14]*y*y*y;
- }
- else throw std::invalid_argument("received unknown polynomial degree");
- }
- else if (dim == 3)
- {
- const double z = xvector[2];
- if (polynomialdegree == "taylor4")
- {
- grad_p[0] = c[1] + 2*c[2]*x + 3*c[3]*x*x + 4*c[4]*x*x*x + c[6]*y + 2*c[7]*x*y + 3*c[8]*x*x*y + c[10]*y*y + 2*c[11]*x*y*y + c[13]*y*y*y
- + c[16]*z + 2*c[17]*x*z + 3*c[18]*x*x*z + c[20]*y*z + 2*c[21]*x*y*z + c[23]*y*y*z + c[26]*z*z + 2*c[27]*x*z*z + c[29]*y*z*z + c[32]*z*z*z;
- grad_p[1] = c[5] + c[6]*x + c[7]*x*x + c[8]*x*x*x + 2*c[9]*y + 2*c[10]*x*y + 2*c[11]*x*x*y + 3*c[12]*y*y + 3*c[13]*y*y*x + 4*c[14]*y*y*y
- + c[19]*z + c[20]*x*z + c[21]*x*x*z + 2*c[22]*y*z + 2*c[23]*x*y*z + 3*c[24]*y*y*z + c[28]*z*z + c[29]*x*z*z + 2*c[30]*y*z*z + c[33]*z*z*z;
- grad_p[2] = c[15] + c[16]*x + c[17]*x*x + c[18]*x*x*x + c[19]*y + c[20]*x*y + c[21]*x*x*y + c[22]*y*y + c[23]*x*y*y + c[24]*y*y*y
- + 2*c[25]*z + 2*c[26]*x*z + 2*c[27]*x*x*z + 2*c[28]*y*z + 2*c[29]*x*y*z + 2*c[30]*y*y*z + 3*c[31]*z*z + 3*c[32]*x*z*z + 3*c[33]*y*z*z + 4*c[34]*z*z*z;
- }
- }
- else throw std::invalid_argument("received unknown input dimension. Spatial variable needs to be either 2D or 3D.");
-
- return(grad_p);
-
- }
-
- inline EMatrix<double, Eigen::Dynamic, Eigen::Dynamic> get_H_p(EMatrix<double, Eigen::Dynamic, 1> xvector, EMatrix<double, Eigen::Dynamic, 1> c, std::string polynomialdegree)
- {
- const double x = xvector[0];
- const double y = xvector[1];
- EMatrix<double, Eigen::Dynamic, Eigen::Dynamic> H_p(dim_r, dim_r);
-
- if (dim == 2)
- {
- if (polynomialdegree == "quadratic")
- {
- H_p(0, 0) = 2*c[2];
- H_p(0, 1) = c[4];
- H_p(1, 0) = H_p(0, 1);
- H_p(1, 1) = 2*c[5];
- }
- else if (polynomialdegree == "bicubic")
- {
- H_p(0, 0) = 2*c[2] + 6*c[3]*x + 2*c[6]*y + 6*c[7]*x*y + 2*c[10]*y*y + 6*c[11]*y*y*x + 2*c[14]*y*y*y + 6*c[15]*y*y*y*x;
- H_p(0, 1) = c[5] + 2*c[6]*x + 3*c[7]*x*x + 2*c[9]*y + 4*c[10]*x*y + 6*c[11]*x*x*y + 3*c[13]*y*y + 6*c[14]*x*y*y + 9*c[15]*x*x*y*y;
- H_p(1, 0) = H_p(0, 1);
- H_p(1, 1) = 2*c[8] + 2*c[9]*x + 2*c[10]*x*x + 2*c[11]*x*x*x + 6*c[12]*y + 6*c[13]*x*y + 6*c[14]*x*x*y + 6*c[15]*x*x*x*y;
- }
- else if (polynomialdegree == "taylor4")
- {
- H_p(0, 0) = 2*c[2] + 6*c[3]*x + 12*c[4]*x*x + 2*c[7]*y + 6*c[8]*x*y + 2*c[11]*y*y;
- H_p(0, 1) = c[6] + 2*c[7]*x + 3*c[8]*x*x + 2*c[10]*y +4*c[11]*x*y + 3*c[13]*y*y;
- H_p(1, 0) = H_p(0, 1);
- H_p(1, 1) = 2*c[9] + 2*c[10]*x + 2*c[11]*x*x + 6*c[12]*y + 6*c[13]*x*y + 12*c[14]*y*y;
- }
- else throw std::invalid_argument("received unknown polynomial degree");
- }
- else if (dim == 3)
- {
- const double z = xvector[2];
- if (polynomialdegree == "taylor4")
- {
- H_p(0, 0) = 2*c[2] + 6*c[3]*x + 12*c[4]*x*x + 2*c[7]*y + 6*c[8]*x*y + 2*c[11]*y*y
- + 2*c[17]*z + 6*c[18]*x*z + 2*c[21]*y*z + 2*c[27]*z*z;
- H_p(1, 1) = 2*c[9] + 2*c[10]*x + 2*c[11]*x*x + 6*c[12]*y + 6*c[13]*x*y + 12*c[14]*y*y
- + 2*c[22]*z + 2*c[23]*x*z + 6*c[24]*y*z + 2*c[30]*z*z;
- H_p(2, 2) = 2*c[25] + 2*c[26]*x + 2*c[27]*x*x + 2*c[28]*y + 2*c[29]*x*y + 2*c[30]*y*y
- + 6*c[31]*z + 6*c[32]*x*z + 6*c[33]*y*z + 12*c[34]*z*z;
-
- H_p(0, 1) = (c[6] + 2*c[7]*x + 3*c[8]*x*x + 2*c[10]*y +4*c[11]*x*y + 3*c[13]*y*y
- + c[20]*z + 2*c[21]*x*z + 2*c[23]*y*z + c[29]*z*z);
- H_p(0, 2) = (c[16] + 2*c[17]*x + 3*c[18]*x*x + c[20]*y + 2*c[21]*x*y + c[23]*y*y
- + 2*c[26]*z + 4*c[27]*x*z + 2*c[29]*y*z + 3*c[32]*z*z);
- H_p(1, 2) = (c[19] + c[20]*x + c[21]*x*x + 2*c[22]*y + 2*c[23]*x*y + 3*c[24]*y*y
- + 2*c[28]*z + 2*c[29]*x*z + 4*c[30]*y*z + 3*c[33]*z*z);
- H_p(1, 0) = H_p(0, 1);
- H_p(2, 0) = H_p(0, 2);
- H_p(2, 1) = H_p(1, 2);
- }
- }
- else throw std::invalid_argument("received unknown input dimension. Spatial variable needs to be either 2D or 3D.");
-
- return(H_p);
-
- }
-
// minterface
template<typename PolyType>
inline double get_p_minter(EMatrix<double, Eigen::Dynamic, 1> xvector, PolyType model)
@@ -891,4 +599,3 @@ private:
};
-
diff --git a/src/level_set/particle_cp/pcp_unit_tests.cpp b/src/level_set/particle_cp/pcp_unit_tests.cpp
index 7155b2f8baab13d4f1d9ed86d1cd7a9ce2576d77..4e7d8569512080cb4e01d422205e18c98488ffe5 100644
--- a/src/level_set/particle_cp/pcp_unit_tests.cpp
+++ b/src/level_set/particle_cp/pcp_unit_tests.cpp
@@ -124,7 +124,7 @@ BOOST_AUTO_TEST_CASE( ellipsoid )
static constexpr unsigned int num_coeffs = minter_lp_degree_one_num_coeffs(3, poly_order);
- particle_cp_redistancing<particles, minter_polynomial, sdf, cp, normal, curvature, num_coeffs> pcprdist(vd, rdistoptions);
+ 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();