/*
* Copyright 2009-2020 The VOTCA Development Team (http://www.votca.org)
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*/
#define BOOST_TEST_MAIN
#define BOOST_TEST_MODULE gaussian_quadratures_test
// Standard includes
#include <fstream>
// Third party includes
#include <boost/test/unit_test.hpp>
// Local VOTCA includes
#include "votca/xtp/quadrature_factory.h"
// VOTCA includes
#include <votca/tools/eigenio_matrixmarket.h>
using namespace votca::xtp;
using namespace votca;
// defines a Gaussian as integration test
// result should by ~sqrt(pi)
class FunctionEvaluation {
public:
FunctionEvaluation(){};
double operator()(Index j, double point, bool symmetry) const {
double factor = 0.0;
// this is only here to staisfy the unused variable warning
if (j < 100) {
if (symmetry) {
factor = 2.0;
} else {
factor = 1.0;
}
}
return factor * exp(-std::pow(point, 2));
}
};
BOOST_AUTO_TEST_SUITE(gaussian_quadratures_test)
BOOST_AUTO_TEST_CASE(gauss_legendre) {
QuadratureFactory::RegisterAll();
std::unique_ptr<GaussianQuadratureBase> _gq =
Quadratures().Create("legendre");
std::vector<int> orders{8, 10, 12, 14, 16, 18, 20, 40, 100};
FunctionEvaluation f = FunctionEvaluation();
Eigen::VectorXd integrals(9);
for (Index i = 0; i < 9; i++) {
_gq->configure(orders[i]);
integrals(i) = _gq->Integrate(f);
}
Eigen::VectorXd integrals_ref =
votca::tools::EigenIO_MatrixMarket::ReadVector(
std::string(XTP_TEST_DATA_FOLDER) +
"/gaussian_quadratures/gauss_legendre.mm");
bool check_integral = integrals.isApprox(integrals_ref, 1e-10);
if (!check_integral) {
std::cout << "Gauss-Legendre" << std::endl;
std::cout << integrals << std::endl;
std::cout << "Gauss-Legendre ref" << std::endl;
std::cout << integrals_ref << std::endl;
}
BOOST_CHECK_EQUAL(check_integral, true);
}
BOOST_AUTO_TEST_CASE(modified_gauss_legendre) {
QuadratureFactory::RegisterAll();
std::unique_ptr<GaussianQuadratureBase> _gq =
std::unique_ptr<GaussianQuadratureBase>(
Quadratures().Create("modified_legendre"));
std::vector<int> orders{8, 10, 12, 14, 16, 18, 20, 40, 100};
FunctionEvaluation f = FunctionEvaluation();
Eigen::VectorXd integrals(9);
for (Index i = 0; i < 9; i++) {
_gq->configure(orders[i]);
integrals(i) = _gq->Integrate(f);
}
Eigen::VectorXd integrals_ref =
votca::tools::EigenIO_MatrixMarket::ReadVector(
std::string(XTP_TEST_DATA_FOLDER) +
"/gaussian_quadratures/modified_gauss_legendre.mm");
bool check_integral = integrals.isApprox(integrals_ref, 1e-10);
if (!check_integral) {
std::cout << "modified Gauss-Legendre" << std::endl;
std::cout << integrals << std::endl;
std::cout << "modified Gauss-Legendre ref" << std::endl;
std::cout << integrals_ref << std::endl;
}
BOOST_CHECK_EQUAL(check_integral, true);
}
BOOST_AUTO_TEST_CASE(gauss_laguerre) {
QuadratureFactory::RegisterAll();
std::unique_ptr<GaussianQuadratureBase> _gq =
std::unique_ptr<GaussianQuadratureBase>(Quadratures().Create("laguerre"));
std::vector<int> orders{8, 10, 12, 14, 16, 18, 20, 40, 100};
FunctionEvaluation f = FunctionEvaluation();
Eigen::VectorXd integrals(9);
for (Index i = 0; i < 9; i++) {
_gq->configure(orders[i]);
integrals(i) = _gq->Integrate(f);
}
Eigen::VectorXd integrals_ref =
votca::tools::EigenIO_MatrixMarket::ReadVector(
std::string(XTP_TEST_DATA_FOLDER) +
"/gaussian_quadratures/gauss_laguerre.mm");
bool check_integral = integrals.isApprox(integrals_ref, 1e-10);
if (!check_integral) {
std::cout << "Gauss-Laguerre" << std::endl;
std::cout << integrals << std::endl;
std::cout << "Gauss-Laguerre ref" << std::endl;
std::cout << integrals_ref << std::endl;
}
BOOST_CHECK_EQUAL(check_integral, true);
}
BOOST_AUTO_TEST_CASE(gauss_hermite) {
QuadratureFactory::RegisterAll();
std::unique_ptr<GaussianQuadratureBase> _gq =
std::unique_ptr<GaussianQuadratureBase>(Quadratures().Create("hermite"));
std::vector<int> orders{8, 10, 12, 14, 16, 18, 20, 40, 100};
FunctionEvaluation f = FunctionEvaluation();
Eigen::VectorXd integrals(9);
for (Index i = 0; i < 9; i++) {
_gq->configure(orders[i]);
integrals(i) = _gq->Integrate(f);
}
Eigen::VectorXd integrals_ref =
votca::tools::EigenIO_MatrixMarket::ReadVector(
std::string(XTP_TEST_DATA_FOLDER) +
"/gaussian_quadratures/gauss_hermite.mm");
bool check_integral = integrals.isApprox(integrals_ref, 1e-10);
if (!check_integral) {
std::cout << "Gauss-Hermite" << std::endl;
std::cout << integrals << std::endl;
std::cout << "Gauss-Hermite ref" << std::endl;
std::cout << integrals_ref << std::endl;
}
BOOST_CHECK_EQUAL(check_integral, true);
}
BOOST_AUTO_TEST_SUITE_END()