OGS 6.1.0-1721-g6382411ad
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  _is_internal_node[i] = false;
39  }
40 
41  bool operator()(MeshLib::Node const& node) const
42  {
43  if (!_is_internal_node[node.getID()] || !_mesh.isBaseNode(node.getID()))
44  return false;
45 
46  unsigned n_connected_fracture_elements = 0;
47  for (MeshLib::Element const* e : node.getElements())
48  if (e->getDimension() == _fracture_element_dim)
49  n_connected_fracture_elements++;
50  assert(n_connected_fracture_elements > 0);
51 
52  return (n_connected_fracture_elements == 1);
53  }
54 
55 private:
57  unsigned const _fracture_element_dim;
58  std::vector<bool> _is_internal_node;
59 };
60 
62  MeshLib::Mesh const& mesh,
63  std::vector<int> const& vec_fracture_mat_IDs,
64  std::vector<std::vector<MeshLib::Node*>> const& vec_fracture_nodes,
65  std::vector<std::vector<MeshLib::Element*>>& intersected_fracture_elements,
66  std::vector<std::pair<std::size_t, std::vector<int>>>&
67  vec_branch_nodeID_matIDs,
68  std::vector<std::pair<std::size_t, std::vector<int>>>&
69  vec_junction_nodeID_matIDs)
70 {
71  auto const n_fractures = vec_fracture_mat_IDs.size();
72  std::map<unsigned, unsigned> matID_to_fid;
73  for (unsigned i = 0; i < n_fractures; i++)
74  matID_to_fid[vec_fracture_mat_IDs[i]] = i;
75 
76  // make a vector all fracture nodes
77  std::vector<std::size_t> all_fracture_nodes;
78  for (auto& vec : vec_fracture_nodes)
79  for (auto* node : vec)
80  all_fracture_nodes.push_back(node->getID());
81 
82  // create a table of a node id and connected material IDs
83  std::map<std::size_t, std::vector<std::size_t>> frac_nodeID_to_matIDs;
84  for (unsigned i = 0; i < n_fractures; i++)
85  for (auto* node : vec_fracture_nodes[i])
86  frac_nodeID_to_matIDs[node->getID()].push_back(
87  vec_fracture_mat_IDs[i]);
88 
89  auto const* const opt_material_ids = MeshLib::materialIDs(mesh);
90 
91  // find branch/junction nodes which connect to multiple fractures
92  intersected_fracture_elements.resize(n_fractures);
93  for (auto entry : frac_nodeID_to_matIDs)
94  {
95  auto nodeID = entry.first;
96  auto const* node = mesh.getNode(entry.first);
97  auto const& matIDs = entry.second;
98  if (matIDs.size() < 2)
99  continue; // no intersection
100 
101  std::vector<MeshLib::Element*> conn_fracture_elements;
102  {
103  for (auto const* e : node->getElements())
104  if (e->getDimension() == (mesh.getDimension() - 1))
105  conn_fracture_elements.push_back(
106  const_cast<MeshLib::Element*>(e));
107  }
108 
109  std::map<int, int> vec_matID_counts;
110  {
111  for (auto matid : matIDs)
112  vec_matID_counts[matid] = 0;
113 
114  for (auto const* e : conn_fracture_elements)
115  {
116  auto matid = (*opt_material_ids)[e->getID()];
117  vec_matID_counts[matid]++;
118  }
119  }
120 
121  for (auto matid : matIDs)
122  {
123  auto fid = matID_to_fid[matid];
124  for (auto* e : conn_fracture_elements)
125  {
126  auto e_matid = (*opt_material_ids)[e->getID()];
127  if (matID_to_fid[e_matid] != fid) // only slave elements
128  intersected_fracture_elements[fid].push_back(e);
129  }
130  }
131 
132  bool isBranch = false;
133  {
134  for (auto entry : vec_matID_counts)
135  {
136  auto count = entry.second;
137  if (count % 2 == 1)
138  {
139  isBranch = true;
140  break;
141  }
142  }
143  }
144 
145  if (isBranch)
146  {
147  std::vector<int> matIDs(2);
148  for (auto entry : vec_matID_counts)
149  {
150  auto matid = entry.first;
151  auto count = entry.second;
152  if (count % 2 == 0)
153  {
154  matIDs[0] = matid; // master
155  }
156  else
157  {
158  matIDs[1] = matid; // slave
159  }
160  }
161  vec_branch_nodeID_matIDs.push_back(std::make_pair(nodeID, matIDs));
162  }
163  else
164  {
165  std::vector<int> matIDs(2);
166  matIDs[0] = std::min(vec_matID_counts.begin()->first,
167  vec_matID_counts.rbegin()->first);
168  matIDs[1] = std::max(vec_matID_counts.begin()->first,
169  vec_matID_counts.rbegin()->first);
170  vec_junction_nodeID_matIDs.push_back(
171  std::make_pair(nodeID, matIDs));
172  }
173  }
174 
175  for (auto& eles : intersected_fracture_elements)
177 
178  DBUG("-> found %d branchs and %d junctions",
179  vec_branch_nodeID_matIDs.size(), vec_junction_nodeID_matIDs.size());
180 }
181 
182 } // namespace
183 
185  MeshLib::Mesh const& mesh,
186  std::vector<MeshLib::Element*>& vec_matrix_elements,
187  std::vector<int>& vec_fracture_mat_IDs,
188  std::vector<std::vector<MeshLib::Element*>>& vec_fracture_elements,
189  std::vector<std::vector<MeshLib::Element*>>& vec_fracture_matrix_elements,
190  std::vector<std::vector<MeshLib::Node*>>& vec_fracture_nodes,
191  std::vector<std::pair<std::size_t, std::vector<int>>>&
192  vec_branch_nodeID_matIDs,
193  std::vector<std::pair<std::size_t, std::vector<int>>>&
194  vec_junction_nodeID_matIDs)
195 {
196  IsCrackTip isCrackTip(mesh);
197 
198  // get vectors of matrix elements and fracture elements
199  vec_matrix_elements.reserve(mesh.getNumberOfElements());
200  std::vector<MeshLib::Element*> all_fracture_elements;
201  for (MeshLib::Element* e : mesh.getElements())
202  {
203  if (e->getDimension() == mesh.getDimension())
204  vec_matrix_elements.push_back(e);
205  else
206  all_fracture_elements.push_back(e);
207  }
208  DBUG("-> found total %d matrix elements and %d fracture elements",
209  vec_matrix_elements.size(), all_fracture_elements.size());
210 
211  // get fracture material IDs
212  auto const material_ids = materialIDs(mesh);
213  if (!material_ids)
214  {
215  OGS_FATAL("Could not access MaterialIDs property from mesh.");
216  }
217  for (MeshLib::Element* e : all_fracture_elements)
218  vec_fracture_mat_IDs.push_back((*material_ids)[e->getID()]);
219  BaseLib::makeVectorUnique(vec_fracture_mat_IDs);
220  DBUG("-> found %d fracture material groups", vec_fracture_mat_IDs.size());
221 
222  // create a vector of fracture elements for each material
223  vec_fracture_elements.resize(vec_fracture_mat_IDs.size());
224  for (unsigned frac_id = 0; frac_id < vec_fracture_mat_IDs.size(); frac_id++)
225  {
226  const auto frac_mat_id = vec_fracture_mat_IDs[frac_id];
227  std::vector<MeshLib::Element*>& vec_elements =
228  vec_fracture_elements[frac_id];
229  std::copy_if(all_fracture_elements.begin(), all_fracture_elements.end(),
230  std::back_inserter(vec_elements),
231  [&](MeshLib::Element* e) {
232  return (*material_ids)[e->getID()] == frac_mat_id;
233  });
234  DBUG("-> found %d elements on the fracture %d", vec_elements.size(),
235  frac_id);
236  }
237 
238  // get a vector of fracture nodes for each material
239  vec_fracture_nodes.resize(vec_fracture_mat_IDs.size());
240  for (unsigned frac_id = 0; frac_id < vec_fracture_mat_IDs.size(); frac_id++)
241  {
242  std::vector<MeshLib::Node*>& vec_nodes = vec_fracture_nodes[frac_id];
243  for (MeshLib::Element* e : vec_fracture_elements[frac_id])
244  {
245  for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
246  {
247  if (isCrackTip(*e->getNode(i)))
248  continue;
249  vec_nodes.push_back(const_cast<MeshLib::Node*>(e->getNode(i)));
250  }
251  }
253  vec_nodes, [](MeshLib::Node* node1, MeshLib::Node* node2) {
254  return node1->getID() < node2->getID();
255  });
256  DBUG("-> found %d nodes on the fracture %d", vec_nodes.size(), frac_id);
257  }
258 
259  // find branch/junction nodes which connect to multiple fractures
260  std::vector<std::vector<MeshLib::Element*>> intersected_fracture_elements;
261  findFracutreIntersections(mesh, vec_fracture_mat_IDs, vec_fracture_nodes,
262  intersected_fracture_elements,
263  vec_branch_nodeID_matIDs,
264  vec_junction_nodeID_matIDs);
265 
266  // create a vector fracture elements and connected matrix elements,
267  // which are passed to a DoF table
268  for (unsigned fid = 0; fid < vec_fracture_elements.size(); fid++)
269  {
270  auto const& fracture_elements = vec_fracture_elements[fid];
271  std::vector<MeshLib::Element*> vec_ele;
272  // first, collect matrix elements
273  for (MeshLib::Element* e : fracture_elements)
274  {
275  // it is sufficient to iterate over base nodes, because they are
276  // already connected to all neighbours
277  for (unsigned i = 0; i < e->getNumberOfBaseNodes(); i++)
278  {
279  MeshLib::Node const* node = e->getNode(i);
280  if (isCrackTip(*node))
281  continue;
282  for (unsigned j = 0; j < node->getNumberOfElements(); j++)
283  {
284  // only matrix elements
285  if (node->getElement(j)->getDimension() <
286  mesh.getDimension())
287  continue;
288  vec_ele.push_back(
289  const_cast<MeshLib::Element*>(node->getElement(j)));
290  }
291  }
292  }
294  vec_ele, [](MeshLib::Element* e1, MeshLib::Element* e2) {
295  return e1->getID() < e2->getID();
296  });
297 
298  // second, append fracture elements
299  std::copy(fracture_elements.begin(), fracture_elements.end(),
300  std::back_inserter(vec_ele));
301  // thirdly, append intersected fracture elements
302  std::copy(intersected_fracture_elements[fid].begin(),
303  intersected_fracture_elements[fid].end(),
304  std::back_inserter(vec_ele));
305 
306  vec_fracture_matrix_elements.push_back(vec_ele);
307  }
308 }
309 
310 } // namespace LIE
311 } // 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:65
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:184
const Element * getElement(std::size_t idx) const
Get an element the node is part of.
Definition: Node.h:62
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:63
std::size_t getNumberOfElements() const
Get number of elements the node is part of.
Definition: Node.h:68
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:352
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:73
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:71
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:61