...
 
Commits (2)
  • bianucci's avatar
    Fix symbolic derivative computation in Monomial and restructure code to make... · 0030cd72
    bianucci authored
    Fix symbolic derivative computation in Monomial and restructure code to make it suitable to final integration into DCPSE.
    
    1) The derivative computation in Monomial was totally wrong as it did not account for the scalar factor. Now it has been fixed and the tests modified accordingly.
    2) The interface of Support has been changed to make it more suitable for the expected usage in the integrated DCPSE computation. Interfaces of the other classes interacting with Support has been changed accordingly.
    0030cd72
  • bianucci's avatar
    Add DC-PSE feature. · 12121488
    bianucci authored
    12121488
......@@ -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 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/Support.hpp DCPSE/tests/Support_unit_tests.cpp DCPSE/DcpseRhs.hpp DCPSE/tests/DcpseRhs_unit_tests.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)
if ( CMAKE_COMPILER_IS_GNUCC )
target_compile_options(pdata PRIVATE "-Wno-deprecated-declarations")
......
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
......@@ -8,7 +8,7 @@
#include "../../openfpm_numerics/src/DMatrix/EMatrix.hpp"
#include "MonomialBasis.hpp"
template<unsigned int dim, typename T, typename MatrixType>
template<unsigned int dim>
class DcpseRhs
{
private:
......@@ -18,29 +18,30 @@ private:
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, typename T, typename MatrixType>
DcpseRhs<dim, T, MatrixType>::DcpseRhs(const MonomialBasis<dim> &monomialBasis,
const Point<dim, unsigned int> &differentialSignature)
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
} else
{
sign = -1;
}
}
template<unsigned int dim, typename T, typename MatrixType>
MatrixType &DcpseRhs<dim, T, MatrixType>::getVector(MatrixType &b)
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);
......@@ -48,7 +49,7 @@ MatrixType &DcpseRhs<dim, T, MatrixType>::getVector(MatrixType &b)
for (unsigned int i = 0; i < derivatives.size(); ++i)
{
const Monomial<dim> dm = derivatives.getElement(i);
b(i,0) = sign * dm.evaluate(Point<2, double>({0,0}));
b(i, 0) = sign * dm.evaluate(Point<dim, T>(0));
}
return b;
}
......
......@@ -17,17 +17,18 @@ class Monomial
private:
unsigned int sum = 0;
Point<dim, unsigned int> exponents;
unsigned int scalar = 1;
public:
Monomial();
explicit Monomial(const Point<dim, unsigned int> &other);
explicit Monomial(const Point<dim, unsigned int> &other, unsigned int scalar = 1);
explicit Monomial(const Point<dim, long int> &other);
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(other.exponents) {}
Monomial(const Monomial<dim> &other);
Monomial<dim> &operator=(const Monomial<dim> &other);
......@@ -51,7 +52,7 @@ public:
friend std::basic_ostream<charT, traits> &
operator<<(std::basic_ostream<charT, traits> &lhs, Monomial<dim> const &rhs)
{
return lhs << rhs.exponents.toString();
return lhs << rhs.scalar << " : " << rhs.exponents.toString();
}
private:
......@@ -68,14 +69,15 @@ Monomial<dim>::Monomial()
}
template<unsigned int dim>
Monomial<dim>::Monomial(const Point<dim, unsigned int> &other)
: exponents(other)
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)
Monomial<dim>::Monomial(const Point<dim, long int> &other, unsigned int scalar)
: scalar(scalar)
{
for (size_t i = 0; i < other.nvals; ++i)
{
......@@ -88,11 +90,18 @@ 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;
}
......@@ -130,14 +139,14 @@ template<unsigned int dim>
bool Monomial<dim>::operator==
(const Monomial<dim> &other) const
{
return exponents == other.exponents;
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 = 1;
T res = scalar;
for (unsigned int i = 0; i < dim; ++i)
{
res *= pow(x.value(i), getExponent(i));
......@@ -148,14 +157,19 @@ T Monomial<dim>::evaluate(const Point<dim, T> x) const
template<unsigned int dim>
Monomial<dim> Monomial<dim>::getDerivative(const Point<dim, unsigned int> differentialOrder) const
{
Point<dim, unsigned int> res(exponents);
unsigned int s = scalar;
Point<dim, unsigned int> e(exponents);
for (unsigned int i = 0; i < dim; ++i)
{
res.get(i) = static_cast<unsigned int>(
std::max(static_cast<int>(res.value(i)) - static_cast<int>(differentialOrder.value(i)), 0)
);
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(res);
return Monomial(e, s);
}
template<unsigned int dim>
......
//
// Created by tommaso on 25/03/19.
// Created by tommaso on 29/03/19.
//
#ifndef OPENFPM_PDATA_SUPPORT_HPP
#define OPENFPM_PDATA_SUPPORT_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>
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:
vector_dist<dim, T, Prop> &domain;
CellList<dim, T, Mem_fast<HeapMemory, local_index_>> cellList;
const Point<dim, unsigned int> differentialOrder;
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(vector_dist<dim, T, Prop> &domain, Point<dim, unsigned int> differentialOrder, T rCut);
Support() {};
Support(vector_dist<dim, T, Prop> &domain, unsigned int differentialOrder[dim], T rCut);
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)) {}
std::vector<Point<dim, T>>
getSupport(vector_dist_iterator itPoint, unsigned int requiredSize);
Support(const Support<dim, T, Prop> &other);
private:
size_t getCellLinId(const grid_key_dx<dim> &cellKey);
size_t size();
size_t getNumElementsInCell(const grid_key_dx<dim> &cellKey);
const Point<dim, T> getReferencePoint() const;
size_t getNumElementsInSetOfCells(const std::set<grid_key_dx<dim>> &set);
const size_t getReferencePointKey() const;
void enlargeSetOfCellsUntilSize(std::set<grid_key_dx<dim>> &set, unsigned int requiredSize);
const std::vector<size_t> &getKeys() const;
std::vector<Point<dim, T>> getPointsInSetOfCells(std::set<grid_key_dx<dim>> set);
const std::vector<Point<dim, T>> &getOffsets() const;
bool isCellKeyInBounds(grid_key_dx<dim> key);
private:
std::vector<Point<dim, T>>
computeOffsets(const size_t referencePoint, const std::vector<size_t> &keys);
};
// Method definitions below
template<unsigned int dim, typename T, typename Prop>
Support<dim, T, Prop>::Support(vector_dist<dim, T, Prop> &domain, const Point<dim, unsigned int> differentialOrder,
T rCut) : domain(domain), differentialOrder(differentialOrder)
std::vector<Point<dim, T>>
Support<dim, T, Prop>::computeOffsets(const size_t referencePoint, const std::vector<size_t> &keys)
{
cellList = domain.template getCellList<CellList<dim, T, Mem_fast<HeapMemory, local_index_>>>(rCut);
}
template<unsigned int dim, typename T, typename Prop>
std::vector<Point<dim, T>> Support<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);
// Now return all the points from the support into a vector
return getPointsInSetOfCells(supportCells);
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>
size_t Support<dim, T, Prop>::getNumElementsInCell(const grid_key_dx<dim> &cellKey)
const Point<dim, T> Support<dim, T, Prop>::getReferencePoint() const
{
const size_t curCellId = getCellLinId(cellKey);
size_t numElements = cellList.getNelements(curCellId);
return numElements;
return domain.getPos(referencePointKey);
}
template<unsigned int dim, typename T, typename Prop>
size_t Support<dim, T, Prop>::getNumElementsInSetOfCells(const std::set<grid_key_dx<dim>> &set)
const std::vector<Point<dim, T>> &Support<dim, T, Prop>::getOffsets() const
{
size_t tot = 0;
for (const auto cell : set)
{
tot += getNumElementsInCell(cell);
}
return tot;
return offsets;
}
template<unsigned int dim, typename T, typename Prop>
void Support<dim, T, Prop>::enlargeSetOfCellsUntilSize(std::set<grid_key_dx<dim>> &set, unsigned int requiredSize)
size_t Support<dim, T, Prop>::size()
{
while (getNumElementsInSetOfCells(set) < requiredSize)
{
auto tmpSet = set;
for (const auto el : tmpSet)
{
for (unsigned int i = 0; i < dim; ++i)
{
if (differentialOrder.value(i) != 0)
{
// TODO: here check not to go out of bound with the grid keys!
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);
}
}
}
}
}
return offsets.size();
}
template<unsigned int dim, typename T, typename Prop>
size_t Support<dim, T, Prop>::getCellLinId(const grid_key_dx<dim> &cellKey)
{
mem_id id = cellList.getGrid().LinId(cellKey);
return static_cast<size_t>(id);
}
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>
std::vector<Point<dim, T>> Support<dim, T, Prop>::getPointsInSetOfCells(std::set<grid_key_dx<dim>> set)
const size_t Support<dim, T, Prop>::getReferencePointKey() const
{
std::vector<Point<dim, 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)
{
auto el = cellList.get(cellLinId, k);
Point<dim, T> pos = domain.getPos(el);
points.push_back(pos);
}
}
return points;
return referencePointKey;
}
template<unsigned int dim, typename T, typename Prop>
Support<dim, T, Prop>::Support(vector_dist<dim, T, Prop> &domain, unsigned int *differentialOrder, T rCut)
: Support(domain, Point<dim, unsigned int>(differentialOrder), rCut) {}
template<unsigned int dim, typename T, typename Prop>
bool Support<dim, T, Prop>::isCellKeyInBounds(grid_key_dx<dim> key)
const std::vector<size_t> &Support<dim, T, Prop>::getKeys() const
{
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;
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
......@@ -8,13 +8,14 @@
#include "../../openfpm_numerics/src/DMatrix/EMatrix.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;
std::vector<Point<dim, T>> offsets;
const std::vector<Point<dim, T>> offsets;
const MonomialBasis<dim> monomialBasis;
T eps;
......@@ -22,41 +23,72 @@ 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:
void computeOffsets(const std::vector<Point<dim, T>> &neighbours);
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(Point<dim, T> &x);
static T computeAbsSum(const Point<dim, T> &x);
void initialize();
};
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),
monomialBasis(monomialBasis)
void Vandermonde<dim, T, MatrixType>::initialize()
{
// First check that the number of points given is enough for building the Vandermonde matrix
if (neighbours.size() < monomialBasis.size())
if (offsets.size() < monomialBasis.size())
{
ACTION_ON_ERROR(std::length_error("Not enough neighbour points passed for Vandermonde matrix construction!"));
}
// Compute the offsets from point to all neighbours (and store them)
computeOffsets(neighbours);
// Compute eps for this point
computeEps(2);
}
template<unsigned int dim, typename T, typename MatrixType>
void Vandermonde<dim, T, MatrixType>::computeOffsets(const std::vector<Point<dim, T>> &neighbours)
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>
......@@ -86,14 +118,20 @@ void Vandermonde<dim, T, MatrixType>::computeEps(T factor)
}
template<unsigned int dim, typename T, typename MatrixType>
T Vandermonde<dim, T, MatrixType>::computeAbsSum(Point<dim, T> &x)
T Vandermonde<dim, T, MatrixType>::computeAbsSum(const Point<dim, T> &x)
{
T absSum = 0;
for (const auto elem : x.asArray())
for (unsigned int i = 0; i < dim; ++i)
{
absSum += abs(elem);
absSum += abs(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
......@@ -12,39 +12,27 @@ template <unsigned int dim, typename T>
class VandermondeRowBuilder
{
private:
MonomialBasis<dim> monomialBasis;
const MonomialBasis<dim> monomialBasis;
public:
VandermondeRowBuilder(const MonomialBasis<dim> &monomialBasis) : monomialBasis(monomialBasis) {}
void buildRow(EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> &M, unsigned int row, Point<dim, T> x, T eps);
private:
T computePower(Monomial<dim> mbe, Point<dim, T> x);
template <typename MatrixType>
void buildRow(MatrixType &M, unsigned int row, Point<dim, T> x, T eps);
};
template<unsigned int dim, typename T>
void VandermondeRowBuilder<dim, T>::buildRow(EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> &M, unsigned int row, Point<dim, T> x, T eps)
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())
{
auto & mbe = monomialBasis.getElement(col);
M(row, col) = computePower(mbe, x);
M(row, col) /= pow(eps, mbe.order());
Monomial<dim> m = monomialBasis.getElement(col);
M(row, col) = m.evaluate(x);
M(row, col) /= pow(eps, m.order());
++col;
}
}
template<unsigned int dim, typename T>
T VandermondeRowBuilder<dim, T>::computePower(Monomial<dim> mbe, Point<dim, T> x)
{
T res = 1;
for (int i = 0; i < dim; ++i)
{
res *= pow(x.value(i), mbe.getExponent(i));
}
return res;
}
#endif //OPENFPM_PDATA_VANDERMONDEROW_HPP
......@@ -19,8 +19,8 @@ BOOST_AUTO_TEST_SUITE(DcpseRhs_tests)
std::cout << mb << std::endl;
EMatrix<double, Eigen::Dynamic, 1> b(mb.size());
DcpseRhs<2, double, EMatrix<double, Eigen::Dynamic, 1>> rhs(mb, m);
rhs.getVector(b);
DcpseRhs<2> rhs(mb, m);
rhs.getVector<double>(b);
// for (unsigned int i = 0; i < mb.size(); ++i)
// {
// std::cout << b(i) << std::endl;
......@@ -28,6 +28,7 @@ BOOST_AUTO_TEST_SUITE(DcpseRhs_tests)
// 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);
......@@ -45,8 +46,8 @@ BOOST_AUTO_TEST_SUITE(DcpseRhs_tests)
std::cout << mb << std::endl;
EMatrix<double, Eigen::Dynamic, Eigen::Dynamic> b(mb.size(), 1);
DcpseRhs<2, double, EMatrix<double, Eigen::Dynamic, Eigen::Dynamic>> rhs(mb, m);
rhs.getVector(b);
DcpseRhs<2> rhs(mb, m);
rhs.getVector<double>(b);
// for (unsigned int i = 0; i < mb.size(); ++i)
// {
// std::cout << b(i) << std::endl;
......@@ -54,6 +55,34 @@ BOOST_AUTO_TEST_SUITE(DcpseRhs_tests)
// 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);
......
......@@ -77,38 +77,38 @@ BOOST_AUTO_TEST_SUITE(MonomialBasis_tests)
{
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}));
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}));
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}));
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}));
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}));
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}));
Monomial<3> expected(Point<3, unsigned int>({6, 6, 0}), 0);
BOOST_REQUIRE_EQUAL(derivative, expected);
}
}
......
......@@ -7,17 +7,16 @@
#include <boost/test/unit_test.hpp>
#include "Vector/vector_dist.hpp"
#include <Space/Shape/Point.hpp>
#include "DCPSE/Support.hpp"
#include "DCPSE/SupportBuilder.hpp"
BOOST_AUTO_TEST_SUITE(Support_tests)
BOOST_AUTO_TEST_CASE(Support_2D_1_0_2spacing_test)
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});
// 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);
......@@ -58,18 +57,18 @@ BOOST_AUTO_TEST_SUITE(Support_tests)
// 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 Support
// Now that domain is built and populated, let's test SupportBuilder
// We use (0,0) as initial point
Support<2, double, aggregate<double>> support(domain, {1,0}, 2*spacing[0]);
auto supportPoints = support.getSupport(itPoint, 6);
// for (const auto &pt : supportPoints)
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 << pt.toString() << std::endl;
// std::cout << off.toString() << std::endl;
// }
BOOST_REQUIRE_GE(supportPoints.size(), 6);
BOOST_REQUIRE_GE(support.size(), 6);
}
BOOST_AUTO_TEST_CASE(Support_2D_2_2_2spacing_test)
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;
......@@ -118,15 +117,21 @@ BOOST_AUTO_TEST_SUITE(Support_tests)
// 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 Support
// Now that domain is built and populated, let's test SupportBuilder
// We use (0,0) as initial point
Support<2, double, aggregate<double>> support(domain, {2,2}, 2*spacing[0]);
auto supportPoints = support.getSupport(itPoint, 20);
// for (const auto &pt : supportPoints)
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
......@@ -82,6 +82,7 @@ BOOST_AUTO_TEST_SUITE(Vandermonde_tests)
neighbours.push_back({-1,1});
neighbours.push_back({0,1});
neighbours.push_back({0,-1});
// ...and get the matrix V
Vandermonde<2, double, EMatrix<double, Eigen::Dynamic, Eigen::Dynamic>> vandermonde(x, neighbours, mb);
vandermonde.getMatrix(V);
......