FESADev/tests/unit/model/shell_geometry_test.cpp

#include "fesa/model/shell_geometry.h"

#include <gtest/gtest.h>

#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <limits>
#include <string>
#include <vector>

#include "fesa/math/vector3.h"

namespace {

using Vector3 = std::array<double, 3>;

fesa::Node Node(fesa::EntityIndex label, Vector3 coordinates) {
  return {{"Shell-Instance", static_cast<std::int64_t>(label),
           std::to_string(label)},
          coordinates,
          {"shell-geometry.inp", static_cast<std::size_t>(label + 1U)}};
}

fesa::Mitc4ShellDefinition Element(
    fesa::EntityIndex label, std::array<fesa::EntityIndex, 4> node_indices) {
  return {{"Shell-Instance", static_cast<std::int64_t>(label),
           std::to_string(label)},
          fesa::ShellSourceElementType::kS4,
          node_indices,
          0U,
          0U,
          {"shell-geometry.inp", static_cast<std::size_t>(100U + label)}};
}

std::vector<fesa::ShellSection> Sections(double thickness = 0.2) {
  return {{"Section", thickness, 0U, {"shell-geometry.inp", 90U}}};
}

double Dot(const Vector3& left, const Vector3& right) {
  return left[0] * right[0] + left[1] * right[1] + left[2] * right[2];
}

Vector3 Cross(const Vector3& left, const Vector3& right) {
  return {left[1] * right[2] - left[2] * right[1],
          left[2] * right[0] - left[0] * right[2],
          left[0] * right[1] - left[1] * right[0]};
}

double Norm(const Vector3& value) { return std::sqrt(Dot(value, value)); }

void ExpectVectorNear(const Vector3& actual, const Vector3& expected,
                      double tolerance = 1.0e-12) {
  for (std::size_t component = 0U; component < actual.size(); ++component) {
    EXPECT_NEAR(actual[component], expected[component], tolerance);
  }
}

void ExpectRightHandedFrame(const fesa::ShellNodeInitialFrame& frame) {
  EXPECT_NEAR(Norm(frame.director), 1.0, 1.0e-12);
  EXPECT_NEAR(Norm(frame.tangent_a), 1.0, 1.0e-12);
  EXPECT_NEAR(Norm(frame.tangent_b), 1.0, 1.0e-12);
  EXPECT_NEAR(Dot(frame.director, frame.tangent_a), 0.0, 1.0e-12);
  EXPECT_NEAR(Dot(frame.director, frame.tangent_b), 0.0, 1.0e-12);
  EXPECT_NEAR(Dot(frame.tangent_a, frame.tangent_b), 0.0, 1.0e-12);
  ExpectVectorNear(Cross(frame.tangent_a, frame.tangent_b), frame.director);
}

const fesa::ShellNodeInitialFrame& FrameFor(const fesa::ShellGeometry& geometry,
                                            fesa::EntityIndex node_index) {
  const auto found =
      std::find_if(geometry.nodal_frames.begin(), geometry.nodal_frames.end(),
                   [node_index](const fesa::ShellNodeInitialFrame& frame) {
                     return frame.node_index == node_index;
                   });
  EXPECT_NE(found, geometry.nodal_frames.end());
  return *found;
}

void ExpectFailureCode(const fesa::Result<fesa::ShellGeometry>& result,
                       const std::string& code) {
  ASSERT_FALSE(result.HasValue());
  EXPECT_EQ(result.GetStatus().Category(), fesa::FailureCategory::kModel);
  ASSERT_EQ(result.GetStatus().Diagnostics().size(), 1U);
  EXPECT_EQ(result.GetStatus().Diagnostics()[0].code, code);
}

}  // namespace

// MITC4-GEO-001
TEST(Mitc4Geometry,
     BuildsDeterministicFramesForPlanarRotatedAndWarpedElements) {
  const std::vector<fesa::Node> planar_nodes{
      Node(0U, {0.0, 0.0, 0.0}), Node(1U, {1.0, 0.0, 0.0}),
      Node(2U, {1.0, 1.0, 0.0}), Node(3U, {0.0, 1.0, 0.0})};
  auto planar = fesa::PreprocessShellGeometry(
      planar_nodes, {Element(10U, {0U, 1U, 2U, 3U})}, Sections());

  ASSERT_TRUE(planar.HasValue());
  ASSERT_EQ(planar.Value().element_data.size(), 1U);
  ExpectVectorNear(planar.Value().element_data[0].normal_candidate,
                   {0.0, 0.0, 1.0});
  EXPECT_NEAR(planar.Value().element_data[0].surface_area_weight, 1.0, 1.0e-12);
  ASSERT_EQ(planar.Value().nodal_frames.size(), 4U);
  for (const auto& frame : planar.Value().nodal_frames) {
    ExpectVectorNear(frame.director, {0.0, 0.0, 1.0});
    ExpectVectorNear(frame.tangent_a, {1.0, 0.0, 0.0});
    ExpectVectorNear(frame.tangent_b, {0.0, 1.0, 0.0});
    ExpectRightHandedFrame(frame);
  }

  const std::vector<fesa::Node> rotated_nodes{
      Node(0U, {0.0, 0.0, 0.0}), Node(1U, {0.0, 1.0, 0.0}),
      Node(2U, {0.0, 1.0, 1.0}), Node(3U, {0.0, 0.0, 1.0})};
  auto rotated = fesa::PreprocessShellGeometry(
      rotated_nodes, {Element(11U, {0U, 1U, 2U, 3U})}, Sections());

  ASSERT_TRUE(rotated.HasValue());
  const auto& rotated_frame = FrameFor(rotated.Value(), 0U);
  ExpectVectorNear(rotated_frame.director, {1.0, 0.0, 0.0});
  ExpectVectorNear(rotated_frame.tangent_a, {0.0, 1.0, 0.0});
  ExpectVectorNear(rotated_frame.tangent_b, {0.0, 0.0, 1.0});
  ExpectRightHandedFrame(rotated_frame);

  const std::vector<fesa::Node> warped_nodes{
      Node(0U, {0.0, 0.0, 0.0}), Node(1U, {2.0, 0.0, 0.0}),
      Node(2U, {2.0, 1.0, 0.2}), Node(3U, {0.0, 1.0, 0.0})};
  auto warped = fesa::PreprocessShellGeometry(
      warped_nodes, {Element(12U, {0U, 1U, 2U, 3U})}, Sections());

  ASSERT_TRUE(warped.HasValue());
  EXPECT_GT(warped.Value().element_data[0].surface_area_weight, 2.0);
  for (const auto& frame : warped.Value().nodal_frames) {
    ExpectRightHandedFrame(frame);
  }
}

TEST(Mitc4Geometry, PreservesExactWarpedResultsAcrossVector3Migration) {
  const std::vector<fesa::Node> warped_nodes{
      Node(0U, {0.0, 0.0, 0.0}), Node(1U, {2.0, 0.0, 0.0}),
      Node(2U, {2.0, 1.0, 0.2}), Node(3U, {0.0, 1.0, 0.0})};
  const auto warped = fesa::PreprocessShellGeometry(
      warped_nodes, {Element(12U, {0U, 1U, 2U, 3U})}, Sections());
  ASSERT_TRUE(warped.HasValue());
  const auto& element = warped.Value().element_data[0U];
  const auto& frame = FrameFor(warped.Value(), 0U);
  const fesa::Vector3 normal{element.normal_candidate};
  const fesa::Vector3 director{frame.director};
  const fesa::Vector3 tangent_a{frame.tangent_a};
  const fesa::Vector3 tangent_b{frame.tangent_b};
  EXPECT_DOUBLE_EQ(normal.X(), -0x1.97105218e28c2p-5);
  EXPECT_DOUBLE_EQ(normal.Y(), -0x1.97105218e28c2p-4);
  EXPECT_DOUBLE_EQ(normal.Z(), 0x1.fcd4669f1b2f2p-1);
  EXPECT_DOUBLE_EQ(element.surface_area_weight, 0x1.021ebe8f40622p+1);
  EXPECT_DOUBLE_EQ(director.X(), -0x1.97105218e28c2p-5);
  EXPECT_DOUBLE_EQ(director.Y(), -0x1.97105218e28c2p-4);
  EXPECT_DOUBLE_EQ(director.Z(), 0x1.fcd4669f1b2f2p-1);
  EXPECT_DOUBLE_EQ(tangent_a.X(), 0x1.ff5e152c2d2abp-1);
  EXPECT_DOUBLE_EQ(tangent_a.Y(), -0x1.4408ec773d925p-8);
  EXPECT_DOUBLE_EQ(tangent_a.Z(), 0x1.950b27950cf6dp-5);
  EXPECT_DOUBLE_EQ(tangent_b.X(), 0.0);
  EXPECT_DOUBLE_EQ(tangent_b.Y(), 0x1.fd7583bc82e28p-1);
  EXPECT_DOUBLE_EQ(tangent_b.Z(), 0x1.9791363068b53p-4);
}

// MITC4-GEO-002
TEST(Mitc4Geometry, AreaWeightsSharedDirectorsInStableSourceIdentityOrder) {
  const std::vector<fesa::Node> nodes{
      Node(0U, {0.0, 0.0, 0.0}), Node(1U, {1.0, 0.0, 0.0}),
      Node(2U, {1.0, 1.0, 0.0}), Node(3U, {0.0, 1.0, 0.0}),
      Node(4U, {2.0, 0.0, 1.0}), Node(5U, {2.0, 1.0, 1.0})};
  const auto flat = Element(10U, {0U, 1U, 2U, 3U});
  const auto tilted = Element(20U, {1U, 4U, 5U, 2U});

  auto first = fesa::PreprocessShellGeometry(nodes, {tilted, flat}, Sections());
  auto second =
      fesa::PreprocessShellGeometry(nodes, {flat, tilted}, Sections());

  ASSERT_TRUE(first.HasValue());
  ASSERT_TRUE(second.HasValue());
  const Vector3 expected_shared_director{-1.0 / std::sqrt(5.0), 0.0,
                                         2.0 / std::sqrt(5.0)};
  for (const auto sharedNode : {1U, 2U}) {
    const auto& firstFrame = FrameFor(first.Value(), sharedNode);
    const auto& secondFrame = FrameFor(second.Value(), sharedNode);
    ExpectVectorNear(firstFrame.director, expected_shared_director);
    ExpectVectorNear(firstFrame.director, secondFrame.director, 0.0);
    ExpectVectorNear(firstFrame.tangent_a, {0.0, 1.0, 0.0});
    ExpectRightHandedFrame(firstFrame);
  }
}

// MITC4-GEO-003
TEST(Mitc4Geometry, RejectsInvalidSurfaceJacobianAndIncidentOrientationCases) {
  const auto valid_element = Element(10U, {0U, 1U, 2U, 3U});

  ExpectFailureCode(fesa::PreprocessShellGeometry(
                        {Node(0U, {0.0, 0.0, 0.0}), Node(1U, {0.0, 0.0, 0.0}),
                         Node(2U, {1.0, 1.0, 0.0}), Node(3U, {0.0, 1.0, 0.0})},
                        {valid_element}, Sections()),
                    "invalid-shell-geometry");

  ExpectFailureCode(fesa::PreprocessShellGeometry(
                        {Node(0U, {0.0, 0.0, 0.0}), Node(1U, {1.0, 1.0, 0.0}),
                         Node(2U, {0.0, 1.0, 0.0}), Node(3U, {1.0, 0.0, 0.0})},
                        {valid_element}, Sections()),
                    "invalid-shell-geometry");

  ExpectFailureCode(fesa::PreprocessShellGeometry(
                        {Node(0U, {0.0, 0.0, 0.0}), Node(1U, {1.0, 0.0, 0.0}),
                         Node(2U, {2.0, 0.0, 0.0}), Node(3U, {3.0, 0.0, 0.0})},
                        {valid_element}, Sections()),
                    "invalid-shell-geometry");

  ExpectFailureCode(
      fesa::PreprocessShellGeometry(
          {Node(0U, {0.0, 0.0, 0.0}), Node(1U, {1.0, 0.0, 0.0}),
           Node(2U, {0.05, 0.05, 0.0}), Node(3U, {0.0, 1.0, 0.0})},
          {valid_element}, Sections()),
      "invalid-shell-geometry");

  ExpectFailureCode(
      fesa::PreprocessShellGeometry(
          {Node(0U, {0.0, 0.0, 0.0}), Node(1U, {1.0, 0.0, 0.0}),
           Node(2U, {1.0, std::numeric_limits<double>::quiet_NaN(), 0.0}),
           Node(3U, {0.0, 1.0, 0.0})},
          {valid_element}, Sections()),
      "invalid-shell-geometry");

  ExpectFailureCode(fesa::PreprocessShellGeometry(
                        {Node(0U, {0.0, 0.0, 0.0}), Node(1U, {1.0, 0.0, 0.0}),
                         Node(2U, {1.0, 1.0, 0.0}), Node(3U, {0.0, 1.0, 0.0})},
                        {valid_element}, Sections(0.0)),
                    "invalid-shell-jacobian");

  const std::vector<fesa::Node> opposed_nodes{
      Node(0U, {0.0, 0.0, 0.0}), Node(1U, {1.0, 0.0, 0.0}),
      Node(2U, {1.0, 1.0, 0.0}), Node(3U, {0.0, 1.0, 0.0})};
  ExpectFailureCode(
      fesa::PreprocessShellGeometry(
          opposed_nodes,
          {Element(10U, {0U, 1U, 2U, 3U}), Element(20U, {0U, 3U, 2U, 1U})},
          Sections()),
      "opposed-incident-normal");
}

// MITC4-GEO-004
TEST(Mitc4Geometry, ExposesTheCompleteRequiredValidationPointInventory) {
  const auto& points = fesa::ShellGeometryValidationPoints();
  ASSERT_EQ(points.size(), 17U);
  EXPECT_EQ(std::count_if(points.begin(), points.end(),
                          [](const auto& point) {
                            return point.kind ==
                                   fesa::ShellGeometryPointKind::kCenter;
                          }),
            1);
  EXPECT_EQ(std::count_if(points.begin(), points.end(),
                          [](const auto& point) {
                            return point.kind ==
                                   fesa::ShellGeometryPointKind::kStiffness;
                          }),
            8);
  EXPECT_EQ(std::count_if(points.begin(), points.end(),
                          [](const auto& point) {
                            return point.kind ==
                                   fesa::ShellGeometryPointKind::kTying;
                          }),
            4);
  EXPECT_EQ(std::count_if(points.begin(), points.end(),
                          [](const auto& point) {
                            return point.kind ==
                                   fesa::ShellGeometryPointKind::kRecovery;
                          }),
            4);

  EXPECT_EQ(points.front().natural_coordinates, (Vector3{0.0, 0.0, 0.0}));
  const double gauss = 1.0 / std::sqrt(3.0);
  EXPECT_EQ(points[1].natural_coordinates, (Vector3{-gauss, -gauss, -gauss}));
  EXPECT_EQ(points[8].natural_coordinates, (Vector3{-gauss, gauss, gauss}));
  EXPECT_EQ(points[9].natural_coordinates, (Vector3{0.0, -1.0, 0.0}));
  EXPECT_EQ(points[12].natural_coordinates, (Vector3{1.0, 0.0, 0.0}));
  EXPECT_EQ(points[13].natural_coordinates, (Vector3{-gauss, -gauss, 0.0}));
  EXPECT_EQ(points[16].natural_coordinates, (Vector3{-gauss, gauss, 0.0}));
}