OGS
ComputeNaturalCoordsAlgorithm.h
Go to the documentation of this file.
1
10#pragma once
11
12#include <vtkCell.h>
13#include <vtkLine.h>
14#include <vtkPointData.h>
15
16#include <algorithm>
17
19#include "FindCellsForPoints.h"
23
25{
27{
28 Eigen::MatrixXd natural_coords;
29 Eigen::MatrixXd real_coords;
30 Eigen::VectorXd initial_anchor_stress;
31 Eigen::VectorXd maximum_anchor_stress;
32 Eigen::VectorXd residual_anchor_stress;
34 Eigen::VectorXd anchor_stiffness;
35 Eigen::VectorX<vtkIdType> bulk_element_ids;
36 Eigen::VectorX<vtkIdType> point_cloud_node_ids;
37 bool success;
38};
39
41{
42public:
43 void append(Eigen::Vector3d const& natural_coords,
44 Eigen::Vector3d const& real_coords,
45 vtkIdType bulk_element_id,
46 vtkIdType point_cloud_node_id)
47 {
48 natural_coordss.push_back(natural_coords);
49 real_coordss.push_back(real_coords);
50 bulk_element_ids.push_back(bulk_element_id);
51 point_cloud_node_ids.push_back(point_cloud_node_id);
52 }
53
54 void fail() { success = false; }
55
57 Eigen::VectorXd const& initial_anchor_stress,
58 Eigen::VectorXd const& maximum_anchor_stress,
59 Eigen::VectorXd const& residual_anchor_stress,
60 Eigen::VectorXd const& anchor_cross_sectional_area,
61 Eigen::VectorXd const& anchor_stiffness)
62 {
63 auto bei = Eigen::Map<Eigen::VectorX<vtkIdType>>(
64 bulk_element_ids.data(), bulk_element_ids.size());
65 auto pcni = Eigen::Map<Eigen::VectorX<vtkIdType>>(
67 return {conv(natural_coordss),
69 initial_anchor_stress,
70 maximum_anchor_stress,
71 residual_anchor_stress,
72 anchor_cross_sectional_area,
73 anchor_stiffness,
74 bei,
75 pcni,
76 success};
77 }
78
79private:
80 static Eigen::MatrixXd conv(std::vector<Eigen::Vector3d> const& in)
81 {
82 Eigen::MatrixXd out(in.size(), 3);
83 for (std::size_t i = 0; i < in.size(); ++i)
84 {
85 out.row(i).noalias() = in[i].transpose();
86 }
87 return out;
88 }
89
90 std::vector<Eigen::Vector3d> natural_coordss;
91 std::vector<Eigen::Vector3d> real_coordss;
92 std::vector<vtkIdType> bulk_element_ids;
93 std::vector<vtkIdType> point_cloud_node_ids;
94 bool success = true;
95};
96
98 vtkUnstructuredGrid* const bulk_mesh,
99 ComputeNaturalCoordsResult const& json_data, double const tolerance,
100 int const max_iter)
101{
102 // Check inputs ------------------------------------------------------------
103
104 if (json_data.real_coords.cols() != 3)
105 {
106 OGS_FATAL("Wrong number of input coordinates");
107 }
108
109 // General definitions -----------------------------------------------------
110 using SolverRegistry =
113
114 double const real_coords_tolerance = tolerance;
115
116 // Initialization ----------------------------------------------------------
117 FindCellsForPoints findCellsForPoints;
118 findCellsForPoints.initialize(bulk_mesh);
119
121
122 // Do computation for each point ------------------------------------------
123 for (Eigen::Index rc_idx = 0; rc_idx < json_data.real_coords.rows();
124 ++rc_idx)
125 {
126 Eigen::RowVector3d const real_coords_single_point =
127 json_data.real_coords.row(rc_idx);
128
129 // Find cells for point
130 auto const filtered_bulk_mesh =
131 findCellsForPoints.find(real_coords_single_point, tolerance);
132
133 // Check multiplicity
134 int const actual_multiplicity = filtered_bulk_mesh->GetNumberOfCells();
135 if (actual_multiplicity == 0)
136 {
137 OGS_FATAL(
138 "Did not find any cell for the point {} (coordinates: {})",
139 rc_idx, real_coords_single_point);
140 }
141 assert(actual_multiplicity > 0);
142 if (actual_multiplicity != 1)
143 {
144 INFO(
145 "Found more then one cell (namely {}) for point #{} "
146 "(coordinates: {})",
147 actual_multiplicity, rc_idx, real_coords_single_point);
148 }
149
150 // Convert to OGS
151 std::unique_ptr<MeshLib::Mesh> ogs_mesh(
153 filtered_bulk_mesh, "filtered_bulk_mesh"));
154
155 auto const* bulk_element_ids = MeshLib::bulkElementIDs(*ogs_mesh);
156
157 auto const& elements = ogs_mesh->getElements();
158
159 // Compute natural coordinates
160 // Found a solution, no need to test other elements ignoring
161 // multiplicity.
162 constexpr std::size_t elt_idx = 0;
163 auto const& element = *elements[elt_idx];
164 auto const bulk_element_id = bulk_element_ids->getComponent(elt_idx, 0);
165
166 auto const& solver = SolverRegistry::getFor(element);
167 auto const opt_natural_coords = solver.solve(
168 element, real_coords_single_point, max_iter, real_coords_tolerance);
169
170 result.append(opt_natural_coords, real_coords_single_point,
171 bulk_element_id, rc_idx);
172 }
173
174 return result.finished(
175 json_data.initial_anchor_stress, json_data.maximum_anchor_stress,
177 json_data.anchor_stiffness);
178}
179
180template <typename T>
181void addPointData(vtkUnstructuredGrid* grid, std::string const& name,
182 Eigen::MatrixX<T> const& data)
183{
184 INFO("converting {}", name);
185 vtkNew<vtkAOSDataArrayTemplate<T>> array;
186 array->SetName(name.c_str());
187 array->SetNumberOfComponents(data.cols());
188 array->SetNumberOfTuples(grid->GetNumberOfPoints());
189
190 if (grid->GetNumberOfPoints() != data.rows())
191 {
192 OGS_FATAL(
193 "Got {} rows in the table but expected {} rows, same as number of "
194 "points in the grid.",
195 data.rows(), grid->GetNumberOfPoints());
196 }
197 for (Eigen::Index i = 0; i < data.rows(); ++i)
198 {
199 // copy to contiguous storage
200 Eigen::RowVectorX<T> const row = data.row(i);
201
202 array->SetTypedTuple(i, row.data());
203 DBUG("> {}", row);
204 }
205
206 grid->GetPointData()->AddArray(array);
207}
208
209template <typename T>
210void addCellData(vtkUnstructuredGrid* grid, std::string const& name,
211 Eigen::MatrixX<T> const& data)
212{
213 INFO("converting {}", name);
214 vtkNew<vtkAOSDataArrayTemplate<T>> array;
215 array->SetName(name.c_str());
216 array->SetNumberOfComponents(data.cols());
217 array->SetNumberOfTuples(grid->GetNumberOfCells());
218
219 if (grid->GetNumberOfCells() != data.rows())
220 {
221 OGS_FATAL(
222 "Got {} rows in the table but expected {} rows, same as number of "
223 "cells in the grid.",
224 data.rows(), grid->GetNumberOfCells());
225 }
226 for (Eigen::Index i = 0; i < data.rows(); ++i)
227 {
228 // copy to contiguous storage
229 Eigen::RowVectorX<T> const row = data.row(i);
230
231 array->SetTypedTuple(i, row.data());
232 DBUG("> {}", row);
233 }
234
235 grid->GetCellData()->AddArray(array);
236}
237
244vtkSmartPointer<vtkUnstructuredGrid> toVTKGrid(
245 ComputeNaturalCoordsResult const& result)
246{
247 INFO("converting points");
248 auto const n_pts = result.natural_coords.rows();
249 if (n_pts == 0)
250 {
251 OGS_FATAL(
252 "No point was found in the mesh. Please check whether all anchors "
253 "are within the model region.");
254 }
255 else if (n_pts == 1)
256 {
257 OGS_FATAL(
258 "Only one point was found in the mesh. Please check whether all "
259 "anchors are within the model region.");
260 }
261 else if (n_pts % 2 != 0)
262 {
263 // number should be even if multiplicity is ignored
264 OGS_FATAL(
265 "Number of points is not even. Anchors are believed to consist of "
266 "a start and an end point.");
267 }
268
269 vtkNew<vtkUnstructuredGrid> grid;
270
271 // Points ------------------------------------------------------------------
272 vtkNew<vtkPoints> points;
273 points->SetNumberOfPoints(n_pts);
274
275 auto const& css = result.real_coords;
276 for (Eigen::Index i = 0; i < css.rows(); ++i)
277 {
278 points->SetPoint(i, css(i, 0), css(i, 1), css(i, 2));
279
280 if ((i % 2) == 1)
281 {
282 vtkNew<vtkLine> line;
283 line->GetPointIds()->SetId(0, i - 1);
284 line->GetPointIds()->SetId(1, i);
285 grid->InsertNextCell(line->GetCellType(), line->GetPointIds());
286 }
287 }
288
289 grid->SetPoints(points);
290
291 // Point data --------------------------------------------------------------
292 addPointData<double>(grid, "natural_coordinates", result.natural_coords);
293 addPointData<std::size_t>(grid, "bulk_element_ids",
294 result.bulk_element_ids.cast<std::size_t>());
295 addPointData<std::size_t>(grid, "point_cloud_node_ids",
296 result.point_cloud_node_ids.cast<std::size_t>());
297 addCellData<double>(grid, "initial_anchor_stress",
298 result.initial_anchor_stress);
299 addCellData<double>(grid, "maximum_anchor_stress",
300 result.maximum_anchor_stress);
301 addCellData<double>(grid, "residual_anchor_stress",
303 addCellData<double>(grid, "anchor_cross_sectional_area",
305 addCellData<double>(grid, "anchor_stiffness", result.anchor_stiffness);
306 return grid;
307}
308} // namespace ApplicationUtils
#define OGS_FATAL(...)
Definition Error.h:26
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Definition of the VtkMeshConverter class.
ComputeNaturalCoordsResult finished(Eigen::VectorXd const &initial_anchor_stress, Eigen::VectorXd const &maximum_anchor_stress, Eigen::VectorXd const &residual_anchor_stress, Eigen::VectorXd const &anchor_cross_sectional_area, Eigen::VectorXd const &anchor_stiffness)
static Eigen::MatrixXd conv(std::vector< Eigen::Vector3d > const &in)
void append(Eigen::Vector3d const &natural_coords, Eigen::Vector3d const &real_coords, vtkIdType bulk_element_id, vtkIdType point_cloud_node_id)
vtkSmartPointer< vtkUnstructuredGrid > find(Eigen::Vector3d const &coords, double const tolerance)
void initialize(vtkUnstructuredGrid *const bulk_mesh)
static MeshLib::Mesh * convertUnstructuredGrid(vtkUnstructuredGrid *grid, bool const compute_element_neighbors=false, std::string const &mesh_name="vtkUnstructuredGrid")
Converts a vtkUnstructuredGrid object to a Mesh.
void addPointData(vtkUnstructuredGrid *grid, std::string const &name, Eigen::MatrixX< T > const &data)
void addCellData(vtkUnstructuredGrid *grid, std::string const &name, Eigen::MatrixX< T > const &data)
vtkSmartPointer< vtkUnstructuredGrid > toVTKGrid(ComputeNaturalCoordsResult const &result)
ComputeNaturalCoordsResult computeNaturalCoords(vtkUnstructuredGrid *const bulk_mesh, ComputeNaturalCoordsResult const &json_data, double const tolerance, int const max_iter)
PropertyVector< std::size_t > const * bulkElementIDs(Mesh const &mesh)
Definition Mesh.cpp:301