Skip to content
Snippets Groups Projects
Commit cf94da13 authored by Sachin Krishnan T V's avatar Sachin Krishnan T V
Browse files

Add PolyLevelset class

parent 78491429
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/CMakeLists.txt
+ 1
− 0
View file @ cf94da13
......@@ -76,6 +76,7 @@ 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
......
src/regression/poly_levelset.hpp 0 → 100644
+ 110
− 0
View file @ cf94da13
/*
* 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 vector_type, typename MatType = EMatrixXd, typename VecType = EVectorXd>
class PolyLevelset
{
minter::LevelsetPoly<spatial_dim, MatType, VecType> *model;
public:
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;
}
// TODO: Make the return types more generic
double eval(Point<vector_type::dims, typename vector_type::stype> 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);
}
double deriv(Point<vector_type::dims, typename vector_type::stype> pos, \
Point<vector_type::dims, int> 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);
}
Point<vector_type::dims, typename vector_type::stype> estimate_normals_at(Point<vector_type::dims, typename vector_type::stype> pos)
{
int dim = pos.dims;
MatType point(1,dim);
for(int j = 0;j < dim;++j)
point(0,j) = pos.get(j);
Point<vector_type::dims, typename vector_type::stype> 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;
}
double estimate_mean_curvature_at(Point<vector_type::dims, typename vector_type::stype> 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
src/regression/poly_levelset_test.cpp 0 → 100644
+ 90
− 0
View file @ cf94da13
/*
* 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,float> 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,float> g(0.01);
using vectorType = vector_dist<3,float, 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 = new PolyLevelset<3, vectorType>(vd, 1e-4);
double max_err = -1.0;
auto it2 = vd.getDomainIterator();
while (it2.isNext())
{
auto key = it2.get();
vd.template getProp<mean_curvature>(key) = model->estimate_mean_curvature_at(vd.getPos(key));
// 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()
src/regression/regression.hpp
+ 8
− 153
View file @ cf94da13
......@@ -47,46 +47,7 @@ public:
}
};
// template<int spatial_dim, typename MatType, typename VecType>
// class OpenFPMPolyModel
// {
// public:
// // minter::PolyModel<spatial_dim, typename MatType::BaseMatrix, typename VecType::BaseMatrix> *model = nullptr;
// minter::PolyModel<spatial_dim, MatType, VecType> *model = nullptr;
// // OpenFPMPolyModel(minter::PolyModel<spatial_dim, typename MatType::BaseMatrix, typename VecType::BaseMatrix> *mdl) : model(mdl) {}
// OpenFPMPolyModel(minter::PolyModel<spatial_dim, MatType, VecType> *mdl) : model(mdl) {}
// // TODO: Make the return types more generic
// double eval(Point<spatial_dim, double> 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);
// }
// double deriv(Point<spatial_dim, double> pos, \
// Point<spatial_dim, int> 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);
// }
// };
template<int spatial_dim, unsigned int prp_id, typename vector_type, typename MatType, typename VecType>
template<int spatial_dim, unsigned int prp_id, typename vector_type, typename MatType = EMatrixXd, typename VecType = EVectorXd>
class RegressionModel
{
......@@ -95,7 +56,7 @@ public:
minter::PolyModel<spatial_dim, MatType, VecType> *deriv_model[spatial_dim];
template<typename dom_type>
RegressionModel(vector_type &vd, dom_type &domain, unsigned int poly_degree, float lp_degree)
RegressionModel(vector_type &vd, dom_type &domain, unsigned int poly_degree, float lp_degree = 2.0)
{
int num_particles = domain->getNumParticles();
int dim = vector_type::dims;
......@@ -111,21 +72,17 @@ public:
values(i) = vd.template getProp<prp_id>(keys.get(i));
}
// std::cout << boost::core::demangle(typeid(decltype(obj)).name()) << '\n';
// construct polynomial model (degree 4)
// auto mdl = new minter::PolyModel<spatial_dim, typename MatType::BaseMatrix, typename VecType::BaseMatrix>(points, values, poly_degree, lp_degree);
// construct polynomial model
model = new minter::PolyModel<spatial_dim, MatType, VecType>(points, values, poly_degree, lp_degree);
// model = new minter::PolyModel<spatial_dim, MatType, VecType>(mdl);
// 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
RegressionModel(vector_type &vd, unsigned int poly_degree)
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;
......@@ -147,11 +104,8 @@ public:
++i;
}
// construct polynomial model (degree 4)
//auto mdl = new minter::PolyModel<spatial_dim, typename MatType::BaseMatrix, typename VecType::BaseMatrix>(points, values, poly_degree, 2.0);
// auto mdl = new minter::PolyModel<spatial_dim, MatType, VecType>(points, values, poly_degree, 2.0);
model = new minter::PolyModel<spatial_dim, MatType, VecType>(points, values, poly_degree, 2.0);
// 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;
......@@ -182,7 +136,6 @@ public:
int poly_degree = 1;
double error = -1.0;
// minter::PolyModel<spatial_dim, typename MatType::BaseMatrix, typename VecType::BaseMatrix> *mdl = nullptr;
minter::PolyModel<spatial_dim, MatType, VecType> *mdl = nullptr;
do
......@@ -192,7 +145,6 @@ public:
delete mdl;
// construct polynomial model
// mdl = new minter::PolyModel<spatial_dim, typename MatType::BaseMatrix, typename VecType::BaseMatrix>(points, values, poly_degree, 2.0);
mdl = new minter::PolyModel<spatial_dim, MatType, VecType>(points, values, poly_degree, 2.0);
// check if linf_error is within the tolerance
......@@ -208,8 +160,8 @@ public:
}while(error > tolerance);
// model = new OpenFPMPolyModel<spatial_dim, MatType, VecType>(mdl);
model = mdl;
for(i = 0;i < dim;++i)
deriv_model[i] = nullptr;
......@@ -228,7 +180,6 @@ public:
}
}
// TODO: Make the return types more generic
double eval(Point<vector_type::dims, typename vector_type::stype> pos)
{
int dim = pos.dims;
......@@ -261,7 +212,6 @@ public:
{
std::vector<int> ord(spatial_dim, 0);
ord[i] = 1;
// deriv_model[i] = new OpenFPMPolyModel<spatial_dim, MatType, VecType>(model->model->derivative(ord));
deriv_model[i] = model->derivative(ord);
}
}
......@@ -270,11 +220,6 @@ public:
{
Point<vector_type::dims, typename vector_type::stype> res;
// int dim = pos.dims;
// typename MatType::BaseMatrix point(1,dim);
// for(int j = 0;j < dim;++j)
// point(0,j) = pos.get(j);
if(!deriv_model[0])
compute_grad();
......@@ -286,94 +231,4 @@ public:
};
template<int spatial_dim, typename vector_type, typename MatType, typename VecType>
class LevelsetPoly
{
minter::LevelsetPoly<spatial_dim, MatType, VecType> *model;
public:
LevelsetPoly(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 (degree 4)
model = new minter::LevelsetPoly<spatial_dim, MatType, VecType>(points, tol);
}
~LevelsetPoly()
{
if(model)
delete model;
}
// TODO: Make the return types more generic
double eval(Point<vector_type::dims, typename vector_type::stype> 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);
}
double deriv(Point<vector_type::dims, typename vector_type::stype> pos, \
Point<vector_type::dims, int> 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);
}
Point<vector_type::dims, typename vector_type::stype> estimate_normals_at(Point<vector_type::dims, typename vector_type::stype> pos)
{
int dim = pos.dims;
MatType point(1,dim);
for(int j = 0;j < dim;++j)
point(0,j) = pos.get(j);
Point<vector_type::dims, typename vector_type::stype> 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;
}
double estimate_mean_curvatures_at(Point<vector_type::dims, typename vector_type::stype> 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 /* REGRESSION_HPP_ */
\ No newline at end of file
src/regression/regression_test.cpp
+ 7
− 13
View file @ cf94da13
......@@ -5,10 +5,6 @@
*
*/
#ifndef OPENFPM_NUMERICS_REGRESSION_TEST_HPP_
#define OPENFPM_NUMERICS_REGRESSION_TEST_HPP_
#define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>
#include "regression.hpp"
......@@ -28,7 +24,7 @@ BOOST_AUTO_TEST_CASE ( Regression_domain_initialization )
Ghost<2,float> g(0.01);
using vectorType = vector_dist<2,float, aggregate<double> >;
vectorType vd(4096,domain,bc,g);
vectorType vd(2048,domain,bc,g);
// the scalar is the element at position 0 in the aggregate
const int scalar = 0;
......@@ -54,7 +50,7 @@ BOOST_AUTO_TEST_CASE ( Regression_domain_initialization )
std::cout<<"Initialized domain with "<<dom->getNumParticles()<<" particles.\n";
auto model = new RegressionModel<2, 0, vectorType, EMatrixXd, EVectorXd>(vd, dom, 6, 2.0);
auto model = new 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)
......@@ -69,7 +65,7 @@ BOOST_AUTO_TEST_CASE ( Regression_domain_initialization )
}
}
double tolerance = 1e-7;
double tolerance = 1e-5;
bool check;
if (std::abs(max_err) < tolerance)
check = true;
......@@ -97,7 +93,7 @@ BOOST_AUTO_TEST_CASE ( Regression_without_domain_initialization)
Ghost<2,float> g(0.01);
using vectorType = vector_dist<2,float, aggregate<double> >;
vectorType vd(4096,domain,bc,g);
vectorType vd(2048,domain,bc,g);
// the scalar is the element at position 0 in the aggregate
const int scalar = 0;
......@@ -117,7 +113,7 @@ BOOST_AUTO_TEST_CASE ( Regression_without_domain_initialization)
}
vd.map();
auto model = new RegressionModel<2, 0, vectorType, EMatrixXd, EVectorXd>(vd, static_cast<unsigned int>(10));
auto model = new RegressionModel<2, 0, vectorType>(vd, static_cast<unsigned int>(10));
double max_err = -1.0;
for(double x = 0.75; x < 0.85;x+=0.01)
......@@ -132,7 +128,7 @@ BOOST_AUTO_TEST_CASE ( Regression_without_domain_initialization)
}
}
double tolerance = 1e-7;
double tolerance = 1e-5;
bool check;
if (std::abs(max_err) < tolerance)
check = true;
......@@ -147,6 +143,4 @@ BOOST_AUTO_TEST_CASE ( Regression_without_domain_initialization)
}
BOOST_AUTO_TEST_SUITE_END()
#endif /* OPENFPM_NUMERICS_REGRESSION_TEST_HPP_ */
BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment