diff --git a/src/beamme/core/base_mesh_item.py b/src/beamme/core/base_mesh_item.py
index c87db9fb..63daa604 100644
--- a/src/beamme/core/base_mesh_item.py
+++ b/src/beamme/core/base_mesh_item.py
@@ -20,15 +20,14 @@
 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 # THE SOFTWARE.
 """This module implements the class that will be used as the base for all items
-that are in a mesh."""
-
-from typing import Optional as _Optional
+that are in a mesh, except for nodes and elements."""
 
 
 class BaseMeshItem:
-    """Base class for all objects that are related to a mesh."""
+    """Base class for boundary conditions, functions, geometry sets and
+    materials."""
 
-    def __init__(self, data: _Optional[dict] = None):
+    def __init__(self, data: dict | None = None):
         """Create the base object.
 
         Args:
diff --git a/src/beamme/core/conf.py b/src/beamme/core/conf.py
index c4b579a2..7d7d8c29 100644
--- a/src/beamme/core/conf.py
+++ b/src/beamme/core/conf.py
@@ -81,6 +81,14 @@ class CouplingDofType(_Enum):
     joint = _auto()
 
 
+class ElementType(_Enum):
+    """Enum for element types."""
+
+    beam = _auto()
+    nurbs = _auto()
+    solid = _auto()
+
+
 class DoubleNodes(_Enum):
     """Enum for handing double nodes in Neumann conditions."""
 
@@ -93,17 +101,10 @@ class BeamMe(object):
 
     def __init__(self):
         self.set_default_values()
-
-        # Geometry types.
         self.geo = Geometry
-
-        # Boundary conditions types.
         self.bc = BoundaryCondition
-
-        # Coupling types.
         self.coupling_dof = CouplingDofType
-
-        # Handling of multiple nodes in Neumann bcs.
+        self.element_type = ElementType
         self.double_nodes = DoubleNodes
 
     def set_default_values(self):
diff --git a/src/beamme/core/element.py b/src/beamme/core/element.py
index 32a1b090..3998e929 100644
--- a/src/beamme/core/element.py
+++ b/src/beamme/core/element.py
@@ -21,14 +21,13 @@
 # THE SOFTWARE.
 """This module implements the class that represents one element in the Mesh."""
 
-from beamme.core.base_mesh_item import BaseMeshItem as _BaseMeshItem
 
-
-class Element(_BaseMeshItem):
+class Element:
     """A base class for an FEM element in the mesh."""
 
-    def __init__(self, nodes=None, material=None, **kwargs):
-        super().__init__(**kwargs)
+    def __init__(self, nodes=None, material=None):
+        # Global index of this item in a mesh.
+        self.i_global: None | int = None
 
         # List of nodes that are connected to the element.
         if nodes is None:
diff --git a/src/beamme/core/node.py b/src/beamme/core/node.py
index fd3053ae..6bf26643 100644
--- a/src/beamme/core/node.py
+++ b/src/beamme/core/node.py
@@ -23,15 +23,15 @@
 
 import numpy as _np
 
-from beamme.core.base_mesh_item import BaseMeshItem as _BaseMeshItem
 from beamme.core.rotation import Rotation as _Rotation
 
 
-class Node(_BaseMeshItem):
+class Node:
     """This object represents one node in the mesh."""
 
-    def __init__(self, coordinates, *, is_middle_node=False, **kwargs):
-        super().__init__(**kwargs)
+    def __init__(self, coordinates, *, is_middle_node=False):
+        # Global index of this item in a mesh.
+        self.i_global: None | int = None
 
         # Coordinates of this node.
         self.coordinates = _np.array(coordinates, dtype=float)
diff --git a/src/beamme/core/nurbs_patch.py b/src/beamme/core/nurbs_patch.py
index c3f83cfb..93918e8e 100644
--- a/src/beamme/core/nurbs_patch.py
+++ b/src/beamme/core/nurbs_patch.py
@@ -39,15 +39,8 @@ class NURBSPatch(_Element):
     # A list of valid material types for this element
     valid_materials = [_MaterialSolidBase]
 
-    def __init__(
-        self,
-        knot_vectors,
-        polynomial_orders,
-        material=None,
-        nodes=None,
-        data=None,
-    ):
-        super().__init__(nodes=nodes, material=material, data=data)
+    def __init__(self, knot_vectors, polynomial_orders, material=None, nodes=None):
+        super().__init__(nodes=nodes, material=material)
 
         # Knot vectors
         self.knot_vectors = knot_vectors
@@ -55,7 +48,8 @@ def __init__(
         # Polynomial degrees
         self.polynomial_orders = polynomial_orders
 
-        # Set numbers for elements
+        # Global indices
+        self.i_global_start = None
         self.i_nurbs_patch = None
 
     def get_nurbs_dimension(self) -> int:
diff --git a/src/beamme/four_c/element_beam.py b/src/beamme/four_c/element_beam.py
index c9bbbd6e..39b08ebb 100644
--- a/src/beamme/four_c/element_beam.py
+++ b/src/beamme/four_c/element_beam.py
@@ -22,11 +22,13 @@
 """This file implements beam elements for 4C."""
 
 import warnings as _warnings
+from typing import Any as _Any
 
 import numpy as _np
 
 from beamme.core.conf import bme as _bme
 from beamme.core.element_beam import Beam as _Beam
+from beamme.core.element_beam import Beam2 as _Beam2
 from beamme.core.element_beam import generate_beam_class as _generate_beam_class
 from beamme.four_c.four_c_types import (
     BeamKirchhoffConstraintType as _BeamKirchhoffConstraintType,
@@ -35,6 +37,9 @@
     BeamKirchhoffParametrizationType as _BeamKirchhoffParametrizationType,
 )
 from beamme.four_c.four_c_types import BeamType as _BeamType
+from beamme.four_c.input_file_dump_item import (
+    get_four_c_legacy_dict as _get_four_c_legacy_dict,
+)
 from beamme.four_c.input_file_mappings import (
     INPUT_FILE_MAPPINGS as _INPUT_FILE_MAPPINGS,
 )
@@ -56,38 +61,35 @@ def dump_four_c_beam_to_list(self) -> dict:
     # Check the material.
     self._check_material()
 
-    # Gather the element data for the input file.
-    element_data = type(self).four_c_element_data | self.data
-    element_data["MAT"] = self.material
-    if type(self).four_c_triads:
-        element_data["TRIADS"] = [
+    # Get the legacy dictionary for FourCIPP.
+    node_ordering = _INPUT_FILE_MAPPINGS["beam_n_nodes_to_four_c_ordering"][
+        len(self.nodes)
+    ]
+    legacy_dict = _get_four_c_legacy_dict(self, node_ordering=node_ordering)
+    dump_triads = {
+        _BeamType.reissner: True,
+        _BeamType.kirchhoff: True,
+        _BeamType.euler_bernoulli: False,
+    }
+    if dump_triads[type(self).beam_type]:
+        legacy_dict["data"]["TRIADS"] = [
             item
-            for i in _INPUT_FILE_MAPPINGS["n_nodes_to_node_ordering"][len(self.nodes)]
+            for i in node_ordering
             for item in self.nodes[i].rotation.get_rotation_vector()
         ]
 
-    return {
-        "id": self.i_global + 1,
-        "cell": {
-            "type": _INPUT_FILE_MAPPINGS["n_nodes_to_cell_type"][len(self.nodes)],
-            "connectivity": [
-                self.nodes[i]
-                for i in _INPUT_FILE_MAPPINGS["n_nodes_to_node_ordering"][
-                    len(self.nodes)
-                ]
-            ],
-        },
-        "data": element_data,
-    }
+    return legacy_dict
 
 
 def get_four_c_reissner_beam(n_nodes: int, is_hermite_centerline: bool) -> type[_Beam]:
     """Return a Simo-Reissner beam for 4C."""
 
-    four_c_element_data = {
-        "type": _INPUT_FILE_MAPPINGS["beam_types"][_BeamType.reissner],
-        "HERMITE_CENTERLINE": is_hermite_centerline,
-    }
+    four_c_type = _INPUT_FILE_MAPPINGS["beam_type_to_four_c_type"][_BeamType.reissner]
+    four_c_cell = _INPUT_FILE_MAPPINGS["element_type_and_n_nodes_to_four_c_cell"][
+        _bme.element_type.beam, n_nodes
+    ]
+    element_technology = {"HERMITE_CENTERLINE": is_hermite_centerline}
+
     if is_hermite_centerline:
         coupling_fix_dict = {"NUMDOF": 9, "ONOFF": [1, 1, 1, 1, 1, 1, 0, 0, 0]}
         coupling_joint_dict = {"NUMDOF": 9, "ONOFF": [1, 1, 1, 0, 0, 0, 0, 0, 0]}
@@ -99,9 +101,9 @@ def get_four_c_reissner_beam(n_nodes: int, is_hermite_centerline: bool) -> type[
         "BeamFourCSimoReissner",
         (_generate_beam_class(n_nodes),),
         {
-            "four_c_beam_type": _BeamType.reissner,
-            "four_c_triads": True,
-            "four_c_element_data": four_c_element_data,
+            "element_type": _bme.element_type.beam,
+            "beam_type": _BeamType.reissner,
+            "data": {four_c_type: {four_c_cell: element_technology}},
             "valid_materials": [_MaterialReissner, _MaterialReissnerElastoplastic],
             "coupling_fix_dict": coupling_fix_dict,
             "coupling_joint_dict": coupling_joint_dict,
@@ -117,14 +119,6 @@ def get_four_c_kirchhoff_beam(
 ) -> type[_Beam]:
     """Return a Kirchhoff-Love beam for 4C."""
 
-    # Set the parameters for this beam.
-    four_c_element_data = {
-        "type": _INPUT_FILE_MAPPINGS["beam_types"][_BeamType.kirchhoff],
-        "CONSTRAINT": constraint.name,
-        "PARAMETRIZATION": parametrization.name,
-        "USE_FAD": is_fad,
-    }
-
     # Show warning when not using rotvec.
     if not parametrization == _BeamKirchhoffParametrizationType.rot:
         _warnings.warn(
@@ -132,16 +126,29 @@ def get_four_c_kirchhoff_beam(
             " applying the boundary conditions and couplings."
         )
 
+    n_nodes = 3
+    four_c_type = _INPUT_FILE_MAPPINGS["beam_type_to_four_c_type"][_BeamType.kirchhoff]
+    four_c_cell = _INPUT_FILE_MAPPINGS["element_type_and_n_nodes_to_four_c_cell"][
+        _bme.element_type.beam, n_nodes
+    ]
+
+    element_technology = {
+        "CONSTRAINT": constraint.name,
+        "PARAMETRIZATION": parametrization.name,
+        "USE_FAD": is_fad,
+    }
+
     coupling_fix_dict = {"NUMDOF": 7, "ONOFF": [1, 1, 1, 1, 1, 1, 0]}
     coupling_joint_dict = {"NUMDOF": 7, "ONOFF": [1, 1, 1, 0, 0, 0, 0]}
 
     return type(
         "BeamFourCKirchhoffLove",
-        (_generate_beam_class(3),),
+        (_generate_beam_class(n_nodes),),
         {
-            "four_c_beam_type": _BeamType.kirchhoff,
-            "four_c_triads": True,
-            "four_c_element_data": four_c_element_data,
+            "element_type": _bme.element_type.beam,
+            "beam_type": _BeamType.kirchhoff,
+            "kirchhoff_parametrization": parametrization,
+            "data": {four_c_type: {four_c_cell: element_technology}},
             "valid_materials": [_MaterialKirchhoff],
             "coupling_fix_dict": coupling_fix_dict,
             "coupling_joint_dict": coupling_joint_dict,
@@ -150,13 +157,17 @@ def get_four_c_kirchhoff_beam(
     )
 
 
-class BeamFourCEulerBernoulli(_generate_beam_class(2)):  # type: ignore[misc]
+class BeamFourCEulerBernoulli(_Beam2):
     """Represents a Euler Bernoulli beam element."""
 
-    four_c_beam_type = _BeamType.euler_bernoulli
-    four_c_triads = False
-    four_c_element_data = {
-        "type": _INPUT_FILE_MAPPINGS["beam_types"][_BeamType.euler_bernoulli]
+    element_type = _bme.element_type.beam
+    beam_type = _BeamType.euler_bernoulli
+    data: dict[str, dict[str, _Any]] = {
+        _INPUT_FILE_MAPPINGS["beam_type_to_four_c_type"][_BeamType.euler_bernoulli]: {
+            _INPUT_FILE_MAPPINGS["element_type_and_n_nodes_to_four_c_cell"][
+                _bme.element_type.beam, len(_Beam2.nodes_create)
+            ]: {}
+        }
     }
 
     valid_materials = [_MaterialEulerBernoulli]
@@ -171,8 +182,8 @@ def dump_to_list(self):
         # the start point to the end point.
         if not self.nodes[0].rotation == self.nodes[1].rotation:
             raise ValueError(
-                "The two nodal rotations in Euler Bernoulli beams must be the same, i.e. the beam "
-                "has to be straight!"
+                "The two nodal rotations in Euler Bernoulli beams must be the same,"
+                "i.e., the beam has to be straight!"
             )
         direction = self.nodes[1].coordinates - self.nodes[0].coordinates
         t1 = self.nodes[0].rotation * [1, 0, 0]
diff --git a/src/beamme/four_c/element_solid.py b/src/beamme/four_c/element_solid.py
new file mode 100644
index 00000000..3a68ddd5
--- /dev/null
+++ b/src/beamme/four_c/element_solid.py
@@ -0,0 +1,120 @@
+# The MIT License (MIT)
+#
+# Copyright (c) 2018-2026 BeamMe Authors
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+"""This file implements solid elements for 4C."""
+
+import copy as _copy
+
+from beamme.core.conf import bme as _bme
+from beamme.core.element import Element as _Element
+from beamme.core.element_volume import VolumeElement as _VolumeElement
+from beamme.core.element_volume import VolumeHEX8 as _VolumeHEX8
+from beamme.core.element_volume import VolumeHEX20 as _VolumeHEX20
+from beamme.core.element_volume import VolumeHEX27 as _VolumeHEX27
+from beamme.core.element_volume import VolumeTET4 as _VolumeTET4
+from beamme.core.element_volume import VolumeTET10 as _VolumeTET10
+from beamme.core.element_volume import VolumeWEDGE6 as _VolumeWEDGE6
+from beamme.core.nurbs_patch import NURBSSurface as _NURBSSurface
+from beamme.core.nurbs_patch import NURBSVolume as _NURBSVolume
+from beamme.four_c.input_file_mappings import (
+    INPUT_FILE_MAPPINGS as _INPUT_FILE_MAPPINGS,
+)
+
+
+def get_four_c_solid(
+    solid_type_string: str,
+    n_nodes: int | None,
+    *,
+    element_technology: dict | None = None,
+) -> type[_Element]:
+    """Return a type that defines a solid element block for 4C solid elements.
+
+    Args:
+        solid_type_string: The string defining the type of solid element. Supported
+            values are:
+            - "nurbs_2d"
+            - "nurbs_3d"
+            - "nurbs_shell"
+            - "solid"
+            - "rigid_sphere"
+        n_nodes: The number of nodes of the solid element.
+        element_technology: Dictionary specifying element technologies. Will be deep copied.
+
+    Returns:
+        A type that defines a solid element block for 4C solid elements.
+    """
+
+    if element_technology is None:
+        element_technology = {}
+    else:
+        element_technology = _copy.deepcopy(element_technology)
+
+    base_type: type[_Element]
+
+    if solid_type_string.startswith("nurbs"):
+        if solid_type_string.endswith("2d"):
+            base_type = _NURBSSurface
+        elif solid_type_string.endswith("3d"):
+            base_type = _NURBSVolume
+        elif solid_type_string.endswith("shell"):
+            base_type = _NURBSSurface
+        else:
+            raise ValueError(f"Unsupported NURBS type {solid_type_string}!")
+        element_type = _bme.element_type.nurbs
+    elif solid_type_string == "rigid_sphere":
+        base_type = _VolumeElement
+        element_type = _bme.element_type.solid
+    elif solid_type_string == "solid":
+        match n_nodes:
+            case 8:
+                base_type = _VolumeHEX8
+            case 20:
+                base_type = _VolumeHEX20
+            case 27:
+                base_type = _VolumeHEX27
+            case 4:
+                base_type = _VolumeTET4
+            case 10:
+                base_type = _VolumeTET10
+            case 6:
+                base_type = _VolumeWEDGE6
+            case _:
+                raise ValueError(
+                    f"Unsupported number of nodes {n_nodes} for solid element!"
+                )
+        element_type = _bme.element_type.solid
+    else:
+        raise ValueError(f"Unsupported solid type {solid_type_string}!")
+
+    four_c_cell = _INPUT_FILE_MAPPINGS["element_type_and_n_nodes_to_four_c_cell"][
+        element_type, n_nodes
+    ]
+    four_c_type = _INPUT_FILE_MAPPINGS["solid_type_to_four_c_type"][solid_type_string]
+
+    solid_type = type(
+        "FourCSolidElementType",
+        (base_type,),
+        {
+            "element_type": element_type,
+            "data": {four_c_type: {four_c_cell: element_technology}},
+        },
+    )
+    return solid_type
diff --git a/src/beamme/four_c/element_volume.py b/src/beamme/four_c/element_volume.py
deleted file mode 100644
index c01d8717..00000000
--- a/src/beamme/four_c/element_volume.py
+++ /dev/null
@@ -1,32 +0,0 @@
-# The MIT License (MIT)
-#
-# Copyright (c) 2018-2026 BeamMe Authors
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in
-# all copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-# THE SOFTWARE.
-"""This file implements volume elements for 4C."""
-
-from beamme.core.element_volume import VolumeElement as _VolumeElement
-
-
-class SolidRigidSphere(_VolumeElement):
-    """A rigid sphere solid element."""
-
-    def __init__(self, **kwargs):
-        """Initialize solid sphere object."""
-        _VolumeElement.__init__(self, **kwargs)
diff --git a/src/beamme/four_c/input_file.py b/src/beamme/four_c/input_file.py
index 3c3f2658..5550a57b 100644
--- a/src/beamme/four_c/input_file.py
+++ b/src/beamme/four_c/input_file.py
@@ -316,12 +316,14 @@ def add_mesh_to_input_file(self, mesh: _Mesh) -> None:
         nurbs_count = 0
 
         for element in mesh.elements:
-            element.i_global = i
             if isinstance(element, _NURBSPatch):
+                element.i_global_start = i
                 element.i_nurbs_patch = nurbs_count
                 i += element.get_number_of_elements()
                 nurbs_count += 1
                 continue
+            else:
+                element.i_global = i
             i += 1
 
         #   Materials: Get a list of all materials in the mesh,
diff --git a/src/beamme/four_c/input_file_dump_item.py b/src/beamme/four_c/input_file_dump_item.py
index 4976bd60..443c22bc 100644
--- a/src/beamme/four_c/input_file_dump_item.py
+++ b/src/beamme/four_c/input_file_dump_item.py
@@ -26,9 +26,9 @@
 from beamme.core.boundary_condition import BoundaryCondition as _BoundaryCondition
 from beamme.core.conf import bme as _bme
 from beamme.core.coupling import Coupling as _Coupling
+from beamme.core.element import Element as _Element
 from beamme.core.element_volume import VolumeElement as _VolumeElement
-from beamme.core.geometry_set import GeometrySet as _GeometrySet
-from beamme.core.geometry_set import GeometrySetNodes as _GeometrySetNodes
+from beamme.core.geometry_set import GeometrySetBase as _GeometrySetBase
 from beamme.core.node import ControlPoint as _ControlPoint
 from beamme.core.node import Node as _Node
 from beamme.core.nurbs_patch import NURBSPatch as _NURBSPatch
@@ -41,6 +41,51 @@
 )
 
 
+def get_four_c_legacy_dict(
+    element: _Element, *, node_ordering: list[int] | None = None
+) -> dict[str, _Any]:
+    """Get the 4C legacy data dict for the given element.
+
+    Args:
+        element: The element for which the legacy data dict should be returned.
+        node_ordering: The ordering of nodes for the element in 4C.
+
+    Returns:
+        A dict with the 4C legacy data for the given element (including material entry).
+    """
+
+    element_type_data = getattr(type(element), "data", {})
+    if len(element_type_data) == 1:
+        four_c_type, four_c_cell_dict = next(iter(element_type_data.items()))
+    else:
+        raise ValueError("Expected exactly one entry in type dictionary")
+    if len(four_c_cell_dict) == 1:
+        four_c_cell, four_c_data = next(iter(four_c_cell_dict.items()))
+    else:
+        raise ValueError("Expected exactly one entry in cell dictionary")
+
+    if node_ordering is None:
+        connectivity = element.nodes
+    else:
+        connectivity = [element.nodes[i] for i in node_ordering]
+
+    if element.material is None:
+        material_dict = {}
+    else:
+        material_dict = {"MAT": element.material}
+
+    if element.i_global is not None:
+        four_c_id = element.i_global + 1
+    else:
+        four_c_id = None
+
+    return {
+        "id": four_c_id,
+        "cell": {"type": four_c_cell, "connectivity": connectivity},
+        "data": {"type": four_c_type, **material_dict, **four_c_data},
+    }
+
+
 def dump_node(node):
     """Return the representation of a node in the 4C input file."""
 
@@ -62,29 +107,7 @@ def dump_node(node):
 
 def dump_solid_element(solid_element):
     """Return a dict with the items representing the given solid element."""
-
-    if "MAT" in solid_element.data:
-        raise ValueError(
-            f"Element {solid_element.i_global} has a MAT entry in its data, this "
-            "is not supported, the materials have to be assigned via the material "
-            "attribute of the element."
-        )
-
-    return {
-        "id": solid_element.i_global + 1,
-        "cell": {
-            "type": _INPUT_FILE_MAPPINGS["element_type_to_four_c_string"][
-                type(solid_element)
-            ],
-            "connectivity": solid_element.nodes,
-        },
-        "data": solid_element.data
-        | (
-            {"MAT": solid_element.material}
-            if solid_element.material is not None
-            else {}
-        ),
-    }
+    return get_four_c_legacy_dict(solid_element)
 
 
 def dump_coupling(coupling):
@@ -106,11 +129,9 @@ def dump_coupling(coupling):
                 f"Expected a single connected type of beam elements, got {element_types}"
             )
         element_type = element_types.pop()
-        if element_type.four_c_beam_type is _BeamType.kirchhoff:
+        if element_type.beam_type is _BeamType.kirchhoff:
             unique_parametrization_flags = {
-                _BeamKirchhoffParametrizationType[
-                    type(element).four_c_element_data["PARAMETRIZATION"]
-                ]
+                type(element).kirchhoff_parametrization
                 for element in connected_elements
             }
             if (
@@ -205,10 +226,12 @@ def dump_nurbs_patch_elements(nurbs_patch: _NURBSPatch) -> list[dict[str, _Any]]
     """Return a list with all the element definitions contained in this
     patch."""
 
-    if nurbs_patch.i_global is None:
+    if nurbs_patch.i_global is not None:
         raise ValueError(
-            "i_global is not set, make sure that the NURBS patch is added to the mesh"
+            f"For NURBS patches, i_global shall not be set, got {nurbs_patch.i_global}."
         )
+    if nurbs_patch.i_nurbs_patch is None:
+        raise ValueError("Expected i_nurbs_patch to be set for NURBS patch.")
 
     # Check the material
     nurbs_patch._check_material()
@@ -218,25 +241,10 @@ def dump_nurbs_patch_elements(nurbs_patch: _NURBSPatch) -> list[dict[str, _Any]]
 
     for knot_span in nurbs_patch.get_knot_span_iterator():
         element_cps_ids = nurbs_patch.get_ids_ctrlpts(*knot_span)
-        connectivity = [nurbs_patch.nodes[i] for i in element_cps_ids]
-        num_cp = len(connectivity)
-
-        patch_elements.append(
-            {
-                "id": nurbs_patch.i_global + j + 1,
-                "cell": {
-                    "type": f"NURBS{num_cp}",
-                    "connectivity": connectivity,
-                },
-                "data": {
-                    "type": _INPUT_FILE_MAPPINGS["nurbs_type_to_default_four_c_type"][
-                        type(nurbs_patch)
-                    ],
-                    "MAT": nurbs_patch.material,
-                    **(nurbs_patch.data if nurbs_patch.data else {}),
-                },
-            }
-        )
+        legacy_dict = get_four_c_legacy_dict(nurbs_patch, node_ordering=element_cps_ids)
+        # We need to overwrite the element ID here, as one patch contains multiple elements.
+        legacy_dict["id"] = nurbs_patch.i_global_start + j + 1
+        patch_elements.append(legacy_dict)
         j += 1
 
     return patch_elements
@@ -250,7 +258,7 @@ def dump_item_to_list(dumped_list, item) -> None:
         dumped_list.append(dump_node(item))
     elif isinstance(item, _VolumeElement):
         dumped_list.append(dump_solid_element(item))
-    elif isinstance(item, _GeometrySet) or isinstance(item, _GeometrySetNodes):
+    elif isinstance(item, _GeometrySetBase):
         dumped_list.extend(dump_geometry_set(item))
     elif isinstance(item, _NURBSPatch):
         dumped_list.extend(dump_nurbs_patch_elements(item))
diff --git a/src/beamme/four_c/input_file_mappings.py b/src/beamme/four_c/input_file_mappings.py
index 86891b83..49472a65 100644
--- a/src/beamme/four_c/input_file_mappings.py
+++ b/src/beamme/four_c/input_file_mappings.py
@@ -25,35 +25,55 @@
 from typing import Any as _Any
 
 from beamme.core.conf import bme as _bme
-from beamme.core.element_volume import (
-    VolumeHEX8 as _VolumeHEX8,
-)
-from beamme.core.element_volume import (
-    VolumeHEX20 as _VolumeHEX20,
-)
-from beamme.core.element_volume import (
-    VolumeHEX27 as _VolumeHEX27,
-)
-from beamme.core.element_volume import (
-    VolumeTET4 as _VolumeTET4,
-)
-from beamme.core.element_volume import (
-    VolumeTET10 as _VolumeTET10,
-)
-from beamme.core.element_volume import (
-    VolumeWEDGE6 as _VolumeWEDGE6,
-)
-from beamme.core.nurbs_patch import NURBSSurface as _NURBSSurface
-from beamme.core.nurbs_patch import NURBSVolume as _NURBSVolume
-from beamme.four_c.element_volume import SolidRigidSphere as _SolidRigidSphere
 from beamme.four_c.four_c_types import BeamType as _BeamType
+from beamme.utils.data_structures import (
+    create_inverse_mapping as _create_inverse_mapping,
+)
 
 INPUT_FILE_MAPPINGS: dict[str, _Any] = {}
-INPUT_FILE_MAPPINGS["beam_types"] = {
+INPUT_FILE_MAPPINGS["beam_type_to_four_c_type"] = {
     _BeamType.reissner: "BEAM3R",
     _BeamType.kirchhoff: "BEAM3K",
     _BeamType.euler_bernoulli: "BEAM3EB",
 }
+INPUT_FILE_MAPPINGS["solid_type_to_four_c_type"] = {
+    "nurbs_2d": "WALLNURBS",
+    "nurbs_3d": "SOLID",
+    "nurbs_shell": "SHELL_KIRCHHOFF_LOVE_NURBS",
+    "solid": "SOLID",
+    "rigid_sphere": "RIGIDSPHERE",
+}
+INPUT_FILE_MAPPINGS["four_c_type_to_solid_type"] = {
+    "SOLID": "solid",
+    "RIGIDSPHERE": "rigid_sphere",
+}
+INPUT_FILE_MAPPINGS["element_type_and_n_nodes_to_four_c_cell"] = {
+    (_bme.element_type.beam, 2): "LINE2",
+    (_bme.element_type.beam, 3): "LINE3",
+    (_bme.element_type.beam, 4): "LINE4",
+    (_bme.element_type.beam, 5): "LINE5",
+    (_bme.element_type.nurbs, 9): "NURBS9",
+    (_bme.element_type.nurbs, 27): "NURBS27",
+    (_bme.element_type.solid, 8): "HEX8",
+    (_bme.element_type.solid, 20): "HEX20",
+    (_bme.element_type.solid, 27): "HEX27",
+    (_bme.element_type.solid, 4): "TET4",
+    (_bme.element_type.solid, 10): "TET10",
+    (_bme.element_type.solid, 6): "WEDGE6",
+    (_bme.element_type.solid, 1): "POINT1",
+}
+INPUT_FILE_MAPPINGS["geometry_sets_geometry_to_entry_name"] = {
+    _bme.geo.point: "DNODE",
+    _bme.geo.line: "DLINE",
+    _bme.geo.surface: "DSURFACE",
+    _bme.geo.volume: "DVOL",
+}
+INPUT_FILE_MAPPINGS["beam_n_nodes_to_four_c_ordering"] = {
+    2: [0, 1],
+    3: [0, 2, 1],
+    4: [0, 3, 1, 2],
+    5: [0, 4, 1, 2, 3],
+}
 INPUT_FILE_MAPPINGS["boundary_conditions"] = {
     (_bme.bc.dirichlet, _bme.geo.point): "DESIGN POINT DIRICH CONDITIONS",
     (_bme.bc.dirichlet, _bme.geo.line): "DESIGN LINE DIRICH CONDITIONS",
@@ -113,50 +133,14 @@
         _bme.geo.surface,
     ): "DESIGN SURF MORTAR CONTACT CONDITIONS 3D",
 }
-INPUT_FILE_MAPPINGS["element_type_to_four_c_string"] = {
-    _VolumeHEX8: "HEX8",
-    _VolumeHEX20: "HEX20",
-    _VolumeHEX27: "HEX27",
-    _VolumeTET4: "TET4",
-    _VolumeTET10: "TET10",
-    _VolumeWEDGE6: "WEDGE6",
-    _SolidRigidSphere: "POINT1",
-}
-INPUT_FILE_MAPPINGS["element_four_c_string_to_type"] = {
-    value: key
-    for key, value in INPUT_FILE_MAPPINGS["element_type_to_four_c_string"].items()
-}
 INPUT_FILE_MAPPINGS["geometry_sets_geometry_to_condition_name"] = {
     _bme.geo.point: "DNODE-NODE TOPOLOGY",
     _bme.geo.line: "DLINE-NODE TOPOLOGY",
     _bme.geo.surface: "DSURF-NODE TOPOLOGY",
     _bme.geo.volume: "DVOL-NODE TOPOLOGY",
 }
-INPUT_FILE_MAPPINGS["geometry_sets_condition_to_geometry_name"] = {
-    value: key
-    for key, value in INPUT_FILE_MAPPINGS[
-        "geometry_sets_geometry_to_condition_name"
-    ].items()
-}
-INPUT_FILE_MAPPINGS["geometry_sets_geometry_to_entry_name"] = {
-    _bme.geo.point: "DNODE",
-    _bme.geo.line: "DLINE",
-    _bme.geo.surface: "DSURFACE",
-    _bme.geo.volume: "DVOL",
-}
-INPUT_FILE_MAPPINGS["n_nodes_to_cell_type"] = {
-    2: "LINE2",
-    3: "LINE3",
-    4: "LINE4",
-    5: "LINE5",
-}
-INPUT_FILE_MAPPINGS["n_nodes_to_node_ordering"] = {
-    2: [0, 1],
-    3: [0, 2, 1],
-    4: [0, 3, 1, 2],
-    5: [0, 4, 1, 2, 3],
-}
-INPUT_FILE_MAPPINGS["nurbs_type_to_default_four_c_type"] = {
-    _NURBSSurface: "WALLNURBS",
-    _NURBSVolume: "SOLID",
-}
+INPUT_FILE_MAPPINGS["geometry_sets_condition_to_geometry_name"] = (
+    _create_inverse_mapping(
+        INPUT_FILE_MAPPINGS["geometry_sets_geometry_to_condition_name"]
+    )
+)
diff --git a/src/beamme/four_c/model_importer.py b/src/beamme/four_c/model_importer.py
index 144277e8..fb33bf4a 100644
--- a/src/beamme/four_c/model_importer.py
+++ b/src/beamme/four_c/model_importer.py
@@ -34,11 +34,15 @@
 from beamme.core.geometry_set import GeometrySetNodes as _GeometrySetNodes
 from beamme.core.mesh import Mesh as _Mesh
 from beamme.core.node import Node as _Node
+from beamme.four_c.element_solid import get_four_c_solid as _get_four_c_solid
 from beamme.four_c.input_file import InputFile as _InputFile
 from beamme.four_c.input_file_mappings import (
     INPUT_FILE_MAPPINGS as _INPUT_FILE_MAPPINGS,
 )
 from beamme.four_c.material import MaterialSolid as _MaterialSolid
+from beamme.utils.data_structures import (
+    compare_nested_dicts_or_lists as _compare_nested_dicts_or_lists,
+)
 from beamme.utils.environment import cubitpy_is_available as _cubitpy_is_available
 
 if _cubitpy_is_available():
@@ -159,21 +163,35 @@ def _extract_mesh_sections(input_file: _InputFile) -> _Tuple[_InputFile, _Mesh]:
     mesh.nodes = [_Node(node["COORD"]) for node in _pop_section("NODE COORDS")]
 
     # extract elements
+    element_types: dict[type, dict] = {}
     for input_element in _pop_section("STRUCTURE ELEMENTS"):
-        if (
-            input_element["cell"]["type"]
-            not in _INPUT_FILE_MAPPINGS["element_four_c_string_to_type"]
-        ):
-            raise TypeError(
-                f"Could not create a BeamMe element for `{input_element['data']['type']}` `{input_element['cell']['type']}`!"
+        # At this point `input_element` contains the full element data. We can not
+        # compare this directly with the other block data dictionaries, we first
+        # need to extract connectivity, cell ID and material information from
+        # `input_element`.
+        input_element.pop("id")
+        connectivity = input_element["cell"].pop("connectivity")
+        four_c_type = input_element["data"].pop("type")
+        material_id = input_element["data"].pop("MAT", None)
+
+        # Loop over already found block element types and check if the current
+        # element matches one of those. If not, create a new block element type
+        # for this element.
+        for element_type, type_data in element_types.items():
+            if _compare_nested_dicts_or_lists(type_data, input_element):
+                break
+        else:
+            element_type = _get_four_c_solid(
+                _INPUT_FILE_MAPPINGS["four_c_type_to_solid_type"][four_c_type],
+                n_nodes=len(connectivity),
+                element_technology=input_element["data"],
             )
-        nodes = [mesh.nodes[i - 1] for i in input_element["cell"]["connectivity"]]
-        element_class = _INPUT_FILE_MAPPINGS["element_four_c_string_to_type"][
-            input_element["cell"]["type"]
-        ]
-        element = element_class(nodes=nodes, data=input_element["data"])
-        if "MAT" in element.data:
-            element.material = material_id_map[element.data.pop("MAT")]
+            element_types[element_type] = input_element
+
+        nodes = [mesh.nodes[i - 1] for i in connectivity]
+        element = element_type(nodes=nodes)
+        if material_id is not None:
+            element.material = material_id_map[material_id]
         mesh.elements.append(element)
 
     # extract geometry sets
diff --git a/src/beamme/four_c/solid_shell_thickness_direction.py b/src/beamme/four_c/solid_shell_thickness_direction.py
index 3c89ac6a..6b093b91 100644
--- a/src/beamme/four_c/solid_shell_thickness_direction.py
+++ b/src/beamme/four_c/solid_shell_thickness_direction.py
@@ -197,11 +197,14 @@ def set_solid_shell_thickness_direction(
         raise ValueError("Expected a non empty element list")
 
     for element in elements:
-        if (
-            isinstance(element, _VolumeHEX8)
-            and isinstance(element.data, dict)
-            and element.data.get("TECH") == "shell_eas_ans"
-        ):
+        is_hex8_solid_shell = (
+            getattr(type(element), "data", {})
+            .get("SOLID", {})
+            .get("HEX8", {})
+            .get("TECH", "")
+            == "shell_eas_ans"
+        )
+        if is_hex8_solid_shell:
             # Get the element center and the Jacobian at the center
             (
                 reference_position_center,
diff --git a/src/beamme/mesh_creation_functions/nurbs_generic.py b/src/beamme/mesh_creation_functions/nurbs_generic.py
index bc54c4d4..5c40a108 100644
--- a/src/beamme/mesh_creation_functions/nurbs_generic.py
+++ b/src/beamme/mesh_creation_functions/nurbs_generic.py
@@ -22,7 +22,6 @@
 """Generic function used to create NURBS meshes."""
 
 import itertools as _itertools
-from typing import Type as _Type
 
 import numpy as _np
 
@@ -35,20 +34,38 @@
 from beamme.core.nurbs_patch import NURBSVolume as _NURBSVolume
 
 
+def _check_nurbs_dimension_and_element_type(
+    nurbs_dimension: int, element_type: type
+) -> None:
+    """Check if the element type is compatible with the NURBS dimension.
+
+    Args:
+        nurbs_dimension: The dimension of the NURBS patch (2 for surface, 3 for volume).
+        element_type: The type of element to be created.
+
+    Raises:
+        ValueError: If the element type is not compatible with the NURBS dimension.
+    """
+    if nurbs_dimension == 2 and not issubclass(element_type, _NURBSSurface):
+        raise ValueError(
+            "Error, expected element type to be a NURBSSurface for a NURBS surface!"
+        )
+    elif nurbs_dimension == 3 and not issubclass(element_type, _NURBSVolume):
+        raise ValueError(
+            "Error, expected element type to be a NURBSVolume for a NURBS volume!"
+        )
+
+
 def add_splinepy_nurbs_to_mesh(
-    mesh: _Mesh,
-    splinepy_obj,
-    *,
-    material=None,
-    data: dict | None = None,
+    mesh: _Mesh, element_type: type, splinepy_obj, *, material=None
 ) -> _GeometryName:
     """Add a splinepy NURBS to the mesh.
 
     Args:
         mesh: Mesh that the created NURBS geometry will be added to.
+        element_type: The type of element to be created.
         splinepy_obj (splinepy object): NURBS geometry created using splinepy.
         material (Material): Material for this geometry.
-        data: General element data, e.g., material, formulation, ...
 
     Returns:
         GeometryName:
@@ -85,24 +102,16 @@ def add_splinepy_nurbs_to_mesh(
         )
     ]
 
-    # Fill element
-    manifold_dim = len(splinepy_obj.knot_vectors)
-    nurbs_object: _Type[_NURBSSurface] | _Type[_NURBSVolume]
-    if manifold_dim == 2:
-        nurbs_object = _NURBSSurface
-    elif manifold_dim == 3:
-        nurbs_object = _NURBSVolume
-    else:
-        raise NotImplementedError(
-            "Error, not implemented for a NURBS {}!".format(type(splinepy_obj))
-        )
+    # Create elements
+    _check_nurbs_dimension_and_element_type(
+        len(splinepy_obj.knot_vectors), element_type
+    )
 
-    element = nurbs_object(
+    element = element_type(
         [_np.asarray(knot_vector) for knot_vector in splinepy_obj.knot_vectors],
         _np.asarray(splinepy_obj.degrees),
         nodes=control_points,
         material=material,
-        data=data,
     )
 
     # Add element and control points to the mesh
@@ -116,19 +125,15 @@ def add_splinepy_nurbs_to_mesh(
 
 
 def add_geomdl_nurbs_to_mesh(
-    mesh: _Mesh,
-    geomdl_obj,
-    *,
-    material=None,
-    data: dict | None = None,
+    mesh: _Mesh, element_type: type, geomdl_obj, *, material=None
 ) -> _GeometryName:
     """Add a geomdl NURBS to the mesh.
 
     Args:
         mesh: Mesh that the created NURBS geometry will be added to.
+        element_type: The type of element to be created.
         geomdl_obj (geomdl object): NURBS geometry created using geomdl.
         material (Material): Material for this geometry.
-        data: General element data, e.g., material, formulation, ...
 
     Returns:
         GeometryName:
@@ -166,24 +171,13 @@ def add_geomdl_nurbs_to_mesh(
             )
         )
 
-    # Fill element
-    manifold_dim = len(geomdl_obj.knotvector)
-    nurbs_object: _Type[_NURBSSurface] | _Type[_NURBSVolume]
-    if manifold_dim == 2:
-        nurbs_object = _NURBSSurface
-    elif manifold_dim == 3:
-        nurbs_object = _NURBSVolume
-    else:
-        raise NotImplementedError(
-            "Error, not implemented for a NURBS {}!".format(type(geomdl_obj))
-        )
-
-    element = nurbs_object(
+    # Create elements
+    _check_nurbs_dimension_and_element_type(len(geomdl_obj.knotvector), element_type)
+    element = element_type(
         geomdl_obj.knotvector,
         geomdl_obj.degree,
         nodes=control_points,
         material=material,
-        data=data,
     )
 
     # Add element and control points to the mesh
diff --git a/src/beamme/utils/data_structures.py b/src/beamme/utils/data_structures.py
new file mode 100644
index 00000000..2e04c4fb
--- /dev/null
+++ b/src/beamme/utils/data_structures.py
@@ -0,0 +1,101 @@
+# The MIT License (MIT)
+#
+# Copyright (c) 2018-2026 BeamMe Authors
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+"""Helper functions for data structure related functionality."""
+
+from typing import Any as _Any
+from typing import Callable as _Callable
+
+from fourcipp.utils.dict_utils import (
+    compare_nested_dicts_or_lists as _compare_nested_dicts_or_lists,
+)
+
+
+def create_inverse_mapping(mapping: dict[_Any, _Any]) -> dict[_Any, _Any]:
+    """Create an inverse mapping from a given mapping.
+
+    Args:
+        mapping: The mapping to create the inverse mapping from.
+
+    Returns:
+        The inverse mapping.
+    """
+    # Check that the mapping is invertible, i.e., that all values are unique.
+    unique_values = set(mapping.values())
+    if len(unique_values) != len(mapping):
+        unique_values = set()
+        duplicate_values = set()
+        for value in mapping.values():
+            if value in unique_values:
+                duplicate_values.add(value)
+            else:
+                unique_values.add(value)
+        raise ValueError(
+            "The mapping is not invertible, values are not unique. "
+            f"Non-unique values: {duplicate_values}"
+        )
+    return {value: key for key, value in mapping.items()}
+
+
+def compare_nested_dicts_or_lists(
+    obj: _Any,
+    reference_obj: _Any,
+    allow_int_vs_float_comparison: bool = False,
+    rtol: float = 1.0e-5,
+    atol: float = 1.0e-8,
+    equal_nan: bool = False,
+    custom_compare: _Callable | None = None,
+) -> bool:
+    """Recursively compare two nested dictionaries or lists.
+
+    This function is taken from FourCIPP and modified such that no assertion error
+    is raised if the objects are not equal but instead a boolean is returned.
+
+    To compare custom python objects, a `custom_compare` callable can be provided which:
+        - Returns nothing/`None` if the objects where not compared within `custom_compare`
+        - Returns `True` if the objects are seen as equal
+        - Raises AssertionError if the objects are not equal
+
+    Args:
+        obj: Object for comparison
+        reference_obj: Reference object
+        allow_int_vs_float_comparison: Allow a tolerance based comparison between int and
+            float
+        rtol: The relative tolerance parameter for numpy.isclose
+        atol: The absolute tolerance parameter for numpy.isclose
+        equal_nan: Whether to compare NaN's as equal for numpy.isclose
+        custom_compare: Callable to compare objects within this nested framework
+
+    Returns:
+        True if the dictionaries are equal
+    """
+    try:
+        return _compare_nested_dicts_or_lists(
+            obj,
+            reference_obj,
+            allow_int_vs_float_comparison,
+            rtol,
+            atol,
+            equal_nan,
+            custom_compare,
+        )
+    except AssertionError:
+        return False
diff --git a/tests/beamme/mesh_creation_functions/test_beamme_mesh_creation_functions_nurbs_generic.py b/tests/beamme/mesh_creation_functions/test_beamme_mesh_creation_functions_nurbs_generic.py
index 2144504b..64393de7 100644
--- a/tests/beamme/mesh_creation_functions/test_beamme_mesh_creation_functions_nurbs_generic.py
+++ b/tests/beamme/mesh_creation_functions/test_beamme_mesh_creation_functions_nurbs_generic.py
@@ -24,6 +24,7 @@
 import pytest
 
 from beamme.core.mesh import Mesh
+from beamme.four_c.element_solid import get_four_c_solid
 from beamme.mesh_creation_functions.nurbs_generic import (
     add_geomdl_nurbs_to_mesh,
     create_geometry_sets,
@@ -35,10 +36,12 @@
 
 
 @pytest.mark.parametrize(
-    ("nurbs_patch", "reference_values"),
+    ("nurbs_patch", "solid_type_string", "n_nodes", "reference_values"),
     [
         (
             create_nurbs_flat_plate_2d(1, 2, n_ele_u=1, n_ele_v=2),
+            "nurbs_2d",
+            9,
             {
                 "vertex_u_min_v_min": [0],
                 "vertex_u_min_v_max": [9],
@@ -57,6 +60,8 @@
         ),
         (
             create_nurbs_brick(1, 2, 3, n_ele_u=1, n_ele_v=2, n_ele_w=3),
+            "nurbs_3d",
+            27,
             {
                 "vertex_u_min_v_min_w_min": [0],
                 "vertex_u_min_v_min_w_max": [48],
@@ -209,14 +214,17 @@
     ],
 )
 def test_beamme_mesh_creation_functions_nurbs_generic_sets(
-    nurbs_patch, reference_values
+    nurbs_patch, solid_type_string, n_nodes, reference_values
 ):
     """Test that the add NURBS to mesh functionality returns the correct
     geometry sets."""
 
     # Add the nurbs to a mesh
     mesh = Mesh()
-    add_geomdl_nurbs_to_mesh(mesh, nurbs_patch)
+    element_type = get_four_c_solid(
+        solid_type_string=solid_type_string, n_nodes=n_nodes
+    )
+    add_geomdl_nurbs_to_mesh(mesh, element_type, nurbs_patch)
     nurbs_patch = mesh.elements[0]
 
     # Create the geometry sets for this patch
diff --git a/tests/beamme/four_c/test_beamme_four_c_input_file.py b/tests/beamme/utils/test_beamme_utils_data_structures.py
similarity index 64%
rename from tests/beamme/four_c/test_beamme_four_c_input_file.py
rename to tests/beamme/utils/test_beamme_utils_data_structures.py
index d1f28fd6..7fd28793 100644
--- a/tests/beamme/four_c/test_beamme_four_c_input_file.py
+++ b/tests/beamme/utils/test_beamme_utils_data_structures.py
@@ -19,25 +19,25 @@
 # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 # THE SOFTWARE.
-"""Unit tests for the 4C input file."""
+"""This script is used to unit test the data structure utility functions."""
 
 import pytest
 
-from beamme.core.element_volume import VolumeElement
-from beamme.four_c.input_file_dump_item import dump_solid_element
+from beamme.utils.data_structures import create_inverse_mapping
 
 
-def test_beamme_four_c_input_file_material_in_element_data():
-    """Check that material information can not be included in solid element
-    data."""
+def test_beamme_utils_data_structures_create_inverse_mapping():
+    """Test the create_inverse_mapping function."""
 
-    element = VolumeElement(
-        nodes=[0, 1, 2, 3],
-        data={"MAT": 1},
-    )
+    # Test with a simple mapping.
+    mapping = {1: "a", 2: "b", 3: "c"}
+    inverse_mapping = create_inverse_mapping(mapping)
+    assert inverse_mapping == {"a": 1, "b": 2, "c": 3}
+
+    # Test with a mapping that has duplicate values.
+    mapping = {1: "a", 2: "b", 3: "a"}
     with pytest.raises(
         ValueError,
-        match="Element None has a MAT entry in its data, this is not supported, "
-        "the materials have to be assigned via the material attribute of the element.",
+        match="The mapping is not invertible, values are not unique. Non-unique values: {'a'}",
     ):
-        dump_solid_element(element)
+        create_inverse_mapping(mapping)
diff --git a/tests/conftest_test_object_generators.py b/tests/conftest_test_object_generators.py
index 4e76c8fb..31253d71 100644
--- a/tests/conftest_test_object_generators.py
+++ b/tests/conftest_test_object_generators.py
@@ -28,8 +28,10 @@
 import pytest
 import splinepy
 
+from beamme.core.element import Element
 from beamme.core.material import MaterialBeamBase
 from beamme.cosserat_curve.cosserat_curve import CosseratCurve
+from beamme.four_c.element_solid import get_four_c_solid
 from beamme.four_c.material import (
     MaterialReissner,
     MaterialSolid,
@@ -167,45 +169,53 @@ def _get_default_test_solid_material(
 
 
 @pytest.fixture(scope="function")
-def get_default_test_solid_element_description() -> Callable:
-    """Return a function to create a default solid element description for
-    testing purposes.
+def get_default_test_solid_element() -> Callable:
+    """Return a function to create a default solid element type for testing
+    purposes.
 
     Returns:
-        A function that creates a default solid element description.
+        A function that creates a default solid element type.
     """
 
-    def _get_default_test_solid_element_description(
-        element_type: str = "2d_solid",
-    ):
-        """Return a default solid element description for testing purposes.
+    def _get_default_test_solid_element(element_type: str) -> type[Element]:
+        """Return a default solid element type for testing purposes.
 
         Args:
             element_type: The type of solid element to return.
 
         Returns:
-            A dictionary containing the solid element description parameters.
+            A type defining a solid element for testing purposes.
         """
 
-        if element_type == "2d_solid":
-            return {
-                "KINEM": "nonlinear",
-                "EAS": "none",
-                "THICK": 1.0,
-                "STRESS_STRAIN": "plane_strain",
-                "GP": [3, 3],
-            }
+        if element_type == "nurbs_2d":
+            return get_four_c_solid(
+                element_type,
+                n_nodes=9,
+                element_technology={
+                    "KINEM": "nonlinear",
+                    "EAS": "none",
+                    "THICK": 1.0,
+                    "STRESS_STRAIN": "plane_strain",
+                    "GP": [3, 3],
+                },
+            )
 
-        elif element_type == "2d_shell":
-            return {"type": "SHELL_KIRCHHOFF_LOVE_NURBS", "GP": [3, 3]}
+        elif element_type == "nurbs_3d":
+            return get_four_c_solid(
+                element_type, n_nodes=27, element_technology={"KINEM": "nonlinear"}
+            )
 
-        elif element_type == "3d_solid":
-            return {"KINEM": "nonlinear"}
+        elif element_type == "nurbs_shell":
+            return get_four_c_solid(
+                element_type,
+                n_nodes=9,
+                element_technology={"GP": [3, 3]},
+            )
 
         else:
             raise ValueError(f"Unknown solid element type: {element_type}")
 
-    return _get_default_test_solid_element_description
+    return _get_default_test_solid_element
 
 
 @pytest.fixture(scope="function")
diff --git a/tests/integration/test_integration_four_c.py b/tests/integration/test_integration_four_c.py
index 8ba761fb..5c70cefc 100644
--- a/tests/integration/test_integration_four_c.py
+++ b/tests/integration/test_integration_four_c.py
@@ -36,6 +36,7 @@
     Beam3rHerm2Line3,
     get_four_c_reissner_beam,
 )
+from beamme.four_c.element_solid import get_four_c_solid
 from beamme.four_c.header_functions import (
     add_result_description,
     set_beam_to_solid_meshtying,
@@ -223,8 +224,13 @@ def test_integration_four_c_nurbs_import(
     volume_set = GeometrySetNodes(geometry_type=bme.geo.volume)
     fix_set = GeometrySetNodes(geometry_type=bme.geo.surface)
     for i in range(3):
+        element_type = get_four_c_solid(
+            solid_type_string="nurbs_3d",
+            n_nodes=27,
+            element_technology=element_description,
+        )
         patch_set = add_splinepy_nurbs_to_mesh(
-            mesh, extruded, material=mat, data=element_description
+            mesh, element_type, extruded, material=mat
         )
         volume_set = volume_set + patch_set["vol"]
         fix_set = fix_set + patch_set["surf_w_min"]
diff --git a/tests/integration/test_integration_mesh_creation_functions_nurbs.py b/tests/integration/test_integration_mesh_creation_functions_nurbs.py
index 2056545a..4573f1a4 100644
--- a/tests/integration/test_integration_mesh_creation_functions_nurbs.py
+++ b/tests/integration/test_integration_mesh_creation_functions_nurbs.py
@@ -26,6 +26,7 @@
 
 from beamme.core.mesh import Mesh
 from beamme.core.rotation import Rotation
+from beamme.four_c.element_solid import get_four_c_solid
 from beamme.mesh_creation_functions.nurbs_generic import (
     add_geomdl_nurbs_to_mesh,
     add_splinepy_nurbs_to_mesh,
@@ -46,7 +47,7 @@
 
 
 def test_integration_mesh_creation_functions_nurbs_hollow_cylinder_segment_2d(
-    get_default_test_solid_element_description,
+    get_default_test_solid_element,
     get_default_test_solid_material,
     assert_results_close,
     get_corresponding_reference_file_path,
@@ -62,15 +63,12 @@ def test_integration_mesh_creation_functions_nurbs_hollow_cylinder_segment_2d(
     mesh = Mesh()
 
     # Create patch set
-    element_description = get_default_test_solid_element_description(
-        element_type="2d_solid"
-    )
-
+    element_type = get_default_test_solid_element("nurbs_2d")
     patch_set = add_geomdl_nurbs_to_mesh(
         mesh,
+        element_type,
         surf_obj,
         material=get_default_test_solid_material(material_type="st_venant_kirchhoff"),
-        data=element_description,
     )
 
     mesh.add(patch_set)
@@ -80,7 +78,7 @@ def test_integration_mesh_creation_functions_nurbs_hollow_cylinder_segment_2d(
 
 
 def test_integration_mesh_creation_functions_nurbs_flat_plate_2d(
-    get_default_test_solid_element_description,
+    get_default_test_solid_element,
     get_default_test_solid_material,
     assert_results_close,
     get_corresponding_reference_file_path,
@@ -97,12 +95,8 @@ def test_integration_mesh_creation_functions_nurbs_flat_plate_2d(
     mat = get_default_test_solid_material(material_type="2d_shell")
 
     # Create patch set
-    element_description = get_default_test_solid_element_description(
-        element_type="2d_shell"
-    )
-    patch_set = add_geomdl_nurbs_to_mesh(
-        mesh, surf_obj, material=mat, data=element_description
-    )
+    element_type = get_default_test_solid_element("nurbs_shell")
+    patch_set = add_geomdl_nurbs_to_mesh(mesh, element_type, surf_obj, material=mat)
 
     mesh.add(patch_set)
 
@@ -111,7 +105,7 @@ def test_integration_mesh_creation_functions_nurbs_flat_plate_2d(
 
 
 def test_integration_mesh_creation_functions_nurbs_flat_plate_2d_splinepy(
-    get_default_test_solid_element_description,
+    get_default_test_solid_element,
     get_default_test_solid_material,
     assert_results_close,
     get_corresponding_reference_file_path,
@@ -131,12 +125,8 @@ def test_integration_mesh_creation_functions_nurbs_flat_plate_2d_splinepy(
     # Create the shell mesh
     mesh = Mesh()
     mat = get_default_test_solid_material(material_type="2d_shell")
-    element_description = get_default_test_solid_element_description(
-        element_type="2d_shell"
-    )
-    patch_set = add_splinepy_nurbs_to_mesh(
-        mesh, surf_obj, material=mat, data=element_description
-    )
+    element_type = get_default_test_solid_element("nurbs_shell")
+    patch_set = add_splinepy_nurbs_to_mesh(mesh, element_type, surf_obj, material=mat)
     mesh.add(patch_set)
     assert_results_close(
         get_corresponding_reference_file_path(
@@ -147,7 +137,7 @@ def test_integration_mesh_creation_functions_nurbs_flat_plate_2d_splinepy(
 
 
 def test_integration_mesh_creation_functions_nurbs_flat_plate_2d_splinepy_copy(
-    get_default_test_solid_element_description,
+    get_default_test_solid_element,
     get_default_test_solid_material,
     assert_results_close,
     get_corresponding_reference_file_path,
@@ -162,10 +152,8 @@ def test_integration_mesh_creation_functions_nurbs_flat_plate_2d_splinepy_copy(
     # Create mesh
     mesh = Mesh()
     mat = get_default_test_solid_material(material_type="2d_shell")
-    element_description = get_default_test_solid_element_description(
-        element_type="2d_shell"
-    )
-    add_splinepy_nurbs_to_mesh(mesh, surf_obj, material=mat, data=element_description)
+    element_type = get_default_test_solid_element("nurbs_shell")
+    add_splinepy_nurbs_to_mesh(mesh, element_type, surf_obj, material=mat)
 
     mesh_copy = mesh.copy()
     mesh_copy.translate([3, 0, 0])
@@ -176,7 +164,7 @@ def test_integration_mesh_creation_functions_nurbs_flat_plate_2d_splinepy_copy(
 
 
 def test_integration_mesh_creation_functions_nurbs_brick(
-    get_default_test_solid_element_description,
+    get_default_test_solid_element,
     get_default_test_solid_material,
     assert_results_close,
     get_corresponding_reference_file_path,
@@ -190,11 +178,12 @@ def test_integration_mesh_creation_functions_nurbs_brick(
     mesh = Mesh()
 
     # Create patch set
+    element_type = get_default_test_solid_element("nurbs_3d")
     patch_set = add_geomdl_nurbs_to_mesh(
         mesh,
+        element_type,
         vol_obj,
         material=get_default_test_solid_material(material_type="st_venant_kirchhoff"),
-        data=get_default_test_solid_element_description(element_type="3d_solid"),
     )
 
     mesh.add(patch_set)
@@ -204,7 +193,7 @@ def test_integration_mesh_creation_functions_nurbs_brick(
 
 
 def test_integration_mesh_creation_functions_nurbs_brick_splinepy(
-    get_default_test_solid_element_description,
+    get_default_test_solid_element,
     get_default_test_solid_material,
     assert_results_close,
     get_corresponding_reference_file_path,
@@ -227,11 +216,12 @@ def test_integration_mesh_creation_functions_nurbs_brick_splinepy(
     mesh = Mesh()
 
     # Create patch set
+    element_type = get_default_test_solid_element("nurbs_3d")
     patch_set = add_splinepy_nurbs_to_mesh(
         mesh,
+        element_type,
         vol_obj,
         material=get_default_test_solid_material(material_type="st_venant_kirchhoff"),
-        data=get_default_test_solid_element_description(element_type="3d_solid"),
     )
 
     mesh.add(patch_set)
@@ -246,7 +236,7 @@ def test_integration_mesh_creation_functions_nurbs_brick_splinepy(
 
 
 def test_integration_mesh_creation_functions_nurbs_rotation_nurbs_surface(
-    get_default_test_solid_element_description,
+    get_default_test_solid_element,
     get_default_test_solid_material,
     assert_results_close,
     get_corresponding_reference_file_path,
@@ -262,15 +252,12 @@ def test_integration_mesh_creation_functions_nurbs_rotation_nurbs_surface(
     mesh = Mesh()
 
     # Create patch set
-    element_description = get_default_test_solid_element_description(
-        element_type="2d_solid"
-    )
-
+    element_type = get_default_test_solid_element("nurbs_2d")
     patch_set = add_geomdl_nurbs_to_mesh(
         mesh,
+        element_type,
         surf_obj,
         material=get_default_test_solid_material(material_type="st_venant_kirchhoff"),
-        data=element_description,
     )
 
     mesh.add(patch_set)
@@ -282,7 +269,7 @@ def test_integration_mesh_creation_functions_nurbs_rotation_nurbs_surface(
 
 
 def test_integration_mesh_creation_functions_nurbs_translate_nurbs_surface(
-    get_default_test_solid_element_description,
+    get_default_test_solid_element,
     get_default_test_solid_material,
     assert_results_close,
     get_corresponding_reference_file_path,
@@ -296,16 +283,12 @@ def test_integration_mesh_creation_functions_nurbs_translate_nurbs_surface(
     mesh = Mesh()
 
     # Create patch set
-
-    element_description = get_default_test_solid_element_description(
-        element_type="2d_solid"
-    )
-
+    element_type = get_default_test_solid_element("nurbs_2d")
     patch_set = add_geomdl_nurbs_to_mesh(
         mesh,
+        element_type,
         surf_obj,
         material=get_default_test_solid_material(material_type="st_venant_kirchhoff"),
-        data=element_description,
     )
 
     mesh.add(patch_set)
@@ -317,7 +300,7 @@ def test_integration_mesh_creation_functions_nurbs_translate_nurbs_surface(
 
 
 def test_integration_mesh_creation_functions_nurbs_cylindrical_shell_sector(
-    get_default_test_solid_element_description,
+    get_default_test_solid_element,
     get_default_test_solid_material,
     assert_results_close,
     get_corresponding_reference_file_path,
@@ -333,11 +316,12 @@ def test_integration_mesh_creation_functions_nurbs_cylindrical_shell_sector(
     mesh = Mesh()
 
     # Create patch set
+    element_type = get_default_test_solid_element("nurbs_2d")
     patch_set = add_geomdl_nurbs_to_mesh(
         mesh,
+        element_type,
         surf_obj,
         material=get_default_test_solid_material(material_type="st_venant_kirchhoff"),
-        data=get_default_test_solid_element_description(element_type="2d_solid"),
     )
 
     mesh.add(patch_set)
@@ -347,7 +331,7 @@ def test_integration_mesh_creation_functions_nurbs_cylindrical_shell_sector(
 
 
 def test_integration_mesh_creation_functions_nurbs_couple_nurbs_meshes(
-    get_default_test_solid_element_description,
+    get_default_test_solid_element,
     get_default_test_solid_material,
     assert_results_close,
     get_corresponding_reference_file_path,
@@ -364,14 +348,8 @@ def test_integration_mesh_creation_functions_nurbs_couple_nurbs_meshes(
 
     # Create first patch set
     mat = get_default_test_solid_material(material_type="st_venant_kirchhoff")
-
-    element_description = get_default_test_solid_element_description(
-        element_type="2d_solid"
-    )
-
-    patch_set_1 = add_geomdl_nurbs_to_mesh(
-        mesh, surf_obj_1, material=mat, data=element_description
-    )
+    element_type = get_default_test_solid_element("nurbs_2d")
+    patch_set_1 = add_geomdl_nurbs_to_mesh(mesh, element_type, surf_obj_1, material=mat)
 
     mesh.add(patch_set_1)
 
@@ -382,9 +360,7 @@ def test_integration_mesh_creation_functions_nurbs_couple_nurbs_meshes(
         0.65, 1.46, np.pi / 3, n_ele_u=3, n_ele_v=2
     )
 
-    patch_set_2 = add_geomdl_nurbs_to_mesh(
-        mesh, surf_obj_2, material=mat, data=element_description
-    )
+    patch_set_2 = add_geomdl_nurbs_to_mesh(mesh, element_type, surf_obj_2, material=mat)
 
     mesh.add(patch_set_2)
 
@@ -395,7 +371,7 @@ def test_integration_mesh_creation_functions_nurbs_couple_nurbs_meshes(
 
 
 def test_integration_mesh_creation_functions_nurbs_sphere_surface(
-    get_default_test_solid_element_description,
+    get_default_test_solid_element,
     get_default_test_solid_material,
     assert_results_close,
     get_corresponding_reference_file_path,
@@ -409,15 +385,12 @@ def test_integration_mesh_creation_functions_nurbs_sphere_surface(
     surf_obj = create_nurbs_sphere_surface(1, n_ele_u=3, n_ele_v=2)
 
     # Create first patch set
-    element_description = get_default_test_solid_element_description(
-        element_type="2d_solid"
-    )
-
+    element_type = get_default_test_solid_element("nurbs_2d")
     patch_set = add_geomdl_nurbs_to_mesh(
         mesh,
+        element_type,
         surf_obj,
         material=get_default_test_solid_material(material_type="st_venant_kirchhoff"),
-        data=element_description,
     )
 
     mesh.add(patch_set)
@@ -427,7 +400,7 @@ def test_integration_mesh_creation_functions_nurbs_sphere_surface(
 
 
 def test_integration_mesh_creation_functions_nurbs_string_types(
-    get_default_test_solid_element_description,
+    get_default_test_solid_element,
     get_default_test_solid_material,
     assert_results_close,
     get_corresponding_reference_file_path,
@@ -442,11 +415,12 @@ def test_integration_mesh_creation_functions_nurbs_string_types(
     surf_obj = create_nurbs_flat_plate_2d(1, 3, n_ele_u=3, n_ele_v=2)
 
     # Create first patch set
+    element_type = get_default_test_solid_element("nurbs_2d")
     patch_set = add_geomdl_nurbs_to_mesh(
         mesh,
+        element_type,
         surf_obj,
         material=get_default_test_solid_material(material_type="st_venant_kirchhoff"),
-        data=get_default_test_solid_element_description(element_type="2d_solid"),
     )
 
     mesh.add(patch_set)
@@ -456,7 +430,7 @@ def test_integration_mesh_creation_functions_nurbs_string_types(
 
 
 def test_integration_mesh_creation_functions_nurbs_hemisphere_surface(
-    get_default_test_solid_element_description,
+    get_default_test_solid_element,
     get_default_test_solid_material,
     assert_results_close,
     get_corresponding_reference_file_path,
@@ -471,17 +445,15 @@ def test_integration_mesh_creation_functions_nurbs_hemisphere_surface(
 
     # Create first patch set
     mat = get_default_test_solid_material(material_type="st_venant_kirchhoff")
-    element_description = get_default_test_solid_element_description(
-        element_type="2d_solid"
-    )
+    element_type = get_default_test_solid_element("nurbs_2d")
 
     # Add the patch sets of every surface section of the hemisphere to the input file
     for surf in surfs:
         patch_set = add_geomdl_nurbs_to_mesh(
             mesh,
+            element_type,
             surf,
             material=mat,
-            data=element_description,
         )
 
         mesh.add(patch_set)
@@ -491,7 +463,7 @@ def test_integration_mesh_creation_functions_nurbs_hemisphere_surface(
 
 
 def test_integration_mesh_creation_functions_nurbs_torus_surface(
-    get_default_test_solid_element_description,
+    get_default_test_solid_element,
     get_default_test_solid_material,
     assert_results_close,
     get_corresponding_reference_file_path,
@@ -506,17 +478,15 @@ def test_integration_mesh_creation_functions_nurbs_torus_surface(
 
     # Define element description
     mat = get_default_test_solid_material(material_type="st_venant_kirchhoff")
-    element_description = get_default_test_solid_element_description(
-        element_type="2d_solid"
-    )
+    element_type = get_default_test_solid_element("nurbs_2d")
 
     # Add the patch sets of every surface section of the torus to the input file
     for surf in surfs:
         patch_set = add_geomdl_nurbs_to_mesh(
             mesh,
+            element_type,
             surf,
             material=mat,
-            data=element_description,
         )
 
         mesh.add(patch_set)
@@ -543,7 +513,8 @@ def test_integration_mesh_creation_functions_nurbs_empty_knot_spans(
     # Create mesh
     mesh = Mesh()
     mat = get_default_test_solid_material(material_type="st_venant_kirchhoff")
-    patch_set = add_splinepy_nurbs_to_mesh(mesh, pipe, material=mat)
+    element_type = get_four_c_solid("nurbs_3d", n_nodes=27)
+    patch_set = add_splinepy_nurbs_to_mesh(mesh, element_type, pipe, material=mat)
     mesh.add(patch_set)
     mesh.couple_nodes(reuse_matching_nodes=True)