1 // Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC and
2 // other Axom Project Developers. See the top-level LICENSE file for details.
3 //
4 // SPDX-License-Identifier: (BSD-3-Clause)
5
6 #include "axom/mint/mesh/internal/MeshHelpers.hpp"
7 #include "axom/mint/mesh/Mesh.hpp"
8 #include "axom/core/utilities/Utilities.hpp"
9
10 #include <algorithm>
11 #include <unordered_map>
12 #include <vector>
13
14 namespace axom
15 {
16 namespace mint
17 {
18 namespace internal
19 {
20 //------------------------------------------------------------------------------
21 // IMPLEMENTATION OF FREE METHODS
22 //------------------------------------------------------------------------------
join_ints_into_string(int count,IndexType * values,char sep)23 std::string join_ints_into_string(int count, IndexType* values, char sep)
24 {
25 std::stringstream joined;
26
27 for(int i = 0; i < count; ++i)
28 {
29 if(i > 0)
30 {
31 joined << sep;
32 }
33 joined << values[i];
34 }
35
36 return joined.str();
37 }
38
39 //------------------------------------------------------------------------------
make_face_key(int count,IndexType * values,char sep)40 std::string make_face_key(int count, IndexType* values, char sep)
41 {
42 std::vector<IndexType> locvalues(values, values + count);
43 std::sort(locvalues.begin(), locvalues.end());
44 return join_ints_into_string(count, locvalues.data(), sep);
45 }
46
47 //------------------------------------------------------------------------------
otherSide(IndexType * f2c,IndexType thisSide)48 IndexType otherSide(IndexType* f2c, IndexType thisSide)
49 {
50 return (f2c[0] != thisSide) ? f2c[0] : f2c[1];
51 }
52
53 //------------------------------------------------------------------------------
54 struct FaceTypeCellsNodes
55 {
FaceTypeCellsNodesaxom::mint::internal::FaceTypeCellsNodes56 FaceTypeCellsNodes() : facetype(UNDEFINED_CELL) { }
57
FaceTypeCellsNodesaxom::mint::internal::FaceTypeCellsNodes58 FaceTypeCellsNodes(CellType ftype,
59 std::vector<IndexType>& fcells,
60 std::vector<IndexType>& fnodes)
61 : facetype(ftype)
62 , facecells(fcells)
63 , facenodes(fnodes)
64 { }
65
66 CellType facetype;
67 std::vector<IndexType> facecells;
68 std::vector<IndexType> facenodes;
69 };
70
71 //------------------------------------------------------------------------------
initFaces(Mesh * mesh,IndexType & facecount,IndexType * & f2c,IndexType * & c2f,IndexType * & c2n,IndexType * & c2foffsets,IndexType * & f2n,IndexType * & f2noffsets,CellType * & f2ntypes)72 bool initFaces(Mesh* mesh,
73 IndexType& facecount,
74 IndexType*& f2c,
75 IndexType*& c2f,
76 IndexType*& c2n,
77 IndexType*& c2foffsets,
78 IndexType*& f2n,
79 IndexType*& f2noffsets,
80 CellType*& f2ntypes)
81 {
82 bool success = true;
83
84 facecount = 0;
85 f2c = nullptr;
86 c2f = nullptr;
87 c2n = nullptr;
88 c2foffsets = nullptr;
89 f2n = nullptr;
90 f2noffsets = nullptr;
91 f2ntypes = nullptr;
92
93 using IDtoKeyType = std::unordered_map<IndexType, std::vector<IndexType>>;
94 using FaceBuilderType = std::map<std::vector<IndexType>,
95 FaceTypeCellsNodes,
96 utilities::LexiComparator<IndexType>>;
97 FaceBuilderType workface;
98
99 // Iterate over each cell.
100 const IndexType cellcount = mesh->getNumberOfCells();
101 for(int c = 0; c < cellcount; ++c)
102 {
103 // Step 1. For every cell, get the nodes.
104 IndexType nodes[MAX_CELL_NODES];
105 mesh->getCellNodeIDs(c, nodes);
106 const CellType celltype = mesh->getCellType(c);
107 const CellInfo thisCell = getCellInfo(celltype);
108
109 int base = 0;
110 for(int f = 0; f < thisCell.num_faces; ++f)
111 {
112 // Step 2. The cell nodes will be in "VTK Order," as specified by
113 // https://www.vtk.org/wp-content/uploads/2015/04/file-formats.pdf.
114 // This routine selects the face nodes as listed in the registered
115 // cell info to make sure that face normals point outward.
116 //
117 // For every face, sort the node IDs and join them together in a
118 // string with '.' delimiters. This is the key for an associative
119 // array whose value is a list of cell IDs.
120 int num_face_nodes = thisCell.face_nodecount[f];
121 std::vector<IndexType> inorder_facenodes(num_face_nodes);
122 const IndexType* face_node_offsets = thisCell.face_nodes + base;
123 for(int fn = 0; fn < num_face_nodes; ++fn)
124 {
125 inorder_facenodes[fn] = nodes[face_node_offsets[fn]];
126 }
127 base += num_face_nodes;
128
129 std::vector<IndexType> sorted_facenodes(inorder_facenodes);
130 std::sort(sorted_facenodes.begin(), sorted_facenodes.end());
131
132 FaceTypeCellsNodes cells_of_face;
133 if(workface.count(sorted_facenodes) > 0)
134 {
135 cells_of_face = workface[sorted_facenodes];
136 }
137 else
138 {
139 cells_of_face.facetype = thisCell.face_types[f];
140 cells_of_face.facenodes = inorder_facenodes;
141 }
142 cells_of_face.facecells.push_back(c);
143 workface[sorted_facenodes] = cells_of_face;
144 }
145 }
146
147 // Step 3. For each face, record its incident cells.
148 // Here we use ConnectivityArray::reserve() and then append() each
149 // face. This means the kth inserted face gets ID k.
150 IndexType faceID = 0;
151 IndexType faceNodeTotal = 0;
152 IDtoKeyType keys;
153 f2c = new IndexType[2 * workface.size()];
154 for(FaceBuilderType::value_type v : workface)
155 {
156 FaceTypeCellsNodes theFace = v.second;
157 faceNodeTotal += static_cast<IndexType>(theFace.facenodes.size());
158 int faceCellCount = static_cast<int>(theFace.facecells.size());
159 IndexType* faceCells = theFace.facecells.data();
160
161 if(faceCellCount < 1 || faceCellCount > 2)
162 {
163 success = false;
164 }
165 else
166 {
167 f2c[2 * faceID] = faceCells[0];
168 f2c[2 * faceID + 1] = (faceCellCount > 1 ? faceCells[1] : -1);
169 }
170
171 keys[faceID] = v.first;
172
173 faceID += 1;
174 }
175
176 // If we have any face with less than one or more than two incident cells,
177 // clean up and return failure. We won't do any more work here.
178 if(!success)
179 {
180 delete[] f2c;
181 f2c = nullptr;
182 return success;
183 }
184
185 // Record how many faces we have in this mesh.
186 facecount = faceID;
187
188 // Now that we have a count of all the face-nodes, and we know we have no
189 // faces with more than two nodes, record the face-node relations.
190 f2n = new IndexType[faceNodeTotal];
191 f2noffsets = new IndexType[facecount + 1];
192 f2ntypes = new CellType[facecount];
193 int faceNodeOffset = 0;
194 for(int fidx = 0; fidx < facecount; ++fidx)
195 {
196 FaceTypeCellsNodes theFace = workface[keys[fidx]];
197 std::vector<IndexType>& faceNodes = theFace.facenodes;
198 std::copy(faceNodes.begin(), faceNodes.end(), f2n + faceNodeOffset);
199 f2noffsets[fidx] = faceNodeOffset;
200 f2ntypes[fidx] = theFace.facetype;
201 faceNodeOffset += static_cast<int>(faceNodes.size());
202 }
203 f2noffsets[facecount] = faceNodeOffset;
204
205 // Step 4. Now that we have face IDs, record cell-to-face relation.
206 typedef std::unordered_map<IndexType, std::vector<IndexType>> CellFaceBuilderType;
207
208 CellFaceBuilderType cell_to_face;
209 int cellFaceCount = 0;
210 for(int f = 0; f < facecount; ++f)
211 {
212 for(IndexType c : workface[keys[f]].facecells)
213 {
214 std::vector<IndexType> cell_faces;
215 if(cell_to_face.count(c) > 0)
216 {
217 cell_faces = cell_to_face[c];
218 }
219 cellFaceCount += 1;
220 cell_faces.push_back(f);
221 cell_to_face[c] = cell_faces;
222 }
223 }
224
225 // Step 4b. Put cell-to-face relation into output arrays.
226 c2f = new IndexType[cellFaceCount];
227 c2n = new IndexType[cellFaceCount];
228 c2foffsets = new IndexType[cellcount + 1];
229 cellFaceCount = 0;
230
231 for(IndexType cellID = 0; cellID < cellcount; ++cellID)
232 {
233 // For the current cell ID, set its type and value offset
234 // cellTypes[cellID] = m_cell_connectivity->getIDType(cellID);
235 c2foffsets[cellID] = cellFaceCount;
236
237 if(cell_to_face.count(cellID) > 0)
238 {
239 // If we have faces for the current cell (which we certainly
240 // *SHOULD* have), copy them in.
241 std::vector<IndexType>& theCells = cell_to_face[cellID];
242 IndexType* faceIDs = theCells.data();
243 int thisCellFaceCount = static_cast<int>(theCells.size());
244
245 // Copy in values
246 for(int f = 0; f < thisCellFaceCount; ++f)
247 {
248 c2f[cellFaceCount + f] = faceIDs[f];
249 c2n[cellFaceCount + f] = otherSide(&f2c[2 * faceIDs[f]], cellID);
250 }
251
252 // Maintain offset
253 cellFaceCount += thisCellFaceCount;
254 }
255 }
256 c2foffsets[cellcount] = cellFaceCount;
257
258 return success;
259 }
260
261 } /* namespace internal */
262 } /* namespace mint */
263 } /* namespace axom */
264