Merge pull request #2235 from norihiro-w/lie-intersects-pr

[L-IE] supporting fracture intersections in mechanics

Applications/Utils/PostProcessing/postLIE.cpp

@@ -42,12 +42,21 @@ void postVTU(std::string const& int_vtu_filename,
     std::vector<std::vector<MeshLib::Element*>> vec_fracture_elements;
     std::vector<std::vector<MeshLib::Element*>> vec_fracture_matrix_elements;
     std::vector<std::vector<MeshLib::Node*>> vec_fracture_nodes;
+    std::vector<std::pair<std::size_t, std::vector<int>>>
+        vec_branch_nodeID_matIDs;
+    std::vector<std::pair<std::size_t, std::vector<int>>>
+        vec_junction_nodeID_matIDs;
     ProcessLib::LIE::getFractureMatrixDataInMesh(
         *mesh, vec_matrix_elements, vec_fracture_mat_IDs, vec_fracture_elements,
-        vec_fracture_matrix_elements, vec_fracture_nodes);
-
-    ProcessLib::LIE::PostProcessTool post(*mesh, vec_fracture_nodes,
-                                          vec_fracture_matrix_elements);
+        vec_fracture_matrix_elements, vec_fracture_nodes,
+        vec_branch_nodeID_matIDs, vec_junction_nodeID_matIDs);
+
+    ProcessLib::LIE::PostProcessTool post(*mesh,
+                                          vec_fracture_mat_IDs,
+                                          vec_fracture_nodes,
+                                          vec_fracture_matrix_elements,
+                                          vec_branch_nodeID_matIDs,
+                                          vec_junction_nodeID_matIDs);
 
     // create a new VTU file
     INFO("create %s", out_vtu_filename.c_str());

BaseLib/Algorithm.h

@@ -216,4 +216,12 @@ void uniquePushBack(Container& container,
         container.end())
         container.push_back(element);
 }
+
+template <typename Container, typename ValueType>
+inline bool contains(Container const& container, ValueType const& element)
+{
+    return (std::find(container.begin(), container.end(), element) !=
+            container.end());
+}
+
 }  // namespace BaseLib

ProcessLib/LIE/Common/BranchProperty.h

@@ -0,0 +1,31 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <Eigen/Eigen>
+
+namespace ProcessLib
+{
+namespace LIE
+{
+struct BranchProperty final
+{
+    Eigen::Vector3d 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;
+
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+};
+
+}  // namespace LIE
+}  // namespace ProcessLib

ProcessLib/LIE/Common/FractureProperty.h

@@ -9,8 +9,12 @@
 
 #pragma once
 
+#include <memory>
+
 #include <Eigen/Eigen>
 
+#include "BranchProperty.h"
+#include "JunctionProperty.h"
 #include "Utils.h"
 
 namespace MeshLib
@@ -36,6 +40,8 @@ struct FractureProperty
     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;
 
     virtual ~FractureProperty() = default;
 };
@@ -54,11 +60,43 @@ inline void setFractureProperty(unsigned dim, MeshLib::Element const& e,
     // 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.getNode(0)->getCoords()[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)
+{
+    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;
+    // set a normal vector from the master to the slave fracture
+    Eigen::Vector3d branch_vector =
+        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;
+    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/JunctionProperty.h

@@ -0,0 +1,31 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <array>
+
+#include <Eigen/Eigen>
+
+namespace ProcessLib
+{
+namespace LIE
+{
+struct JunctionProperty final
+{
+    Eigen::Vector3d coords;
+    int junction_id;
+    int node_id;
+    std::array<int, 2> fracture_IDs;
+
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+};
+
+}  // namespace LIE
+}  // namespace ProcessLib

ProcessLib/LIE/Common/LevelSetFunction.cpp

@@ -10,7 +10,11 @@
 
 #include <boost/math/special_functions/sign.hpp>
 
+#include "BaseLib/Algorithm.h"
+
+#include "BranchProperty.h"
 #include "FractureProperty.h"
+#include "JunctionProperty.h"
 
 namespace
 {
@@ -26,11 +30,107 @@ namespace ProcessLib
 {
 namespace LIE
 {
-double calculateLevelSetFunction(FractureProperty const& frac, double const* x_)
+double levelset_fracture(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)
+{
+    return boost::math::sign(
+        branch.normal_vector_branch.dot(x - branch.coords));
+}
+
+std::vector<double> u_global_enrichments(
+    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)
+{
+    // 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));
+
+    std::vector<double> enrichments(frac_props.size() + junction_props.size());
+    // fractures possibly with branches
+    for (std::size_t i = 0; i < frac_props.size(); i++)
+    {
+        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));
+        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]);
+        double enrich = levelsets[fid1] * levelsets[fid2];
+        enrichments[i + frac_props.size()] = enrich;
+    }
+
+    return enrichments;
+}
+
+std::vector<double> du_global_enrichments(
+    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)
 {
-    Eigen::Map<Eigen::Vector3d const> x(x_, 3);
-    return Heaviside(
-        boost::math::sign(frac.normal_vector.dot(x - frac.point_on_fracture)));
+    auto this_frac_local_index = fracID_to_local.at(this_frac_id);
+    auto const& this_frac = *frac_props[this_frac_local_index];
+    // 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));
+
+    std::vector<double> enrichments(frac_props.size() + junction_props.size());
+    enrichments[this_frac_local_index] = 1.0;
+
+    // fractures possibly with branches
+    if (frac_props.size() > 1)
+    {
+        for (auto const& branch : this_frac.branches_master)
+        {
+            if (branch->master_fracture_ID != this_frac.fracture_id)
+                continue;
+
+            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);
+            double enrich = singned * levelsets[slave_fid];
+            enrichments[slave_fid] = enrich;
+        }
+    }
+
+    // junctions
+    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))
+            continue;
+
+        auto another_frac_id =
+            (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;
+    }
+
+    return enrichments;
 }
 
 }  // namespace LIE

ProcessLib/LIE/Common/LevelSetFunction.h

@@ -8,20 +8,51 @@
 
 #pragma once
 
+#include <unordered_map>
+#include <vector>
+
+#include <Eigen/Eigen>
+
 namespace ProcessLib
 {
 namespace LIE
 {
 struct FractureProperty;
+struct JunctionProperty;
+
+/// calculate the enrichment function for displacements at the given point
+/// Remarks:
+/// * branch/junction intersections of two fractures are supported in 2D
+///
+/// @param frac_props fracture properties
+/// @param junction_props junction properties
+/// @param fracID_to_local a mapping table from a fracture ID to a local index
+/// in frac_props
+/// @param x evaluating point coordiates
+/// @return a vector of enrichment values for displacements
+std::vector<double> u_global_enrichments(
+    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);
 
-/// calculate the level set function
-/// \f$ \psi(\mathbf{x}) = H(|\mathbf{x} - \mathbf{x_d}|
-/// \mathrm{sign}[\mathbf{n_d} \cdot (\mathbf{x}-\mathbf{x_d}]) \f$ where
-/// \f$H(u)\f$ is the Heaviside step function, \f$\mathbf{x_d}\f$ is a point on
-/// the fracture plane, and \f$\mathbf{n_d}\f$ is the normal vector of a
-/// fracture plane
-double calculateLevelSetFunction(FractureProperty const& fracture_property,
-                                 double const* x);
+/// calculate the enrichment function for fracture relative displacements
+/// Remarks:
+/// * branch/junction intersections of two fractures are supported in 2D
+///
+/// @param this_fracID the fracture ID
+/// @param frac_props fracture properties
+/// @param junction_props junction properties
+/// @param fracID_to_local a mapping table from a fracture ID to a local index
+/// 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::size_t this_fracID,
+    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);
 
 }  // namespace LIE
 }  // namespace ProcessLib