tboudreaux/MeanField
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Export enthalpy nodal values on limiting diagnostic elements

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-17 08:40:24 -04:00
parent 9c975422c1
commit 4280c532e3
1 changed files with 27 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

27
experiments/plateau_pressure_obstruction_diagnostics.hpp
View File

@@ -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<int>(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<int>(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<double>::quiet_NaN()) << ','
<< h + dh << '\n';
}
const int extraOrder = static_cast<int>(measures.at(prefix + "eos_extra_order"));
const auto query = pressure ? EnthalpyField::template make_query<field::Enthalpy::Form::PressureForce>(
quadrature::QuadratureRole::discretization, transformation->OrderW(), std::array<int, 1>{extraOrder},