| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // Copyright (C) 2025 Tiago V. L. Amorim (@tiagovla) | ||
| 2 | // | ||
| 3 | // This file is part of oiseau (https://github.com/tiagovla/oiseau) | ||
| 4 | // | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | |||
| 7 | #include "oiseau/dg/dg_space.hpp" | ||
| 8 | |||
| 9 | #include <array> | ||
| 10 | #include <cstddef> | ||
| 11 | #include <span> | ||
| 12 | #include <stdexcept> | ||
| 13 | #include <vector> | ||
| 14 | #include <xtensor-blas/xlinalg.hpp> | ||
| 15 | #include <xtensor/containers/xadapt.hpp> | ||
| 16 | #include <xtensor/containers/xbuffer_adaptor.hpp> | ||
| 17 | #include <xtensor/io/xio.hpp> | ||
| 18 | #include <xtensor/views/xslice.hpp> | ||
| 19 | #include <xtensor/views/xview.hpp> | ||
| 20 | |||
| 21 | #include "oiseau/dg/nodal/element.hpp" | ||
| 22 | #include "oiseau/dg/nodal/ref_element.hpp" | ||
| 23 | #include "oiseau/mesh/cell.hpp" | ||
| 24 | #include "oiseau/mesh/mesh.hpp" | ||
| 25 | |||
| 26 | namespace oiseau::dg { | ||
| 27 | |||
| 28 | ✗ | DGSpace::DGSpace(const mesh::Mesh& mesh, const std::vector<unsigned>& orders) | |
| 29 | ✗ | : m_mesh(mesh), m_orders(orders) { | |
| 30 | ✗ | m_elements.reserve(orders.size()); | |
| 31 | |||
| 32 | ✗ | auto topology = mesh.topology(); | |
| 33 | ✗ | auto geometry = mesh.geometry(); | |
| 34 | ✗ | auto cell_types = topology.cell_types(); | |
| 35 | |||
| 36 | ✗ | std::array<std::size_t, 2> shape = {geometry.x().size() / geometry.dim(), geometry.dim()}; | |
| 37 | ✗ | auto nodes = xt::adapt(geometry.x().data(), geometry.x().size(), xt::no_ownership(), shape); | |
| 38 | |||
| 39 | ✗ | for (std::size_t i = 0; i < cell_types.size(); ++i) { | |
| 40 | ✗ | const auto& cell_type = cell_types[i]; | |
| 41 | ✗ | mesh::CellKind kind = cell_type->kind(); | |
| 42 | |||
| 43 | nodal::RefElementType ref_type; | ||
| 44 | ✗ | switch (kind) { | |
| 45 | ✗ | case mesh::CellKind::Triangle: | |
| 46 | ✗ | ref_type = nodal::RefElementType::Triangle; | |
| 47 | ✗ | break; | |
| 48 | ✗ | case mesh::CellKind::Quadrilateral: | |
| 49 | ✗ | ref_type = nodal::RefElementType::Quadrilateral; | |
| 50 | // FIX: wrong quadrilateral ordering | ||
| 51 | ✗ | break; | |
| 52 | ✗ | default: | |
| 53 | ✗ | throw std::runtime_error("Unsupported cell type"); | |
| 54 | } | ||
| 55 | |||
| 56 | ✗ | auto interp_elem = nodal::get_ref_element(ref_type, 1); | |
| 57 | ✗ | auto ref_elem = nodal::get_ref_element(ref_type, orders[i]); | |
| 58 | ✗ | auto x_view = xt::view(nodes, xt::keep(topology.conn()[i]), xt::all()); | |
| 59 | |||
| 60 | ✗ | auto inv_v = xt::linalg::inv(interp_elem->v()); | |
| 61 | ✗ | auto v = interp_elem->vandermonde(ref_elem->r()); | |
| 62 | ✗ | auto interp_x = xt::linalg::dot(xt::linalg::dot(v, inv_v), x_view); | |
| 63 | |||
| 64 | ✗ | m_elements.emplace_back(ref_elem, interp_x); | |
| 65 | ✗ | } | |
| 66 | // TODO(tiagovla): clean up this mess, introduce proper api | ||
| 67 | ✗ | } | |
| 68 | ✗ | std::span<const nodal::Element> DGSpace::elements() const { return {m_elements}; } | |
| 69 | ✗ | std::span<const unsigned> DGSpace::orders() const { return {m_orders}; } | |
| 70 | |||
| 71 | } // namespace oiseau::dg | ||
| 72 |