math.h

// Copyright (C) 2021 Igor Baratta
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier:    LGPL-3.0-or-later
 
#pragma once
 
#include <xtensor/xarray.hpp>
#include <xtensor/xfixed.hpp>
#include <xtensor/xtensor.hpp>
 
 
namespace basix::math
{
 
template <typename U, typename V>
xt::xtensor<typename U::value_type, 2> outer(const U& u, const V& v)
{
  xt::xtensor<typename U::value_type, 2> results({u.size(), v.size()});
  for (std::size_t i = 0; i < u.size(); i++)
    for (std::size_t j = 0; j < u.size(); j++)
      results(i, j) = u(i) * v(j);
 
  return results;
}
 
template <typename U, typename V>
xt::xtensor_fixed<typename U::value_type, xt::xshape<3>> cross(const U& u,
                                                               const V& v)
{
  assert(u.size() == 3);
  assert(v.size() == 3);
  return {u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2],
          u[0] * v[1] - u[1] * v[0]};
}
 
std::pair<xt::xtensor<double, 1>,
          xt::xtensor<double, 2, xt::layout_type::column_major>>
eigh(const xt::xtensor<double, 2>& A);
 
xt::xarray<double, xt::layout_type::column_major>
solve(const xt::xtensor<double, 2>& A, const xt::xarray<double>& B);
 
bool is_singular(const xt::xtensor<double, 2>& A);
 
void dot(const xt::xtensor<double, 2>& A, const xt::xtensor<double, 2>& B,
         xt::xtensor<double, 2>& C);
 
xt::xtensor<double, 2> dot(const xt::xtensor<double, 2>& A,
                           const xt::xtensor<double, 2>& B);
 
} // namespace basix::math