OGS
MeshUtils.cpp
Go to the documentation of this file.
1
10#include "MeshUtils.h"
11
12#include "BaseLib/Algorithm.h"
14#include "MeshLib/Mesh.h"
16#include "MeshLib/Node.h"
17
18namespace ProcessLib
19{
20namespace LIE
21{
22namespace
23{
24// A class to check whether a node is located on a crack tip with
25// the following conditions:
26// - the number of connected fracture elements is one
27// - the node is not located on a domain boundary
29{
30public:
31 explicit IsCrackTip(MeshLib::Mesh const& mesh)
32 : _mesh(mesh), _fracture_element_dim(_mesh.getDimension() - 1)
33 {
34 _is_internal_node.resize(_mesh.getNumberOfNodes(), true);
35
36 MeshLib::NodeSearch nodeSearch(mesh);
37 nodeSearch.searchBoundaryNodes();
38 for (auto i : nodeSearch.getSearchedNodeIDs())
39 {
40 _is_internal_node[i] = false;
41 }
42 }
43
44 bool operator()(MeshLib::Node const& node) const
45 {
46 if (!_is_internal_node[node.getID()] ||
47 !isBaseNode(node, _mesh.getElementsConnectedToNode(node)))
48 {
49 return false;
50 }
51
52 auto const elements_connected_to_node =
53 _mesh.getElementsConnectedToNode(node);
54 auto const n_connected_fracture_elements =
55 count_if(cbegin(elements_connected_to_node),
56 cend(elements_connected_to_node),
57 [&](auto* const e)
58 { return e->getDimension() == _fracture_element_dim; });
59 assert(n_connected_fracture_elements > 0);
60
61 return (n_connected_fracture_elements == 1);
62 }
63
64private:
66 unsigned const _fracture_element_dim;
67 std::vector<bool> _is_internal_node;
68};
69
71 MeshLib::Mesh const& mesh, std::vector<int> const& vec_fracture_mat_IDs,
72 std::vector<std::vector<MeshLib::Node*>> const& vec_fracture_nodes,
73 std::vector<std::vector<MeshLib::Element*>>& intersected_fracture_elements,
74 std::vector<std::pair<std::size_t, std::vector<int>>>&
75 vec_branch_nodeID_matIDs,
76 std::vector<std::pair<std::size_t, std::vector<int>>>&
77 vec_junction_nodeID_matIDs)
78{
79 auto const n_fractures = vec_fracture_mat_IDs.size();
80 std::map<unsigned, unsigned> matID_to_fid;
81 for (unsigned i = 0; i < n_fractures; i++)
82 {
83 matID_to_fid[vec_fracture_mat_IDs[i]] = i;
84 }
85
86 // make a vector all fracture nodes
87 std::vector<std::size_t> all_fracture_nodes;
88 for (auto& vec : vec_fracture_nodes)
89 {
90 transform(cbegin(vec), cend(vec), back_inserter(all_fracture_nodes),
91 [](auto* const n) { return n->getID(); });
92 }
93
94 // create a table of a node id and connected material IDs
95 std::map<std::size_t, std::vector<std::size_t>> frac_nodeID_to_matIDs;
96 for (unsigned i = 0; i < n_fractures; i++)
97 {
98 for (auto* node : vec_fracture_nodes[i])
99 {
100 frac_nodeID_to_matIDs[node->getID()].push_back(
101 vec_fracture_mat_IDs[i]);
102 }
103 }
104
105 auto const* const opt_material_ids = MeshLib::materialIDs(mesh);
106
107 // find branch/junction nodes which connect to multiple fractures
108 intersected_fracture_elements.resize(n_fractures);
109 for (auto frac_nodeID_to_matID : frac_nodeID_to_matIDs)
110 {
111 auto nodeID = frac_nodeID_to_matID.first;
112 auto const* node = mesh.getNode(frac_nodeID_to_matID.first);
113 auto const& matIDs = frac_nodeID_to_matID.second;
114 if (matIDs.size() < 2)
115 {
116 continue; // no intersection
117 }
118
119 auto const elements_connected_to_node =
120 mesh.getElementsConnectedToNode(*node);
121 std::vector<MeshLib::Element*> conn_fracture_elements;
122 {
123 for (auto const* e : elements_connected_to_node)
124 {
125 if (e->getDimension() == (mesh.getDimension() - 1))
126 {
127 conn_fracture_elements.push_back(
128 const_cast<MeshLib::Element*>(e));
129 }
130 }
131 }
132
133 std::map<int, int> vec_matID_counts;
134 {
135 for (auto matid : matIDs)
136 {
137 vec_matID_counts[matid] = 0;
138 }
139
140 for (auto const* e : conn_fracture_elements)
141 {
142 auto matid = (*opt_material_ids)[e->getID()];
143 vec_matID_counts[matid]++;
144 }
145 }
146
147 for (auto matid : matIDs)
148 {
149 auto fid = matID_to_fid[matid];
150 for (auto* e : conn_fracture_elements)
151 {
152 auto e_matid = (*opt_material_ids)[e->getID()];
153 if (matID_to_fid[e_matid] != fid)
154 { // only slave elements
155 intersected_fracture_elements[fid].push_back(e);
156 }
157 }
158 }
159
160 bool isBranch = false;
161 {
162 for (auto vec_matID_count : vec_matID_counts)
163 {
164 auto count = vec_matID_count.second;
165 if (count % 2 == 1)
166 {
167 isBranch = true;
168 break;
169 }
170 }
171 }
172
173 if (isBranch)
174 {
175 std::vector<int> branch_matIDs(2);
176 for (auto vec_matID_count : vec_matID_counts)
177 {
178 auto matid = vec_matID_count.first;
179 auto count = vec_matID_count.second;
180 if (count % 2 == 0)
181 {
182 branch_matIDs[0] = matid; // master
183 }
184 else
185 {
186 branch_matIDs[1] = matid; // slave
187 }
188 }
189 vec_branch_nodeID_matIDs.emplace_back(nodeID, branch_matIDs);
190 }
191 else
192 {
193 std::vector<int> junction_matIDs(2);
194 junction_matIDs[0] = std::min(vec_matID_counts.begin()->first,
195 vec_matID_counts.rbegin()->first);
196 junction_matIDs[1] = std::max(vec_matID_counts.begin()->first,
197 vec_matID_counts.rbegin()->first);
198 vec_junction_nodeID_matIDs.emplace_back(nodeID, junction_matIDs);
199 }
200 }
201
202 for (auto& elements : intersected_fracture_elements)
203 {
205 }
206
207 DBUG("-> found {:d} branches and {:d} junctions",
208 vec_branch_nodeID_matIDs.size(), vec_junction_nodeID_matIDs.size());
209}
210
211} // namespace
212
214 MeshLib::Mesh const& mesh,
215 std::vector<MeshLib::Element*>& vec_matrix_elements,
216 std::vector<int>& vec_fracture_mat_IDs,
217 std::vector<std::vector<MeshLib::Element*>>& vec_fracture_elements,
218 std::vector<std::vector<MeshLib::Element*>>& vec_fracture_matrix_elements,
219 std::vector<std::vector<MeshLib::Node*>>& vec_fracture_nodes,
220 std::vector<std::pair<std::size_t, std::vector<int>>>&
221 vec_branch_nodeID_matIDs,
222 std::vector<std::pair<std::size_t, std::vector<int>>>&
223 vec_junction_nodeID_matIDs)
224{
225 IsCrackTip isCrackTip(mesh);
226
227 // get vectors of matrix elements and fracture elements
228 vec_matrix_elements.reserve(mesh.getNumberOfElements());
229 std::vector<MeshLib::Element*> all_fracture_elements;
230 for (MeshLib::Element* e : mesh.getElements())
231 {
232 if (e->getDimension() == mesh.getDimension())
233 {
234 vec_matrix_elements.push_back(e);
235 }
236 else
237 {
238 all_fracture_elements.push_back(e);
239 }
240 }
241 DBUG("-> found total {:d} matrix elements and {:d} fracture elements",
242 vec_matrix_elements.size(), all_fracture_elements.size());
243
244 // get fracture material IDs
245 auto const material_ids = materialIDs(mesh);
246 if (!material_ids)
247 {
248 OGS_FATAL("Could not access MaterialIDs property from mesh.");
249 }
250 transform(cbegin(all_fracture_elements), cend(all_fracture_elements),
251 back_inserter(vec_fracture_mat_IDs),
252 [&material_ids](auto const* e)
253 { return (*material_ids)[e->getID()]; });
254
255 BaseLib::makeVectorUnique(vec_fracture_mat_IDs);
256 DBUG("-> found {:d} fracture material groups", vec_fracture_mat_IDs.size());
257
258 // create a vector of fracture elements for each material
259 vec_fracture_elements.resize(vec_fracture_mat_IDs.size());
260 for (unsigned frac_id = 0; frac_id < vec_fracture_mat_IDs.size(); frac_id++)
261 {
262 const auto frac_mat_id = vec_fracture_mat_IDs[frac_id];
263 std::vector<MeshLib::Element*>& vec_elements =
264 vec_fracture_elements[frac_id];
265 std::copy_if(all_fracture_elements.begin(), all_fracture_elements.end(),
266 std::back_inserter(vec_elements),
267 [&](MeshLib::Element const* const e)
268 { return (*material_ids)[e->getID()] == frac_mat_id; });
269 DBUG("-> found {:d} elements on the fracture {:d}", vec_elements.size(),
270 frac_id);
271 }
272
273 // get a vector of fracture nodes for each material
274 vec_fracture_nodes.resize(vec_fracture_mat_IDs.size());
275 for (unsigned frac_id = 0; frac_id < vec_fracture_mat_IDs.size(); frac_id++)
276 {
277 std::vector<MeshLib::Node*>& vec_nodes = vec_fracture_nodes[frac_id];
278 for (MeshLib::Element const* const e : vec_fracture_elements[frac_id])
279 {
280 for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
281 {
282 if (isCrackTip(*e->getNode(i)))
283 {
284 continue;
285 }
286 vec_nodes.push_back(const_cast<MeshLib::Node*>(e->getNode(i)));
287 }
288 }
291 DBUG("-> found {:d} nodes on the fracture {:d}", vec_nodes.size(),
292 frac_id);
293 }
294
295 // find branch/junction nodes which connect to multiple fractures
296 std::vector<std::vector<MeshLib::Element*>> intersected_fracture_elements;
297 findFracutreIntersections(mesh, vec_fracture_mat_IDs, vec_fracture_nodes,
298 intersected_fracture_elements,
299 vec_branch_nodeID_matIDs,
300 vec_junction_nodeID_matIDs);
301
302 // create a vector fracture elements and connected matrix elements,
303 // which are passed to a DoF table
304 for (unsigned fid = 0; fid < vec_fracture_elements.size(); fid++)
305 {
306 auto const& fracture_elements = vec_fracture_elements[fid];
307 std::vector<MeshLib::Element*> vec_ele;
308 // first, collect matrix elements
309 for (MeshLib::Element const* const e : fracture_elements)
310 {
311 // it is sufficient to iterate over base nodes, because they are
312 // already connected to all neighbours
313 for (unsigned i = 0; i < e->getNumberOfBaseNodes(); i++)
314 {
315 MeshLib::Node const* node = e->getNode(i);
316 if (isCrackTip(*node))
317 {
318 continue;
319 }
320 auto const& elements_connected_to_node =
321 mesh.getElementsConnectedToNode(*node);
322 for (unsigned j = 0; j < elements_connected_to_node.size(); j++)
323 {
324 // only matrix elements
325 if (elements_connected_to_node[j]->getDimension() <
326 mesh.getDimension())
327 {
328 continue;
329 }
330 vec_ele.push_back(const_cast<MeshLib::Element*>(
331 elements_connected_to_node[j]));
332 }
333 }
334 }
337
338 // second, append fracture elements
339 std::copy(fracture_elements.begin(), fracture_elements.end(),
340 std::back_inserter(vec_ele));
341 // thirdly, append intersected fracture elements
342 std::copy(intersected_fracture_elements[fid].begin(),
343 intersected_fracture_elements[fid].end(),
344 std::back_inserter(vec_ele));
345
346 vec_fracture_matrix_elements.push_back(vec_ele);
347 }
348}
349
350} // namespace LIE
351} // namespace ProcessLib
Definition of the Element class.
#define OGS_FATAL(...)
Definition Error.h:26
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Definition of the Mesh class.
Definition of the Node class.
std::size_t getID() const
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:109
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition Mesh.h:88
const Node * getNode(std::size_t idx) const
Get the node with the given index.
Definition Mesh.h:91
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition Mesh.cpp:256
std::size_t getNumberOfElements() const
Get the number of elements.
Definition Mesh.h:97
Node search class.
Definition NodeSearch.h:25
const std::vector< std::size_t > & getSearchedNodeIDs() const
return marked node IDs
Definition NodeSearch.h:30
std::size_t searchBoundaryNodes()
Marks all boundary nodes.
void makeVectorUnique(std::vector< T > &v)
Definition Algorithm.h:176
bool idsComparator(T const a, T const b)
Definition Mesh.h:206
PropertyVector< int > const * materialIDs(Mesh const &mesh)
Definition Mesh.cpp:268
void findFracutreIntersections(MeshLib::Mesh const &mesh, std::vector< int > const &vec_fracture_mat_IDs, std::vector< std::vector< MeshLib::Node * > > const &vec_fracture_nodes, std::vector< std::vector< MeshLib::Element * > > &intersected_fracture_elements, std::vector< std::pair< std::size_t, std::vector< int > > > &vec_branch_nodeID_matIDs, std::vector< std::pair< std::size_t, std::vector< int > > > &vec_junction_nodeID_matIDs)
Definition MeshUtils.cpp:70
void getFractureMatrixDataInMesh(MeshLib::Mesh const &mesh, std::vector< MeshLib::Element * > &vec_matrix_elements, std::vector< int > &vec_fracture_mat_IDs, std::vector< std::vector< MeshLib::Element * > > &vec_fracture_elements, std::vector< std::vector< MeshLib::Element * > > &vec_fracture_matrix_elements, std::vector< std::vector< MeshLib::Node * > > &vec_fracture_nodes, std::vector< std::pair< std::size_t, std::vector< int > > > &vec_branch_nodeID_matIDs, std::vector< std::pair< std::size_t, std::vector< int > > > &vec_junction_nodeID_matIDs)