LIE intersect (continued).

BaseLib/Algorithm.h

@@ -12,6 +12,7 @@
 #pragma once
 
 #include <algorithm>
+#include <boost/optional.hpp>
 #include <cassert>
 #include <typeindex>
 #include <typeinfo>
@@ -217,11 +218,40 @@ void uniquePushBack(Container& container,
         container.push_back(element);
 }
 
-template <typename Container, typename ValueType>
-inline bool contains(Container const& container, ValueType const& element)
+template <typename Container>
+bool contains(Container const& container,
+              typename Container::value_type const& element)
+{
+    return std::find(container.begin(), container.end(), element) !=
+           container.end();
+}
+
+template <typename Container>
+boost::optional<typename Container::value_type> findFirstNotEqualElement(
+    Container const& container, typename Container::value_type const& element)
 {
-    return (std::find(container.begin(), container.end(), element) !=
-            container.end());
+    auto const it =
+        std::find_if_not(container.begin(), container.end(),
+                         [&element](typename Container::value_type const& e) {
+                             return e == element;
+                         });
+    return it == container.end() ? boost::none : boost::make_optional(*it);
 }
 
+/// Returns the index of first element in container or, if the element is not
+/// found a std::size_t maximum value.
+///
+/// The maximum value of std::size_t is chosen, because such an index cannot
+/// exist in a container; the maximum index is std::size_t::max-1.
+template <typename Container>
+std::size_t findIndex(Container const& container,
+                      typename Container::value_type const& element)
+{
+    auto const it = std::find(container.begin(), container.end(), element);
+    if (it == container.end())
+    {
+        return std::numeric_limits<std::size_t>::max();
+    }
+    return std::distance(container.begin(), it);
+}
 }  // namespace BaseLib

ProcessLib/LIE/Common/BranchProperty.h

@@ -11,18 +11,30 @@
 
 #include <Eigen/Eigen>
 
+#include "MeshLib/Node.h"
+
 namespace ProcessLib
 {
 namespace LIE
 {
 struct BranchProperty final
 {
-    Eigen::Vector3d coords;
+    BranchProperty(MeshLib::Node const& branchNode,
+                   int const master_fracture_id_,
+                   int const slave_fracture_id_)
+        : coords{branchNode.getCoords()},
+          node_id{branchNode.getID()},
+          master_fracture_id{master_fracture_id_},
+          slave_fracture_id{slave_fracture_id_}
+    {
+    }
+
+    Eigen::Vector3d const coords;
     // unit vector normal to the master fracture in a direction to the slave
     Eigen::Vector3d normal_vector_branch;
-    int node_id;
-    int master_fracture_ID;
-    int slave_fracture_ID;
+    std::size_t const node_id;
+    int const master_fracture_id;
+    int const slave_fracture_id;
 
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW
 };

ProcessLib/LIE/Common/FractureProperty.h

@@ -39,64 +39,65 @@ struct FractureProperty
     /// Rotation matrix from global to local coordinates
     Eigen::MatrixXd R;
     /// Initial aperture
-    ProcessLib::Parameter<double> const* aperture0 = nullptr;
-    std::vector<std::unique_ptr<BranchProperty>> branches_master;
-    std::vector<std::unique_ptr<BranchProperty>> branches_slave;
+    ProcessLib::Parameter<double> const& aperture0;
+    std::vector<BranchProperty> branches_master;
+    std::vector<BranchProperty> branches_slave;
+
+    FractureProperty(int const fracture_id_, int const material_id,
+                     ProcessLib::Parameter<double> const& initial_aperture)
+        : fracture_id(fracture_id_),
+          mat_id(material_id),
+          aperture0(initial_aperture)
+    {
+    }
 
     virtual ~FractureProperty() = default;
 };
 
 struct FracturePropertyHM : public FractureProperty
 {
-    ProcessLib::Parameter<double> const* specific_storage = nullptr;
-    ProcessLib::Parameter<double> const* biot_coefficient = nullptr;
+    FracturePropertyHM(int const fracture_id_, int const material_id,
+                       ProcessLib::Parameter<double> const& initial_aperture,
+                       ProcessLib::Parameter<double> const& specific_storage_,
+                       ProcessLib::Parameter<double> const& biot_coefficient_)
+        : FractureProperty(fracture_id_, material_id, initial_aperture),
+          specific_storage(specific_storage_),
+          biot_coefficient(biot_coefficient_)
+    {
+    }
+    ProcessLib::Parameter<double> const& specific_storage;
+    ProcessLib::Parameter<double> const& biot_coefficient;
 };
 
 /// configure fracture property based on a fracture element assuming
 /// a fracture is a straight line/flat plane
-inline void setFractureProperty(unsigned dim, MeshLib::Element const& e,
+inline void setFractureProperty(int const dim, MeshLib::Element const& e,
                                 FractureProperty& frac_prop)
 {
     // 1st node is used but using other node is also possible, because
     // a fracture is not curving
     for (int j = 0; j < 3; j++)
+    {
         frac_prop.point_on_fracture[j] = e.getCenterOfGravity().getCoords()[j];
+    }
     computeNormalVector(e, dim, frac_prop.normal_vector);
     frac_prop.R.resize(dim, dim);
     computeRotationMatrix(e, frac_prop.normal_vector, dim, frac_prop.R);
 }
 
-inline BranchProperty* createBranchProperty(MeshLib::Node const& branchNode,
-                                            FractureProperty const& master_frac,
-                                            FractureProperty const& slave_frac)
+inline BranchProperty createBranchProperty(MeshLib::Node const& branchNode,
+                                           FractureProperty const& master_frac,
+                                           FractureProperty const& slave_frac)
 {
-    BranchProperty* branch = new BranchProperty();
-    branch->node_id = branchNode.getID();
-    branch->coords = Eigen::Vector3d(branchNode.getCoords());
-    branch->master_fracture_ID = master_frac.fracture_id;
-    branch->slave_fracture_ID = slave_frac.fracture_id;
+    BranchProperty branch{branchNode, master_frac.fracture_id,
+                          slave_frac.fracture_id};
+
     // set a normal vector from the master to the slave fracture
     Eigen::Vector3d branch_vector =
-        slave_frac.point_on_fracture - branch->coords;
+        slave_frac.point_on_fracture - branch.coords;
     double sign = (branch_vector.dot(master_frac.normal_vector) < 0) ? -1 : 1;
-    branch->normal_vector_branch = sign * master_frac.normal_vector;
+    branch.normal_vector_branch = sign * master_frac.normal_vector;
     return branch;
 }
-
-inline JunctionProperty* createJunctionProperty(
-    int junction_id,
-    MeshLib::Node const& junctionNode,
-    std::vector<int>
-        frac_ids)
-{
-    JunctionProperty* junction = new JunctionProperty();
-    junction->junction_id = junction_id;
-    junction->node_id = junctionNode.getID();
-    junction->coords = Eigen::Vector3d(junctionNode.getCoords());
-    for (int j = 0; j < 2; j++)
-        junction->fracture_IDs[j] = frac_ids[j];
-    return junction;
-}
-
 }  // namespace LIE
 }  // namespace ProcessLib

ProcessLib/LIE/Common/HMatrixUtils.h

@@ -58,7 +58,9 @@ void computeHMatrix(N_Type const& N, HMatrixType& H)
     H.setZero();
 
     for (unsigned j = 0; j < DisplacementDim; j++)
+    {
         H.block(j, j * NPOINTS, 1, NPOINTS) = N;
+    }
 }
 
 }  // namespace ProcessLib

ProcessLib/LIE/Common/JunctionProperty.h

@@ -9,20 +9,30 @@
 
 #pragma once
 
+#include <Eigen/Eigen>
 #include <array>
 
-#include <Eigen/Eigen>
+#include "MeshLib/Node.h"
 
 namespace ProcessLib
 {
 namespace LIE
 {
 struct JunctionProperty final
 {
-    Eigen::Vector3d coords;
-    int junction_id;
-    int node_id;
-    std::array<int, 2> fracture_IDs;
+    JunctionProperty(int const junction_id_,
+                     MeshLib::Node const& junctionNode,
+                     std::array<int, 2> const fracture_ids_)
+        : coords{junctionNode.getCoords()},
+          node_id{junctionNode.getID()},
+          fracture_ids{fracture_ids_},
+          junction_id{junction_id_}
+    {
+    }
+    Eigen::Vector3d const coords;
+    std::size_t const node_id;
+    std::array<int, 2> const fracture_ids;
+    int const junction_id;
 
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW
 };

ProcessLib/LIE/Common/LevelSetFunction.cpp

@@ -30,19 +30,19 @@ namespace ProcessLib
 {
 namespace LIE
 {
-double levelset_fracture(FractureProperty const& frac, Eigen::Vector3d const& x)
+double levelsetFracture(FractureProperty const& frac, Eigen::Vector3d const& x)
 {
     return boost::math::sign(
         frac.normal_vector.dot(x - frac.point_on_fracture));
 }
 
-double levelset_branch(BranchProperty const& branch, Eigen::Vector3d const& x)
+double levelsetBranch(BranchProperty const& branch, Eigen::Vector3d const& x)
 {
     return boost::math::sign(
         branch.normal_vector_branch.dot(x - branch.coords));
 }
 
-std::vector<double> u_global_enrichments(
+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,
@@ -51,7 +51,7 @@ std::vector<double> u_global_enrichments(
     // pre-calculate levelsets for all fractures
     std::vector<double> levelsets(frac_props.size());
     for (std::size_t i = 0; i < frac_props.size(); i++)
-        levelsets[i] = Heaviside(levelset_fracture(*frac_props[i], x));
+        levelsets[i] = Heaviside(levelsetFracture(*frac_props[i], x));
 
     std::vector<double> enrichments(frac_props.size() + junction_props.size());
     // fractures possibly with branches
@@ -60,24 +60,24 @@ std::vector<double> u_global_enrichments(
         auto const* frac = frac_props[i];
         double enrich = levelsets[i];
         for (std::size_t j = 0; j < frac->branches_slave.size(); j++)
-            enrich *= Heaviside(levelset_branch(*frac->branches_slave[j], x));
+            enrich *= Heaviside(levelsetBranch(frac->branches_slave[j], x));
         enrichments[i] = enrich;
     }
 
     // junctions
     for (std::size_t i = 0; i < junction_props.size(); i++)
     {
         auto const* junction = junction_props[i];
-        auto fid1 = fracID_to_local.at(junction->fracture_IDs[0]);
-        auto fid2 = fracID_to_local.at(junction->fracture_IDs[1]);
+        auto fid1 = fracID_to_local.at(junction->fracture_ids[0]);
+        auto fid2 = fracID_to_local.at(junction->fracture_ids[1]);
         double enrich = levelsets[fid1] * levelsets[fid2];
         enrichments[i + frac_props.size()] = enrich;
     }
 
     return enrichments;
 }
 
-std::vector<double> du_global_enrichments(
+std::vector<double> duGlobalEnrichments(
     std::size_t this_frac_id,
     std::vector<FractureProperty*> const& frac_props,
     std::vector<JunctionProperty*> const& junction_props,
@@ -89,7 +89,7 @@ std::vector<double> du_global_enrichments(
     // pre-calculate levelsets for all fractures
     std::vector<double> levelsets(frac_props.size());
     for (std::size_t i = 0; i < frac_props.size(); i++)
-        levelsets[i] = Heaviside(levelset_fracture(*frac_props[i], x));
+        levelsets[i] = Heaviside(levelsetFracture(*frac_props[i], x));
 
     std::vector<double> enrichments(frac_props.size() + junction_props.size());
     enrichments[this_frac_local_index] = 1.0;
@@ -99,16 +99,16 @@ std::vector<double> du_global_enrichments(
     {
         for (auto const& branch : this_frac.branches_master)
         {
-            if (branch->master_fracture_ID != this_frac.fracture_id)
+            if (branch.master_fracture_id != this_frac.fracture_id)
                 continue;
 
-            if (fracID_to_local.find(branch->slave_fracture_ID) ==
+            if (fracID_to_local.find(branch.slave_fracture_id) ==
                 fracID_to_local.end())
                 continue;
 
             double singned = boost::math::sign(
-                this_frac.normal_vector.dot(branch->normal_vector_branch));
-            auto slave_fid = fracID_to_local.at(branch->slave_fracture_ID);
+                this_frac.normal_vector.dot(branch.normal_vector_branch));
+            auto slave_fid = fracID_to_local.at(branch.slave_fracture_id);
             double enrich = singned * levelsets[slave_fid];
             enrichments[slave_fid] = enrich;
         }
@@ -118,13 +118,13 @@ std::vector<double> du_global_enrichments(
     for (unsigned i = 0; i < junction_props.size(); i++)
     {
         auto const* junction = junction_props[i];
-        if (!BaseLib::contains(junction->fracture_IDs, this_frac.fracture_id))
+        if (!BaseLib::contains(junction->fracture_ids, this_frac.fracture_id))
             continue;
 
         auto another_frac_id =
-            (junction->fracture_IDs[0] == this_frac.fracture_id)
-                ? junction->fracture_IDs[1]
-                : junction->fracture_IDs[0];
+            (junction->fracture_ids[0] == this_frac.fracture_id)
+                ? junction->fracture_ids[1]
+                : junction->fracture_ids[0];
         auto fid = fracID_to_local.at(another_frac_id);
         double enrich = levelsets[fid];
         enrichments[i + frac_props.size()] = enrich;

ProcessLib/LIE/Common/LevelSetFunction.h

@@ -30,7 +30,7 @@ struct JunctionProperty;
 /// in frac_props
 /// @param x evaluating point coordiates
 /// @return a vector of enrichment values for displacements
-std::vector<double> u_global_enrichments(
+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,
@@ -47,7 +47,7 @@ std::vector<double> u_global_enrichments(
 /// in frac_props
 /// @param x evaluating point coordiates
 /// @return a vector of enrichment values for fracture relative displacements
-std::vector<double> du_global_enrichments(
+std::vector<double> duGlobalEnrichments(
     std::size_t this_fracID,
     std::vector<FractureProperty*> const& frac_props,
     std::vector<JunctionProperty*> const& junction_props,