diff --git a/ProcessLib/ProcessOutput.cpp b/ProcessLib/ProcessOutput.cpp index 65b483fd4e1871af80b571a314ddc52e3199f562..18de2526442ff8b2d16f11bd6c743934a9f295f9 100644 --- a/ProcessLib/ProcessOutput.cpp +++ b/ProcessLib/ProcessOutput.cpp @@ -12,6 +12,98 @@ #include "MeshLib/IO/VtkIO/VtuInterface.h" #include "NumLib/DOF/LocalToGlobalIndexMap.h" +static void addSecondaryVariableNodes( + double const t, + GlobalVector const& x, + NumLib::LocalToGlobalIndexMap const& dof_table, + ProcessLib::SecondaryVariable const& var, + std::string const& output_name, + MeshLib::Mesh& mesh) +{ + DBUG(" secondary variable %s", output_name.c_str()); + + auto& nodal_values_mesh = *MeshLib::getOrCreateMeshProperty<double>( + mesh, output_name, MeshLib::MeshItemType::Node, + var.fcts.num_components); + if (nodal_values_mesh.size() != + mesh.getNumberOfNodes() * var.fcts.num_components) + { + OGS_FATAL( + "Nodal property `%s' does not have the right number of " + "components. Expected: %d, actual: %d", + output_name.c_str(), + mesh.getNumberOfNodes() * var.fcts.num_components, + nodal_values_mesh.size()); + } + + std::unique_ptr<GlobalVector> result_cache; + auto const& nodal_values = + var.fcts.eval_field(t, x, dof_table, result_cache); + if (nodal_values_mesh.size() != + static_cast<std::size_t>(nodal_values.size())) + { + OGS_FATAL( + "Secondary variable `%s' did not evaluate to the right " + "number of components. Expected: %d, actual: %d.", + var.name.c_str(), nodal_values_mesh.size(), nodal_values.size()); + } + + // Copy result + for (GlobalIndexType i = 0; i < nodal_values.size(); ++i) + { + assert(!std::isnan(nodal_values[i])); + nodal_values_mesh[i] = nodal_values[i]; + } +} + +static void addSecondaryVariableResiduals( + double const t, + GlobalVector const& x, + NumLib::LocalToGlobalIndexMap const& dof_table, + ProcessLib::SecondaryVariable const& var, + std::string const& output_name, + MeshLib::Mesh& mesh) +{ + if (!var.fcts.eval_residuals) + return; + + DBUG(" secondary variable %s residual", output_name.c_str()); + auto const& property_name_res = output_name + "_residual"; + + auto& residuals_mesh = *MeshLib::getOrCreateMeshProperty<double>( + mesh, property_name_res, MeshLib::MeshItemType::Cell, + var.fcts.num_components); + if (residuals_mesh.size() != + mesh.getNumberOfElements() * var.fcts.num_components) + { + OGS_FATAL( + "Cell property `%s' does not have the right number of " + "components. Expected: %d, actual: %d", + property_name_res.c_str(), + mesh.getNumberOfElements() * var.fcts.num_components, + residuals_mesh.size()); + } + + std::unique_ptr<GlobalVector> result_cache; + auto const& residuals = + var.fcts.eval_residuals(t, x, dof_table, result_cache); + if (residuals_mesh.size() != static_cast<std::size_t>(residuals.size())) + { + OGS_FATAL( + "Thee residual of secondary variable `%s' did not evaluate " + "to the right number of components. Expected: %d, actual: " + "%d.", + var.name.c_str(), residuals_mesh.size(), residuals.size()); + } + + // Copy result + for (GlobalIndexType i = 0; i < residuals.size(); ++i) + { + assert(!std::isnan(residuals[i])); + residuals_mesh[i] = residuals[i]; + } +} + namespace ProcessLib { void processOutputData( @@ -101,87 +193,8 @@ void processOutputData( } #ifndef USE_PETSC - auto add_secondary_var = [&](SecondaryVariable const& var, - std::string const& output_name) { - { - DBUG(" secondary variable %s", output_name.c_str()); - - auto& nodal_values_mesh = *MeshLib::getOrCreateMeshProperty<double>( - mesh, output_name, MeshLib::MeshItemType::Node, - var.fcts.num_components); - if (nodal_values_mesh.size() != - mesh.getNumberOfNodes() * var.fcts.num_components) - { - OGS_FATAL( - "Nodal property `%s' does not have the right number of " - "components. Expected: %d, actual: %d", - output_name.c_str(), - mesh.getNumberOfNodes() * var.fcts.num_components, - nodal_values_mesh.size()); - } - - std::unique_ptr<GlobalVector> result_cache; - auto const& nodal_values = - var.fcts.eval_field(t, x, dof_table, result_cache); - if (nodal_values_mesh.size() != - static_cast<std::size_t>(nodal_values.size())) - { - OGS_FATAL( - "Secondary variable `%s' did not evaluate to the right " - "number of components. Expected: %d, actual: %d.", - var.name.c_str(), nodal_values_mesh.size(), - nodal_values.size()); - } - - // Copy result - for (GlobalIndexType i = 0; i < nodal_values.size(); ++i) - { - assert(!std::isnan(nodal_values[i])); - nodal_values_mesh[i] = nodal_values[i]; - } - } - - if (process_output.output_residuals && var.fcts.eval_residuals) - { - DBUG(" secondary variable %s residual", output_name.c_str()); - auto const& property_name_res = output_name + "_residual"; - - auto& residuals_mesh = *MeshLib::getOrCreateMeshProperty<double>( - mesh, property_name_res, MeshLib::MeshItemType::Cell, - var.fcts.num_components); - if (residuals_mesh.size() != - mesh.getNumberOfElements() * var.fcts.num_components) - { - OGS_FATAL( - "Cell property `%s' does not have the right number of " - "components. Expected: %d, actual: %d", - property_name_res.c_str(), - mesh.getNumberOfElements() * var.fcts.num_components, - residuals_mesh.size()); - } - - std::unique_ptr<GlobalVector> result_cache; - auto const& residuals = - var.fcts.eval_residuals(t, x, dof_table, result_cache); - if (residuals_mesh.size() != - static_cast<std::size_t>(residuals.size())) - { - OGS_FATAL( - "Thee residual of secondary variable `%s' did not evaluate " - "to the right number of components. Expected: %d, actual: " - "%d.", - var.name.c_str(), residuals_mesh.size(), residuals.size()); - } - - // Copy result - for (GlobalIndexType i = 0; i < residuals.size(); ++i) - { - assert(!std::isnan(residuals[i])); - residuals_mesh[i] = residuals[i]; - } - } - }; + // Secondary variables output for (auto const& external_variable_name : output_variables) { if (!already_output.insert(external_variable_name).second) { @@ -189,8 +202,16 @@ void processOutputData( continue; } - add_secondary_var(secondary_variables.get(external_variable_name), - external_variable_name); + addSecondaryVariableNodes( + t, x, dof_table, secondary_variables.get(external_variable_name), + external_variable_name, mesh); + if (process_output.output_residuals) + { + addSecondaryVariableResiduals( + t, x, dof_table, + secondary_variables.get(external_variable_name), + external_variable_name, mesh); + } } #else (void)secondary_variables;