OGS 6.1.0-1699-ge946d4c5f
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  continue;
109 
110  auto const eid = org_e->getID();
111  // keep original if the element has levelset=0
112  if ((*prop_levelset)[eid] == 0)
113  continue;
114 
115  // replace fracture nodes with duplicated ones
116  MeshLib::Element* e = new_eles[eid];
117  for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
118  {
119  const auto node_id = e->getNodeIndex(i);
120  if (!includesNodeID(vec_fracture_nodes, node_id))
121  continue;
122 
123  // list of duplicated node IDs
124  auto itr = _map_dup_newNodeIDs.find(node_id);
125  if (itr == _map_dup_newNodeIDs.end())
126  continue;
127  const auto& dup_newNodeIDs = itr->second;
128 
129  // choose new node id
130  unsigned new_node_id = 0;
131  if (dup_newNodeIDs.size() == 1)
132  {
133  // non-intersected nodes
134  new_node_id = dup_newNodeIDs[0];
135  }
136  else if (dup_newNodeIDs.size() == 2)
137  {
138  // branch nodes
139  const auto& br_matids = getMaterialIdsForNode(
140  vec_branch_nodeID_matIDs, node_id);
141  auto const frac2_matid = BaseLib::findFirstNotEqualElement(
142  br_matids, frac_matid);
143  assert(frac2_matid);
144  auto const frac2_id =
145  BaseLib::findIndex(vec_fracture_mat_IDs, *frac2_matid);
146  assert(frac2_id != std::numeric_limits<std::size_t>::max());
147  auto prop_levelset2 =
148  org_mesh.getProperties().getPropertyVector<double>(
149  "levelset" + std::to_string(frac2_id + 1));
150  std::size_t pos = 0;
151  if ((*prop_levelset2)[eid] == 0)
152  {
153  // index of this fracture
154  pos = BaseLib::findIndex(br_matids, frac_matid);
155  }
156  else if ((*prop_levelset2)[eid] == 1)
157  {
158  // index of the other fracture
159  pos = BaseLib::findIndex(br_matids, *frac2_matid);
160  }
161  assert(pos != std::numeric_limits<std::size_t>::max());
162  new_node_id = dup_newNodeIDs[pos];
163  }
164  else
165  {
166  // junction nodes
167  const auto& jct_matids = getMaterialIdsForNode(
168  vec_junction_nodeID_matIDs, node_id);
169  auto const frac2_matid = BaseLib::findFirstNotEqualElement(
170  jct_matids, frac_matid);
171  assert(frac2_matid);
172  auto const frac2_id =
173  BaseLib::findIndex(vec_fracture_mat_IDs, *frac2_matid);
174  assert(frac2_id != std::numeric_limits<std::size_t>::max());
175  auto prop_levelset2 =
176  org_mesh.getProperties().getPropertyVector<double>(
177  "levelset" + std::to_string(frac2_id + 1));
178 
179  //
180  if ((*prop_levelset2)[eid] == 0)
181  {
182  // index of this frac
183  auto const pos =
184  BaseLib::findIndex(jct_matids, frac_matid);
185  assert(pos != std::numeric_limits<std::size_t>::max());
186  new_node_id = dup_newNodeIDs[pos];
187  }
188  else if ((*prop_levelset2)[eid] == 1)
189  {
190  // set the last duplicated node
191  new_node_id = dup_newNodeIDs.back();
192  }
193  }
194 
195  // replace node
196  e->setNode(i, new_nodes[new_node_id]);
197  }
198  }
199  }
200 
201  // new mesh
202  _output_mesh = std::make_unique<MeshLib::Mesh>(org_mesh.getName(),
203  new_nodes, new_eles);
204  createProperties<int>();
205  createProperties<double>();
206  copyProperties<int>();
207  copyProperties<double>();
208  calculateTotalDisplacement(vec_vec_fracture_nodes.size(),
209  vec_junction_nodeID_matIDs.size());
210 }
211 
212 template <typename T>
214 {
215  MeshLib::Properties const& src_properties = _org_mesh.getProperties();
216  for (auto name : src_properties.getPropertyVectorNames())
217  {
218  if (!src_properties.existsPropertyVector<T>(name))
219  continue;
220  auto const* src_prop = src_properties.getPropertyVector<T>(name);
221 
222  auto const n_src_comp = src_prop->getNumberOfComponents();
223  // convert 2D vector to 3D. Otherwise Paraview Calculator filter does
224  // not recognize it as a vector
225  auto const n_dest_comp = (n_src_comp == 2) ? 3 : n_src_comp;
226 
227  auto new_prop = MeshLib::getOrCreateMeshProperty<T>(
228  *_output_mesh, name, src_prop->getMeshItemType(), n_dest_comp);
229 
230  if (src_prop->getMeshItemType() == MeshLib::MeshItemType::Node)
231  {
232  assert(new_prop->size() ==
233  _output_mesh->getNumberOfNodes() * n_dest_comp);
234  (void)(new_prop); // to avoid compilation warning.
235  }
236  else if (src_prop->getMeshItemType() == MeshLib::MeshItemType::Cell)
237  {
238  assert(new_prop->size() ==
239  _output_mesh->getNumberOfElements() * n_dest_comp);
240  }
241  else
242  {
243  WARN(
244  "Property '%s' cannot be created because its mesh item type is "
245  "not supported.",
246  name.c_str());
247  }
248  }
249 }
250 
251 template <typename T>
253 {
254  MeshLib::Properties const& src_properties = _org_mesh.getProperties();
255  for (auto name : src_properties.getPropertyVectorNames())
256  {
257  if (!src_properties.existsPropertyVector<T>(name))
258  continue;
259  auto const* src_prop = src_properties.getPropertyVector<T>(name);
260  auto* dest_prop =
261  _output_mesh->getProperties().getPropertyVector<T>(name);
262 
263  if (src_prop->getMeshItemType() == MeshLib::MeshItemType::Node)
264  {
265  auto const n_src_comp = src_prop->getNumberOfComponents();
266  auto const n_dest_comp = dest_prop->getNumberOfComponents();
267  // copy existing
268  for (unsigned i = 0; i < _org_mesh.getNumberOfNodes(); i++)
269  {
270  for (int j = 0; j < n_src_comp; j++)
271  (*dest_prop)[i * n_dest_comp + j] =
272  (*src_prop)[i * n_src_comp + j];
273  // set zero for components not existing in the original
274  for (int j = n_src_comp; j < n_dest_comp; j++)
275  (*dest_prop)[i * n_dest_comp + j] = 0;
276  }
277  // copy duplicated
278  for (auto itr : _map_dup_newNodeIDs)
279  {
280  for (int j = 0; j < n_dest_comp; j++)
281  for (unsigned k = 0; k < itr.second.size(); k++)
282  (*dest_prop)[itr.second[k] * n_dest_comp + j] =
283  (*dest_prop)[itr.first * n_dest_comp + j];
284  }
285  }
286  else if (src_prop->getMeshItemType() == MeshLib::MeshItemType::Cell)
287  {
288  std::copy(src_prop->begin(), src_prop->end(), dest_prop->begin());
289  }
290  else
291  {
292  WARN(
293  "Property '%s' cannot be created because its mesh item type is "
294  "not supported.",
295  name.c_str());
296  }
297  }
298 }
299 
300 void PostProcessTool::calculateTotalDisplacement(unsigned const n_fractures,
301  unsigned const n_junctions)
302 {
303  auto const& u = *_output_mesh->getProperties().getPropertyVector<double>(
304  "displacement");
305  auto const n_u_comp = u.getNumberOfComponents();
306  assert(u.size() == _output_mesh->getNodes().size() * 3);
307  auto& total_u =
308  *_output_mesh->getProperties().createNewPropertyVector<double>(
309  "u", MeshLib::MeshItemType::Node, n_u_comp);
310  total_u.resize(u.size());
311  for (unsigned i = 0; i < _output_mesh->getNodes().size(); i++)
312  {
313  for (int j = 0; j < n_u_comp; j++)
314  total_u[i * n_u_comp + j] = u[i * n_u_comp + j];
315  }
316 
317  for (unsigned enrich_id = 0; enrich_id < n_fractures + n_junctions;
318  enrich_id++)
319  {
320  // nodal value of levelset
321  std::vector<double> nodal_levelset(_output_mesh->getNodes().size(),
322  0.0);
323  auto const& ele_levelset =
324  *_output_mesh->getProperties().getPropertyVector<double>(
325  "levelset" + std::to_string(enrich_id + 1));
326  for (MeshLib::Element const* e : _output_mesh->getElements())
327  {
328  if (e->getDimension() != _output_mesh->getDimension())
329  continue;
330  const double e_levelset = ele_levelset[e->getID()];
331  if (e_levelset == 0)
332  continue;
333 
334  for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
335  nodal_levelset[e->getNodeIndex(i)] = e_levelset;
336  }
337 
338  // update total displacements
339  auto const& g =
340  *_output_mesh->getProperties().getPropertyVector<double>(
341  "displacement_jump" + std::to_string(enrich_id + 1));
342  for (unsigned i = 0; i < _output_mesh->getNodes().size(); i++)
343  {
344  for (int j = 0; j < n_u_comp; j++)
345  total_u[i * n_u_comp + j] +=
346  nodal_levelset[i] * g[i * n_u_comp + j];
347  }
348  }
349 }
350 
351 } // namespace LIE
352 } // 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:131
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:156
boost::optional< typename Container::value_type > findFirstNotEqualElement(Container const &container, typename Container::value_type const &element)
Definition: Algorithm.h:230
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:247
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
unsigned getNodeIndex(unsigned i) const
Definition: Element.cpp:164
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:71
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:300
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:36
Definition of the Element class.