CreateOrthotropicEmbeddedFracturePermeability.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#include <Eigen/Geometry>
 
#include "BaseLib/ConfigTree.h"
#include "OrthotropicEmbeddedFracturePermeability.h"
#include "ParameterLib/Utils.h"
 
namespace MaterialPropertyLib
{
std::unique_ptr<Property> createOrthotropicEmbeddedFracturePermeability(
    int const geometry_dimension, BaseLib::ConfigTree const& config,
    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters)
{
    if ((geometry_dimension != 2) && (geometry_dimension != 3))
    {
        OGS_FATAL(
            "The OrthotropicEmbeddedFracturePermeability is implemented only "
            "for 2D or 3D problems");
    }

    config.checkConfigParameter("type",
                                "OrthotropicEmbeddedFracturePermeability");
 
    // Second access for storage.
    auto property_name = config.peekConfigParameter<std::string>("name");
 
    DBUG("Create OrthotropicEmbeddedFracturePermeability medium property");
 
    auto const a_i =
        config.getConfigParameter<std::vector<double>>("mean_frac_distances");
    if (a_i.size() != 3)
    {
        OGS_FATAL(
            "The size of the mean fracture distances vector must be 3, but is "
            "{}.",
            a_i.size());
    }
 
    auto const e_i0 =
        config.getConfigParameter<std::vector<double>>("threshold_strains");
    if (e_i0.size() != 3)
    {
        OGS_FATAL(
            "The size of the mean threshold strains vector must be 3, but is "
            "{}.",
            e_i0.size());
    }
 
    auto const n =
        config.getConfigParameter<std::vector<double>>("fracture_normals");
    if (n.size() != 6)
    {
        OGS_FATAL(
            "The size of the fracture normals vector must be 6, but is {}.",
            n.size());
    }
    Eigen::Vector3d const n1 = Eigen::Vector3d({n[0], n[1], n[2]}).normalized();
    Eigen::Vector3d const n2 = Eigen::Vector3d({n[3], n[4], n[5]}).normalized();
 
    if (n1.dot(n2) > std::numeric_limits<double>::epsilon())
    {
        OGS_FATAL(
            "The given fracture normals are not orthogonal. Please provide two "
            "orthogonal fracture normals");
    }
 
    Eigen::Matrix3d const n_i =
        (Eigen::Matrix3d() << n1, n2, n1.cross(n2)).finished();
 
    std::string const intrinsic_permeability_param_name =
        config.getConfigParameter<std::string>("intrinsic_permeability");
 
    auto const& k = ParameterLib::findParameter<double>(
        intrinsic_permeability_param_name, parameters, 0, nullptr);
 
    std::string const fracture_rotation_xy_param_name =
        config.getConfigParameter<std::string>("fracture_rotation_xy");
 
    auto const& phi_xy = ParameterLib::findParameter<double>(
        fracture_rotation_xy_param_name, parameters, 0, nullptr);
 
    std::string const fracture_rotation_yz_param_name =
        config.getConfigParameter<std::string>("fracture_rotation_yz");
 
    auto const& phi_yz = ParameterLib::findParameter<double>(
        fracture_rotation_yz_param_name, parameters, 0, nullptr);
 
    auto const jf =
        config.getConfigParameter<double>("jacobian_factor", 0.);
 
    if (geometry_dimension == 2)
    {
        return std::make_unique<OrthotropicEmbeddedFracturePermeability<2>>(
            std::move(property_name), a_i, e_i0, n_i, k, phi_xy, phi_yz, jf);
    }
    return std::make_unique<OrthotropicEmbeddedFracturePermeability<3>>(
        std::move(property_name), a_i, e_i0, n_i, k, phi_xy, phi_yz, jf);
}
}  // namespace MaterialPropertyLib