OGS
ComputeIntersections.cpp
Go to the documentation of this file.
1
9
11
12#include <spdlog/spdlog.h>
13#include <vtkCellData.h>
14#include <vtkCellLocator.h>
15#include <vtkExtractEdges.h>
16#include <vtkIdTypeArray.h>
17#include <vtkLine.h>
18
19#include "BaseLib/Logging.h"
20
21int GetGridDimension(vtkUnstructuredGrid* grid)
22{
23 vtkIdType numCells = grid->GetNumberOfCells();
24 int maxDim = 0;
25
26 vtkSmartPointer<vtkGenericCell> cell =
27 vtkSmartPointer<vtkGenericCell>::New();
28
29 for (vtkIdType i = 0; i < numCells; ++i)
30 {
31 grid->GetCell(i, cell);
32 int dim = cell->GetCellDimension();
33 if (dim > maxDim)
34 {
35 maxDim = dim;
36 }
37 }
38 return maxDim;
39}
40
41bool isPointClose(const Eigen::Vector3d& a, const Eigen::Vector3d& b,
42 double tol)
43{
44 return (a - b).squaredNorm() <= tol * tol;
45}
46
47// Comparator to sort intersections by t value
49{
50 return a.t < b.t;
51}
52
53std::vector<IntersectionResult> findOrderedIntersections(
54 vtkUnstructuredGrid* grid,
55 Eigen::Vector3d const& p0,
56 Eigen::Vector3d const& p1,
57 const double free_fraction,
58 double const tol)
59{
60 std::vector<IntersectionResult> intersections;
61
62 intersections.push_back(IntersectionResult(p0, 0.0));
63 intersections.push_back(IntersectionResult(p1, 1.0));
64 Eigen::Vector3d transition_point = p0 + (p1 - p0) * free_fraction;
65 if (!isPointClose(p0, transition_point, tol) &&
66 !isPointClose(p1, transition_point, tol))
67 {
68 INFO("Adding transition point at t = {:.2f} for free fraction {:.2f}",
69 free_fraction, free_fraction);
70 intersections.push_back(
71 IntersectionResult(transition_point, free_fraction));
72 }
73 // For cell intersection
74 vtkSmartPointer<vtkGenericCell> cell =
75 vtkSmartPointer<vtkGenericCell>::New();
76
77 vtkNew<vtkIdList> filteredCellIds;
78 vtkNew<vtkCellLocator> locator;
79 vtkSmartPointer<vtkDataSet> cells_to_be_intersected;
80 if (GetGridDimension(grid) < 3)
81 {
82 vtkSmartPointer<vtkExtractEdges> extractEdges =
83 vtkSmartPointer<vtkExtractEdges>::New();
84 extractEdges->SetInputData(grid);
85 extractEdges->Update();
86
87 cells_to_be_intersected = extractEdges->GetOutput();
88 }
89 else
90 {
91 cells_to_be_intersected = grid;
92 }
93 locator->SetDataSet(cells_to_be_intersected);
94 locator->BuildLocator();
95 locator->FindCellsAlongLine(p0.data(), p1.data(), tol, filteredCellIds);
96 for (vtkIdType i = 0; i < filteredCellIds->GetNumberOfIds(); ++i)
97 {
98 vtkIdType cellId = filteredCellIds->GetId(i);
99 cells_to_be_intersected->GetCell(cellId, cell);
100
101 double t;
102 double x[3];
103 double pcoords[3];
104 int subId;
105
106 int intersected = cell->IntersectWithLine(p0.data(), p1.data(), tol, t,
107 x, pcoords, subId);
108
109 Eigen::Vector3d x_vec(x[0], x[1], x[2]);
110
111 if (!intersected)
112 {
113 continue;
114 }
115 bool isDuplicate = false;
116 for (const auto& existing : intersections)
117 {
118 if (isPointClose(existing.point, x_vec, tol))
119 {
120 isDuplicate = true;
121 break;
122 }
123 }
124 if (!isDuplicate && (t < tol || t > free_fraction))
125 {
126 IntersectionResult result;
127 result.t = t;
128 result.point = x_vec;
129 intersections.push_back(result);
130 }
131 }
132 std::sort(intersections.begin(), intersections.end(), compareByT);
133 return intersections;
134}
135
136std::vector<std::vector<IntersectionResult>> getOrderedAnchorCoords(
137 vtkUnstructuredGrid* grid,
138 Eigen::MatrixX3d const& realcoords,
139 Eigen::VectorXd const& free_fraction,
140 double const tol)
141{
142 std::vector<std::vector<IntersectionResult>> ordered_intersections(
143 realcoords.rows() / 2);
144 assert(realcoords.rows() % 2 == 0);
145 assert(free_fraction.size() == realcoords.rows() / 2);
146 for (int i = 0; i < realcoords.rows(); i += 2)
147 {
148 ordered_intersections[i / 2] = findOrderedIntersections(
149 grid, realcoords.row(i), realcoords.row(i + 1),
150 free_fraction(i / 2), tol);
151 }
152 return ordered_intersections;
153}
154
156 std::vector<std::vector<IntersectionResult>> const& anchor_coords,
157 AU::ComputeNaturalCoordsResult const& original_anchor_data)
158{
159 std::vector<IntersectionResult> all_intersections;
160 std::vector<int> anchorids;
161 for (auto&& [anchor_id, intersections] :
162 ranges::views::enumerate(anchor_coords))
163 {
164 for (auto&& [index, result] : ranges::views::enumerate(intersections))
165 {
166 if (index > 0 && index < intersections.size() - 1)
167 {
168 anchorids.push_back(anchor_id);
169 all_intersections.push_back(result);
170 }
171 anchorids.push_back(anchor_id);
172 all_intersections.push_back(result);
173 }
174 }
175 assert(all_intersections.size() == anchorids.size());
176 assert(all_intersections.size() % 2 == 0);
177 Eigen::MatrixXd realcoords(anchorids.size(), 3);
178 auto const number_of_anchors = anchorids.size() / 2;
179 Eigen::VectorXd initial_stress(number_of_anchors);
180 Eigen::VectorXd maximum_stress(number_of_anchors);
181 Eigen::VectorXd residual_stress(number_of_anchors);
182 Eigen::VectorXd cross_sectional_area(number_of_anchors);
183 Eigen::VectorXd stiffness(number_of_anchors);
184 for (std::size_t i = 0; i < all_intersections.size(); ++i)
185 {
186 realcoords.row(i).noalias() = all_intersections[i].point;
187 if (i % 2 == 0)
188 {
189 initial_stress(i / 2) =
190 original_anchor_data.initial_anchor_stress(anchorids[i]);
191 maximum_stress(i / 2) =
192 original_anchor_data.maximum_anchor_stress(anchorids[i]);
193 residual_stress(i / 2) =
194 original_anchor_data.residual_anchor_stress(anchorids[i]);
195 cross_sectional_area(i / 2) =
196 original_anchor_data.anchor_cross_sectional_area(anchorids[i]);
197 stiffness(i / 2) =
198 original_anchor_data.anchor_stiffness(anchorids[i]);
199 }
200 }
202 anchor_data.real_coords = realcoords;
203 anchor_data.initial_anchor_stress = initial_stress;
204 anchor_data.maximum_anchor_stress = maximum_stress;
205 anchor_data.residual_anchor_stress = residual_stress;
206 anchor_data.anchor_cross_sectional_area = cross_sectional_area;
207 anchor_data.anchor_stiffness = stiffness;
208 return anchor_data;
209}
bool compareByT(const IntersectionResult &a, const IntersectionResult &b)
bool isPointClose(const Eigen::Vector3d &a, const Eigen::Vector3d &b, double tol)
std::vector< std::vector< IntersectionResult > > getOrderedAnchorCoords(vtkUnstructuredGrid *grid, Eigen::MatrixX3d const &realcoords, Eigen::VectorXd const &free_fraction, double const tol)
Finds intersection points of a line segment with the cells of a vtkUnstructuredGrid....
int GetGridDimension(vtkUnstructuredGrid *grid)
AU::ComputeNaturalCoordsResult setPhysicalPropertiesForIntersectionPoints(std::vector< std::vector< IntersectionResult > > const &anchor_coords, AU::ComputeNaturalCoordsResult const &original_anchor_data)
fills the physical properties of the intersection points based on the original anchor data.
std::vector< IntersectionResult > findOrderedIntersections(vtkUnstructuredGrid *grid, Eigen::Vector3d const &p0, Eigen::Vector3d const &p1, const double free_fraction, double const tol)
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:36
Result of an intersection of a line with a cell.