24065262673

Committed 06 Apr 2026 04:56PM UTC coverage: 75.039%. Remained the same

Build # 24065262673

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #1014 from wdhawkins/ff_fixes

Numerous bug fixes for field function and interpolation routines

Coverage Stats

48 of 71 new or added lines in 7 files covered. (67.61%)

152 existing lines in 9 files now uncovered.

20984 of 27964 relevant lines covered (75.04%)

66658029.72 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

31.18

/framework/field_functions/field_function_grid_based.cc

// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
// SPDX-License-Identifier: MIT

#include "framework/field_functions/field_function_grid_based.h"
#include "framework/math/spatial_discretization/finite_element/piecewise_linear/piecewise_linear_continuous.h"
#include "framework/math/spatial_discretization/finite_element/piecewise_linear/piecewise_linear_discontinuous.h"
#include "framework/math/spatial_discretization/spatial_discretization.h"
#include "framework/mesh/mesh.h"
#include "framework/mesh/mesh_continuum/mesh_continuum.h"
#include "framework/mesh/mesh_continuum/grid_vtk_utils.h"
#include "framework/object_factory.h"
#include "framework/logging/log.h"
#include "framework/utils/error.h"
#include "framework/runtime.h"
#include <limits>
#include <petsc.h>
#include <vtkUnstructuredGrid.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <vtkDoubleArray.h>

namespace opensn
{

InputParameters
FieldFunctionGridBased::GetInputParameters()
{
  InputParameters params = FieldFunction::GetInputParameters();

  params.AddOptionalParameter(
    "discretization", "PWLC", "The spatial discretization type to be used");
  params.AddOptionalParameter(
    "coord_sys", "cartesian", "Coordinate system to apply to element mappings");
  params.AddOptionalParameter("quadrature_order",
                              0,
                              "If supplied, will overwrite the default for the "
                              "specific discretization-coordinate system combination.");

  params.AddOptionalParameter(
    "initial_value", 0.0, "The initial value to assign to the field function");

  params.ConstrainParameterRange("discretization", AllowableRangeList::New({"PWLC", "PWLD"}));
  params.ConstrainParameterRange(
    "coord_sys", AllowableRangeList::New({"cartesian", "cylindrical", "spherical"}));

  return params;
}

FieldFunctionGridBased::FieldFunctionGridBased(const InputParameters& params)
  : FieldFunction(params),
    discretization_(MakeSpatialDiscretization(params)),
    ghosted_field_vector_(MakeFieldVector(*discretization_, GetUnknownManager())),
    local_grid_bounding_box_(params.GetSharedPtrParam<MeshContinuum>("mesh")->GetLocalBoundingBox())
{
  ghosted_field_vector_->Set(params.GetParamValue<double>("initial_value"));
}

FieldFunctionGridBased::FieldFunctionGridBased(
  std::string name, std::shared_ptr<SpatialDiscretization>& discretization_ptr, Unknown unknown)
  : FieldFunction(std::move(name), std::move(unknown)),
    discretization_(discretization_ptr),
    ghosted_field_vector_(MakeFieldVector(*discretization_, GetUnknownManager())),
    local_grid_bounding_box_(discretization_->GetGrid()->GetLocalBoundingBox())
{
}

FieldFunctionGridBased::FieldFunctionGridBased(
  std::string name,
  std::shared_ptr<SpatialDiscretization>& discretization_ptr,
  Unknown unknown,
  const std::vector<double>& field_vector)
  : FieldFunction(std::move(name), std::move(unknown)),
    discretization_(discretization_ptr),
    ghosted_field_vector_(MakeFieldVector(*discretization_, GetUnknownManager())),
    local_grid_bounding_box_(discretization_->GetGrid()->GetLocalBoundingBox())
{
  OpenSnInvalidArgumentIf(field_vector.size() != ghosted_field_vector_->GetLocalSize(),
                          "Constructor called with incompatible size field vector.");

  ghosted_field_vector_->Set(field_vector);
}

FieldFunctionGridBased::FieldFunctionGridBased(
  std::string name,
  std::shared_ptr<SpatialDiscretization>& discretization_ptr,
  Unknown unknown,
  double field_value)
  : FieldFunction(std::move(name), std::move(unknown)),
    discretization_(discretization_ptr),
    ghosted_field_vector_(MakeFieldVector(*discretization_, GetUnknownManager())),
    local_grid_bounding_box_(discretization_->GetGrid()->GetLocalBoundingBox())
{
  ghosted_field_vector_->Set(field_value);
}

const SpatialDiscretization&
FieldFunctionGridBased::GetSpatialDiscretization() const
{
  return *discretization_;
}

std::vector<double>&
FieldFunctionGridBased::GetLocalFieldVector()
{
  return ghosted_field_vector_->GetLocalSTLData();
}

const std::vector<double>&
FieldFunctionGridBased::GetLocalFieldVector() const
{
  return ghosted_field_vector_->GetLocalSTLData();
}

std::vector<double>
FieldFunctionGridBased::GetGhostedFieldVector() const
{
  return ghosted_field_vector_->GetLocalSTLData();
}

void
FieldFunctionGridBased::UpdateFieldVector(const std::vector<double>& field_vector)
{
  OpenSnInvalidArgumentIf(field_vector.size() < ghosted_field_vector_->GetLocalSize(),
                          "Attempted update with a vector of insufficient size.");

  ghosted_field_vector_->Set(field_vector);
  ghosted_field_vector_->CommunicateGhostEntries();
}

void
FieldFunctionGridBased::UpdateFieldVector(const Vec& field_vector)
{
  ghosted_field_vector_->CopyLocalValues(field_vector);

  ghosted_field_vector_->CommunicateGhostEntries();
}

std::vector<double>
FieldFunctionGridBased::GetPointValue(const Vector3& point) const
{
  const auto& uk_man = GetUnknownManager();
  const size_t num_components = uk_man.GetTotalUnknownStructureSize();

  size_t local_num_point_hits = 0;
  std::vector<double> local_point_value(num_components, 0.0);

  const auto& xyz_min = local_grid_bounding_box_.first;
  const auto& xyz_max = local_grid_bounding_box_.second;

  const auto xmin = xyz_min.x;
  const auto ymin = xyz_min.y;
  const auto zmin = xyz_min.z;

  const auto xmax = xyz_max.x;
  const auto ymax = xyz_max.y;
  const auto zmax = xyz_max.z;

  const auto& field_vector = *ghosted_field_vector_;

  if (point.x >= xmin and point.x <= xmax and point.y >= ymin and point.y <= ymax and
      point.z >= zmin and point.z <= zmax)
  {
    const auto& grid = discretization_->GetGrid();
    for (const auto& cell : grid->local_cells)
    {
      if (grid->CheckPointInsideCell(cell, point))
      {
        const auto& cell_mapping = discretization_->GetCellMapping(cell);
        Vector<double> shape_values;
        cell_mapping.ShapeValues(point, shape_values);

        local_num_point_hits += 1;

        const auto num_nodes = cell_mapping.GetNumNodes();
        for (size_t c = 0; c < num_components; ++c)
        {
          for (size_t j = 0; j < num_nodes; ++j)
          {
            const auto dof_map = discretization_->MapDOFLocal(cell, j, uk_man, 0, c);
            const double dof_value = field_vector[dof_map];

            local_point_value[c] += dof_value * shape_values(j);
          } // for node i
        } // for component c
      } // if inside cell
    } // for cell
  } // if in bounding box

  // Communicate number of point hits
  size_t global_num_point_hits = 0;
  mpi_comm.all_reduce(local_num_point_hits, global_num_point_hits, mpi::op::sum<size_t>());
  OpenSnLogicalErrorIf(global_num_point_hits == 0,
                       "FieldFunctionGridBased::GetPointValue: No cell identified containing "
                       "the specified point.");

  std::vector<double> global_point_value(num_components, 0.0);
  OpenSnLogicalErrorIf(num_components > static_cast<size_t>(std::numeric_limits<int>::max()),
                       "FieldFunctionGridBased::GetPointValue: Number of components exceeds "
                       "MPI reduction count limit.");
  const auto point_value_size = static_cast<int>(num_components);
  mpi_comm.all_reduce(
    local_point_value.data(), point_value_size, global_point_value.data(), mpi::op::sum<double>());

  Scale(global_point_value, 1.0 / static_cast<double>(global_num_point_hits));

  return global_point_value;
}

double
FieldFunctionGridBased::Evaluate(const Cell& cell, const Vector3& position, int component) const
{
  const auto& field_vector = *ghosted_field_vector_;

  const auto& cell_mapping = discretization_->GetCellMapping(cell);

  Vector<double> shape_values;
  cell_mapping.ShapeValues(position, shape_values);

  double value = 0.0;
  const size_t num_nodes = cell_mapping.GetNumNodes();
  for (size_t j = 0; j < num_nodes; ++j)
  {
    const auto dof_map = discretization_->MapDOFLocal(cell, j, GetUnknownManager(), 0, component);
    value += field_vector[dof_map] * shape_values(j);
  }

  return value;
}

void
FieldFunctionGridBased::ExportMultipleToPVTU(
  const std::string& file_base_name,
  const std::vector<std::shared_ptr<const FieldFunctionGridBased>>& ff_list)
{
  const std::string fname = "FieldFunctionGridBased::ExportMultipleToPVTU";
  log.Log() << "Exporting field functions to PVTU with file base \"" << file_base_name << "\"";

  if (ff_list.empty())
    throw std::logic_error(fname + ": Cannot be used with empty field-function"
                                   " list");

  // Setup master and check slaves
  const auto& master_ff_ptr = ff_list.front();
  const auto& master_ff = *master_ff_ptr;

  for (const auto& ff_ptr : ff_list)
    if (ff_ptr != master_ff_ptr)
      if (ff_ptr->discretization_->GetGrid() != master_ff_ptr->discretization_->GetGrid())
        throw std::logic_error(fname +
                               ": Cannot be used with field functions based on different grids.");

  // Get grid
  const auto& grid = master_ff.discretization_->GetGrid();

  auto ugrid = PrepareVtkUnstructuredGrid(grid);

  // Upload cell/point data
  auto* point_data = ugrid->GetPointData();
  for (const auto& ff_ptr : ff_list)
  {
    const auto field_vector = ff_ptr->GetGhostedFieldVector();

    const auto& uk_man = ff_ptr->GetUnknownManager();
    const auto& unknown = ff_ptr->GetUnknown();
    const auto& sdm = ff_ptr->discretization_;
    const size_t num_comps = unknown.GetNumComponents();

    for (uint c = 0; c < num_comps; ++c)
    {
      std::string component_name = ff_ptr->GetName() + unknown.name;
      if (num_comps > 1)
        component_name += unknown.component_names[c];

      vtkNew<vtkDoubleArray> point_array;

      point_array->SetName(component_name.c_str());

      // Populate the array here
      for (const auto& cell : grid->local_cells)
      {
        const size_t num_nodes = sdm->GetCellNumNodes(cell);

        if (num_nodes == cell.vertex_ids.size())
        {
          for (size_t n = 0; n < num_nodes; ++n)
          {
            const auto nmap = sdm->MapDOFLocal(cell, n, uk_man, 0, c);

            const double field_value = field_vector[nmap];

            point_array->InsertNextValue(field_value);
          } // for node
        }
        else
        {
          double node_average = 0.0;
          for (size_t n = 0; n < num_nodes; ++n)
          {
            const auto nmap = sdm->MapDOFLocal(cell, n, uk_man, 0, c);

            const double field_value = field_vector[nmap];
            node_average += field_value;
          } // for node
          node_average /= static_cast<double>(num_nodes);
          for (std::size_t n = 0; n < cell.vertex_ids.size(); ++n)
          {
            point_array->InsertNextValue(node_average);
          } // for vertex
        }

      } // for cell

      point_data->AddArray(point_array);
    } // for component
  } // for ff_ptr

  WritePVTUFiles(ugrid, file_base_name);

  log.Log() << "Done exporting field functions to PVTU.";
  opensn::mpi_comm.barrier();
}

std::shared_ptr<SpatialDiscretization>
FieldFunctionGridBased::MakeSpatialDiscretization(const InputParameters& params)
{
  const auto grid = params.GetSharedPtrParam<MeshContinuum>("mesh");
  const auto sdm_type = params.GetParamValue<std::string>("discretization");

  std::string cs = "cartesian";
  if (params.IsParameterValid("coord_sys"))
    cs = params.GetParamValue<std::string>("coord_sys");

  QuadratureOrder q_order = QuadratureOrder::SECOND;

  if (params.IsParameterValid("quadrature_order"))
  {
    const auto max_order = static_cast<uint32_t>(QuadratureOrder::FORTYTHIRD);
    const auto q_order_int = params.GetParamValue<uint32_t>("quadrature_order");
    OpenSnInvalidArgumentIf(q_order_int > max_order,
                            "Invalid quadrature point order " + std::to_string(q_order_int));
    q_order = static_cast<QuadratureOrder>(q_order_int);
  }
  else // Defaulted
  {
    if (cs == "cartesian")
      q_order = QuadratureOrder::SECOND;
    else if (cs == "cylindrical")
      q_order = QuadratureOrder::THIRD;
    else if (cs == "spherical")
      q_order = QuadratureOrder::FOURTH;
    else
      throw std::logic_error("Unrecognized coordinate system.");
  }

  if (sdm_type == "PWLC")
    return PieceWiseLinearContinuous::New(grid, q_order);
  else if (sdm_type == "PWLD")
    return PieceWiseLinearDiscontinuous::New(grid, q_order);

  // If not returned by now
  OpenSnInvalidArgument("Unsupported discretization \"" + sdm_type + "\"");
}

std::unique_ptr<GhostedParallelSTLVector>
FieldFunctionGridBased::MakeFieldVector(const SpatialDiscretization& discretization,
                                        const UnknownManager& uk_man)
{
  auto field = std::make_unique<GhostedParallelSTLVector>(discretization.GetNumLocalDOFs(uk_man),
                                                          discretization.GetNumGlobalDOFs(uk_man),
                                                          discretization.GetGhostDOFIndices(uk_man),
                                                          mpi_comm);

  return field;
}

} // namespace opensn

1	// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
2	// SPDX-License-Identifier: MIT
3
4	#include "framework/field_functions/field_function_grid_based.h"
5	#include "framework/math/spatial_discretization/finite_element/piecewise_linear/piecewise_linear_continuous.h"
6	#include "framework/math/spatial_discretization/finite_element/piecewise_linear/piecewise_linear_discontinuous.h"
7	#include "framework/math/spatial_discretization/spatial_discretization.h"
8	#include "framework/mesh/mesh.h"
9	#include "framework/mesh/mesh_continuum/mesh_continuum.h"
10	#include "framework/mesh/mesh_continuum/grid_vtk_utils.h"
11	#include "framework/object_factory.h"
12	#include "framework/logging/log.h"
13	#include "framework/utils/error.h"
14	#include "framework/runtime.h"
15	#include <limits>
16	#include <petsc.h>
17	#include <vtkUnstructuredGrid.h>
18	#include <vtkCellData.h>
19	#include <vtkPointData.h>
20	#include <vtkDoubleArray.h>
21
22	namespace opensn
23	{
24
25	InputParameters
26	FieldFunctionGridBased::GetInputParameters()	×
27	{
28	InputParameters params = FieldFunction::GetInputParameters();	×
29
30	params.AddOptionalParameter(	×
31	"discretization", "PWLC", "The spatial discretization type to be used");
32	params.AddOptionalParameter(	×
33	"coord_sys", "cartesian", "Coordinate system to apply to element mappings");
34	params.AddOptionalParameter("quadrature_order",	×
35	0,
36	"If supplied, will overwrite the default for the "
37	"specific discretization-coordinate system combination.");
38
39	params.AddOptionalParameter(	×
40	"initial_value", 0.0, "The initial value to assign to the field function");
41
42	params.ConstrainParameterRange("discretization", AllowableRangeList::New({"PWLC", "PWLD"}));	×
43	params.ConstrainParameterRange(	×
44	"coord_sys", AllowableRangeList::New({"cartesian", "cylindrical", "spherical"}));	×
45
46	return params;	×
47	}	×
48
49	FieldFunctionGridBased::FieldFunctionGridBased(const InputParameters& params)	×
50	: FieldFunction(params),
51	discretization_(MakeSpatialDiscretization(params)),	×
52	ghosted_field_vector_(MakeFieldVector(*discretization_, GetUnknownManager())),	×
53	local_grid_bounding_box_(params.GetSharedPtrParam<MeshContinuum>("mesh")->GetLocalBoundingBox())	×
54	{
55	ghosted_field_vector_->Set(params.GetParamValue<double>("initial_value"));	×
56	}	×
57
58	FieldFunctionGridBased::FieldFunctionGridBased(	42,320✔
59	std::string name, std::shared_ptr<SpatialDiscretization>& discretization_ptr, Unknown unknown)	42,320✔
60	: FieldFunction(std::move(name), std::move(unknown)),
61	discretization_(discretization_ptr),	84,640✔
62	ghosted_field_vector_(MakeFieldVector(*discretization_, GetUnknownManager())),	42,320✔
63	local_grid_bounding_box_(discretization_->GetGrid()->GetLocalBoundingBox())	126,960✔
64	{
65	}	42,320✔
66
67	FieldFunctionGridBased::FieldFunctionGridBased(	×
68	std::string name,
69	std::shared_ptr<SpatialDiscretization>& discretization_ptr,
70	Unknown unknown,
71	const std::vector<double>& field_vector)	×
72	: FieldFunction(std::move(name), std::move(unknown)),
73	discretization_(discretization_ptr),	×
74	ghosted_field_vector_(MakeFieldVector(*discretization_, GetUnknownManager())),	×
75	local_grid_bounding_box_(discretization_->GetGrid()->GetLocalBoundingBox())	×
76	{
77	OpenSnInvalidArgumentIf(field_vector.size() != ghosted_field_vector_->GetLocalSize(),	×
78	"Constructor called with incompatible size field vector.");
79
80	ghosted_field_vector_->Set(field_vector);	×
81	}	×
82
83	FieldFunctionGridBased::FieldFunctionGridBased(	×
84	std::string name,
85	std::shared_ptr<SpatialDiscretization>& discretization_ptr,
86	Unknown unknown,
87	double field_value)	×
88	: FieldFunction(std::move(name), std::move(unknown)),
89	discretization_(discretization_ptr),	×
90	ghosted_field_vector_(MakeFieldVector(*discretization_, GetUnknownManager())),	×
91	local_grid_bounding_box_(discretization_->GetGrid()->GetLocalBoundingBox())	×
92	{
93	ghosted_field_vector_->Set(field_value);	×
94	}	×
95
96	const SpatialDiscretization&
97	FieldFunctionGridBased::GetSpatialDiscretization() const	14,680✔
98	{
99	return *discretization_;	14,680✔
100	}
101
102	std::vector<double>&
103	FieldFunctionGridBased::GetLocalFieldVector()	×
104	{
105	return ghosted_field_vector_->GetLocalSTLData();	×
106	}
107
108	const std::vector<double>&
109	FieldFunctionGridBased::GetLocalFieldVector() const	×
110	{
111	return ghosted_field_vector_->GetLocalSTLData();	×
112	}
113
114	std::vector<double>
115	FieldFunctionGridBased::GetGhostedFieldVector() const	13,281✔
116	{
117	return ghosted_field_vector_->GetLocalSTLData();	13,281✔
118	}
119
120	void
121	FieldFunctionGridBased::UpdateFieldVector(const std::vector<double>& field_vector)	42,320✔
122	{
123	OpenSnInvalidArgumentIf(field_vector.size() < ghosted_field_vector_->GetLocalSize(),	42,320✔
124	"Attempted update with a vector of insufficient size.");
125
126	ghosted_field_vector_->Set(field_vector);	42,320✔
127	ghosted_field_vector_->CommunicateGhostEntries();	42,320✔
128	}	42,320✔
129
130	void
131	FieldFunctionGridBased::UpdateFieldVector(const Vec& field_vector)	×
132	{
133	ghosted_field_vector_->CopyLocalValues(field_vector);	×
134
135	ghosted_field_vector_->CommunicateGhostEntries();	×
136	}	×
137
138	std::vector<double>
139	FieldFunctionGridBased::GetPointValue(const Vector3& point) const	×
140	{
141	const auto& uk_man = GetUnknownManager();	×
142	const size_t num_components = uk_man.GetTotalUnknownStructureSize();	×
143
144	size_t local_num_point_hits = 0;	×
145	std::vector<double> local_point_value(num_components, 0.0);	×
146
147	const auto& xyz_min = local_grid_bounding_box_.first;	×
148	const auto& xyz_max = local_grid_bounding_box_.second;	×
149
150	const auto xmin = xyz_min.x;	×
151	const auto ymin = xyz_min.y;	×
152	const auto zmin = xyz_min.z;	×
153
154	const auto xmax = xyz_max.x;	×
155	const auto ymax = xyz_max.y;	×
156	const auto zmax = xyz_max.z;	×
157
158	const auto& field_vector = *ghosted_field_vector_;	×
159
160	if (point.x >= xmin and point.x <= xmax and point.y >= ymin and point.y <= ymax and	×
161	point.z >= zmin and point.z <= zmax)	×
162	{
163	const auto& grid = discretization_->GetGrid();	×
164	for (const auto& cell : grid->local_cells)	×
165	{
166	if (grid->CheckPointInsideCell(cell, point))	×
167	{
168	const auto& cell_mapping = discretization_->GetCellMapping(cell);	×
169	Vector<double> shape_values;	×
170	cell_mapping.ShapeValues(point, shape_values);	×
171
172	local_num_point_hits += 1;	×
173
174	const auto num_nodes = cell_mapping.GetNumNodes();	×
175	for (size_t c = 0; c < num_components; ++c)	×
176	{
177	for (size_t j = 0; j < num_nodes; ++j)	×
178	{
179	const auto dof_map = discretization_->MapDOFLocal(cell, j, uk_man, 0, c);	×
180	const double dof_value = field_vector[dof_map];	×
181
182	local_point_value[c] += dof_value * shape_values(j);	×
183	} // for node i
184	} // for component c
185	} // if inside cell	×
186	} // for cell
187	} // if in bounding box	×
188
189	// Communicate number of point hits
190	size_t global_num_point_hits = 0;	×
191	mpi_comm.all_reduce(local_num_point_hits, global_num_point_hits, mpi::op::sum<size_t>());	×
NEW 192	OpenSnLogicalErrorIf(global_num_point_hits == 0,	×
193	"FieldFunctionGridBased::GetPointValue: No cell identified containing "
194	"the specified point.");
195
196	std::vector<double> global_point_value(num_components, 0.0);	×
NEW 197	OpenSnLogicalErrorIf(num_components > static_cast<size_t>(std::numeric_limits<int>::max()),	×
198	"FieldFunctionGridBased::GetPointValue: Number of components exceeds "
199	"MPI reduction count limit.");
NEW 200	const auto point_value_size = static_cast<int>(num_components);	×
201	mpi_comm.all_reduce(	×
NEW 202	local_point_value.data(), point_value_size, global_point_value.data(), mpi::op::sum<double>());	×
203
204	Scale(global_point_value, 1.0 / static_cast<double>(global_num_point_hits));	×
205
206	return global_point_value;	×
207	}	×
208
209	double
210	FieldFunctionGridBased::Evaluate(const Cell& cell, const Vector3& position, int component) const	×
211	{
212	const auto& field_vector = *ghosted_field_vector_;	×
213
214	const auto& cell_mapping = discretization_->GetCellMapping(cell);	×
215
216	Vector<double> shape_values;	×
217	cell_mapping.ShapeValues(position, shape_values);	×
218
219	double value = 0.0;	×
220	const size_t num_nodes = cell_mapping.GetNumNodes();	×
221	for (size_t j = 0; j < num_nodes; ++j)	×
222	{
223	const auto dof_map = discretization_->MapDOFLocal(cell, j, GetUnknownManager(), 0, component);	×
224	value += field_vector[dof_map] * shape_values(j);	×
225	}
226
227	return value;	×
228	}	×
229
230	void
231	FieldFunctionGridBased::ExportMultipleToPVTU(	58✔
232	const std::string& file_base_name,
233	const std::vector<std::shared_ptr<const FieldFunctionGridBased>>& ff_list)
234	{
235	const std::string fname = "FieldFunctionGridBased::ExportMultipleToPVTU";	58✔
236	log.Log() << "Exporting field functions to PVTU with file base \"" << file_base_name << "\"";	174✔
237
238	if (ff_list.empty())	58✔
239	throw std::logic_error(fname + ": Cannot be used with empty field-function"	×
240	" list");	×
241
242	// Setup master and check slaves
243	const auto& master_ff_ptr = ff_list.front();	58✔
244	const auto& master_ff = *master_ff_ptr;	58✔
245
246	for (const auto& ff_ptr : ff_list)	5,998✔
247	if (ff_ptr != master_ff_ptr)	5,940✔
248	if (ff_ptr->discretization_->GetGrid() != master_ff_ptr->discretization_->GetGrid())	17,646✔
249	throw std::logic_error(fname +	×
250	": Cannot be used with field functions based on different grids.");	×
251
252	// Get grid
253	const auto& grid = master_ff.discretization_->GetGrid();	58✔
254
255	auto ugrid = PrepareVtkUnstructuredGrid(grid);	116✔
256
257	// Upload cell/point data
258	auto* point_data = ugrid->GetPointData();	58✔
259	for (const auto& ff_ptr : ff_list)	5,998✔
260	{
261	const auto field_vector = ff_ptr->GetGhostedFieldVector();	5,940✔
262
263	const auto& uk_man = ff_ptr->GetUnknownManager();	5,940✔
264	const auto& unknown = ff_ptr->GetUnknown();	5,940✔
265	const auto& sdm = ff_ptr->discretization_;	5,940✔
266	const size_t num_comps = unknown.GetNumComponents();	5,940✔
267
268	for (uint c = 0; c < num_comps; ++c)	11,880✔
269	{
270	std::string component_name = ff_ptr->GetName() + unknown.name;	5,940✔
271	if (num_comps > 1)	5,940✔
272	component_name += unknown.component_names[c];	×
273
274	vtkNew<vtkDoubleArray> point_array;	5,940✔
275
276	point_array->SetName(component_name.c_str());	5,940✔
277
278	// Populate the array here
279	for (const auto& cell : grid->local_cells)	1,276,198✔
280	{
281	const size_t num_nodes = sdm->GetCellNumNodes(cell);	1,270,258✔
282
283	if (num_nodes == cell.vertex_ids.size())	1,270,258✔
284	{
285	for (size_t n = 0; n < num_nodes; ++n)	6,238,430✔
286	{
287	const auto nmap = sdm->MapDOFLocal(cell, n, uk_man, 0, c);	4,968,172✔
288
289	const double field_value = field_vector[nmap];	4,968,172✔
290
291	point_array->InsertNextValue(field_value);	4,968,172✔
292	} // for node
293	}
294	else
295	{
296	double node_average = 0.0;
297	for (size_t n = 0; n < num_nodes; ++n)	×
298	{
299	const auto nmap = sdm->MapDOFLocal(cell, n, uk_man, 0, c);	×
300
301	const double field_value = field_vector[nmap];	×
302	node_average += field_value;	×
303	} // for node
304	node_average /= static_cast<double>(num_nodes);	×
305	for (std::size_t n = 0; n < cell.vertex_ids.size(); ++n)	×
306	{
307	point_array->InsertNextValue(node_average);	×
308	} // for vertex
309	}
310
311	} // for cell
312
313	point_data->AddArray(point_array);	5,940✔
314	} // for component	5,940✔
315	} // for ff_ptr	5,940✔
316
317	WritePVTUFiles(ugrid, file_base_name);	58✔
318
319	log.Log() << "Done exporting field functions to PVTU.";	116✔
320	opensn::mpi_comm.barrier();	58✔
321	}	116✔
322
323	std::shared_ptr<SpatialDiscretization>
324	FieldFunctionGridBased::MakeSpatialDiscretization(const InputParameters& params)	×
325	{
326	const auto grid = params.GetSharedPtrParam<MeshContinuum>("mesh");	×
327	const auto sdm_type = params.GetParamValue<std::string>("discretization");	×
328
329	std::string cs = "cartesian";	×
330	if (params.IsParameterValid("coord_sys"))	×
331	cs = params.GetParamValue<std::string>("coord_sys");	×
332
333	QuadratureOrder q_order = QuadratureOrder::SECOND;	×
334
335	if (params.IsParameterValid("quadrature_order"))	×
336	{
337	const auto max_order = static_cast<uint32_t>(QuadratureOrder::FORTYTHIRD);	×
338	const auto q_order_int = params.GetParamValue<uint32_t>("quadrature_order");	×
339	OpenSnInvalidArgumentIf(q_order_int > max_order,	×
340	"Invalid quadrature point order " + std::to_string(q_order_int));
341	q_order = static_cast<QuadratureOrder>(q_order_int);	×
342	}
343	else // Defaulted
344	{
345	if (cs == "cartesian")	×
346	q_order = QuadratureOrder::SECOND;
347	else if (cs == "cylindrical")	×
348	q_order = QuadratureOrder::THIRD;
349	else if (cs == "spherical")	×
350	q_order = QuadratureOrder::FOURTH;
351	else
352	throw std::logic_error("Unrecognized coordinate system.");	×
353	}
354
355	if (sdm_type == "PWLC")	×
356	return PieceWiseLinearContinuous::New(grid, q_order);	×
357	else if (sdm_type == "PWLD")	×
358	return PieceWiseLinearDiscontinuous::New(grid, q_order);	×
359
360	// If not returned by now
361	OpenSnInvalidArgument("Unsupported discretization \"" + sdm_type + "\"");	×
362	}	×
363
364	std::unique_ptr<GhostedParallelSTLVector>
365	FieldFunctionGridBased::MakeFieldVector(const SpatialDiscretization& discretization,	42,320✔
366	const UnknownManager& uk_man)
367	{
368	auto field = std::make_unique<GhostedParallelSTLVector>(discretization.GetNumLocalDOFs(uk_man),	169,280✔
369	discretization.GetNumGlobalDOFs(uk_man),	84,640✔
370	discretization.GetGhostDOFIndices(uk_man),	42,320✔
371	mpi_comm);	42,320✔
372
373	return field;	42,320✔
374	}
375
376	} // namespace opensn

Open-Sn / opensn / 24065262673

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous