Detect inconsistent fiber entry and throw an error

doc/documentation/src/analysis_guide_templates/geometrydefinition.rst.j2

@@ -53,9 +53,9 @@ The nodes are read in by the section ``NODE COORDS``.
 
 Within this section, the nodes are read one per line. There are different types of nodes:
 
-- General nodes (NODE)
-- Control points for nurbs geometry (CP)
-- Nodes with addition fiber information (FNODE)
+- General nodes (`NODE`)
+- Control points for nurbs geometry (`CP`)
+- Nodes with addition fiber information  (`FNODE`)
 
 All node coordinates must be given with three coordinates, even if the problem is 1D or 2D:
 
@@ -64,8 +64,14 @@ All node coordinates must be given with three coordinates, even if the problem i
    NODE COORDS:
    - "NODE <ID> COORD <coord-x> <coord-y> <coord-z> [ROTANGLE <coord-yz> <coord-xz> <coord-xy>]"
    - "CP <ID> COORD <coord-x> <coord-y> <coord-z> <weight>"
-   - "FNODE <ID> COORD <coord-x> <coord-y> <coord-z> [FIBER1|FIBER2|FIBER3|CIR|TAN|RAD|HELIX|TRANS] <further parameters>"
+   - "FNODE <ID> COORD <coord-x> <coord-y> <coord-z> FIBER1 <x y z> [FIBER2 <x y z> [FIBER3 <x y z>]]"
+   - "FNODE <ID> "COORD <coord-x> <coord-y> <coord-z> CIR <x y z> TAN <x y z> HELIX <angle> TRANS <angle>"
 
+For ``FNODE`` in combination with fiber directions, the directions must be numbered consecutively starting with ``FIBER1``.
+For the problemtype `Cardiac_Monodomain`, the fiber directions can be entered alternatively
+using ``CIR``, ``TAN``, ``HELIX``, and ``TRANS``.
+All nodes of an element must provide the same set of parameters, either explicit or cardiac fiber directions.
+Their direction vectors and angle values may vary between nodes.
 
 .. _geometrysets:
 
@@ -177,4 +183,3 @@ while the parameter :math:`\sigma` is defined by the parameter ``RADIUSDISTRIBUT
     PARTICLE DYNAMIC/DEM:
       INITIAL_RADIUS: RadiusFromParticleInput|NormalRadiusDistribution|LogNormalRadiusDistribution
       RADIUSDISTRIBUTION_SIGMA: <variation>   # variation sigma for a radius distribution
-

src/core/fem/src/general/node/4C_fem_general_fiber_node_holder.cpp

@@ -9,6 +9,7 @@
 
 #include "4C_comm_parobject.hpp"
 #include "4C_fem_general_fiber_node.hpp"
+#include "4C_utils_exceptions.hpp"
 
 FOUR_C_NAMESPACE_OPEN
 
@@ -60,7 +61,13 @@ void Core::Nodes::NodalFiberHolder::set_angle(AngleType type, const std::vector<
 const std::vector<double>& Core::Nodes::NodalFiberHolder::get_angle(
     Core::Nodes::AngleType type) const
 {
-  return angles_.at(type);
+  const auto angle = angles_.find(type);
+  if (angle == angles_.end())
+  {
+    FOUR_C_THROW("Nodal fiber data does not contain the requested {} angle.",
+        type == AngleType::Helix ? "HELIX" : "TRANS");
+  }
+  return angle->second;
 }
 
 std::size_t Core::Nodes::NodalFiberHolder::fibers_size() const { return fibers_.size(); }

src/core/fem/src/general/node/4C_fem_general_fiber_node_utils.cpp

@@ -36,6 +36,43 @@ void Core::Nodes::project_fibers_to_gauss_points(const Core::Nodes::Node* const*
       FOUR_C_THROW("At least one node of the element does not provide fibers.");
     }
 
+    if (inode > 0)
+    {
+      if (fiberNodes[inode]->fibers().size() != fiberNodes[0]->fibers().size())
+      {
+        FOUR_C_THROW(
+            "All nodes of an element must define the same number of FIBER directions. Node {} "
+            "defines {}, but node {} defines {}.",
+            fiberNodes[inode]->id(), fiberNodes[inode]->fibers().size(), fiberNodes[0]->id(),
+            fiberNodes[0]->fibers().size());
+      }
+
+      const auto& reference_directions = fiberNodes[0]->coordinate_system_directions();
+      const auto& node_directions = fiberNodes[inode]->coordinate_system_directions();
+      const bool same_coordinate_directions =
+          reference_directions.size() == node_directions.size() &&
+          std::all_of(reference_directions.begin(), reference_directions.end(),
+              [&node_directions](const auto& direction)
+              { return node_directions.contains(direction.first); });
+      if (!same_coordinate_directions)
+      {
+        FOUR_C_THROW(
+            "All nodes of an element must define the same set of coordinate system direction "
+            "types (CIR and TAN).");
+      }
+
+      const auto& reference_angles = fiberNodes[0]->angles();
+      const auto& node_angles = fiberNodes[inode]->angles();
+      const bool same_angles =
+          reference_angles.size() == node_angles.size() &&
+          std::all_of(reference_angles.begin(), reference_angles.end(),
+              [&node_angles](const auto& angle) { return node_angles.contains(angle.first); });
+      if (!same_angles)
+      {
+        FOUR_C_THROW("All nodes of an element must define the same angle types.");
+      }
+    }
+
     for (const auto& pair : fiberNodes[inode]->coordinate_system_directions())
     {
       coordinateSystemDirections[pair.first][inode] = pair.second;
@@ -120,21 +157,6 @@ void Core::Nodes::project_fibers_to_gauss_points(const Core::Nodes::Node* const*
       tan[gp] += -tancir * cir[gp];
       tan[gp] *= 1.0 / Core::LinAlg::norm2(tan[gp]);
     }
-
-    // orthogonalize radial vector, preserve circular and tangential direction
-    if (gpFiberHolder.contains_coordinate_system_direction(CoordinateSystemDirection::Radial))
-    {
-      std::vector<Core::LinAlg::Tensor<double, 3>>& rad =
-          gpFiberHolder.get_coordinate_system_direction_mutual(CoordinateSystemDirection::Radial);
-      for (std::size_t gp = 0; gp < tan.size(); ++gp)
-      {
-        double radcir = rad[gp] * cir[gp];
-        double radtan = rad[gp] * tan[gp];
-        // double
-        rad[gp] += -radcir * cir[gp] - radtan * tan[gp];
-        rad[gp] *= 1.0 / Core::LinAlg::norm2(rad[gp]);
-      }
-    }
   }
 }
 

src/core/io/src/4C_io_meshreader.cpp

@@ -66,6 +66,74 @@ namespace Core::IO::Internal
      */
     std::optional<Core::IO::MeshInput::Mesh<3>> filtered_mesh_on_rank_zero{};
   };
+
+  FiberNodeData read_fiber_node_data(Core::IO::ValueParser& parser)
+  {
+    FiberNodeData data;
+
+    while (!parser.at_end())
+    {
+      const auto next = parser.read<std::string>();
+
+      if (next.starts_with("FIBER"))
+      {
+        const std::string expected = "FIBER" + std::to_string(data.fibers.size() + 1);
+        if (next != expected)
+        {
+          FOUR_C_THROW(
+              "Fiber directions in FNODE must be numbered consecutively. Expected '{}', found "
+              "'{}'.",
+              expected, next);
+        }
+        data.fibers.emplace_back(parser.read<std::array<double, 3>>());
+      }
+      else if (next == "CIR")
+      {
+        data.coordinate_system_directions[Core::Nodes::CoordinateSystemDirection::Circular] =
+            parser.read<std::array<double, 3>>();
+      }
+      else if (next == "TAN")
+      {
+        data.coordinate_system_directions[Core::Nodes::CoordinateSystemDirection::Tangential] =
+            parser.read<std::array<double, 3>>();
+      }
+      else if (next == "HELIX")
+      {
+        data.angles[Core::Nodes::AngleType::Helix] = parser.read<double>();
+      }
+      else if (next == "TRANS")
+      {
+        data.angles[Core::Nodes::AngleType::Transverse] = parser.read<double>();
+      }
+      else
+      {
+        FOUR_C_THROW("Unknown FNODE parameter '{}'.", next);
+      }
+    }
+
+    const bool has_circular = data.coordinate_system_directions.contains(
+        Core::Nodes::CoordinateSystemDirection::Circular);
+    const bool has_tangential = data.coordinate_system_directions.contains(
+        Core::Nodes::CoordinateSystemDirection::Tangential);
+    const bool has_helix = data.angles.contains(Core::Nodes::AngleType::Helix);
+    const bool has_transverse = data.angles.contains(Core::Nodes::AngleType::Transverse);
+    const bool has_cardiac_directions =
+        has_circular || has_tangential || has_helix || has_transverse;
+
+    if (!data.fibers.empty() && has_cardiac_directions)
+    {
+      FOUR_C_THROW(
+          "Fiber directions in FNODE must be defined either by CIR/TAN/HELIX/TRANS or by FIBER1, "
+          "FIBER2, etc., but not by both.");
+    }
+
+    if (has_cardiac_directions && !(has_circular && has_tangential && has_helix && has_transverse))
+    {
+      FOUR_C_THROW("Cardiac fiber directions in FNODE require all of CIR, TAN, HELIX, and TRANS.");
+    }
+
+    return data;
+  }
 }  // namespace Core::IO::Internal
 
 namespace
@@ -453,65 +521,21 @@ namespace
       // this is a special node with additional fiber information
       else if (type == "FNODE")
       {
-        enum class FiberType
-        {
-          Unknown,
-          Angle,
-          Fiber,
-          CosyDirection
-        };
-
-        // read fiber node
-        std::map<Core::Nodes::CoordinateSystemDirection, std::array<double, 3>> cosyDirections;
-        std::vector<std::array<double, 3>> fibers;
-        std::map<Core::Nodes::AngleType, double> angles;
-
         int nodeid = parser.read<int>() - 1;
         parser.consume("COORD");
         auto coords = parser.read<std::vector<double>>(3);
         max_node_id = std::max(max_node_id, nodeid) + 1;
-
-        while (!parser.at_end())
-        {
-          auto next = parser.read<std::string>();
-
-          if (next == "FIBER" + std::to_string(1 + fibers.size()))
-          {
-            fibers.emplace_back(parser.read<std::array<double, 3>>());
-          }
-          else if (next == "CIR")
-          {
-            cosyDirections[Core::Nodes::CoordinateSystemDirection::Circular] =
-                parser.read<std::array<double, 3>>();
-          }
-          else if (next == "TAN")
-          {
-            cosyDirections[Core::Nodes::CoordinateSystemDirection::Tangential] =
-                parser.read<std::array<double, 3>>();
-          }
-          else if (next == "RAD")
-          {
-            cosyDirections[Core::Nodes::CoordinateSystemDirection::Radial] =
-                parser.read<std::array<double, 3>>();
-          }
-          else if (next == "HELIX")
-          {
-            angles[Core::Nodes::AngleType::Helix] = parser.read<double>();
-          }
-          else if (next == "TRANS")
-          {
-            angles[Core::Nodes::AngleType::Transverse] = parser.read<double>();
-          }
-        }
+        auto fiber_node_data = Core::IO::Internal::read_fiber_node_data(parser);
 
         // add fiber information to node
         std::vector<std::shared_ptr<Core::FE::Discretization>> discretizations =
             find_dis_node(element_readers, nodeid);
         for (auto& dis : discretizations)
         {
           sanitize_node_coordinates(nodeid, dis->n_dim(), coords);
-          auto node = std::make_shared<Core::Nodes::FiberNode>(
-              nodeid, coords, cosyDirections, fibers, angles, myrank);
+          auto node = std::make_shared<Core::Nodes::FiberNode>(nodeid, coords,
+              fiber_node_data.coordinate_system_directions, fiber_node_data.fibers,
+              fiber_node_data.angles, myrank);
           dis->add_node(coords, nodeid, node);
         }
       }

src/core/io/src/4C_io_meshreader.hpp

@@ -10,11 +10,13 @@
 
 #include "4C_config.hpp"
 
+#include "4C_fem_general_fiber_node.hpp"
 #include "4C_linalg_graph.hpp"
 #include "4C_rebalance.hpp"
 
 #include <Teuchos_ParameterList.hpp>
 
+#include <array>
 #include <map>
 #include <string>
 #include <vector>
@@ -29,6 +31,7 @@ namespace Core::FE
 namespace Core::IO
 {
   class InputFile;
+  class ValueParser;
 
   namespace MeshInput
   {
@@ -39,7 +42,17 @@ namespace Core::IO
   namespace Internal
   {
     struct MeshReader;
-  }
+
+    struct FiberNodeData
+    {
+      std::map<Core::Nodes::CoordinateSystemDirection, std::array<double, 3>>
+          coordinate_system_directions;
+      std::vector<std::array<double, 3>> fibers;
+      std::map<Core::Nodes::AngleType, double> angles;
+    };
+
+    FiberNodeData read_fiber_node_data(Core::IO::ValueParser& parser);
+  }  // namespace Internal
 
 
   /**

src/solid_ele/4C_solid_ele_fiber.cpp

@@ -14,11 +14,45 @@
 #include <algorithm>
 
 FOUR_C_NAMESPACE_OPEN
+namespace
+{
+  bool has_nonempty_optional_vector_parameter(
+      const Core::IO::InputParameterContainer& input_data, const std::string& name)
+  {
+    const auto* parameter = input_data.get_if<std::optional<std::vector<double>>>(name);
+    return parameter != nullptr && parameter->has_value() && !parameter->value().empty();
+  }
 
+  bool has_nonempty_optional_vector(const std::optional<std::vector<double>>* parameter)
+  {
+    return parameter != nullptr && parameter->has_value() && !parameter->value().empty();
+  }
+
+  void validate_direction_input(const Core::IO::InputParameterContainer& input_data)
+  {
+    const bool has_coordinate_system = has_nonempty_optional_vector_parameter(input_data, "RAD") ||
+                                       has_nonempty_optional_vector_parameter(input_data, "AXI") ||
+                                       has_nonempty_optional_vector_parameter(input_data, "CIR");
+    const bool has_fiber1 = has_nonempty_optional_vector_parameter(input_data, "FIBER1");
+    const bool has_fiber2 = has_nonempty_optional_vector_parameter(input_data, "FIBER2");
+    const bool has_fiber3 = has_nonempty_optional_vector_parameter(input_data, "FIBER3");
+    const bool has_fibers = has_fiber1 || has_fiber2 || has_fiber3;
+
+    FOUR_C_ASSERT_ALWAYS(!(has_coordinate_system && has_fibers),
+        "Material directions must be specified either by RAD/AXI/CIR or by FIBER1, FIBER2, "
+        "FIBER3, but not by both.");
+    FOUR_C_ASSERT_ALWAYS(!has_fiber2 || has_fiber1,
+        "Fiber directions must be numbered consecutively: FIBER2 requires FIBER1.");
+    FOUR_C_ASSERT_ALWAYS(!has_fiber3 || has_fiber2,
+        "Fiber directions must be numbered consecutively: FIBER3 requires FIBER2.");
+  }
+}  // namespace
 
 Discret::Elements::Fibers Discret::Elements::read_fibers(
     const Core::IO::InputParameterContainer& input_data)
 {
+  validate_direction_input(input_data);
+
   Fibers fibers{};
 
   // for now only support the old style fiber input via FIBER1, FIBER2, ... keywords
@@ -27,7 +61,7 @@ Discret::Elements::Fibers Discret::Elements::read_fibers(
     const std::string fiber_name = "FIBER" + std::to_string(i);
 
     const auto* fiber_ptr = input_data.get_if<std::optional<std::vector<double>>>(fiber_name);
-    if (!fiber_ptr || !fiber_ptr->has_value())
+    if (!has_nonempty_optional_vector(fiber_ptr))
     {
       break;
     }
@@ -38,26 +72,26 @@ Discret::Elements::Fibers Discret::Elements::read_fibers(
 
     fibers.element_fibers.emplace_back(tensor);
   }
-
   return fibers;
 }
 
 std::optional<Discret::Elements::CoordinateSystem> Discret::Elements::read_coordinate_system(
     const Core::IO::InputParameterContainer& input_data)
 {
+  validate_direction_input(input_data);
+
   // for now only support the old style fiber input via RAD, AXI, CIR keywords
   const std::array rad_axi_cir = {input_data.get_if<std::optional<std::vector<double>>>("RAD"),
       input_data.get_if<std::optional<std::vector<double>>>("AXI"),
       input_data.get_if<std::optional<std::vector<double>>>("CIR")};
 
-  if (std::ranges::none_of(rad_axi_cir, [](const auto* p) { return p && p->has_value(); }))
+  if (std::ranges::none_of(rad_axi_cir, has_nonempty_optional_vector))
   {
     // no local coordinate system defined
     return std::nullopt;
   }
 
-  FOUR_C_ASSERT_ALWAYS(
-      std::ranges::all_of(rad_axi_cir, [](const auto* p) { return p && p->has_value(); }),
+  FOUR_C_ASSERT_ALWAYS(std::ranges::all_of(rad_axi_cir, has_nonempty_optional_vector),
       "If you specify a coordinate system, you need to define all of RAD, AXI and CIR!");
 
   CoordinateSystem coord_sys{};