diff --git a/include/fesa/elements/mitc4_shell.h b/include/fesa/elements/mitc4_shell.h
index 9b49f38..f0dc0fd 100644
--- a/include/fesa/elements/mitc4_shell.h
+++ b/include/fesa/elements/mitc4_shell.h
@@ -7,6 +7,7 @@
 #include "fesa/core/status.h"
 #include "fesa/math/matrix.h"
 #include "fesa/math/vector.h"
+#include "fesa/math/vector3.h"
 #include "fesa/model/model_types.h"
 
 namespace fesa {
@@ -132,8 +133,6 @@ class Mitc4Shell {
       const Vector& global_element_displacement24) const;
 
  private:
-  using Vector3 = std::array<double, 3>;
-
   /// @brief Stores covariant, reciprocal, frame, and Jacobian data at one
   /// point.
   struct GeometryData {
diff --git a/include/fesa/math/vector3.h b/include/fesa/math/vector3.h
index c27fa6c..ae4e16e 100644
--- a/include/fesa/math/vector3.h
+++ b/include/fesa/math/vector3.h
@@ -21,6 +21,10 @@ class Vector3 {
   constexpr Vector3(double x, double y, double z) noexcept
       : components_{{x, y, z}} {}
 
+  /// @brief Copies components from an existing array-backed carrier.
+  explicit constexpr Vector3(const std::array<double, 3>& components) noexcept
+      : components_{components} {}
+
   /// @brief Returns the first component.
   constexpr double X() const noexcept { return components_[0]; }
 
@@ -36,6 +40,11 @@ class Vector3 {
     return components_[index];
   }
 
+  /// @brief Returns the immutable array-backed component carrier.
+  constexpr const std::array<double, 3>& Components() const noexcept {
+    return components_;
+  }
+
   /// @brief Adds corresponding vector components.
   constexpr Vector3 operator+(const Vector3& rhs) const noexcept {
     return Vector3{X() + rhs.X(), Y() + rhs.Y(), Z() + rhs.Z()};
@@ -51,6 +60,22 @@ class Vector3 {
     return Vector3{X() * scalar, Y() * scalar, Z() * scalar};
   }
 
+  /// @brief Divides every component by a scalar.
+  constexpr Vector3 operator/(double scalar) const noexcept {
+    return Vector3{X() / scalar, Y() / scalar, Z() / scalar};
+  }
+
+  /// @brief Multiplies every component with the scalar as the left operand.
+  friend constexpr Vector3 operator*(double scalar,
+                                     const Vector3& rhs) noexcept {
+    return Vector3{scalar * rhs.X(), scalar * rhs.Y(), scalar * rhs.Z()};
+  }
+
+  /// @brief Compares every component exactly.
+  constexpr bool operator==(const Vector3& rhs) const noexcept {
+    return X() == rhs.X() && Y() == rhs.Y() && Z() == rhs.Z();
+  }
+
   /// @brief Computes the Euclidean dot product with rhs.
   double Dot(const Vector3& rhs) const noexcept {
     return X() * rhs.X() + Y() * rhs.Y() + Z() * rhs.Z();
@@ -63,7 +88,7 @@ class Vector3 {
   }
 
   /// @brief Computes the Euclidean norm.
-  double Norm() const noexcept { return std::sqrt(Dot(*this)); }
+  double Norm() const noexcept { return std::hypot(X(), Y(), Z()); }
 
   /// @brief Returns a unit vector when the norm is usable.
   /// @return Empty when the norm is exactly zero or nonfinite.
diff --git a/src/fesa/elements/euler_beam_3d.cpp b/src/fesa/elements/euler_beam_3d.cpp
index 8d0a2ca..51c87dd 100644
--- a/src/fesa/elements/euler_beam_3d.cpp
+++ b/src/fesa/elements/euler_beam_3d.cpp
@@ -9,6 +9,8 @@
 #include <utility>
 #include <vector>
 
+#include "fesa/math/vector3.h"
+
 namespace fesa {
 namespace {
 
@@ -17,26 +19,6 @@ constexpr std::size_t kGeneralizedComponentCount = 4U;
 constexpr double kGeometryTolerance = 1.0e-12;
 constexpr double kStiffnessInvariantTolerance = 1.0e-12;
 
-using Vector3 = std::array<double, 3>;
-
-double Norm(const Vector3& value) {
-  return std::hypot(value[0], value[1], value[2]);
-}
-
-double Dot(const Vector3& lhs, const Vector3& rhs) {
-  return lhs[0] * rhs[0] + lhs[1] * rhs[1] + lhs[2] * rhs[2];
-}
-
-Vector3 Cross(const Vector3& lhs, const Vector3& rhs) {
-  return {lhs[1] * rhs[2] - lhs[2] * rhs[1], lhs[2] * rhs[0] - lhs[0] * rhs[2],
-          lhs[0] * rhs[1] - lhs[1] * rhs[0]};
-}
-
-bool IsFinite(const Vector3& value) {
-  return std::isfinite(value[0]) && std::isfinite(value[1]) &&
-         std::isfinite(value[2]);
-}
-
 std::string ElementIdentity(const Node& first_node, const Node& second_node) {
   return first_node.source_id.instance_name + ":" +
          first_node.source_id.source_label_text + "-" +
@@ -233,13 +215,13 @@ Result<EulerBeam3D> EulerBeam3D::Create(const Node& first_node,
                                         const GeneralBeamSection& section,
                                         const LinearElasticMaterial& material) {
   const std::string identity = ElementIdentity(first_node, second_node);
-  const Vector3& first = first_node.coordinates;
-  const Vector3& second = second_node.coordinates;
-  const Vector3 delta = {second[0] - first[0], second[1] - first[1],
-                         second[2] - first[2]};
-  const double length = Norm(delta);
-  const double coordinate_scale = (std::max)({1.0, Norm(first), Norm(second)});
-  if (!IsFinite(first) || !IsFinite(second) || !IsFinite(delta) ||
+  const Vector3 first{first_node.coordinates};
+  const Vector3 second{second_node.coordinates};
+  const Vector3 delta = second - first;
+  const double length = delta.Norm();
+  const double coordinate_scale =
+      (std::max)({1.0, first.Norm(), second.Norm()});
+  if (!first.IsFinite() || !second.IsFinite() || !delta.IsFinite() ||
       !std::isfinite(length) || !std::isfinite(coordinate_scale) ||
       !(length > kGeometryTolerance * coordinate_scale)) {
     return ModelFailure(
@@ -247,15 +229,13 @@ Result<EulerBeam3D> EulerBeam3D::Create(const Node& first_node,
         "Beam length must exceed the scale-aware geometry threshold.");
   }
 
-  const Vector3 ex = {delta[0] / length, delta[1] / length, delta[2] / length};
-  const Vector3& guide = section.first_axis;
-  const double guide_norm = Norm(guide);
-  const double guide_projection = Dot(guide, ex);
-  const Vector3 ey_trial = {guide[0] - guide_projection * ex[0],
-                            guide[1] - guide_projection * ex[1],
-                            guide[2] - guide_projection * ex[2]};
-  const double ey_trial_norm = Norm(ey_trial);
-  if (!IsFinite(guide) || !std::isfinite(guide_norm) || !IsFinite(ey_trial) ||
+  const Vector3 ex = delta / length;
+  const Vector3 guide{section.first_axis};
+  const double guide_norm = guide.Norm();
+  const double guide_projection = guide.Dot(ex);
+  const Vector3 ey_trial = guide - guide_projection * ex;
+  const double ey_trial_norm = ey_trial.Norm();
+  if (!guide.IsFinite() || !std::isfinite(guide_norm) || !ey_trial.IsFinite() ||
       !std::isfinite(ey_trial_norm) ||
       !(ey_trial_norm > kGeometryTolerance * (std::max)(1.0, guide_norm))) {
     return ModelFailure(
@@ -323,9 +303,8 @@ Result<EulerBeam3D> EulerBeam3D::Create(const Node& first_node,
         "Derived beam stiffness coefficients must be finite and positive.");
   }
 
-  const Vector3 ey = {ey_trial[0] / ey_trial_norm, ey_trial[1] / ey_trial_norm,
-                      ey_trial[2] / ey_trial_norm};
-  const Vector3 ez = Cross(ex, ey);
+  const Vector3 ey = ey_trial / ey_trial_norm;
+  const Vector3 ez = ex.Cross(ey);
   // Rows map global vectors to the approved right-handed local (ex,ey,ez)
   // basis.
   const std::array<double, 9> rotation = {ex[0], ex[1], ex[2], ey[0], ey[1],
diff --git a/src/fesa/elements/mitc4_shell.cpp b/src/fesa/elements/mitc4_shell.cpp
index d9634b4..3b28f8c 100644
--- a/src/fesa/elements/mitc4_shell.cpp
+++ b/src/fesa/elements/mitc4_shell.cpp
@@ -10,8 +10,6 @@
 namespace fesa {
 namespace {
 
-using Vector3 = std::array<double, 3>;
-
 constexpr std::size_t kNodeCount = 4U;
 constexpr std::size_t kPhysicalDofsPerNode = 5U;
 constexpr std::size_t kGlobalDofsPerNode = 6U;
@@ -20,46 +18,11 @@ constexpr std::size_t kGlobalDofCount = 24U;
 constexpr double kShearCorrection = 5.0 / 6.0;
 constexpr double kFrameTolerance = 1.0e-12;
 
-Vector3 Add(const Vector3& left, const Vector3& right) {
-  return {left[0] + right[0], left[1] + right[1], left[2] + right[2]};
-}
-
-Vector3 Subtract(const Vector3& left, const Vector3& right) {
-  return {left[0] - right[0], left[1] - right[1], left[2] - right[2]};
-}
-
-Vector3 Scale(double factor, const Vector3& value) {
-  return {factor * value[0], factor * value[1], factor * value[2]};
-}
-
-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::hypot(value[0], value[1], value[2]);
-}
-
-bool IsFinite(const Vector3& value) {
-  return std::all_of(value.begin(), value.end(),
-                     [](double component) { return std::isfinite(component); });
-}
-
-Vector3 Normalized(const Vector3& value) {
-  return Scale(1.0 / Norm(value), value);
-}
-
 Vector3 WeightedSum(const std::array<double, kNodeCount>& weights,
                     const std::array<Vector3, kNodeCount>& values) {
   Vector3 result{};
   for (std::size_t node = 0U; node < kNodeCount; ++node) {
-    result = Add(result, Scale(weights[node], values[node]));
+    result = result + weights[node] * values[node];
   }
   return result;
 }
@@ -68,7 +31,7 @@ Vector3 DerivativeSum(const std::array<double, kNodeCount>& derivatives,
                       const std::array<Vector3, kNodeCount>& values) {
   std::array<Vector3, kNodeCount> relative{};
   for (std::size_t node = 0U; node < kNodeCount; ++node) {
-    relative[node] = Subtract(values[node], values[0]);
+    relative[node] = values[node] - values[0];
   }
   return WeightedSum(derivatives, relative);
 }
@@ -86,18 +49,18 @@ std::array<Vector3, kNodeCount> NodalTangentsA(
   std::array<Vector3, kNodeCount> tangents{};
   for (std::size_t node = 0U; node < kNodeCount; ++node) {
     std::size_t selected = 0U;
-    double alignment = std::abs(Dot(global_axes[0], directors[node]));
+    double alignment = std::abs(global_axes[0].Dot(directors[node]));
     for (std::size_t axis = 1U; axis < global_axes.size(); ++axis) {
-      const double candidate =
-          std::abs(Dot(global_axes[axis], directors[node]));
+      const double candidate = std::abs(global_axes[axis].Dot(directors[node]));
       if (candidate < alignment) {
         selected = axis;
         alignment = candidate;
       }
     }
-    tangents[node] = Normalized(Subtract(
-        global_axes[selected],
-        Scale(Dot(global_axes[selected], directors[node]), directors[node])));
+    const Vector3 tangent_candidate =
+        global_axes[selected] -
+        global_axes[selected].Dot(directors[node]) * directors[node];
+    tangents[node] = (1.0 / tangent_candidate.Norm()) * tangent_candidate;
   }
   return tangents;
 }
@@ -108,7 +71,7 @@ std::array<Vector3, kNodeCount> NodalTangentsB(
     const std::array<Vector3, kNodeCount>& tangent_a) {
   std::array<Vector3, kNodeCount> tangents{};
   for (std::size_t node = 0U; node < kNodeCount; ++node) {
-    tangents[node] = Cross(directors[node], tangent_a[node]);
+    tangents[node] = directors[node].Cross(tangent_a[node]);
   }
   return tangents;
 }
@@ -144,9 +107,8 @@ std::array<std::array<double, 3>, 3> CovariantStrainColumn(
   std::array<std::array<double, 3>, 3> strain{};
   for (std::size_t first = 0U; first < 3U; ++first) {
     for (std::size_t second = 0U; second < 3U; ++second) {
-      strain[first][second] =
-          0.5 * (Dot(covariant[first], derivatives[second]) +
-                 Dot(covariant[second], derivatives[first]));
+      strain[first][second] = 0.5 * (covariant[first].Dot(derivatives[second]) +
+                                     covariant[second].Dot(derivatives[first]));
     }
   }
   // Thickness stretch is excluded from the five-component shell law.
@@ -190,11 +152,13 @@ double FrameComponent(const Vector3& left,
 std::array<double, 5> LocalEngineeringComponents(
     const std::array<std::array<double, 3>, 3>& tensor,
     const Mitc4LocalFrame& frame) {
-  return {FrameComponent(frame.e1, tensor, frame.e1),
-          FrameComponent(frame.e2, tensor, frame.e2),
-          2.0 * FrameComponent(frame.e1, tensor, frame.e2),
-          2.0 * FrameComponent(frame.e1, tensor, frame.e3),
-          2.0 * FrameComponent(frame.e2, tensor, frame.e3)};
+  const Vector3 e1{frame.e1};
+  const Vector3 e2{frame.e2};
+  const Vector3 e3{frame.e3};
+  return {FrameComponent(e1, tensor, e1), FrameComponent(e2, tensor, e2),
+          2.0 * FrameComponent(e1, tensor, e2),
+          2.0 * FrameComponent(e1, tensor, e3),
+          2.0 * FrameComponent(e2, tensor, e3)};
 }
 
 Matrix ScaledMatrix(const Matrix& source, double factor) {
@@ -287,8 +251,8 @@ Result<Mitc4Shell> Mitc4Shell::Create(
 
   std::array<Vector3, kNodeCount> coordinates{};
   for (std::size_t node = 0U; node < kNodeCount; ++node) {
-    coordinates[node] = nodes[node]->coordinates;
-    if (!IsFinite(coordinates[node])) {
+    coordinates[node] = Vector3{nodes[node]->coordinates};
+    if (!coordinates[node].IsFinite()) {
       return ModelFailure("invalid-shell-geometry", nodes[node]->location,
                           identity, "MITC4 node coordinates must be finite.");
     }
@@ -301,9 +265,12 @@ Result<Mitc4Shell> Mitc4Shell::Create(
     }
   }
 
+  const std::array<Vector3, kNodeCount> directors{
+      Vector3{initial_directors[0U]}, Vector3{initial_directors[1U]},
+      Vector3{initial_directors[2U]}, Vector3{initial_directors[3U]}};
   for (std::size_t node = 0U; node < kNodeCount; ++node) {
-    const double director_norm = Norm(initial_directors[node]);
-    if (!IsFinite(initial_directors[node]) || !std::isfinite(director_norm) ||
+    const double director_norm = directors[node].Norm();
+    if (!directors[node].IsFinite() || !std::isfinite(director_norm) ||
         std::abs(director_norm - 1.0) > kFrameTolerance) {
       return ModelFailure(
           "invalid-shell-director", nodes[node]->location, identity,
@@ -329,28 +296,28 @@ Result<Mitc4Shell> Mitc4Shell::Create(
       DerivativeSum(center_shape.xi_derivatives, coordinates);
   const Vector3 center_eta =
       DerivativeSum(center_shape.eta_derivatives, coordinates);
-  const Vector3 center_area = Cross(center_xi, center_eta);
-  const double center_measure = Norm(center_area);
-  if (!IsFinite(center_area) || !std::isfinite(center_measure) ||
+  const Vector3 center_area = center_xi.Cross(center_eta);
+  const double center_measure = center_area.Norm();
+  if (!center_area.IsFinite() || !std::isfinite(center_measure) ||
       !(center_measure > 0.0)) {
     return ModelFailure(
         "invalid-shell-geometry", nodes[0]->location, identity,
         "MITC4 center surface basis must be finite and nonzero.");
   }
-  const Vector3 normal_candidate = Scale(1.0 / center_measure, center_area);
-  if (std::any_of(initial_directors.begin(), initial_directors.end(),
+  const Vector3 normal_candidate = (1.0 / center_measure) * center_area;
+  if (std::any_of(directors.begin(), directors.end(),
                   [&normal_candidate](const Vector3& director) {
-                    return !(Dot(normal_candidate, director) > 0.0);
+                    return !(normal_candidate.Dot(director) > 0.0);
                   })) {
     return ModelFailure(
         "invalid-shell-director", nodes[0]->location, identity,
         "MITC4 directors must follow the source-order positive face.");
   }
 
-  const auto tangent_a = NodalTangentsA(initial_directors);
-  const auto tangent_b = NodalTangentsB(initial_directors, tangent_a);
+  const auto tangent_a = NodalTangentsA(directors);
+  const auto tangent_b = NodalTangentsB(directors, tangent_a);
   Mitc4Shell shell{coordinates,
-                   initial_directors,
+                   directors,
                    tangent_a,
                    tangent_b,
                    normal_candidate,
@@ -805,48 +772,47 @@ bool Mitc4Shell::EvaluateGeometry(double xi, double eta, double zeta,
   const Vector3 director_value = WeightedSum(shape.values, directors_);
   const double half_thickness = 0.5 * thickness_;
 
-  result.covariant[0] =
-      Add(midsurface_xi, Scale(half_thickness * zeta, director_xi));
-  result.covariant[1] =
-      Add(midsurface_eta, Scale(half_thickness * zeta, director_eta));
-  result.covariant[2] = Scale(half_thickness, director_value);
+  result.covariant[0] = midsurface_xi + half_thickness * zeta * director_xi;
+  result.covariant[1] = midsurface_eta + half_thickness * zeta * director_eta;
+  result.covariant[2] = half_thickness * director_value;
   result.jacobian =
-      Dot(result.covariant[0], Cross(result.covariant[1], result.covariant[2]));
-  if (!IsFinite(result.covariant[0]) || !IsFinite(result.covariant[1]) ||
-      !IsFinite(result.covariant[2]) || !std::isfinite(result.jacobian) ||
+      result.covariant[0].Dot(result.covariant[1].Cross(result.covariant[2]));
+  if (!result.covariant[0].IsFinite() || !result.covariant[1].IsFinite() ||
+      !result.covariant[2].IsFinite() || !std::isfinite(result.jacobian) ||
       !(result.jacobian > 0.0)) {
     return false;
   }
 
-  result.reciprocal[0] = Scale(1.0 / result.jacobian,
-                               Cross(result.covariant[1], result.covariant[2]));
-  result.reciprocal[1] = Scale(1.0 / result.jacobian,
-                               Cross(result.covariant[2], result.covariant[0]));
-  result.reciprocal[2] = Scale(1.0 / result.jacobian,
-                               Cross(result.covariant[0], result.covariant[1]));
+  result.reciprocal[0] =
+      (1.0 / result.jacobian) * result.covariant[1].Cross(result.covariant[2]);
+  result.reciprocal[1] =
+      (1.0 / result.jacobian) * result.covariant[2].Cross(result.covariant[0]);
+  result.reciprocal[2] =
+      (1.0 / result.jacobian) * result.covariant[0].Cross(result.covariant[1]);
 
-  const Vector3 area = Cross(midsurface_xi, midsurface_eta);
-  const double director_norm = Norm(director_value);
-  if (!IsFinite(area) || !IsFinite(director_value) ||
+  const Vector3 area = midsurface_xi.Cross(midsurface_eta);
+  const double director_norm = director_value.Norm();
+  if (!area.IsFinite() || !director_value.IsFinite() ||
       !std::isfinite(director_norm) || !(director_norm > 0.0) ||
-      !(Dot(area, normal_candidate_) > 0.0)) {
+      !(area.Dot(normal_candidate_) > 0.0)) {
     return false;
   }
-  result.frame.e3 = Scale(1.0 / director_norm, director_value);
-  if (!(Dot(area, result.frame.e3) > 0.0)) {
+  const Vector3 e3 = (1.0 / director_norm) * director_value;
+  result.frame.e3 = e3.Components();
+  if (!(area.Dot(e3) > 0.0)) {
     return false;
   }
-  const Vector3 e1_candidate =
-      Subtract(midsurface_xi,
-               Scale(Dot(midsurface_xi, result.frame.e3), result.frame.e3));
-  const double e1_norm = Norm(e1_candidate);
-  if (!IsFinite(e1_candidate) || !std::isfinite(e1_norm) || !(e1_norm > 0.0)) {
+  const Vector3 e1_candidate = midsurface_xi - midsurface_xi.Dot(e3) * e3;
+  const double e1_norm = e1_candidate.Norm();
+  if (!e1_candidate.IsFinite() || !std::isfinite(e1_norm) || !(e1_norm > 0.0)) {
     return false;
   }
-  result.frame.e1 = Scale(1.0 / e1_norm, e1_candidate);
-  result.frame.e2 = Cross(result.frame.e3, result.frame.e1);
-  return IsFinite(result.reciprocal[0]) && IsFinite(result.reciprocal[1]) &&
-         IsFinite(result.reciprocal[2]) && IsFinite(result.frame.e2);
+  const Vector3 e1 = (1.0 / e1_norm) * e1_candidate;
+  const Vector3 e2 = e3.Cross(e1);
+  result.frame.e1 = e1.Components();
+  result.frame.e2 = e2.Components();
+  return result.reciprocal[0].IsFinite() && result.reciprocal[1].IsFinite() &&
+         result.reciprocal[2].IsFinite() && e2.IsFinite();
 }
 
 std::array<std::array<Vector3, 3>, 20> Mitc4Shell::BasisDerivatives(
@@ -861,24 +827,24 @@ std::array<std::array<Vector3, 3>, 20> Mitc4Shell::BasisDerivatives(
     const std::size_t offset = node * kPhysicalDofsPerNode;
     for (std::size_t component = 0U; component < 3U; ++component) {
       derivatives[offset + component][0] =
-          Scale(shape.xi_derivatives[node], global_axes[component]);
+          shape.xi_derivatives[node] * global_axes[component];
       derivatives[offset + component][1] =
-          Scale(shape.eta_derivatives[node], global_axes[component]);
+          shape.eta_derivatives[node] * global_axes[component];
     }
 
-    const Vector3 alpha_direction = Scale(-half_thickness, tangent_b_[node]);
+    const Vector3 alpha_direction = -half_thickness * tangent_b_[node];
     derivatives[offset + 3U][0] =
-        Scale(zeta * shape.xi_derivatives[node], alpha_direction);
+        zeta * shape.xi_derivatives[node] * alpha_direction;
     derivatives[offset + 3U][1] =
-        Scale(zeta * shape.eta_derivatives[node], alpha_direction);
-    derivatives[offset + 3U][2] = Scale(shape.values[node], alpha_direction);
+        zeta * shape.eta_derivatives[node] * alpha_direction;
+    derivatives[offset + 3U][2] = shape.values[node] * alpha_direction;
 
-    const Vector3 beta_direction = Scale(half_thickness, tangent_a_[node]);
+    const Vector3 beta_direction = half_thickness * tangent_a_[node];
     derivatives[offset + 4U][0] =
-        Scale(zeta * shape.xi_derivatives[node], beta_direction);
+        zeta * shape.xi_derivatives[node] * beta_direction;
     derivatives[offset + 4U][1] =
-        Scale(zeta * shape.eta_derivatives[node], beta_direction);
-    derivatives[offset + 4U][2] = Scale(shape.values[node], beta_direction);
+        zeta * shape.eta_derivatives[node] * beta_direction;
+    derivatives[offset + 4U][2] = shape.values[node] * beta_direction;
   }
   return derivatives;
 }
diff --git a/src/fesa/model/shell_geometry.cpp b/src/fesa/model/shell_geometry.cpp
index 1277b41..9e152fe 100644
--- a/src/fesa/model/shell_geometry.cpp
+++ b/src/fesa/model/shell_geometry.cpp
@@ -5,11 +5,11 @@
 #include <tuple>
 #include <utility>
 
+#include "fesa/math/vector3.h"
+
 namespace fesa {
 namespace {
 
-using Vector3 = std::array<double, 3>;
-
 constexpr std::array<double, 4> kXiSigns{-1.0, 1.0, 1.0, -1.0};
 constexpr std::array<double, 4> kEtaSigns{-1.0, -1.0, 1.0, 1.0};
 
@@ -25,37 +25,6 @@ struct ElementWork {
   double area_weight;
 };
 
-Vector3 Add(const Vector3& left, const Vector3& right) {
-  return {left[0] + right[0], left[1] + right[1], left[2] + right[2]};
-}
-
-Vector3 Subtract(const Vector3& left, const Vector3& right) {
-  return {left[0] - right[0], left[1] - right[1], left[2] - right[2]};
-}
-
-Vector3 Scale(double factor, const Vector3& value) {
-  return {factor * value[0], factor * value[1], factor * value[2]};
-}
-
-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::hypot(value[0], value[1], value[2]);
-}
-
-bool IsFinite(const Vector3& value) {
-  return std::all_of(value.begin(), value.end(),
-                     [](double component) { return std::isfinite(component); });
-}
-
 ShapeData ShapeDataAt(double xi, double eta) {
   ShapeData data{};
   for (std::size_t node = 0U; node < 4U; ++node) {
@@ -73,7 +42,7 @@ Vector3 WeightedSum(const std::array<double, 4>& weights,
                     const std::array<Vector3, 4>& values) {
   Vector3 result{};
   for (std::size_t node = 0U; node < values.size(); ++node) {
-    result = Add(result, Scale(weights[node], values[node]));
+    result = result + weights[node] * values[node];
   }
   return result;
 }
@@ -85,7 +54,7 @@ Vector3 DerivativeSum(const std::array<double, 4>& derivatives,
   // numerical cancellation without changing the covariant tangent.
   std::array<Vector3, 4> relative{};
   for (std::size_t node = 0U; node < coordinates.size(); ++node) {
-    relative[node] = Subtract(coordinates[node], coordinates[0]);
+    relative[node] = coordinates[node] - coordinates[0];
   }
   return WeightedSum(derivatives, relative);
 }
@@ -106,7 +75,7 @@ bool SameCoordinates(const Vector3& left, const Vector3& right) {
 
 double Orientation(const Vector3& first, const Vector3& second,
                    const Vector3& third, const Vector3& normal) {
-  return Dot(Cross(Subtract(second, first), Subtract(third, first)), normal);
+  return (second - first).Cross(third - first).Dot(normal);
 }
 
 bool HasOppositeSigns(double first, double second) {
@@ -197,8 +166,8 @@ Result<ShellGeometry> PreprocessShellGeometry(
             "Shell geometry references an unavailable internal node.");
       }
       current.coordinates[local_node] =
-          nodes[element.node_indices[local_node]].coordinates;
-      if (!IsFinite(current.coordinates[local_node])) {
+          Vector3{nodes[element.node_indices[local_node]].coordinates};
+      if (!current.coordinates[local_node].IsFinite()) {
         return GeometryFailure(
             "invalid-shell-geometry", element.location, "ELEMENT",
             element.source_id.source_label_text,
@@ -222,17 +191,17 @@ Result<ShellGeometry> PreprocessShellGeometry(
         DerivativeSum(center.xi_derivatives, current.coordinates);
     const Vector3 center_eta =
         DerivativeSum(center.eta_derivatives, current.coordinates);
-    const Vector3 center_cross = Cross(center_xi, center_eta);
-    const double center_measure = Norm(center_cross);
-    if (!IsFinite(center_xi) || !IsFinite(center_eta) ||
-        !IsFinite(center_cross) || !std::isfinite(center_measure) ||
+    const Vector3 center_cross = center_xi.Cross(center_eta);
+    const double center_measure = center_cross.Norm();
+    if (!center_xi.IsFinite() || !center_eta.IsFinite() ||
+        !center_cross.IsFinite() || !std::isfinite(center_measure) ||
         !(center_measure > 0.0)) {
       return GeometryFailure(
           "invalid-shell-geometry", element.location, "ELEMENT",
           element.source_id.source_label_text,
           "Shell center has no finite nonzero normal candidate.");
     }
-    current.normal = Scale(1.0 / center_measure, center_cross);
+    current.normal = (1.0 / center_measure) * center_cross;
 
     if (SegmentsProperlyIntersect(
             current.coordinates[0], current.coordinates[1],
@@ -253,11 +222,11 @@ Result<ShellGeometry> PreprocessShellGeometry(
           DerivativeSum(shape.xi_derivatives, current.coordinates);
       const Vector3 tangent_eta =
           DerivativeSum(shape.eta_derivatives, current.coordinates);
-      const Vector3 area_vector = Cross(tangent_xi, tangent_eta);
-      const double measure = Norm(area_vector);
-      if (!IsFinite(tangent_xi) || !IsFinite(tangent_eta) ||
-          !IsFinite(area_vector) || !std::isfinite(measure) ||
-          !(measure > 0.0) || !(Dot(area_vector, current.normal) > 0.0)) {
+      const Vector3 area_vector = tangent_xi.Cross(tangent_eta);
+      const double measure = area_vector.Norm();
+      if (!tangent_xi.IsFinite() || !tangent_eta.IsFinite() ||
+          !area_vector.IsFinite() || !std::isfinite(measure) ||
+          !(measure > 0.0) || !(area_vector.Dot(current.normal) > 0.0)) {
         return GeometryFailure("invalid-shell-geometry", element.location,
                                "ELEMENT", element.source_id.source_label_text,
                                "Shell surface is zero-area or locally reversed "
@@ -272,8 +241,8 @@ Result<ShellGeometry> PreprocessShellGeometry(
     }
     current.area_weight = area_weight;
     work.push_back(current);
-    geometry.element_data.push_back(
-        {static_cast<EntityIndex>(element_index), current.normal, area_weight});
+    geometry.element_data.push_back({static_cast<EntityIndex>(element_index),
+                                     current.normal.Components(), area_weight});
   }
 
   std::vector<std::vector<std::size_t>> incident(nodes.size());
@@ -301,8 +270,8 @@ Result<ShellGeometry> PreprocessShellGeometry(
     for (std::size_t first = 0U; first < node_incident.size(); ++first) {
       for (std::size_t second = first + 1U; second < node_incident.size();
            ++second) {
-        const double pair_dot = Dot(work[node_incident[first]].normal,
-                                    work[node_incident[second]].normal);
+        const double pair_dot = work[node_incident[first]].normal.Dot(
+            work[node_incident[second]].normal);
         if (!std::isfinite(pair_dot) || !(pair_dot > 0.0)) {
           return GeometryFailure("opposed-incident-normal",
                                  nodes[node_index].location, "NODE",
@@ -320,52 +289,53 @@ Result<ShellGeometry> PreprocessShellGeometry(
     }
     Vector3 director_sum{};
     for (const std::size_t element_index : node_incident) {
-      director_sum = Add(director_sum,
-                         Scale(work[element_index].area_weight / maximum_weight,
-                               work[element_index].normal));
+      director_sum =
+          director_sum + (work[element_index].area_weight / maximum_weight) *
+                             work[element_index].normal;
     }
-    const double director_norm = Norm(director_sum);
-    if (!IsFinite(director_sum) || !std::isfinite(director_norm) ||
+    const double director_norm = director_sum.Norm();
+    if (!director_sum.IsFinite() || !std::isfinite(director_norm) ||
         !(director_norm > 0.0)) {
       return GeometryFailure(
           "invalid-shell-director", nodes[node_index].location, "NODE",
           nodes[node_index].source_id.source_label_text,
           "Area-weighted shell director is nonfinite or zero.");
     }
-    const Vector3 director = Scale(1.0 / director_norm, director_sum);
+    const Vector3 director = (1.0 / director_norm) * director_sum;
 
     const std::array<Vector3, 3> global_axes{
         Vector3{1.0, 0.0, 0.0}, Vector3{0.0, 1.0, 0.0}, Vector3{0.0, 0.0, 1.0}};
     std::size_t selected_axis = 0U;
-    double selected_alignment = std::abs(Dot(global_axes[0], director));
+    double selected_alignment = std::abs(global_axes[0].Dot(director));
     for (std::size_t axis = 1U; axis < global_axes.size(); ++axis) {
-      const double alignment = std::abs(Dot(global_axes[axis], director));
+      const double alignment = std::abs(global_axes[axis].Dot(director));
       if (alignment < selected_alignment) {
         selected_alignment = alignment;
         selected_axis = axis;
       }
     }
     const Vector3 tangent_candidate =
-        Subtract(global_axes[selected_axis],
-                 Scale(Dot(global_axes[selected_axis], director), director));
-    const double tangent_norm = Norm(tangent_candidate);
-    if (!IsFinite(tangent_candidate) || !std::isfinite(tangent_norm) ||
+        global_axes[selected_axis] -
+        global_axes[selected_axis].Dot(director) * director;
+    const double tangent_norm = tangent_candidate.Norm();
+    if (!tangent_candidate.IsFinite() || !std::isfinite(tangent_norm) ||
         !(tangent_norm > 0.0)) {
       return GeometryFailure(
           "invalid-shell-director", nodes[node_index].location, "NODE",
           nodes[node_index].source_id.source_label_text,
           "Least-aligned-axis tangent frame construction failed.");
     }
-    const Vector3 tangent_a = Scale(1.0 / tangent_norm, tangent_candidate);
-    const Vector3 tangent_b = Cross(director, tangent_a);
-    if (!IsFinite(tangent_b) || !(Norm(tangent_b) > 0.0)) {
+    const Vector3 tangent_a = (1.0 / tangent_norm) * tangent_candidate;
+    const Vector3 tangent_b = director.Cross(tangent_a);
+    if (!tangent_b.IsFinite() || !(tangent_b.Norm() > 0.0)) {
       return GeometryFailure(
           "invalid-shell-director", nodes[node_index].location, "NODE",
           nodes[node_index].source_id.source_label_text,
           "Right-handed shell tangent frame construction failed.");
     }
     geometry.nodal_frames.push_back(
-        {static_cast<EntityIndex>(node_index), director, tangent_a, tangent_b});
+        {static_cast<EntityIndex>(node_index), director.Components(),
+         tangent_a.Components(), tangent_b.Components()});
   }
 
   // Build lookup only after frame storage is complete so later code never
@@ -392,7 +362,7 @@ Result<ShellGeometry> PreprocessShellGeometry(
                                "ELEMENT", element.source_id.source_label_text,
                                "Shell element is missing a nodal director.");
       }
-      directors[local_node] = frame->director;
+      directors[local_node] = Vector3{frame->director};
     }
 
     for (const auto& point : ShellGeometryValidationPoints()) {
@@ -404,12 +374,12 @@ Result<ShellGeometry> PreprocessShellGeometry(
           DerivativeSum(shape.xi_derivatives, work[element_index].coordinates);
       const Vector3 midsurface_eta =
           DerivativeSum(shape.eta_derivatives, work[element_index].coordinates);
-      const Vector3 area_vector = Cross(midsurface_xi, midsurface_eta);
-      const double surface_measure = Norm(area_vector);
-      if (!IsFinite(midsurface_xi) || !IsFinite(midsurface_eta) ||
-          !IsFinite(area_vector) || !std::isfinite(surface_measure) ||
+      const Vector3 area_vector = midsurface_xi.Cross(midsurface_eta);
+      const double surface_measure = area_vector.Norm();
+      if (!midsurface_xi.IsFinite() || !midsurface_eta.IsFinite() ||
+          !area_vector.IsFinite() || !std::isfinite(surface_measure) ||
           !(surface_measure > 0.0) ||
-          !(Dot(area_vector, work[element_index].normal) > 0.0)) {
+          !(area_vector.Dot(work[element_index].normal) > 0.0)) {
         return GeometryFailure("invalid-shell-geometry", element.location,
                                "ELEMENT", element.source_id.source_label_text,
                                "Shell surface basis is nonfinite, zero, or "
@@ -421,14 +391,14 @@ Result<ShellGeometry> PreprocessShellGeometry(
           WeightedSum(shape.eta_derivatives, directors);
       const Vector3 director_value = WeightedSum(shape.values, directors);
       const Vector3 covariant_xi =
-          Add(midsurface_xi, Scale(0.5 * thickness * zeta, director_xi));
+          midsurface_xi + 0.5 * thickness * zeta * director_xi;
       const Vector3 covariant_eta =
-          Add(midsurface_eta, Scale(0.5 * thickness * zeta, director_eta));
-      const Vector3 covariant_zeta = Scale(0.5 * thickness, director_value);
+          midsurface_eta + 0.5 * thickness * zeta * director_eta;
+      const Vector3 covariant_zeta = 0.5 * thickness * director_value;
       const double jacobian =
-          Dot(covariant_xi, Cross(covariant_eta, covariant_zeta));
-      if (!IsFinite(covariant_xi) || !IsFinite(covariant_eta) ||
-          !IsFinite(covariant_zeta) || !std::isfinite(jacobian) ||
+          covariant_xi.Dot(covariant_eta.Cross(covariant_zeta));
+      if (!covariant_xi.IsFinite() || !covariant_eta.IsFinite() ||
+          !covariant_zeta.IsFinite() || !std::isfinite(jacobian) ||
           !(jacobian > 0.0)) {
         return GeometryFailure(
             "invalid-shell-jacobian", element.location, "ELEMENT",
@@ -436,13 +406,13 @@ Result<ShellGeometry> PreprocessShellGeometry(
             "Shell Jacobian is nonfinite or nonpositive at a required point.");
       }
       const Vector3 reciprocal_xi =
-          Scale(1.0 / jacobian, Cross(covariant_eta, covariant_zeta));
+          (1.0 / jacobian) * covariant_eta.Cross(covariant_zeta);
       const Vector3 reciprocal_eta =
-          Scale(1.0 / jacobian, Cross(covariant_zeta, covariant_xi));
+          (1.0 / jacobian) * covariant_zeta.Cross(covariant_xi);
       const Vector3 reciprocal_zeta =
-          Scale(1.0 / jacobian, Cross(covariant_xi, covariant_eta));
-      if (!IsFinite(reciprocal_xi) || !IsFinite(reciprocal_eta) ||
-          !IsFinite(reciprocal_zeta)) {
+          (1.0 / jacobian) * covariant_xi.Cross(covariant_eta);
+      if (!reciprocal_xi.IsFinite() || !reciprocal_eta.IsFinite() ||
+          !reciprocal_zeta.IsFinite()) {
         return GeometryFailure(
             "invalid-shell-jacobian", element.location, "ELEMENT",
             element.source_id.source_label_text,
diff --git a/tests/unit/elements/euler_beam_3d_test.cpp b/tests/unit/elements/euler_beam_3d_test.cpp
index 65eb965..d47d8d6 100644
--- a/tests/unit/elements/euler_beam_3d_test.cpp
+++ b/tests/unit/elements/euler_beam_3d_test.cpp
@@ -12,6 +12,8 @@
 #include <utility>
 #include <vector>
 
+#include "fesa/math/vector3.h"
+
 namespace fesa {
 namespace {
 
@@ -720,6 +722,40 @@ TEST(EulerBeam3D, RotatedTransformPreservesWorkAndEnergy) {
               (local_displacement[6U] - local_displacement[0U]) / 3.0, 1.0e-14);
 }
 
+TEST(EulerBeam3D, PreservesExactRotatedResultsAcrossVector3Migration) {
+  const auto section = MakeSection({-2.0, 2.0, 0.0});
+  const auto material = MakeMaterial();
+  const Node first_node = MakeNode({1.0, -2.0, 0.5}, 1U);
+  const Node second_node = MakeNode({3.0, 0.0, 1.5}, 2U);
+  const Vector3 delta =
+      Vector3(second_node.coordinates) - Vector3(first_node.coordinates);
+  EXPECT_DOUBLE_EQ(delta.Norm(), 3.0);
+
+  const auto beam = RequireBeam(first_node, second_node, section, material);
+  const Matrix global = beam.GlobalStiffness();
+  Vector displacement{kElementDofCount};
+  for (std::size_t index = 0U; index < displacement.Size(); ++index) {
+    displacement[index] = 0.01 * static_cast<double>(index + 1U) - 0.04;
+  }
+  const BeamRecovery recovery = beam.Recover(displacement);
+
+  EXPECT_DOUBLE_EQ(global(0U, 0U), 0x1.65f135da12f68p+28);
+  EXPECT_DOUBLE_EQ(global(0U, 1U), 0x1.6261c084bda12p+28);
+  EXPECT_DOUBLE_EQ(global(0U, 2U), 0x1.630ca684bda13p+27);
+  EXPECT_DOUBLE_EQ(global(4U, 4U), 0x1.068e359b59b58p+22);
+  EXPECT_DOUBLE_EQ(global(5U, 11U), 0x1.2d14a7ee7ee7bp+22);
+  EXPECT_DOUBLE_EQ(recovery.gauss_generalized_strains[0U][0U],
+                   0x1.1111111111110p-5);
+  EXPECT_DOUBLE_EQ(recovery.gauss_generalized_strains[0U][1U],
+                   0x1.1111111111111p-5);
+  EXPECT_DOUBLE_EQ(recovery.gauss_generalized_strains[0U][2U],
+                   -0x1.382425a7d2473p-6);
+  EXPECT_DOUBLE_EQ(recovery.gauss_generalized_strains[0U][3U],
+                   -0x1.a9389137b051dp-6);
+  EXPECT_DOUBLE_EQ(recovery.endpoint_section_resultants[1U][2U],
+                   0x1.525c94a87359fp+17);
+}
+
 TEST(EulerBeam3D, ConstantLineLoadMatchesAllSignedComponents) {
   const double length = 4.0;
   const ConstantLocalLineLoad load{2.5, -3.0, 5.5, -7.0};
diff --git a/tests/unit/elements/mitc4_shell_test.cpp b/tests/unit/elements/mitc4_shell_test.cpp
index 05b3955..214500c 100644
--- a/tests/unit/elements/mitc4_shell_test.cpp
+++ b/tests/unit/elements/mitc4_shell_test.cpp
@@ -12,6 +12,8 @@
 #include <string>
 #include <vector>
 
+#include "fesa/math/vector3.h"
+
 namespace {
 
 using Vector3 = std::array<double, 3>;
@@ -588,6 +590,64 @@ TEST(Mitc4ShellKernel,
   EXPECT_DOUBLE_EQ(repeated.drilling_stiffness, stiffness.drilling_stiffness);
 }
 
+TEST(Mitc4ShellKernel,
+     PreservesExactStiffnessRecoveryAndPatchAcrossVector3Migration) {
+  const auto nodes = PlanarNodes();
+  const auto shell_candidate = fesa::Mitc4Shell::Create(
+      NodePointers(nodes), Directors(), Section(), Material());
+  ASSERT_TRUE(shell_candidate.HasValue());
+  const auto& shell = shell_candidate.Value();
+  const auto stiffness_candidate = shell.Stiffness();
+  ASSERT_TRUE(stiffness_candidate.HasValue());
+  const auto& stiffness = stiffness_candidate.Value();
+  const fesa::Vector3 frame_e3{shell.LocalFrame(0.0, 0.0).e3};
+  EXPECT_DOUBLE_EQ(frame_e3.X(), 0.0);
+  EXPECT_DOUBLE_EQ(frame_e3.Y(), 0.0);
+  EXPECT_DOUBLE_EQ(frame_e3.Z(), 1.0);
+
+  constexpr std::array<double, 8> kGeneralized{0.1,   -0.05, 0.2, 0.3,
+                                               -0.15, 0.25,  0.4, -0.3};
+  fesa::Vector global_field{24U};
+  for (std::size_t node_index = 0U; node_index < nodes.size(); ++node_index) {
+    const double x = nodes[node_index].coordinates[0];
+    const double y = nodes[node_index].coordinates[1];
+    const std::size_t offset = 6U * node_index;
+    global_field[offset] = kGeneralized[0] * x + 0.5 * kGeneralized[2] * y;
+    global_field[offset + 1U] = kGeneralized[1] * y + 0.5 * kGeneralized[2] * x;
+    global_field[offset + 2U] = kGeneralized[6] * x + kGeneralized[7] * y -
+                                0.5 * kGeneralized[5] * x * y;
+    global_field[offset + 3U] =
+        -kGeneralized[4] * y - 0.5 * kGeneralized[5] * x;
+    global_field[offset + 4U] = kGeneralized[3] * x + 0.5 * kGeneralized[5] * y;
+  }
+  const auto recovery_candidate = shell.RecoverPhysical(global_field);
+  ASSERT_TRUE(recovery_candidate.HasValue());
+  const auto& recovery = recovery_candidate.Value();
+
+  std::array<std::array<double, 5>, 4> e11_values{};
+  for (std::size_t node_index = 0U; node_index < nodes.size(); ++node_index) {
+    e11_values[node_index][0U] = 0.2 * nodes[node_index].coordinates[0U];
+  }
+  const double gauss = 1.0 / std::sqrt(3.0);
+  const fesa::Vector patch = shell.StrainDisplacement20(gauss, -gauss, gauss)
+                                 .Multiply(PhysicalField(e11_values));
+
+  EXPECT_DOUBLE_EQ(stiffness.physical_local20(0U, 0U), 0x1.d555555555554p+6);
+  EXPECT_DOUBLE_EQ(stiffness.physical_local20(0U, 4U), 0.0);
+  EXPECT_DOUBLE_EQ(stiffness.physical_global24(2U, 2U), 0x1.aaaaaaaaaaaadp+5);
+  EXPECT_DOUBLE_EQ(stiffness.drilling_stiffness, 0x1.0d6cffc5beeb5p-4);
+  EXPECT_DOUBLE_EQ(recovery.strain_energy, 0x1.20a3d70a3d70bp+6);
+  constexpr std::array<double, 8> kExpectedStrain{
+      0x1.999999999999ap-4, -0x1.999999999999bp-5, 0x1.9999999999998p-3,
+      0x1.3333333333335p-2, -0x1.3333333333335p-3, 0x1.0000000000000p-2,
+      0x1.999999999999cp-2, -0x1.3333333333334p-2};
+  EXPECT_EQ(recovery.points[0U].generalized_strain, kExpectedStrain);
+  EXPECT_DOUBLE_EQ(patch[0U], 0x1.9999999999999p-3);
+  for (std::size_t component = 1U; component < patch.Size(); ++component) {
+    EXPECT_DOUBLE_EQ(patch[component], 0.0);
+  }
+}
+
 // MITC4-KERNEL-002
 TEST(Mitc4ShellKernel,
      PreservesPhysicalEnergyUnderTwentyToTwentyFourCongruence) {
diff --git a/tests/unit/math/vector3_test.cpp b/tests/unit/math/vector3_test.cpp
index 17593ed..bd51337 100644
--- a/tests/unit/math/vector3_test.cpp
+++ b/tests/unit/math/vector3_test.cpp
@@ -2,6 +2,8 @@
 
 #include <gtest/gtest.h>
 
+#include <array>
+#include <cmath>
 #include <limits>
 #include <optional>
 
@@ -19,6 +21,7 @@ TEST(Vector3, StoresComponentsWithExactArithmetic) {
   constexpr Vector3 sum = lhs + rhs;
   constexpr Vector3 difference = lhs - rhs;
   constexpr Vector3 scaled = lhs * -2.0;
+  constexpr Vector3 divided = lhs / -2.0;
 
   EXPECT_DOUBLE_EQ(lhs.X(), 1.0);
   EXPECT_DOUBLE_EQ(lhs.Y(), -2.0);
@@ -35,6 +38,22 @@ TEST(Vector3, StoresComponentsWithExactArithmetic) {
   EXPECT_DOUBLE_EQ(scaled[0], -2.0);
   EXPECT_DOUBLE_EQ(scaled[1], 4.0);
   EXPECT_DOUBLE_EQ(scaled[2], -6.0);
+  EXPECT_DOUBLE_EQ(divided[0], -0.5);
+  EXPECT_DOUBLE_EQ(divided[1], 1.0);
+  EXPECT_DOUBLE_EQ(divided[2], -1.5);
+}
+
+TEST(Vector3, ConvertsArrayAndPreservesScalarLeftEvaluationOrder) {
+  constexpr std::array<double, 3> kComponents{1.25, -2.5, 5.0};
+  constexpr Vector3 value{kComponents};
+  constexpr Vector3 scaled = -2.0 * value;
+
+  EXPECT_TRUE(value == Vector3(1.25, -2.5, 5.0));
+  EXPECT_FALSE(value == Vector3(1.25, -2.5, 4.0));
+  EXPECT_EQ(value.Components(), kComponents);
+  EXPECT_DOUBLE_EQ(scaled.X(), -2.5);
+  EXPECT_DOUBLE_EQ(scaled.Y(), 5.0);
+  EXPECT_DOUBLE_EQ(scaled.Z(), -10.0);
 }
 
 TEST(Vector3, ComputesDotProduct) {
@@ -61,6 +80,9 @@ TEST(Vector3, UsesRightHandedCrossProductOrientation) {
 
 TEST(Vector3, ComputesEuclideanNorm) {
   EXPECT_DOUBLE_EQ(Vector3(2.0, -3.0, 6.0).Norm(), 7.0);
+
+  const double expected = std::hypot(1.0e308, 1.0e308, 0.0);
+  EXPECT_DOUBLE_EQ(Vector3(1.0e308, 1.0e308, 0.0).Norm(), expected);
 }
 
 TEST(Vector3, NormalizesNonzeroVectors) {
diff --git a/tests/unit/model/shell_geometry_test.cpp b/tests/unit/model/shell_geometry_test.cpp
index bd563b6..41ab691 100644
--- a/tests/unit/model/shell_geometry_test.cpp
+++ b/tests/unit/model/shell_geometry_test.cpp
@@ -10,6 +10,8 @@
 #include <string>
 #include <vector>
 
+#include "fesa/math/vector3.h"
+
 namespace {
 
 using Vector3 = std::array<double, 3>;
@@ -134,6 +136,34 @@ TEST(Mitc4Geometry,
   }
 }
 
+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{