diff --git a/BaseLib/Algorithm.h b/BaseLib/Algorithm.h index 57baff5af97..22b74c4ae03 100644 --- a/BaseLib/Algorithm.h +++ b/BaseLib/Algorithm.h @@ -12,6 +12,7 @@ #pragma once #include +#include #include #include #include @@ -217,11 +218,40 @@ void uniquePushBack(Container& container, container.push_back(element); } -template -inline bool contains(Container const& container, ValueType const& element) +template +bool contains(Container const& container, + typename Container::value_type const& element) +{ + return std::find(container.begin(), container.end(), element) != + container.end(); +} + +template +boost::optional 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 +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::max(); + } + return std::distance(container.begin(), it); +} } // namespace BaseLib diff --git a/ProcessLib/LIE/Common/BranchProperty.h b/ProcessLib/LIE/Common/BranchProperty.h index 4558f5bc52d..df37c4c5961 100644 --- a/ProcessLib/LIE/Common/BranchProperty.h +++ b/ProcessLib/LIE/Common/BranchProperty.h @@ -11,18 +11,30 @@ #include +#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 }; diff --git a/ProcessLib/LIE/Common/FractureProperty.h b/ProcessLib/LIE/Common/FractureProperty.h index a205aa068cb..6b7ee503921 100644 --- a/ProcessLib/LIE/Common/FractureProperty.h +++ b/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 const* aperture0 = nullptr; - std::vector> branches_master; - std::vector> branches_slave; + ProcessLib::Parameter const& aperture0; + std::vector branches_master; + std::vector branches_slave; + + FractureProperty(int const fracture_id_, int const material_id, + ProcessLib::Parameter const& initial_aperture) + : fracture_id(fracture_id_), + mat_id(material_id), + aperture0(initial_aperture) + { + } virtual ~FractureProperty() = default; }; struct FracturePropertyHM : public FractureProperty { - ProcessLib::Parameter const* specific_storage = nullptr; - ProcessLib::Parameter const* biot_coefficient = nullptr; + FracturePropertyHM(int const fracture_id_, int const material_id, + ProcessLib::Parameter const& initial_aperture, + ProcessLib::Parameter const& specific_storage_, + ProcessLib::Parameter const& biot_coefficient_) + : FractureProperty(fracture_id_, material_id, initial_aperture), + specific_storage(specific_storage_), + biot_coefficient(biot_coefficient_) + { + } + ProcessLib::Parameter const& specific_storage; + ProcessLib::Parameter 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 - 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 diff --git a/ProcessLib/LIE/Common/HMatrixUtils.h b/ProcessLib/LIE/Common/HMatrixUtils.h index fde4653a64b..6e9839fdecb 100644 --- a/ProcessLib/LIE/Common/HMatrixUtils.h +++ b/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 diff --git a/ProcessLib/LIE/Common/JunctionProperty.h b/ProcessLib/LIE/Common/JunctionProperty.h index 7d88ddc2ca2..c1b13110b69 100644 --- a/ProcessLib/LIE/Common/JunctionProperty.h +++ b/ProcessLib/LIE/Common/JunctionProperty.h @@ -9,9 +9,10 @@ #pragma once +#include #include -#include +#include "MeshLib/Node.h" namespace ProcessLib { @@ -19,10 +20,19 @@ namespace LIE { struct JunctionProperty final { - Eigen::Vector3d coords; - int junction_id; - int node_id; - std::array fracture_IDs; + JunctionProperty(int const junction_id_, + MeshLib::Node const& junctionNode, + std::array 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 const fracture_ids; + int const junction_id; EIGEN_MAKE_ALIGNED_OPERATOR_NEW }; diff --git a/ProcessLib/LIE/Common/LevelSetFunction.cpp b/ProcessLib/LIE/Common/LevelSetFunction.cpp index fa675fc7ee3..99fa1864bdf 100644 --- a/ProcessLib/LIE/Common/LevelSetFunction.cpp +++ b/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 u_global_enrichments( +std::vector uGlobalEnrichments( std::vector const& frac_props, std::vector const& junction_props, std::unordered_map const& fracID_to_local, @@ -51,7 +51,7 @@ std::vector u_global_enrichments( // pre-calculate levelsets for all fractures std::vector 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 enrichments(frac_props.size() + junction_props.size()); // fractures possibly with branches @@ -60,7 +60,7 @@ std::vector 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; } @@ -68,8 +68,8 @@ std::vector u_global_enrichments( 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; } @@ -77,7 +77,7 @@ std::vector u_global_enrichments( return enrichments; } -std::vector du_global_enrichments( +std::vector duGlobalEnrichments( std::size_t this_frac_id, std::vector const& frac_props, std::vector const& junction_props, @@ -89,7 +89,7 @@ std::vector du_global_enrichments( // pre-calculate levelsets for all fractures std::vector 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 enrichments(frac_props.size() + junction_props.size()); enrichments[this_frac_local_index] = 1.0; @@ -99,16 +99,16 @@ std::vector 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 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; diff --git a/ProcessLib/LIE/Common/LevelSetFunction.h b/ProcessLib/LIE/Common/LevelSetFunction.h index 0e090ad2bac..6fd8c97e5d8 100644 --- a/ProcessLib/LIE/Common/LevelSetFunction.h +++ b/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 u_global_enrichments( +std::vector uGlobalEnrichments( std::vector const& frac_props, std::vector const& junction_props, std::unordered_map const& fracID_to_local, @@ -47,7 +47,7 @@ std::vector u_global_enrichments( /// in frac_props /// @param x evaluating point coordiates /// @return a vector of enrichment values for fracture relative displacements -std::vector du_global_enrichments( +std::vector duGlobalEnrichments( std::size_t this_fracID, std::vector const& frac_props, std::vector const& junction_props, diff --git a/ProcessLib/LIE/Common/PostUtils.cpp b/ProcessLib/LIE/Common/PostUtils.cpp index 2cc727e897f..2d3ecd3498c 100644 --- a/ProcessLib/LIE/Common/PostUtils.cpp +++ b/ProcessLib/LIE/Common/PostUtils.cpp @@ -40,30 +40,6 @@ std::vector const& getMaterialIdsForNode( assert(itr != vec_nodeID_matIDs.end()); return itr->second; }; - -int getfrac2matid(std::vector const& fracmatids, int frac1matid) -{ - auto itr_mat2 = - std::find_if(fracmatids.begin(), fracmatids.end(), - [&](int matid) { return matid != frac1matid; }); - assert(itr_mat2 != fracmatids.end()); - return *itr_mat2; -}; - -int matid2fracid(std::vector const& vec, int matid) -{ - auto itr = std::find(vec.begin(), vec.end(), matid); - assert(itr != vec.end()); - return itr - vec.begin(); -}; - -unsigned getpos_in_ids(std::vector const& ids, int id) -{ - auto itr = std::find(ids.begin(), ids.end(), id); - assert(itr != ids.end()); - auto const pos = itr - ids.begin(); - return pos; -} } // namespace PostProcessTool::PostProcessTool( @@ -162,23 +138,27 @@ PostProcessTool::PostProcessTool( // branch nodes const auto& br_matids = getMaterialIdsForNode( vec_branch_nodeID_matIDs, node_id); - auto frac2_matid = getfrac2matid(br_matids, frac_matid); + auto const frac2_matid = BaseLib::findFirstNotEqualElement( + br_matids, frac_matid); + assert(frac2_matid); auto const frac2_id = - matid2fracid(vec_fracture_mat_IDs, frac2_matid); + BaseLib::findIndex(vec_fracture_mat_IDs, *frac2_matid); + assert(frac2_id != std::numeric_limits::max()); auto prop_levelset2 = org_mesh.getProperties().getPropertyVector( "levelset" + std::to_string(frac2_id + 1)); - unsigned pos = 0; + std::size_t pos = 0; if ((*prop_levelset2)[eid] == 0) { - // index of this frac - pos = getpos_in_ids(br_matids, frac_matid); + // index of this fracture + pos = BaseLib::findIndex(br_matids, frac_matid); } else if ((*prop_levelset2)[eid] == 1) { - // index of the other frac - pos = getpos_in_ids(br_matids, frac2_matid); + // index of the other fracture + pos = BaseLib::findIndex(br_matids, *frac2_matid); } + assert(pos != std::numeric_limits::max()); new_node_id = dup_newNodeIDs[pos]; } else @@ -186,9 +166,12 @@ PostProcessTool::PostProcessTool( // junction nodes const auto& jct_matids = getMaterialIdsForNode( vec_junction_nodeID_matIDs, node_id); - auto frac2_matid = getfrac2matid(jct_matids, frac_matid); + auto const frac2_matid = BaseLib::findFirstNotEqualElement( + jct_matids, frac_matid); + assert(frac2_matid); auto const frac2_id = - matid2fracid(vec_fracture_mat_IDs, frac2_matid); + BaseLib::findIndex(vec_fracture_mat_IDs, *frac2_matid); + assert(frac2_id != std::numeric_limits::max()); auto prop_levelset2 = org_mesh.getProperties().getPropertyVector( "levelset" + std::to_string(frac2_id + 1)); @@ -197,7 +180,9 @@ PostProcessTool::PostProcessTool( if ((*prop_levelset2)[eid] == 0) { // index of this frac - auto const pos = getpos_in_ids(jct_matids, frac_matid); + auto const pos = + BaseLib::findIndex(jct_matids, frac_matid); + assert(pos != std::numeric_limits::max()); new_node_id = dup_newNodeIDs[pos]; } else if ((*prop_levelset2)[eid] == 1) diff --git a/ProcessLib/LIE/Common/Utils.h b/ProcessLib/LIE/Common/Utils.h index bcd8090d839..92c7826272c 100644 --- a/ProcessLib/LIE/Common/Utils.h +++ b/ProcessLib/LIE/Common/Utils.h @@ -51,7 +51,9 @@ MathLib::Point3d computePhysicalCoordinates( { MeshLib::Node const& node = *e.getNode(i); for (unsigned j = 0; j < 3; j++) + { pt[j] += shape[i] * node[j]; + } } return pt; } diff --git a/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp index 67c8e20bb18..71c6d8964dd 100644 --- a/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp +++ b/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp @@ -97,8 +97,10 @@ std::unique_ptr createHydroMechanicsProcess( if (!use_monolithic_scheme) { if (pv_name == "pressure") + { p_process_variables.emplace_back( const_cast(*variable)); + } else { u_process_variables.emplace_back( @@ -113,7 +115,9 @@ std::unique_ptr createHydroMechanicsProcess( } if (p_u_process_variables.size() > 3 || u_process_variables.size() > 2) + { OGS_FATAL("Currently only one displacement jump is supported"); + } if (!use_monolithic_scheme) { @@ -121,7 +125,9 @@ std::unique_ptr createHydroMechanicsProcess( process_variables.push_back(std::move(u_process_variables)); } else + { process_variables.push_back(std::move(p_u_process_variables)); + } auto solid_constitutive_relations = MaterialLib::Solids::createConstitutiveRelations(parameters, @@ -251,26 +257,26 @@ std::unique_ptr createHydroMechanicsProcess( { auto& fracture_properties_config = *opt_fracture_properties_config; - frac_prop = std::make_unique(); - frac_prop->mat_id = + frac_prop = std::make_unique( + 0 /*fracture_id*/, //! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS_WITH_LIE__fracture_properties__material_id} - fracture_properties_config.getConfigParameter("material_id"); - frac_prop->aperture0 = &ProcessLib::findParameter( - //! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS_WITH_LIE__fracture_properties__initial_aperture} - fracture_properties_config, "initial_aperture", parameters, 1); - if (frac_prop->aperture0->isTimeDependent()) + fracture_properties_config.getConfigParameter("material_id"), + ProcessLib::findParameter( + //! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS_WITH_LIE__fracture_properties__initial_aperture} + fracture_properties_config, "initial_aperture", parameters, 1), + ProcessLib::findParameter( + //! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS_WITH_LIE__fracture_properties__specific_storage} + fracture_properties_config, "specific_storage", parameters, 1), + ProcessLib::findParameter( + //! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS_WITH_LIE__fracture_properties__biot_coefficient} + fracture_properties_config, "biot_coefficient", parameters, 1)); + if (frac_prop->aperture0.isTimeDependent()) { OGS_FATAL( "The initial aperture parameter '%s' must not be " "time-dependent.", - frac_prop->aperture0->name.c_str()); + frac_prop->aperture0.name.c_str()); } - frac_prop->specific_storage = &ProcessLib::findParameter( - //! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS_WITH_LIE__fracture_properties__specific_storage} - fracture_properties_config, "specific_storage", parameters, 1); - frac_prop->biot_coefficient = &ProcessLib::findParameter( - //! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS_WITH_LIE__fracture_properties__biot_coefficient} - fracture_properties_config, "biot_coefficient", parameters, 1); } // initial effective stress in matrix diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp index 0b1461e29e0..0c0397e9ae6 100644 --- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp +++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp @@ -83,7 +83,7 @@ HydroMechanicsProcess::HydroMechanicsProcess( // set fracture property assuming a fracture forms a straight line setFractureProperty(GlobalDim, *_vec_fracture_elements[0], - *_process_data.fracture_property.get()); + *_process_data.fracture_property); } // @@ -108,7 +108,9 @@ HydroMechanicsProcess::HydroMechanicsProcess( if (std::find(vec_fracture_mat_IDs.begin(), vec_fracture_mat_IDs.end(), matID) == vec_fracture_mat_IDs.end()) + { vec_p_inactive_matIDs.push_back(matID); + } } _process_data.p_element_status = std::make_unique(&mesh, @@ -300,13 +302,15 @@ void HydroMechanicsProcess::initializeConcreteProcess( for (MeshLib::Element const* e : _mesh.getElements()) { if (e->getDimension() < GlobalDim) + { continue; + } std::vector fracture_props( {_process_data.fracture_property.get()}); std::vector junction_props; std::unordered_map fracID_to_local({{0, 0}}); - std::vector levelsets = u_global_enrichments( + std::vector levelsets = uGlobalEnrichments( fracture_props, junction_props, fracID_to_local, Eigen::Vector3d(e->getCenterOfGravity().getCoords())); (*mesh_prop_levelset)[e->getID()] = levelsets[0]; @@ -332,16 +336,20 @@ void HydroMechanicsProcess::initializeConcreteProcess( { OGS_FATAL("Could not access MaterialIDs property from mesh."); } - auto frac = _process_data.fracture_property.get(); + auto const& frac = _process_data.fracture_property; for (MeshLib::Element const* e : _mesh.getElements()) { if (e->getDimension() == GlobalDim) + { continue; + } if ((*mesh_prop_matid)[e->getID()] != frac->mat_id) + { continue; + } ProcessLib::SpatialPosition x; x.setElementID(e->getID()); - (*mesh_prop_b)[e->getID()] = (*frac->aperture0)(0, x)[0]; + (*mesh_prop_b)[e->getID()] = frac->aperture0(0, x)[0]; } _process_data.mesh_prop_b = mesh_prop_b; @@ -492,7 +500,7 @@ void HydroMechanicsProcess::computeSecondaryVariableConcrete( // compute nodal w and aperture auto const& R = _process_data.fracture_property->R; - auto* const b0 = _process_data.fracture_property->aperture0; + auto const& b0 = _process_data.fracture_property->aperture0; MeshLib::PropertyVector& vec_w = *_process_data.mesh_prop_nodal_w; MeshLib::PropertyVector& vec_b = *_process_data.mesh_prop_nodal_b; @@ -500,12 +508,14 @@ void HydroMechanicsProcess::computeSecondaryVariableConcrete( double const w_n) { // skip aperture computation for element-wise defined b0 because there // are jumps on the nodes between the element's values. - if (dynamic_cast const*>(b0)) + if (dynamic_cast const*>(&b0)) + { return std::numeric_limits::quiet_NaN(); + } ProcessLib::SpatialPosition x; x.setNodeID(node_id); - return w_n + (*b0)(/*time independent*/ 0, x)[0]; + return w_n + b0(/*time independent*/ 0, x)[0]; }; Eigen::VectorXd g(GlobalDim), w(GlobalDim); @@ -514,11 +524,15 @@ void HydroMechanicsProcess::computeSecondaryVariableConcrete( auto const node_id = node->getID(); g.setZero(); for (int k = 0; k < GlobalDim; k++) + { g[k] = mesh_prop_g[node_id * GlobalDim + k]; + } w.noalias() = R * g; for (int k = 0; k < GlobalDim; k++) + { vec_w[node_id * GlobalDim + k] = w[k]; + } vec_b[node_id] = compute_nodal_aperture(node_id, w[GlobalDim - 1]); } diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h index 4b8d396744e..a8275f288a7 100644 --- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h +++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h @@ -66,7 +66,7 @@ HydroMechanicsLocalAssemblerFracturegetNodalValuesOnElement( + aperture0_node_values = frac_prop.aperture0.getNodalValuesOnElement( e, /*time independent*/ 0); SpatialPosition x_position; @@ -213,9 +213,9 @@ void HydroMechanicsLocalAssemblerFracture( i, i * displacement_size / GlobalDim) .noalias() = ip_data.N_u; + } ip_data.N_p = sm_p.N; ip_data.dNdx_p = sm_p.dNdx; @@ -243,7 +245,9 @@ void HydroMechanicsLocalAssemblerMatrix C; std::tie(sigma_eff, state, C) = std::move(*solution); @@ -308,7 +312,9 @@ void HydroMechanicsLocalAssemblerMatrix(local_x).segment(pressure_index, pressure_size); if (_process_data.deactivate_matrix_in_flow) + { setPressureOfInactiveNodes(t, p); + } auto u = local_x.segment(displacement_index, displacement_size); computeSecondaryVariableConcreteWithBlockVectors(t, p, u); @@ -361,7 +367,9 @@ void HydroMechanicsLocalAssemblerMatrix C; std::tie(sigma_eff, state, C) = std::move(*solution); @@ -419,8 +427,10 @@ void HydroMechanicsLocalAssemblerMatrixisActiveNode(_element.getNode(i))) + { continue; + } x_position.setNodeID(_element.getNodeIndex(i)); auto const p0 = (*_process_data.p0)(t, x_position)[0]; p[i] = p0; @@ -459,7 +471,9 @@ void HydroMechanicsLocalAssemblerMatrix< { // only inactive nodes if (_process_data.p_element_status->isActiveNode(_element.getNode(i))) + { continue; + } p_dot[i] = 0; } } diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix.h index 35485a3e775..760339a4dda 100644 --- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix.h +++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix.h @@ -52,7 +52,9 @@ class HydroMechanicsLocalAssemblerMatrix double const /*delta_t*/) override { for (auto& data : _ip_data) + { data.pushBackState(); + } } Eigen::Map getShapeMatrix( diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrixNearFracture-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrixNearFracture-impl.h index 4ba25c8e435..e1cec75f5d3 100644 --- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrixNearFracture-impl.h +++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrixNearFracture-impl.h @@ -79,7 +79,7 @@ void HydroMechanicsLocalAssemblerMatrixNearFracture< // levelset value of the element // remark: this assumes the levelset function is uniform within an element - std::vector levelsets = u_global_enrichments( + std::vector levelsets = uGlobalEnrichments( _fracture_props, _junction_props, _fracID_to_local, _e_center_coords); double const ele_levelset = levelsets[0]; // single fracture @@ -137,12 +137,14 @@ void HydroMechanicsLocalAssemblerMatrixNearFracture< auto p = const_cast(local_x).segment(pressure_index, pressure_size); if (_process_data.deactivate_matrix_in_flow) + { Base::setPressureOfInactiveNodes(t, p); + } auto u = local_x.segment(displacement_index, displacement_size); // levelset value of the element // remark: this assumes the levelset function is uniform within an element - std::vector levelsets = u_global_enrichments( + std::vector levelsets = uGlobalEnrichments( _fracture_props, _junction_props, _fracID_to_local, _e_center_coords); double const ele_levelset = levelsets[0]; // single fracture diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/LocalDataInitializer.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/LocalDataInitializer.h index d6f9b2e4e68..cf342feecfa 100644 --- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/LocalDataInitializer.h +++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/LocalDataInitializer.h @@ -141,10 +141,12 @@ class LocalDataInitializer final : _dof_table(dof_table) { if (shapefunction_order != 2) + { OGS_FATAL( "The given shape function order %d is not supported.\nOnly " "shape functions of order 2 are supported.", shapefunction_order); + } // /// Quads and Hexahedra /////////////////////////////////// #if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 && \ @@ -205,7 +207,9 @@ class LocalDataInitializer final std::vector involved_varIDs; // including deactived elements involved_varIDs.reserve(varIDs.size() + 1); if (isPressureDeactivated) + { involved_varIDs.push_back(0); // always pressure come in + } involved_varIDs.insert(involved_varIDs.end(), varIDs.begin(), varIDs.end()); @@ -219,8 +223,10 @@ class LocalDataInitializer final mesh_item.getNumberOfBaseNodes()); // pressure auto const max_varID = *std::max_element(varIDs.begin(), varIDs.end()); for (int i = 1; i < max_varID + 1; i++) + { vec_n_element_nodes.push_back( mesh_item.getNumberOfNodes()); // displacements + } unsigned local_id = 0; unsigned dof_id = 0; @@ -242,7 +248,9 @@ class LocalDataInitializer final auto global_index = _dof_table.getGlobalIndex(l, var_id, var_comp_id); if (global_index != NumLib::MeshComponentMap::nop) + { dofIndex_to_localIndex[dof_id++] = local_id; + } local_id++; } } diff --git a/ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.cpp b/ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.cpp index 91026a36618..bc0f18daee8 100644 --- a/ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.cpp +++ b/ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.cpp @@ -60,7 +60,9 @@ std::unique_ptr createSmallDeformationProcess( "'displacement_junctionN'"); } if (pv_name.find("displacement_jump") == 0) + { n_var_du++; + } auto variable = std::find_if(variables.cbegin(), variables.cend(), [&pv_name](ProcessVariable const& v) { @@ -146,25 +148,22 @@ std::unique_ptr createSmallDeformationProcess( } // Fracture properties - std::vector> vec_fracture_property; + std::vector fracture_properties; for ( auto fracture_properties_config : //! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION_WITH_LIE__fracture_properties} config.getConfigSubtreeList("fracture_properties")) { - auto& para_b0 = ProcessLib::findParameter( - //! \ogs_file_param_special{prj__processes__process__SMALL_DEFORMATION_WITH_LIE__fracture_properties__initial_aperture} - fracture_properties_config, "initial_aperture", parameters, 1); - auto frac_prop(new FractureProperty()); - frac_prop->fracture_id = vec_fracture_property.size(); - frac_prop->mat_id = + fracture_properties.emplace_back( + fracture_properties.size(), //! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION_WITH_LIE__fracture_properties__material_id} - fracture_properties_config.getConfigParameter("material_id"); - frac_prop->aperture0 = ¶_b0; - vec_fracture_property.emplace_back(frac_prop); + fracture_properties_config.getConfigParameter("material_id"), + ProcessLib::findParameter( + //! \ogs_file_param_special{prj__processes__process__SMALL_DEFORMATION_WITH_LIE__fracture_properties__initial_aperture} + fracture_properties_config, "initial_aperture", parameters, 1)); } - if (n_var_du < vec_fracture_property.size()) + if (n_var_du < fracture_properties.size()) { OGS_FATAL( "The number of displacement jumps and the number of " @@ -180,7 +179,7 @@ std::unique_ptr createSmallDeformationProcess( SmallDeformationProcessData process_data( materialIDs(mesh), std::move(solid_constitutive_relations), - std::move(fracture_model), std::move(vec_fracture_property), + std::move(fracture_model), std::move(fracture_properties), reference_temperature); SecondaryVariableCollection secondary_variables; diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h index 70478e9a465..3855fc685ca 100644 --- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h +++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h @@ -57,18 +57,16 @@ SmallDeformationLocalAssemblerFractureaperture0)(0, x_position)[0]; + ip_data._aperture0 = _fracture_property->aperture0(0, x_position)[0]; ip_data._aperture_prev = ip_data._aperture0; _secondary_data.N[ip] = sm.N; @@ -177,12 +175,6 @@ void SmallDeformationLocalAssemblerFracture< SpatialPosition x_position; x_position.setElementID(_element.getID()); - std::size_t this_frac_local_index = 0; - for (; this_frac_local_index < _fracture_props.size(); - this_frac_local_index++) - if (_fracture_props[this_frac_local_index] == _fracture_property) - break; - for (unsigned ip = 0; ip < n_integration_points; ip++) { x_position.setIntegrationPoint(ip); @@ -201,7 +193,7 @@ void SmallDeformationLocalAssemblerFracture< Eigen::Vector3d const ip_physical_coords( computePhysicalCoordinates(_element, N).getCoords()); - std::vector const levelsets(du_global_enrichments( + std::vector const levelsets(duGlobalEnrichments( _fracture_property->fracture_id, _fracture_props, _junction_props, _fracID_to_local, ip_physical_coords)); @@ -283,12 +275,6 @@ void SmallDeformationLocalAssemblerFracture const levelsets(du_global_enrichments( + std::vector const levelsets(duGlobalEnrichments( _fracture_property->fracture_id, _fracture_props, _junction_props, _fracID_to_local, ip_physical_coords)); @@ -326,10 +312,12 @@ void SmallDeformationLocalAssemblerFracture const levelsets( - u_global_enrichments(_fracture_props, _junction_props, - _fracID_to_local, ip_physical_coords)); + uGlobalEnrichments(_fracture_props, _junction_props, + _fracID_to_local, ip_physical_coords)); // u = u^hat + sum_i(enrich^br_i(x) * [u]_i) + sum_i(enrich^junc_i(x) * // [u]_i) diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp index e2c8332823f..dc2c9f5232c 100644 --- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp +++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp @@ -55,13 +55,13 @@ SmallDeformationProcess::SmallDeformationProcess( vec_junction_nodeID_matIDs); if (_vec_fracture_mat_IDs.size() != - _process_data._vec_fracture_property.size()) + _process_data.fracture_properties.size()) { OGS_FATAL( "The number of the given fracture properties (%d) are not " "consistent" " with the number of fracture groups in a mesh (%d).", - _process_data._vec_fracture_property.size(), + _process_data.fracture_properties.size(), _vec_fracture_mat_IDs.size()); } @@ -88,13 +88,13 @@ SmallDeformationProcess::SmallDeformationProcess( } // set fracture property - for (auto& fracture_prop : _process_data._vec_fracture_property) + for (auto& fracture_prop : _process_data.fracture_properties) { // based on the 1st element assuming a fracture forms a straight line setFractureProperty( DisplacementDim, - *_vec_fracture_elements[fracture_prop->fracture_id][0], - *fracture_prop.get()); + *_vec_fracture_elements[fracture_prop.fracture_id][0], + fracture_prop); } // set branches @@ -103,22 +103,19 @@ SmallDeformationProcess::SmallDeformationProcess( auto master_matId = vec_branch_nodeID_matIDs[i].second[0]; auto slave_matId = vec_branch_nodeID_matIDs[i].second[1]; auto& master_frac = - *_process_data._vec_fracture_property - [_process_data._map_materialID_to_fractureID[master_matId]]; + _process_data.fracture_properties + [_process_data._map_materialID_to_fractureID[master_matId]]; auto& slave_frac = - *_process_data._vec_fracture_property - [_process_data._map_materialID_to_fractureID[slave_matId]]; + _process_data.fracture_properties + [_process_data._map_materialID_to_fractureID[slave_matId]]; - auto* branch = createBranchProperty( + master_frac.branches_master.push_back(createBranchProperty( *mesh.getNode(vec_branch_nodeID_matIDs[i].first), master_frac, - slave_frac); + slave_frac)); - master_frac.branches_master.emplace_back(branch); - - auto* branch2 = createBranchProperty( + slave_frac.branches_slave.push_back(createBranchProperty( *mesh.getNode(vec_branch_nodeID_matIDs[i].first), master_frac, - slave_frac); - slave_frac.branches_slave.emplace_back(branch2); + slave_frac)); } // set junctions @@ -129,14 +126,15 @@ SmallDeformationProcess::SmallDeformationProcess( } for (std::size_t i = 0; i < vec_junction_nodeID_matIDs.size(); i++) { - std::vector fracIDs; - for (auto matid : vec_junction_nodeID_matIDs[i].second) - fracIDs.push_back( - _process_data._map_materialID_to_fractureID[matid]); - auto* junction = createJunctionProperty( - i, *mesh.getNode(vec_junction_nodeID_matIDs[i].first), fracIDs); - - _process_data._vec_junction_property.emplace_back(junction); + auto const& material_ids = vec_junction_nodeID_matIDs[i].second; + assert(material_ids.size() == 2); + std::array fracture_ids{ + _process_data._map_materialID_to_fractureID[material_ids[0]], + _process_data._map_materialID_to_fractureID[material_ids[1]]}; + + _process_data.junction_properties.emplace_back( + i, *mesh.getNode(vec_junction_nodeID_matIDs[i].first), + fracture_ids); } // create a table of connected junction IDs for each element @@ -418,7 +416,9 @@ void SmallDeformationProcess::initializeConcreteProcess( for (MeshLib::Element const* e : _mesh.getElements()) { if (e->getDimension() < DisplacementDim) + { continue; + } Eigen::Vector3d const pt(e->getCenterOfGravity().getCoords()); std::vector e_fracture_props; @@ -427,8 +427,7 @@ void SmallDeformationProcess::initializeConcreteProcess( for (auto fid : _process_data._vec_ele_connected_fractureIDs[e->getID()]) { - e_fracture_props.push_back( - _process_data._vec_fracture_property[fid].get()); + e_fracture_props.push_back(&_process_data.fracture_properties[fid]); e_fracID_to_local.insert({fid, tmpi++}); } std::vector e_junction_props; @@ -437,11 +436,10 @@ void SmallDeformationProcess::initializeConcreteProcess( for (auto fid : _process_data._vec_ele_connected_junctionIDs[e->getID()]) { - e_junction_props.push_back( - _process_data._vec_junction_property[fid].get()); + e_junction_props.push_back(&_process_data.junction_properties[fid]); e_juncID_to_local.insert({fid, tmpi++}); } - std::vector const levelsets(u_global_enrichments( + std::vector const levelsets(uGlobalEnrichments( e_fracture_props, e_junction_props, e_fracID_to_local, pt)); for (unsigned i = 0; i < e_fracture_props.size(); i++) @@ -460,7 +458,7 @@ void SmallDeformationProcess::initializeConcreteProcess( const_cast(mesh), "levelset" + std::to_string(e_junction_props[i]->junction_id + 1 + - _process_data._vec_fracture_property.size()), + _process_data.fracture_properties.size()), MeshLib::MeshItemType::Cell, 1); mesh_prop_levelset->resize(mesh.getNumberOfElements()); (*mesh_prop_levelset)[e->getID()] = @@ -485,7 +483,7 @@ void SmallDeformationProcess::initializeConcreteProcess( MeshLib::MeshItemType::Cell, 1); mesh_prop_b->resize(mesh.getNumberOfElements()); auto const& mesh_prop_matid = *_process_data._mesh_prop_materialIDs; - for (auto const& fracture_prop : _process_data._vec_fracture_property) + for (auto const& fracture_prop : _process_data.fracture_properties) { for (MeshLib::Element const* e : _mesh.getElements()) { @@ -493,13 +491,13 @@ void SmallDeformationProcess::initializeConcreteProcess( { continue; } - if (mesh_prop_matid[e->getID()] != fracture_prop->mat_id) + if (mesh_prop_matid[e->getID()] != fracture_prop.mat_id) { continue; } ProcessLib::SpatialPosition x; x.setElementID(e->getID()); - (*mesh_prop_b)[e->getID()] = (*fracture_prop->aperture0)(0, x)[0]; + (*mesh_prop_b)[e->getID()] = fracture_prop.aperture0(0, x)[0]; } } _process_data._mesh_prop_b = mesh_prop_b; diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcessData.h b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcessData.h index b7e369e03d2..4f63353afd6 100644 --- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcessData.h +++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcessData.h @@ -42,12 +42,12 @@ struct SmallDeformationProcessData std::unique_ptr< MaterialLib::Fracture::FractureModelBase>&& fracture_model, - std::vector>&& vec_fracture_prop, + std::vector&& fracture_properties_, double const reference_temperature) : material_ids(material_ids_), solid_materials{std::move(solid_materials_)}, _fracture_model{std::move(fracture_model)}, - _vec_fracture_property(std::move(vec_fracture_prop)), + fracture_properties(std::move(fracture_properties_)), _reference_temperature(reference_temperature) { } @@ -73,8 +73,8 @@ struct SmallDeformationProcessData std::unique_ptr> _fracture_model; - std::vector> _vec_fracture_property; - std::vector> _vec_junction_property; + std::vector fracture_properties; + std::vector junction_properties; MeshLib::PropertyVector const* _mesh_prop_materialIDs = nullptr; std::vector _map_materialID_to_fractureID;