OGS
LevelSetFunction.cpp
Go to the documentation of this file.
1 
10 #include "LevelSetFunction.h"
11 
12 #include <boost/math/special_functions/sign.hpp>
13 #include <numeric>
14 
15 #include "BaseLib/Algorithm.h"
16 #include "BranchProperty.h"
17 #include "FractureProperty.h"
18 #include "JunctionProperty.h"
19 
20 namespace
21 {
22 // Heaviside step function
23 inline double Heaviside(bool const v)
24 {
25  return v ? 0.5 : -0.5;
26 }
27 
28 } // namespace
29 
30 namespace ProcessLib
31 {
32 namespace LIE
33 {
34 bool levelsetFracture(FractureProperty const& frac, Eigen::Vector3d const& x)
35 {
36  return frac.normal_vector.dot(x - frac.point_on_fracture) > 0;
37 }
38 
39 bool levelsetBranch(BranchProperty const& branch, Eigen::Vector3d const& x)
40 {
41  return branch.normal_vector_branch.dot(x - branch.coords) > 0;
42 }
43 
44 std::vector<double> uGlobalEnrichments(
45  std::vector<FractureProperty*> const& frac_props,
46  std::vector<JunctionProperty*> const& junction_props,
47  std::unordered_map<int, int> const& fracID_to_local,
48  Eigen::Vector3d const& x)
49 {
50  // pre-calculate levelsets for all fractures
51  std::vector<bool> levelsets(frac_props.size());
52  for (std::size_t i = 0; i < frac_props.size(); i++)
53  {
54  levelsets[i] = levelsetFracture(*frac_props[i], x);
55  }
56 
57  std::vector<double> enrichments(frac_props.size() + junction_props.size());
58  // fractures possibly with branches
59  for (std::size_t i = 0; i < frac_props.size(); i++)
60  {
61  auto const* frac = frac_props[i];
62  enrichments[i] = Heaviside(
63  std::accumulate(cbegin(frac->branches_slave),
64  cend(frac->branches_slave), levelsets[i],
65  [&](bool const enrich, auto const& branch)
66  { return enrich & levelsetBranch(branch, x); }));
67  }
68 
69  // junctions
70  for (std::size_t i = 0; i < junction_props.size(); i++)
71  {
72  auto const* junction = junction_props[i];
73  auto fid1 = fracID_to_local.at(junction->fracture_ids[0]);
74  auto fid2 = fracID_to_local.at(junction->fracture_ids[1]);
75  bool const enrich = levelsets[fid1] && levelsets[fid2];
76  enrichments[i + frac_props.size()] = Heaviside(enrich);
77  }
78 
79  return enrichments;
80 }
81 
82 std::vector<double> duGlobalEnrichments(
83  std::size_t this_frac_id,
84  std::vector<FractureProperty*> const& frac_props,
85  std::vector<JunctionProperty*> const& junction_props,
86  std::unordered_map<int, int> const& fracID_to_local,
87  Eigen::Vector3d const& x)
88 {
89  auto this_frac_local_index = fracID_to_local.at(this_frac_id);
90  auto const& this_frac = *frac_props[this_frac_local_index];
91  // pre-calculate levelsets for all fractures
92  std::vector<bool> levelsets(frac_props.size());
93  for (std::size_t i = 0; i < frac_props.size(); i++)
94  {
95  levelsets[i] = levelsetFracture(*frac_props[i], x);
96  }
97 
98  std::vector<double> enrichments(frac_props.size() + junction_props.size());
99  enrichments[this_frac_local_index] = 1.0;
100 
101  // fractures possibly with branches
102  if (frac_props.size() > 1)
103  {
104  for (auto const& branch : this_frac.branches_master)
105  {
106  if (branch.master_fracture_id != this_frac.fracture_id)
107  {
108  continue;
109  }
110 
111  if (fracID_to_local.find(branch.slave_fracture_id) ==
112  fracID_to_local.end())
113  {
114  continue;
115  }
116 
117  double sign = boost::math::sign(
118  this_frac.normal_vector.dot(branch.normal_vector_branch));
119  auto slave_fid = fracID_to_local.at(branch.slave_fracture_id);
120  double const enrich = levelsets[slave_fid] ? 1. : 0.;
121  enrichments[slave_fid] = sign * enrich;
122  }
123  }
124 
125  // junctions
126  for (unsigned i = 0; i < junction_props.size(); i++)
127  {
128  auto const* junction = junction_props[i];
129  if (!BaseLib::contains(junction->fracture_ids, this_frac.fracture_id))
130  {
131  continue;
132  }
133 
134  auto another_frac_id =
135  (junction->fracture_ids[0] == this_frac.fracture_id)
136  ? junction->fracture_ids[1]
137  : junction->fracture_ids[0];
138  auto fid = fracID_to_local.at(another_frac_id);
139  double const enrich = levelsets[fid] ? 1. : 0.;
140  enrichments[i + frac_props.size()] = enrich;
141  }
142 
143  return enrichments;
144 }
145 
146 } // namespace LIE
147 } // namespace ProcessLib
bool contains(Container const &container, typename Container::value_type const &element)
Definition: Algorithm.h:256
static const double v
bool levelsetFracture(FractureProperty const &frac, Eigen::Vector3d const &x)
std::vector< double > duGlobalEnrichments(std::size_t this_frac_id, std::vector< FractureProperty * > const &frac_props, std::vector< JunctionProperty * > const &junction_props, std::unordered_map< int, int > const &fracID_to_local, Eigen::Vector3d const &x)
std::vector< double > uGlobalEnrichments(std::vector< FractureProperty * > const &frac_props, std::vector< JunctionProperty * > const &junction_props, std::unordered_map< int, int > const &fracID_to_local, Eigen::Vector3d const &x)
bool levelsetBranch(BranchProperty const &branch, Eigen::Vector3d const &x)
Eigen::Vector3d const coords
Eigen::Vector3d normal_vector_branch