Merge pull request #1455 from norihiro-w/sd-lie-common

[SD-LIE] Common tools

ProcessLib/CMakeLists.txt

@@ -16,6 +16,8 @@ APPEND_SOURCE_FILES(SOURCES GroundwaterFlow)
 
 APPEND_SOURCE_FILES(SOURCES SmallDeformation)
 
+APPEND_SOURCE_FILES(SOURCES SmallDeformationWithLIE/Common)
+
 add_subdirectory(TES)
 APPEND_SOURCE_FILES(SOURCES TES)
 

ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp

@@ -9,6 +9,7 @@
 
 #include "CreateGroundwaterFlowProcess.h"
 
+#include "BaseLib/FileTools.h"
 #include "ProcessLib/Utils/ParseSecondaryVariables.h"
 #include "ProcessLib/Utils/ProcessUtils.h"
 #include "ProcessLib/CalculateSurfaceFlux/ParseCalculateSurfaceFluxData.h"

ProcessLib/SmallDeformationWithLIE/Common/FractureProperty.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_FRACTUREPROPERTY_H_
+#define PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_FRACTUREPROPERTY_H_
+
+#include <Eigen/Eigen>
+
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Node.h"
+
+#include "ProcessLib/Parameter/Parameter.h"
+
+#include "FractureProperty.h"
+#include "Utils.h"
+
+
+namespace ProcessLib
+{
+namespace SmallDeformationWithLIE
+{
+
+struct FractureProperty
+{
+    int mat_id = 0;
+    Eigen::Vector3d point_on_fracture;
+    Eigen::Vector3d normal_vector;
+    /// Rotation matrix from global to local coordinates
+    Eigen::MatrixXd R;
+    /// Initial aperture
+    ProcessLib::Parameter<double> const* aperture0 = nullptr;
+};
+
+/// 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,
+                                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.getNode(0)->getCoords()[j];
+    computeNormalVector(e, frac_prop.normal_vector);
+    frac_prop.R.resize(dim, dim);
+    computeRotationMatrix(frac_prop.normal_vector, dim, frac_prop.R);
+}
+
+}  // namespace SmallDeformationWithLIE
+}  // namespace ProcessLib
+
+#endif // PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_FRACTUREPROPERTY_H_

ProcessLib/SmallDeformationWithLIE/Common/HMatrixUtils.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_HMATRIXPOLICYTYPE_H_
+#define PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_HMATRIXPOLICYTYPE_H_
+
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+
+namespace ProcessLib
+{
+
+/// An implementation of H-Matrix policy using same matrix and vector types
+/// (fixed size or dynamic) as in the ShapeMatrixPolicyType.
+template <typename ShapeFunction, unsigned DisplacementDim>
+class HMatrixPolicyType
+{
+private:
+    /// Reusing the ShapeMatrixPolicy vector type.
+    template <int N>
+    using VectorType =
+        typename ShapeMatrixPolicyType<ShapeFunction,
+                                       DisplacementDim>::template VectorType<N>;
+
+    /// Reusing the ShapeMatrixPolicy matrix type.
+    template <int N, int M>
+    using MatrixType = typename ShapeMatrixPolicyType<
+        ShapeFunction, DisplacementDim>::template MatrixType<N, M>;
+
+    // Dimensions of specific H-matrix for n-points and displacement dimension.
+    static int const _number_of_dof = ShapeFunction::NPOINTS * DisplacementDim;
+
+public:
+    using StiffnessMatrixType = MatrixType<_number_of_dof, _number_of_dof>;
+    using NodalForceVectorType = VectorType<_number_of_dof>;
+
+    using HMatrixType = MatrixType<DisplacementDim, _number_of_dof>;
+
+    using ConstitutiveMatrixType = MatrixType<DisplacementDim, DisplacementDim>;
+    using ForceVectorType = VectorType<DisplacementDim>;
+};
+
+
+/// Fills a H-matrix based on given shape function
+template <int DisplacementDim, int NPOINTS, typename N_Type,
+          typename HMatrixType>
+void computeHMatrix(N_Type const& N, HMatrixType& H)
+{
+    static_assert(1 < DisplacementDim && DisplacementDim <= 3,
+                  "LinearHMatrix::computeHMatrix: DisplacementDim must be in "
+                  "range (1,3].");
+
+    H.setZero();
+
+    for (unsigned j=0; j<DisplacementDim; j++)
+        H.block(j, j*NPOINTS, 1, NPOINTS) = N;
+}
+
+}  // namespace ProcessLib
+
+#endif // PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_HMATRIXPOLICYTYPE_H_

ProcessLib/SmallDeformationWithLIE/Common/LevelSetFunction.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#include "LevelSetFunction.h"
+
+#ifndef Q_MOC_RUN // to avoid Qt4 bug, https://bugreports.qt.io/browse/QTBUG-22829
+#include <boost/math/special_functions/sign.hpp>
+#endif
+
+#include "FractureProperty.h"
+
+namespace
+{
+// Heaviside step function
+inline double Heaviside(double v)
+{
+    return (v < 0.0) ? 0.0 : 1.0;
+}
+
+} // no named namespace
+
+namespace ProcessLib
+{
+namespace SmallDeformationWithLIE
+{
+
+double calculateLevelSetFunction(
+        FractureProperty const& frac, double const* x_)
+{
+    Eigen::Map<Eigen::Vector3d const> x(x_, 3);
+    return Heaviside(
+                boost::math::sign(
+                    frac.normal_vector.dot(x - frac.point_on_fracture)));
+}
+
+} // SmallDeformationWithLIE
+} // ProcessLib

ProcessLib/SmallDeformationWithLIE/Common/LevelSetFunction.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#ifndef PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_LEVELSETFUNCTION_H_
+#define PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_LEVELSETFUNCTION_H_
+
+namespace ProcessLib
+{
+namespace SmallDeformationWithLIE
+{
+
+struct FractureProperty;
+
+/// 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);
+
+} // SmallDeformationWithLIE
+} // ProcessLib
+
+#endif // PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_LEVELSETFUNCTION_H_

ProcessLib/SmallDeformationWithLIE/Common/MeshUtils.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#include "MeshUtils.h"
+
+
+namespace ProcessLib
+{
+namespace SmallDeformationWithLIE
+{
+
+void getFractureMatrixDataInMesh(
+        MeshLib::Mesh const& mesh,
+        std::vector<MeshLib::Element*>& vec_matrix_elements,
+        std::vector<MeshLib::Element*>& vec_fracture_elements,
+        std::vector<MeshLib::Element*>& vec_fracture_matrix_elements,
+        std::vector<MeshLib::Node*>& vec_fracture_nodes
+        )
+{
+    // get vectors of matrix elements and fracture elements
+    vec_matrix_elements.reserve(mesh.getNumberOfElements());
+    for (MeshLib::Element* e : mesh.getElements())
+    {
+        if (e->getDimension() == mesh.getDimension())
+            vec_matrix_elements.push_back(e);
+        else
+            vec_fracture_elements.push_back(e);
+    }
+    DBUG("-> found total %d matrix elements and %d fracture elements",
+         vec_matrix_elements.size(), vec_fracture_elements.size());
+
+    // get a vector of fracture nodes
+    for (MeshLib::Element* e : vec_fracture_elements)
+    {
+        for (unsigned i=0; i<e->getNumberOfNodes(); i++)
+        {
+            vec_fracture_nodes.push_back(const_cast<MeshLib::Node*>(e->getNode(i)));
+        }
+    }
+    std::sort(vec_fracture_nodes.begin(), vec_fracture_nodes.end(),
+        [](MeshLib::Node* node1, MeshLib::Node* node2) { return (node1->getID() < node2->getID()); }
+        );
+    vec_fracture_nodes.erase(
+                std::unique(vec_fracture_nodes.begin(), vec_fracture_nodes.end()),
+                vec_fracture_nodes.end());
+    DBUG("-> found %d nodes on the fracture", vec_fracture_nodes.size());
+
+    // create a vector fracture elements and connected matrix elements,
+    // which are passed to a DoF table
+    // first, collect matrix elements
+    for (MeshLib::Element *e : vec_fracture_elements)
+    {
+        for (unsigned i=0; i<e->getNumberOfBaseNodes(); i++)
+        {
+            MeshLib::Node const* node = e->getNode(i);
+            for (unsigned j=0; j<node->getNumberOfElements(); j++)
+            {
+                // only matrix elements
+                if (node->getElement(j)->getDimension() == mesh.getDimension()-1)
+                    continue;
+                vec_fracture_matrix_elements.push_back(const_cast<MeshLib::Element*>(node->getElement(j)));
+            }
+        }
+    }
+    std::sort(vec_fracture_matrix_elements.begin(), vec_fracture_matrix_elements.end(),
+        [](MeshLib::Element* p1, MeshLib::Element* p2) { return (p1->getID() < p2->getID()); }
+        );
+    vec_fracture_matrix_elements.erase(
+                std::unique(vec_fracture_matrix_elements.begin(), vec_fracture_matrix_elements.end()),
+                vec_fracture_matrix_elements.end());
+
+    // second, append fracture elements
+    vec_fracture_matrix_elements.insert(
+                vec_fracture_matrix_elements.end(),
+                vec_fracture_elements.begin(), vec_fracture_elements.end());
+}
+
+}  // namespace SmallDeformationWithLIE
+}  // namespace ProcessLib

ProcessLib/SmallDeformationWithLIE/Common/MeshUtils.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#ifndef PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_MESHUTILS_H_
+#define PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_MESHUTILS_H_
+
+#include <vector>
+
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/Node.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationWithLIE
+{
+
+/**
+ * get data about fracture and matrix elements/nodes from a mesh
+ *
+ * @param mesh  A mesh which includes fracture elements, i.e. lower-dimensional elements.
+ * It is assumed that elements forming a fracture have a distinct material ID.
+ * @param vec_matrix_elements  a vector of matrix elements
+ * @param vec_fracture_elements  a vector of fracture elements
+ * @param vec_fracture_matrix_elements  a vector of fracture elements and matrix elements connecting to the fracture
+ * @param vec_fracture_nodes  a vector of fracture element nodes
+ */
+void getFractureMatrixDataInMesh(
+        MeshLib::Mesh const& mesh,
+        std::vector<MeshLib::Element*>& vec_matrix_elements,
+        std::vector<MeshLib::Element*>& vec_fracture_elements,
+        std::vector<MeshLib::Element*>& vec_fracture_matrix_elements,
+        std::vector<MeshLib::Node*>& vec_fracture_nodes
+        );
+
+}  // namespace SmallDeformationWithLIE
+}  // namespace ProcessLib
+
+#endif // PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_MESHUTILS_H_

ProcessLib/SmallDeformationWithLIE/Common/Utils.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#include "Utils.h"
+
+#include <cmath>
+
+#ifndef Q_MOC_RUN // to avoid Qt4 bug, https://bugreports.qt.io/browse/QTBUG-22829
+#include <boost/math/constants/constants.hpp>
+#endif
+
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Node.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationWithLIE
+{
+
+namespace
+{
+
+// Computed normal vector is oriented in the left direction of the given line element
+// such that computeRotationMatrix2D() returns the indentity matrix for line elements
+// parallel to a vector (1,0,0)
+void computeNormalVector2D(MeshLib::Element const& e, Eigen::Vector3d & normal_vector)
+{
+    assert(e.getGeomType() == MeshLib::MeshElemType::LINE);
+
+    auto v1 = (*e.getNode(1)) - (*e.getNode(0));
+    normal_vector[0] = -v1[1];
+    normal_vector[1] = v1[0];
+    normal_vector[2] = 0.0;
+    normal_vector.normalize();
+}
+
+// Compute a rotation matrix from global coordinates to local coordinates whose y axis
+// should be same as the given normal vector
+void computeRotationMatrix2D(Eigen::Vector3d const& n, Eigen::MatrixXd& matR)
+{
+    matR(0,0) = n[1];
+    matR(0,1) = -n[0];
+    matR(1,0) = n[0];
+    matR(1,1) = n[1];
+}
+
+} // no named namespace
+
+
+void computeNormalVector(MeshLib::Element const& e, Eigen::Vector3d & normal_vector)
+{
+    if (e.getDimension() == 1)
+        computeNormalVector2D(e, normal_vector);
+    else
+        OGS_FATAL("2D elements are not supported in computeNormalVector()");
+}
+
+void computeRotationMatrix(Eigen::Vector3d const& normal_vector, int dim, Eigen::MatrixXd &matR)
+{
+    if (dim==2)
+        computeRotationMatrix2D(normal_vector, matR);
+    else
+        OGS_FATAL("2D elements are not supported in computeNormalVector()");
+}
+
+}  // namespace SmallDeformationWithLIE
+}  // namespace ProcessLib

ProcessLib/SmallDeformationWithLIE/Common/Utils.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#ifndef PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_UTILS_H_
+#define PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_UTILS_H_
+
+#include <Eigen/Eigen>
+
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Node.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationWithLIE
+{
+
+/// compute a normal vector of the given element
+void computeNormalVector(MeshLib::Element const& e, Eigen::Vector3d & normal_vector);
+
+/// compute a rotation matrix from global to local coordinates based on the given normal vector
+void computeRotationMatrix(Eigen::Vector3d const& normal_vector, int dim, Eigen::MatrixXd &matR);
+
+/// compute physical coordinates from the given shape vector, i.e. from the natural coordinates
+template <typename Derived>
+MathLib::Point3d computePhysicalCoordinates(MeshLib::Element const&e, Eigen::MatrixBase<Derived> const& shape)
+{
+    MathLib::Point3d pt;
+    for (unsigned i=0; i<e.getNumberOfNodes(); i++)
+    {
+        MeshLib::Node const& node = *e.getNode(i);
+        for (unsigned j=0; j<3; j++)
+            pt[j] += shape[i]*node[j];
+    }
+    return pt;
+}
+
+}  // namespace SmallDeformationWithLIE
+}  // namespace ProcessLib
+
+#endif

Tests/ProcessLib/TestLIE.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <cmath>
+
+#include <gtest/gtest.h>
+
+#include <Eigen/Eigen>
+
+#include "MeshLib/Elements/Line.h"
+#include "MeshLib/Mesh.h"
+
+#include "ProcessLib/SmallDeformationWithLIE/Common/Utils.h"
+
+
+namespace
+{
+
+std::unique_ptr<MeshLib::Mesh> createLine(
+    std::array<double, 3> const& a, std::array<double, 3> const& b)
+{
+    MeshLib::Node** nodes = new MeshLib::Node*[2];
+    nodes[0] = new MeshLib::Node(a);
+    nodes[1] = new MeshLib::Node(b);
+    MeshLib::Element* e = new MeshLib::Line(nodes);
+
+    return std::unique_ptr<MeshLib::Mesh>(
+                new MeshLib::Mesh("",
+                                  std::vector<MeshLib::Node*>{nodes[0], nodes[1]},
+                                  std::vector<MeshLib::Element*>{e})
+                );
+}
+
+std::unique_ptr<MeshLib::Mesh> createX()
+{
+    return createLine({{-1.0, 0.0, 0.0}}, {{1.0,  0.0, 0.0}});
+}
+
+std::unique_ptr<MeshLib::Mesh> createY()
+{
+    return createLine({{0.0, -1.0, 0.0}}, {{0.0,  1.0, 0.0}});
+}
+
+std::unique_ptr<MeshLib::Mesh> createXY()
+{
+    // 45degree inclined
+    return createLine({{0.0, 0.0, 0.0}}, {{2./sqrt(2), 2./sqrt(2), 0.0}});
+}
+
+const double eps = std::numeric_limits<double>::epsilon();
+
+}
+
+TEST(LIE, rotationMatrixX)
+{
+    auto msh(createX());
+    auto e(msh->getElement(0));
+    Eigen::Vector3d nv;
+    ProcessLib::SmallDeformationWithLIE::computeNormalVector(*e, nv);
+    ASSERT_EQ(0., nv[0]);
+    ASSERT_EQ(1., nv[1]);
+    ASSERT_EQ(0., nv[2]);
+
+    Eigen::MatrixXd R(2,2);
+    ProcessLib::SmallDeformationWithLIE::computeRotationMatrix(nv, 2, R);
+
+    ASSERT_NEAR(1., R(0,0), eps);
+    ASSERT_NEAR(0., R(0,1), eps);
+    ASSERT_NEAR(0., R(1,0), eps);
+    ASSERT_NEAR(1., R(1,1), eps);
+}
+
+TEST(LIE, rotationMatrixY)
+{
+    auto msh(createY());
+    auto e(msh->getElement(0));
+    Eigen::Vector3d nv;
+    ProcessLib::SmallDeformationWithLIE::computeNormalVector(*e, nv);
+    ASSERT_EQ(-1., nv[0]);
+    ASSERT_EQ(0., nv[1]);
+    ASSERT_EQ(0., nv[2]);
+
+    Eigen::MatrixXd R(2,2);
+    ProcessLib::SmallDeformationWithLIE::computeRotationMatrix(nv, 2, R);
+
+    ASSERT_NEAR(0., R(0,0), eps);
+    ASSERT_NEAR(1., R(0,1), eps);
+    ASSERT_NEAR(-1., R(1,0), eps);
+    ASSERT_NEAR(0., R(1,1), eps);
+}
+
+TEST(LIE, rotationMatrixXY)
+{
+    auto msh(createXY());
+    auto e(msh->getElement(0));
+    Eigen::Vector3d nv;
+    ProcessLib::SmallDeformationWithLIE::computeNormalVector(*e, nv);
+    ASSERT_NEAR(-1./sqrt(2), nv[0], eps);
+    ASSERT_NEAR(1./sqrt(2), nv[1], eps);
+    ASSERT_EQ(0., nv[2]);
+
+    Eigen::MatrixXd R(2,2);
+    ProcessLib::SmallDeformationWithLIE::computeRotationMatrix(nv, 2, R);
+
+    ASSERT_NEAR(1./sqrt(2), R(0,0), eps);
+    ASSERT_NEAR(1./sqrt(2), R(0,1), eps);
+    ASSERT_NEAR(-1./sqrt(2), R(1,0), eps);
+    ASSERT_NEAR(1./sqrt(2), R(1,1), eps);
+}