DiscreteFunctionVectorIntegrator.cpp

#include <scheme/DiscreteFunctionVectorIntegrator.hpp>

#include <language/utils/IntegrateCellArray.hpp>
#include <mesh/MeshCellZone.hpp>
#include <mesh/MeshTraits.hpp>
#include <scheme/DiscreteFunctionP0Vector.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
#include <utils/Exceptions.hpp>

template <MeshConcept MeshType, typename DataType>
DiscreteFunctionVariant
DiscreteFunctionVectorIntegrator::_integrateOnZoneList() const
{
  Assert(m_zone_list.size() > 0, "no zone list provided");
  constexpr size_t Dimension = MeshType::Dimension;

  std::shared_ptr p_mesh = m_mesh->get<MeshType>();

  CellValue<bool> is_in_zone{p_mesh->connectivity()};
  is_in_zone.fill(false);

  size_t number_of_cells = 0;
  for (const auto& zone : m_zone_list) {
    auto mesh_cell_zone   = getMeshCellZone(*p_mesh, *zone);
    const auto& cell_list = mesh_cell_zone.cellList();
    for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
      const CellId cell_id = cell_list[i_cell];
      if (is_in_zone[cell_id]) {
        std::ostringstream os;
        os << "cell " << cell_id << " (number " << p_mesh->connectivity().cellNumber()[cell_id]
           << ") is present multiple times in zone list";
        throw NormalError(os.str());
      }
      ++number_of_cells;
      is_in_zone[cell_id] = true;
    }
  }

  Array<CellId> zone_cell_list{number_of_cells};
  {
    size_t i_cell = 0;
    for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
      if (is_in_zone[cell_id]) {
        zone_cell_list[i_cell++] = cell_id;
      }
    }
  }

  Table<const DataType> table =
    IntegrateCellArray<DataType(TinyVector<Dimension>)>::template integrate(m_function_id_list,
                                                                            *m_quadrature_descriptor, *p_mesh,
                                                                            zone_cell_list);

  CellArray<DataType> cell_array{p_mesh->connectivity(), m_function_id_list.size()};
  if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
    cell_array.fill(zero);
  } else {
    cell_array.fill(0);
  }

  parallel_for(
    zone_cell_list.size(), PUGS_LAMBDA(const size_t i_cell) {
      for (size_t i = 0; i < table.numberOfColumns(); ++i) {
        cell_array[zone_cell_list[i_cell]][i] = table[i_cell][i];
      }
    });

  return DiscreteFunctionP0Vector<DataType>(m_mesh, cell_array);
}

template <MeshConcept MeshType, typename DataType>
DiscreteFunctionVariant
DiscreteFunctionVectorIntegrator::_integrateGlobally() const
{
  Assert(m_zone_list.size() == 0, "invalid call when zones are defined");

  std::shared_ptr mesh       = m_mesh->get<MeshType>();
  constexpr size_t Dimension = MeshType::Dimension;

  return DiscreteFunctionP0Vector<
    DataType>(m_mesh,
              IntegrateCellArray<DataType(TinyVector<Dimension>)>::template integrate(m_function_id_list,
                                                                                      *m_quadrature_descriptor, *mesh));
}

template <MeshConcept MeshType, typename DataType>
DiscreteFunctionVariant
DiscreteFunctionVectorIntegrator::_integrate() const
{
  if (m_zone_list.size() == 0) {
    return this->_integrateGlobally<MeshType, DataType>();
  } else {
    return this->_integrateOnZoneList<MeshType, DataType>();
  }
}

template <MeshConcept MeshType>
DiscreteFunctionVariant
DiscreteFunctionVectorIntegrator::_integrate() const
{
  for (const auto& function_id : m_function_id_list) {
    const auto& function_descriptor = function_id.descriptor();
    Assert(function_descriptor.domainMappingNode().children[1]->m_data_type == ASTNodeDataType::typename_t);
    const ASTNodeDataType& data_type = function_descriptor.domainMappingNode().children[1]->m_data_type.contentType();

    switch (data_type) {
    case ASTNodeDataType::bool_t:
    case ASTNodeDataType::unsigned_int_t:
    case ASTNodeDataType::int_t:
    case ASTNodeDataType::double_t: {
      break;
    }
    default: {
      std::ostringstream os;
      os << "vector functions require scalar value type.\n"
         << "Invalid value type: " << rang::fgB::red << dataTypeName(data_type) << rang::style::reset;
      throw NormalError(os.str());
    }
    }
  }
  return this->_integrate<MeshType, double>();
}

DiscreteFunctionVariant
DiscreteFunctionVectorIntegrator::integrate() const
{
  if (m_discrete_function_descriptor->type() != DiscreteFunctionType::P0Vector) {
    throw NormalError("invalid discrete function type for vector integration");
  }

  return std::visit(
    [&](auto&& mesh) {
      using MeshType = mesh_type_t<decltype(mesh)>;

      if constexpr (is_polygonal_mesh_v<MeshType>) {
        return this->_integrate<MeshType>();
      } else {
        throw NormalError("unexpected mesh type");
      }
    },
    m_mesh->variant());
}