//
// Created by cflin on 6/12/23.
//

#include "ThermalLinearFEA.h"

namespace da {
    namespace sha {
        namespace top {
            Eigen::MatrixXd ThermalLinearFEA::computeKe() {
                Eigen::Vector2d GP(-1.0 / sqrt(3.0), 1.0 / sqrt(3.0));
                Eigen::Vector2d GW(1.0, 1.0);
                Eigen::Matrix<double, 8, 8> heatKe;
                heatKe.setZero();

                Eigen::Matrix<double, 8, 3> tmp{
                        {-a, -b, -c},
                        {a,  -b, -c},
                        {a,  b,  -c},
                        {-a, b,  -c},
                        {-a, -b, c},
                        {a,  -b, c},
                        {a,  b,  c},
                        {-a, b,  c}
                };

                for (int i = 0; i < 2; ++i) {
                    for (int j = 0; j < 2; ++j) {
                        for (int k = 0; k < 2; ++k) {
                            double x = GP(i), y = GP(j), z = GP(k);
                            Eigen::RowVector<double, 8> dNx{-(1.0 - y) * (1.0 - z) / 8.0,
                                                            (1.0 - y) * (1.0 - z) / 8.0,
                                                            (1.0 + y) * (1.0 - z) / 8.0,
                                                            -(1.0 + y) * (1.0 - z) / 8.0,
                                                            -(1.0 - y) * (1.0 + z) / 8.0,
                                                            (1.0 - y) * (1.0 + z) / 8.0,
                                                            (1.0 + y) * (1.0 + z) / 8.0,
                                                            -(1.0 + y) * (1.0 + z) / 8.0};
                            Eigen::RowVector<double, 8> dNy{-(1.0 - x) * (1.0 - z) / 8.0,
                                                            -(1.0 + x) * (1.0 - z) / 8.0,
                                                            (1.0 + x) * (1.0 - z) / 8.0,
                                                            (1.0 - x) * (1.0 - z) / 8.0,
                                                            -(1.0 - x) * (1.0 + z) / 8.0,
                                                            -(1.0 + x) * (1.0 + z) / 8.0,
                                                            (1.0 + x) * (1.0 + z) / 8.0,
                                                            (1.0 - x) * (1.0 + z) / 8.0};
                            Eigen::RowVector<double, 8> dNz{-(1.0 - x) * (1.0 - y) / 8.0,
                                                            -(1.0 + x) * (1.0 - y) / 8.0,
                                                            -(1.0 + x) * (1.0 + y) / 8.0,
                                                            -(1.0 - x) * (1.0 + y) / 8.0,
                                                            (1.0 - x) * (1.0 - y) / 8.0,
                                                            (1.0 + x) * (1.0 - y) / 8.0,
                                                            (1.0 + x) * (1.0 + y) / 8.0,
                                                            (1.0 - x) * (1.0 + y) / 8.0};

                            Eigen::Matrix<double, 3, 8> tmp1;
                            tmp1(0, Eigen::all) = dNx;
                            tmp1(1, Eigen::all) = dNy;
                            tmp1(2, Eigen::all) = dNz;
                            Eigen::Matrix3d J = tmp1 * tmp;
                            Eigen::Matrix3d JInv = J.inverse();
                            double JDet = J.determinant();
                            heatKe = heatKe + GW(i) * GW(j) * GW(k) * JDet * sp_material_->thermal_conductivity *
                                              (dNx.transpose() * dNx + dNy.transpose() * dNy +
                                               dNz.transpose() * dNz);
                        }
                    }
                }
                return heatKe;
            }

            Eigen::MatrixXd ThermalLinearFEA::computeBe(const Eigen::Vector3d &P) {
                double x = P(0), y = P(1), z = P(2);
                assert(x >= -1.0 && x <= 1.0 && y >= -1.0 && y <= 1.0 && z >= -1.0 && z <= 1.0);
                Eigen::Matrix<double, 3, 8> dNi;
                dNi(0,0) = 0.125 *  - 1 * (1.0 - y) * (1.0 - z);
                dNi(0,1) = 0.125 *    1 * (1.0 - y) * (1.0 - z);
                dNi(0,2) = 0.125 *    1 * (1.0 - y) * (1.0 - z);
                dNi(0,3) = 0.125 *  - 1 * (1.0 + y) * (1.0 - z);
                dNi(0,4) = 0.125 *  - 1 * (1.0 + y) * (1.0 - z);
                dNi(0,5) = 0.125 *    1 * (1.0 - y) * (1.0 + z);
                dNi(0,6) = 0.125 *    1 * (1.0 - y) * (1.0 + z);
                dNi(0,7) = 0.125 *  - 1 * (1.0 + y) * (1.0 + z);
                dNi(1,0) = 0.125 * (1.0 - x) *  - 1 * (1.0 - z);
                dNi(1,1) = 0.125 * (1.0 + x) *  - 1 * (1.0 - z);
                dNi(1,2) = 0.125 * (1.0 + x) *  - 1 * (1.0 - z);
                dNi(1,3) = 0.125 * (1.0 - x) *  + 1 * (1.0 - z);
                dNi(1,4) = 0.125 * (1.0 - x) *  + 1 * (1.0 - z);
                dNi(1,5) = 0.125 * (1.0 + x) *  - 1 * (1.0 + z);
                dNi(1,6) = 0.125 * (1.0 + x) *  - 1 * (1.0 + z);
                dNi(1,7) = 0.125 * (1.0 - x) *  + 1 * (1.0 + z);
                dNi(2,0) = 0.125 * (1.0 - x) * (1.0 - y) *  - 1;
                dNi(2,1) = 0.125 * (1.0 + x) * (1.0 - y) *  - 1;
                dNi(2,2) = 0.125 * (1.0 + x) * (1.0 + y) *  - 1;
                dNi(2,3) = 0.125 * (1.0 - x) * (1.0 + y) *  - 1;
                dNi(2,4) = 0.125 * (1.0 - x) * (1.0 - y) *  + 1;
                dNi(2,5) = 0.125 * (1.0 + x) * (1.0 - y) *  + 1;
                dNi(2,6) = 0.125 * (1.0 + x) * (1.0 + y) *  + 1;
                dNi(2,7) = 0.125 * (1.0 - x) * (1.0 + y) *  + 1;
                return dNi;
            }

            Eigen::MatrixXd ThermalLinearFEA::computeN(const Eigen::Vector3d &P) {
                Eigen::RowVector<double, 8> Ni;
                Ni.setZero();
                double x = P(0), y = P(1), z = P(2);
                assert(x >= -1.0 && x <= 1.0 && y >= -1.0 && y <= 1.0 && z >= -1.0 && z <= 1.0);
                Eigen::RowVector<double, 8> tmp;
                Ni(0) = 0.125 * (1.0 - x) * (1.0 - y) * (1.0 - z);
                Ni(1) = 0.125 * (1.0 + x) * (1.0 - y) * (1.0 - z);
                Ni(2) = 0.125 * (1.0 + x) * (1.0 + y) * (1.0 - z);
                Ni(3) = 0.125 * (1.0 - x) * (1.0 + y) * (1.0 - z);
                Ni(4) = 0.125 * (1.0 - x) * (1.0 - y) * (1.0 + z);
                Ni(5) = 0.125 * (1.0 + x) * (1.0 - y) * (1.0 + z);
                Ni(6) = 0.125 * (1.0 + x) * (1.0 + y) * (1.0 + z);
                Ni(7) = 0.125 * (1.0 - x) * (1.0 + y) * (1.0 + z);
                return Ni;
            }
        } // da
    } // sha
} // top