From 6d9e0ef2429bbf73a76804691de756f545477f37 Mon Sep 17 00:00:00 2001
From: lschulze <lschulze@mpi-cbg.de>
Date: Mon, 6 Dec 2021 19:16:47 +0100
Subject: [PATCH] Add draft of quadratic interpolation. Note that
MonomialBasis.hpp needs to be set correctly.
---
src/level_set/particle_cp/particle_cp.hpp | 38 +++++++++++------------
1 file changed, 18 insertions(+), 20 deletions(-)
diff --git a/src/level_set/particle_cp/particle_cp.hpp b/src/level_set/particle_cp/particle_cp.hpp
index 45982cfa..cd8642bc 100644
--- a/src/level_set/particle_cp/particle_cp.hpp
+++ b/src/level_set/particle_cp/particle_cp.hpp
@@ -131,7 +131,6 @@ private:
double r2 = norm2(dr);
if ((sqrt(r2) < (1.3*redistOptions.H)) && !vd_in.template getProp<vd_in_close_part>(bkey) && (sgn_a != sgn_b)) isclose = 1;
-
//int isclose = (abs(vd_in.template getProp<vd_in_sdf>(akey)) < (redistOptions.H/2.0 + 1e-8));
if (r2 < r_cutoff2)
@@ -188,7 +187,7 @@ private:
void interpolate_sdf_field()
{
- int verbose = 1;
+ int verbose = 0;
auto NN_s = vd_s.getCellList(sqrt(r_cutoff2));
vd_s.updateCellList(NN_s);
auto part = vd_s.getDomainIterator();
@@ -217,9 +216,9 @@ private:
unsigned int degrees[2];
degrees[0] = 1;
-<<<<<<< Updated upstream
degrees[1] = 1;
MonomialBasis<2> m(degrees, 1);
+ //MonomialBasis<2> m(4);
if (redistOptions.polynomial_degree == "quadratic")
{
@@ -233,12 +232,6 @@ private:
//MonomialBasis<2> m(4);
}
-=======
- degrees[1] = 2;
- // order limit 4 corresponds to bicubic basis functions m(4)
- //MonomialBasis<2> m(degrees,1);
- MonomialBasis<2> m(3);
->>>>>>> Stashed changes
// std::cout<<"m size"<<m.size()<<std::endl;
VandermondeRowBuilder<2, double> vrb(m);
@@ -261,8 +254,8 @@ private:
// debug given data
//std::cout<<std::setprecision(15)<<xb[0]<<"\t"<<xb[1]<<"\t"<<vd.getProp<sdf>(b)<<std::endl;
- xb[0] = 0.1;
- xb[1] = 0.5;
+ //xb[0] = 0.1;
+ //xb[1] = 0.5;
// Fill phi-vector from the right hand side
phi[neib] = vd_s.template getProp<vd_s_sdf>(b);
@@ -273,19 +266,13 @@ private:
}
EMatrix<double, Eigen::Dynamic, 1> c(m.size(), 1);
- std::cout<<"testinghere,..."<<std::endl;
+
c = V.completeOrthogonalDecomposition().solve(phi);
-<<<<<<< Updated upstream
-=======
- std::cout<<"testinghere,..."<<std::endl;
->>>>>>> Stashed changes
for (int k = 0; k<m.size(); k++)
{
- std::cout<<"testinghere,..."<<k<<std::endl;
vd_s.template getProp<interpol_coeff>(a)[k] = c[k];
}
- std::cout<<"testinghere,..."<<std::endl;
if (verbose)
{
@@ -384,6 +371,10 @@ private:
val = vd_s.template getProp<vd_s_sdf>(b_min);
x[0] = vd_s.getPos(b_min)[0];
x[1] = vd_s.getPos(b_min)[1];
+
+ EMatrix<double, Eigen::Dynamic, 1> x00(2,1);
+ x00[0] = x[0];
+ x00[1] = x[1];
// take the interpolation polynomial of the particle closest to the surface
for (int k = 0 ; k < msize ; k++)
{
@@ -391,7 +382,7 @@ private:
c[k] = vd_s.template getProp<interpol_coeff>(b_min)[k];
}
- //if (a.getKey() == 78) verbose = 1;
+ if (a.getKey() == 32290) verbose = 1;
if(verbose)
{
@@ -465,16 +456,23 @@ private:
{
std::cout<<"x_final: "<<x[0]<<", "<<x[1]<<std::endl;
std::cout<<"p(x_final) :"<<get_p(x, c, redistOptions.polynomial_degree)<<std::endl;
- std::cout<<"\nabla p(x_final)"<<get_dpdx(x, c, redistOptions.polynomial_degree)<<", "<<get_dpdy(x, c, redistOptions.polynomial_degree)<<std::endl;
+ std::cout<<"nabla p(x_final)"<<get_dpdx(x, c, redistOptions.polynomial_degree)<<", "<<get_dpdy(x, c, redistOptions.polynomial_degree)<<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;
verbose = 0;
+ if (((abs(x00[0] - x[0]))>sqrt(r_cutoff2))||((abs(x00[1] - x[1]))>sqrt(r_cutoff2)))
+ {
+ std::cout<<"straying out of local neighborhood.."<<std::endl;
+ std::cout<<"computed sdf: "<<vd_in.template getProp<vd_in_sdf>(a)<<std::endl;
+ std::cout<<"analytical value: "<<vd_in.template getProp<8>(a)<<", diff: "<<abs(vd_in.template getProp<vd_in_sdf>(a) - vd_in.template getProp<8>(a))<<std::endl;
+ }
++part;
}
}
+ //todo: template these
inline double get_p(EMatrix<double, Eigen::Dynamic, 1> xvector, EMatrix<double, Eigen::Dynamic, 1> c, std::string polynomialdegree)
{
const double x = xvector[0];
--
GitLab