OGS
ComputeIntersections.cpp File Reference
#include "ComputeIntersections.h"
#include <spdlog/spdlog.h>
#include <vtkCellData.h>
#include <vtkCellLocator.h>
#include <vtkExtractEdges.h>
#include <vtkIdTypeArray.h>
#include <vtkLine.h>
#include "BaseLib/Logging.h"
Include dependency graph for ComputeIntersections.cpp:

Go to the source code of this file.

Functions

int GetGridDimension (vtkUnstructuredGrid *grid)
bool isPointClose (const Eigen::Vector3d &a, const Eigen::Vector3d &b, double tol)
bool compareByT (const IntersectionResult &a, const IntersectionResult &b)
std::vector< IntersectionResultfindOrderedIntersections (vtkUnstructuredGrid *grid, Eigen::Vector3d const &p0, Eigen::Vector3d const &p1, const double free_fraction, double const 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. The line segment is defined by consecutive entries in realcoords.
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.

Function Documentation

◆ compareByT()

bool compareByT ( const IntersectionResult & a,
const IntersectionResult & b )

Definition at line 48 of file ComputeIntersections.cpp.

49{
50 return a.t < b.t;
51}

References IntersectionResult::t.

Referenced by findOrderedIntersections().

◆ findOrderedIntersections()

std::vector< IntersectionResult > findOrderedIntersections ( vtkUnstructuredGrid * grid,
Eigen::Vector3d const & p0,
Eigen::Vector3d const & p1,
const double free_fraction,
double const tol )

Definition at line 53 of file ComputeIntersections.cpp.

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}
bool compareByT(const IntersectionResult &a, const IntersectionResult &b)
bool isPointClose(const Eigen::Vector3d &a, const Eigen::Vector3d &b, double tol)
int GetGridDimension(vtkUnstructuredGrid *grid)
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:36
Result of an intersection of a line with a cell.

References compareByT(), GetGridDimension(), INFO(), isPointClose(), IntersectionResult::point, and IntersectionResult::t.

Referenced by getOrderedAnchorCoords().

◆ GetGridDimension()

int GetGridDimension ( vtkUnstructuredGrid * grid)

Definition at line 21 of file ComputeIntersections.cpp.

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}

Referenced by findOrderedIntersections().

◆ getOrderedAnchorCoords()

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. The line segment is defined by consecutive entries in realcoords.

Parameters
gridThe vtkUnstructuredGrid to intersect with.
realcoordsoriginal coordinates of the anchor start and end points.
free_fractionThe fraction of the line segment to consider for intersections.
tolTolerance for considering points as identical.
Returns
a vector of vectors containing the intersection results which need to be converted using setPhysicalPropertiesForIntersectionPoints() to fill in physical properties and to obtain the final anchor format

Definition at line 136 of file ComputeIntersections.cpp.

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}
std::vector< IntersectionResult > findOrderedIntersections(vtkUnstructuredGrid *grid, Eigen::Vector3d const &p0, Eigen::Vector3d const &p1, const double free_fraction, double const tol)

References findOrderedIntersections().

Referenced by main().

◆ isPointClose()

bool isPointClose ( const Eigen::Vector3d & a,
const Eigen::Vector3d & b,
double tol )

Definition at line 41 of file ComputeIntersections.cpp.

43{
44 return (a - b).squaredNorm() <= tol * tol;
45}

Referenced by findOrderedIntersections().

◆ setPhysicalPropertiesForIntersectionPoints()

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.

Parameters
anchor_coordsThe intersection points of the anchors with the bulk mesh.
original_anchor_dataThe original anchor data read from the JSON file.
Returns
struct containing the anchor data which will be stored in the unstructured grid

Definition at line 155 of file ComputeIntersections.cpp.

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}

References ApplicationUtils::ComputeNaturalCoordsResult::anchor_cross_sectional_area, ApplicationUtils::ComputeNaturalCoordsResult::anchor_stiffness, ApplicationUtils::ComputeNaturalCoordsResult::initial_anchor_stress, ApplicationUtils::ComputeNaturalCoordsResult::maximum_anchor_stress, ApplicationUtils::ComputeNaturalCoordsResult::real_coords, and ApplicationUtils::ComputeNaturalCoordsResult::residual_anchor_stress.

Referenced by main().