GridOneDim.h

Go to the documentation of this file.

#ifndef GRID_ONE_DIM_H
#define GRID_ONE_DIM_H
 
#include <algorithm>
#include <memory>
#include <string>
#include <vector>
 
#include "BoundaryConditionRegistry.h"
#include "CellRegistry.h"
#include "CsvParser.h"
#include "GenericBoundaryCondition.h"
#include "GenericElement.h"
#include "GridBase.h"
#include "IniParser.h"
#include "ShapeBase.h"
#include "SimulationClassBase.h"
 
template <typename T>
class GridOneDim : public GridBase<T>
{
    public:
        void createGrid() override;
 
        T calculateCellLength();
 
        GridOneDim(SimulationClassBase<T>* sim_ptr, std::shared_ptr<ShapeBase<T>> mesh_ptr);
        ~GridOneDim() = default;
 
    private:
        void setupBoundaryCell(const std::vector<int>& node_numbers,
                               std::string_view boundary_condition_type,
                               const std::map<std::string_view, std::vector<int>>& boundary_points,
                               int current_facet, T areas, const std::string& face_prefix);
        void incrementVector(std::vector<int>& vec);
        int getNumberOfFacets() const;
        std::vector<T> getGridArray(int number_of_elements);
};
 
template <typename T>
GridOneDim<T>::GridOneDim(SimulationClassBase<T>* sim_ptr, std::shared_ptr<ShapeBase<T>> mesh_ptr)
    : GridBase<T>{sim_ptr, mesh_ptr}
{
    createGrid();
}
 
template <typename T>
void GridOneDim<T>::incrementVector(std::vector<int>& vec)
{
    for (auto& val : vec)
    {
        val++;
    }
}
 
template <typename T>
T GridOneDim<T>::calculateCellLength()
{
    bool depth_present{false};
    bool length_present{false};
    double cl_depth{};
    double cl_length{};
    try
    {
        int numerical_layers{this->sim->m_simulation_config.getIntParameters("numerical_layers")};
        double max_depth{this->sim->m_simulation_config.getDoubleParameters("max_depth")};
        cl_depth = max_depth / static_cast<double>(numerical_layers);
        depth_present = true;
    }
    catch (const BadInput& e)
    {
    }
    try
    {
        cl_length = this->sim->m_simulation_config.getDoubleParameters("cell_length");
        length_present = true;
    }
    catch (const BadInput& e)
    {
    }
 
    if (depth_present && !length_present && cl_depth >= 0)
    {
        return cl_depth;
    }
    if (length_present && !depth_present && cl_length >= 0)
    {
        return cl_length;
    }
    else if (!length_present && !depth_present)
    {
        throw std::runtime_error(
            "Could not find the entries 'cell_length' or 'max_depth' within the ini file. Please "
            "supply one of these entries!");
    }
    else
    {
        const double epsilon{1e-8};
        if (std::abs(cl_length - cl_depth) <= epsilon
            || std::abs(cl_length - cl_depth)
                   <= epsilon * std::max(std::abs(cl_depth), std::abs(cl_length)))
        {
            return cl_length;
        }
        throw std::runtime_error(
            "Both entries 'cell_length' and 'max_depth' are present within the ini file and result "
            "in different lengths. Please only supply one of these entries or have them produce "
            "the same result!");
    }
}
 
template <typename T>
std::vector<T> GridOneDim<T>::getGridArray(int number_of_elements)
{
    try
    {
        CsvParser grid_data{this->sim->m_simulation_config.getStringParameters("grid")};
        std::vector<T> depth{grid_data.getCsvContents()[0]};
        std::vector<T> length_arr;
        for (int i{1}; i < depth.size(); i++)
        {
            length_arr.push_back(depth[i] - depth[i - 1]);
        }
        // Overwrite the numerical_layers parameter with the one given by the grid file
        if (this->sim->m_simulation_config.getIntParameters("numerical_layers") != depth.size())
        {
            std::cerr
                << "[Warning] Missmatch between number of numerical layers within provided file "
                   "and given value in the ini! The value in the provided file will be used!\n";
            this->sim->m_simulation_config.forceValueOverwrite("numerical_layers",
                                                               std::to_string(depth.size()));
        }
        return std::move(length_arr);
    }
    catch (const BadInput&)
    {
        std::vector<T> length_arr(number_of_elements, calculateCellLength());
        length_arr[0] = length_arr[0] / 2;
        return length_arr;
    }
}
 
template <typename T>
void GridOneDim<T>::setupBoundaryCell(
    const std::vector<int>& node_numbers, std::string_view boundary_condition_type,
    const std::map<std::string_view, std::vector<int>>& boundary_points, int current_facet, T areas,
    const std::string& face_prefix)
{
    std::vector<std::vector<int>> connected_points({{node_numbers[0], node_numbers[0] + 1}});
    std::map<std::string_view, std::unique_ptr<GenericBoundaryCondition<T>>> boundary_condition;
    boundary_condition[boundary_condition_type] = createBoundaryCondition<T>(
        boundary_condition_type, this->sim, boundary_points.at(boundary_condition_type),
        connected_points, current_facet, areas, face_prefix);
    this->all_cells[this->all_cells.size() - 1]->setBoundaryConditions(boundary_condition);
}
 
template <typename T>
int GridOneDim<T>::getNumberOfFacets() const
{
#ifdef MPI_PRESENT
    return this->sim->m_simulation_config.getIntParameters("number_of_facets");
#else
    return static_cast<int>(this->m_poly_mesh->getNumberOfFacets());
#endif
}
 
template <typename T>
void GridOneDim<T>::createGrid()
{
    ElementProperties<T> props{};
    const std::string cell_type{this->sim->m_simulation_config.getStringParameters("cell_type")};
    const std::string top_boundary_condition_type{
        this->sim->m_simulation_config.getStringParameters("top_boundary_condition_type")};
    const std::string bottom_boundary_condition_type{
        this->sim->m_simulation_config.getStringParameters("bottom_boundary_condition_type")};
    int number_of_elements{this->sim->m_simulation_config.getIntParameters("numerical_layers") - 1};
    int nr_of_facets{getNumberOfFacets()};
    std::vector length_arr{getGridArray(number_of_elements)};
 
    // This is ugly, please revise
    std::valarray<double> cell_length(number_of_elements);
    for (int i{0}; i < number_of_elements; i++)
    {
        cell_length[i] = length_arr[i];
    }
    this->sim->m_field_map.insert({"CellLength", cell_length});
 
    this->grid_size =
        nr_of_facets * this->sim->m_simulation_config.getIntParameters("numerical_layers");
 
    // this->sim->set_matrix_size(nr_of_facets *
    // this->sim->m_simulation_config.getIntParameters("numerical_layers"));
    std::vector<T> areas{this->m_poly_mesh->getAreas()};
    std::vector<int> node_numbers = {0, 1};
    int num_cells{0};
    std::vector<int> boundary_vector = {-1};
    std::map<std::string_view, std::vector<int>> boundary_points_top{
        {top_boundary_condition_type, boundary_vector}};
    std::map<std::string_view, std::vector<int>> boundary_points_bottom{
        {bottom_boundary_condition_type, boundary_vector}};
    int current_facet = 0;
 
    for (int i{0}; i < nr_of_facets * (number_of_elements); i++)
    {
        if (node_numbers[1] > nr_of_facets * (number_of_elements + 1))
        {
            throw std::runtime_error(
                static_cast<std::string>("Node number out of bounds in OneDimGrid. The node number "
                                         "exceeds the requested number"));
        }
        props.element_length = length_arr[i % number_of_elements];
        if (i % number_of_elements == 0)
        {
            current_facet = i / number_of_elements;
            props.element_area = areas[current_facet];
            if (i > 0)
            {
                incrementVector(node_numbers);  // The extra increase is needed to separate the
                                                // facets in the matrix
            }
            boundary_points_top.at(top_boundary_condition_type)[0] = node_numbers[0];
            this->all_cells.push_back({createCell<T>(cell_type, true, boundary_points_top,
                                                     node_numbers, this->sim, num_cells, props)});
            setupBoundaryCell(node_numbers, top_boundary_condition_type, boundary_points_top,
                              current_facet, areas[current_facet], "top_");
        }
        else if (i % number_of_elements == number_of_elements - 1
                 && bottom_boundary_condition_type != "None")
        {
            boundary_points_bottom.at(bottom_boundary_condition_type)[0] = node_numbers[1];
            this->all_cells.push_back({createCell<T>(cell_type, true, boundary_points_bottom,
                                                     node_numbers, this->sim, num_cells, props)});
            setupBoundaryCell(node_numbers, bottom_boundary_condition_type, boundary_points_bottom,
                              current_facet, areas[current_facet], "bottom_");
        }
        else
        {
            // boundary_points.at(boundary_condition_type)[0] = node_numbers[0]; // This value is
            // never used for internal elements
            this->all_cells.push_back({createCell<T>(cell_type, false, boundary_points_top,
                                                     node_numbers, this->sim, num_cells, props)});
        }
 
        num_cells++;
        incrementVector(node_numbers);
    }
    // std::cout << "Number interior cells: " << this->compute_cells.size() << ". Number boundary
    // cells: " << this->boundary_cells.size() << "\n";
}
 
#endif