Commit ec577798 authored by bianucci's avatar bianucci

Add RHS constuction for DCPSE and fix build issue when EIGEN package not properly found at first.

The DcpseRhs class provides machinery for getting the RHS to be used in the linear system associated with DCPSE computation. Analogously to Vandermonde, the user has to initialize the class with the MonomialBasis and differential signature in use and then pass an existing vector object to be filled to the getVector() method.
The fix in CMakeLists.txt is to allow MacOs and CentOS builds to be successful: in these cases the EIGEN package may not be found at first and this case is now supported.
parent 15adac4e
......@@ -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)
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)
if ( CMAKE_COMPILER_IS_GNUCC )
target_compile_options(pdata PRIVATE "-Wno-deprecated-declarations")
......@@ -49,7 +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})
target_include_directories (pdata PUBLIC ${EIGEN3_INCLUDE_DIR})
if(EIGEN3_FOUND)
target_include_directories (pdata PUBLIC ${EIGEN3_INCLUDE_DIR})
endif()
target_include_directories (pdata PUBLIC ${Boost_INCLUDE_DIRS})
......
//
// Created by tommaso on 28/03/19.
//
#ifndef OPENFPM_PDATA_DCPSERHS_HPP
#define OPENFPM_PDATA_DCPSERHS_HPP
#include "../../openfpm_numerics/src/DMatrix/EMatrix.hpp"
#include "MonomialBasis.hpp"
template<unsigned int dim, typename T, typename MatrixType>
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);
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)
: 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, typename T, typename MatrixType>
MatrixType &DcpseRhs<dim, T, MatrixType>::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<2, double>({0,0}));
}
return b;
}
#endif //OPENFPM_PDATA_DCPSERHS_HPP
......@@ -42,6 +42,9 @@ public:
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>
......@@ -155,5 +158,12 @@ Monomial<dim> Monomial<dim>::getDerivative(const Point<dim, unsigned int> differ
return Monomial(res);
}
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
......@@ -17,10 +17,12 @@ private:
std::vector<Monomial<dim>> basis;
public:
MonomialBasis(std::vector<unsigned int> degrees, unsigned int convergenceOrder);
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);
......@@ -29,11 +31,13 @@ public:
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);
MonomialBasis<dim> getDerivative(Point<dim, unsigned int> differentialOrder) const;
bool operator==(const MonomialBasis &other) const;
......@@ -57,7 +61,7 @@ private:
//// Definitions below
template<unsigned int dim>
MonomialBasis<dim>::MonomialBasis(std::vector<unsigned int> degrees, unsigned int convergenceOrder)
MonomialBasis<dim>::MonomialBasis(const std::vector<unsigned int> &degrees, unsigned int convergenceOrder)
{
generateBasis(degrees, convergenceOrder);
}
......@@ -85,6 +89,12 @@ 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)
{
......@@ -137,10 +147,10 @@ const std::vector<Monomial<dim>> &MonomialBasis<dim>::getElements() const
}
template<unsigned int dim>
MonomialBasis<dim> MonomialBasis<dim>::getDerivative(const Point<dim, unsigned int> differentialOrder)
MonomialBasis<dim> MonomialBasis<dim>::getDerivative(const Point<dim, unsigned int> differentialOrder) const
{
std::vector<Monomial<dim>> derivatives;
for (const auto & monomial : getElements())
for (const auto &monomial : getElements())
{
derivatives.push_back(monomial.getDerivative(differentialOrder));
}
......@@ -153,4 +163,8 @@ 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 28/03/19.
//
#define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>
#include <DCPSE/MonomialBasis.hpp>
#include <DCPSE/DcpseRhs.hpp>
#include "../../openfpm_numerics/src/DMatrix/EMatrix.hpp"
BOOST_AUTO_TEST_SUITE(DcpseRhs_tests)
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, double, EMatrix<double, Eigen::Dynamic, 1>> rhs(mb, m);
rhs.getVector(b);
// for (unsigned int i = 0; i < mb.size(); ++i)
// {
// std::cout << b(i) << std::endl;
// }
// Validation
const MonomialBasis<2> &Dmb = mb.getDerivative(m);
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, double, EMatrix<double, Eigen::Dynamic, Eigen::Dynamic>> rhs(mb, m);
rhs.getVector(b);
// for (unsigned int i = 0; i < mb.size(); ++i)
// {
// std::cout << b(i) << std::endl;
// }
// Validation
const MonomialBasis<2> &Dmb = mb.getDerivative(m);
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_SUITE_END()
\ No newline at end of file
Markdown is supported
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