diff --git a/include/fesa/assembly/surface_pressure_assembler.h b/include/fesa/assembly/surface_pressure_assembler.h
new file mode 100644
index 0000000..7d86721
--- /dev/null
+++ b/include/fesa/assembly/surface_pressure_assembler.h
@@ -0,0 +1,42 @@
+#ifndef FESA_ASSEMBLY_SURFACE_PRESSURE_ASSEMBLER_H_
+#define FESA_ASSEMBLY_SURFACE_PRESSURE_ASSEMBLER_H_
+
+#include <vector>
+
+#include "fesa/analysis/analysis_model.h"
+#include "fesa/assembly/parallel_for.h"
+#include "fesa/core/diagnostic.h"
+#include "fesa/core/status.h"
+#include "fesa/fem/dof_manager.h"
+#include "fesa/loads/shell_surface_pressure.h"
+#include "fesa/math/sparse_matrix.h"
+#include "fesa/math/vector.h"
+
+namespace fesa {
+
+/// @brief Stores one complete full-space surface-pressure assembly candidate.
+struct SurfacePressureAssemblyCandidate {
+  Vector full_external_force;
+  SparseMatrix full_external_tangent;
+  std::vector<Diagnostic> diagnostics;
+};
+
+/// @brief Assembles active surface-pressure records into stable full DOFs.
+/// @param model Non-owning active model whose Domain outlives this call.
+/// @param dofs Owner of the stable shell scatter and full-DOF ordering.
+/// @param full_generalized_coordinates Full `[u1,u2,u3,phi1,phi2,phi3]`
+/// coordinates used only for current-configuration evaluation.
+/// @param load_factor Applied exactly once to every source pressure record.
+/// @param configuration Selects reference or current surface geometry.
+/// @param parallel_for Runs only index-owned local evaluator candidates.
+/// @return A fully validated force/tangent candidate or a structured failure.
+/// @note Global vector and CSR reduction occur serially after all local work.
+Result<SurfacePressureAssemblyCandidate> AssembleSurfacePressure(
+    const AnalysisModel& model, const DofManager& dofs,
+    const Vector& full_generalized_coordinates, double load_factor,
+    SurfacePressureConfiguration configuration,
+    const ParallelFor& parallel_for);
+
+}  // namespace fesa
+
+#endif  // FESA_ASSEMBLY_SURFACE_PRESSURE_ASSEMBLER_H_
diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt
index db367c6..bd64328 100644
--- a/src/fesa/CMakeLists.txt
+++ b/src/fesa/CMakeLists.txt
@@ -15,6 +15,7 @@ add_library(
     assembly/load_assembler.cpp
     assembly/parallel_for.cpp
     assembly/sparse_assembler.cpp
+    assembly/surface_pressure_assembler.cpp
     build_info.cpp
     constraints/essential_constraint_policy.cpp
     constraints/prescribed_displacement.cpp
diff --git a/src/fesa/analysis/linear_static_procedure.cpp b/src/fesa/analysis/linear_static_procedure.cpp
index cbe74df..ed1aafa 100644
--- a/src/fesa/analysis/linear_static_procedure.cpp
+++ b/src/fesa/analysis/linear_static_procedure.cpp
@@ -1,6 +1,8 @@
 #include "fesa/analysis/linear_static_procedure.h"
 
+#include <cmath>
 #include <memory>
+#include <string>
 #include <utility>
 #include <vector>
 
@@ -8,6 +10,7 @@
 #include "fesa/assembly/load_assembler.h"
 #include "fesa/assembly/parallel_for.h"
 #include "fesa/assembly/sparse_assembler.h"
+#include "fesa/assembly/surface_pressure_assembler.h"
 #include "fesa/constraints/essential_constraint_policy.h"
 #include "fesa/elements/element.h"
 #include "fesa/elements/element_factory.h"
@@ -67,12 +70,67 @@ Result<LinearStaticExecutionResult> LinearStaticProcedureExecutor::Execute(
     return Result<LinearStaticExecutionResult>::Failure(factorized);
   }
 
-  Result<Vector> external =
+  Result<Vector> nodal_external =
       LoadAssembler::AssembleFullNodalLoad(model.Value(), dofs);
-  if (!external.HasValue()) {
-    return Result<LinearStaticExecutionResult>::Failure(external.GetStatus());
+  if (!nodal_external.HasValue()) {
+    return Result<LinearStaticExecutionResult>::Failure(
+        nodal_external.GetStatus());
   }
-  state.ExternalForce() = std::move(external.Value());
+  Vector reference_coordinates{dofs.FullDofCount()};
+  Result<SurfacePressureAssemblyCandidate> pressure_external =
+      AssembleSurfacePressure(model.Value(), dofs, reference_coordinates, 1.0,
+                              SurfacePressureConfiguration::kReference,
+                              parallel_for_);
+  if (!pressure_external.HasValue()) {
+    return Result<LinearStaticExecutionResult>::Failure(
+        pressure_external.GetStatus());
+  }
+  if (pressure_external.Value().full_external_force.Size() !=
+          nodal_external.Value().Size() ||
+      pressure_external.Value().full_external_tangent.Rows() !=
+          dofs.FullDofCount() ||
+      pressure_external.Value().full_external_tangent.Columns() !=
+          dofs.FullDofCount()) {
+    return Result<LinearStaticExecutionResult>::Failure(Status::Failure(
+        FailureCategory::kModel,
+        {{Severity::kError,
+          "invalid-pressure-dimensions",
+          {domain.SourcePath(), 0U},
+          "DSLOAD",
+          domain.SourceContentIdentity(),
+          "Pressure assembly must return the stable full-DOF dimensions."}}));
+  }
+  for (const double value :
+       pressure_external.Value().full_external_tangent.Values()) {
+    if (value != 0.0) {
+      return Result<LinearStaticExecutionResult>::Failure(Status::Failure(
+          FailureCategory::kModel,
+          {{Severity::kError,
+            "nonzero-linear-pressure-tangent",
+            {domain.SourcePath(), 0U},
+            "DSLOAD",
+            domain.SourceContentIdentity(),
+            "Reference-configuration linear pressure requires zero "
+            "tangent."}}));
+    }
+  }
+  for (std::size_t index = 0U; index < nodal_external.Value().Size(); ++index) {
+    const double combined =
+        nodal_external.Value()[index] +
+        pressure_external.Value().full_external_force[index];
+    if (!std::isfinite(combined)) {
+      return Result<LinearStaticExecutionResult>::Failure(Status::Failure(
+          FailureCategory::kSolver,
+          {{Severity::kError,
+            "nonfinite-load-accumulation",
+            {domain.SourcePath(), 0U},
+            "DSLOAD",
+            std::to_string(index),
+            "Nodal and surface pressure load accumulation is nonfinite."}}));
+    }
+    nodal_external.Value()[index] = combined;
+  }
+  state.ExternalForce() = std::move(nodal_external.Value());
 
   Result<Vector> rhs = LoadAssembler::EffectiveFreeRhs(
       state.ExternalForce(), partitions.Value().k_fc, dofs.PrescribedValues(),
@@ -97,6 +155,9 @@ Result<LinearStaticExecutionResult> LinearStaticProcedureExecutor::Execute(
   }
 
   std::vector<Diagnostic> diagnostics = domain.Warnings();
+  diagnostics.insert(diagnostics.end(),
+                     pressure_external.Value().diagnostics.begin(),
+                     pressure_external.Value().diagnostics.end());
   SortDiagnostics(diagnostics);
   return Result<LinearStaticExecutionResult>::Success(
       {std::move(state), std::move(diagnostics)});
diff --git a/src/fesa/assembly/surface_pressure_assembler.cpp b/src/fesa/assembly/surface_pressure_assembler.cpp
new file mode 100644
index 0000000..aafe1a6
--- /dev/null
+++ b/src/fesa/assembly/surface_pressure_assembler.cpp
@@ -0,0 +1,291 @@
+#include "fesa/assembly/surface_pressure_assembler.h"
+
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <optional>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "fesa/elements/element_definition.h"
+#include "fesa/model/domain.h"
+
+namespace fesa {
+namespace {
+
+constexpr std::size_t kNodesPerFacet = 4U;
+constexpr std::size_t kDofsPerNode = 6U;
+constexpr std::size_t kFacetDofCount = kNodesPerFacet * kDofsPerNode;
+
+struct PressureTask {
+  ShellSurfacePressureInput input;
+  std::array<std::size_t, kFacetDofCount> scatter;
+};
+
+struct LocalPressureCandidate {
+  ShellSurfacePressureContribution contribution;
+};
+
+std::string Identity(const SurfacePressureLoadDefinition& load,
+                     const ElementSurfaceFacetDefinition* const facet) {
+  std::string identity =
+      "surface_raw=" + load.surface_name_raw + "; surface_canonical=" +
+      load.surface_name_canonical + "; load_type=" + load.load_type +
+      "; load_ordinal=" + std::to_string(load.load_ordinal);
+  if (facet != nullptr) {
+    identity += "; part_raw=" + facet->part_name_raw +
+                "; instance_raw=" + facet->instance_name_raw +
+                "; element_label=" + facet->source_element_label_raw +
+                "; source_element_type=" + facet->source_element_type +
+                "; side=" +
+                (facet->side == ShellSurfaceSide::kSpos ? "SPOS" : "SNEG") +
+                "; member_ordinal=" +
+                std::to_string(facet->first_member_ordinal);
+  }
+  return identity;
+}
+
+Status Failure(const FailureCategory category, const SourceLocation& location,
+               const std::string& code, const std::string& identity,
+               const std::string& message) {
+  return Status::Failure(category,
+                         {{Severity::kError, code, location, "DSLOAD",
+                           identity, message}});
+}
+
+Status ValidateFullCoordinates(const Domain& domain, const DofManager& dofs,
+                               const Vector& coordinates,
+                               const double load_factor) {
+  if (!std::isfinite(load_factor)) {
+    return Failure(FailureCategory::kSolver, {domain.SourcePath(), 0U},
+                   "nonfinite-pressure-contribution", "load_factor",
+                   "The pressure load factor must be finite.");
+  }
+  if (domain.Nodes().size() >
+      (std::numeric_limits<std::size_t>::max)() / kDofsPerNode) {
+    return Failure(FailureCategory::kModel, {domain.SourcePath(), 0U},
+                   "invalid-pressure-dimensions",
+                   domain.SourceContentIdentity(),
+                   "The semantic node count cannot form a full DOF vector.");
+  }
+  const std::size_t expected_count = domain.Nodes().size() * kDofsPerNode;
+  if (dofs.FullDofCount() != expected_count ||
+      coordinates.Size() != expected_count) {
+    return Failure(FailureCategory::kModel, {domain.SourcePath(), 0U},
+                   "invalid-pressure-dimensions",
+                   domain.SourceContentIdentity(),
+                   "Pressure assembly requires the stable full-DOF dimension.");
+  }
+  const Status dof_status = dofs.ValidateInvariants();
+  if (!dof_status.IsOk()) {
+    return dof_status;
+  }
+  for (std::size_t index = 0U; index < coordinates.Size(); ++index) {
+    if (!std::isfinite(coordinates[index])) {
+      return Failure(FailureCategory::kSolver, {domain.SourcePath(), 0U},
+                     "nonfinite-pressure-contribution", std::to_string(index),
+                     "The full generalized coordinates must be finite.");
+    }
+  }
+  return Status::Ok();
+}
+
+Result<PressureTask> CreateTask(
+    const Domain& domain, const DofManager& dofs,
+    const SurfacePressureLoadDefinition& load,
+    const ElementSurfaceFacetDefinition& facet,
+    const Vector& full_generalized_coordinates, const double load_factor,
+    const SurfacePressureConfiguration configuration) {
+  const std::string identity = Identity(load, &facet);
+  if (facet.element_index >= domain.Elements().Size()) {
+    return Result<PressureTask>::Failure(Failure(
+        FailureCategory::kModel, load.data_location, "invalid-pressure-facet",
+        identity, "The pressure facet element is outside the Domain."));
+  }
+  const ElementDefinition& element = domain.Elements()[facet.element_index];
+  if (element.Kind() != ElementDefinitionKind::kMitc4Shell ||
+      (facet.source_element_type != "S4" &&
+       facet.source_element_type != "S4R") ||
+      element.SourceElementType() != facet.source_element_type ||
+      element.NodeIndices().size() != kNodesPerFacet) {
+    return Result<PressureTask>::Failure(Failure(
+        FailureCategory::kModel, load.data_location, "invalid-pressure-facet",
+        identity, "Pressure requires one matching four-node S4 or S4R facet."));
+  }
+
+  PressureTask task{};
+  task.input.side = facet.side;
+  task.input.pressure_magnitude = load.magnitude;
+  task.input.load_factor = load_factor;
+  task.input.configuration = configuration;
+  task.input.location = load.data_location;
+  task.input.entity_identity = identity;
+  try {
+    task.scatter = dofs.ShellElementScatter(facet.element_index);
+  } catch (const std::out_of_range&) {
+    return Result<PressureTask>::Failure(
+        Failure(FailureCategory::kModel, load.data_location,
+                "invalid-pressure-scatter", identity,
+                "The pressure facet does not have a 24-DOF shell scatter."));
+  }
+
+  for (std::size_t local_node = 0U; local_node < kNodesPerFacet;
+       ++local_node) {
+    const EntityIndex node_index = element.NodeIndices()[local_node];
+    if (node_index >= domain.Nodes().size()) {
+      return Result<PressureTask>::Failure(Failure(
+          FailureCategory::kModel, load.data_location, "invalid-pressure-facet",
+          identity, "The pressure facet node is outside the Domain."));
+    }
+    task.input.reference_coordinates[local_node] =
+        Vector3{domain.Nodes()[node_index].coordinates};
+    std::array<double, 3> displacement{};
+    for (std::size_t component = 0U; component < displacement.size();
+         ++component) {
+      const std::size_t full_dof =
+          task.scatter[local_node * kDofsPerNode + component];
+      if (full_dof >= full_generalized_coordinates.Size()) {
+        return Result<PressureTask>::Failure(
+            Failure(FailureCategory::kModel, load.data_location,
+                    "invalid-pressure-scatter", identity,
+                    "The pressure facet scatter is outside the full vector."));
+      }
+      displacement[component] = full_generalized_coordinates[full_dof];
+    }
+    task.input.translational_displacements[local_node] = Vector3{displacement};
+  }
+  return Result<PressureTask>::Success(std::move(task));
+}
+
+}  // namespace
+
+Result<SurfacePressureAssemblyCandidate> AssembleSurfacePressure(
+    const AnalysisModel& model, const DofManager& dofs,
+    const Vector& full_generalized_coordinates, const double load_factor,
+    const SurfacePressureConfiguration configuration,
+    const ParallelFor& parallel_for) {
+  const Domain& domain = model.GetDomain();
+  const Status coordinate_status = ValidateFullCoordinates(
+      domain, dofs, full_generalized_coordinates, load_factor);
+  if (!coordinate_status.IsOk()) {
+    return Result<SurfacePressureAssemblyCandidate>::Failure(coordinate_status);
+  }
+
+  const auto& loads = model.ActiveSurfacePressureLoads();
+  const auto& surfaces = domain.ElementSurfaces();
+  std::vector<PressureTask> tasks;
+  for (std::size_t load_index = 0U; load_index < loads.size(); ++load_index) {
+    const auto& load = loads[load_index];
+    if (load.load_ordinal != load_index ||
+        load.surface_index >= surfaces.size()) {
+      return Result<SurfacePressureAssemblyCandidate>::Failure(Failure(
+          FailureCategory::kModel, load.data_location, "invalid-pressure-order",
+          Identity(load, nullptr),
+          "Pressure records must retain complete stable source order."));
+    }
+    const auto& surface = surfaces[load.surface_index];
+    if (load.load_type != "P" ||
+        load.surface_name_canonical != surface.surface_name_canonical) {
+      return Result<SurfacePressureAssemblyCandidate>::Failure(Failure(
+          FailureCategory::kModel, load.data_location,
+          "invalid-pressure-target",
+          Identity(load, nullptr),
+          "Pressure record type and resolved surface identity must agree."));
+    }
+    if (surface.facets.empty()) {
+      return Result<SurfacePressureAssemblyCandidate>::Failure(Failure(
+          FailureCategory::kModel, load.data_location, "empty-surface",
+          Identity(load, nullptr),
+          "Pressure assembly requires a nonempty resolved surface."));
+    }
+    for (std::size_t facet_index = 0U; facet_index < surface.facets.size();
+         ++facet_index) {
+      const auto& facet = surface.facets[facet_index];
+      if (facet.facet_ordinal != facet_index) {
+        return Result<SurfacePressureAssemblyCandidate>::Failure(Failure(
+            FailureCategory::kModel, load.data_location,
+            "invalid-pressure-order",
+            Identity(load, &facet),
+            "Pressure facets must retain their stable resolved order."));
+      }
+      auto task = CreateTask(domain, dofs, load, facet,
+                             full_generalized_coordinates, load_factor,
+                             configuration);
+      if (!task.HasValue()) {
+        return Result<SurfacePressureAssemblyCandidate>::Failure(
+            task.GetStatus());
+      }
+      tasks.push_back(std::move(task.Value()));
+    }
+  }
+
+  std::vector<LocalPressureCandidate> local_candidates(tasks.size());
+  std::vector<std::optional<Status>> local_failures(tasks.size());
+  parallel_for.Execute(tasks.size(),
+                       [&tasks, &local_candidates,
+                        &local_failures](const std::size_t index) {
+    const auto result = EvaluateShellSurfacePressure(tasks[index].input);
+    if (!result.HasValue()) {
+      local_failures[index] = result.GetStatus();
+      return;
+    }
+    local_candidates[index].contribution = result.Value();
+  });
+  for (std::size_t index = 0U; index < tasks.size(); ++index) {
+    if (local_failures[index].has_value()) {
+      return Result<SurfacePressureAssemblyCandidate>::Failure(
+          std::move(*local_failures[index]));
+    }
+  }
+
+  Vector full_force{dofs.FullDofCount()};
+  std::vector<CooContribution> tangent_contributions;
+  tangent_contributions.reserve(tasks.size() * kFacetDofCount * kFacetDofCount);
+  for (std::size_t task_index = 0U; task_index < tasks.size(); ++task_index) {
+    const auto& task = tasks[task_index];
+    const auto& contribution = local_candidates[task_index].contribution;
+    for (std::size_t local_row = 0U; local_row < kFacetDofCount; ++local_row) {
+      const std::size_t full_row = task.scatter[local_row];
+      const double force = contribution.force[local_row];
+      const double accumulated = full_force[full_row] + force;
+      if (!std::isfinite(force) || !std::isfinite(accumulated)) {
+        return Result<SurfacePressureAssemblyCandidate>::Failure(Failure(
+            FailureCategory::kSolver, task.input.location,
+            "nonfinite-pressure-accumulation", task.input.entity_identity,
+            "Stable pressure force reduction produced a nonfinite value."));
+      }
+      full_force[full_row] = accumulated;
+      for (std::size_t local_column = 0U; local_column < kFacetDofCount;
+           ++local_column) {
+        const double tangent = contribution.external_tangent[
+            local_row * kFacetDofCount + local_column];
+        if (!std::isfinite(tangent)) {
+          return Result<SurfacePressureAssemblyCandidate>::Failure(Failure(
+              FailureCategory::kSolver, task.input.location,
+              "nonfinite-pressure-tangent", task.input.entity_identity,
+              "A pressure local tangent contains a nonfinite value."));
+        }
+        tangent_contributions.push_back(
+            {full_row,
+             task.scatter[local_column],
+             tangent,
+             task_index,
+             local_row * kFacetDofCount + local_column});
+      }
+    }
+  }
+  auto tangent = SparseMatrix::FromCoo(
+      dofs.FullDofCount(), dofs.FullDofCount(),
+      std::move(tangent_contributions), dofs.GetSparsePattern());
+  if (!tangent.HasValue()) {
+    return Result<SurfacePressureAssemblyCandidate>::Failure(
+        tangent.GetStatus());
+  }
+  return Result<SurfacePressureAssemblyCandidate>::Success(
+      {std::move(full_force), std::move(tangent.Value()), {}});
+}
+
+}  // namespace fesa
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 22685cf..fe178c3 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -11,6 +11,7 @@ add_executable(
     unit/analysis/newton_iteration_log_test.cpp
     unit/assembly/parallel_for_test.cpp
     unit/assembly/load_assembler_test.cpp
+    unit/assembly/surface_pressure_assembler_test.cpp
     unit/assembly/sparse_assembler_test.cpp
     unit/app/output_paths_test.cpp
     unit/constraints/boundary_condition_test.cpp
@@ -74,8 +75,10 @@ target_link_libraries(
 add_executable(
     fesa_integration_tests
     integration/analysis/linear_static_analysis_test.cpp
+    integration/analysis/shell_surface_pressure_linear_test.cpp
     integration/analysis/nonlinear_static_procedure_test.cpp
     integration/app/fesa_application_test.cpp
+    integration/app/shell_surface_pressure_application_test.cpp
 )
 
 target_link_libraries(
diff --git a/tests/integration/analysis/shell_surface_pressure_linear_test.cpp b/tests/integration/analysis/shell_surface_pressure_linear_test.cpp
new file mode 100644
index 0000000..8912b87
--- /dev/null
+++ b/tests/integration/analysis/shell_surface_pressure_linear_test.cpp
@@ -0,0 +1,272 @@
+#include "fesa/analysis/linear_static_analysis.h"
+
+#include <gtest/gtest.h>
+
+#include <array>
+#include <atomic>
+#include <chrono>
+#include <cstddef>
+#include <cstdint>
+#include <filesystem>
+#include <fstream>
+#include <memory>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <hdf5.h>
+
+#include "fesa/analysis/analysis.h"
+#include "fesa/app/fesa_application.h"
+#include "fesa/assembly/parallel_for.h"
+#include "fesa/results/results_writer.h"
+#include "fesa/solvers/linear/mkl_pardiso_solver.h"
+
+namespace {
+
+class TempDirectory {
+ public:
+  explicit TempDirectory(const std::string& label) {
+    static std::atomic<unsigned long long> sequence{0U};
+    const auto tick =
+        std::chrono::steady_clock::now().time_since_epoch().count();
+    path_ = std::filesystem::temp_directory_path() /
+            ("fesa-surface-pressure-linear-" + label + "-" +
+             std::to_string(tick) + "-" +
+             std::to_string(sequence.fetch_add(1U)));
+    std::error_code error;
+    if (!std::filesystem::create_directory(path_, error) || error) {
+      throw std::runtime_error{"Unable to create surface pressure fixture."};
+    }
+  }
+
+  TempDirectory(const TempDirectory&) = delete;
+  TempDirectory& operator=(const TempDirectory&) = delete;
+
+  ~TempDirectory() {
+    std::error_code ignored;
+    std::filesystem::remove_all(path_, ignored);
+  }
+
+  const std::filesystem::path& Path() const noexcept { return path_; }
+
+ private:
+  std::filesystem::path path_;
+};
+
+void WriteText(const std::filesystem::path& path, const std::string& text) {
+  std::ofstream stream{path, std::ios::binary | std::ios::trunc};
+  stream.write(text.data(), static_cast<std::streamsize>(text.size()));
+  if (!stream) {
+    throw std::runtime_error{"Unable to write surface pressure input."};
+  }
+}
+
+std::string ShellDeck(const std::string& source_type,
+                      const std::string& connectivity,
+                      const bool prescribed_displacement,
+                      const bool with_pressure) {
+  std::string boundary;
+  if (prescribed_displacement) {
+    boundary = R"inp(*Boundary
+N1, 1, 6
+N2, 1, 2
+N2, 3, 3, 5.0
+N2, 4, 6
+N3, 1, 6
+N4, 1, 6
+)inp";
+  } else {
+    boundary = "*Boundary\nAll, 1, 6\n";
+  }
+  const std::string surface = with_pressure
+                                  ? "*Surface, TYPE=ELEMENT, NAME=Pressure\n"
+                                    "Shell-1.1, SPOS\n"
+                                  : "";
+  const std::string dsload = with_pressure
+                                 ? "*Dsload\nPressure, P, 4.0\n"
+                                 : "";
+  return std::string{R"inp(*Part, name=ShellPart
+*Node
+1, 0., 0., 0.
+2, 1., 0., 0.
+3, 1., 1., 0.
+4, 0., 1., 0.
+*Element, type=)inp"} +
+         source_type + "\n1, " + connectivity + R"inp(
+*Elset, elset=ShellSet
+1
+*Shell Section, elset=ShellSet, material=Steel
+0.1
+*End Part
+*Assembly, name=Assembly
+*Instance, name=Shell-1, part=ShellPart
+*End Instance
+)inp" +
+         surface + R"inp(*Nset, nset=N1, instance=Shell-1
+1
+*Nset, nset=N2, instance=Shell-1
+2
+*Nset, nset=N3, instance=Shell-1
+3
+*Nset, nset=N4, instance=Shell-1
+4
+*Nset, nset=All, instance=Shell-1
+1, 2, 3, 4
+*End Assembly
+*Material, name=Steel
+*Elastic
+1000., 0.25
+)inp" +
+         boundary + R"inp(*Step, name=Load, nlgeom=NO
+*Static
+0.1, 1., 0.01, 1.
+)inp" +
+         dsload + "*End Step\n";
+}
+
+class CapturingResultsWriter final : public fesa::ResultsWriter {
+ public:
+  fesa::Status Write(const std::filesystem::path&, const fesa::Domain&,
+                     const fesa::AnalysisState& state,
+                     const std::vector<fesa::Diagnostic>&) override {
+    state_ = std::make_unique<fesa::AnalysisState>(state);
+    return fesa::Status::Ok();
+  }
+
+  const fesa::AnalysisState& State() const {
+    if (!state_) {
+      throw std::logic_error{"No linear state was captured."};
+    }
+    return *state_;
+  }
+
+ private:
+  std::unique_ptr<fesa::AnalysisState> state_;
+};
+
+fesa::AnalysisState RunLinear(const std::filesystem::path& input,
+                              const std::filesystem::path& output) {
+  fesa::SerialParallelFor serial;
+  fesa::MklPardisoSolver solver;
+  CapturingResultsWriter writer;
+  std::unique_ptr<fesa::Analysis> analysis =
+      std::make_unique<fesa::LinearStaticAnalysis>(serial, solver, writer);
+  const fesa::Status status = analysis->Run({input, output});
+  if (!status.IsOk()) {
+    throw std::runtime_error{status.Diagnostics().empty()
+                                 ? "Linear pressure analysis failed."
+                                 : status.Diagnostics().front().code};
+  }
+  return writer.State();
+}
+
+std::array<double, 3> SumTranslation(const fesa::Vector& values) {
+  std::array<double, 3> sum{};
+  for (std::size_t node = 0U; node < values.Size() / 6U; ++node) {
+    for (std::size_t component = 0U; component < 3U; ++component) {
+      sum[component] += values[node * 6U + component];
+    }
+  }
+  return sum;
+}
+
+}  // namespace
+
+TEST(LinearShellSurfacePressure,
+     UsesReferenceNormalAndAreaAfterPrescribedDisplacement) {
+  TempDirectory directory{"reference-configuration"};
+  const auto input = directory.Path() / "prescribed.inp";
+  WriteText(input, ShellDeck("S4", "1, 2, 3, 4", true, true));
+
+  const fesa::AnalysisState state =
+      RunLinear(input, directory.Path() / "captured.h5");
+
+  EXPECT_DOUBLE_EQ(state.Displacement()[8U], 5.0);
+  EXPECT_EQ(SumTranslation(state.ExternalForce()),
+            (std::array<double, 3>{0.0, 0.0, -4.0}));
+}
+
+TEST(LinearShellSurfacePressure,
+     AssemblesBeforeProjectionAndIncludesConstrainedReaction) {
+  TempDirectory directory{"reaction"};
+  const auto input = directory.Path() / "reaction.inp";
+  WriteText(input, ShellDeck("S4", "1, 2, 3, 4", false, true));
+
+  const fesa::AnalysisState state =
+      RunLinear(input, directory.Path() / "captured.h5");
+
+  EXPECT_EQ(SumTranslation(state.ExternalForce()),
+            (std::array<double, 3>{0.0, 0.0, -4.0}));
+  EXPECT_EQ(SumTranslation(state.Reaction()),
+            (std::array<double, 3>{0.0, 0.0, 4.0}));
+}
+
+TEST(LinearShellSurfacePressure,
+     SatisfiesFreeResidualGlobalResultantAndMoment) {
+  TempDirectory directory{"equilibrium"};
+  const auto input = directory.Path() / "equilibrium.inp";
+  WriteText(input, ShellDeck("S4", "1, 2, 3, 4", false, true));
+
+  const fesa::AnalysisState state =
+      RunLinear(input, directory.Path() / "captured.h5");
+
+  for (const double value : state.Equilibrium()) {
+    EXPECT_NEAR(value, 0.0, 1.0e-8);
+  }
+  EXPECT_NEAR(state.VerificationMetrics()[0U], 0.0, 1.0e-10);
+}
+
+TEST(LinearShellSurfacePressure, PreservesS4S4rMeshAndCyclicOrderInvariance) {
+  TempDirectory directory{"source-invariance"};
+  const auto s4_input = directory.Path() / "s4.inp";
+  const auto s4r_input = directory.Path() / "s4r.inp";
+  const auto cyclic_input = directory.Path() / "cyclic.inp";
+  WriteText(s4_input, ShellDeck("S4", "1, 2, 3, 4", false, true));
+  WriteText(s4r_input, ShellDeck("S4R", "1, 2, 3, 4", false, true));
+  WriteText(cyclic_input, ShellDeck("S4", "2, 3, 4, 1", false, true));
+
+  const fesa::AnalysisState s4 =
+      RunLinear(s4_input, directory.Path() / "s4.h5");
+  const fesa::AnalysisState s4r =
+      RunLinear(s4r_input, directory.Path() / "s4r.h5");
+  const fesa::AnalysisState cyclic =
+      RunLinear(cyclic_input, directory.Path() / "cyclic.h5");
+
+  ASSERT_EQ(s4.ExternalForce().Size(), s4r.ExternalForce().Size());
+  ASSERT_EQ(s4.ExternalForce().Size(), cyclic.ExternalForce().Size());
+  for (std::size_t index = 0U; index < s4.ExternalForce().Size(); ++index) {
+    EXPECT_DOUBLE_EQ(s4.ExternalForce()[index], s4r.ExternalForce()[index]);
+    EXPECT_DOUBLE_EQ(s4.ExternalForce()[index], cyclic.ExternalForce()[index]);
+    EXPECT_DOUBLE_EQ(s4.Reaction()[index], s4r.Reaction()[index]);
+    EXPECT_DOUBLE_EQ(s4.Reaction()[index], cyclic.Reaction()[index]);
+  }
+}
+
+TEST(LinearShellSurfacePressure, WritesSchemaZeroWithoutPressureDataset) {
+  TempDirectory directory{"schema-zero"};
+  const auto input = directory.Path() / "pressure.inp";
+  const auto output = directory.Path() / "pressure_results.h5";
+  WriteText(input, ShellDeck("S4", "1, 2, 3, 4", false, true));
+
+  EXPECT_EQ(fesa::FesaApplication{}.Run({input.string()}), 0);
+  ASSERT_GT(H5Fis_hdf5(output.string().c_str()), 0);
+  const hid_t file =
+      H5Fopen(output.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+  ASSERT_GE(file, 0);
+  const hid_t metadata = H5Gopen2(file, "/metadata", H5P_DEFAULT);
+  ASSERT_GE(metadata, 0);
+  const hid_t schema = H5Aopen(metadata, "schema_version", H5P_DEFAULT);
+  ASSERT_GE(schema, 0);
+  std::uint64_t schema_version = 0U;
+  ASSERT_GE(H5Aread(schema, H5T_NATIVE_UINT64, &schema_version), 0);
+  EXPECT_EQ(schema_version, 0U);
+  EXPECT_EQ(H5Lexists(file, "/steps/Step-1/frames/0/surface_pressure",
+                       H5P_DEFAULT),
+            0);
+  EXPECT_EQ(H5Lexists(file, "/pressure", H5P_DEFAULT), 0);
+  EXPECT_GE(H5Aclose(schema), 0);
+  EXPECT_GE(H5Gclose(metadata), 0);
+  EXPECT_GE(H5Fclose(file), 0);
+}
diff --git a/tests/integration/app/shell_surface_pressure_application_test.cpp b/tests/integration/app/shell_surface_pressure_application_test.cpp
new file mode 100644
index 0000000..e270098
--- /dev/null
+++ b/tests/integration/app/shell_surface_pressure_application_test.cpp
@@ -0,0 +1,145 @@
+#include "fesa/app/fesa_application.h"
+
+#include <gtest/gtest.h>
+
+#include <atomic>
+#include <chrono>
+#include <filesystem>
+#include <fstream>
+#include <iterator>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+#include <hdf5.h>
+
+namespace {
+
+class TempDirectory {
+ public:
+  explicit TempDirectory(const std::string& label) {
+    static std::atomic<unsigned long long> sequence{0U};
+    const auto tick =
+        std::chrono::steady_clock::now().time_since_epoch().count();
+    path_ = std::filesystem::temp_directory_path() /
+            ("fesa-surface-pressure-app-" + label + "-" +
+             std::to_string(tick) + "-" +
+             std::to_string(sequence.fetch_add(1U)));
+    std::error_code error;
+    if (!std::filesystem::create_directory(path_, error) || error) {
+      throw std::runtime_error{"Unable to create application fixture."};
+    }
+  }
+
+  TempDirectory(const TempDirectory&) = delete;
+  TempDirectory& operator=(const TempDirectory&) = delete;
+
+  ~TempDirectory() {
+    std::error_code ignored;
+    std::filesystem::remove_all(path_, ignored);
+  }
+
+  const std::filesystem::path& Path() const noexcept { return path_; }
+
+ private:
+  std::filesystem::path path_;
+};
+
+void WriteText(const std::filesystem::path& path, const std::string& text) {
+  std::ofstream stream{path, std::ios::binary | std::ios::trunc};
+  stream.write(text.data(), static_cast<std::streamsize>(text.size()));
+  if (!stream) {
+    throw std::runtime_error{"Unable to write application input."};
+  }
+}
+
+std::vector<char> ReadBytes(const std::filesystem::path& path) {
+  std::ifstream stream{path, std::ios::binary};
+  return {std::istreambuf_iterator<char>{stream},
+          std::istreambuf_iterator<char>{}};
+}
+
+std::string ShellDeck(const bool with_pressure, const bool collapsed) {
+  const std::string third_node = collapsed ? "3, 1., 0., 0.\n"
+                                           : "3, 1., 1., 0.\n";
+  const std::string surface = with_pressure
+                                  ? "*Surface, TYPE=ELEMENT, NAME=Pressure\n"
+                                    "Shell-1.1, SPOS\n"
+                                  : "";
+  const std::string dsload = with_pressure
+                                 ? "*Dsload\nPressure, P, 4.0\n"
+                                 : "";
+  return std::string{R"inp(*Part, name=ShellPart
+*Node
+1, 0., 0., 0.
+2, 1., 0., 0.
+)inp"} +
+         third_node + R"inp(4, 0., 1., 0.
+*Element, type=S4
+1, 1, 2, 3, 4
+*Elset, elset=ShellSet
+1
+*Shell Section, elset=ShellSet, material=Steel
+0.1
+*End Part
+*Assembly, name=Assembly
+*Instance, name=Shell-1, part=ShellPart
+*End Instance
+)inp" +
+         surface + R"inp(*Nset, nset=All, instance=Shell-1
+1, 2, 3, 4
+*End Assembly
+*Material, name=Steel
+*Elastic
+1000., 0.25
+*Boundary
+All, 1, 6
+*Step, name=Load, nlgeom=NO
+*Static
+0.1, 1., 0.01, 1.
+)inp" +
+         dsload + "*End Step\n";
+}
+
+}  // namespace
+
+TEST(ShellSurfacePressureApplication,
+     LinearFailurePreservesExistingFinalAndNoCandidate) {
+  TempDirectory directory{"preserve-final"};
+  const auto input = directory.Path() / "invalid-pressure.inp";
+  const auto output = directory.Path() / "invalid-pressure_results.h5";
+  WriteText(input, ShellDeck(true, true));
+  const std::vector<char> previous{'p', 'r', 'i', 'o', 'r'};
+  {
+    std::ofstream stream{output, std::ios::binary | std::ios::trunc};
+    stream.write(previous.data(),
+                 static_cast<std::streamsize>(previous.size()));
+  }
+
+  const int exit_code = fesa::FesaApplication{}.Run({input.string()});
+
+  EXPECT_NE(exit_code, 0);
+  EXPECT_EQ(ReadBytes(output), previous);
+}
+
+TEST(ShellSurfacePressureApplication,
+     InputWithoutPressureRetainsExistingLinearBehavior) {
+  TempDirectory directory{"no-pressure"};
+  const auto input = directory.Path() / "no-pressure.inp";
+  const auto output = directory.Path() / "no-pressure_results.h5";
+  WriteText(input, ShellDeck(false, false));
+
+  EXPECT_EQ(fesa::FesaApplication{}.Run({input.string()}), 0);
+  ASSERT_TRUE(std::filesystem::exists(output));
+  ASSERT_GT(H5Fis_hdf5(output.string().c_str()), 0);
+  const hid_t file =
+      H5Fopen(output.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+  ASSERT_GE(file, 0);
+  EXPECT_GT(H5Lexists(file, "/steps/Step-1/frames/0/nodal/displacement",
+                       H5P_DEFAULT),
+            0);
+  EXPECT_EQ(H5Lexists(file, "/steps/Step-1/frames/0/surface_pressure",
+                       H5P_DEFAULT),
+            0);
+  EXPECT_GE(H5Fclose(file), 0);
+}
diff --git a/tests/unit/assembly/surface_pressure_assembler_test.cpp b/tests/unit/assembly/surface_pressure_assembler_test.cpp
new file mode 100644
index 0000000..b9515e6
--- /dev/null
+++ b/tests/unit/assembly/surface_pressure_assembler_test.cpp
@@ -0,0 +1,195 @@
+#include "fesa/assembly/surface_pressure_assembler.h"
+
+#include <gtest/gtest.h>
+
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <filesystem>
+#include <memory>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "fesa/analysis/analysis_model.h"
+#include "fesa/assembly/parallel_for.h"
+#include "fesa/fem/dof_manager.h"
+#include "fesa/model/domain.h"
+
+namespace {
+
+struct PressureFixture {
+  std::unique_ptr<fesa::Domain> domain;
+  std::unique_ptr<fesa::AnalysisModel> model;
+  std::unique_ptr<fesa::DofManager> dofs;
+};
+
+PressureFixture MakeFixture(const std::vector<double>& magnitudes) {
+  const std::filesystem::path source{"surface-pressure-assembly.inp"};
+  fesa::ModelDefinition definition{};
+  definition.source_path = source;
+  definition.source_content_identity = "fnv1a64:surfacepressure";
+  definition.nodes = {
+      {{"Shell-1", 1, "1"}, {0.0, 0.0, 0.0}, {source, 2U}},
+      {{"Shell-1", 2, "2"}, {1.0, 0.0, 0.0}, {source, 3U}},
+      {{"Shell-1", 3, "3"}, {1.0, 1.0, 0.0}, {source, 4U}},
+      {{"Shell-1", 4, "4"}, {0.0, 1.0, 0.0}, {source, 5U}}};
+  definition.materials = {{"Steel", 1000.0, 0.25, {source, 6U}}};
+  definition.shell_sections = {{"Shell", 0.1, 0U, {source, 7U}}};
+  definition.shell_elements = {{{"Shell-1", 1, "1"},
+                                fesa::ShellSourceElementType::kS4,
+                                {0U, 1U, 2U, 3U},
+                                0U,
+                                0U,
+                                {source, 8U}}};
+  for (std::size_t node = 0U; node < definition.nodes.size(); ++node) {
+    definition.shell_node_initial_frames.push_back(
+        {static_cast<fesa::EntityIndex>(node),
+         {0.0, 0.0, 1.0},
+         {1.0, 0.0, 0.0},
+         {0.0, 1.0, 0.0}});
+  }
+  definition.element_surfaces = {
+      {"Pressure-Surface",
+       "PRESSURE-SURFACE",
+       {{0U,
+         0U,
+         "ShellPart",
+         "SHELLPART",
+         "Shell-1",
+         "SHELL-1",
+         1,
+         "1",
+         "S4",
+         fesa::ShellSurfaceSide::kSpos,
+         {source, 12U},
+         0U,
+         fesa::SurfaceMemberKind::kDirect,
+         "",
+         ""}},
+       {source, 11U},
+       0U}};
+
+  std::vector<fesa::SurfacePressureLoadDefinition> pressure_loads;
+  pressure_loads.reserve(magnitudes.size());
+  for (std::size_t ordinal = 0U; ordinal < magnitudes.size(); ++ordinal) {
+    pressure_loads.push_back({"Step-1",
+                              "STEP-1",
+                              0U,
+                              "Pressure-Surface",
+                              "PRESSURE-SURFACE",
+                              "P",
+                              magnitudes[ordinal],
+                              std::to_string(magnitudes[ordinal]),
+                              {source, 20U + ordinal},
+                              {source, 21U + ordinal},
+                              ordinal,
+                              ordinal});
+  }
+  definition.steps.emplace_back(
+      "Step-1", std::vector<fesa::PrescribedDisplacementDefinition>{},
+      std::vector<fesa::NodalLoad>{}, std::move(pressure_loads),
+      fesa::StaticProcedureKind::kLinearSmallDisplacement, false, 0.1, 1.0,
+      0.01, 1.0, fesa::SourceLocation{source, 19U});
+
+  auto domain_result = fesa::Domain::Create(std::move(definition));
+  if (!domain_result.HasValue()) {
+    throw std::runtime_error{"Unable to construct pressure Domain fixture."};
+  }
+  auto domain =
+      std::make_unique<fesa::Domain>(std::move(domain_result.Value()));
+  auto model_result = fesa::AnalysisModel::Create(*domain);
+  if (!model_result.HasValue()) {
+    throw std::runtime_error{"Unable to construct pressure AnalysisModel."};
+  }
+  auto model =
+      std::make_unique<fesa::AnalysisModel>(std::move(model_result.Value()));
+  auto dofs_result = fesa::DofManager::Create(*model);
+  if (!dofs_result.HasValue()) {
+    throw std::runtime_error{"Unable to construct pressure DofManager."};
+  }
+  auto dofs =
+      std::make_unique<fesa::DofManager>(std::move(dofs_result.Value()));
+  return {std::move(domain), std::move(model), std::move(dofs)};
+}
+
+std::vector<double> Values(const fesa::Vector& vector) {
+  return {vector.Data(), vector.Data() + vector.Size()};
+}
+
+}  // namespace
+
+TEST(SurfacePressureAssembly, AssemblesFullDofsInPressureFacetAndRecordOrder) {
+  auto fixture = MakeFixture({2.0});
+  fesa::Vector coordinates{fixture.dofs->FullDofCount()};
+  coordinates[2U] = 11.0;
+  fesa::SerialParallelFor serial;
+
+  const auto result = fesa::AssembleSurfacePressure(
+      *fixture.model, *fixture.dofs, coordinates, 1.0,
+      fesa::SurfacePressureConfiguration::kReference, serial);
+
+  ASSERT_TRUE(result.HasValue());
+  EXPECT_TRUE(result.Value().diagnostics.empty());
+  ASSERT_EQ(result.Value().full_external_force.Size(), 24U);
+  for (std::size_t node = 0U; node < 4U; ++node) {
+    const std::size_t offset = node * 6U;
+    EXPECT_DOUBLE_EQ(result.Value().full_external_force[offset], 0.0);
+    EXPECT_DOUBLE_EQ(result.Value().full_external_force[offset + 1U], 0.0);
+    EXPECT_DOUBLE_EQ(result.Value().full_external_force[offset + 2U], -0.5);
+    EXPECT_DOUBLE_EQ(result.Value().full_external_force[offset + 3U], 0.0);
+    EXPECT_DOUBLE_EQ(result.Value().full_external_force[offset + 4U], 0.0);
+    EXPECT_DOUBLE_EQ(result.Value().full_external_force[offset + 5U], 0.0);
+  }
+  EXPECT_EQ(result.Value().full_external_tangent.Rows(), 24U);
+  EXPECT_EQ(result.Value().full_external_tangent.Columns(), 24U);
+  for (const double value : result.Value().full_external_tangent.Values()) {
+    EXPECT_EQ(value, 0.0);
+  }
+}
+
+TEST(SurfacePressureAssembly,
+     AddsRepeatedRecordsAndAppliesUnitFactorExactlyOnce) {
+  auto fixture = MakeFixture({6.0, 2.0});
+  fesa::Vector coordinates{fixture.dofs->FullDofCount()};
+  fesa::SerialParallelFor serial;
+
+  const auto result = fesa::AssembleSurfacePressure(
+      *fixture.model, *fixture.dofs, coordinates, 0.5,
+      fesa::SurfacePressureConfiguration::kReference, serial);
+
+  ASSERT_TRUE(result.HasValue());
+  for (std::size_t node = 0U; node < 4U; ++node) {
+    EXPECT_DOUBLE_EQ(result.Value().full_external_force[node * 6U + 2U],
+                     -1.0);
+  }
+}
+
+TEST(SurfacePressureAssembly,
+     IsBitwiseStableForSerialAndParallelLocalEvaluation) {
+  auto fixture = MakeFixture({1.0e16, -1.0e16, 3.5});
+  fesa::Vector coordinates{fixture.dofs->FullDofCount()};
+  fesa::SerialParallelFor serial;
+  fesa::TbbParallelFor parallel;
+
+  const auto serial_result = fesa::AssembleSurfacePressure(
+      *fixture.model, *fixture.dofs, coordinates, 1.0,
+      fesa::SurfacePressureConfiguration::kReference, serial);
+  const auto parallel_result = fesa::AssembleSurfacePressure(
+      *fixture.model, *fixture.dofs, coordinates, 1.0,
+      fesa::SurfacePressureConfiguration::kReference, parallel);
+
+  ASSERT_TRUE(serial_result.HasValue());
+  ASSERT_TRUE(parallel_result.HasValue());
+  EXPECT_EQ(Values(serial_result.Value().full_external_force),
+            Values(parallel_result.Value().full_external_force));
+  EXPECT_EQ(serial_result.Value().full_external_tangent.RowOffsets(),
+            parallel_result.Value().full_external_tangent.RowOffsets());
+  EXPECT_EQ(serial_result.Value().full_external_tangent.ColumnIndices(),
+            parallel_result.Value().full_external_tangent.ColumnIndices());
+  EXPECT_EQ(serial_result.Value().full_external_tangent.Values(),
+            parallel_result.Value().full_external_tangent.Values());
+  EXPECT_TRUE(serial_result.Value().diagnostics.empty());
+  EXPECT_TRUE(parallel_result.Value().diagnostics.empty());
+}
diff --git a/tests/unit/io/hdf5/hdf5_results_writer_test.cpp b/tests/unit/io/hdf5/hdf5_results_writer_test.cpp
index 2cef031..5250cd4 100644
--- a/tests/unit/io/hdf5/hdf5_results_writer_test.cpp
+++ b/tests/unit/io/hdf5/hdf5_results_writer_test.cpp
@@ -989,6 +989,9 @@ TEST(Hdf5ResultsWriter, WritesExactSchemaShapesAttributesAndIdentity) {
                       "/steps/Step-1/frames/0/element/transverse_shear_stress",
                       H5P_DEFAULT),
             0);
+  EXPECT_EQ(H5Lexists(file.Get(), "/steps/Step-1/frames/0/surface_pressure",
+                      H5P_DEFAULT),
+            0);
 }
 
 TEST(Hdf5ResultsWriter, WritesMandatoryOutputsDespiteOutputRequests) {