/*=========================================================================

  Module:    V_GaussIntegration.cpp

  Copyright 2003,2006,2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
  Under the terms of Contract DE-NA0003525 with NTESS,
  the U.S. Government retains certain rights in this software.

  See LICENSE for details.

=========================================================================*/

/*
 *
 * GaussIntegration.cpp performs Gauss Integrations
 *
 * This file is part of VERDICT
 *
 */

#include "V_GaussIntegration.hpp"
#include "verdict.h"

namespace VERDICT_NAMESPACE
{
void GaussIntegration::initialize(int n, int m, int dim, int tri)
{
  numberGaussPoints = n;
  numberNodes = m;
  numberDims = dim;

  if (tri == 1)
  // triangular element
  {
    if (numberDims == 2)
    {
      totalNumberGaussPts = numberGaussPoints;
    }
    else if (numberDims == 3)
    {
      totalNumberGaussPts = numberGaussPoints;
    }
  }
  else if (tri == 0)
  {
    if (numberDims == 2)
    {
      totalNumberGaussPts = numberGaussPoints * numberGaussPoints;
    }
    else if (numberDims == 3)
    {
      totalNumberGaussPts = numberGaussPoints * numberGaussPoints * numberGaussPoints;
    }
  }
}

void GaussIntegration::get_shape_func(double shape_function[], double dndy1_at_gauss_pts[],
  double dndy2_at_gauss_pts[], double gauss_weight[])
{
  int i, j;
  for (i = 0; i < totalNumberGaussPts; i++)
  {
    for (j = 0; j < numberNodes; j++)
    {
      shape_function[i * maxNumberNodes + j] = shapeFunction[i][j];
      dndy1_at_gauss_pts[i * maxNumberNodes + j] = dndy1GaussPts[i][j];
      dndy2_at_gauss_pts[i * maxNumberNodes + j] = dndy2GaussPts[i][j];
    }
  }

  for (i = 0; i < totalNumberGaussPts; i++)
  {
    gauss_weight[i] = totalGaussWeight[i];
  }
}

void GaussIntegration::get_shape_func(double shape_function[], double dndy1_at_gauss_pts[],
  double dndy2_at_gauss_pts[], double dndy3_at_gauss_pts[], double gauss_weight[])
{
  int i, j;
  for (i = 0; i < totalNumberGaussPts; i++)
  {
    for (j = 0; j < numberNodes; j++)
    {
      shape_function[i * maxNumberNodes + j] = shapeFunction[i][j];
      dndy1_at_gauss_pts[i * maxNumberNodes + j] = dndy1GaussPts[i][j];
      dndy2_at_gauss_pts[i * maxNumberNodes + j] = dndy2GaussPts[i][j];
      dndy3_at_gauss_pts[i * maxNumberNodes + j] = dndy3GaussPts[i][j];
    }
  }

  for (i = 0; i < totalNumberGaussPts; i++)
  {
    gauss_weight[i] = totalGaussWeight[i];
  }
}

void GaussIntegration::get_gauss_pts_and_weight()
{

  switch (numberGaussPoints)
  {
    case 1:
      gaussPointY[0] = 0.0;
      gaussWeight[0] = 2.0;
      break;
    case 2:
      gaussPointY[0] = -0.577350269189626;
      gaussPointY[1] = 0.577350269189626;
      gaussWeight[0] = 1.0;
      gaussWeight[1] = 1.0;
      break;
    case 3:
      gaussPointY[0] = -0.774596669241483;
      gaussPointY[1] = 0.0;
      gaussPointY[2] = 0.774596669241483;
      gaussWeight[0] = 0.555555555555555;
      gaussWeight[1] = 0.888888888888889;
      gaussWeight[2] = 0.555555555555555;
      break;
  }
}

void GaussIntegration::calculate_shape_function_2d_quad()
{
  int ife = 0, i, j;
  double y1, y2;
  get_gauss_pts_and_weight();

  switch (numberNodes)
  {
    case 4:
      for (i = 0; i < numberGaussPoints; i++)
      {
        for (j = 0; j < numberGaussPoints; j++)
        {
          y1 = gaussPointY[i];
          y2 = gaussPointY[j];
          shapeFunction[ife][0] = 0.25 * (1 - y1) * (1 - y2);
          shapeFunction[ife][1] = 0.25 * (1 + y1) * (1 - y2);
          shapeFunction[ife][2] = 0.25 * (1 + y1) * (1 + y2);
          shapeFunction[ife][3] = 0.25 * (1 - y1) * (1 + y2);

          dndy1GaussPts[ife][0] = -0.25 * (1 - y2);
          dndy1GaussPts[ife][1] = 0.25 * (1 - y2);
          dndy1GaussPts[ife][2] = 0.25 * (1 + y2);
          dndy1GaussPts[ife][3] = -0.25 * (1 + y2);

          dndy2GaussPts[ife][0] = -0.25 * (1 - y1);
          dndy2GaussPts[ife][1] = -0.25 * (1 + y1);
          dndy2GaussPts[ife][2] = 0.25 * (1 + y1);
          dndy2GaussPts[ife][3] = 0.25 * (1 - y1);

          totalGaussWeight[ife] = gaussWeight[i] * gaussWeight[j];
          ife++;
        }
      }
      break;
    case 8:
      for (i = 0; i < numberGaussPoints; i++)
      {
        for (j = 0; j < numberGaussPoints; j++)
        {
          y1 = gaussPointY[i];
          y2 = gaussPointY[j];
          shapeFunction[ife][0] = 0.25 * (1 - y1) * (1 - y2) * (-y1 - y2 - 1);
          shapeFunction[ife][1] = 0.25 * (1 + y1) * (1 - y2) * (y1 - y2 - 1);
          shapeFunction[ife][2] = 0.25 * (1 + y1) * (1 + y2) * (y1 + y2 - 1);
          shapeFunction[ife][3] = 0.25 * (1 - y1) * (1 + y2) * (-y1 + y2 - 1);
          shapeFunction[ife][4] = 0.5 * (1 - y1 * y1) * (1 - y2);
          shapeFunction[ife][5] = 0.5 * (1 - y2 * y2) * (1 + y1);
          shapeFunction[ife][6] = 0.5 * (1 - y1 * y1) * (1 + y2);
          shapeFunction[ife][7] = 0.5 * (1 - y2 * y2) * (1 - y1);

          dndy1GaussPts[ife][0] = 0.25 * (1 - y2) * (2.0 * y1 + y2);
          dndy1GaussPts[ife][1] = 0.25 * (1 - y2) * (2.0 * y1 - y2);
          dndy1GaussPts[ife][2] = 0.25 * (1 + y2) * (2.0 * y1 + y2);
          dndy1GaussPts[ife][3] = 0.25 * (1 + y2) * (2.0 * y1 - y2);

          dndy1GaussPts[ife][4] = -y1 * (1 - y2);
          dndy1GaussPts[ife][5] = 0.5 * (1 - y2 * y2);
          dndy1GaussPts[ife][6] = -y1 * (1 + y2);
          dndy1GaussPts[ife][7] = -0.5 * (1 - y2 * y2);

          dndy2GaussPts[ife][0] = 0.25 * (1 - y1) * (2.0 * y2 + y1);
          dndy2GaussPts[ife][1] = 0.25 * (1 + y1) * (2.0 * y2 - y1);
          dndy2GaussPts[ife][2] = 0.25 * (1 + y1) * (2.0 * y2 + y1);
          dndy2GaussPts[ife][3] = 0.25 * (1 - y1) * (2.0 * y2 - y1);

          dndy2GaussPts[ife][4] = -0.5 * (1 - y1 * y1);
          dndy2GaussPts[ife][5] = -y2 * (1 + y1);
          dndy2GaussPts[ife][6] = 0.5 * (1 - y1 * y1);
          dndy2GaussPts[ife][7] = -y2 * (1 - y1);

          totalGaussWeight[ife] = gaussWeight[i] * gaussWeight[j];
          ife++;
        }
      }
      break;
  }
}

void GaussIntegration::calculate_shape_function_3d_hex()
{
  int ife = 0, i, j, k, node_id;
  double y1, y2, y3, sign_node_y1, sign_node_y2, sign_node_y3;
  double y1_term, y2_term, y3_term, y123_temp;

  get_gauss_pts_and_weight();

  switch (numberNodes)
  {
    case 8:
      for (i = 0; i < numberGaussPoints; i++)
      {
        for (j = 0; j < numberGaussPoints; j++)
        {
          for (k = 0; k < numberGaussPoints; k++)
          {
            y1 = gaussPointY[i];
            y2 = gaussPointY[j];
            y3 = gaussPointY[k];

            for (node_id = 0; node_id < numberNodes; node_id++)
            {
              get_signs_for_node_local_coord_hex(node_id, sign_node_y1, sign_node_y2, sign_node_y3);

              y1_term = 1 + sign_node_y1 * y1;
              y2_term = 1 + sign_node_y2 * y2;
              y3_term = 1 + sign_node_y3 * y3;

              shapeFunction[ife][node_id] = 0.125 * y1_term * y2_term * y3_term;
              dndy1GaussPts[ife][node_id] = 0.125 * sign_node_y1 * y2_term * y3_term;
              dndy2GaussPts[ife][node_id] = 0.125 * sign_node_y2 * y1_term * y3_term;
              dndy3GaussPts[ife][node_id] = 0.125 * sign_node_y3 * y1_term * y2_term;
            }
            totalGaussWeight[ife] = gaussWeight[i] * gaussWeight[j] * gaussWeight[k];
            ife++;
          }
        }
      }
      break;
    case 20:
      for (i = 0; i < numberGaussPoints; i++)
      {
        for (j = 0; j < numberGaussPoints; j++)
        {
          for (k = 0; k < numberGaussPoints; k++)
          {
            y1 = gaussPointY[i];
            y2 = gaussPointY[j];
            y3 = gaussPointY[k];

            for (node_id = 0; node_id < numberNodes; node_id++)
            {
              get_signs_for_node_local_coord_hex(node_id, sign_node_y1, sign_node_y2, sign_node_y3);

              y1_term = 1 + sign_node_y1 * y1;
              y2_term = 1 + sign_node_y2 * y2;
              y3_term = 1 + sign_node_y3 * y3;
              y123_temp = sign_node_y1 * y1 + sign_node_y2 * y2 + sign_node_y3 * y3 - 2.;

              switch (node_id)
              {
                case 0:
                case 1:
                case 2:
                case 3:
                case 4:
                case 5:
                case 6:
                case 7:
                {
                  shapeFunction[ife][node_id] = 0.125 * y1_term * y2_term * y3_term * y123_temp;
                  dndy1GaussPts[ife][node_id] =
                    0.125 * sign_node_y1 * y123_temp * y2_term * y3_term +
                    0.125 * y1_term * y2_term * y3_term * sign_node_y1;
                  dndy2GaussPts[ife][node_id] =
                    0.125 * sign_node_y2 * y1_term * y3_term * y123_temp +
                    0.125 * y1_term * y2_term * y3_term * sign_node_y2;
                  dndy3GaussPts[ife][node_id] =
                    0.125 * sign_node_y3 * y1_term * y2_term * y123_temp +
                    0.125 * y1_term * y2_term * y3_term * sign_node_y3;
                  break;
                }
                case 8:
                case 10:
                case 16:
                case 18:
                {
                  shapeFunction[ife][node_id] = 0.25 * (1 - y1 * y1) * y2_term * y3_term;
                  dndy1GaussPts[ife][node_id] = -0.5 * y1 * y2_term * y3_term;
                  dndy2GaussPts[ife][node_id] = 0.25 * (1 - y1 * y1) * sign_node_y2 * y3_term;
                  dndy3GaussPts[ife][node_id] = 0.25 * (1 - y1 * y1) * y2_term * sign_node_y3;
                  break;
                }
                case 9:
                case 11:
                case 17:
                case 19:
                {
                  shapeFunction[ife][node_id] = 0.25 * (1 - y2 * y2) * y1_term * y3_term;
                  dndy1GaussPts[ife][node_id] = 0.25 * (1 - y2 * y2) * sign_node_y1 * y3_term;
                  dndy2GaussPts[ife][node_id] = -0.5 * y2 * y1_term * y3_term;
                  dndy3GaussPts[ife][node_id] = 0.25 * (1 - y2 * y2) * y1_term * sign_node_y3;
                  break;
                }
                case 12:
                case 13:
                case 14:
                case 15:
                {
                  shapeFunction[ife][node_id] = 0.25 * (1 - y3 * y3) * y1_term * y2_term;
                  dndy1GaussPts[ife][node_id] = 0.25 * (1 - y3 * y3) * sign_node_y1 * y2_term;
                  dndy2GaussPts[ife][node_id] = 0.25 * (1 - y3 * y3) * y1_term * sign_node_y2;
                  dndy3GaussPts[ife][node_id] = -0.5 * y3 * y1_term * y2_term;
                  break;
                }
              }
            }
            totalGaussWeight[ife] = gaussWeight[i] * gaussWeight[j] * gaussWeight[k];
            ife++;
          }
        }
      }
      break;
  }
}

void GaussIntegration::calculate_derivative_at_nodes(
  double dndy1_at_nodes[][maxNumberNodes], double dndy2_at_nodes[][maxNumberNodes])
{
  double y1 = 0., y2 = 0.;
  int i;
  for (i = 0; i < numberNodes; i++)
  {
    switch (i)
    {
      case 0:
        y1 = -1.;
        y2 = -1.;
        break;

      case 1:
        y1 = 1.;
        y2 = -1.;
        break;

      case 2:
        y1 = 1.;
        y2 = 1.;
        break;

      case 3:
        y1 = -1.;
        y2 = 1.;
        break;

      // mid-side nodes if there is any
      case 4:
        y1 = 0.;
        y2 = -1.;
        break;

      case 5:
        y1 = 1.;
        y2 = 0.;
        break;

      case 6:
        y1 = 0.;
        y2 = 1.;
        break;

      case 7:
        y1 = -1.;
        y2 = 0.;
        break;
    }

    switch (numberNodes)
    {
      case 4:
        // dn_i/dy1 evaluated at node i
        dndy1_at_nodes[i][0] = -0.25 * (1 - y2);
        dndy1_at_nodes[i][1] = 0.25 * (1 - y2);
        dndy1_at_nodes[i][2] = 0.25 * (1 + y2);
        dndy1_at_nodes[i][3] = -0.25 * (1 + y2);

        // dn_i/dy2 evaluated at node i
        dndy2_at_nodes[i][0] = -0.25 * (1 - y1);
        dndy2_at_nodes[i][1] = -0.25 * (1 + y1);
        dndy2_at_nodes[i][2] = 0.25 * (1 + y1);
        dndy2_at_nodes[i][3] = 0.25 * (1 - y1);
        break;

      case 8:
        dndy1_at_nodes[i][0] = 0.25 * (1 - y2) * (2.0 * y1 + y2);
        dndy1_at_nodes[i][1] = 0.25 * (1 - y2) * (2.0 * y1 - y2);
        dndy1_at_nodes[i][2] = 0.25 * (1 + y2) * (2.0 * y1 + y2);
        dndy1_at_nodes[i][3] = 0.25 * (1 + y2) * (2.0 * y1 - y2);

        dndy1_at_nodes[i][4] = -y1 * (1 - y2);
        dndy1_at_nodes[i][5] = 0.5 * (1 - y2 * y2);
        dndy1_at_nodes[i][6] = -y1 * (1 + y2);
        dndy1_at_nodes[i][7] = -0.5 * (1 - y2 * y2);

        dndy2_at_nodes[i][0] = 0.25 * (1 - y1) * (2.0 * y2 + y1);
        dndy2_at_nodes[i][1] = 0.25 * (1 + y1) * (2.0 * y2 - y1);
        dndy2_at_nodes[i][2] = 0.25 * (1 + y1) * (2.0 * y2 + y1);
        dndy2_at_nodes[i][3] = 0.25 * (1 - y1) * (2.0 * y2 - y1);

        dndy2_at_nodes[i][4] = -0.5 * (1 - y1 * y1);
        dndy2_at_nodes[i][5] = -y2 * (1 + y1);
        dndy2_at_nodes[i][6] = 0.5 * (1 - y1 * y1);
        dndy2_at_nodes[i][7] = -y2 * (1 - y1);
        break;
    }
  }
}

void GaussIntegration::calculate_derivative_at_nodes_3d(double dndy1_at_nodes[][maxNumberNodes],
  double dndy2_at_nodes[][maxNumberNodes], double dndy3_at_nodes[][maxNumberNodes])
{
  double y1, y2, y3, sign_node_y1, sign_node_y2, sign_node_y3;
  double y1_term, y2_term, y3_term, y123_temp;
  int node_id, node_id_2;
  for (node_id = 0; node_id < numberNodes; node_id++)
  {
    get_signs_for_node_local_coord_hex(node_id, y1, y2, y3);

    switch (numberNodes)
    {
      case 8:
        for (node_id_2 = 0; node_id_2 < numberNodes; node_id_2++)
        {
          get_signs_for_node_local_coord_hex(node_id_2, sign_node_y1, sign_node_y2, sign_node_y3);
          y1_term = 1 + sign_node_y1 * y1;
          y2_term = 1 + sign_node_y2 * y2;
          y3_term = 1 + sign_node_y3 * y3;

          dndy1_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y1 * y2_term * y3_term;

          dndy2_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y2 * y1_term * y3_term;

          dndy3_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y3 * y1_term * y2_term;
        }
        break;
      case 20:
        for (node_id_2 = 0; node_id_2 < numberNodes; node_id_2++)
        {
          get_signs_for_node_local_coord_hex(node_id_2, sign_node_y1, sign_node_y2, sign_node_y3);

          y1_term = 1 + sign_node_y1 * y1;
          y2_term = 1 + sign_node_y2 * y2;
          y3_term = 1 + sign_node_y3 * y3;
          y123_temp = sign_node_y1 * y1 + sign_node_y2 * y2 + sign_node_y3 * y3 - 2.;
          switch (node_id_2)
          {
            case 0:
            case 1:
            case 2:
            case 3:
            case 4:
            case 5:
            case 6:
            case 7:
            {
              dndy1_at_nodes[node_id][node_id_2] =
                0.125 * sign_node_y1 * y2_term * y3_term * y123_temp +
                0.125 * y1_term * y2_term * y3_term * sign_node_y1;
              dndy2_at_nodes[node_id][node_id_2] =
                0.125 * sign_node_y2 * y1_term * y3_term * y123_temp +
                0.125 * y1_term * y2_term * y3_term * sign_node_y2;
              dndy3_at_nodes[node_id][node_id_2] =
                0.125 * sign_node_y3 * y1_term * y2_term * y123_temp +
                0.125 * y1_term * y2_term * y3_term * sign_node_y3;
              break;
            }
            case 8:
            case 10:
            case 16:
            case 18:
            {
              dndy1_at_nodes[node_id][node_id_2] = -0.5 * y1 * y2_term * y3_term;
              dndy2_at_nodes[node_id][node_id_2] = 0.25 * (1 - y1 * y1) * sign_node_y2 * y3_term;
              dndy3_at_nodes[node_id][node_id_2] = 0.25 * (1 - y1 * y1) * y2_term * sign_node_y3;
              break;
            }
            case 9:
            case 11:
            case 17:
            case 19:
            {
              dndy1_at_nodes[node_id][node_id_2] = 0.25 * (1 - y2 * y2) * sign_node_y1 * y3_term;
              dndy2_at_nodes[node_id][node_id_2] = -0.5 * y2 * y1_term * y3_term;
              dndy3_at_nodes[node_id][node_id_2] = 0.25 * (1 - y2 * y2) * y1_term * sign_node_y3;
              break;
            }
            case 12:
            case 13:
            case 14:
            case 15:
            {
              dndy1_at_nodes[node_id][node_id_2] = 0.25 * (1 - y3 * y3) * sign_node_y1 * y2_term;
              dndy2_at_nodes[node_id][node_id_2] = 0.25 * (1 - y3 * y3) * y1_term * sign_node_y2;
              dndy3_at_nodes[node_id][node_id_2] = -0.5 * y3 * y1_term * y2_term;
              break;
            }
          }
        }
        break;
    }
  }
}

void GaussIntegration::get_signs_for_node_local_coord_hex(
  int node_id, double& sign_node_y1, double& sign_node_y2, double& sign_node_y3)
{
  switch (node_id)
  {
    case 0:
      sign_node_y1 = -1.;
      sign_node_y2 = -1.;
      sign_node_y3 = -1.;
      break;
    case 1:
      sign_node_y1 = 1.;
      sign_node_y2 = -1.;
      sign_node_y3 = -1.;
      break;
    case 2:
      sign_node_y1 = 1.;
      sign_node_y2 = 1.;
      sign_node_y3 = -1.;
      break;
    case 3:
      sign_node_y1 = -1.;
      sign_node_y2 = 1.;
      sign_node_y3 = -1.;
      break;
    case 4:
      sign_node_y1 = -1.;
      sign_node_y2 = -1.;
      sign_node_y3 = 1.;
      break;
    case 5:
      sign_node_y1 = 1.;
      sign_node_y2 = -1.;
      sign_node_y3 = 1.;
      break;
    case 6:
      sign_node_y1 = 1.;
      sign_node_y2 = 1.;
      sign_node_y3 = 1.;
      break;
    case 7:
      sign_node_y1 = -1.;
      sign_node_y2 = 1.;
      sign_node_y3 = 1.;
      break;
    case 8:
      sign_node_y1 = 0;
      sign_node_y2 = -1.;
      sign_node_y3 = -1.;
      break;
    case 9:
      sign_node_y1 = 1.;
      sign_node_y2 = 0;
      sign_node_y3 = -1.;
      break;
    case 10:
      sign_node_y1 = 0;
      sign_node_y2 = 1.;
      sign_node_y3 = -1.;
      break;
    case 11:
      sign_node_y1 = -1.;
      sign_node_y2 = 0.;
      sign_node_y3 = -1.;
      break;
    case 12:
      sign_node_y1 = -1.;
      sign_node_y2 = -1.;
      sign_node_y3 = 0.;
      break;
    case 13:
      sign_node_y1 = 1.;
      sign_node_y2 = -1.;
      sign_node_y3 = 0.;
      break;
    case 14:
      sign_node_y1 = 1.;
      sign_node_y2 = 1.;
      sign_node_y3 = 0.;
      break;
    case 15:
      sign_node_y1 = -1.;
      sign_node_y2 = 1.;
      sign_node_y3 = 0.;
      break;
    case 16:
      sign_node_y1 = 0;
      sign_node_y2 = -1.;
      sign_node_y3 = 1.;
      break;
    case 17:
      sign_node_y1 = 1.;
      sign_node_y2 = 0;
      sign_node_y3 = 1.;
      break;
    case 18:
      sign_node_y1 = 0;
      sign_node_y2 = 1.;
      sign_node_y3 = 1.;
      break;
    case 19:
      sign_node_y1 = -1.;
      sign_node_y2 = 0.;
      sign_node_y3 = 1.;
      break;
    default:
      // Should not be possible to get here, but if we do, at least results will be consistent, not
      // random
      sign_node_y1 = 0.;
      sign_node_y2 = 0.;
      sign_node_y3 = 0.;
      break;
  }
}

void GaussIntegration::get_tri_rule_pts_and_weight()
{
  // get triangular rule integration points and weight
  switch (numberGaussPoints)
  {
    case 6:
      y1Area[0] = 0.09157621;
      y2Area[0] = 0.09157621;

      y1Area[1] = 0.09157621;
      y2Area[1] = 0.8168476;

      y1Area[2] = 0.8168476;
      y2Area[2] = 0.09157621;

      y1Area[3] = 0.4459485;
      y2Area[3] = 0.4459485;

      y1Area[4] = 0.4459485;
      y2Area[4] = 0.1081030;

      y1Area[5] = 0.1081030;
      y2Area[5] = 0.4459485;

      int i;
      for (i = 0; i < 3; i++)
      {
        totalGaussWeight[i] = 0.06348067;
      }
      for (i = 3; i < 6; i++)
      {
        totalGaussWeight[i] = 0.1289694;
      }
      break;
  }
}

void GaussIntegration::calculate_shape_function_2d_tri()
{
  int ife;
  double y1, y2, y3;
  get_tri_rule_pts_and_weight();

  for (ife = 0; ife < totalNumberGaussPts; ife++)
  {
    y1 = y1Area[ife];
    y2 = y2Area[ife];
    y3 = 1.0 - y1 - y2;

    shapeFunction[ife][0] = y1 * (2. * y1 - 1.);
    shapeFunction[ife][1] = y2 * (2. * y2 - 1.);
    shapeFunction[ife][2] = y3 * (2. * y3 - 1.);

    shapeFunction[ife][3] = 4. * y1 * y2;
    shapeFunction[ife][4] = 4. * y2 * y3;
    shapeFunction[ife][5] = 4. * y1 * y3;

    dndy1GaussPts[ife][0] = 4 * y1 - 1.;
    dndy1GaussPts[ife][1] = 0;
    dndy1GaussPts[ife][2] = 1 - 4. * y3;

    dndy1GaussPts[ife][3] = 4. * y2;
    dndy1GaussPts[ife][4] = -4. * y2;
    dndy1GaussPts[ife][5] = 4. * (1 - 2 * y1 - y2);

    dndy2GaussPts[ife][0] = 0.0;
    dndy2GaussPts[ife][1] = 4. * y2 - 1.;
    dndy2GaussPts[ife][2] = 1 - 4. * y3;

    dndy2GaussPts[ife][3] = 4. * y1;
    dndy2GaussPts[ife][4] = 4. * (1 - y1 - 2. * y2);
    dndy2GaussPts[ife][5] = -4. * y1;
  }
}

void GaussIntegration::calculate_derivative_at_nodes_2d_tri(
  double dndy1_at_nodes[][maxNumberNodes], double dndy2_at_nodes[][maxNumberNodes])
{
  double y1 = 0., y2 = 0., y3;
  int i;
  for (i = 0; i < numberNodes; i++)
  {
    switch (i)
    {
      case 0:
        y1 = 1.;
        y2 = 0.;
        break;
      case 1:
        y1 = 0.;
        y2 = 1.;
        break;
      case 2:
        y1 = 0.;
        y2 = 0.;
        break;
      case 3:
        y1 = 0.5;
        y2 = 0.5;
        break;
      case 4:
        y1 = 0.;
        y2 = 0.5;
        break;
      case 5:
        y1 = 0.5;
        y2 = 0.0;
        break;
    }

    y3 = 1. - y1 - y2;

    dndy1_at_nodes[i][0] = 4 * y1 - 1.;
    dndy1_at_nodes[i][1] = 0;
    dndy1_at_nodes[i][2] = 1 - 4. * y3;

    dndy1_at_nodes[i][3] = 4. * y2;
    dndy1_at_nodes[i][4] = -4. * y2;
    dndy1_at_nodes[i][5] = 4. * (1 - 2 * y1 - y2);

    dndy2_at_nodes[i][0] = 0.0;
    dndy2_at_nodes[i][1] = 4. * y2 - 1.;
    dndy2_at_nodes[i][2] = 1 - 4. * y3;

    dndy2_at_nodes[i][3] = 4. * y1;
    dndy2_at_nodes[i][4] = 4. * (1 - y1 - 2. * y2);
    dndy2_at_nodes[i][5] = -4. * y1;
  }
}

void GaussIntegration::get_tet_rule_pts_and_weight()
{
  // get tetrahedron rule integration points and weight
  double a, b;
  switch (numberGaussPoints)
  {
    case 1:
      // 1 integration point formula, degree of precision 1
      y1Volume[0] = 0.25;
      y2Volume[0] = 0.25;
      y3Volume[0] = 0.25;
      y4Volume[0] = 0.25;
      totalGaussWeight[0] = 1.;
      break;
    case 4:
      // 4 integration points formula, degree of precision 2
      a = 0.58541020;
      b = 0.13819660;

      y1Volume[0] = a;
      y2Volume[0] = b;
      y3Volume[0] = b;
      y4Volume[0] = b;

      y1Volume[1] = b;
      y2Volume[1] = a;
      y3Volume[1] = b;
      y4Volume[1] = b;

      y1Volume[2] = b;
      y2Volume[2] = b;
      y3Volume[2] = a;
      y4Volume[2] = b;

      y1Volume[3] = b;
      y2Volume[3] = b;
      y3Volume[3] = b;
      y4Volume[3] = a;

      int i;
      for (i = 0; i < 4; i++)
      {
        totalGaussWeight[i] = 0.25;
      }
      break;
  }
}

void GaussIntegration::calculate_shape_function_3d_tet()
{
  int ife;
  double y1, y2, y3, y4;
  get_tet_rule_pts_and_weight();

  switch (numberNodes)
  {
    case 10: // 10 nodes quadratic tet
    {
      for (ife = 0; ife < totalNumberGaussPts; ife++)
      {
        // y1,y2,y3,y4 are the volume coordinates
        y1 = y1Volume[ife];
        y2 = y2Volume[ife];
        y3 = y3Volume[ife];
        y4 = y4Volume[ife];

        // shape function is the same as in ABAQUS
        // it is different from that in all the FEA book
        // in which node is the first node
        // here at node 1 y4=1
        shapeFunction[ife][0] = y4 * (2. * y4 - 1.);
        shapeFunction[ife][1] = y1 * (2. * y1 - 1.);
        shapeFunction[ife][2] = y2 * (2. * y2 - 1.);
        shapeFunction[ife][3] = y3 * (2. * y3 - 1.);

        shapeFunction[ife][4] = 4. * y1 * y4;
        shapeFunction[ife][5] = 4. * y1 * y2;
        shapeFunction[ife][6] = 4. * y2 * y4;
        shapeFunction[ife][7] = 4. * y3 * y4;
        shapeFunction[ife][8] = 4. * y1 * y3;
        shapeFunction[ife][9] = 4. * y2 * y3;

        dndy1GaussPts[ife][0] = 1 - 4 * y4;
        dndy1GaussPts[ife][1] = 4 * y1 - 1.;
        dndy1GaussPts[ife][2] = 0;
        dndy1GaussPts[ife][3] = 0;

        dndy1GaussPts[ife][4] = 4. * (y4 - y1);
        dndy1GaussPts[ife][5] = 4. * y2;
        dndy1GaussPts[ife][6] = -4. * y2;
        dndy1GaussPts[ife][7] = -4. * y3;
        dndy1GaussPts[ife][8] = 4. * y3;
        dndy1GaussPts[ife][9] = 0;

        dndy2GaussPts[ife][0] = 1 - 4 * y4;
        dndy2GaussPts[ife][1] = 0;
        dndy2GaussPts[ife][2] = 4. * y2 - 1.;
        dndy2GaussPts[ife][3] = 0;

        dndy2GaussPts[ife][4] = -4. * y1;
        dndy2GaussPts[ife][5] = 4. * y1;
        dndy2GaussPts[ife][6] = 4. * (y4 - y2);
        dndy2GaussPts[ife][7] = -4. * y3;
        dndy2GaussPts[ife][8] = 0.;
        dndy2GaussPts[ife][9] = 4. * y3;

        dndy3GaussPts[ife][0] = 1 - 4 * y4;
        dndy3GaussPts[ife][1] = 0;
        dndy3GaussPts[ife][2] = 0;
        dndy3GaussPts[ife][3] = 4. * y3 - 1.;

        dndy3GaussPts[ife][4] = -4. * y1;
        dndy3GaussPts[ife][5] = 0;
        dndy3GaussPts[ife][6] = -4. * y2;
        dndy3GaussPts[ife][7] = 4. * (y4 - y3);
        dndy3GaussPts[ife][8] = 4. * y1;
        dndy3GaussPts[ife][9] = 4. * y2;
      }
      break;
    }
    case 4: // four node linear tet for debug purpose
    {
      for (ife = 0; ife < totalNumberGaussPts; ife++)
      {
        y1 = y1Volume[ife];
        y2 = y2Volume[ife];
        y3 = y3Volume[ife];
        y4 = y4Volume[ife];

        shapeFunction[ife][0] = y4;
        shapeFunction[ife][1] = y1;
        shapeFunction[ife][2] = y2;
        shapeFunction[ife][3] = y3;

        dndy1GaussPts[ife][0] = -1.;
        dndy1GaussPts[ife][1] = 1;
        dndy1GaussPts[ife][2] = 0;
        dndy1GaussPts[ife][3] = 0;

        dndy2GaussPts[ife][0] = -1.;
        dndy2GaussPts[ife][1] = 0;
        dndy2GaussPts[ife][2] = 1;
        dndy2GaussPts[ife][3] = 0;

        dndy3GaussPts[ife][0] = -1.;
        dndy3GaussPts[ife][1] = 0;
        dndy3GaussPts[ife][2] = 0;
        dndy3GaussPts[ife][3] = 1;
      }
      break;
    }
  }
}

void GaussIntegration::calculate_derivative_at_nodes_3d_tet(double dndy1_at_nodes[][maxNumberNodes],
  double dndy2_at_nodes[][maxNumberNodes], double dndy3_at_nodes[][maxNumberNodes])
{
  double y1, y2, y3, y4;
  int i;

  switch (numberNodes)
  {
    case 10:
    {
      for (i = 0; i < numberNodes; i++)
      {
        get_node_local_coord_tet(i, y1, y2, y3, y4);

        dndy1_at_nodes[i][0] = 1 - 4 * y4;
        dndy1_at_nodes[i][1] = 4 * y1 - 1.;
        dndy1_at_nodes[i][2] = 0;
        dndy1_at_nodes[i][3] = 0;

        dndy1_at_nodes[i][4] = 4. * (y4 - y1);
        dndy1_at_nodes[i][5] = 4. * y2;
        dndy1_at_nodes[i][6] = -4. * y2;
        dndy1_at_nodes[i][7] = -4. * y3;
        dndy1_at_nodes[i][8] = 4. * y3;
        dndy1_at_nodes[i][9] = 0;

        dndy2_at_nodes[i][0] = 1 - 4 * y4;
        dndy2_at_nodes[i][1] = 0;
        dndy2_at_nodes[i][2] = 4. * y2 - 1.;
        dndy2_at_nodes[i][3] = 0;
        dndy2_at_nodes[i][4] = -4. * y1;
        dndy2_at_nodes[i][5] = 4. * y1;
        dndy2_at_nodes[i][6] = 4. * (y4 - y2);
        dndy2_at_nodes[i][7] = -4. * y3;
        dndy2_at_nodes[i][8] = 0.;
        dndy2_at_nodes[i][9] = 4. * y3;

        dndy3_at_nodes[i][0] = 1 - 4 * y4;
        dndy3_at_nodes[i][1] = 0;
        dndy3_at_nodes[i][2] = 0;
        dndy3_at_nodes[i][3] = 4. * y3 - 1.;

        dndy3_at_nodes[i][4] = -4. * y1;
        dndy3_at_nodes[i][5] = 0;
        dndy3_at_nodes[i][6] = -4. * y2;
        dndy3_at_nodes[i][7] = 4. * (y4 - y3);
        dndy3_at_nodes[i][8] = 4. * y1;
        dndy3_at_nodes[i][9] = 4. * y2;
      }
      break;
    }
    case 4:
    {
      for (i = 0; i < numberNodes; i++)
      {
        get_node_local_coord_tet(i, y1, y2, y3, y4);
        dndy1_at_nodes[i][0] = -1.;
        dndy1_at_nodes[i][1] = 1;
        dndy1_at_nodes[i][2] = 0;
        dndy1_at_nodes[i][3] = 0;

        dndy2_at_nodes[i][0] = -1.;
        dndy2_at_nodes[i][1] = 0;
        dndy2_at_nodes[i][2] = 1;
        dndy2_at_nodes[i][3] = 0;

        dndy3_at_nodes[i][0] = -1.;
        dndy3_at_nodes[i][1] = 0;
        dndy3_at_nodes[i][2] = 0;
        dndy3_at_nodes[i][3] = 1;
      }
      break;
    }
  }
}

void GaussIntegration::get_node_local_coord_tet(
  int node_id, double& y1, double& y2, double& y3, double& y4)
{
  switch (node_id)
  {
    case 0:
      y1 = 0.;
      y2 = 0.;
      y3 = 0.;
      y4 = 1.;
      break;
    case 1:
      y1 = 1.;
      y2 = 0.;
      y3 = 0.;
      y4 = 0.;
      break;
    case 2:
      y1 = 0.;
      y2 = 1.;
      y3 = 0.;
      y4 = 0.;
      break;
    case 3:
      y1 = 0.;
      y2 = 0.;
      y3 = 1.;
      y4 = 0.;
      break;
    case 4:
      y1 = 0.5;
      y2 = 0.;
      y3 = 0.;
      y4 = 0.5;
      break;
    case 5:
      y1 = 0.5;
      y2 = 0.5;
      y3 = 0.;
      y4 = 0.;
      break;
    case 6:
      y1 = 0.;
      y2 = 0.5;
      y3 = 0.;
      y4 = 0.5;
      break;
    case 7:
      y1 = 0.;
      y2 = 0.0;
      y3 = 0.5;
      y4 = 0.5;
      break;
    case 8:
      y1 = 0.5;
      y2 = 0.;
      y3 = 0.5;
      y4 = 0.0;
      break;
    case 9:
      y1 = 0.;
      y2 = 0.5;
      y3 = 0.5;
      y4 = 0.;
      break;
    default:
      // Should not be possible to get here, but if we do, at least results will be consistent, not
      // random
      y1 = 0.;
      y2 = 0.;
      y3 = 0.;
      y4 = 0.;
      break;
  }
}
} // namespace verdict