@@ -44,6 +44,153 @@ std::size_t Splitter::splitSurfaces(Mesh& mesh) {
4444 return totalNewQuads;
4545}
4646
47+ std::size_t Splitter::splitLines (Mesh& mesh) {
48+ std::size_t totalNewLines = 0 ;
49+
50+ for (GroupId g = 0 ; g < mesh.groups .size (); g++) {
51+ std::vector<Element> newElements;
52+
53+ for (ElementId e = 0 ; e < mesh.groups [g].elements .size (); e++) {
54+ const Element& elem = mesh.groups [g].elements [e];
55+ if (elem.type == Element::Type::Line) {
56+ // Split this line into unit grid lines
57+ std::map<Coordinate, CoordinateId> coordMap;
58+
59+ // Build initial coord map from existing coordinates
60+ for (CoordinateId i = 0 ; i < static_cast <CoordinateId>(mesh.coordinates .size ()); ++i) {
61+ coordMap[mesh.coordinates [i]] = i;
62+ }
63+
64+ std::vector<Element> splitLines = splitLine_ (
65+ elem, mesh.coordinates , mesh.grid , coordMap);
66+
67+ // Add new coordinates from coordMap
68+ for (const auto & [coord, id] : coordMap) {
69+ if (id >= mesh.coordinates .size ()) {
70+ mesh.coordinates .push_back (coord);
71+ }
72+ }
73+
74+ newElements.insert (newElements.end (), splitLines.begin (), splitLines.end ());
75+ totalNewLines += splitLines.size ();
76+ } else {
77+ newElements.push_back (elem);
78+ }
79+ }
80+
81+ mesh.groups [g].elements = std::move (newElements);
82+ }
83+
84+ return totalNewLines;
85+ }
86+
87+ std::vector<Element> Splitter::splitLine_ (
88+ const Element& line,
89+ const std::vector<Coordinate>& coords,
90+ const Grid& grid,
91+ std::map<Coordinate, CoordinateId>& coordMap) {
92+ std::vector<Element> unitLines;
93+
94+ // Get the axis direction of the line
95+ Axis lineAxis = getLineAxis_ (line, coords);
96+
97+ // Get the grid cell bounds
98+ auto [minCell, maxCell] = getLineBounds_ (line, coords);
99+
100+ // Determine the fixed coordinate axes (the two axes perpendicular to lineAxis)
101+ Axis axis1 = (lineAxis + 1 ) % 3 ;
102+ Axis axis2 = (lineAxis + 2 ) % 3 ;
103+
104+ // Get the fixed coordinate values from minCell
105+ CellDir fixedCoord1 = minCell (axis1);
106+ CellDir fixedCoord2 = minCell (axis2);
107+
108+ // Generate unit lines for each cell along the line axis
109+ for (CellDir i = minCell (lineAxis); i < maxCell (lineAxis); i++) {
110+ Element unitLine = createUnitLine_ (
111+ fixedCoord1, fixedCoord2, lineAxis, i, grid, coordMap);
112+ unitLines.push_back (unitLine);
113+ }
114+
115+ return unitLines;
116+ }
117+
118+ std::pair<Cell, Cell> Splitter::getLineBounds_ (
119+ const Element& line,
120+ const std::vector<Coordinate>& coords) {
121+ Cell minCell = utils::GridTools::toCell (coords[line.vertices [0 ]]);
122+ Cell maxCell = utils::GridTools::toCell (coords[line.vertices [1 ]]);
123+
124+ // Ensure minCell <= maxCell for all axes
125+ for (Axis d = 0 ; d < 3 ; d++) {
126+ if (minCell (d) > maxCell (d)) {
127+ std::swap (minCell (d), maxCell (d));
128+ }
129+ }
130+
131+ return {minCell, maxCell};
132+ }
133+
134+ Axis Splitter::getLineAxis_ (
135+ const Element& line,
136+ const std::vector<Coordinate>& coords) {
137+ // Find which axis has different coordinate values (the line direction)
138+ Cell cell0 = utils::GridTools::toCell (coords[line.vertices [0 ]]);
139+ Cell cell1 = utils::GridTools::toCell (coords[line.vertices [1 ]]);
140+
141+ for (Axis d = 0 ; d < 3 ; d++) {
142+ if (cell0 (d) != cell1 (d)) {
143+ return d;
144+ }
145+ }
146+
147+ // Fallback (should not happen for valid lines)
148+ return 0 ;
149+ }
150+
151+ Element Splitter::createUnitLine_ (
152+ CellDir fixedCoord1,
153+ CellDir fixedCoord2,
154+ Axis lineAxis,
155+ CellDir cell,
156+ const Grid& grid,
157+ std::map<Coordinate, CoordinateId>& coordMap) {
158+ Axis axis1 = (lineAxis + 1 ) % 3 ;
159+ Axis axis2 = (lineAxis + 2 ) % 3 ;
160+
161+ // Create 2 coordinates for the unit line
162+ std::array<Coordinate, 2 > endpoints;
163+ endpoints[0 ] = Coordinate ({0 , 0 , 0 });
164+ endpoints[1 ] = Coordinate ({0 , 0 , 0 });
165+
166+ // Set coordinates for each endpoint
167+ endpoints[0 ](lineAxis) = grid[lineAxis][cell];
168+ endpoints[0 ](axis1) = grid[axis1][fixedCoord1];
169+ endpoints[0 ](axis2) = grid[axis2][fixedCoord2];
170+
171+ endpoints[1 ](lineAxis) = grid[lineAxis][cell + 1 ];
172+ endpoints[1 ](axis1) = grid[axis1][fixedCoord1];
173+ endpoints[1 ](axis2) = grid[axis2][fixedCoord2];
174+
175+ // Get or create coordinate IDs
176+ std::array<CoordinateId, 2 > vids;
177+ for (int i = 0 ; i < 2 ; i++) {
178+ auto it = coordMap.find (endpoints[i]);
179+ if (it != coordMap.end ()) {
180+ vids[i] = it->second ;
181+ } else {
182+ coordMap[endpoints[i]] = static_cast <CoordinateId>(coordMap.size ());
183+ vids[i] = coordMap[endpoints[i]];
184+ }
185+ }
186+
187+ Element unitLine;
188+ unitLine.type = Element::Type::Line;
189+ unitLine.vertices = {vids[0 ], vids[1 ]};
190+
191+ return unitLine;
192+ }
193+
47194std::vector<Element> Splitter::splitSurface_ (
48195 const Element& surface,
49196 const std::vector<Coordinate>& coords,
0 commit comments