/**
 * \file
 *
 * \copyright
 * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 */

#pragma once

#include <cmath>

namespace ProcessLib
{
namespace Deformation
{
/// Fills a G-matrix based on given shape function dN/dx values.
template <int DisplacementDim,
          int NPOINTS,
          typename N_Type,
          typename DNDX_Type,
          typename GMatrixType>
void computeGMatrix(DNDX_Type const& dNdx,
                    GMatrixType& g_matrix,
                    const bool is_axially_symmetric,
                    N_Type const& N,
                    const double radius)
{
    static_assert(0 < DisplacementDim && DisplacementDim <= 3,
                  "LinearGMatrix::computeGMatrix: DisplacementDim must be in "
                  "range [1,3].");

    g_matrix.setZero();

    switch (DisplacementDim)
    {
        case 3:
            // The gradient coordinates are organized in the following order:
            // (1,1), (1,2), (1,3)
            // (2,1), (2,2), (2,3)
            // (3,1), (3,2), (3,3)
            for (int d = 0; d < DisplacementDim; ++d)
            {
                for (int i = 0; i < NPOINTS; ++i)
                {
                    g_matrix(d + 0 * DisplacementDim, i + 0 * NPOINTS) =
                        dNdx(d, i);
                    g_matrix(d + 1 * DisplacementDim, i + 1 * NPOINTS) =
                        dNdx(d, i);
                    g_matrix(d + 2 * DisplacementDim, i + 2 * NPOINTS) =
                        dNdx(d, i);
                }
            }
            break;
        case 2:
            // The gradient coordinates are organized in the following order:
            // (1,1), (1,2)
            // (2,1), (2,2)
            // (3,3)
            for (int d = 0; d < DisplacementDim; ++d)
            {
                for (int i = 0; i < NPOINTS; ++i)
                {
                    g_matrix(d, i) = dNdx(d, i);
                    g_matrix(d + DisplacementDim, i + NPOINTS) = dNdx(d, i);
                }
            }
            if (is_axially_symmetric)
            {
                for (int i = 0; i < NPOINTS; ++i)
                {
                    g_matrix(4, i) = N[i] / radius;
                }
            }
            break;
        default:
            break;
    }
}

}  // namespace Deformation
}  // namespace ProcessLib