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#include <range/v3/algorithm/contains.hpp>
15
16#include "BranchProperty.h"
17#include "FractureProperty.h"
18#include "JunctionProperty.h"
19
20namespace
21{
22// Heaviside step function
23inline double Heaviside(bool const v)
24{
25 return v ? 0.5 : -0.5;
26}
27
28} // namespace
29
30namespace ProcessLib
31{
32namespace LIE
33{
34bool levelsetFracture(FractureProperty const& frac, Eigen::Vector3d const& x)
35{
36 return frac.normal_vector.dot(x - frac.point_on_fracture) > 0;
37}
38
39bool levelsetBranch(BranchProperty const& branch, Eigen::Vector3d const& x)
40{
41 return branch.normal_vector_branch.dot(x - branch.coords) > 0;
42}
43
44std::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
82std::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 (!ranges::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 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