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 
18 namespace ProcessLib
19 {
20 namespace LIE
21 {
22 namespace
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 {
30 public:
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 
64 private:
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  {
204  BaseLib::makeVectorUnique(elements);
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* 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* 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  }
289  BaseLib::makeVectorUnique(vec_nodes,
290  [](MeshLib::Node* node1, MeshLib::Node* node2)
291  { return node1->getID() < node2->getID(); });
292  DBUG("-> found {:d} nodes on the fracture {:d}", vec_nodes.size(),
293  frac_id);
294  }
295 
296  // find branch/junction nodes which connect to multiple fractures
297  std::vector<std::vector<MeshLib::Element*>> intersected_fracture_elements;
298  findFracutreIntersections(mesh, vec_fracture_mat_IDs, vec_fracture_nodes,
299  intersected_fracture_elements,
300  vec_branch_nodeID_matIDs,
301  vec_junction_nodeID_matIDs);
302 
303  // create a vector fracture elements and connected matrix elements,
304  // which are passed to a DoF table
305  for (unsigned fid = 0; fid < vec_fracture_elements.size(); fid++)
306  {
307  auto const& fracture_elements = vec_fracture_elements[fid];
308  std::vector<MeshLib::Element*> vec_ele;
309  // first, collect matrix elements
310  for (MeshLib::Element* e : fracture_elements)
311  {
312  // it is sufficient to iterate over base nodes, because they are
313  // already connected to all neighbours
314  for (unsigned i = 0; i < e->getNumberOfBaseNodes(); i++)
315  {
316  MeshLib::Node const* node = e->getNode(i);
317  if (isCrackTip(*node))
318  {
319  continue;
320  }
321  auto const& elements_connected_to_node =
322  mesh.getElementsConnectedToNode(*node);
323  for (unsigned j = 0; j < elements_connected_to_node.size(); j++)
324  {
325  // only matrix elements
326  if (elements_connected_to_node[j]->getDimension() <
327  mesh.getDimension())
328  {
329  continue;
330  }
331  vec_ele.push_back(const_cast<MeshLib::Element*>(
332  elements_connected_to_node[j]));
333  }
334  }
335  }
338  { return e1->getID() < e2->getID(); });
339 
340  // second, append fracture elements
341  std::copy(fracture_elements.begin(), fracture_elements.end(),
342  std::back_inserter(vec_ele));
343  // thirdly, append intersected fracture elements
344  std::copy(intersected_fracture_elements[fid].begin(),
345  intersected_fracture_elements[fid].end(),
346  std::back_inserter(vec_ele));
347 
348  vec_fracture_matrix_elements.push_back(vec_ele);
349  }
350 }
351 
352 } // namespace LIE
353 } // namespace ProcessLib
Definition of the Element class.
#define OGS_FATAL(...)
Definition: Error.h:26
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:27
Definition of the Mesh class.
Definition of the Node class.
std::size_t getID() const
Definition: Point3dWithID.h:62
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:82
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:71
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:98
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition: Mesh.cpp:232
const Node * getNode(std::size_t idx) const
Get the node with the given index.
Definition: Mesh.h:74
std::size_t getNumberOfElements() const
Get the number of elements.
Definition: Mesh.h:86
Node search class.
Definition: NodeSearch.h:25
std::size_t searchBoundaryNodes()
Marks all boundary nodes.
Definition: NodeSearch.cpp:76
const std::vector< std::size_t > & getSearchedNodeIDs() const
return marked node IDs
Definition: NodeSearch.h:30
void makeVectorUnique(std::vector< T > &v)
Definition: Algorithm.h:209
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:37
PropertyVector< int > const * materialIDs(Mesh const &mesh)
Definition: Mesh.cpp:258
bool isBaseNode(Node const &node, std::vector< Element const * > const &elements_connected_to_node)
Definition: Mesh.cpp:370
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)
Definition: MeshUtils.cpp:213
unsigned getDimension(MeshLib::MeshElemType eleType)