OGS 6.2.1-97-g73d1aeda3
PostUtils.cpp
Go to the documentation of this file.
1 
9 #include "PostUtils.h"
10 
13 #include "MeshLib/Node.h"
14 
15 namespace ProcessLib
16 {
17 namespace LIE
18 {
19 namespace
20 {
21 bool includesNodeID(std::vector<MeshLib::Node*> const& vec_nodes,
22  std::size_t node_id)
23 {
24  auto itr2 = std::find_if(
25  vec_nodes.begin(), vec_nodes.end(),
26  [&](MeshLib::Node const* node) { return node->getID() == node_id; });
27  return (itr2 != vec_nodes.end());
28 }
29 
30 std::vector<int> const& getMaterialIdsForNode(
31  std::vector<std::pair<std::size_t, std::vector<int>>> const&
32  vec_nodeID_matIDs,
33  std::size_t nodeID)
34 {
35  auto itr = std::find_if(
36  vec_nodeID_matIDs.begin(), vec_nodeID_matIDs.end(),
37  [&](std::pair<std::size_t, std::vector<int>> const& entry) {
38  return entry.first == nodeID;
39  });
40  assert(itr != vec_nodeID_matIDs.end());
41  return itr->second;
42 };
43 } // namespace
44 
46  MeshLib::Mesh const& org_mesh,
47  std::vector<int> const& vec_fracture_mat_IDs,
48  std::vector<std::vector<MeshLib::Node*>> const& vec_vec_fracture_nodes,
49  std::vector<std::vector<MeshLib::Element*>> const&
50  vec_vec_fracture_matrix_elements,
51  std::vector<std::pair<std::size_t, std::vector<int>>> const&
52  vec_branch_nodeID_matIDs,
53  std::vector<std::pair<std::size_t, std::vector<int>>> const&
54  vec_junction_nodeID_matIDs)
55  : _org_mesh(org_mesh)
56 {
57  if (!org_mesh.getProperties().hasPropertyVector("displacement") ||
58  !org_mesh.getProperties().hasPropertyVector("displacement_jump1") ||
59  !org_mesh.getProperties().hasPropertyVector("levelset1"))
60  {
61  OGS_FATAL("The given mesh does not have relevant properties");
62  }
63 
64  // clone nodes and elements
65  std::vector<MeshLib::Node*> new_nodes(
66  MeshLib::copyNodeVector(org_mesh.getNodes()));
67  std::vector<MeshLib::Element*> new_eles(
68  MeshLib::copyElementVector(org_mesh.getElements(), new_nodes));
69 
70  // duplicate fracture nodes (two dup. nodes created at a branch)
71  for (auto const& vec_fracture_nodes : vec_vec_fracture_nodes)
72  {
73  for (auto const* org_node : vec_fracture_nodes)
74  {
75  auto duplicated_node =
76  new MeshLib::Node(org_node->getCoords(), new_nodes.size());
77  new_nodes.push_back(duplicated_node);
78  _map_dup_newNodeIDs[org_node->getID()].push_back(
79  duplicated_node->getID());
80  }
81  }
82  // at a junction, generate one more duplicated node (total 4 nodes)
83  for (auto& entry : vec_junction_nodeID_matIDs)
84  {
85  auto* org_node = org_mesh.getNode(entry.first);
86  auto duplicated_node =
87  new MeshLib::Node(org_node->getCoords(), new_nodes.size());
88  new_nodes.push_back(duplicated_node);
89  _map_dup_newNodeIDs[org_node->getID()].push_back(
90  duplicated_node->getID());
91  }
92 
93  // split elements using the new duplicated nodes
94  for (unsigned fracture_id = 0;
95  fracture_id < vec_vec_fracture_matrix_elements.size();
96  fracture_id++)
97  {
98  auto const frac_matid = vec_fracture_mat_IDs[fracture_id];
99  auto const& vec_fracture_matrix_elements =
100  vec_vec_fracture_matrix_elements[fracture_id];
101  auto const& vec_fracture_nodes = vec_vec_fracture_nodes[fracture_id];
102  auto prop_levelset = org_mesh.getProperties().getPropertyVector<double>(
103  "levelset" + std::to_string(fracture_id + 1));
104  for (auto const* org_e : vec_fracture_matrix_elements)
105  {
106  // only matrix elements
107  if (org_e->getDimension() != org_mesh.getDimension())
108  {
109  continue;
110  }
111 
112  auto const eid = org_e->getID();
113  // keep original if the element has levelset=0
114  if ((*prop_levelset)[eid] == 0)
115  {
116  continue;
117  }
118 
119  // replace fracture nodes with duplicated ones
120  MeshLib::Element* e = new_eles[eid];
121  for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
122  {
123  const auto node_id = e->getNodeIndex(i);
124  if (!includesNodeID(vec_fracture_nodes, node_id))
125  {
126  continue;
127  }
128 
129  // list of duplicated node IDs
130  auto itr = _map_dup_newNodeIDs.find(node_id);
131  if (itr == _map_dup_newNodeIDs.end())
132  {
133  continue;
134  }
135  const auto& dup_newNodeIDs = itr->second;
136 
137  // choose new node id
138  unsigned new_node_id = 0;
139  if (dup_newNodeIDs.size() == 1)
140  {
141  // non-intersected nodes
142  new_node_id = dup_newNodeIDs[0];
143  }
144  else if (dup_newNodeIDs.size() == 2)
145  {
146  // branch nodes
147  const auto& br_matids = getMaterialIdsForNode(
148  vec_branch_nodeID_matIDs, node_id);
149  auto const frac2_matid = BaseLib::findFirstNotEqualElement(
150  br_matids, frac_matid);
151  assert(frac2_matid);
152  auto const frac2_id =
153  BaseLib::findIndex(vec_fracture_mat_IDs, *frac2_matid);
154  assert(frac2_id != std::numeric_limits<std::size_t>::max());
155  auto prop_levelset2 =
156  org_mesh.getProperties().getPropertyVector<double>(
157  "levelset" + std::to_string(frac2_id + 1));
158  std::size_t pos = 0;
159  if ((*prop_levelset2)[eid] == 0)
160  {
161  // index of this fracture
162  pos = BaseLib::findIndex(br_matids, frac_matid);
163  }
164  else if ((*prop_levelset2)[eid] == 1)
165  {
166  // index of the other fracture
167  pos = BaseLib::findIndex(br_matids, *frac2_matid);
168  }
169  assert(pos != std::numeric_limits<std::size_t>::max());
170  new_node_id = dup_newNodeIDs[pos];
171  }
172  else
173  {
174  // junction nodes
175  const auto& jct_matids = getMaterialIdsForNode(
176  vec_junction_nodeID_matIDs, node_id);
177  auto const frac2_matid = BaseLib::findFirstNotEqualElement(
178  jct_matids, frac_matid);
179  assert(frac2_matid);
180  auto const frac2_id =
181  BaseLib::findIndex(vec_fracture_mat_IDs, *frac2_matid);
182  assert(frac2_id != std::numeric_limits<std::size_t>::max());
183  auto prop_levelset2 =
184  org_mesh.getProperties().getPropertyVector<double>(
185  "levelset" + std::to_string(frac2_id + 1));
186 
187  //
188  if ((*prop_levelset2)[eid] == 0)
189  {
190  // index of this frac
191  auto const pos =
192  BaseLib::findIndex(jct_matids, frac_matid);
193  assert(pos != std::numeric_limits<std::size_t>::max());
194  new_node_id = dup_newNodeIDs[pos];
195  }
196  else if ((*prop_levelset2)[eid] == 1)
197  {
198  // set the last duplicated node
199  new_node_id = dup_newNodeIDs.back();
200  }
201  }
202 
203  // replace node
204  e->setNode(i, new_nodes[new_node_id]);
205  }
206  }
207  }
208 
209  // new mesh
210  _output_mesh = std::make_unique<MeshLib::Mesh>(org_mesh.getName(),
211  new_nodes, new_eles);
212  createProperties<int>();
213  createProperties<double>();
214  copyProperties<int>();
215  copyProperties<double>();
216  calculateTotalDisplacement(vec_vec_fracture_nodes.size(),
217  vec_junction_nodeID_matIDs.size());
218 }
219 
220 template <typename T>
222 {
223  MeshLib::Properties const& src_properties = _org_mesh.getProperties();
224  for (auto name : src_properties.getPropertyVectorNames())
225  {
226  if (!src_properties.existsPropertyVector<T>(name))
227  {
228  continue;
229  }
230  auto const* src_prop = src_properties.getPropertyVector<T>(name);
231 
232  auto const n_src_comp = src_prop->getNumberOfComponents();
233  // convert 2D vector to 3D. Otherwise Paraview Calculator filter does
234  // not recognize it as a vector
235  auto const n_dest_comp = (n_src_comp == 2) ? 3 : n_src_comp;
236 
237  auto new_prop = MeshLib::getOrCreateMeshProperty<T>(
238  *_output_mesh, name, src_prop->getMeshItemType(), n_dest_comp);
239 
240  if (src_prop->getMeshItemType() == MeshLib::MeshItemType::Node)
241  {
242  assert(new_prop->size() ==
243  _output_mesh->getNumberOfNodes() * n_dest_comp);
244  (void)(new_prop); // to avoid compilation warning.
245  }
246  else if (src_prop->getMeshItemType() == MeshLib::MeshItemType::Cell)
247  {
248  assert(new_prop->size() ==
249  _output_mesh->getNumberOfElements() * n_dest_comp);
250  }
251  else
252  {
253  WARN(
254  "Property '%s' cannot be created because its mesh item type is "
255  "not supported.",
256  name.c_str());
257  }
258  }
259 }
260 
261 template <typename T>
263 {
264  MeshLib::Properties const& src_properties = _org_mesh.getProperties();
265  for (auto name : src_properties.getPropertyVectorNames())
266  {
267  if (!src_properties.existsPropertyVector<T>(name))
268  {
269  continue;
270  }
271  auto const* src_prop = src_properties.getPropertyVector<T>(name);
272  auto* dest_prop =
273  _output_mesh->getProperties().getPropertyVector<T>(name);
274 
275  if (src_prop->getMeshItemType() == MeshLib::MeshItemType::Node)
276  {
277  auto const n_src_comp = src_prop->getNumberOfComponents();
278  auto const n_dest_comp = dest_prop->getNumberOfComponents();
279  // copy existing
280  for (unsigned i = 0; i < _org_mesh.getNumberOfNodes(); i++)
281  {
282  for (int j = 0; j < n_src_comp; j++)
283  {
284  (*dest_prop)[i * n_dest_comp + j] =
285  (*src_prop)[i * n_src_comp + j];
286  }
287  // set zero for components not existing in the original
288  for (int j = n_src_comp; j < n_dest_comp; j++)
289  {
290  (*dest_prop)[i * n_dest_comp + j] = 0;
291  }
292  }
293  // copy duplicated
294  for (auto itr : _map_dup_newNodeIDs)
295  {
296  for (int j = 0; j < n_dest_comp; j++)
297  {
298  for (unsigned k = 0; k < itr.second.size(); k++)
299  {
300  (*dest_prop)[itr.second[k] * n_dest_comp + j] =
301  (*dest_prop)[itr.first * n_dest_comp + j];
302  }
303  }
304  }
305  }
306  else if (src_prop->getMeshItemType() == MeshLib::MeshItemType::Cell)
307  {
308  std::copy(src_prop->begin(), src_prop->end(), dest_prop->begin());
309  }
310  else
311  {
312  WARN(
313  "Property '%s' cannot be created because its mesh item type is "
314  "not supported.",
315  name.c_str());
316  }
317  }
318 }
319 
320 void PostProcessTool::calculateTotalDisplacement(unsigned const n_fractures,
321  unsigned const n_junctions)
322 {
323  auto const& u = *_output_mesh->getProperties().getPropertyVector<double>(
324  "displacement");
325  auto const n_u_comp = u.getNumberOfComponents();
326  assert(u.size() == _output_mesh->getNodes().size() * 3);
327  auto& total_u =
328  *_output_mesh->getProperties().createNewPropertyVector<double>(
329  "u", MeshLib::MeshItemType::Node, n_u_comp);
330  total_u.resize(u.size());
331  for (unsigned i = 0; i < _output_mesh->getNodes().size(); i++)
332  {
333  for (int j = 0; j < n_u_comp; j++)
334  {
335  total_u[i * n_u_comp + j] = u[i * n_u_comp + j];
336  }
337  }
338 
339  for (unsigned enrich_id = 0; enrich_id < n_fractures + n_junctions;
340  enrich_id++)
341  {
342  // nodal value of levelset
343  std::vector<double> nodal_levelset(_output_mesh->getNodes().size(),
344  0.0);
345  auto const& ele_levelset =
346  *_output_mesh->getProperties().getPropertyVector<double>(
347  "levelset" + std::to_string(enrich_id + 1));
348  for (MeshLib::Element const* e : _output_mesh->getElements())
349  {
350  if (e->getDimension() != _output_mesh->getDimension())
351  {
352  continue;
353  }
354  const double e_levelset = ele_levelset[e->getID()];
355  if (e_levelset == 0)
356  {
357  continue;
358  }
359 
360  for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
361  {
362  nodal_levelset[e->getNodeIndex(i)] = e_levelset;
363  }
364  }
365 
366  // update total displacements
367  auto const& g =
368  *_output_mesh->getProperties().getPropertyVector<double>(
369  "displacement_jump" + std::to_string(enrich_id + 1));
370  for (unsigned i = 0; i < _output_mesh->getNodes().size(); i++)
371  {
372  for (int j = 0; j < n_u_comp; j++)
373  {
374  total_u[i * n_u_comp + j] +=
375  nodal_levelset[i] * g[i * n_u_comp + j];
376  }
377  }
378  }
379 }
380 
381 } // namespace LIE
382 } // namespace ProcessLib
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition: Mesh.h:99
std::vector< MeshLib::Element * > copyElementVector(const std::vector< MeshLib::Element *> &elements, const std::vector< MeshLib::Node *> &nodes)
const Node * getNode(std::size_t idx) const
Get the node with the given index.
Definition: Mesh.h:84
Definition of the Node class.
std::vector< std::string > getPropertyVectorNames() const
Definition: Properties.cpp:37
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition: Mesh.h:105
MeshLib::Properties & getProperties()
Definition: Mesh.h:134
MeshLib::Mesh const & _org_mesh
Definition: PostUtils.h:53
PostProcessTool(MeshLib::Mesh const &org_mesh, std::vector< int > const &vec_fracture_mat_IDs, std::vector< std::vector< MeshLib::Node *>> const &vec_vec_fracture_nodes, std::vector< std::vector< MeshLib::Element *>> const &vec_vec_fracture_matrix_elements, std::vector< std::pair< std::size_t, std::vector< int >>> const &vec_branch_nodeID_matIDs, std::vector< std::pair< std::size_t, std::vector< int >>> const &vec_junction_nodeID_matIDs)
Definition: PostUtils.cpp:45
bool includesNodeID(std::vector< MeshLib::Node *> const &vec_nodes, std::size_t node_id)
Definition: PostUtils.cpp:21
void setNode(unsigned idx, Node *node)
Definition: Element.cpp:172
boost::optional< typename Container::value_type > findFirstNotEqualElement(Container const &container, typename Container::value_type const &element)
Definition: Algorithm.h:257
std::map< std::size_t, std::vector< std::size_t > > _map_dup_newNodeIDs
Definition: PostUtils.h:55
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties...
Definition: Properties.h:37
const std::string getName() const
Get name of the mesh.
Definition: Mesh.h:102
std::unique_ptr< MeshLib::Mesh > _output_mesh
Definition: PostUtils.h:54
std::size_t findIndex(Container const &container, typename Container::value_type const &element)
Definition: Algorithm.h:274
std::size_t getNodeIndex(unsigned i) const
Definition: Element.cpp:182
Definition of Duplicate functions.
std::vector< MeshLib::Node * > copyNodeVector(const std::vector< MeshLib::Node *> &nodes)
Creates a deep copy of a Node vector.
PropertyVector< T > const * getPropertyVector(std::string const &name) const
Definition: Properties.h:119
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:108
bool existsPropertyVector(std::string const &name) const
Definition: Properties.h:79
std::vector< int > const & getMaterialIdsForNode(std::vector< std::pair< std::size_t, std::vector< int >>> const &vec_nodeID_matIDs, std::size_t nodeID)
Definition: PostUtils.cpp:30
bool hasPropertyVector(std::string const &name) const
Definition: Properties.cpp:32
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:81
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
virtual unsigned getNumberOfNodes() const =0
Returns the number of all nodes including both linear and nonlinear nodes.
void calculateTotalDisplacement(unsigned const n_fractures, unsigned const n_junctions)
Definition: PostUtils.cpp:320
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:36
Definition of the Element class.