OGS 6.2.0-97-g4a610c866
MeshUtils.cpp
Go to the documentation of this file.
1 
9 #include "MeshUtils.h"
10 
11 #include "BaseLib/Algorithm.h"
13 #include "MeshLib/Mesh.h"
15 #include "MeshLib/Node.h"
16 
17 namespace ProcessLib
18 {
19 namespace LIE
20 {
21 namespace
22 {
23 // A class to check whether a node is located on a crack tip with
24 // the following conditions:
25 // - the number of connected fracture elements is one
26 // - the node is not located on a domain boundary
28 {
29 public:
30  explicit IsCrackTip(MeshLib::Mesh const& mesh)
31  : _mesh(mesh), _fracture_element_dim(mesh.getDimension() - 1)
32  {
33  _is_internal_node.resize(mesh.getNumberOfNodes(), true);
34 
35  MeshLib::NodeSearch nodeSearch(mesh);
36  nodeSearch.searchBoundaryNodes();
37  for (auto i : nodeSearch.getSearchedNodeIDs())
38  {
39  _is_internal_node[i] = false;
40  }
41  }
42 
43  bool operator()(MeshLib::Node const& node) const
44  {
45  if (!_is_internal_node[node.getID()] || !_mesh.isBaseNode(node.getID()))
46  {
47  return false;
48  }
49 
50  unsigned n_connected_fracture_elements = 0;
51  for (MeshLib::Element const* e : node.getElements())
52  {
53  if (e->getDimension() == _fracture_element_dim)
54  {
55  n_connected_fracture_elements++;
56  }
57  }
58  assert(n_connected_fracture_elements > 0);
59 
60  return (n_connected_fracture_elements == 1);
61  }
62 
63 private:
65  unsigned const _fracture_element_dim;
66  std::vector<bool> _is_internal_node;
67 };
68 
70  MeshLib::Mesh const& mesh,
71  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  for (auto* node : vec)
91  {
92  all_fracture_nodes.push_back(node->getID());
93  }
94  }
95 
96  // create a table of a node id and connected material IDs
97  std::map<std::size_t, std::vector<std::size_t>> frac_nodeID_to_matIDs;
98  for (unsigned i = 0; i < n_fractures; i++)
99  {
100  for (auto* node : vec_fracture_nodes[i])
101  {
102  frac_nodeID_to_matIDs[node->getID()].push_back(
103  vec_fracture_mat_IDs[i]);
104  }
105  }
106 
107  auto const* const opt_material_ids = MeshLib::materialIDs(mesh);
108 
109  // find branch/junction nodes which connect to multiple fractures
110  intersected_fracture_elements.resize(n_fractures);
111  for (auto entry : frac_nodeID_to_matIDs)
112  {
113  auto nodeID = entry.first;
114  auto const* node = mesh.getNode(entry.first);
115  auto const& matIDs = entry.second;
116  if (matIDs.size() < 2)
117  {
118  continue; // no intersection
119  }
120 
121  std::vector<MeshLib::Element*> conn_fracture_elements;
122  {
123  for (auto const* e : node->getElements())
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 entry : vec_matID_counts)
163  {
164  auto count = entry.second;
165  if (count % 2 == 1)
166  {
167  isBranch = true;
168  break;
169  }
170  }
171  }
172 
173  if (isBranch)
174  {
175  std::vector<int> matIDs(2);
176  for (auto entry : vec_matID_counts)
177  {
178  auto matid = entry.first;
179  auto count = entry.second;
180  if (count % 2 == 0)
181  {
182  matIDs[0] = matid; // master
183  }
184  else
185  {
186  matIDs[1] = matid; // slave
187  }
188  }
189  vec_branch_nodeID_matIDs.emplace_back(nodeID, matIDs);
190  }
191  else
192  {
193  std::vector<int> matIDs(2);
194  matIDs[0] = std::min(vec_matID_counts.begin()->first,
195  vec_matID_counts.rbegin()->first);
196  matIDs[1] = std::max(vec_matID_counts.begin()->first,
197  vec_matID_counts.rbegin()->first);
198  vec_junction_nodeID_matIDs.emplace_back(nodeID, matIDs);
199  }
200  }
201 
202  for (auto& eles : intersected_fracture_elements)
203  {
205  }
206 
207  DBUG("-> found %d branchs 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  for (MeshLib::Element* e : all_fracture_elements)
251  {
252  vec_fracture_mat_IDs.push_back((*material_ids)[e->getID()]);
253  }
254  BaseLib::makeVectorUnique(vec_fracture_mat_IDs);
255  DBUG("-> found %d fracture material groups", vec_fracture_mat_IDs.size());
256 
257  // create a vector of fracture elements for each material
258  vec_fracture_elements.resize(vec_fracture_mat_IDs.size());
259  for (unsigned frac_id = 0; frac_id < vec_fracture_mat_IDs.size(); frac_id++)
260  {
261  const auto frac_mat_id = vec_fracture_mat_IDs[frac_id];
262  std::vector<MeshLib::Element*>& vec_elements =
263  vec_fracture_elements[frac_id];
264  std::copy_if(all_fracture_elements.begin(), all_fracture_elements.end(),
265  std::back_inserter(vec_elements),
266  [&](MeshLib::Element* e) {
267  return (*material_ids)[e->getID()] == frac_mat_id;
268  });
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  }
290  vec_nodes, [](MeshLib::Node* node1, MeshLib::Node* node2) {
291  return node1->getID() < node2->getID();
292  });
293  DBUG("-> found %d nodes on the fracture %d", vec_nodes.size(), 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  for (unsigned j = 0; j < node->getNumberOfElements(); j++)
322  {
323  // only matrix elements
324  if (node->getElement(j)->getDimension() <
325  mesh.getDimension())
326  {
327  continue;
328  }
329  vec_ele.push_back(
330  const_cast<MeshLib::Element*>(node->getElement(j)));
331  }
332  }
333  }
335  vec_ele, [](MeshLib::Element* e1, MeshLib::Element* e2) {
336  return e1->getID() < e2->getID();
337  });
338 
339  // second, append fracture elements
340  std::copy(fracture_elements.begin(), fracture_elements.end(),
341  std::back_inserter(vec_ele));
342  // thirdly, append intersected fracture elements
343  std::copy(intersected_fracture_elements[fid].begin(),
344  intersected_fracture_elements[fid].end(),
345  std::back_inserter(vec_ele));
346 
347  vec_fracture_matrix_elements.push_back(vec_ele);
348  }
349 }
350 
351 } // namespace LIE
352 } // namespace ProcessLib
Node search class.
Definition: NodeSearch.h:23
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition: Mesh.h:99
const Node * getNode(std::size_t idx) const
Get the node with the given index.
Definition: Mesh.h:84
const std::vector< Element * > & getElements() const
Get all elements the node is part of.
Definition: Node.h:67
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
const Element * getElement(std::size_t idx) const
Get an element the node is part of.
Definition: Node.h:64
Definition of the Node class.
void makeVectorUnique(std::vector< T > &v)
Definition: Algorithm.h:177
virtual unsigned getDimension() const =0
Get dimension of the mesh element.
Definition of the Mesh class.
std::size_t getID() const
Definition: Point3dWithID.h:64
std::size_t getNumberOfElements() const
Get number of elements the node is part of.
Definition: Node.h:70
const std::vector< std::size_t > & getSearchedNodeIDs() const
return marked node IDs
Definition: NodeSearch.h:29
PropertyVector< int > const * materialIDs(Mesh const &mesh)
Definition: Mesh.cpp:403
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:108
std::size_t searchBoundaryNodes()
Marks all boundary nodes.
Definition: NodeSearch.cpp:79
unsigned getDimension(MeshLib::MeshElemType eleType)
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:81
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:90
std::size_t getNumberOfElements() const
Get the number of elements.
Definition: Mesh.h:96
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:36
Definition of the Element class.
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:69