14#include <vtkPointData.h>
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)
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)
63 auto bei = Eigen::Map<Eigen::VectorX<vtkIdType>>(
65 auto pcni = Eigen::Map<Eigen::VectorX<vtkIdType>>(
69 initial_anchor_stress,
70 maximum_anchor_stress,
71 residual_anchor_stress,
72 anchor_cross_sectional_area,
80 static Eigen::MatrixXd
conv(std::vector<Eigen::Vector3d>
const& in)
82 Eigen::MatrixXd out(in.size(), 3);
83 for (std::size_t i = 0; i < in.size(); ++i)
85 out.row(i).noalias() = in[i].transpose();
98 vtkUnstructuredGrid*
const bulk_mesh,
106 OGS_FATAL(
"Wrong number of input coordinates");
110 using SolverRegistry =
114 double const real_coords_tolerance = tolerance;
123 for (Eigen::Index rc_idx = 0; rc_idx < json_data.
real_coords.rows();
126 Eigen::RowVector3d
const real_coords_single_point =
130 auto const filtered_bulk_mesh =
131 findCellsForPoints.
find(real_coords_single_point, tolerance);
134 int const actual_multiplicity = filtered_bulk_mesh->GetNumberOfCells();
135 if (actual_multiplicity == 0)
138 "Did not find any cell for the point {} (coordinates: {})",
139 rc_idx, real_coords_single_point);
141 assert(actual_multiplicity > 0);
142 if (actual_multiplicity != 1)
145 "Found more then one cell (namely {}) for point #{} "
147 actual_multiplicity, rc_idx, real_coords_single_point);
151 std::unique_ptr<MeshLib::Mesh> ogs_mesh(
153 filtered_bulk_mesh,
"filtered_bulk_mesh"));
157 auto const& elements = ogs_mesh->getElements();
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);
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);
170 result.
append(opt_natural_coords, real_coords_single_point,
171 bulk_element_id, rc_idx);
182 Eigen::MatrixX<T>
const& data)
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());
190 if (grid->GetNumberOfPoints() != data.rows())
193 "Got {} rows in the table but expected {} rows, same as number of "
194 "points in the grid.",
195 data.rows(), grid->GetNumberOfPoints());
197 for (Eigen::Index i = 0; i < data.rows(); ++i)
200 Eigen::RowVectorX<T>
const row = data.row(i);
202 array->SetTypedTuple(i, row.data());
206 grid->GetPointData()->AddArray(array);
210void addCellData(vtkUnstructuredGrid* grid, std::string
const& name,
211 Eigen::MatrixX<T>
const& data)
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());
219 if (grid->GetNumberOfCells() != data.rows())
222 "Got {} rows in the table but expected {} rows, same as number of "
223 "cells in the grid.",
224 data.rows(), grid->GetNumberOfCells());
226 for (Eigen::Index i = 0; i < data.rows(); ++i)
229 Eigen::RowVectorX<T>
const row = data.row(i);
231 array->SetTypedTuple(i, row.data());
235 grid->GetCellData()->AddArray(array);
247 INFO(
"converting points");
252 "No point was found in the mesh. Please check whether all anchors "
253 "are within the model region.");
258 "Only one point was found in the mesh. Please check whether all "
259 "anchors are within the model region.");
261 else if (n_pts % 2 != 0)
265 "Number of points is not even. Anchors are believed to consist of "
266 "a start and an end point.");
269 vtkNew<vtkUnstructuredGrid> grid;
272 vtkNew<vtkPoints> points;
273 points->SetNumberOfPoints(n_pts);
276 for (Eigen::Index i = 0; i < css.rows(); ++i)
278 points->SetPoint(i, css(i, 0), css(i, 1), css(i, 2));
282 vtkNew<vtkLine> line;
283 line->GetPointIds()->SetId(0, i - 1);
284 line->GetPointIds()->SetId(1, i);
285 grid->InsertNextCell(line->GetCellType(), line->GetPointIds());
289 grid->SetPoints(points);
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition of the VtkMeshConverter class.
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)
Eigen::VectorXd anchor_cross_sectional_area
Eigen::VectorXd maximum_anchor_stress
Eigen::VectorXd anchor_stiffness
Eigen::VectorX< vtkIdType > bulk_element_ids
Eigen::MatrixXd real_coords
Eigen::VectorXd residual_anchor_stress
Eigen::VectorXd initial_anchor_stress
Eigen::MatrixXd natural_coords
Eigen::VectorX< vtkIdType > point_cloud_node_ids