Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/third-party/expolit/include/expolit/integration.hh
1fc7644bad52features/netcdf-mpi

Last Change
Raw File

integration.hh
View Options

/*
* Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Expolit is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Expolit is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
* more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Expolit. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef INTEGRATION_HH
#define INTEGRATION_HH
#include "differentiation.hh"
#include "exponential.hh"
#include "litteral.hh"
#include "polynomial.hh"
#include "types.hh"
#include <thrust/pair.h>
namespace expolit {
namespace detail {
template <std::size_t i>
struct for_integ {
template <typename CT, std::size_t n>
static constexpr void integ(std::array<CT, n+1>& res,
const std::array<CT, n>& a) {
std::get<i>(res) = a[i-1] / static_cast<CT>(i);
for_integ<i-1>::integ(res, a);
}
};
template <>
struct for_integ<1> {
template <typename CT, std::size_t n>
static constexpr void integ(std::array<CT, n+1>& res,
const std::array<CT, n>& a) {
std::get<1>(res) = a[0];
}
};
} // namespace detail
/// Integrate polynomial
template <typename CT, UInt order>
constexpr auto integrate(const Polynomial<CT, order>& p) {
Polynomial<CT, order + 1> ip;
detail::for_integ<ip.terms-1>::integ(ip.coeffs, p.coeffs);
return ip;
}
/// Integrate litteral
template <typename Tag>
constexpr auto integrate(const Litteral<Tag>& q) {
constexpr Polynomial<Int, 1> id({0, 1});
return q * id;
}
/// Integrate constant * expression
template <typename CT, typename Derived>
constexpr auto integrate(const Product<Constant<CT>, Derived>& p) {
return p.operands.first * integrate(p.operands.second);
}
/// Integrate expression * constant
template <typename CT, typename Derived>
constexpr auto integrate(const Product<Derived, Constant<CT>>& p) {
return integrate(p.commute());
}
/// Integrate litteral * expression
template <typename Tag, typename Derived>
constexpr auto integrate(const Product<Litteral<Tag>, Derived>& p) {
return p.operands.first * integrate(p.operands.second);
}
/// Integrate expression * litteral
template <typename Tag, typename Derived>
constexpr auto integrate(const Product<Derived, Litteral<Tag>>& p) {
return integrate(p.commute());
}
/// Integrate litteral / expression
template <typename Tag, typename Derived>
constexpr auto integrate(const Division<Derived, Litteral<Tag>>& p) {
return integrate(p.operands.first) / p.operands.second;
}
/// Integrate exponential of order 1 polynomial
template <typename CT>
constexpr auto integrate(const Exponential<Polynomial<CT, 1>>& e) {
return Constant<CT>({1 / e.expression.coeffs.back()}) * e;
}
/// Integrate exponential of litteral * order 1 polynomial
template <typename Tag, typename CT>
constexpr auto
integrate(const Exponential<Product<Litteral<Tag>, Polynomial<CT, 1>>>& e) {
return e * Constant<CT>({1 / e.expression.operands.second.coeffs.back()}) /
e.expression.operands.first;
}
/// Integrate exponential of order 1 polynomial * litteral
template <typename Tag, typename CT>
constexpr auto
integrate(const Exponential<Product<Polynomial<CT, 1>, Litteral<Tag>>>& e) {
return integrate(e.op(e.expression.commute()));
}
/// Integration by parts of expression * polynomial
template <typename Derived, typename CT, UInt order>
constexpr auto integrate(const Product<Derived, Polynomial<CT, order>>& p) {
auto u = integrate(p.operands.first);
return u * p.operands.second -
integrate(u * differentiate(p.operands.second));
}
/// Integration by parts of polynomial * expression
template <typename Derived, typename CT, UInt order>
constexpr auto integrate(const Product<Polynomial<CT, order>, Derived>& p) {
return integrate(p.commute());
}
/// Integration of sum
template <typename Der1, typename Der2>
constexpr auto integrate(const Sum<Der1, Der2>& s) {
return integrate(s.operands.first) + integrate(s.operands.second);
}
/// Compute definite integral
template <typename T, typename Derived>
constexpr auto definite_integral(const std::pair<T, T>& bounds,
const Expression<Derived>& exp) {
auto integ = integrate(exp.downcast());
return integ(bounds.second) - integ(bounds.first);
}
/// Compute definite integral
template <typename T, typename Derived>
constexpr auto definite_integral(const thrust::pair<T, T>& bounds,
const Expression<Derived>& exp) {
auto integ = integrate(exp.downcast());
return integ(bounds.second) - integ(bounds.first);
}
} // namespace expolit
#endif
c4science · Help