Skip to content
Snippets Groups Projects

Add files for high order integration with quadratures

3 files
+ 428
0
Compare changes
  • Side-by-side
  • Inline

Files

+ 166
0
#ifndef EDGE_INTEGRATOR_HPP
#define EDGE_INTEGRATOR_HPP
#include <analysis/IQuadratureDescriptor.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/LineTransformation.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/ItemId.hpp>
#include <utils/Array.hpp>
#include <functional>
class EdgeIntegrator
{
private:
template <size_t Dimension>
class AllEdges
{
private:
const Connectivity<Dimension>& m_connectivity;
public:
using index_type = EdgeId;
PUGS_INLINE
EdgeId
edgeId(const EdgeId edge_id) const
{
return edge_id;
}
PUGS_INLINE
size_t
size() const
{
return m_connectivity.numberOfEdges();
}
PUGS_INLINE
AllEdges(const Connectivity<Dimension>& connectivity) : m_connectivity{connectivity} {}
};
template <template <typename T> typename ArrayT>
class EdgeList
{
private:
const ArrayT<EdgeId>& m_edge_list;
public:
using index_type = typename ArrayT<EdgeId>::index_type;
PUGS_INLINE
EdgeId
edgeId(const index_type index) const
{
return m_edge_list[index];
}
PUGS_INLINE
size_t
size() const
{
return m_edge_list.size();
}
PUGS_INLINE
EdgeList(const ArrayT<EdgeId>& edge_list) : m_edge_list{edge_list} {}
};
template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
static PUGS_INLINE void
_integrateTo(std::function<OutputType(InputType)> function,
const IQuadratureDescriptor& quadrature_descriptor,
const MeshType& mesh,
const ListTypeT& edge_list,
OutputArrayT& value)
{
constexpr size_t Dimension = MeshType::Dimension;
static_assert(Dimension == 3);
static_assert(std::is_same_v<TinyVector<Dimension>, std::decay_t<InputType>>, "invalid input data type");
static_assert(std::is_same_v<std::remove_const_t<typename OutputArrayT::data_type>, OutputType>,
"invalid output data type");
using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
if constexpr (std::is_arithmetic_v<OutputType>) {
value.fill(0);
} else if constexpr (is_tiny_vector_v<OutputType> or is_tiny_matrix_v<OutputType>) {
value.fill(zero);
} else {
static_assert(std::is_same_v<OutputType, double>, "unexpected output type");
}
const auto& connectivity = mesh.connectivity();
const auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
const auto xr = mesh.xr();
const auto qf = QuadratureManager::instance().getLineFormula(quadrature_descriptor);
parallel_for(edge_list.size(), [=, &qf, &edge_list](typename ListTypeT::index_type index) {
const EdgeId edge_id = edge_list.edgeId(index);
const auto edge_to_node = edge_to_node_matrix[edge_id];
Assert(edge_to_node.size() == 2);
const LineTransformation<3> T(xr[edge_to_node[0]], xr[edge_to_node[1]]);
for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) {
const auto xi = qf.point(i_point);
value[index] += qf.weight(i_point) * T.velocityNorm() * function(T(xi));
}
});
}
public:
template <typename MeshType, typename OutputArrayT, typename OutputType, typename InputType>
static PUGS_INLINE void
integrateTo(const std::function<OutputType(InputType)>& function,
const IQuadratureDescriptor& quadrature_descriptor,
const MeshType& mesh,
OutputArrayT& value)
{
constexpr size_t Dimension = MeshType::Dimension;
_integrateTo<MeshType, OutputArrayT>(function, quadrature_descriptor, mesh,
AllEdges<Dimension>{mesh.connectivity()}, value);
}
template <typename MeshType, typename OutputArrayT, typename FunctionType>
static PUGS_INLINE void
integrateTo(const FunctionType& function,
const IQuadratureDescriptor& quadrature_descriptor,
const MeshType& mesh,
OutputArrayT& value)
{
integrateTo(std::function(function), quadrature_descriptor, mesh, value);
}
template <typename MeshType, typename OutputType, typename InputType, template <typename DataType> typename ArrayT>
static PUGS_INLINE ArrayT<OutputType>
integrate(const std::function<OutputType(InputType)>& function,
const IQuadratureDescriptor& quadrature_descriptor,
const MeshType& mesh,
const ArrayT<EdgeId>& edge_list)
{
ArrayT<OutputType> value(size(edge_list));
_integrateTo<MeshType, ArrayT<OutputType>>(function, quadrature_descriptor, mesh, EdgeList{edge_list}, value);
return value;
}
template <typename MeshType, typename FunctionType, template <typename DataType> typename ArrayT>
static PUGS_INLINE auto
integrate(const FunctionType& function,
const IQuadratureDescriptor& quadrature_descriptor,
const MeshType& mesh,
const ArrayT<EdgeId>& edge_list)
{
return integrate(std::function(function), quadrature_descriptor, mesh, edge_list);
}
};
#endif // EDGE_INTEGRATOR_HPP
Loading