...
 
Commits (18)
......@@ -36,6 +36,14 @@
ehthumbs.db
Thumbs.db
# CLion IDE-related
.idea/
cmake-build-debug/
cmake-build-release/
# Script-generated build folder
build/
###### Other
AUTHORS
......@@ -47,6 +55,7 @@ README
**/src/Makefile
./Makefile
Makefile.in
Makefile
config.status
configure
numerics
......
......@@ -156,7 +156,10 @@ else()
endif()
if(EIGEN3_FOUND)
set(DEFINE_HAVE_EIGEN "#define HAVE_EIGEN")
message("-- EIGEN3 found.")
set(DEFINE_HAVE_EIGEN "#define HAVE_EIGEN")
else()
message("-- EIGEN3 not found!")
endif()
if(LIBHILBERT_FOUND)
......
openfpm_data @ 0c3cdf5a
Subproject commit da7bebf0e5101b3e0e85e42695d755d843f7c418
Subproject commit 0c3cdf5a817b05578d2eb919683bc64e1a75efdb
......@@ -9,7 +9,7 @@ else()
set(CUDA_SOURCES)
endif()
add_executable(pdata ${OPENFPM_INIT_FILE} ${CUDA_SOURCES} main.cpp Debug/debug_test.cpp Grid/tests/grid_dist_id_HDF5_chckpnt_restart_test.cpp Grid/tests/grid_dist_id_unit_test.cpp Grid/tests/staggered_grid_dist_unit_test.cpp Vector/tests/vector_dist_cell_list_tests.cpp Vector/tests/vector_dist_complex_prp_unit_test.cpp Vector/tests/vector_dist_HDF5_chckpnt_restart_test.cpp Vector/tests/vector_dist_MP_unit_tests.cpp Vector/tests/vector_dist_NN_tests.cpp Vector/tests/vector_dist_unit_test.cpp pdata_performance.cpp Decomposition/tests/CartDecomposition_unit_test.cpp Decomposition/tests/shift_vect_converter_tests.cpp Vector/performance/vector_dist_performance_util.cpp lib/pdata.cpp test_multiple_o.cpp ../openfpm_devices/src/memory/HeapMemory.cpp ../openfpm_devices/src/memory/PtrMemory.cpp ../openfpm_vcluster/src/VCluster/VCluster.cpp ../openfpm_devices/src/Memleak_check.cpp)
add_executable(pdata ${OPENFPM_INIT_FILE} ${CUDA_SOURCES} main.cpp Debug/debug_test.cpp Grid/tests/grid_dist_id_HDF5_chckpnt_restart_test.cpp Grid/tests/grid_dist_id_unit_test.cpp Grid/tests/staggered_grid_dist_unit_test.cpp Vector/tests/vector_dist_cell_list_tests.cpp Vector/tests/vector_dist_complex_prp_unit_test.cpp Vector/tests/vector_dist_HDF5_chckpnt_restart_test.cpp Vector/tests/vector_dist_MP_unit_tests.cpp Vector/tests/vector_dist_NN_tests.cpp Vector/tests/vector_dist_unit_test.cpp pdata_performance.cpp Decomposition/tests/CartDecomposition_unit_test.cpp Decomposition/tests/shift_vect_converter_tests.cpp Vector/performance/vector_dist_performance_util.cpp lib/pdata.cpp test_multiple_o.cpp ../openfpm_devices/src/memory/HeapMemory.cpp ../openfpm_devices/src/memory/PtrMemory.cpp ../openfpm_vcluster/src/VCluster/VCluster.cpp ../openfpm_devices/src/Memleak_check.cpp DCPSE/MonomialBasis.hpp DCPSE/Monomial.hpp DCPSE/tests/MonomialBasis_unit_tests.cpp DCPSE/Vandermonde.hpp DCPSE/VandermondeRowBuilder.hpp DCPSE/tests/Vandermonde_unit_tests.cpp DCPSE/SupportBuilder.hpp DCPSE/tests/Support_unit_tests.cpp DCPSE/DcpseRhs.hpp DCPSE/tests/DcpseRhs_unit_tests.cpp DCPSE/Support.hpp DCPSE/DcpseDiagonalScalingMatrix.hpp DCPSE/Dcpse.hpp DCPSE/tests/Dcpse_unit_tests.cpp)
if ( CMAKE_COMPILER_IS_GNUCC )
target_compile_options(pdata PRIVATE "-Wno-deprecated-declarations")
......@@ -49,6 +49,9 @@ target_include_directories (pdata PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/config)
target_include_directories (pdata PUBLIC ${PETSC_INCLUDES})
target_include_directories (pdata PUBLIC ${HDF5_ROOT}/include)
target_include_directories (pdata PUBLIC ${LIBHILBERT_INCLUDE_DIRS})
if(EIGEN3_FOUND)
target_include_directories (pdata PUBLIC ${EIGEN3_INCLUDE_DIR})
endif()
target_include_directories (pdata PUBLIC ${Boost_INCLUDE_DIRS})
......@@ -63,7 +66,13 @@ if (TEST_COVERAGE)
target_link_libraries(pdata -lgcov --coverage)
endif()
# Debug!
# Hack found at https://github.com/LLNL/scr/issues/130#issuecomment-402815952
IF(MPI_CXX_FOUND)
INCLUDE_DIRECTORIES(${MPI_CXX_INCLUDE_PATH})
# LIST(APPEND SCR_EXTERNAL_LIBS ${MPI_CXX_LIBRARIES})
target_link_libraries(pdata ${MPI_CXX_LIBRARIES})
ENDIF(MPI_CXX_FOUND)
target_include_directories (ofpm_pdata PUBLIC ${CUDA_INCLUDE_DIRS})
target_include_directories (ofpm_pdata PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
......
This diff is collapsed.
//
// Created by tommaso on 29/03/19.
//
#ifndef OPENFPM_PDATA_DCPSEDIAGONALSCALINGMATRIX_HPP
#define OPENFPM_PDATA_DCPSEDIAGONALSCALINGMATRIX_HPP
#include "MonomialBasis.hpp"
#include "Support.hpp"
template <unsigned int dim>
class DcpseDiagonalScalingMatrix
{
private:
const MonomialBasis<dim> monomialBasis;
public:
DcpseDiagonalScalingMatrix(const MonomialBasis<dim> &monomialBasis) : monomialBasis(monomialBasis) {}
template <typename T, typename MatrixType, typename Prop>
void buildMatrix(MatrixType &M, Support<dim, T, Prop> support, T eps);
};
template<unsigned int dim>
template<typename T, typename MatrixType, typename Prop>
void DcpseDiagonalScalingMatrix<dim>::buildMatrix(MatrixType &M, Support<dim, T, Prop> support, T eps)
{
// Check that all the dimension constraints are met
assert(support.size() > monomialBasis.size());
assert(M.rows() == support.size());
assert(M.cols() == support.size());
// Fill the diagonal matrix
M.setZero(); // Make sure the rest of the matrix is zero!
int i = 0;
for (const auto& pt : support.getOffsets())
{
M(i,i) = exp(- norm2(pt) / (2.0 * eps * eps));
++i;
}
}
#endif //OPENFPM_PDATA_DCPSEDIAGONALSCALINGMATRIX_HPP
//
// Created by tommaso on 28/03/19.
//
#ifndef OPENFPM_PDATA_DCPSERHS_HPP
#define OPENFPM_PDATA_DCPSERHS_HPP
#include "MonomialBasis.hpp"
template<unsigned int dim>
class DcpseRhs
{
private:
const Point<dim, unsigned int> differentialSignature;
const MonomialBasis<dim> derivatives;
int sign;
public:
DcpseRhs(const MonomialBasis<dim> &monomialBasis, const Point<dim, unsigned int> &differentialSignature);
template<typename T, typename MatrixType>
MatrixType &getVector(MatrixType &b);
};
// Definitions below
template<unsigned int dim>
DcpseRhs<dim>::DcpseRhs(const MonomialBasis<dim> &monomialBasis,
const Point<dim, unsigned int> &differentialSignature)
: differentialSignature(differentialSignature), derivatives(monomialBasis.getDerivative(differentialSignature))
{
unsigned int order = (Monomial<dim>(differentialSignature)).order();
if (order % 2 == 0)
{
sign = 1;
} else
{
sign = -1;
}
}
template<unsigned int dim>
template<typename T, typename MatrixType>
MatrixType &DcpseRhs<dim>::getVector(MatrixType &b)
{
// The given vector V should have the right dimensions
assert(b.cols() == 1);
assert(b.rows() == derivatives.size());
for (unsigned int i = 0; i < derivatives.size(); ++i)
{
const Monomial<dim> dm = derivatives.getElement(i);
b(i, 0) = sign * dm.evaluate(Point<dim, T>(0));
}
return b;
}
#endif //OPENFPM_PDATA_DCPSERHS_HPP
//
// Created by tommaso on 21/03/19.
//
#ifndef OPENFPM_PDATA_MONOMIALBASISELEMENT_H
#define OPENFPM_PDATA_MONOMIALBASISELEMENT_H
#include <vector>
#include <ostream>
#include "Space/Shape/Point.hpp"
template<unsigned int dim>
class Monomial
{
// The base idea is that this should behave like a pair <unsigned int, std::vector<unsigned int>>
private:
unsigned int sum = 0;
Point<dim, unsigned int> exponents;
unsigned int scalar = 1;
public:
Monomial();
explicit Monomial(const Point<dim, unsigned int> &other, unsigned int scalar = 1);
explicit Monomial(const Point<dim, long int> &other, unsigned int scalar = 1);
explicit Monomial(const unsigned int other[dim]);
Monomial(const Monomial<dim> &other);
Monomial<dim> &operator=(const Monomial<dim> &other);
bool operator==(const Monomial<dim> &other) const;
unsigned int order() const;
unsigned int getExponent(unsigned int i) const;
void setExponent(unsigned int i, unsigned int value);
template<typename T>
T evaluate(const Point<dim, T> x) const;
template<typename T>
T evaluate(const T x[dim]) const;
Monomial<dim> getDerivative(const Point<dim, unsigned int> differentialOrder) const;
template<typename charT, typename traits>
friend std::basic_ostream<charT, traits> &
operator<<(std::basic_ostream<charT, traits> &lhs, Monomial<dim> const &rhs)
{
return lhs << rhs.scalar << " : " << rhs.exponents.toString();
}
private:
void updateSum();
};
////// Definitions below
template<unsigned int dim>
Monomial<dim>::Monomial()
{
exponents.zero();
sum = 0;
}
template<unsigned int dim>
Monomial<dim>::Monomial(const Point<dim, unsigned int> &other, unsigned int scalar)
: exponents(other), scalar(scalar)
{
updateSum();
}
template<unsigned int dim>
Monomial<dim>::Monomial(const Point<dim, long int> &other, unsigned int scalar)
: scalar(scalar)
{
for (size_t i = 0; i < other.nvals; ++i)
{
exponents.get(i) = other.value(i);
}
updateSum();
}
template<unsigned int dim>
Monomial<dim>::Monomial(const unsigned int other[dim])
: Monomial(Point<3, unsigned int>(other)) {}
template<unsigned int dim>
Monomial<dim>::Monomial(const Monomial<dim> &other)
: exponents(other.exponents),
sum(other.sum),
scalar(other.scalar) {}
template<unsigned int dim>
Monomial<dim> &Monomial<dim>::operator=(const Monomial<dim> &other)
{
exponents = other.exponents;
sum = other.sum;
scalar = other.scalar;
return *this;
}
template<unsigned int dim>
void Monomial<dim>::updateSum()
{
unsigned int partialSum = 0;
for (unsigned int i = 0; i < dim; ++i)
{
partialSum += exponents.value(i);
}
sum = partialSum;
}
template<unsigned int dim>
unsigned int Monomial<dim>::order() const
{
return sum;
}
template<unsigned int dim>
unsigned int Monomial<dim>::getExponent(unsigned int i) const
{
return exponents.value(i);
}
template<unsigned int dim>
void Monomial<dim>::setExponent(unsigned int i, unsigned int value)
{
exponents.get(i) = value;
updateSum();
}
template<unsigned int dim>
bool Monomial<dim>::operator==
(const Monomial<dim> &other) const
{
return (exponents == other.exponents) && (scalar == other.scalar);
}
template<unsigned int dim>
template<typename T>
T Monomial<dim>::evaluate(const Point<dim, T> x) const
{
T res = scalar;
for (unsigned int i = 0; i < dim; ++i)
{
res *= pow(x.value(i), getExponent(i));
}
return res;
}
template<unsigned int dim>
Monomial<dim> Monomial<dim>::getDerivative(const Point<dim, unsigned int> differentialOrder) const
{
unsigned int s = scalar;
Point<dim, unsigned int> e(exponents);
for (unsigned int i = 0; i < dim; ++i)
{
unsigned int origExp = e.value(i);
int targetExp = static_cast<int>(origExp) - static_cast<int>(differentialOrder.value(i));
for (int k = origExp; k > targetExp && k >= 0; --k)
{
s *= k;
}
e.get(i) = static_cast<unsigned int>(std::max(targetExp, 0));
}
return Monomial(e, s);
}
template<unsigned int dim>
template<typename T>
T Monomial<dim>::evaluate(const T x[dim]) const
{
return evaluate(Point<dim, T>(x));
}
#endif //OPENFPM_PDATA_MONOMIALBASISELEMENT_H
//
// Created by tommaso on 20/03/19.
//
#ifndef OPENFPM_PDATA_MONOMIALBASIS_H
#define OPENFPM_PDATA_MONOMIALBASIS_H
#include <vector>
#include <Grid/grid_sm.hpp>
#include <Grid/iterators/grid_key_dx_iterator_sub_bc.hpp>
#include "Monomial.hpp"
template<unsigned int dim>
class MonomialBasis
{
private:
std::vector<Monomial<dim>> basis;
public:
MonomialBasis(const std::vector<unsigned int> &degrees, unsigned int convergenceOrder);
MonomialBasis(unsigned int degrees[dim], unsigned int convergenceOrder);
// explicit MonomialBasis(Point<dim, unsigned int> degrees, unsigned int convergenceOrder);
explicit MonomialBasis(const std::vector<Monomial<dim>> &basis) : basis(basis) {}
MonomialBasis(const MonomialBasis &other);
MonomialBasis &operator=(const MonomialBasis &other);
unsigned int size() const;
const Monomial<dim> &getElement(unsigned int i) const;
Monomial<dim> &getElement(unsigned int i);
const std::vector<Monomial<dim>> &getElements() const;
MonomialBasis<dim> getDerivative(Point<dim, unsigned int> differentialOrder) const;
bool operator==(const MonomialBasis &other) const;
template<typename charT, typename traits>
friend std::basic_ostream<charT, traits> &
operator<<(std::basic_ostream<charT, traits> &lhs, MonomialBasis<dim> const &rhs)
{
lhs << "MonomialBasis: size=" << rhs.size() << ", elements={ ";
for (const auto &el : rhs.getElements())
{
lhs << "(" << el << ") ";
}
lhs << "}" << std::endl;
return lhs;
}
private:
void generateBasis(std::vector<unsigned int> m, unsigned int r);
};
//// Definitions below
template<unsigned int dim>
MonomialBasis<dim>::MonomialBasis(const std::vector<unsigned int> &degrees, unsigned int convergenceOrder)
{
generateBasis(degrees, convergenceOrder);
}
template<unsigned int dim>
MonomialBasis<dim>::MonomialBasis(unsigned int *degrees, unsigned int convergenceOrder)
: MonomialBasis(std::vector<unsigned int>(degrees, degrees + dim), convergenceOrder) {}
template<unsigned int dim>
MonomialBasis<dim>::MonomialBasis(const MonomialBasis &other)
{
basis = other.basis; // Here it works because both std::vector and Monomial perform a deep copy.
}
template<unsigned int dim>
MonomialBasis<dim> &MonomialBasis<dim>::operator=(const MonomialBasis &other)
{
basis = other.basis; // Here it works because both std::vector and Monomial perform a deep copy.
return *this;
}
template<unsigned int dim>
unsigned int MonomialBasis<dim>::size() const
{
return basis.size();
}
template<unsigned int dim>
const Monomial<dim> &MonomialBasis<dim>::getElement(unsigned int i) const
{
return basis[i];
}
template<unsigned int dim>
Monomial<dim> &MonomialBasis<dim>::getElement(unsigned int i)
{
return basis[i];
}
template<unsigned int dim>
void MonomialBasis<dim>::generateBasis(std::vector<unsigned int> m, unsigned int r)
{
// Compute the vector of actual dimensions to iterate over
// NOTE: each index can go up to sum(m)+r
unsigned int mSum = std::accumulate(m.begin(), m.end(), 0U);
unsigned int orderLimit = mSum + r;
size_t dimensions[dim];
std::fill(dimensions, dimensions + dim, orderLimit);
// Now initialize grid with appropriate size, then start-stop points and boundary conditions for the iterator
grid_sm<dim, void> grid(dimensions);
long int startV[dim] = {}; // 0-initialized
grid_key_dx<dim, long int> start(startV);
grid_key_dx<dim, long int> stop(dimensions);
size_t bc[dim];
std::fill(bc, bc + dim, NON_PERIODIC);
grid_key_dx_iterator_sub_bc<dim> it(grid, start, stop, bc);
// Finally compute alpha_min
// unsigned char alphaMin = static_cast<unsigned char>(!(mSum % 2)); // if mSum is even, alpha_min must be 1
unsigned char alphaMin = 0; // we want to always have 1 in the basis
while (it.isNext())
{
Point<dim, long int> p = it.get().get_k();
Monomial<dim> candidateBasisElement(p);
// Filter out the elements which don't fullfil the theoretical condition for being in the vandermonde matrix
if (candidateBasisElement.order() < orderLimit && candidateBasisElement.order() >= alphaMin)
{
basis.push_back(candidateBasisElement);
}
++it;
}
}
template<unsigned int dim>
const std::vector<Monomial<dim>> &MonomialBasis<dim>::getElements() const
{
return basis;
}
template<unsigned int dim>
MonomialBasis<dim> MonomialBasis<dim>::getDerivative(const Point<dim, unsigned int> differentialOrder) const
{
std::vector<Monomial<dim>> derivatives;
for (const auto &monomial : getElements())
{
derivatives.push_back(monomial.getDerivative(differentialOrder));
}
return MonomialBasis<dim>(derivatives);
}
template<unsigned int dim>
bool MonomialBasis<dim>::operator==(const MonomialBasis &other) const
{
return basis == other.basis;
}
//template<unsigned int dim>
//MonomialBasis<dim>::MonomialBasis(Point<dim, unsigned int> degrees, unsigned int convergenceOrder)
// : MonomialBasis(degrees.asArray(), convergenceOrder) {}
#endif //OPENFPM_PDATA_MONOMIALBASIS_H
//
// Created by tommaso on 29/03/19.
//
#ifndef OPENFPM_PDATA_SUPPORT_HPP
#define OPENFPM_PDATA_SUPPORT_HPP
#include <Space/Shape/Point.hpp>
#include <Vector/vector_dist.hpp>
template<unsigned int dim, typename T, typename Prop>
class Support
{
// This class is basically a wrapper around a point and a set of offsets.
// Offsets are stored as they are the data that is mostly required in DCPSE, so we
// pre-compute and store them, while the positions can be re-computed everytime is
// necessary (it should almost never be the case) (todo: check if this is required)
private:
const vector_dist<dim, T, Prop> &domain;
const size_t referencePointKey;
const std::vector<size_t> keys;
const std::vector<Point<dim, T>> offsets;
public:
Support() {};
Support(const vector_dist<dim, T, Prop> &domain, const size_t &referencePoint, const std::vector<size_t> &keys)
: domain(domain),
referencePointKey(referencePoint),
keys(keys),
offsets(computeOffsets(referencePoint, keys)) {}
Support(const Support<dim, T, Prop> &other);
size_t size();
const Point<dim, T> getReferencePoint() const;
const size_t getReferencePointKey() const;
const std::vector<size_t> &getKeys() const;
const std::vector<Point<dim, T>> &getOffsets() const;
private:
std::vector<Point<dim, T>>
computeOffsets(const size_t referencePoint, const std::vector<size_t> &keys);
};
template<unsigned int dim, typename T, typename Prop>
std::vector<Point<dim, T>>
Support<dim, T, Prop>::computeOffsets(const size_t referencePoint, const std::vector<size_t> &keys)
{
std::vector<Point<dim, T>> offsets;
for (auto &otherK : keys)
{
Point<dim, T> curOffset(domain.getPos(referencePoint));
curOffset -= domain.getPos(otherK);
offsets.push_back(curOffset);
}
return offsets;
}
template<unsigned int dim, typename T, typename Prop>
const Point<dim, T> Support<dim, T, Prop>::getReferencePoint() const
{
return Point<dim, T>(domain.getPos(referencePointKey));
}
template<unsigned int dim, typename T, typename Prop>
const std::vector<Point<dim, T>> &Support<dim, T, Prop>::getOffsets() const
{
return offsets;
}
template<unsigned int dim, typename T, typename Prop>
size_t Support<dim, T, Prop>::size()
{
return offsets.size();
}
template<unsigned int dim, typename T, typename Prop>
Support<dim, T, Prop>::Support(const Support<dim, T, Prop> &other)
: domain(other.domain),
referencePointKey(other.referencePointKey),
keys(other.keys),
offsets(other.offsets) {}
template<unsigned int dim, typename T, typename Prop>
const size_t Support<dim, T, Prop>::getReferencePointKey() const
{
return referencePointKey;
}
template<unsigned int dim, typename T, typename Prop>
const std::vector<size_t> &Support<dim, T, Prop>::getKeys() const
{
return keys;
}
#endif //OPENFPM_PDATA_SUPPORT_HPP
//
// Created by tommaso on 25/03/19.
//
#ifndef OPENFPM_PDATA_SUPPORTBUILDER_HPP
#define OPENFPM_PDATA_SUPPORTBUILDER_HPP
// This is to automatically get the right support (set of particles) for doing DCPSE on a given particle.
// todo: This could be enhanced to support an algorithm for choosing the support in a smart way (to lower condition num)
#include <Space/Shape/Point.hpp>
#include <Vector/vector_dist.hpp>
#include "Support.hpp"
template<unsigned int dim, typename T, typename Prop>
class SupportBuilder
{
private:
vector_dist<dim, T, Prop> &domain;
CellList<dim, T, Mem_fast<HeapMemory, local_index_>> cellList;
const Point<dim, unsigned int> differentialSignature;
public:
SupportBuilder(vector_dist<dim, T, Prop> &domain, Point<dim, unsigned int> differentialSignature, T rCut);
SupportBuilder(vector_dist<dim, T, Prop> &domain, unsigned int differentialSignature[dim], T rCut);
Support<dim, T, Prop> getSupport(vector_dist_iterator itPoint, unsigned int requiredSize);
private:
size_t getCellLinId(const grid_key_dx<dim> &cellKey);
size_t getNumElementsInCell(const grid_key_dx<dim> &cellKey);
size_t getNumElementsInSetOfCells(const std::set<grid_key_dx<dim>> &set);
void enlargeSetOfCellsUntilSize(std::set<grid_key_dx<dim>> &set, unsigned int requiredSize);
std::vector<size_t> getPointsInSetOfCells(std::set<grid_key_dx<dim>> set);
bool isCellKeyInBounds(grid_key_dx<dim> key);
};
// Method definitions below
template<unsigned int dim, typename T, typename Prop>
SupportBuilder<dim, T, Prop>::SupportBuilder(vector_dist<dim, T, Prop> &domain, const Point<dim, unsigned int> differentialSignature,
T rCut) : domain(domain), differentialSignature(differentialSignature)
{
cellList = domain.template getCellList<CellList<dim, T, Mem_fast<HeapMemory, local_index_>>>(rCut);
}
template<unsigned int dim, typename T, typename Prop>
Support<dim, T, Prop> SupportBuilder<dim, T, Prop>::getSupport(vector_dist_iterator itPoint, unsigned int requiredSize)
{
// Get spatial position from point iterator
vect_dist_key_dx p = itPoint.get();
Point<dim, T> pos = domain.getPos(p.getKey());
// Get cell containing current point and add it to the set of cell keys
grid_key_dx<dim> curCellKey = cellList.getCellGrid(pos); // Here get the key of the cell where the current point is
std::set<grid_key_dx<dim>> supportCells;
supportCells.insert(curCellKey);
// Make sure to consider a set of cells providing enough points for the support
enlargeSetOfCellsUntilSize(supportCells, requiredSize + 1); // NOTE: this +1 is because we then remove the point itself
// Now return all the points from the support into a vector
std::vector<size_t> supportKeys = getPointsInSetOfCells(supportCells);
std::remove(supportKeys.begin(), supportKeys.end(), p.getKey());
return Support<dim, T, Prop>(domain, p.getKey(), supportKeys);
}
template<unsigned int dim, typename T, typename Prop>
size_t SupportBuilder<dim, T, Prop>::getNumElementsInCell(const grid_key_dx<dim> &cellKey)
{
const size_t curCellId = getCellLinId(cellKey);
size_t numElements = cellList.getNelements(curCellId);
return numElements;
}
template<unsigned int dim, typename T, typename Prop>
size_t SupportBuilder<dim, T, Prop>::getNumElementsInSetOfCells(const std::set<grid_key_dx<dim>> &set)
{
size_t tot = 0;
for (const auto cell : set)
{
tot += getNumElementsInCell(cell);
}
return tot;
}
template<unsigned int dim, typename T, typename Prop>
void SupportBuilder<dim, T, Prop>::enlargeSetOfCellsUntilSize(std::set<grid_key_dx<dim>> &set, unsigned int requiredSize)
{
while (getNumElementsInSetOfCells(set) < requiredSize)
{
auto tmpSet = set;
for (const auto el : tmpSet)
{
for (unsigned int i = 0; i < dim; ++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);
}
}
}
}
}
template<unsigned int dim, typename T, typename Prop>
size_t SupportBuilder<dim, T, Prop>::getCellLinId(const grid_key_dx<dim> &cellKey)
{
mem_id id = cellList.getGrid().LinId(cellKey);
return static_cast<size_t>(id);
}
template<unsigned int dim, typename T, typename Prop>
std::vector<size_t> SupportBuilder<dim, T, Prop>::getPointsInSetOfCells(std::set<grid_key_dx<dim>> set)
{
std::vector<size_t> points;
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<dim, T> pos = domain.getPos(el);
points.push_back(el);
}
}
return points;
}
template<unsigned int dim, typename T, typename Prop>
SupportBuilder<dim, T, Prop>::SupportBuilder(vector_dist<dim, T, Prop> &domain, unsigned int *differentialSignature, T rCut)
: SupportBuilder(domain, Point<dim, unsigned int>(differentialSignature), rCut) {}
template<unsigned int dim, typename T, typename Prop>
bool SupportBuilder<dim, T, Prop>::isCellKeyInBounds(grid_key_dx<dim> key)
{
const size_t *cellGridSize = cellList.getGrid().getSize();
for (size_t i = 0; i < dim; ++i)
{
if (key.value(i) < 0 || key.value(i) >= cellGridSize[i])
{
return false;
}
}
return true;
}
#endif //OPENFPM_PDATA_SUPPORTBUILDER_HPP
//
// Created by tommaso on 21/03/19.
//
#ifndef OPENFPM_PDATA_VANDERMONDE_HPP
#define OPENFPM_PDATA_VANDERMONDE_HPP
#include "MonomialBasis.hpp"
#include "VandermondeRowBuilder.hpp"
#include "Support.hpp"
template<unsigned int dim, typename T, typename MatrixType>
class Vandermonde
{
private:
const Point<dim, T> point;
const std::vector<Point<dim, T>> offsets;
const MonomialBasis<dim> monomialBasis;
T eps;
public:
Vandermonde(const Point<dim, T> &point, const std::vector<Point<dim, T>> &neighbours,
const MonomialBasis<dim> &monomialBasis);
template<typename Prop>
Vandermonde(const Support<dim, T, Prop> &support,
const MonomialBasis<dim> &monomialBasis);
MatrixType &getMatrix(MatrixType &M);
T getEps();
private:
static std::vector<Point<dim, T>>
computeOffsets(const Point<dim, T> &point, const std::vector<Point<dim, T>> &neighbours);
void computeEps(T factor);
static T computeAbsSum(const Point<dim, T> &x);
void initialize();
};
template<unsigned int dim, typename T, typename MatrixType>
void Vandermonde<dim, T, MatrixType>::initialize()
{
// First check that the number of points given is enough for building the Vandermonde matrix
if (offsets.size() < monomialBasis.size())
{
ACTION_ON_ERROR(std::length_error("Not enough neighbour points passed for Vandermonde matrix construction!"));
}
// Compute eps for this point
computeEps(2);
}
template<unsigned int dim, typename T, typename MatrixType>
Vandermonde<dim, T, MatrixType>::Vandermonde(const Point<dim, T> &point, const std::vector<Point<dim, T>> &neighbours,
const MonomialBasis<dim> &monomialBasis)
: point(point),
offsets(computeOffsets(point, neighbours)),
monomialBasis(monomialBasis)
{
initialize();
}
template<unsigned int dim, typename T, typename MatrixType>
template<typename Prop>
Vandermonde<dim, T, MatrixType>::Vandermonde(const Support<dim, T, Prop> &support,
const MonomialBasis<dim> &monomialBasis)
: point(support.getReferencePoint()),
offsets(support.getOffsets()),
monomialBasis(monomialBasis)
{
initialize();
}
template<unsigned int dim, typename T, typename MatrixType>
std::vector<Point<dim, T>>
Vandermonde<dim, T, MatrixType>::computeOffsets(const Point<dim, T> &point,
const std::vector<Point<dim, T>> &neighbours)
{
std::vector<Point<dim, T>> offsets;
for (auto &other : neighbours)
{
Point<dim, T> curOffset(point);
curOffset -= other;
offsets.push_back(curOffset);
}
return offsets;
}
template<unsigned int dim, typename T, typename MatrixType>
MatrixType &Vandermonde<dim, T, MatrixType>::getMatrix(MatrixType &M)
{
// Build the Vandermonde matrix, row-by-row
VandermondeRowBuilder<dim, T> vrb(monomialBasis);
unsigned int row = 0;
for (auto &offset : offsets)
{
vrb.buildRow(M, row, offset, eps);
++row;
}
return M;
}
template<unsigned int dim, typename T, typename MatrixType>
void Vandermonde<dim, T, MatrixType>::computeEps(T factor)
{
T avgNeighbourSpacing = 0;
for (auto &offset : offsets)
{
avgNeighbourSpacing += computeAbsSum(offset);
}
avgNeighbourSpacing /= offsets.size();
eps = factor * avgNeighbourSpacing;
assert(eps != 0);
}
template<unsigned int dim, typename T, typename MatrixType>
T Vandermonde<dim, T, MatrixType>::computeAbsSum(const Point<dim, T> &x)
{
T absSum = 0;
for (unsigned int i = 0; i < dim; ++i)
{
absSum += fabs(x.value(i));
}
return absSum;
}
template<unsigned int dim, typename T, typename MatrixType>
T Vandermonde<dim, T, MatrixType>::getEps()
{
return eps;
}
#endif //OPENFPM_PDATA_VANDERMONDE_HPP
//
// Created by tommaso on 22/03/19.
//
#ifndef OPENFPM_PDATA_VANDERMONDEROW_HPP
#define OPENFPM_PDATA_VANDERMONDEROW_HPP
#include "MonomialBasis.hpp"
template <unsigned int dim, typename T>
class VandermondeRowBuilder
{
private:
const MonomialBasis<dim> monomialBasis;
public:
VandermondeRowBuilder(const MonomialBasis<dim> &monomialBasis) : monomialBasis(monomialBasis) {}
template <typename MatrixType>
void buildRow(MatrixType &M, unsigned int row, Point<dim, T> x, T eps);
};
template<unsigned int dim, typename T>
template <typename MatrixType>
void VandermondeRowBuilder<dim, T>::buildRow(MatrixType &M, unsigned int row, Point<dim, T> x, T eps)
{
unsigned int col = 0;
for (auto& basisElement : monomialBasis.getElements())
{
Monomial<dim> m = monomialBasis.getElement(col);
M(row, col) = m.evaluate(x);
M(row, col) /= pow(eps, m.order());
++col;
}
}
#endif //OPENFPM_PDATA_VANDERMONDEROW_HPP
//
// Created by tommaso on 28/03/19.
//
#define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>
#include <DCPSE/DcpseRhs.hpp>
#include "../../openfpm_numerics/src/DMatrix/EMatrix.hpp"
BOOST_AUTO_TEST_SUITE(DcpseRhs_tests)
// If EIGEN is not present, EMatrix is not available and we don't need to build this test
#ifdef HAVE_EIGEN
BOOST_AUTO_TEST_CASE(DcpseRhs_dx_test)
{
unsigned int r = 2;
Point<2, unsigned int> m({1,0});
MonomialBasis<2> mb(m.asArray(), r);
std::cout << mb << std::endl;
EMatrix<double, Eigen::Dynamic, 1> b(mb.size());
DcpseRhs<2> rhs(mb, m);
rhs.getVector<double>(b);
// for (unsigned int i = 0; i < mb.size(); ++i)
// {
// std::cout << b(i) << std::endl;
// }
// Validation
const MonomialBasis<2> &Dmb = mb.getDerivative(m);
std::cout << Dmb << std::endl;
auto p0 = Point<2, double>({0,0});
BOOST_REQUIRE_CLOSE(b(0), -Dmb.getElement(0).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(1), -Dmb.getElement(1).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(2), -Dmb.getElement(2).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(3), -Dmb.getElement(3).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(4), -Dmb.getElement(4).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(5), -Dmb.getElement(5).evaluate(p0), 1e-16);
}
BOOST_AUTO_TEST_CASE(DcpseRhs_dxdy_test)
{
unsigned int r = 2;
Point<2, unsigned int> m({1,1});
MonomialBasis<2> mb(m.asArray(), r);
std::cout << mb << std::endl;
EMatrix<double, Eigen::Dynamic, Eigen::Dynamic> b(mb.size(), 1);
DcpseRhs<2> rhs(mb, m);
rhs.getVector<double>(b);
// for (unsigned int i = 0; i < mb.size(); ++i)
// {
// std::cout << b(i) << std::endl;
// }
// Validation
const MonomialBasis<2> &Dmb = mb.getDerivative(m);
std::cout << Dmb << std::endl;
auto p0 = Point<2, double>({0,0});
BOOST_REQUIRE_CLOSE(b(0), Dmb.getElement(0).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(1), Dmb.getElement(1).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(2), Dmb.getElement(2).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(3), Dmb.getElement(3).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(4), Dmb.getElement(4).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(5), Dmb.getElement(5).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(6), Dmb.getElement(6).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(7), Dmb.getElement(7).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(8), Dmb.getElement(8).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(9), Dmb.getElement(9).evaluate(p0), 1e-16);
}
BOOST_AUTO_TEST_CASE(DcpseRhs_laplacian_test)
{
unsigned int r = 2;
Point<2, unsigned int> m({2,2});
MonomialBasis<2> mb(m.asArray(), r);
std::cout << mb << std::endl;
EMatrix<double, Eigen::Dynamic, Eigen::Dynamic> b(mb.size(), 1);
DcpseRhs<2> rhs(mb, m);
rhs.getVector<double>(b);
// Validation
const MonomialBasis<2> &Dmb = mb.getDerivative(m);
std::cout << Dmb << std::endl;
auto p0 = Point<2, double>({0,0});
BOOST_REQUIRE_CLOSE(b(0), Dmb.getElement(0).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(1), Dmb.getElement(1).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(2), Dmb.getElement(2).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(3), Dmb.getElement(3).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(4), Dmb.getElement(4).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(5), Dmb.getElement(5).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(6), Dmb.getElement(6).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(7), Dmb.getElement(7).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(8), Dmb.getElement(8).evaluate(p0), 1e-16);
BOOST_REQUIRE_CLOSE(b(9), Dmb.getElement(9).evaluate(p0), 1e-16);
}
#endif // HAVE_EIGEN
BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file
This diff is collapsed.
//
// Created by tommaso on 21/03/19.
//
#define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>
#include <Space/Shape/Point.hpp>
#include <DCPSE/MonomialBasis.hpp>
BOOST_AUTO_TEST_SUITE(MonomialBasis_tests)
BOOST_AUTO_TEST_CASE(Monomial_3D_test)
{
// Testing constructor from a Point
Point<3, unsigned int> p({1, 2, 3});
Monomial<3> mbe(p);
BOOST_REQUIRE_EQUAL(mbe.order(), 6);
BOOST_REQUIRE_EQUAL(mbe.getExponent(0), 1);
BOOST_REQUIRE_EQUAL(mbe.getExponent(1), 2);
BOOST_REQUIRE_EQUAL(mbe.getExponent(2), 3);
// Testing copy constructor
Monomial<3> mbe2(mbe);
BOOST_REQUIRE_EQUAL(mbe2.order(), 6);
BOOST_REQUIRE_EQUAL(mbe2.getExponent(0), 1);
BOOST_REQUIRE_EQUAL(mbe2.getExponent(1), 2);
BOOST_REQUIRE_EQUAL(mbe2.getExponent(2), 3);
// Testing copy assignment operator
Monomial<3> mbe3;
mbe3 = mbe;
BOOST_REQUIRE_EQUAL(mbe3.order(), 6);
BOOST_REQUIRE_EQUAL(mbe3.getExponent(0), 1);
BOOST_REQUIRE_EQUAL(mbe3.getExponent(1), 2);
BOOST_REQUIRE_EQUAL(mbe3.getExponent(2), 3);
// Testing default constructor
Monomial<3> mbe0;
BOOST_REQUIRE_EQUAL(mbe0.order(), 0);
BOOST_REQUIRE_EQUAL(mbe0.getExponent(0), 0);
BOOST_REQUIRE_EQUAL(mbe0.getExponent(1), 0);
BOOST_REQUIRE_EQUAL(mbe0.getExponent(2), 0);
// Testing one-by-one set values
mbe0.setExponent(0, 5);
mbe0.setExponent(1, 5);
mbe0.setExponent(2, 5);
BOOST_REQUIRE_EQUAL(mbe0.order(), 15);
BOOST_REQUIRE_EQUAL(mbe0.getExponent(0), 5);
BOOST_REQUIRE_EQUAL(mbe0.getExponent(1), 5);
BOOST_REQUIRE_EQUAL(mbe0.getExponent(2), 5);
}
BOOST_AUTO_TEST_CASE(Monomial_evaluate_test)
{
// Test with uniform exponent and different coordinates
{
Monomial<3> monomial(Point<3, unsigned int>({1, 1, 1}));
double res = monomial.evaluate(Point<3, double>({5, 3, 2}));
double expected = 5 * 3 * 2;
BOOST_REQUIRE_CLOSE(res, expected, 1e-6);
}
// Test with different exponents and uniform coordinates
{
Monomial<3> monomial(Point<3, unsigned int>({5, 3, 2}));
double res = monomial.evaluate(Point<3, double>({2, 2, 2}));
double expected = pow(2,5) * pow(2,3) * pow(2,2);
BOOST_REQUIRE_CLOSE(res, expected, 1e-6);
}
}
BOOST_AUTO_TEST_CASE(Monomial_derivative_test)
{
// Test differentiation along 3 dimensions independently
{
Monomial<3> monomial(Point<3, unsigned int>({6, 6, 6}));
Monomial<3> derivative = monomial.getDerivative(Point<3, double>({3, 0, 0}));
Monomial<3> expected(Point<3, unsigned int>({3, 6, 6}), 6*5*4);
BOOST_REQUIRE_EQUAL(derivative, expected);
}
{
Monomial<3> monomial(Point<3, unsigned int>({6, 6, 6}));
Monomial<3> derivative = monomial.getDerivative(Point<3, double>({0, 3, 0}));
Monomial<3> expected(Point<3, unsigned int>({6, 3, 6}), 6*5*4);
BOOST_REQUIRE_EQUAL(derivative, expected);
}
{
Monomial<3> monomial(Point<3, unsigned int>({6, 6, 6}));
Monomial<3> derivative = monomial.getDerivative(Point<3, double>({0, 0, 3}));
Monomial<3> expected(Point<3, unsigned int>({6, 6, 3}), 6*5*4);
BOOST_REQUIRE_EQUAL(derivative, expected);
}
// Test for correctly differentiating beyond degree
{
Monomial<3> monomial(Point<3, unsigned int>({6, 6, 6}));
Monomial<3> derivative = monomial.getDerivative(Point<3, double>({7, 0, 0}));
Monomial<3> expected(Point<3, unsigned int>({0, 6, 6}), 0);
BOOST_REQUIRE_EQUAL(derivative, expected);
}
{
Monomial<3> monomial(Point<3, unsigned int>({6, 6, 6}));
Monomial<3> derivative = monomial.getDerivative(Point<3, double>({0, 7, 0}));
Monomial<3> expected(Point<3, unsigned int>({6, 0, 6}), 0);
BOOST_REQUIRE_EQUAL(derivative, expected);
}
{
Monomial<3> monomial(Point<3, unsigned int>({6, 6, 6}));
Monomial<3> derivative = monomial.getDerivative(Point<3, double>({0, 0, 7}));
Monomial<3> expected(Point<3, unsigned int>({6, 6, 0}), 0);
BOOST_REQUIRE_EQUAL(derivative, expected);
}
}
BOOST_AUTO_TEST_CASE(MonomialBasis_2D_test)
{
{
MonomialBasis<2> mb({1, 0}, 2);
// Check that the basis contains all the elements it is supposed to contain
std::vector<Monomial<2>> control;
control.emplace_back(Point<2, unsigned int>({0, 0}));
control.emplace_back(Point<2, unsigned int>({1, 0}));
control.emplace_back(Point<2, unsigned int>({0, 1}));
control.emplace_back(Point<2, unsigned int>({1, 1}));
control.emplace_back(Point<2, unsigned int>({2, 0}));
control.emplace_back(Point<2, unsigned int>({0, 2}));
auto first = mb.getElements().begin();
auto last = mb.getElements().end();
int foundElements = 0;
for (const auto &item : control)
{
auto pos = std::find(first, last, item);
foundElements += (pos != last);
}
BOOST_REQUIRE_EQUAL(foundElements, control.size());
}
{
MonomialBasis<2> mb({1, 1}, 2);
// std::cout << mb << std::endl;
// std::cout << mb.size() << std::endl;
// Check that the basis contains all the elements it is supposed to contain
std::vector<Monomial<2>> control;
control.emplace_back(Point<2, unsigned int>({0, 0}));
control.emplace_back(Point<2, unsigned int>({1, 0}));
control.emplace_back(Point<2, unsigned int>({0, 1}));
control.emplace_back(Point<2, unsigned int>({1, 1}));
control.emplace_back(Point<2, unsigned int>({2, 0}));
control.emplace_back(Point<2, unsigned int>({0, 2}));
control.emplace_back(Point<2, unsigned int>({2, 1}));
control.emplace_back(Point<2, unsigned int>({1, 2}));
control.emplace_back(Point<2, unsigned int>({3, 0}));
control.emplace_back(Point<2, unsigned int>({0, 3}));
BOOST_REQUIRE_EQUAL(mb.size(), control.size());
auto first = mb.getElements().begin();
auto last = mb.getElements().end();
int foundElements = 0;
for (const auto &item : control)
{
auto pos = std::find(first, last, item);
foundElements += (pos != last);
}
BOOST_REQUIRE_EQUAL(foundElements, control.size());
}
}
BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file
//
// Created by tommaso on 26/03/19.
//
#define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>
#include "Vector/vector_dist.hpp"
#include <Space/Shape/Point.hpp>
#include "DCPSE/SupportBuilder.hpp"
BOOST_AUTO_TEST_SUITE(Support_tests)
BOOST_AUTO_TEST_CASE(SupportBuilder_2D_1_0_2spacing_test)
{
// Here build some easy domain and get some points around a given one
size_t edgeSemiSize = 100;
const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
Box<2, double> box({-1.1, -1.1}, {1.1, 1.1});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing[2];
spacing[0] = 1.0 / (sz[0] - 1);
spacing[1] = 1.0 / (sz[1] - 1);
Ghost<2, double> ghost(0.1);
vector_dist<2, double, aggregate<double>> domain(0, box, bc, ghost);
auto it = domain.getGridIterator(sz);
size_t pointId = 0;
size_t counter = 0;
double minNormOne = 999;
while (it.isNext())
{
domain.add();
auto key = it.get();
mem_id k0 = key.get(0) - edgeSemiSize;
double x = k0 * spacing[0];
domain.getLastPos()[0] = x;
mem_id k1 = key.get(1) - edgeSemiSize;
double y = k1 * spacing[1];
domain.getLastPos()[1] = y;
domain.template getLastProp<0>() = 0.0;
++counter;
++it;
}
// Now get iterator to point of interest
auto itPoint = domain.getIterator();
size_t foo = (sqrt(counter) + counter) / 2;
for (int i = 0; i < foo; ++i)
{
++itPoint;
}
// Get spatial position from point iterator
vect_dist_key_dx p = itPoint.get();
const auto pos = domain.getPos(p.getKey());
std::cout << "p=(" << pos[0] << "," << pos[1] << ")" << std::endl;
// BOOST_REQUIRE_CLOSE(pos[0], 0, 1e-16);
// BOOST_REQUIRE_CLOSE(pos[1], 0, 1e-16);
// Now that domain is built and populated, let's test SupportBuilder
// We use (0,0) as initial point
SupportBuilder<2, double, aggregate<double>> supportBuilder(domain, {1,0}, 2*spacing[0]);
auto support = supportBuilder.getSupport(itPoint, 6);
// for (const auto &off : support.getOffsets())
// {
// std::cout << off.toString() << std::endl;
// }
BOOST_REQUIRE_GE(support.size(), 6);
}
BOOST_AUTO_TEST_CASE(SupportBuilder_2D_2_2_2spacing_test)
{
// Here build some easy domain and get some points around a given one
size_t edgeSemiSize = 25;
const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
Box<2, double> box({-1.2, -1.2}, {1.2, 1.2});
// Box<2, double> innerDomain({-1.0, -1.0}, {1.0, 1.0});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing[2];
spacing[0] = 1.0 / (sz[0] - 1);
spacing[1] = 1.0 / (sz[1] - 1);
Ghost<2, double> ghost(3.0 * spacing[0]);
vector_dist<2, double, aggregate<double>> domain(0, box, bc, ghost);
auto it = domain.getGridIterator(sz);
size_t pointId = 0;
size_t counter = 0;
while (it.isNext())
{
domain.add();
auto key = it.get();
mem_id k0 = key.get(0) - edgeSemiSize;
double x = k0 * spacing[0];
domain.getLastPos()[0] = x;
mem_id k1 = key.get(1) - edgeSemiSize;
double y = k1 * spacing[1];
domain.getLastPos()[1] = y;
domain.template getLastProp<0>() = 0.0;
if (abs(domain.getLastPos()[0]) + abs(domain.getLastPos()[1]) < 1e-16) // i.e. if we are at (0,0)
{
pointId = counter;
}
++counter;
++it;
}
// Now get iterator to point of interest
auto itPoint = domain.getDomainIterator();
size_t foo = (sqrt(counter) + counter) / 2;
for (int i = 0; i < foo; ++i)
{
++itPoint;
}
// Get spatial position from point iterator
vect_dist_key_dx p = itPoint.get();
const auto pos = domain.getPos(p.getKey());
std::cout << "p=(" << pos[0] << "," << pos[1] << ")" << std::endl;
// BOOST_REQUIRE_CLOSE(pos[0], 0, 1e-16);
// BOOST_REQUIRE_CLOSE(pos[1], 0, 1e-16);
// Now that domain is built and populated, let's test SupportBuilder
// We use (0,0) as initial point
SupportBuilder<2, double, aggregate<double>> supportBuilder(domain, {2,2}, 2*spacing[0]);
auto supportPoints = supportBuilder.getSupport(itPoint, 20);
// for (const auto &k : supportPoints)
// {
// Point<2, double> pt = domain.getPos(k);
// std::cout << pt.toString() << std::endl;
// }
BOOST_REQUIRE_GE(supportPoints.size(), 20);
}
// BOOST_AUTO_TEST_CASE(Support_CopyConstructor_test)
// {
//
// }
BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file
//
// Created by tommaso on 22/03/19.
//
#define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>
#include <DCPSE/MonomialBasis.hpp>
#include <DCPSE/VandermondeRowBuilder.hpp>
#include <DCPSE/Vandermonde.hpp>
#include "../../openfpm_numerics/src/DMatrix/EMatrix.hpp"
BOOST_AUTO_TEST_SUITE(Vandermonde_tests)
// If EIGEN is not present, EMatrix is not available and we don't need to build this test
#ifdef HAVE_EIGEN
BOOST_AUTO_TEST_CASE(VandermondeRowBuilder_AllOnes_test)
{
MonomialBasis<2> mb({1, 0}, 2);
EMatrix<double, Eigen::Dynamic, Eigen::Dynamic> row(1, mb.size());
Point<2, double> x({2, 2});
double eps = 2;
VandermondeRowBuilder<2, double> vrb(mb);
vrb.buildRow(row, 0, x, eps);
// For the way the row has been constructed, it should be composed of only 1s
bool isRowAllOnes = true;
for (int i = 0; i < mb.size(); ++i)
{
isRowAllOnes = isRowAllOnes && (row(0, i) == 1);
}
BOOST_REQUIRE(isRowAllOnes);
}
BOOST_AUTO_TEST_CASE(VandermondeRowBuilder_OneZero_test)
{
MonomialBasis<2> mb({1, 0}, 2);
EMatrix<double, Eigen::Dynamic, Eigen::Dynamic> row(1, mb.size());
Point<2, double> x({1, 0});
double eps = 1;
VandermondeRowBuilder<2, double> vrb(mb);
vrb.buildRow(row, 0, x, eps);
// For the way the row has been constructed, it should be composed of only 1s
bool areValuesOk = true;
for (int i = 0; i < mb.size(); ++i)
{
bool isThereY = mb.getElement(i).getExponent(1) > 0;
bool curCheck = (row(0, i) == !isThereY);
areValuesOk = areValuesOk && curCheck;
}
BOOST_REQUIRE(areValuesOk);
}