BCSinusoidalTemperature.h

Go to the documentation of this file.

#ifndef BOUNDARY_CONDITION_SINUSOIDAL_TEMPERATURE_H
#define BOUNDARY_CONDITION_SINUSOIDAL_TEMPERATURE_H
 
#include <cmath>
#include <memory>
#include <numbers>
#include <string>
#include <valarray>
#include <vector>
 
#include "GenericBoundaryCondition.h"
#include "SimulationClassBase.h"
 
template <typename T>
class BCSinusoidalTemperature : public GenericBoundaryCondition<T>
{
    private:
        T m_temperature_amplitude{};
        T m_base_temperature{};
        T m_period{};
        T m_dt{};
        int m_numerical_layer{};
        std::vector<std::reference_wrapper<T>> m_capacitance_extra_val{};
        std::vector<std::reference_wrapper<T>> m_stiffness_extra_val{};
        std::vector<std::reference_wrapper<T>> m_forcing_extra_val{};
 
    public:
        void writeBoundaryCondition() override;
        void getReferenceToElements(GenericMatrix<T>& capacitance_matrix,
                                    GenericMatrix<T>& stiffness_matrix,
                                    std::valarray<T>& forcing_vector) override;
 
        BCSinusoidalTemperature(SimulationClassBase<T>* sim, std::vector<int> boundary_points,
                                std::vector<std::vector<int>> connected_points, int facet_nr,
                                T area, std::string face_prefix);
        ~BCSinusoidalTemperature() override = default;
};
 
template <typename T>
BCSinusoidalTemperature<T>::BCSinusoidalTemperature(SimulationClassBase<T>* sim,
                                                    std::vector<int> boundary_points,
                                                    std::vector<std::vector<int>> connected_points,
                                                    int facet_nr, T area, std::string face_prefix)
    : GenericBoundaryCondition<T>{sim, boundary_points, connected_points, facet_nr, face_prefix},
      m_base_temperature(
          this->m_sim->m_simulation_config.getDoubleParameters("initial_temperature")),
      m_numerical_layer(this->m_sim->m_simulation_config.getIntParameters("numerical_layers"))
{
    try
    {
        m_temperature_amplitude = this->m_sim->m_simulation_config.getDoubleParameters(
            this->m_face_indicator_prefix + "temperature_amplitude");
        m_period = this->m_sim->m_simulation_config.getDoubleParameters(
            this->m_face_indicator_prefix + "period");
    }
    catch (const BadInput&)
    {
        std::cerr << "[BCSinusoidalTemperature]: Couldn't find amplitude and/or period value in "
                     "simulation "
                     "config under: "
                  << this->m_face_indicator_prefix << "temperature_amplitude / "
                  << this->m_face_indicator_prefix << "period.";
        throw BadInput(static_cast<std::string>("[BCSinusoidalTemperature]: Couldn't find "
                                                "amplitude and/or period value in simulation "
                                                "config under: ")
                       + this->m_face_indicator_prefix
                       + static_cast<std::string>("temperature_amplitude / ")
                       + this->m_face_indicator_prefix + static_cast<std::string>("period."));
    }
    try
    {
        m_dt = this->m_sim->m_simulation_config.getDoubleParameters("time_delta");
    }
    catch (const BadInput& e)
    {
        std::cerr << "[BCSinusoidalTemperature]: " << e.what() << '\n';
        throw BadInput(static_cast<std::string>("[BCSinusoidalTemperature]: ")
                       + static_cast<std::string>(e.what()));
    }
}
 
template <typename T>
void BCSinusoidalTemperature<T>::getReferenceToElements(GenericMatrix<T>& capacitance_matrix,
                                                        GenericMatrix<T>& stiffness_matrix,
                                                        std::valarray<T>& forcing_vector)
{
    for (int i{0}; i < this->m_boundary_points.size(); i++)
    {
        std::vector<std::reference_wrapper<T>> temp_cap_mat{};
        std::vector<std::reference_wrapper<T>> temp_stiff_mat{};
        for (int j{0}; j < this->m_connected_points[i].size(); j++)
        {
            temp_cap_mat.emplace_back(
                capacitance_matrix(this->m_boundary_points[i], this->m_connected_points[i][j]));
            temp_stiff_mat.emplace_back(
                stiffness_matrix(this->m_boundary_points[i], this->m_connected_points[i][j]));
        }
        this->m_capacitance_matrix_elements.emplace_back(temp_cap_mat);
        this->m_stiffness_matrix_elements.emplace_back(temp_stiff_mat);
        // These extra values are needed for dirichlet boundary conditions in the FEM scheme
        this->m_forcing_vector_elements.emplace_back(forcing_vector[this->m_boundary_points[i]]);
        m_capacitance_extra_val.emplace_back(
            capacitance_matrix(this->m_boundary_points[i] + 1, this->m_boundary_points[i]));
        m_stiffness_extra_val.emplace_back(
            stiffness_matrix(this->m_boundary_points[i] + 1, this->m_boundary_points[i]));
        m_forcing_extra_val.emplace_back(forcing_vector[this->m_boundary_points[i] + 1]);
    }
}
 
template <typename T>
void BCSinusoidalTemperature<T>::writeBoundaryCondition()
{
    for (int i{0}; i < this->m_boundary_points.size(); i++)
    {
        for (auto& val : this->m_capacitance_matrix_elements[i])
        {
            val *= 0.0;
        }
        for (auto& val : this->m_stiffness_matrix_elements[i])
        {
            val *= 0.0;
        }
        // Set almost everything to 0, so that the top equation in the matrix looks like T_n+1 * 1 =
        // T_new.
        this->m_stiffness_matrix_elements[i][0] += 1.0;
        this->m_forcing_vector_elements[i] *= 0.0;
        T st_mat_val{m_stiffness_extra_val[i]};
        T ca_mat_val{m_capacitance_extra_val[i]};
        m_stiffness_extra_val[i] *= 0.0;
        m_capacitance_extra_val[i] *= 0.0;
        // Calculate the temperature for this and the previous time step. The first is used for the
        // top equation and the other for the second equation
        T temperature_previous_ts{
            m_base_temperature
            - m_temperature_amplitude
                  * std::sin(2 * std::numbers::pi * (this->m_sim->time_step - 1) * m_dt
                             / m_period)};
        T temperature{
            m_base_temperature
            - m_temperature_amplitude
                  * std::sin(2 * std::numbers::pi * (this->m_sim->time_step * m_dt) / m_period)};
        this->m_forcing_vector_elements[i] += temperature;
        m_forcing_extra_val[i] -= (st_mat_val * temperature + ca_mat_val * temperature / m_dt
                                   - ca_mat_val * temperature_previous_ts / m_dt);
        // this->m_sim->forcing_vector[1] -= (st_mat_val * temperature * dt + 0.0 * ca_mat_val *
        // temperature);
    }
}
 
#endif