From 4280c532e3620f587a346269fc549384e04e2ba0 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 08:40:24 -0400 Subject: [PATCH] Export enthalpy nodal values on limiting diagnostic elements --- ...ateau_pressure_obstruction_diagnostics.hpp | 27 +++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/experiments/plateau_pressure_obstruction_diagnostics.hpp b/experiments/plateau_pressure_obstruction_diagnostics.hpp index 7e19854..153dada 100644 --- a/experiments/plateau_pressure_obstruction_diagnostics.hpp +++ b/experiments/plateau_pressure_obstruction_diagnostics.hpp @@ -365,6 +365,17 @@ void ProbeSurfaceTrace(State &s, const mean_field::fem::FEM &f, double n, mfem::ParGridFunction hBoundary(f.enthalpyFes.get()), dhBoundary(f.enthalpyFes.get()); adapter.scatter(boundaryH, hBoundary); adapter.scatter(boundaryDh, dhBoundary); + mfem::Vector fullH(size), fullDh(size); + for (int i = 0; i < size; ++i) { + fullH(i) = physicalBase(offset + i); + fullDh(i) = physicalDirection(offset + i); + } + mfem::ParGridFunction hFull(f.enthalpyFes.get()), dhFull(f.enthalpyFes.get()); + adapter.scatter(fullH, hFull); + adapter.scatter(fullDh, dhFull); + auto nodes = detail::Open(output / "surface_trace_nodes.csv"); + nodes << "family,element,attribute,geometry_type,fe_order,local_node,node_count,ip_x,ip_y,ip_z," + "h,dh,alpha_union,h_plus_alpha_union_dh,h_plus_dh\n"; auto limits = detail::Open(output / "surface_trace_limits.csv"); limits << "family,element,point,attribute,geometry_type,transformation_order_w,eos_extra_order,rule_points," "ip_x,ip_y,ip_z,alpha_h,full_h,full_dh,boundary_h,boundary_dh,interior_h,interior_dh," @@ -375,6 +386,22 @@ void ProbeSurfaceTrace(State &s, const mean_field::fem::FEM &f, double n, const int point = static_cast(measures.at(prefix + "limiting_point")); if (element < 0) continue; // No decreasing sample: the summary retains the infinite bound. auto *transformation = f.mesh->GetElementTransformation(element); + const auto *finiteElement = f.enthalpyFes->GetFE(element); + const auto &elementNodes = finiteElement->GetNodes(); + const double alphaUnion = measures.at("alpha_thermodynamic"); + // Evaluate the scattered physical fields at FE nodes; do not infer + // nodal values by indexing potentially signed/oriented local DOFs. + // A shared limiting element is exported once per family for identity. + for (int node = 0; node < elementNodes.GetNPoints(); ++node) { + const auto &ip = elementNodes.IntPoint(node); + const double h = hFull.GetValue(element, ip), dh = dhFull.GetValue(element, ip); + nodes << (pressure ? "pressure" : "closure") << ',' << element << ',' + << f.mesh->GetAttribute(element) << ',' << static_cast(transformation->GetGeometryType()) << ',' + << finiteElement->GetOrder() << ',' << node << ',' << elementNodes.GetNPoints() << ',' + << ip.x << ',' << ip.y << ',' << ip.z << ',' << h << ',' << dh << ',' << alphaUnion << ',' + << (std::isfinite(alphaUnion) ? h + alphaUnion * dh : std::numeric_limits::quiet_NaN()) << ',' + << h + dh << '\n'; + } const int extraOrder = static_cast(measures.at(prefix + "eos_extra_order")); const auto query = pressure ? EnthalpyField::template make_query( quadrature::QuadratureRole::discretization, transformation->OrderW(), std::array{extraOrder},