FESA/src/fesa/elements/beam/beam3d2.cpp

#include <fesa/elements/beam/beam3d2.hpp>

#include <array>
#include <cmath>
#include <cstddef>
#include <optional>
#include <string>
#include <utility>

#include <fesa/fem/gauss_rule.hpp>
#include <fesa/fem/line2_shape.hpp>

namespace fesa {
namespace {

using StrainMatrix =
    std::array<std::array<double, 12>, 6>;

BeamKernelResult error_result(std::string code, std::string message) {
    BeamKernelResult result;
    result.diagnostics.push_back({
        DiagnosticStage::model,
        Severity::error,
        std::move(code),
        std::move(message),
        std::nullopt,
    });
    return result;
}

bool is_positive_finite(const double value) {
    return std::isfinite(value) && value > 0.0;
}

std::optional<BeamKernelResult> validate_properties(
    const Beam3D2Input& input) {
    if (!std::isfinite(input.material.young) ||
        !std::isfinite(input.material.poisson)) {
        return error_result(
            "model.nonfinite_value",
            "Beam kernel requires finite elastic constants.");
    }
    if (input.material.young <= 0.0 ||
        input.material.poisson <= -1.0 ||
        input.material.poisson >= 0.5) {
        return error_result(
            "model.invalid_material",
            "Beam kernel requires E > 0 and -1 < nu < 0.5.");
    }

    if (!is_positive_finite(input.section.area) ||
        !is_positive_finite(input.section.iy) ||
        !is_positive_finite(input.section.iz) ||
        !is_positive_finite(input.section.torsion_j) ||
        !is_positive_finite(input.section.shear_area_y) ||
        !is_positive_finite(input.section.shear_area_z)) {
        return error_result(
            "model.invalid_section",
            "Beam kernel requires positive finite section properties.");
    }
    return std::nullopt;
}

StrainMatrix strain_matrix(
    const double xi,
    const double jacobian) {
    StrainMatrix strain{};
    const auto shape = line2_shape(xi);
    const auto natural_derivative = line2_shape_derivative(xi);

    for (std::size_t node = 0; node < shape.size(); ++node) {
        const std::size_t offset = node * 6;
        const double derivative = natural_derivative[node] / jacobian;
        strain[0][offset] = derivative;
        strain[1][offset + 1] = derivative;
        strain[1][offset + 5] = -shape[node];
        strain[2][offset + 2] = derivative;
        strain[2][offset + 4] = shape[node];
        strain[3][offset + 3] = derivative;
        strain[4][offset + 4] = derivative;
        strain[5][offset + 5] = derivative;
    }
    return strain;
}

template <std::size_t ComponentCount>
void integrate_components(
    Matrix12& stiffness,
    const std::span<const GaussPoint1D> rule,
    const std::array<std::size_t, ComponentCount>& components,
    const std::array<double, 6>& constitutive,
    const double jacobian) {
    for (const GaussPoint1D& point : rule) {
        const StrainMatrix strain = strain_matrix(point.xi, jacobian);
        const double integration_weight = jacobian * point.weight;
        for (std::size_t row = 0; row < stiffness.size(); ++row) {
            for (std::size_t column = row;
                 column < stiffness[row].size();
                 ++column) {
                double entry = 0.0;
                for (const std::size_t component : components) {
                    entry +=
                        strain[component][row] *
                        constitutive[component] *
                        strain[component][column];
                }
                stiffness[row][column] +=
                    entry * integration_weight;
            }
        }
    }
}

void mirror_upper_triangle(Matrix12& matrix) {
    for (std::size_t row = 0; row < matrix.size(); ++row) {
        for (std::size_t column = row + 1;
             column < matrix[row].size();
             ++column) {
            matrix[column][row] = matrix[row][column];
        }
    }
}

Matrix12 transform_stiffness(
    const Matrix12& local,
    const Matrix12& transformation) {
    Matrix12 global{};
    for (std::size_t row = 0; row < global.size(); ++row) {
        for (std::size_t column = row;
             column < global[row].size();
             ++column) {
            double entry = 0.0;
            for (std::size_t local_row = 0;
                 local_row < local.size();
                 ++local_row) {
                for (std::size_t local_column = 0;
                     local_column < local[local_row].size();
                     ++local_column) {
                    entry +=
                        transformation[local_row][row] *
                        local[local_row][local_column] *
                        transformation[local_column][column];
                }
            }
            global[row][column] = entry;
        }
    }
    mirror_upper_triangle(global);
    return global;
}

bool is_finite(const Matrix12& matrix) {
    for (const auto& row : matrix) {
        for (const double value : row) {
            if (!std::isfinite(value)) {
                return false;
            }
        }
    }
    return true;
}

}  // namespace

BeamKernelResult compute_beam3d2(const Beam3D2Input& input) {
    if (const auto validation = validate_properties(input);
        validation.has_value()) {
        return *validation;
    }

    BeamFrameResult frame_result = make_beam_frame(
        input.coordinates[0],
        input.coordinates[1],
        input.section.orientation);
    if (!frame_result.frame.has_value()) {
        return {std::nullopt, std::move(frame_result.diagnostics)};
    }

    const Vec3 axis{
        input.coordinates[1].x - input.coordinates[0].x,
        input.coordinates[1].y - input.coordinates[0].y,
        input.coordinates[1].z - input.coordinates[0].z,
    };
    const double length = std::hypot(axis.x, axis.y, axis.z);
    const double jacobian = length / 2.0;
    if (!std::isfinite(jacobian) || jacobian <= 0.0) {
        return error_result(
            "model.zero_length_element",
            "Beam kernel requires a representable positive Jacobian.");
    }

    const double shear_modulus =
        input.material.young / (2.0 * (1.0 + input.material.poisson));
    const std::array<double, 6> constitutive{
        input.material.young * input.section.area,
        shear_modulus * input.section.shear_area_y,
        shear_modulus * input.section.shear_area_z,
        shear_modulus * input.section.torsion_j,
        input.material.young * input.section.iy,
        input.material.young * input.section.iz,
    };
    for (const double value : constitutive) {
        if (!is_positive_finite(value)) {
            return error_result(
                "model.nonfinite_value",
                "Beam constitutive stiffness is not finite and positive.");
        }
    }

    Matrix12 local{};
    integrate_components(
        local,
        gauss_rule_1d(2),
        std::array<std::size_t, 4>{0, 3, 4, 5},
        constitutive,
        jacobian);
    integrate_components(
        local,
        gauss_rule_1d(1),
        std::array<std::size_t, 2>{1, 2},
        constitutive,
        jacobian);
    mirror_upper_triangle(local);

    const Matrix12 transformation =
        beam_transformation(*frame_result.frame);
    Matrix12 global = transform_stiffness(local, transformation);
    if (!is_finite(local) || !is_finite(global)) {
        return error_result(
            "model.nonfinite_value",
            "Beam stiffness contains a nonfinite value.");
    }

    return {
        Beam3D2Contribution{
            std::move(local),
            std::move(global),
            *frame_result.frame,
        },
        {},
    };
}

}  // namespace fesa