From 2221c611cda8581e60bfd4e32272d9932a8f9e71 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Tue, 2 Dec 2025 11:09:28 +0100 Subject: [PATCH 01/18] remove unused imports --- .../tests/test_link_constraint.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/applications/StructuralMechanicsApplication/tests/test_link_constraint.py b/applications/StructuralMechanicsApplication/tests/test_link_constraint.py index 0222a7bf60f3..dcf7cd68cc5b 100644 --- a/applications/StructuralMechanicsApplication/tests/test_link_constraint.py +++ b/applications/StructuralMechanicsApplication/tests/test_link_constraint.py @@ -6,9 +6,6 @@ # --- STD Imports --- import pathlib -import contextlib -import os -from typing import Generator import math From 5bab33a44f5e7722ae77c6493a478a1bf8725fd1 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Tue, 2 Dec 2025 11:09:54 +0100 Subject: [PATCH 02/18] add process for pretension modeling --- .../insert_pretension_process.cpp | 507 ++++++++++++++++++ .../insert_pretension_process.hpp | 46 ++ 2 files changed, 553 insertions(+) create mode 100644 applications/StructuralMechanicsApplication/custom_processes/insert_pretension_process.cpp create mode 100644 applications/StructuralMechanicsApplication/custom_processes/insert_pretension_process.hpp diff --git a/applications/StructuralMechanicsApplication/custom_processes/insert_pretension_process.cpp b/applications/StructuralMechanicsApplication/custom_processes/insert_pretension_process.cpp new file mode 100644 index 000000000000..de62b2cea8dc --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_processes/insert_pretension_process.cpp @@ -0,0 +1,507 @@ +// KRATOS ___| | | | +// \___ \ __| __| | | __| __| | | __| _` | | +// | | | | | ( | | | | ( | | +// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS +// +// License: BSD License +// license: StructuralMechanicsApplication/license.txt +// +// Main authors: Máté Kelemen +// + +// --- Structural Includes --- +#include "custom_processes/insert_pretension_process.hpp" // InsertPretensionProcess + +// --- Core Includes --- +#include "includes/master_slave_constraint.h" // MasterSlaveConstraint +#include "includes/model_part.h" // ModelPart +#include "processes/find_global_nodal_elemental_neighbours_process.h" +#include "utilities/comparison.h" // Comparison +#include "includes/element.h" // Element +//#include "geometries/point_3d.h" // Point3D + +// --- STL Includes --- +#include // std::array +#include // std::stringstream + + +namespace Kratos { + + +//class PretensionElement : public Element { +// +//}; // class PretensionElement + + +struct InsertPretensionProcess::Impl { + ModelPart* mpModelPart; + std::string mImpositionName; + double mPretensionValue; +}; // struct InsertPretensionProcess::Impl + + +InsertPretensionProcess::InsertPretensionProcess(Model& rModel, Parameters Settings) + : Process(), + mpImpl(new Impl) { + Settings.ValidateAndAssignDefaults(this->GetDefaultParameters()); + + KRATOS_TRY + mpImpl->mpModelPart = &rModel.GetModelPart(Settings["model_part_name"].GetString()); + mpImpl->mImpositionName = Settings["imposition"].GetString(); + mpImpl->mPretensionValue = Settings["pretension_value"].GetDouble(); + KRATOS_CATCH("") + + const std::array valid_imposition_names { + std::string {"lagrange_element"}, + std::string {"penalty_element"}, + std::string {"constraint"} + }; + + if (std::find(valid_imposition_names.begin(), + valid_imposition_names.end(), + mpImpl->mImpositionName) == valid_imposition_names.end()) { + std::stringstream message; + message << "Invalid setting for \"imposition\". Options are: "; + for (const auto& r_imposition_name : valid_imposition_names) + message << '"' << r_imposition_name << "\" "; + KRATOS_ERROR << message.str(); + } + + // This process expects the pre-tension surface to be defined + // by a set of Conditions in the input ModelPart. It requires + // at least one Condition. + KRATOS_ERROR_IF(mpImpl->mpModelPart->Conditions().empty()) + << "InsertPretensionProcess requires at least one Condition in the provided ModelPart, " + << "but there aren't any in '" << mpImpl->mpModelPart->Name() << "'."; +} + + +InsertPretensionProcess::InsertPretensionProcess(InsertPretensionProcess&&) noexcept = default; + + +InsertPretensionProcess::~InsertPretensionProcess() = default; + + +struct PretensionSurfacePartition { + array_1d normal; + std::vector positive_side, negative_side; +}; // struct SurfacePartition + + +void FindElementsAdjacentToConditions(ModelPart& rModelPart) { + using NeighborElements = decltype(NEIGHBOUR_ELEMENTS)::Type; + + // Make sure every node on the pre-tension surface is in the model part. + KRATOS_TRY + ModelPart::NodesContainerType nodes; + for (Condition& r_condition : rModelPart.Conditions()) { + auto& r_geometry = r_condition.GetGeometry(); + nodes.insert(r_geometry.begin(), r_geometry.end()); + } // for r_condition in rModelPart.Conditions + rModelPart.Nodes().insert(nodes); + KRATOS_CATCH("") + + KRATOS_TRY + // Find elements containing each node in the model part. + FindGlobalNodalElementalNeighboursProcess(rModelPart.GetRootModelPart()).Execute(); + + block_for_each( + rModelPart.Conditions(), + [] (Condition& r_condition) -> void { + NeighborElements neighbor_elements; + + // Collect all neighbor elements from the condition's nodes. + for (Node& r_node : r_condition.GetGeometry()) { + NeighborElements& r_containing_elements = r_node.GetValue(NEIGHBOUR_ELEMENTS); + for (unsigned i_element=0u; i_element bool { + return rp_other_element->Id() == rp_element->Id(); + }); + if (it_element == neighbor_elements.ptr_end()) { + neighbor_elements.push_back(rp_element); + } + } + } + + // Filter elements that don't contain all nodes of the condition. + neighbor_elements.erase(std::remove_if( + neighbor_elements.begin(), + neighbor_elements.end(), + [&r_condition] (const Element& r_element) -> bool { + for (const Node& r_condition_node : r_condition.GetGeometry()) { + const auto it_node = std::find_if( + r_element.GetGeometry().begin(), + r_element.GetGeometry().end(), + [&r_condition_node] (const Node& r_element_node) -> bool { + return r_element_node.Id() == r_condition_node.Id(); + }); + if (it_node == r_element.GetGeometry().end()) + return true; + } + return false; + } + ), neighbor_elements.end()); + + // Set neighbor elements on the condition. + r_condition.SetValue(NEIGHBOUR_ELEMENTS, neighbor_elements); + }); + KRATOS_CATCH("") +} + + +template +PretensionSurfacePartition PretensionPlanePartition(const ModelPart::ConditionsContainerType::iterator itConditionBegin, + const ModelPart::ConditionsContainerType::iterator itConditionEnd) { + PretensionSurfacePartition pretension_surface; + + // There are three valid possibilities here. + // 1) All surfaces are 0-dimensional (corners) + // => expecting exactly 1 condition + // => expecting exactly 2 connected elements, both + // must have line geometries + // => the normal is not well-defined, so it is chosen + // as the average of the connected elements' directions + // 2) All surfaces are 1-dimensional (edges) + // => expecting at least 1 condition + // => expecting at least 2 connected elements, all of which + // must be 2-dimensional + // => the normal is not well-defined, so it is chosen as the + // average of the rotated endpoints in 2D, in positive + // direction. This assumes coordinates in z-direction + // to vanish. This assumption is checked. + // 3) All surfaces are 2-dimensional (surfaces) + // => expecting at least 1 condition + // => expecting at least 2 connected elements, all + // of which must be 3-dimensional + // => normals are well-defined + // In all cases, at least one connected element is required on + // either side of the average pre-tension plane. + + // Sanity checks. + if constexpr (SurfaceDimension == 0u) { + KRATOS_ERROR_IF_NOT(std::distance(itConditionBegin, itConditionEnd) == 1) + << "Expecting exactly 1 condition on the 0-dimensional pre-tension surface , but got " + << std::distance(itConditionBegin, itConditionEnd) << "."; + } else if constexpr (SurfaceDimension == 1u) { + KRATOS_ERROR_IF_NOT(1 <= std::distance(itConditionBegin, itConditionEnd)) + << "Expecting at least 1 condition on the 1-dimensional pre-tension surface , but got " + << std::distance(itConditionBegin, itConditionEnd) << "."; + } else if constexpr (SurfaceDimension == 2u) { + KRATOS_ERROR_IF_NOT(1 <= std::distance(itConditionBegin, itConditionEnd)) + << "Expecting at least 1 condition on the 2-dimensional pre-tension surface , but got " + << std::distance(itConditionBegin, itConditionEnd) << "."; + } else { + static_assert(float(SurfaceDimension) == 0.5f, "Invalid surface dimension"); + } + + for (auto it_condition=itConditionBegin; it_condition!=itConditionEnd; ++it_condition) { + Condition& r_condition = *it_condition; + + // Sanity checks on the condition. + KRATOS_ERROR_IF_NOT(r_condition.GetGeometry().LocalSpaceDimension() == SurfaceDimension) + << "Expecting all geometries defining the pre-tension surface to be " << SurfaceDimension << " dimensional, " + << "but condition " << r_condition.Id() << " on geometry " << r_condition.GetGeometry().Id() << " (" << r_condition.GetGeometry().Name() << ") " + << "is " << r_condition.GetGeometry().LocalSpaceDimension() << "-dimensional."; + + KRATOS_ERROR_IF_NOT(r_condition.Has(NEIGHBOUR_ELEMENTS)) + << "Expecting condition " << r_condition.Id() << " to have exactly " + << "2 neighboring elements, but found none."; + auto& r_neighbor_elements = r_condition.GetValue(NEIGHBOUR_ELEMENTS); + if (r_neighbor_elements.size() != 2) { + std::stringstream message; + message << "Expecting condition " << r_condition.Id() << " to have exactly " + << "2 neighboring elements, but found " << r_neighbor_elements.size() << " ["; + for (const Element& r_element : r_neighbor_elements) message << r_element.Id() << " "; + message << "]."; + KRATOS_ERROR << message.str(); + } + + // Check whether the 2 connected elements lie on opposite sides of the surface. + const array_1d + first_connected_direction = r_neighbor_elements[0].GetGeometry().Center() - r_condition.GetGeometry().Center(), + second_connected_direction = r_neighbor_elements[1].GetGeometry().Center() - r_condition.GetGeometry().Center(); + const double direction_inner_product = std::inner_product( + first_connected_direction.begin(), + first_connected_direction.end(), + second_connected_direction.begin(), + 0.0); + KRATOS_ERROR_IF_NOT(direction_inner_product < 0.0) + << "Element " << r_neighbor_elements[0].Id() << " and " << r_neighbor_elements[1].Id() << " " + << "do not lie on opposite sides of condition " << r_condition.Id() << " " + << "defining a pre-tension surface."; + + if constexpr (SurfaceDimension == 0u) { + for (unsigned i_element=0u; i_elementGetGeometry()[0].Id(); + const array_1d surface_node_position = itConditionBegin->GetGeometry()[0]; + array_1d opposite_node_position; + if (r_element.GetGeometry()[0].Id() == surface_node_id) { + opposite_node_position = r_element.GetGeometry()[1]; + } else if (r_element.GetGeometry()[1].Id() == surface_node_id) { + opposite_node_position = r_element.GetGeometry()[0]; + } else { + KRATOS_ERROR << "Element " << r_element.Id() << " is assumed to be connected to node " + << surface_node_id << " but is not."; + } + + if (i_element % 2) { + pretension_surface.positive_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); + pretension_surface.normal = opposite_node_position - surface_node_position; + } else { + pretension_surface.negative_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); + } + } // for i_element in range(r_neighbor_elements.size()) + } /*if SurfaceDimension == 0*/ else if constexpr (SurfaceDimension == 1u) { + // Define the line condition's plane. + const array_1d surface_begin = r_condition.GetGeometry()[0], + surface_end = *(r_condition.GetGeometry().end() - 1); + KRATOS_ERROR_IF_NOT(Comparison::Equal(1e-12, 1e-9)(surface_begin[2], 0.0)) + << "Assuming that 1-dimensional pre-tension surfaces lie on the XY plane, but node " + << r_condition.GetGeometry()[0].Id() << " of condition " << r_condition.Id() << " " + << "has a Z-coordinate of " << surface_begin[2] << "."; + KRATOS_ERROR_IF_NOT(Comparison::Equal(1e-12, 1e-9)(surface_end[2], 0.0)) + << "Assuming that 1-dimensional pre-tension surfaces lie on the XY plane, but node " + << r_condition.GetGeometry()[1].Id() << " of condition " << r_condition.Id() << " " + << "has a Z-coordinate of " << surface_end[2] << "."; + const array_1d + surface_normal { + -(surface_end[1] - surface_begin[1]), + surface_end[0] - surface_begin[0], + 0.0}, + surface_center = r_condition.GetGeometry().Center(); + pretension_surface.normal += surface_normal; + + // Sort connected elements. + for (unsigned i_element=0u; i_element element_direction = r_element.GetGeometry().Center() - surface_center; + const double direction_inner_product = std::inner_product( + surface_normal.begin(), + surface_normal.end(), + element_direction.begin(), + 0.0); + if (0 < direction_inner_product) pretension_surface.negative_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); + else pretension_surface.positive_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); + } // for i_element in range(r_neighbor_elements.size()) + + KRATOS_ERROR_IF_NOT(pretension_surface.positive_side.size() == pretension_surface.negative_side.size()) + << "Elements " << r_neighbor_elements[0].Id() << " and " << r_neighbor_elements[1].Id() << " " + << "lie on the same side of condition " << r_condition.Id() << " that defines a pre-tension surface."; + } /*else if SurfaceDimension == 1*/ else if constexpr (SurfaceDimension == 2u) { + // Define the surface's plane. + const array_1d surface_center = r_condition.GetGeometry().Center(); + const array_1d surface_normal = r_condition.GetGeometry().Normal(surface_center); + pretension_surface.normal += surface_normal; + + // Sort connected elements. + for (unsigned i_element=0u; i_element element_direction = r_element.GetGeometry().Center() - surface_center; + const double direction_inner_product = std::inner_product( + surface_normal.begin(), + surface_normal.end(), + element_direction.begin(), + 0.0); + if (0 < direction_inner_product) pretension_surface.negative_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); + else pretension_surface.positive_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); + } // for i_element in range(r_neighbor_elements.size()) + + KRATOS_ERROR_IF_NOT(pretension_surface.positive_side.size() == pretension_surface.negative_side.size()) + << "Elements " << r_neighbor_elements[0].Id() << " and " << r_neighbor_elements[1].Id() << " " + << "lie on the same side of condition " << r_condition.Id() << " that defines a pre-tension surface."; + } /*else if SurfaceDimension == 2*/ + } // for r_condition in r_model_part.Conditions + + return pretension_surface; +} + + +void InsertPretensionProcess::ExecuteBeforeSolutionLoop() { + ModelPart& r_model_part = *mpImpl->mpModelPart; + + // Find elements connected to the pre-tension surface. + // => each node will have a NEIGHBOR_ELEMENTS variable + // in its non-historical storage that lists elements + // containing them. + KRATOS_TRY + FindElementsAdjacentToConditions(r_model_part); + KRATOS_CATCH("") + + // Sort connected elements into 2 categories, lying on either side of + // the pre-tension surface. + PretensionSurfacePartition pretension_surface; + KRATOS_TRY + // Find the dimension of the pre-tension surface, + // and make sure all conditions share it. + KRATOS_ERROR_IF(r_model_part.Conditions().empty()); + const int pre_tension_surface_dimension = r_model_part.Conditions().front().GetGeometry().LocalSpaceDimension(); + switch (pre_tension_surface_dimension) { + case 0u: { + pretension_surface = PretensionPlanePartition<0>(r_model_part.Conditions().begin(), r_model_part.Conditions().end()); + break; + } + case 1u: { + pretension_surface = PretensionPlanePartition<1>(r_model_part.Conditions().begin(), r_model_part.Conditions().end()); + break; + } + case 2u: { + pretension_surface = PretensionPlanePartition<2>(r_model_part.Conditions().begin(), r_model_part.Conditions().end()); + break; + } + default: KRATOS_ERROR << "Invalid pre-tension surface dimension: " << pre_tension_surface_dimension << "."; + } + KRATOS_CATCH("") + + // Duplicate nodes on the pre-tension surface and collect relevant displacement components. + std::unordered_map duplicated_nodes; + KRATOS_TRY + std::unordered_set surface_nodes; + + // Collect nodes on the surface. + for (Condition& r_condition : r_model_part.Conditions()) { + for (Node& r_node : r_condition.GetGeometry()) { + surface_nodes.insert(&r_node); + } // for r_node in r_condition.GetGeometry() + } // for r_condition in r_model_part.Conditions() + + // Find the first available node ID that can safely be + // used to insert new ones. + Node::IndexType node_id = r_model_part.GetRootModelPart().Nodes().back().Id() + 1; + + // Clone surface nodes and associate them with the original ones. + for (Node* p_node : surface_nodes) { + const auto emplace_result = duplicated_nodes.emplace( + p_node->Id(), + p_node->Clone(node_id++)); + r_model_part.AddNode(emplace_result.first->second); + } // for p_node in surface_nodes + + // Replace nodes on the negative side elements with duplicated ones. + for (auto& rp_element : pretension_surface.negative_side) { + for (auto itp_node=rp_element->GetGeometry().ptr_begin(); itp_node!=rp_element->GetGeometry().ptr_end(); ++itp_node) { + const auto it_duplicate_pair = duplicated_nodes.find((*itp_node)->Id()); + if (it_duplicate_pair != duplicated_nodes.end()) { + (*itp_node) = it_duplicate_pair->second; + } + } // for itp_node in rp_element->GetGeometry() + } // for rp_element in pretension_surface.negative_side + KRATOS_CATCH("") + + // Construct the constraint equation. + struct PretensionConstraintData { + std::vector::Pointer> dofs; + std::vector gradient; + double gap; + } constraint_data; + KRATOS_TRY + const std::array*,3> all_displacement_components { + &DISPLACEMENT_X, + &DISPLACEMENT_Y, + &DISPLACEMENT_Z + }; + + Condition::DofsVectorType dofs; + for (Condition& r_condition : r_model_part.Conditions()) { + r_condition.GetDofList(dofs, r_model_part.GetProcessInfo()); + for (Dof::Pointer p_negative_side_dof : dofs) { + const auto it_displacement_component = std::find_if( + all_displacement_components.begin(), + all_displacement_components.end(), + [p_negative_side_dof](const Variable* p_variable) { + return p_variable->Key() == p_negative_side_dof->GetVariable().Key(); + }); + if (it_displacement_component != all_displacement_components.end()) { + // Insert the negative side component. + const std::size_t i_displacement_component = std::distance(all_displacement_components.begin(), it_displacement_component); + + { + const double gradient_component = pretension_surface.normal[i_displacement_component]; + constraint_data.dofs.push_back(p_negative_side_dof); + constraint_data.gradient.push_back(gradient_component); + } + + // Find the duplicate node on the positive side. + { + const Node::IndexType negative_side_node_id = p_negative_side_dof->Id(); + const auto it_positive_side_node = duplicated_nodes.find(negative_side_node_id); + KRATOS_ERROR_IF(it_positive_side_node == duplicated_nodes.end()); + Node& r_positive_side_node = *it_positive_side_node->second; + const auto& r_positive_side_dofs = r_positive_side_node.GetDofs(); + const auto it_positive_side_dof = std::find_if( + r_positive_side_dofs.begin(), + r_positive_side_dofs.end(), + [it_displacement_component](const auto& rp_dof) { + return rp_dof->GetVariable().Key() == (*it_displacement_component)->Key(); + }); + KRATOS_ERROR_IF(it_positive_side_dof == r_positive_side_dofs.end()); + const double gradient_component = -pretension_surface.normal[i_displacement_component]; + constraint_data.dofs.push_back(it_positive_side_dof->get()); + constraint_data.gradient.push_back(gradient_component); + } + } // if p_dof in all_displacement_components + } // for p_negative_side_dof in dofs + } // for r_condition in r_model_part.Conditions + + constraint_data.gap = mpImpl->mPretensionValue; + KRATOS_CATCH("") + + // Insert the constraint. + const std::size_t constraint_id = r_model_part.MasterSlaveConstraints().empty() ? + 1ul : + r_model_part.MasterSlaveConstraints().back().Id() + 1; + + MasterSlaveConstraint::DofPointerVectorType masters, slaves; + MasterSlaveConstraint::MatrixType relation_matrix; + MasterSlaveConstraint::VectorType constants; + + slaves.push_back(constraint_data.dofs.back()); + masters.reserve(constraint_data.dofs.size() - 1); + std::copy( + constraint_data.dofs.begin(), + constraint_data.dofs.end() - 1, + std::back_inserter(masters)); + + relation_matrix.resize(1, constraint_data.dofs.size() - 1, false); + std::transform( + constraint_data.gradient.begin(), + constraint_data.gradient.end() - 1, + relation_matrix.begin2(), + [&constraint_data](double component){return -component / constraint_data.gradient.back();}); + + constants.resize(1, false); + constants[0] = -constraint_data.gap; + + KRATOS_ERROR_IF_NOT(KratosComponents::Has("LinearMasterSlaveConstraint")); + MasterSlaveConstraint::Pointer p_constraint; + + KRATOS_TRY + p_constraint = KratosComponents::Get("LinearMasterSlaveConstraint").Create( + constraint_id, + masters, + slaves, + relation_matrix, + constants); + r_model_part.AddMasterSlaveConstraint(p_constraint); + KRATOS_CATCH("") +} + + +const Parameters InsertPretensionProcess::GetDefaultParameters() const { + return Parameters(R"({ + "model_part_name" : "", + "pretension_value" : 0.0, + "imposition" : "lagrange_element" + })"); +} + + +} // namespace Kratos diff --git a/applications/StructuralMechanicsApplication/custom_processes/insert_pretension_process.hpp b/applications/StructuralMechanicsApplication/custom_processes/insert_pretension_process.hpp new file mode 100644 index 000000000000..cc2a3e48d09f --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_processes/insert_pretension_process.hpp @@ -0,0 +1,46 @@ +// KRATOS ___| | | | +// \___ \ __| __| | | __| __| | | __| _` | | +// | | | | | ( | | | | ( | | +// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS +// +// License: BSD License +// license: StructuralMechanicsApplication/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// --- Core Includes --- +#include "processes/process.h" // Process +#include "containers/model.h" // Model + + +namespace Kratos { + + +class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertPretensionProcess final : public Process { +public: + KRATOS_CLASS_POINTER_DEFINITION(InsertPretensionProcess); + + InsertPretensionProcess(Model& rModel, Parameters Settings); + + InsertPretensionProcess(InsertPretensionProcess&&) noexcept; + + ~InsertPretensionProcess(); + + void ExecuteBeforeSolutionLoop() override; + + const Parameters GetDefaultParameters() const override; + +private: + InsertPretensionProcess(const InsertPretensionProcess&) = delete; + + InsertPretensionProcess& operator=(const InsertPretensionProcess&) = delete; + + struct Impl; + std::unique_ptr mpImpl; +}; // class InsertPretensionProcess + + +} // namespace Kratos From 3241285692dc9503e83744ce591de88517d2d2d7 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Tue, 2 Dec 2025 11:10:17 +0100 Subject: [PATCH 03/18] add a python interface for pretensioning --- .../custom_python/add_custom_processes_to_python.cpp | 7 ++++++- .../python_scripts/insert_pretension_process.py | 12 ++++++++++++ 2 files changed, 18 insertions(+), 1 deletion(-) create mode 100644 applications/StructuralMechanicsApplication/python_scripts/insert_pretension_process.py diff --git a/applications/StructuralMechanicsApplication/custom_python/add_custom_processes_to_python.cpp b/applications/StructuralMechanicsApplication/custom_python/add_custom_processes_to_python.cpp index ccd4059b1692..4ad257ee8f89 100644 --- a/applications/StructuralMechanicsApplication/custom_python/add_custom_processes_to_python.cpp +++ b/applications/StructuralMechanicsApplication/custom_python/add_custom_processes_to_python.cpp @@ -29,6 +29,7 @@ #include "custom_processes/spr_error_process.h" #include "custom_processes/impose_rigid_movement_process.h" #include "custom_processes/impose_z_strain_process.h" +#include "custom_processes/insert_pretension_process.hpp" #include "custom_processes/distribute_load_on_surface_process.h" #include "custom_processes/set_moving_load_process.h" #include "custom_processes/set_cartesian_local_axes_process.h" @@ -99,6 +100,10 @@ void AddCustomProcessesToPython(pybind11::module& m) .def(py::init< ModelPart&, Parameters >()) ; + py::class_(m, "InsertPretensionProcess") + .def(py::init()) + ; + py::class_(m,"DistributeLoadOnSurfaceProcess") .def(py::init()); @@ -113,7 +118,7 @@ void AddCustomProcessesToPython(pybind11::module& m) py::class_(m,"SetSphericalLocalAxesProcess") .def(py::init()); - + py::class_(m,"SetAutomatedInitialVariableProcess") .def(py::init()); } diff --git a/applications/StructuralMechanicsApplication/python_scripts/insert_pretension_process.py b/applications/StructuralMechanicsApplication/python_scripts/insert_pretension_process.py new file mode 100644 index 000000000000..c5ea6edbfec5 --- /dev/null +++ b/applications/StructuralMechanicsApplication/python_scripts/insert_pretension_process.py @@ -0,0 +1,12 @@ +# --- Core Imports --- +import KratosMultiphysics + +# --- Structural Imports --- +from KratosMultiphysics import StructuralMechanicsApplication + + +def Factory(parameters: KratosMultiphysics.Parameters, + model: KratosMultiphysics.Model) -> KratosMultiphysics.Process: + return StructuralMechanicsApplication.InsertPretensionProcess( + model, + parameters["Parameters"]) From 45a3b9abc9f7a8be7cc2702b0aa14859bd899aa2 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Tue, 2 Dec 2025 11:10:33 +0100 Subject: [PATCH 04/18] add test stub for pretensioning --- .../tests/constraints/pretension_1d.json | 92 +++++++++++++++++++ .../tests/constraints/pretension_1d.mdpa | 46 ++++++++++ .../constraints/pretension_1d_materials.json | 20 ++++ .../tests/test_pretension.py | 24 +++++ 4 files changed, 182 insertions(+) create mode 100644 applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.json create mode 100644 applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.mdpa create mode 100644 applications/StructuralMechanicsApplication/tests/constraints/pretension_1d_materials.json create mode 100644 applications/StructuralMechanicsApplication/tests/test_pretension.py diff --git a/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.json b/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.json new file mode 100644 index 000000000000..21c74d374031 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.json @@ -0,0 +1,92 @@ +{ + "analysis_stage": "KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis", + "problem_data": { + "problem_name": "test_pretension_1d", + "parallel_type": "OpenMP", + "echo_level": 1, + "start_time": 0, + "end_time": 1 + }, + "solver_settings": { + "time_stepping": { + "time_step": 0.5 + }, + "solver_type": "Static", + "model_part_name": "root", + "domain_size": 2, + "echo_level": 0, + "analysis_type": "linear", + "move_mesh_flag" : false, + "model_import_settings": { + "input_type": "mdpa", + "input_filename" : "pretension_1d" + }, + "material_import_settings": { + "materials_filename": "pretension_1d_materials.json" + }, + "convergence_criterion": "displacement_criterion", + "displacement_relative_tolerance": 1e-4, + "displacement_absolute_tolerance": 1e-9, + "rotation_dofs" : true + }, + "processes": { + "constraints_process_list": [ + { + "python_module": "assign_vector_variable_process", + "kratos_module": "KratosMultiphysics", + "process_name": "AssignVectorVariableProcess", + "Parameters": { + "model_part_name": "root.dirichlet", + "variable_name": "DISPLACEMENT", + "interval": [0, "End"], + "constrained": [true, true, false], + "value": [0, 0, 0] + } + }, + { + "python_module" : "insert_pretension_process", + "kratos_module" : "KratosMultiphysics.StructuralMechanicsApplication", + "Parameters" : { + "model_part_name" : "root.pretension_surface", + "pretension_value" : 5e-1 + } + } + ], + "loads_process_list": [ + { + "python_module": "assign_vector_by_direction_to_condition_process", + "kratos_module": "KratosMultiphysics", + "process_name": "AssignVectorByDirectionToConditionProcess", + "Parameters": { + "model_part_name": "root.neumann", + "variable_name": "POINT_LOAD", + "interval": [ + 0, "End" + ], + "modulus": 0.0, + "direction": [-2, 7, 0] + } + } + ] + }, + "output_processes": { + "vtk_output": [ + { + "Parameters": { + "folder_name": "vtk", + "model_part_name": "root", + "nodal_data_value_variables": [], + "nodal_solution_step_data_variables": [ + "DISPLACEMENT", + "REACTION" + ], + "output_sub_model_parts": false, + "save_output_files_in_folder": true + }, + "kratos_module": "KratosMultiphysics", + "process_name": "VTKOutputProcess", + "python_module": "vtk_output_process" + } + ] + } +} diff --git a/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.mdpa b/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.mdpa new file mode 100644 index 000000000000..1c4188693421 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.mdpa @@ -0,0 +1,46 @@ +Begin Properties 0 +End Properties + +Begin Nodes + 1 0 0 0 + 2 1 0 0 + 3 2 0 0 + 4 0 2 0 + 5 1 2 0 + 6 2 2 0 + 7 1 1 0 +End Nodes + +Begin Elements CrBeamElement2D2N + 1 0 1 2 + 2 0 2 3 + 3 0 4 5 + 4 0 5 6 + 5 0 2 7 + 6 0 7 5 +End Elements + +Begin Conditions PointLoadCondition2D1N + 1 0 7 +End Conditions + +Begin SubModelPart neumann + Begin SubModelPartConditions + 1 + End SubModelPartConditions +End SubModelPart + +Begin SubModelPart dirichlet + Begin SubModelPartNodes + 1 + 3 + 4 + 6 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart pretension_surface + Begin SubModelPartConditions + 1 + End SubModelPartConditions +End SubModelPart diff --git a/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d_materials.json b/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d_materials.json new file mode 100644 index 000000000000..1c61b4120912 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d_materials.json @@ -0,0 +1,20 @@ +{ + "properties": [ + { + "model_part_name": "root", + "properties_id": 1, + "Material": { + "constitutive_law": { + "name": "BeamConstitutiveLaw" + }, + "Variables": { + "DENSITY": 1e0, + "YOUNG_MODULUS": 1e0, + "POISSON_RATIO": 3e-1, + "CROSS_AREA": 1e-3, + "I33": 1e-5 + } + } + } + ] +} \ No newline at end of file diff --git a/applications/StructuralMechanicsApplication/tests/test_pretension.py b/applications/StructuralMechanicsApplication/tests/test_pretension.py new file mode 100644 index 000000000000..cef9e70f547d --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/test_pretension.py @@ -0,0 +1,24 @@ +# --- Kratos Imports --- +import KratosMultiphysics +import KratosMultiphysics.KratosUnittest +from KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis import StructuralMechanicsAnalysis +from KratosMultiphysics.KratosUnittest import WorkFolderScope + +# --- STD Imports --- +import pathlib + + +class TestPretension(KratosMultiphysics.KratosUnittest.TestCase): + + def test_pretension_1d(self) -> None: + with WorkFolderScope("constraints", pathlib.Path(__file__).absolute()): + # Load config + with open("pretension_1d.json", "r") as project_parameters_file: + parameters = KratosMultiphysics.Parameters(project_parameters_file.read()) + model = KratosMultiphysics.Model() + analysis = StructuralMechanicsAnalysis(model, parameters) + analysis.Run() + + +if __name__ == "__main__": + KratosMultiphysics.KratosUnittest.main() From 2f91a62337c224ae864eba106b8bb65d8b903f0b Mon Sep 17 00:00:00 2001 From: matekelemen Date: Wed, 10 Dec 2025 14:00:21 +0100 Subject: [PATCH 05/18] make the pretension insertion process an operation --- .../CMakeLists.txt | 1 + .../insert_pretension_operation.cpp} | 22 ++++++++-------- .../insert_pretension_operation.hpp} | 20 +++++++-------- .../add_custom_operations_to_python.cpp | 16 ++++++++++++ .../add_custom_operations_to_python.hpp | 16 ++++++++++++ .../add_custom_processes_to_python.cpp | 5 ---- ...tructural_mechanics_python_application.cpp | 2 ++ .../insert_pretension_process.py | 25 +++++++++++++++++-- 8 files changed, 79 insertions(+), 28 deletions(-) rename applications/StructuralMechanicsApplication/{custom_processes/insert_pretension_process.cpp => custom_operations/insert_pretension_operation.cpp} (97%) rename applications/StructuralMechanicsApplication/{custom_processes/insert_pretension_process.hpp => custom_operations/insert_pretension_operation.hpp} (56%) create mode 100644 applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.cpp create mode 100644 applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.hpp diff --git a/applications/StructuralMechanicsApplication/CMakeLists.txt b/applications/StructuralMechanicsApplication/CMakeLists.txt index 7184d6f2e7dc..11d03ec8e3bd 100644 --- a/applications/StructuralMechanicsApplication/CMakeLists.txt +++ b/applications/StructuralMechanicsApplication/CMakeLists.txt @@ -16,6 +16,7 @@ file(GLOB_RECURSE KRATOS_STRUCTURAL_MECHANICS_APPLICATION_CORE ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/*.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_elements/*.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_io/*.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/custom_operations/*.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_processes/*.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_response_functions/*.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_strategies/*.cpp diff --git a/applications/StructuralMechanicsApplication/custom_processes/insert_pretension_process.cpp b/applications/StructuralMechanicsApplication/custom_operations/insert_pretension_operation.cpp similarity index 97% rename from applications/StructuralMechanicsApplication/custom_processes/insert_pretension_process.cpp rename to applications/StructuralMechanicsApplication/custom_operations/insert_pretension_operation.cpp index de62b2cea8dc..08e425fa3a79 100644 --- a/applications/StructuralMechanicsApplication/custom_processes/insert_pretension_process.cpp +++ b/applications/StructuralMechanicsApplication/custom_operations/insert_pretension_operation.cpp @@ -10,7 +10,7 @@ // // --- Structural Includes --- -#include "custom_processes/insert_pretension_process.hpp" // InsertPretensionProcess +#include "custom_operations/insert_pretension_operation.hpp" // InsertPretensionOperation // --- Core Includes --- #include "includes/master_slave_constraint.h" // MasterSlaveConstraint @@ -33,15 +33,15 @@ namespace Kratos { //}; // class PretensionElement -struct InsertPretensionProcess::Impl { +struct InsertPretensionOperation::Impl { ModelPart* mpModelPart; std::string mImpositionName; double mPretensionValue; -}; // struct InsertPretensionProcess::Impl +}; // struct InsertPretensionOperation::Impl -InsertPretensionProcess::InsertPretensionProcess(Model& rModel, Parameters Settings) - : Process(), +InsertPretensionOperation::InsertPretensionOperation(Model& rModel, Parameters Settings) + : Operation(), mpImpl(new Impl) { Settings.ValidateAndAssignDefaults(this->GetDefaultParameters()); @@ -67,19 +67,19 @@ InsertPretensionProcess::InsertPretensionProcess(Model& rModel, Parameters Setti KRATOS_ERROR << message.str(); } - // This process expects the pre-tension surface to be defined + // This operation expects the pre-tension surface to be defined // by a set of Conditions in the input ModelPart. It requires // at least one Condition. KRATOS_ERROR_IF(mpImpl->mpModelPart->Conditions().empty()) - << "InsertPretensionProcess requires at least one Condition in the provided ModelPart, " + << "InsertPretensionOperation requires at least one Condition in the provided ModelPart, " << "but there aren't any in '" << mpImpl->mpModelPart->Name() << "'."; } -InsertPretensionProcess::InsertPretensionProcess(InsertPretensionProcess&&) noexcept = default; +InsertPretensionOperation::InsertPretensionOperation(InsertPretensionOperation&&) noexcept = default; -InsertPretensionProcess::~InsertPretensionProcess() = default; +InsertPretensionOperation::~InsertPretensionOperation() = default; struct PretensionSurfacePartition { @@ -325,7 +325,7 @@ PretensionSurfacePartition PretensionPlanePartition(const ModelPart::ConditionsC } -void InsertPretensionProcess::ExecuteBeforeSolutionLoop() { +void InsertPretensionOperation::Execute() { ModelPart& r_model_part = *mpImpl->mpModelPart; // Find elements connected to the pre-tension surface. @@ -495,7 +495,7 @@ void InsertPretensionProcess::ExecuteBeforeSolutionLoop() { } -const Parameters InsertPretensionProcess::GetDefaultParameters() const { +const Parameters InsertPretensionOperation::GetDefaultParameters() const { return Parameters(R"({ "model_part_name" : "", "pretension_value" : 0.0, diff --git a/applications/StructuralMechanicsApplication/custom_processes/insert_pretension_process.hpp b/applications/StructuralMechanicsApplication/custom_operations/insert_pretension_operation.hpp similarity index 56% rename from applications/StructuralMechanicsApplication/custom_processes/insert_pretension_process.hpp rename to applications/StructuralMechanicsApplication/custom_operations/insert_pretension_operation.hpp index cc2a3e48d09f..73ada50485fa 100644 --- a/applications/StructuralMechanicsApplication/custom_processes/insert_pretension_process.hpp +++ b/applications/StructuralMechanicsApplication/custom_operations/insert_pretension_operation.hpp @@ -12,35 +12,35 @@ #pragma once // --- Core Includes --- -#include "processes/process.h" // Process +#include "operations/operation.h" // Operation #include "containers/model.h" // Model namespace Kratos { -class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertPretensionProcess final : public Process { +class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertPretensionOperation final : public Operation { public: - KRATOS_CLASS_POINTER_DEFINITION(InsertPretensionProcess); + KRATOS_CLASS_POINTER_DEFINITION(InsertPretensionOperation); - InsertPretensionProcess(Model& rModel, Parameters Settings); + InsertPretensionOperation(Model& rModel, Parameters Settings); - InsertPretensionProcess(InsertPretensionProcess&&) noexcept; + InsertPretensionOperation(InsertPretensionOperation&&) noexcept; - ~InsertPretensionProcess(); + ~InsertPretensionOperation(); - void ExecuteBeforeSolutionLoop() override; + void Execute() override; const Parameters GetDefaultParameters() const override; private: - InsertPretensionProcess(const InsertPretensionProcess&) = delete; + InsertPretensionOperation(const InsertPretensionOperation&) = delete; - InsertPretensionProcess& operator=(const InsertPretensionProcess&) = delete; + InsertPretensionOperation& operator=(const InsertPretensionOperation&) = delete; struct Impl; std::unique_ptr mpImpl; -}; // class InsertPretensionProcess +}; // class InsertPretensionOperation } // namespace Kratos diff --git a/applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.cpp b/applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.cpp new file mode 100644 index 000000000000..61551b8c3e4b --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.cpp @@ -0,0 +1,16 @@ +// --- Structural Includes --- +#include "custom_python/add_custom_operations_to_python.hpp" // AddCustomOperationsToPython +#include "custom_operations/insert_pretension_operation.hpp" // InsertPretensionOperation + + +namespace Kratos::Python { + + +void AddCustomOperationsToPython(pybind11::module& rModule) { + pybind11::class_(rModule, "InsertPretensionOperation") + .def(pybind11::init()) + ; +} + + +} // namespace Kratos::Python diff --git a/applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.hpp b/applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.hpp new file mode 100644 index 000000000000..7880d612da93 --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.hpp @@ -0,0 +1,16 @@ +#pragma once + +// --- External Includes --- +#include + +// --- Core Includes --- +#include "includes/define_python.h" + + +namespace Kratos::Python { + + +void AddCustomOperationsToPython(pybind11::module& rModule); + + +} // namespace Kratos::Python diff --git a/applications/StructuralMechanicsApplication/custom_python/add_custom_processes_to_python.cpp b/applications/StructuralMechanicsApplication/custom_python/add_custom_processes_to_python.cpp index 4ad257ee8f89..60a5dbcdcb28 100644 --- a/applications/StructuralMechanicsApplication/custom_python/add_custom_processes_to_python.cpp +++ b/applications/StructuralMechanicsApplication/custom_python/add_custom_processes_to_python.cpp @@ -29,7 +29,6 @@ #include "custom_processes/spr_error_process.h" #include "custom_processes/impose_rigid_movement_process.h" #include "custom_processes/impose_z_strain_process.h" -#include "custom_processes/insert_pretension_process.hpp" #include "custom_processes/distribute_load_on_surface_process.h" #include "custom_processes/set_moving_load_process.h" #include "custom_processes/set_cartesian_local_axes_process.h" @@ -100,10 +99,6 @@ void AddCustomProcessesToPython(pybind11::module& m) .def(py::init< ModelPart&, Parameters >()) ; - py::class_(m, "InsertPretensionProcess") - .def(py::init()) - ; - py::class_(m,"DistributeLoadOnSurfaceProcess") .def(py::init()); diff --git a/applications/StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp b/applications/StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp index 3b6fae63902d..c3536c6e828b 100644 --- a/applications/StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp +++ b/applications/StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp @@ -20,6 +20,7 @@ #include "structural_mechanics_application_variables.h" #include "structural_mechanics_application.h" #include "custom_python/add_custom_strategies_to_python.h" +#include "custom_python/add_custom_operations_to_python.hpp" #include "custom_python/add_custom_processes_to_python.h" #include "custom_python/add_custom_utilities_to_python.h" #include "custom_python/add_custom_constitutive_laws_to_python.h" @@ -39,6 +40,7 @@ PYBIND11_MODULE(KratosStructuralMechanicsApplication,m) ; AddCustomStrategiesToPython(m); + AddCustomOperationsToPython(m); AddCustomProcessesToPython(m); AddCustomUtilitiesToPython(m); AddCustomConstitutiveLawsToPython(m); diff --git a/applications/StructuralMechanicsApplication/python_scripts/insert_pretension_process.py b/applications/StructuralMechanicsApplication/python_scripts/insert_pretension_process.py index c5ea6edbfec5..6b9b27044397 100644 --- a/applications/StructuralMechanicsApplication/python_scripts/insert_pretension_process.py +++ b/applications/StructuralMechanicsApplication/python_scripts/insert_pretension_process.py @@ -2,11 +2,32 @@ import KratosMultiphysics # --- Structural Imports --- -from KratosMultiphysics import StructuralMechanicsApplication +from KratosMultiphysics import StructuralMechanicsApplication as StructuralMechanics + + +class PretensionProcess(KratosMultiphysics.Process): + + def __init__(self, + model: KratosMultiphysics.Model, + parameters: KratosMultiphysics.Parameters) -> None: + parameters.ValidateAndAssignDefaults(self.GetDefaultParameters()) + #self.__model_part = model.GetModelPart(parameters["model_part_name"].GetString()) + self.__insert_pretension_operation = StructuralMechanics.InsertPretensionOperation( + model, + parameters) + + def ExecuteBeforeSolutionLoop(self) -> None: + self.__insert_pretension_operation.Execute() + + def GetDefaultParameters(self) -> KratosMultiphysics.Parameters: + return KratosMultiphysics.Parameters(R"""{ + "model_part_name" : "", + "pretension_value" : 0.0 + }""") def Factory(parameters: KratosMultiphysics.Parameters, model: KratosMultiphysics.Model) -> KratosMultiphysics.Process: - return StructuralMechanicsApplication.InsertPretensionProcess( + return PretensionProcess( model, parameters["Parameters"]) From 8f75dfb52302e8b794ef5e63a2c852e92d4acd33 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Wed, 10 Dec 2025 14:19:52 +0100 Subject: [PATCH 06/18] keep passed linear solver --- .../p_multigrid_builder_and_solver.cpp | 33 ++++++++++++------- .../p_multigrid_builder_and_solver.hpp | 4 +-- 2 files changed, 23 insertions(+), 14 deletions(-) diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.cpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.cpp index 7d98e48cb71b..0c9d0615f52b 100644 --- a/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.cpp +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.cpp @@ -60,6 +60,8 @@ struct PMultigridBuilderAndSolver::Impl PGrid,TUblasDenseSpace> >> mMaybeHierarchy; + typename LinearSolverType::Pointer mpRootGridSolver; + std::unique_ptr mpDiagonalScaling; int mMaxIterations; @@ -79,6 +81,7 @@ struct PMultigridBuilderAndSolver::Impl mpMaybeModelPart(), mpConstraintAssembler(), mMaybeHierarchy(), + mpRootGridSolver(), mpDiagonalScaling(), mMaxIterations(0), mTolerance(std::numeric_limits::max()), @@ -232,12 +235,13 @@ struct PMultigridBuilderAndSolver::Impl { // Prepare and initialize members. KRATOS_TRY - if (mpInterface->GetLinearSolver().AdditionalPhysicalDataIsNeeded()) { - mpInterface->GetLinearSolver().ProvideAdditionalData(rLhs, - rSolution, - rRhs, - mpInterface->GetDofSet(), - rModelPart); + if (mpInterface->GetRootGridSolver().AdditionalPhysicalDataIsNeeded()) { + mpInterface->GetRootGridSolver().ProvideAdditionalData( + rLhs, + rSolution, + rRhs, + mpInterface->GetDofSet(), + rModelPart); } if (mMaybeHierarchy.has_value()) { @@ -980,7 +984,7 @@ void PMultigridBuilderAndSolver::AssignSettings(const Parameters << "\"" << solver_name << "\" is not a valid linear solver name in the registry. " << "Make sure you imported the application it is defined in and that the spelling is correct."; const auto& r_factory = SolverFactoryRegistry::Get(solver_name); - this->mpLinearSystemSolver = r_factory.Create(smoother_settings); + this->mpImpl->mpRootGridSolver = r_factory.Create(smoother_settings); // Construct the coarse hierarchy. const std::string coarse_build_precision = coarse_hierarchy_settings["precision"].Get(); @@ -1019,10 +1023,15 @@ void PMultigridBuilderAndSolver::AssignSettings(const Parameters const std::string solver_name = leaf_solver_settings["solver_type"].Get(); using SolverFactoryRegistry = KratosComponents>; KRATOS_ERROR_IF_NOT(SolverFactoryRegistry::Has(solver_name)) - << "\"" << solver_name << "\" is not a valid linear solver name in the registry. " - << "Make sure you imported the application it is defined in and that the spelling is correct."; + << "'" << solver_name << "' does not name a registered linear solver. Options are:\n" + << []() -> std::string { + std::stringstream message; + for (const auto& r_pair : SolverFactoryRegistry::GetComponents()) + message << "\t" << r_pair.first << "\n"; + return message.str(); + }(); const auto& r_factory = SolverFactoryRegistry::Get(solver_name); - this->mpLinearSystemSolver = r_factory.Create(leaf_solver_settings); + this->mpImpl->mpRootGridSolver = r_factory.Create(leaf_solver_settings); KRATOS_CATCH("") } @@ -1064,7 +1073,7 @@ Parameters PMultigridBuilderAndSolver::GetDefaultParameters() co "solver_type" : "amgcl" }, "constraint_imposition_settings" : { - "method" : "master_slave_elimination" + "method" : "master_slave" }, "coarse_hierarchy_settings" : {} })"); @@ -1088,7 +1097,7 @@ std::size_t PMultigridBuilderAndSolver::GetEquationSystemSize() } template -LinearSolver& PMultigridBuilderAndSolver::GetLinearSolver() noexcept +LinearSolver& PMultigridBuilderAndSolver::GetRootGridSolver() noexcept { return *Interface::mpLinearSystemSolver; } diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.hpp index 55b3251ea75d..1f4336afcef3 100644 --- a/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.hpp +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.hpp @@ -351,7 +351,7 @@ class KRATOS_API(KRATOS_CORE) PMultigridBuilderAndSolver /// "tolerance" : 1e-8, /// "verbosity" : 1, /// "constraint_imposition_settings" : { - /// "method" : "master_slave_elimination" + /// "method" : "master_slave" /// }, /// "coarse_hierarchy_settings" : {} /// "smoother_settings" : { @@ -426,7 +426,7 @@ class KRATOS_API(KRATOS_CORE) PMultigridBuilderAndSolver private: std::size_t GetEquationSystemSize() const noexcept; - LinearSolverType& GetLinearSolver() noexcept; + LinearSolverType& GetRootGridSolver() noexcept; PMultigridBuilderAndSolver(const PMultigridBuilderAndSolver& rOther) = delete; From e0762201cb49c20120907eeefabd7bf3a3cfec35 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Fri, 10 Apr 2026 16:34:56 +0200 Subject: [PATCH 07/18] fix ignored linear solver settings --- .../p_multigrid_builder_and_solver.cpp | 33 ++++++++++--------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.cpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.cpp index 432604031d78..b38363296e92 100644 --- a/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.cpp +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.cpp @@ -1068,6 +1068,7 @@ void PMultigridBuilderAndSolver::AssignSettings(const Parameters << "Make sure you imported the application it is defined in and that the spelling is correct."; const auto& r_factory = SolverFactoryRegistry::Get(solver_name); this->mpImpl->mpRootGridSolver = r_factory.Create(smoother_settings); + Interface::mpLinearSystemSolver = mpImpl->mpRootGridSolver; // Construct the coarse hierarchy. const std::string coarse_build_precision = coarse_hierarchy_settings["precision"].Get(); @@ -1102,19 +1103,21 @@ void PMultigridBuilderAndSolver::AssignSettings(const Parameters KRATOS_CATCH("") } /* if coarse_hierarchy_settings["max_depth"].Get() */ else { KRATOS_TRY - KRATOS_ERROR_IF_NOT(leaf_solver_settings.Has("solver_type")); - const std::string solver_name = leaf_solver_settings["solver_type"].Get(); - using SolverFactoryRegistry = KratosComponents>; - KRATOS_ERROR_IF_NOT(SolverFactoryRegistry::Has(solver_name)) - << "'" << solver_name << "' does not name a registered linear solver. Options are:\n" - << []() -> std::string { - std::stringstream message; - for (const auto& r_pair : SolverFactoryRegistry::GetComponents()) - message << "\t" << r_pair.first << "\n"; - return message.str(); - }(); - const auto& r_factory = SolverFactoryRegistry::Get(solver_name); - this->mpImpl->mpRootGridSolver = r_factory.Create(leaf_solver_settings); + if (leaf_solver_settings.Has("solver_type")) { + const std::string solver_name = leaf_solver_settings["solver_type"].Get(); + using SolverFactoryRegistry = KratosComponents>; + KRATOS_ERROR_IF_NOT(SolverFactoryRegistry::Has(solver_name)) + << "'" << solver_name << "' does not name a registered linear solver. Options are:\n" + << []() -> std::string { + std::stringstream message; + for (const auto& r_pair : SolverFactoryRegistry::GetComponents()) + message << "\t" << r_pair.first << "\n"; + return message.str(); + }(); + const auto& r_factory = SolverFactoryRegistry::Get(solver_name); + this->mpImpl->mpRootGridSolver = r_factory.Create(leaf_solver_settings); + Interface::mpLinearSystemSolver = mpImpl->mpRootGridSolver; + } KRATOS_CATCH("") } @@ -1154,9 +1157,7 @@ Parameters PMultigridBuilderAndSolver::GetDefaultParameters() co "smoother_settings" : { "solver_type" : "" }, - "linear_solver_settings" : { - "solver_type" : "amgcl" - }, + "linear_solver_settings" : {}, "constraint_imposition_settings" : { "method" : "master_slave" }, From 680b9694a8183a06a8c8449224b75fb8512aa978 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Wed, 22 Apr 2026 16:59:54 +0200 Subject: [PATCH 08/18] fix linear mfcs --- .../p_multigrid/linear_multifreedom_constraint.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.cpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.cpp index 5ed7ce5014af..c99584e98e65 100644 --- a/kratos/solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.cpp +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.cpp @@ -41,7 +41,7 @@ LinearMultifreedomConstraint::LinearMultifreedomConstraint(const IndexType Id, << "(" << rRelationMatrix.size1() << "x" << rRelationMatrix.size2() << ")"; // Store input data. - this->SetValue(GEOMETRIC_STIFFNESS_MATRIX, rRelationMatrix); + this->SetValue(CONSTITUTIVE_MATRIX, rRelationMatrix); this->SetValue(INTERNAL_FORCES_VECTOR, rConstraintGaps); } @@ -50,7 +50,7 @@ void LinearMultifreedomConstraint::CalculateLocalSystem(MatrixType& rRelationMat VectorType& rConstraintGaps, const ProcessInfo&) const { - rRelationMatrix = this->GetData().GetValue(GEOMETRIC_STIFFNESS_MATRIX); + rRelationMatrix = this->GetData().GetValue(CONSTITUTIVE_MATRIX); rConstraintGaps = this->GetData().GetValue(INTERNAL_FORCES_VECTOR); } From 30c4d6ffc76fcf4fbb0506038b43a4a0ee06aee8 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Thu, 23 Apr 2026 15:25:04 +0200 Subject: [PATCH 09/18] add generic pre-tensioning system --- .../insert_pre_tension_operation.cpp | 1188 +++++++++++++++++ .../insert_pre_tension_operation.hpp | 101 ++ .../insert_pretension_operation.cpp | 507 ------- .../insert_pretension_operation.hpp | 46 - .../add_custom_operations_to_python.cpp | 21 +- ...ss.py => dirichlet_pre_tension_process.py} | 11 +- .../neumann_pre_tension_process.py | 34 + .../tests/constraints/pre_tensioned_beam.json | 124 ++ ...ension_1d.mdpa => pre_tensioned_beam.mdpa} | 26 +- ...json => pre_tensioned_beam_materials.json} | 10 +- .../tests/constraints/pretension_1d.json | 92 -- ...test_pretension.py => test_pre_tension.py} | 4 +- 12 files changed, 1492 insertions(+), 672 deletions(-) create mode 100644 applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp create mode 100644 applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp delete mode 100644 applications/StructuralMechanicsApplication/custom_operations/insert_pretension_operation.cpp delete mode 100644 applications/StructuralMechanicsApplication/custom_operations/insert_pretension_operation.hpp rename applications/StructuralMechanicsApplication/python_scripts/{insert_pretension_process.py => dirichlet_pre_tension_process.py} (79%) create mode 100644 applications/StructuralMechanicsApplication/python_scripts/neumann_pre_tension_process.py create mode 100644 applications/StructuralMechanicsApplication/tests/constraints/pre_tensioned_beam.json rename applications/StructuralMechanicsApplication/tests/constraints/{pretension_1d.mdpa => pre_tensioned_beam.mdpa} (64%) rename applications/StructuralMechanicsApplication/tests/constraints/{pretension_1d_materials.json => pre_tensioned_beam_materials.json} (59%) delete mode 100644 applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.json rename applications/StructuralMechanicsApplication/tests/{test_pretension.py => test_pre_tension.py} (86%) diff --git a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp new file mode 100644 index 000000000000..218913c3c162 --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp @@ -0,0 +1,1188 @@ +// KRATOS ___| | | | +// \___ \ __| __| | | __| __| | | __| _` | | +// | | | | | ( | | | | ( | | +// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS +// +// License: BSD License +// license: StructuralMechanicsApplication/license.txt +// +// Main authors: Máté Kelemen +// + +// --- Structural Includes --- +#include "custom_operations/insert_pre_tension_operation.hpp" // InsertPreTensionOperation +#include "structural_mechanics_application_variables.h" // POINT_LOAD_X + +// --- Core Includes --- +#include "includes/master_slave_constraint.h" // MasterSlaveConstraint +#include "includes/model_part.h" // ModelPart +#include "processes/find_global_nodal_elemental_neighbours_process.h" +#include "utilities/comparison.h" // Comparison +#include "includes/element.h" // Element +#include "input_output/logger.h" // KRATOS_INFO +#include "solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.hpp" // LinearMultifreedomConstraint +#include "geometries/point_2d.h" // Point2D + +// --- STL Includes --- +#include // std::array +#include // std::stringstream + + +namespace Kratos { + + +struct InsertPreTensionOperation::Impl { + ModelPart* mpModelPart; + double mPretensionValue; + int mVerbosity; +}; // struct InsertPreTensionOperation::Impl + + +InsertPreTensionOperation::InsertPreTensionOperation(Model& rModel, Parameters Settings) + : Operation(), + mpImpl(new Impl) { + Settings.ValidateAndAssignDefaults(this->GetDefaultParameters()); + + KRATOS_TRY + mpImpl->mpModelPart = &rModel.GetModelPart(Settings["model_part_name"].GetString()); + mpImpl->mPretensionValue = Settings["magnitude"].GetDouble(); + mpImpl->mVerbosity = Settings["verbosity"].GetInt(); + KRATOS_CATCH("") + + // This operation expects the pre-tension surface to be defined + // by a set of Geometries in the input ModelPart. It requires + // at least one Geometry. + KRATOS_ERROR_IF(mpImpl->mpModelPart->Geometries().empty()) + << "InsertPreTensionOperation requires at least one Geometry in the provided ModelPart, " + << "but there aren't any in '" << mpImpl->mpModelPart->Name() << "'."; +} + + +InsertPreTensionOperation::InsertPreTensionOperation(InsertPreTensionOperation&&) noexcept = default; + + +InsertPreTensionOperation::~InsertPreTensionOperation() = default; + + +struct PreTensionSurfacePartition { + array_1d normal; + std::vector positive_side, negative_side; +}; // struct SurfacePartition + + +void FindElementsAdjacentToGeometries(ModelPart& rModelPart) { + using NeighborElements = decltype(NEIGHBOUR_ELEMENTS)::Type; + + // Make sure every node on the pre-tension surface is in the model part. + KRATOS_TRY + ModelPart::NodesContainerType nodes; + for (Geometry& r_geometry : rModelPart.Geometries()) { + nodes.insert(r_geometry.begin(), r_geometry.end()); + } // for r_geometry in rModelPart.Conditions + rModelPart.Nodes().insert(nodes); + KRATOS_CATCH("") + + KRATOS_TRY + // Find elements containing each node in the model part. + FindGlobalNodalElementalNeighboursProcess(rModelPart.GetRootModelPart()).Execute(); + + block_for_each( + rModelPart.Geometries(), + [] (Geometry& r_geometry) -> void { + NeighborElements neighbor_elements; + + // Collect all neighbor elements from the geometry's nodes. + for (Node& r_node : r_geometry) { + NeighborElements& r_containing_elements = r_node.GetValue(NEIGHBOUR_ELEMENTS); + for (unsigned i_element=0u; i_element bool { + return rp_other_element->Id() == rp_element->Id(); + }); + if (it_element == neighbor_elements.ptr_end()) { + neighbor_elements.push_back(rp_element); + } + } + } + + // Filter elements that don't contain all nodes of the geometry. + neighbor_elements.erase(std::remove_if( + neighbor_elements.begin(), + neighbor_elements.end(), + [&r_geometry] (const Element& r_element) -> bool { + for (const Node& r_geometry_node : r_geometry) { + const auto it_node = std::find_if( + r_element.GetGeometry().begin(), + r_element.GetGeometry().end(), + [&r_geometry_node] (const Node& r_element_node) -> bool { + return r_element_node.Id() == r_geometry_node.Id(); + }); + if (it_node == r_element.GetGeometry().end()) + return true; + } + return false; + } + ), neighbor_elements.end()); + + // Set neighbor elements on the condition. + r_geometry.SetValue(NEIGHBOUR_ELEMENTS, neighbor_elements); + }); + KRATOS_CATCH("") +} + + +template +PreTensionSurfacePartition PartitionPreTensionSurface( + const ModelPart::GeometryContainerType::iterator itGeometryBegin, + const ModelPart::GeometryContainerType::iterator itGeometryEnd) { + PreTensionSurfacePartition surface_partition; + + // There are three valid possibilities here, separately for split and unsplit cases. + // 1) All surfaces are 0-dimensional (corners) + // => expecting exactly 1 geometry + // - unsplit case + // => expecting exactly 2 connected elements, both + // must have line geometries + // => the normal is not well-defined, so it is chosen + // as the average of the connected elements' directions + // - split case + // => expecting exactly 1 connected element that must have a line geometry + // => the normal is not well-defined, so it is chosen + // as the connected element's direction + // 2) All surfaces are 1-dimensional (edges) + // => expecting at least 1 geometry + // - unsplit case + // => expecting at least 2 connected elements, all of which + // must be 2-dimensional + // => the normal is not well-defined, so it is chosen as the + // average of the rotated endpoints in 2D, in positive + // direction. This assumes coordinates in z-direction vanish. + // This assumption is checked. + // - split case + // => expecting at least 1 connected element, all of which + // must be 2-dimensional + // => the normal is not well-defined, so it is chosen as the + // average of the rotated endpoints in 2D, in positive + // direction. This assumes coordinates in z-direcion vanish. + // This assumption is checked. + // 3) All surfaces are 2-dimensional (surfaces) + // => expecting at least 1 geometry + // => normals are well-defined + // - unsplit case + // => expecting at least 2 connected elements, all + // of which must be 3-dimensional + // - split case + // => expecting at least 1 connected element, all + // of which must be 3-dimensional + // + // In all unsplit cases, at least one connected element is required on + // either side of the average pre-tension plane. + + // Sanity checks. + if constexpr (SurfaceDimension == 0u) { + KRATOS_ERROR_IF_NOT(std::distance(itGeometryBegin, itGeometryEnd) == 1) + << "Expecting exactly 1 geometry on the 0-dimensional pre-tension surface , but got " + << std::distance(itGeometryBegin, itGeometryEnd) << "."; + } else if constexpr (SurfaceDimension == 1u) { + KRATOS_ERROR_IF_NOT(1 <= std::distance(itGeometryBegin, itGeometryEnd)) + << "Expecting at least 1 geometry on the 1-dimensional pre-tension surface , but got " + << std::distance(itGeometryBegin, itGeometryEnd) << "."; + } else if constexpr (SurfaceDimension == 2u) { + KRATOS_ERROR_IF_NOT(1 <= std::distance(itGeometryBegin, itGeometryEnd)) + << "Expecting at least 1 geometry on the 2-dimensional pre-tension surface , but got " + << std::distance(itGeometryBegin, itGeometryEnd) << "."; + } else { + static_assert(float(SurfaceDimension) == 0.5f, "Invalid surface dimension"); + } + + for (auto it_geometry=itGeometryBegin; it_geometry!=itGeometryEnd; ++it_geometry) { + Geometry& r_geometry = *it_geometry; + + // Sanity checks on the geometry. + KRATOS_ERROR_IF_NOT(r_geometry.LocalSpaceDimension() == SurfaceDimension) + << "Expecting all geometries defining the pre-tension surface to be " << SurfaceDimension << " dimensional, " + << "but geometry " << r_geometry.Id() << " on geometry " << r_geometry.Id() << " (" << r_geometry.Name() << ") " + << "is " << r_geometry.LocalSpaceDimension() << "-dimensional."; + + if constexpr (IsSplit) { + KRATOS_ERROR_IF_NOT(r_geometry.Has(NEIGHBOUR_ELEMENTS)) + << "Expecting geometry " << r_geometry.Id() << " to have exactly " + << "1 neighboring element, but found none."; + } else { + KRATOS_ERROR_IF_NOT(r_geometry.Has(NEIGHBOUR_ELEMENTS)) + << "Expecting geometry " << r_geometry.Id() << " to have exactly " + << "2 neighboring elements, but found none."; + } + + auto& r_neighbor_elements = r_geometry.GetValue(NEIGHBOUR_ELEMENTS); + if (r_neighbor_elements.size() != (IsSplit ? 1 : 2)) { + std::stringstream message; + message << "Expecting geometry " << r_geometry.Id() << " to have exactly " + << (IsSplit ? '1' : '2') + << " neighboring elements, but found " << r_neighbor_elements.size() << ": ["; + for (const Element& r_element : r_neighbor_elements) message << r_element.Id() << " "; + message << "]."; + KRATOS_ERROR << message.str(); + } + + // Check whether the 2 connected elements lie on opposite sides of the surface (only in the unsplit case). + const array_1d + first_connected_direction = r_neighbor_elements[0].GetGeometry().Center() - r_geometry.Center(), + second_connected_direction = r_neighbor_elements[IsSplit ? 0 : 1].GetGeometry().Center() - r_geometry.Center(); + const double direction_inner_product = std::inner_product( + first_connected_direction.begin(), + first_connected_direction.end(), + second_connected_direction.begin(), + 0.0); + KRATOS_ERROR_IF_NOT(IsSplit || direction_inner_product < 0.0) + << "Element " << r_neighbor_elements[0].Id() << " and " << r_neighbor_elements[1].Id() << " " + << "do not lie on opposite sides of condition " << r_geometry.Id() << " " + << "defining a pre-tension surface."; + + if constexpr (SurfaceDimension == 0u) { + for (unsigned i_element=0u; i_elementoperator[](0).Id(); + const array_1d surface_node_position = itGeometryBegin->operator[](0); + array_1d opposite_node_position; + if (r_element.GetGeometry()[0].Id() == surface_node_id) { + opposite_node_position = r_element.GetGeometry()[1]; + } else if (r_element.GetGeometry()[1].Id() == surface_node_id) { + opposite_node_position = r_element.GetGeometry()[0]; + } else { + KRATOS_ERROR + << "Element " << r_element.Id() << " is assumed to be connected to node " + << surface_node_id << " but is not."; + } + + if (IsSplit || i_element % 2) { + surface_partition.positive_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); + surface_partition.normal = opposite_node_position - surface_node_position; + } else { + surface_partition.negative_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); + } + } // for i_element in range(r_neighbor_elements.size()) + } /*if SurfaceDimension == 0*/ else if constexpr (SurfaceDimension == 1u) { + // Define the line geometry's plane. + const array_1d + surface_begin = r_geometry[0], + surface_end = *(r_geometry.end() - 1); + KRATOS_ERROR_IF_NOT(Comparison::Equal(1e-12, 1e-9)(surface_begin[2], 0.0)) + << "Assuming that 1-dimensional pre-tension surfaces lie on the XY plane, but node " + << r_geometry[0].Id() << " of geometry " << r_geometry.Id() << " " + << "has a Z-coordinate of " << surface_begin[2] << "."; + KRATOS_ERROR_IF_NOT(Comparison::Equal(1e-12, 1e-9)(surface_end[2], 0.0)) + << "Assuming that 1-dimensional pre-tension surfaces lie on the XY plane, but node " + << r_geometry[1].Id() << " of geometry " << r_geometry.Id() << " " + << "has a Z-coordinate of " << surface_end[2] << "."; + const array_1d + surface_normal { + -(surface_end[1] - surface_begin[1]), + surface_end[0] - surface_begin[0], + 0.0}, + surface_center = r_geometry.Center(); + surface_partition.normal += surface_normal; + + // Sort connected elements. + for (unsigned i_element=0u; i_element element_direction = r_element.GetGeometry().Center() - surface_center; + const double direction_inner_product = std::inner_product( + surface_normal.begin(), + surface_normal.end(), + element_direction.begin(), + 0.0); + if (0 < direction_inner_product) surface_partition.negative_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); + else surface_partition.positive_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); + } // for i_element in range(r_neighbor_elements.size()) + + KRATOS_ERROR_IF_NOT(IsSplit || surface_partition.positive_side.size() == surface_partition.negative_side.size()) + << "Elements " << r_neighbor_elements[0].Id() << " and " << r_neighbor_elements[1].Id() << " " + << "lie on the same side of geometry " << r_geometry.Id() << " that defines a pre-tension surface."; + } /*else if SurfaceDimension == 1*/ else if constexpr (SurfaceDimension == 2u) { + // Define the surface's plane. + const array_1d surface_center = r_geometry.Center(); + const array_1d surface_normal = r_geometry.Normal(surface_center); + surface_partition.normal += surface_normal; + + // Sort connected elements. + for (unsigned i_element=0u; i_element element_direction = r_element.GetGeometry().Center() - surface_center; + const double direction_inner_product = std::inner_product( + surface_normal.begin(), + surface_normal.end(), + element_direction.begin(), + 0.0); + if (0 < direction_inner_product) surface_partition.negative_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); + else surface_partition.positive_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); + } // for i_element in range(r_neighbor_elements.size()) + + KRATOS_ERROR_IF_NOT(IsSplit || surface_partition.positive_side.size() == surface_partition.negative_side.size()) + << "Elements " << r_neighbor_elements[0].Id() << " and " << r_neighbor_elements[1].Id() << " " + << "lie on the same side of condition " << r_geometry.Id() << " that defines a pre-tension surface."; + } /*else if SurfaceDimension == 2*/ + } // for r_condition in r_model_part.Conditions + + // Normalize the surface's normal. + const double norm = std::sqrt(std::accumulate( + surface_partition.normal.begin(), + surface_partition.normal.end(), + 0.0, + [] (double sum, double component) -> double { + return sum + component * component; + })); + KRATOS_ERROR_IF_NOT(norm) << "normal vanished on the pretension surface"; + std::transform( + surface_partition.normal.begin(), + surface_partition.normal.end(), + surface_partition.normal.begin(), + [norm] (double component) -> double { + return component / norm; + }); + + return surface_partition; +} + + +/// @brief Base class for @ref InPlaneRelativeDisplacementConstraint and @ref OutOfPlaneDisplacementConstraint. +class PlaneDisplacementConstraint : public MultifreedomConstraint { +public: + using ProtectedNormal = std::pair< + std::optional>, + LockObject>; + + PlaneDisplacementConstraint() noexcept = default; + + PlaneDisplacementConstraint( + MasterSlaveConstraint::IndexType Id, + const ModelPart::GeometryContainerType& rSurfaceGeometries, + DofPointerVectorType&& rDofs, + const std::vector& rConstraintLabels, + std::shared_ptr pSharedSurfaceNormal, + int Verbosity) + : MultifreedomConstraint( + Id, + std::move(rDofs), + rConstraintLabels), + mVerbosity(Verbosity), + mSurfaceGeometries(rSurfaceGeometries), + mpSharedSurfaceNormal(pSharedSurfaceNormal) { + KRATOS_ERROR_IF(rSurfaceGeometries.empty()); + } + + void InitializeSolutionStep(const ProcessInfo& rProcessInfo) override { + this->InitializeNonLinearIteration(rProcessInfo); + } + + /// @brief Linearize the constraint equation. + void InitializeNonLinearIteration(const ProcessInfo&) override { + KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFileName() << " is virtual"; + } + + void FinalizeNonLinearIteration(const ProcessInfo&) override { + KRATOS_TRY + std::scoped_lock lock(mpSharedSurfaceNormal->second); + mpSharedSurfaceNormal->first.reset(); + KRATOS_CATCH("") + } + + void FinalizeSolutionStep(const ProcessInfo& rProcessInfo) override { + this->FinalizeNonLinearIteration(rProcessInfo); + } + + void CalculateLocalSystem( + MatrixType& rConstraintGradient, + VectorType& rConstraintGaps, + const ProcessInfo&) const override { + rConstraintGradient = this->GetData().GetValue(CONSTITUTIVE_MATRIX); + rConstraintGaps = this->GetData().GetValue(INTERNAL_FORCES_VECTOR); + } + + std::string Info() const override { + return "PlaneDisplacementConstraint"; + } + + void save(Serializer& rSerializer) const override { + KRATOS_TRY + KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, MultifreedomConstraint); + rSerializer.save("SurfaceGeometries", mSurfaceGeometries); + rSerializer.save("Verbosity", mVerbosity); + KRATOS_CATCH("") + } + + void load(Serializer& rDeserializer) override { + KRATOS_TRY + KRATOS_SERIALIZE_LOAD_BASE_CLASS(rDeserializer, MultifreedomConstraint); + rDeserializer.load("SurfaceGeometries", mSurfaceGeometries); + rDeserializer.load("Verbosity", mVerbosity); + KRATOS_CATCH("") + } + +protected: + array_1d GetSurfaceNormal() { + KRATOS_TRY + std::scoped_lock lock(mpSharedSurfaceNormal->second); + if (mpSharedSurfaceNormal->first.has_value()) { + return mpSharedSurfaceNormal->first.value(); + } else { + PreTensionSurfacePartition surface_partition; + switch (mSurfaceGeometries.front().LocalSpaceDimension()) { + case 0: { + surface_partition = PartitionPreTensionSurface<0,true>( + mSurfaceGeometries.begin(), + mSurfaceGeometries.end()); + break; + } + case 1: { + surface_partition = PartitionPreTensionSurface<1,true>( + mSurfaceGeometries.begin(), + mSurfaceGeometries.end()); + break; + } + case 2: { + surface_partition = PartitionPreTensionSurface<2,true>( + mSurfaceGeometries.begin(), + mSurfaceGeometries.end()); + break; + } + default: KRATOS_ERROR << "Invalid pre-tension surface dimension: " << mSurfaceGeometries.front().LocalSpaceDimension() << "."; + } // switch dimension + + KRATOS_INFO_IF(this->Info() + " " + std::to_string(this->Id()), 3 <= mVerbosity) + << "pre-tension surface normal: [" + << surface_partition.normal[0] << ' ' + << surface_partition.normal[1] << ' ' + << surface_partition.normal[2] << "]\n"; + mpSharedSurfaceNormal->first = surface_partition.normal; + return surface_partition.normal; + } + KRATOS_CATCH("") + } + + int mVerbosity; + + ModelPart::GeometryContainerType mSurfaceGeometries; + + /// @details The surface's normal is potentially shared by a set of constraints, + /// since each constraint only restricts one pair of nodes. If the surface + /// contains more than a single pair, computing its normal should not be + /// unnecessarily repeated. + std::shared_ptr>, + LockObject + >> mpSharedSurfaceNormal; +}; // class PlaneDisplacementConstraint + + +/// @brief Forbid the relative in-plane displacement of a pair of nodes. +/// @details The surface's normal is computed in each non-linear iteration +/// as the weighted average of a set of geometries' normals. +class InPlaneRelativeDisplacementConstraint : public PlaneDisplacementConstraint { +public: + using PlaneDisplacementConstraint::PlaneDisplacementConstraint; + + MasterSlaveConstraint::Pointer Clone(IndexType Id) const override { + const auto& r_constraint_labels = this->GetValue(CONSTRAINT_LABELS); + std::vector constraint_labels( + r_constraint_labels.begin(), + r_constraint_labels.end()); + return MasterSlaveConstraint::Pointer(new InPlaneRelativeDisplacementConstraint( + Id, + mSurfaceGeometries, + DofPointerVectorType(this->GetDofs()), + constraint_labels, + mpSharedSurfaceNormal, + mVerbosity)); + } + + void InitializeNonLinearIteration(const ProcessInfo&) override { + KRATOS_TRY + const array_1d normal = this->GetSurfaceNormal(); + + // DoFs are assumed to be stored in a specific layout. + // - contained DoFs are related to two (an original and its duplicate) nodes + // - the same set (n) of variables are referenced + // - the first half of the array contains DISPLACEMENT_X, DISPLACEMENT_Y ... of the original node in that order + // - the second half is identical, but for the duplicate node + auto& r_dofs = this->GetDofs(); + KRATOS_ERROR_IF(r_dofs.size() % 2); + const std::size_t dimension_count = r_dofs.size() / 2; + + Matrix constraint_gradients(dimension_count, 2 * dimension_count); + Vector constraint_gaps, displacements(2 * dimension_count); + + // Build the constraint equations. + for (std::size_t i_dimension=0ul; i_dimensionGetSolutionStepValue(); + displacements[i_dimension + dimension_count] = r_dofs[i_dimension + dimension_count]->GetSolutionStepValue(); + } // for i_dimension in range(dimension_count) + + constraint_gaps = prod(constraint_gradients, displacements); + + // Store constraint gradients and constraint gaps. + this->SetValue(CONSTITUTIVE_MATRIX, constraint_gradients); + this->SetValue(INTERNAL_FORCES_VECTOR, constraint_gaps); + KRATOS_CATCH("") + } + + std::string Info() const override { + return "InPlaneRelativeDisplacementConstraint"; + } +}; // class InPlaneRelativeDisplacementConstraint + + +Dof* FindDofInNode(Node& rNode, const Variable& rVariable) noexcept { + auto& r_dofs = rNode.GetDofs(); + const auto it_dof = std::find_if( + r_dofs.begin(), + r_dofs.end(), + [&rVariable] (const auto& rp_dof) -> bool { + return rp_dof->GetVariable() == rVariable; + }); + return it_dof == r_dofs.end() + ? nullptr + : it_dof->get(); +} + + +void InsertPreTensionOperation::Execute() { + ModelPart& r_model_part = *mpImpl->mpModelPart; + + // Find elements connected to the pre-tension surface. + // => each node will have a NEIGHBOR_ELEMENTS variable + // in its non-historical storage that lists elements + // containing them. + KRATOS_TRY + FindElementsAdjacentToGeometries(r_model_part); + KRATOS_CATCH("") + + // Sort connected elements into 2 categories, lying on either side of + // the pre-tension surface. + PreTensionSurfacePartition surface_partition; + KRATOS_TRY + // Find the dimension of the pre-tension surface, + // and make sure all conditions share it. + KRATOS_ERROR_IF(r_model_part.Geometries().empty()); + const int pre_tension_surface_dimension = r_model_part.Geometries().front().LocalSpaceDimension(); + switch (pre_tension_surface_dimension) { + case 0u: { + surface_partition = PartitionPreTensionSurface<0>( + r_model_part.Geometries().begin(), + r_model_part.Geometries().end()); + break; + } + case 1u: { + surface_partition = PartitionPreTensionSurface<1>( + r_model_part.Geometries().begin(), + r_model_part.Geometries().end()); + break; + } + case 2u: { + surface_partition = PartitionPreTensionSurface<2>( + r_model_part.Geometries().begin(), + r_model_part.Geometries().end()); + break; + } + default: KRATOS_ERROR << "Invalid pre-tension surface dimension: " << pre_tension_surface_dimension << "."; + } + KRATOS_CATCH("") + + if (3 <= mpImpl->mVerbosity) { + KRATOS_INFO(this->Info()) + << "pre-tension surface normal: [" + << surface_partition.normal[0] << ' ' + << surface_partition.normal[1] << ' ' + << surface_partition.normal[2] << "]\n"; + + std::stringstream message; + if (3 <= mpImpl->mVerbosity) { + message << "elements on the negative side: ["; + for (const auto& rp_element : surface_partition.negative_side) + message << rp_element->Id() << ' '; + message << "]\n"; + KRATOS_INFO(this->Info()) << message.view(); + + message = std::stringstream(); + message << "elements on the positive side: ["; + for (const auto& rp_element : surface_partition.positive_side) + message << rp_element->Id() << ' '; + message << "]\n"; + KRATOS_INFO(this->Info()) << message.view(); + } // if 3 <= verbosity + } // if 3 <= verbosity + + // Insert a new control node acting as a master for all + // other nodes on the pre-tension surface. + Node::Pointer p_control_node = mpImpl->mpModelPart->CreateNewNode( + /*Id=*/mpImpl->mpModelPart->GetRootModelPart().Nodes().empty() + ? 1 + : mpImpl->mpModelPart->GetRootModelPart().Nodes().back().Id() + 1, + /*x=*/0.0, + /*y=*/0.0, + /*z=*/0.0); + p_control_node->AddDof(DISPLACEMENT_X); + KRATOS_INFO_IF(this->Info(), 2 <= mpImpl->mVerbosity) + << "insert control node " << p_control_node->Id() << ' ' + << "at [" << p_control_node->X() << ' ' << p_control_node->Y() << ' ' << p_control_node->Z() << "]\n"; + + // Duplicate nodes on the pre-tension surface. + std::unordered_map duplicated_nodes; + KRATOS_TRY + std::unordered_set surface_nodes; + + // Collect nodes on the surface. + for (Geometry& r_geometry : r_model_part.Geometries()) { + for (Node& r_node : r_geometry) { + surface_nodes.insert(&r_node); + } // for r_node in r_geometry + } // for r_geometry in r_model_part.Conditions() + + // Find the first available node ID that can safely be + // used to insert new ones. + Node::IndexType node_id = r_model_part.GetRootModelPart().Nodes().back().Id() + 1; + + // Clone surface nodes and associate them with the original ones. + for (Node* p_node : surface_nodes) { + const auto emplace_result = duplicated_nodes.emplace( + p_node, + p_node->Clone(node_id++)); + KRATOS_INFO_IF(this->Info(), 3 <= mpImpl->mVerbosity) + << "duplicate node " << p_node->Id() << " => " << emplace_result.first->second->Id() << "\n"; + r_model_part.AddNode(emplace_result.first->second); + } // for p_node in surface_nodes + + // Replace nodes on the negative side elements with duplicated ones. + for (auto& rp_element : surface_partition.negative_side) { + for (auto itp_node=rp_element->GetGeometry().ptr_begin(); itp_node!=rp_element->GetGeometry().ptr_end(); ++itp_node) { + const auto it_duplicate_pair = duplicated_nodes.find(&**itp_node); + if (it_duplicate_pair != duplicated_nodes.end()) { + (*itp_node) = it_duplicate_pair->second; + } + } // for itp_node in rp_element->GetGeometry() + } // for rp_element in surface_partition.negative_side + KRATOS_CATCH("") + + // Reset NEIGHBOR_ELEMENTS after duplicating nodes. + KRATOS_TRY + FindElementsAdjacentToGeometries(r_model_part); + KRATOS_CATCH("") + + // Collect DoFs from elements on the positive side of the pre-tension surface. + // This is necessary because elements on the pre-tension surface may not have + // DoFs for all displacement components. In this case, constraints should not + // be inserted for these components. + std::unordered_set*> positive_side_dofs; + { + Element::DofsVectorType dof_buffer; + for (const auto& rp_element : surface_partition.positive_side) { + rp_element->GetDofList(dof_buffer, mpImpl->mpModelPart->GetProcessInfo()); + positive_side_dofs.insert( + dof_buffer.begin(), + dof_buffer.end()); + } + } + + LinearMultifreedomConstraint::IndexType id_constraint = mpImpl->mpModelPart->GetRootModelPart().MasterSlaveConstraints().empty() + ? 1 + : mpImpl->mpModelPart->GetRootModelPart().MasterSlaveConstraints().back().Id() + 1; + + // Forbid relative in-plane displacements for each duplicated node pair. + { + const std::array*,3> all_displacement_components { + &DISPLACEMENT_X, + &DISPLACEMENT_Y, + &DISPLACEMENT_Z}; + + std::vector*,2>> dof_pairs; + auto p_protected_surface_normal = std::make_shared>, + LockObject>>(); + + for (auto& [rp_positive_side_node, rp_negative_side_node] : duplicated_nodes) { + dof_pairs.clear(); + + // Collect relevant DoF pairs. + // dof_pairs will contain [{u_x_positive, u_x_negative}, {u_y_positive, u_y_negative}, ...] + for (const auto& rp_variable : all_displacement_components) { + Dof* p_positive_side_dof = FindDofInNode(*rp_positive_side_node, *rp_variable); + if (positive_side_dofs.find(p_positive_side_dof) == positive_side_dofs.end()) continue; + + Dof* p_negative_side_dof = FindDofInNode(*rp_negative_side_node, *rp_variable); + if (p_positive_side_dof && p_negative_side_dof) + dof_pairs.push_back({p_positive_side_dof, p_negative_side_dof}); + } // for rp_variable in all_displacement_components + + MasterSlaveConstraint::DofPointerVectorType dofs(2 * dof_pairs.size()); + for (std::size_t i_pair=0ul; i_pair constraint_labels(dof_pairs.size()); + std::iota( + constraint_labels.begin(), + constraint_labels.end(), + id_constraint); + MasterSlaveConstraint::Pointer p_constraint(new InPlaneRelativeDisplacementConstraint( + id_constraint, + mpImpl->mpModelPart->Geometries(), + std::move(dofs), + constraint_labels, + p_protected_surface_normal, + mpImpl->mVerbosity)); + + KRATOS_INFO_IF(this->Info(), 3 <= mpImpl->mVerbosity) + << "insert constraint " << p_constraint->Id() << ' ' + << "prohibiting the in-plane relative displacement of nodes " + << rp_positive_side_node->Id() << " and " << rp_negative_side_node->Id() << '\n'; + mpImpl->mpModelPart->AddMasterSlaveConstraint(p_constraint); + + id_constraint += constraint_labels.size(); + } // for id_positive_side_node, rp_negative_side_node : duplicated_nodes + } + + // Tie rotation DoFs between duplicate nodes. + { + const std::array*,3> all_rotation_components { + &ROTATION_X, + &ROTATION_Y, + &ROTATION_Z}; + + for (auto& [rp_positive_side_node, rp_negative_side_node] : duplicated_nodes) { + for (const auto& rp_variable : all_rotation_components) { + Dof* p_positive_side_dof = FindDofInNode(*rp_positive_side_node, *rp_variable); + if (positive_side_dofs.find(p_positive_side_dof) == positive_side_dofs.end()) continue; + Dof* p_negative_side_dof = FindDofInNode(*rp_negative_side_node, *rp_variable); + if (p_positive_side_dof && p_negative_side_dof) { + LinearMultifreedomConstraint::DofPointerVectorType dofs(2); + LinearMultifreedomConstraint::MatrixType constraint_gradient(1, 2); + LinearMultifreedomConstraint::VectorType constraint_gap(1); + dofs[0] = p_positive_side_dof; + dofs[1] = p_negative_side_dof; + constraint_gradient(0, 0) = 1.0; + constraint_gradient(0, 1) = -1.0; + constraint_gap[0] = 0.0; + MasterSlaveConstraint::Pointer p_constraint(new LinearMultifreedomConstraint( + id_constraint, + std::move(dofs), + {id_constraint}, + constraint_gradient, + constraint_gap)); + KRATOS_INFO_IF(this->Info(), 3 <= mpImpl->mVerbosity) + << "tie " << p_positive_side_dof->GetVariable().Name() << ' ' + << "of node " << rp_positive_side_node->Id() << ' ' + << "to " << p_negative_side_dof->GetVariable().Name() << ' ' + << "of node " << rp_negative_side_node->Id() << '\n'; + mpImpl->mpModelPart->AddMasterSlaveConstraint(p_constraint); + ++id_constraint; + } // if p_positive_side_dof && p_negative_side_dof + } // for rp_variable in all_rotation_components + } // for rp_positive_side_node, rp_negative_side_node + } + + // Decide what to do with the control node in derived classes. + this->InsertControlNodeConstraints( + *mpImpl->mpModelPart, + surface_partition.normal, + duplicated_nodes, + p_control_node, + positive_side_dofs); +} + + +const Parameters InsertPreTensionOperation::GetDefaultParameters() const { + return Parameters(R"({ + "model_part_name" : "", + "magnitude" : 0.0, + "verbosity" : 1 + })"); +} + + +std::string InsertPreTensionOperation::Info() const { + return "InsertPreTensionOperation"; +} + + +/// @brief Tie the average [relative] out-of-plane displacement to a common control DoF. +template +class OutOfPlaneDisplacementConstraint : public PlaneDisplacementConstraint { +public: + using PlaneDisplacementConstraint::PlaneDisplacementConstraint; + + MasterSlaveConstraint::Pointer Clone(IndexType Id) const override { + const auto& r_constraint_labels = this->GetValue(CONSTRAINT_LABELS); + std::vector constraint_labels( + r_constraint_labels.begin(), + r_constraint_labels.end()); + return MasterSlaveConstraint::Pointer(new OutOfPlaneDisplacementConstraint( + Id, + mSurfaceGeometries, + DofPointerVectorType(this->GetDofs()), + constraint_labels, + mpSharedSurfaceNormal, + mVerbosity)); + } + + void InitializeNonLinearIteration(const ProcessInfo&) override { + KRATOS_TRY + // Compute the plane's normal. + const array_1d normal = this->GetSurfaceNormal(); + + // DoFs are assumed in a specific layout. + // - DISPLACEMENT_X of node 0 + // - DISPLACEMENT_Y of node 0 + // ... + // - DISPLACEMENT_Z (depending on the plane's dimension) of node n + // - DISPLACEMENT_X of node 0's pair (only if relative) + // - DISPLACEMENT_Y of node 0's pair (only if relative) + // ... + // - DISPLACEMENT_Z (depending on the plane's dimension) of node n's pair (only if relative) + // - DoF of the control node + auto& r_dofs = this->GetDofs(); + KRATOS_ERROR_IF(r_dofs.empty()); + + const std::size_t controlled_dof_count = r_dofs.size() - 1; + const std::size_t positive_side_dof_count = IsRelative + ? controlled_dof_count / 2 + : controlled_dof_count; + + // Find the plane's dimension by checking how many unique variables the DoFs reference. + std::size_t dimension_count = 0ul; + { + std::unordered_set keys; + for (std::size_t i_dof=0ul; i_dofGetVariable().Key()); + } // for i_dof in positive_side_dof_count + dimension_count = keys.size(); + KRATOS_ERROR_IF_NOT(0 < dimension_count && dimension_count < 4) + << "unexpected surface dimension " << dimension_count; + } + + // Build the linearized constraint equation. + Matrix constraint_gradient(1, r_dofs.size()); + Vector constraint_gap(1), displacements(r_dofs.size()); + + for (std::size_t i_dof=0ul; i_dof* p_dof) -> double { + return p_dof->GetSolutionStepValue(); + }); + + constraint_gap += prod(constraint_gradient, displacements); + + // Register the linearized constraint equation. + this->SetValue(CONSTITUTIVE_MATRIX, constraint_gradient); + this->SetValue(INTERNAL_FORCES_VECTOR, constraint_gap); + KRATOS_CATCH("") + } + + std::string Info() const override { + if constexpr (IsRelative) + return "OutOfPlaneRelativeDisplacementConstraint"; + else + return "OutOfPlaneDisplacementConstraint"; + } +}; // class OutOfPlaneDisplacementConstraint + + +void InsertDirichletPreTensionOperation::InsertControlNodeConstraints( + ModelPart& rModelPart, + array_1d SurfaceNormal, + const std::unordered_map rDuplicateNodeMap, + Node::Pointer pControlNode, + const std::unordered_set*> rPositiveSideDofs) const { + KRATOS_TRY + // Insert a constraint that ties the average out-of-plane relative displacement + // to a prescribed value. + const std::array*,3> all_displacement_components { + &DISPLACEMENT_X, + &DISPLACEMENT_Y, + &DISPLACEMENT_Z}; + LinearMultifreedomConstraint::DofPointerVectorType dofs; + + for (const auto& [rp_positive_side_node, rp_negative_side_node] : rDuplicateNodeMap) { + for (std::size_t i_variable=0ul; i_variable& r_variable = *all_displacement_components[i_variable]; + Dof* p_positive_side_dof = FindDofInNode(*rp_positive_side_node, r_variable); + if (rPositiveSideDofs.find(p_positive_side_dof) == rPositiveSideDofs.end()) continue; + + Dof* p_negative_side_dof = FindDofInNode(*rp_negative_side_node, r_variable); + if (p_positive_side_dof && p_negative_side_dof) { + dofs.push_back(p_positive_side_dof); + dofs.push_back(p_negative_side_dof); + } // if p_positive_side_dof && p_negative_side_dof + } // for rp_variable in all_displacement_components + } // for rp_positive_side_node, rp_negative_side_node + + // Add the control node's DoF to the constraint equation. + Dof* p_control_dof = pControlNode->GetDofs()[0].get(); + dofs.push_back(p_control_dof); + + // Construct and register the constraint equation. + const LinearMultifreedomConstraint::IndexType constraint_id = mpImpl->mpModelPart->GetRootModelPart().MasterSlaveConstraints().empty() + ? 1 + : mpImpl->mpModelPart->GetRootModelPart().MasterSlaveConstraints().back().Id() + 1; + auto p_protected_surface_normal = std::make_shared>, + LockObject>>(); + MasterSlaveConstraint::Pointer p_constraint(new OutOfPlaneDisplacementConstraint( + constraint_id, + mpImpl->mpModelPart->Geometries(), + std::move(dofs), + {constraint_id}, + p_protected_surface_normal, + mpImpl->mVerbosity)); + mpImpl->mpModelPart->AddMasterSlaveConstraint(p_constraint); + + // Set a dirichlet condition on the control node's DoF. + p_control_dof->GetSolutionStepValue() = (mpImpl->mPretensionValue - p_control_dof->GetSolutionStepValue()) * rDuplicateNodeMap.size(); + p_control_dof->FixDof(); + KRATOS_CATCH("") +} + + +std::string InsertDirichletPreTensionOperation::Info() const { + return "InsertDirichletPreTensionOperation"; +} + + +class PointLoadCondition1D1N : public Condition { +public: + KRATOS_CLASS_POINTER_DEFINITION(PointLoadCondition1D1N); + + PointLoadCondition1D1N( + Condition::IndexType Id, + Geometry::Pointer pGeometry) + : Condition(Id, pGeometry) + {} + + Condition::Pointer Clone( + Condition::IndexType Id, + const Condition::NodesArrayType& rNodes) const override { + KRATOS_ERROR_IF_NOT(rNodes.size() == 1) + << "PointLoadCondition1D1N::Clone expects 1 node but got " << rNodes.size(); + return Condition::Pointer(new PointLoadCondition1D1N( + Id, + Geometry::Pointer(new Point2D(rNodes)))); + } + + void EquationIdVector( + Condition::EquationIdVectorType& rIndices, + const ProcessInfo&) const override { + rIndices.resize(1); + rIndices[0] = this->GetGeometry()[0].GetDofs()[0]->EquationId(); + } + + void GetDofList( + Condition::DofsVectorType& rDofs, + const ProcessInfo&) const override { + rDofs.resize(1); + rDofs[0] = this->GetGeometry()[0].GetDofs()[0].get(); + } + + void GetValuesVector( + Vector& rOutput, + int Step) const override { + rOutput.resize(1); + rOutput[0] = this->GetGeometry()[0].FastGetSolutionStepValue(DISPLACEMENT_X, Step); + } + + void GetFirstDerivativesVector( + Vector& rOutput, + int Step) const override { + rOutput.resize(1); + rOutput[0] = this->GetGeometry()[0].FastGetSolutionStepValue(VELOCITY_X, Step); + } + + void GetSecondDerivativesVector( + Vector& rOutput, + int Step) const override { + rOutput.resize(1); + rOutput[0] = this->GetGeometry()[0].FastGetSolutionStepValue(ACCELERATION_X, Step); + } + + void CalculateRightHandSide( + Condition::VectorType& rRhs, + const ProcessInfo&) override { + rRhs.resize(1); + rRhs[0] = this->GetValue(POINT_LOAD_X); + } + + void CalculateLocalSystem( + Condition::MatrixType& rLhs, + Condition::VectorType& rRhs, + const ProcessInfo& rProcessInfo) override { + rLhs.resize(1, 1, false); + rLhs(0, 0) = 0.0; + this->CalculateRightHandSide(rRhs, rProcessInfo); + } + + void CalculateMassMatrix( + Condition::MatrixType& rMatrix, + const ProcessInfo&) override { + rMatrix.resize(1, 1, false); + rMatrix(0, 0) = 0.0; + } + + void CalculateDampingMatrix( + Condition::MatrixType& rMatrix, + const ProcessInfo&) override { + rMatrix.resize(1, 1, false); + rMatrix(0, 0) = 0.0; + } +}; // class PointLoadCondition1D1N + + +void InsertNeumannPreTensionOperation::InsertControlNodeConstraints( + ModelPart& rModelPart, + array_1d SurfaceNormal, + const std::unordered_map rDuplicateNodeMap, + Node::Pointer pPositiveSideControlNode, + const std::unordered_set*> rPositiveSideDofs) const { + KRATOS_TRY + const std::array*,3> all_displacement_components { + &DISPLACEMENT_X, + &DISPLACEMENT_Y, + &DISPLACEMENT_Z}; + LinearMultifreedomConstraint::DofPointerVectorType positive_side_dofs, negative_side_dofs; + + for (const auto& [rp_positive_side_node, rp_negative_side_node] : rDuplicateNodeMap) { + for (std::size_t i_variable=0ul; i_variable& r_variable = *all_displacement_components[i_variable]; + Dof* p_positive_side_dof = FindDofInNode(*rp_positive_side_node, r_variable); + if (rPositiveSideDofs.find(p_positive_side_dof) == rPositiveSideDofs.end()) continue; + + Dof* p_negative_side_dof = FindDofInNode(*rp_negative_side_node, r_variable); + if (p_positive_side_dof && p_negative_side_dof) { + positive_side_dofs.push_back(p_positive_side_dof); + negative_side_dofs.push_back(p_negative_side_dof); + } // if p_positive_side_dof && p_negative_side_dof + } // for rp_variable in all_displacement_components + } // for rp_positive_side_node, rp_negative_side_node + + // Insert a new node acting as the control node on the negative side. + Node::Pointer p_negative_side_control_node = pPositiveSideControlNode->Clone( + mpImpl->mpModelPart->Nodes().empty() + ? 1 + : mpImpl->mpModelPart->Nodes().back().Id() + 1); + KRATOS_INFO_IF(this->Info(), 2 <= mpImpl->mVerbosity) + << "insert negative side control node " << p_negative_side_control_node->Id() << ' ' + << "at [" << p_negative_side_control_node->X() << ' ' + << p_negative_side_control_node->Y() << ' ' + << p_negative_side_control_node->Z() << "]\n"; + mpImpl->mpModelPart->AddNode(p_negative_side_control_node); + + // Add control DoFs. + positive_side_dofs.push_back(pPositiveSideControlNode->GetDofs()[0].get()); + negative_side_dofs.push_back(p_negative_side_control_node->GetDofs()[0].get()); + + LinearMultifreedomConstraint::IndexType constraint_id = mpImpl->mpModelPart->GetRootModelPart().MasterSlaveConstraints().empty() + ? 1 + : mpImpl->mpModelPart->GetRootModelPart().MasterSlaveConstraints().back().Id() + 1; + Condition::IndexType condition_id = mpImpl->mpModelPart->GetRootModelPart().Conditions().empty() + ? 1 + : mpImpl->mpModelPart->GetRootModelPart().Conditions().back().Id() + 1; + auto p_protected_surface_normal = std::make_shared>, + LockObject>>(); + + // Insert the positive side constraint and condition. + { + MasterSlaveConstraint::Pointer p_constraint(new OutOfPlaneDisplacementConstraint( + constraint_id, + mpImpl->mpModelPart->Geometries(), + std::move(positive_side_dofs), + {constraint_id}, + p_protected_surface_normal, + mpImpl->mVerbosity)); + Condition::Pointer p_condition(new PointLoadCondition1D1N( + condition_id, + Geometry::Pointer(new Point2D(pPositiveSideControlNode)))); + p_condition->SetValue(POINT_LOAD_X, mpImpl->mPretensionValue); + + KRATOS_INFO_IF(this->Info(), 3 <= mpImpl->mVerbosity) + << "insert constraint " << p_constraint->Id() <<' ' + << "tying the average out-of-plane displacement on the positive side to " + << pPositiveSideControlNode->GetDofs()[0]->GetVariable().Name() << ' ' + << "of node " << pPositiveSideControlNode->Id() << "\n"; + mpImpl->mpModelPart->AddMasterSlaveConstraint(p_constraint); + + KRATOS_INFO_IF(this->Info(), 3 <= mpImpl->mVerbosity) + << "insert condition " << p_condition->Id() <<' ' + << "loading " << pPositiveSideControlNode->GetDofs()[0]->GetVariable().Name() << ' ' + << "of node " << pPositiveSideControlNode->Id() << "\n"; + mpImpl->mpModelPart->AddCondition(p_condition); + + ++constraint_id; + ++condition_id; + } + + // Insert the negative side constraint and condition. + { + MasterSlaveConstraint::Pointer p_constraint(new OutOfPlaneDisplacementConstraint( + constraint_id, + mpImpl->mpModelPart->Geometries(), + std::move(negative_side_dofs), + {constraint_id}, + p_protected_surface_normal, + mpImpl->mVerbosity)); + Condition::Pointer p_condition(new PointLoadCondition1D1N( + condition_id, + Geometry::Pointer(new Point2D(p_negative_side_control_node)))); + p_condition->SetValue(POINT_LOAD_X, mpImpl->mPretensionValue); + + KRATOS_INFO_IF(this->Info(), 3 <= mpImpl->mVerbosity) + << "insert constraint " << p_constraint->Id() <<' ' + << "tying the average out-of-plane displacement on the positive side to " + << p_negative_side_control_node->GetDofs()[0]->GetVariable().Name() << ' ' + << "of node " << p_negative_side_control_node->Id() << "\n"; + mpImpl->mpModelPart->AddMasterSlaveConstraint(p_constraint); + + KRATOS_INFO_IF(this->Info(), 3 <= mpImpl->mVerbosity) + << "insert condition " << p_condition->Id() <<' ' + << "loading " << p_negative_side_control_node->GetDofs()[0]->GetVariable().Name() << ' ' + << "of node " << p_negative_side_control_node->Id() << "\n"; + mpImpl->mpModelPart->AddCondition(p_condition); + + ++constraint_id; + ++condition_id; + } + KRATOS_CATCH("") +} + + +std::string InsertNeumannPreTensionOperation::Info() const { + return "InsertNeumannPreTensionOperation"; +} + + +} // namespace Kratos diff --git a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp new file mode 100644 index 000000000000..ac4ecfa6f5b2 --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp @@ -0,0 +1,101 @@ +// KRATOS ___| | | | +// \___ \ __| __| | | __| __| | | __| _` | | +// | | | | | ( | | | | ( | | +// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS +// +// License: BSD License +// license: StructuralMechanicsApplication/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// --- Core Includes --- +#include "operations/operation.h" // Operation +#include "containers/model.h" // Model + +// --- STL Includes --- +#include +#include + + +namespace Kratos { + + +class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertPreTensionOperation : public Operation { +public: + KRATOS_CLASS_POINTER_DEFINITION(InsertPreTensionOperation); + + InsertPreTensionOperation(Model& rModel, Parameters Settings); + + InsertPreTensionOperation(InsertPreTensionOperation&&) noexcept; + + ~InsertPreTensionOperation(); + + void Execute() override; + + const Parameters GetDefaultParameters() const override; + + virtual std::string Info() const; + +protected: + virtual void InsertControlNodeConstraints( + ModelPart& rModelPart, + array_1d SurfaceNormal, + const std::unordered_map rDuplicateNodeMap, + Node::Pointer pControlNode, + const std::unordered_set*> rPositiveSideDofs) const = 0; + + struct Impl; + std::unique_ptr mpImpl; + +private: + InsertPreTensionOperation(const InsertPreTensionOperation&) = delete; + + InsertPreTensionOperation& operator=(const InsertPreTensionOperation&) = delete; +}; // class InsertPreTensionOperation + + +class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertDirichletPreTensionOperation final : public InsertPreTensionOperation { +public: + KRATOS_CLASS_POINTER_DEFINITION(InsertDirichletPreTensionOperation); + + using InsertPreTensionOperation::InsertPreTensionOperation; + + std::string Info() const override; + +protected: + /// @details Inserts a constraint that ties the average out-of-plane + /// relative displacement of duplicated nodes to a prescribed value. + void InsertControlNodeConstraints( + ModelPart& rModelPart, + array_1d SurfaceNormal, + const std::unordered_map rDuplicateNodeMap, + Node::Pointer pControlNode, + const std::unordered_set*> rPositiveSideDofs) const override; +}; // class InsertDirichletPreTensionOperation + + +class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertNeumannPreTensionOperation final : public InsertPreTensionOperation { +public: + KRATOS_CLASS_POINTER_DEFINITION(InsertNeumannPreTensionOperation); + + using InsertPreTensionOperation::InsertPreTensionOperation; + + std::string Info() const override; + +protected: + /// @details Loads the average out-of-plane displacement with a prescribed + /// value on both sides of the pre-tension surface in opposite + /// directions. + void InsertControlNodeConstraints( + ModelPart& rModelPart, + array_1d SurfaceNormal, + const std::unordered_map rDuplicateNodeMap, + Node::Pointer pControlNode, + const std::unordered_set*> rPositiveSideDofs) const override; +}; // class InsertNeumannPreTensionOperation + + +} // namespace Kratos diff --git a/applications/StructuralMechanicsApplication/custom_operations/insert_pretension_operation.cpp b/applications/StructuralMechanicsApplication/custom_operations/insert_pretension_operation.cpp deleted file mode 100644 index 08e425fa3a79..000000000000 --- a/applications/StructuralMechanicsApplication/custom_operations/insert_pretension_operation.cpp +++ /dev/null @@ -1,507 +0,0 @@ -// KRATOS ___| | | | -// \___ \ __| __| | | __| __| | | __| _` | | -// | | | | | ( | | | | ( | | -// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS -// -// License: BSD License -// license: StructuralMechanicsApplication/license.txt -// -// Main authors: Máté Kelemen -// - -// --- Structural Includes --- -#include "custom_operations/insert_pretension_operation.hpp" // InsertPretensionOperation - -// --- Core Includes --- -#include "includes/master_slave_constraint.h" // MasterSlaveConstraint -#include "includes/model_part.h" // ModelPart -#include "processes/find_global_nodal_elemental_neighbours_process.h" -#include "utilities/comparison.h" // Comparison -#include "includes/element.h" // Element -//#include "geometries/point_3d.h" // Point3D - -// --- STL Includes --- -#include // std::array -#include // std::stringstream - - -namespace Kratos { - - -//class PretensionElement : public Element { -// -//}; // class PretensionElement - - -struct InsertPretensionOperation::Impl { - ModelPart* mpModelPart; - std::string mImpositionName; - double mPretensionValue; -}; // struct InsertPretensionOperation::Impl - - -InsertPretensionOperation::InsertPretensionOperation(Model& rModel, Parameters Settings) - : Operation(), - mpImpl(new Impl) { - Settings.ValidateAndAssignDefaults(this->GetDefaultParameters()); - - KRATOS_TRY - mpImpl->mpModelPart = &rModel.GetModelPart(Settings["model_part_name"].GetString()); - mpImpl->mImpositionName = Settings["imposition"].GetString(); - mpImpl->mPretensionValue = Settings["pretension_value"].GetDouble(); - KRATOS_CATCH("") - - const std::array valid_imposition_names { - std::string {"lagrange_element"}, - std::string {"penalty_element"}, - std::string {"constraint"} - }; - - if (std::find(valid_imposition_names.begin(), - valid_imposition_names.end(), - mpImpl->mImpositionName) == valid_imposition_names.end()) { - std::stringstream message; - message << "Invalid setting for \"imposition\". Options are: "; - for (const auto& r_imposition_name : valid_imposition_names) - message << '"' << r_imposition_name << "\" "; - KRATOS_ERROR << message.str(); - } - - // This operation expects the pre-tension surface to be defined - // by a set of Conditions in the input ModelPart. It requires - // at least one Condition. - KRATOS_ERROR_IF(mpImpl->mpModelPart->Conditions().empty()) - << "InsertPretensionOperation requires at least one Condition in the provided ModelPart, " - << "but there aren't any in '" << mpImpl->mpModelPart->Name() << "'."; -} - - -InsertPretensionOperation::InsertPretensionOperation(InsertPretensionOperation&&) noexcept = default; - - -InsertPretensionOperation::~InsertPretensionOperation() = default; - - -struct PretensionSurfacePartition { - array_1d normal; - std::vector positive_side, negative_side; -}; // struct SurfacePartition - - -void FindElementsAdjacentToConditions(ModelPart& rModelPart) { - using NeighborElements = decltype(NEIGHBOUR_ELEMENTS)::Type; - - // Make sure every node on the pre-tension surface is in the model part. - KRATOS_TRY - ModelPart::NodesContainerType nodes; - for (Condition& r_condition : rModelPart.Conditions()) { - auto& r_geometry = r_condition.GetGeometry(); - nodes.insert(r_geometry.begin(), r_geometry.end()); - } // for r_condition in rModelPart.Conditions - rModelPart.Nodes().insert(nodes); - KRATOS_CATCH("") - - KRATOS_TRY - // Find elements containing each node in the model part. - FindGlobalNodalElementalNeighboursProcess(rModelPart.GetRootModelPart()).Execute(); - - block_for_each( - rModelPart.Conditions(), - [] (Condition& r_condition) -> void { - NeighborElements neighbor_elements; - - // Collect all neighbor elements from the condition's nodes. - for (Node& r_node : r_condition.GetGeometry()) { - NeighborElements& r_containing_elements = r_node.GetValue(NEIGHBOUR_ELEMENTS); - for (unsigned i_element=0u; i_element bool { - return rp_other_element->Id() == rp_element->Id(); - }); - if (it_element == neighbor_elements.ptr_end()) { - neighbor_elements.push_back(rp_element); - } - } - } - - // Filter elements that don't contain all nodes of the condition. - neighbor_elements.erase(std::remove_if( - neighbor_elements.begin(), - neighbor_elements.end(), - [&r_condition] (const Element& r_element) -> bool { - for (const Node& r_condition_node : r_condition.GetGeometry()) { - const auto it_node = std::find_if( - r_element.GetGeometry().begin(), - r_element.GetGeometry().end(), - [&r_condition_node] (const Node& r_element_node) -> bool { - return r_element_node.Id() == r_condition_node.Id(); - }); - if (it_node == r_element.GetGeometry().end()) - return true; - } - return false; - } - ), neighbor_elements.end()); - - // Set neighbor elements on the condition. - r_condition.SetValue(NEIGHBOUR_ELEMENTS, neighbor_elements); - }); - KRATOS_CATCH("") -} - - -template -PretensionSurfacePartition PretensionPlanePartition(const ModelPart::ConditionsContainerType::iterator itConditionBegin, - const ModelPart::ConditionsContainerType::iterator itConditionEnd) { - PretensionSurfacePartition pretension_surface; - - // There are three valid possibilities here. - // 1) All surfaces are 0-dimensional (corners) - // => expecting exactly 1 condition - // => expecting exactly 2 connected elements, both - // must have line geometries - // => the normal is not well-defined, so it is chosen - // as the average of the connected elements' directions - // 2) All surfaces are 1-dimensional (edges) - // => expecting at least 1 condition - // => expecting at least 2 connected elements, all of which - // must be 2-dimensional - // => the normal is not well-defined, so it is chosen as the - // average of the rotated endpoints in 2D, in positive - // direction. This assumes coordinates in z-direction - // to vanish. This assumption is checked. - // 3) All surfaces are 2-dimensional (surfaces) - // => expecting at least 1 condition - // => expecting at least 2 connected elements, all - // of which must be 3-dimensional - // => normals are well-defined - // In all cases, at least one connected element is required on - // either side of the average pre-tension plane. - - // Sanity checks. - if constexpr (SurfaceDimension == 0u) { - KRATOS_ERROR_IF_NOT(std::distance(itConditionBegin, itConditionEnd) == 1) - << "Expecting exactly 1 condition on the 0-dimensional pre-tension surface , but got " - << std::distance(itConditionBegin, itConditionEnd) << "."; - } else if constexpr (SurfaceDimension == 1u) { - KRATOS_ERROR_IF_NOT(1 <= std::distance(itConditionBegin, itConditionEnd)) - << "Expecting at least 1 condition on the 1-dimensional pre-tension surface , but got " - << std::distance(itConditionBegin, itConditionEnd) << "."; - } else if constexpr (SurfaceDimension == 2u) { - KRATOS_ERROR_IF_NOT(1 <= std::distance(itConditionBegin, itConditionEnd)) - << "Expecting at least 1 condition on the 2-dimensional pre-tension surface , but got " - << std::distance(itConditionBegin, itConditionEnd) << "."; - } else { - static_assert(float(SurfaceDimension) == 0.5f, "Invalid surface dimension"); - } - - for (auto it_condition=itConditionBegin; it_condition!=itConditionEnd; ++it_condition) { - Condition& r_condition = *it_condition; - - // Sanity checks on the condition. - KRATOS_ERROR_IF_NOT(r_condition.GetGeometry().LocalSpaceDimension() == SurfaceDimension) - << "Expecting all geometries defining the pre-tension surface to be " << SurfaceDimension << " dimensional, " - << "but condition " << r_condition.Id() << " on geometry " << r_condition.GetGeometry().Id() << " (" << r_condition.GetGeometry().Name() << ") " - << "is " << r_condition.GetGeometry().LocalSpaceDimension() << "-dimensional."; - - KRATOS_ERROR_IF_NOT(r_condition.Has(NEIGHBOUR_ELEMENTS)) - << "Expecting condition " << r_condition.Id() << " to have exactly " - << "2 neighboring elements, but found none."; - auto& r_neighbor_elements = r_condition.GetValue(NEIGHBOUR_ELEMENTS); - if (r_neighbor_elements.size() != 2) { - std::stringstream message; - message << "Expecting condition " << r_condition.Id() << " to have exactly " - << "2 neighboring elements, but found " << r_neighbor_elements.size() << " ["; - for (const Element& r_element : r_neighbor_elements) message << r_element.Id() << " "; - message << "]."; - KRATOS_ERROR << message.str(); - } - - // Check whether the 2 connected elements lie on opposite sides of the surface. - const array_1d - first_connected_direction = r_neighbor_elements[0].GetGeometry().Center() - r_condition.GetGeometry().Center(), - second_connected_direction = r_neighbor_elements[1].GetGeometry().Center() - r_condition.GetGeometry().Center(); - const double direction_inner_product = std::inner_product( - first_connected_direction.begin(), - first_connected_direction.end(), - second_connected_direction.begin(), - 0.0); - KRATOS_ERROR_IF_NOT(direction_inner_product < 0.0) - << "Element " << r_neighbor_elements[0].Id() << " and " << r_neighbor_elements[1].Id() << " " - << "do not lie on opposite sides of condition " << r_condition.Id() << " " - << "defining a pre-tension surface."; - - if constexpr (SurfaceDimension == 0u) { - for (unsigned i_element=0u; i_elementGetGeometry()[0].Id(); - const array_1d surface_node_position = itConditionBegin->GetGeometry()[0]; - array_1d opposite_node_position; - if (r_element.GetGeometry()[0].Id() == surface_node_id) { - opposite_node_position = r_element.GetGeometry()[1]; - } else if (r_element.GetGeometry()[1].Id() == surface_node_id) { - opposite_node_position = r_element.GetGeometry()[0]; - } else { - KRATOS_ERROR << "Element " << r_element.Id() << " is assumed to be connected to node " - << surface_node_id << " but is not."; - } - - if (i_element % 2) { - pretension_surface.positive_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); - pretension_surface.normal = opposite_node_position - surface_node_position; - } else { - pretension_surface.negative_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); - } - } // for i_element in range(r_neighbor_elements.size()) - } /*if SurfaceDimension == 0*/ else if constexpr (SurfaceDimension == 1u) { - // Define the line condition's plane. - const array_1d surface_begin = r_condition.GetGeometry()[0], - surface_end = *(r_condition.GetGeometry().end() - 1); - KRATOS_ERROR_IF_NOT(Comparison::Equal(1e-12, 1e-9)(surface_begin[2], 0.0)) - << "Assuming that 1-dimensional pre-tension surfaces lie on the XY plane, but node " - << r_condition.GetGeometry()[0].Id() << " of condition " << r_condition.Id() << " " - << "has a Z-coordinate of " << surface_begin[2] << "."; - KRATOS_ERROR_IF_NOT(Comparison::Equal(1e-12, 1e-9)(surface_end[2], 0.0)) - << "Assuming that 1-dimensional pre-tension surfaces lie on the XY plane, but node " - << r_condition.GetGeometry()[1].Id() << " of condition " << r_condition.Id() << " " - << "has a Z-coordinate of " << surface_end[2] << "."; - const array_1d - surface_normal { - -(surface_end[1] - surface_begin[1]), - surface_end[0] - surface_begin[0], - 0.0}, - surface_center = r_condition.GetGeometry().Center(); - pretension_surface.normal += surface_normal; - - // Sort connected elements. - for (unsigned i_element=0u; i_element element_direction = r_element.GetGeometry().Center() - surface_center; - const double direction_inner_product = std::inner_product( - surface_normal.begin(), - surface_normal.end(), - element_direction.begin(), - 0.0); - if (0 < direction_inner_product) pretension_surface.negative_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); - else pretension_surface.positive_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); - } // for i_element in range(r_neighbor_elements.size()) - - KRATOS_ERROR_IF_NOT(pretension_surface.positive_side.size() == pretension_surface.negative_side.size()) - << "Elements " << r_neighbor_elements[0].Id() << " and " << r_neighbor_elements[1].Id() << " " - << "lie on the same side of condition " << r_condition.Id() << " that defines a pre-tension surface."; - } /*else if SurfaceDimension == 1*/ else if constexpr (SurfaceDimension == 2u) { - // Define the surface's plane. - const array_1d surface_center = r_condition.GetGeometry().Center(); - const array_1d surface_normal = r_condition.GetGeometry().Normal(surface_center); - pretension_surface.normal += surface_normal; - - // Sort connected elements. - for (unsigned i_element=0u; i_element element_direction = r_element.GetGeometry().Center() - surface_center; - const double direction_inner_product = std::inner_product( - surface_normal.begin(), - surface_normal.end(), - element_direction.begin(), - 0.0); - if (0 < direction_inner_product) pretension_surface.negative_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); - else pretension_surface.positive_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); - } // for i_element in range(r_neighbor_elements.size()) - - KRATOS_ERROR_IF_NOT(pretension_surface.positive_side.size() == pretension_surface.negative_side.size()) - << "Elements " << r_neighbor_elements[0].Id() << " and " << r_neighbor_elements[1].Id() << " " - << "lie on the same side of condition " << r_condition.Id() << " that defines a pre-tension surface."; - } /*else if SurfaceDimension == 2*/ - } // for r_condition in r_model_part.Conditions - - return pretension_surface; -} - - -void InsertPretensionOperation::Execute() { - ModelPart& r_model_part = *mpImpl->mpModelPart; - - // Find elements connected to the pre-tension surface. - // => each node will have a NEIGHBOR_ELEMENTS variable - // in its non-historical storage that lists elements - // containing them. - KRATOS_TRY - FindElementsAdjacentToConditions(r_model_part); - KRATOS_CATCH("") - - // Sort connected elements into 2 categories, lying on either side of - // the pre-tension surface. - PretensionSurfacePartition pretension_surface; - KRATOS_TRY - // Find the dimension of the pre-tension surface, - // and make sure all conditions share it. - KRATOS_ERROR_IF(r_model_part.Conditions().empty()); - const int pre_tension_surface_dimension = r_model_part.Conditions().front().GetGeometry().LocalSpaceDimension(); - switch (pre_tension_surface_dimension) { - case 0u: { - pretension_surface = PretensionPlanePartition<0>(r_model_part.Conditions().begin(), r_model_part.Conditions().end()); - break; - } - case 1u: { - pretension_surface = PretensionPlanePartition<1>(r_model_part.Conditions().begin(), r_model_part.Conditions().end()); - break; - } - case 2u: { - pretension_surface = PretensionPlanePartition<2>(r_model_part.Conditions().begin(), r_model_part.Conditions().end()); - break; - } - default: KRATOS_ERROR << "Invalid pre-tension surface dimension: " << pre_tension_surface_dimension << "."; - } - KRATOS_CATCH("") - - // Duplicate nodes on the pre-tension surface and collect relevant displacement components. - std::unordered_map duplicated_nodes; - KRATOS_TRY - std::unordered_set surface_nodes; - - // Collect nodes on the surface. - for (Condition& r_condition : r_model_part.Conditions()) { - for (Node& r_node : r_condition.GetGeometry()) { - surface_nodes.insert(&r_node); - } // for r_node in r_condition.GetGeometry() - } // for r_condition in r_model_part.Conditions() - - // Find the first available node ID that can safely be - // used to insert new ones. - Node::IndexType node_id = r_model_part.GetRootModelPart().Nodes().back().Id() + 1; - - // Clone surface nodes and associate them with the original ones. - for (Node* p_node : surface_nodes) { - const auto emplace_result = duplicated_nodes.emplace( - p_node->Id(), - p_node->Clone(node_id++)); - r_model_part.AddNode(emplace_result.first->second); - } // for p_node in surface_nodes - - // Replace nodes on the negative side elements with duplicated ones. - for (auto& rp_element : pretension_surface.negative_side) { - for (auto itp_node=rp_element->GetGeometry().ptr_begin(); itp_node!=rp_element->GetGeometry().ptr_end(); ++itp_node) { - const auto it_duplicate_pair = duplicated_nodes.find((*itp_node)->Id()); - if (it_duplicate_pair != duplicated_nodes.end()) { - (*itp_node) = it_duplicate_pair->second; - } - } // for itp_node in rp_element->GetGeometry() - } // for rp_element in pretension_surface.negative_side - KRATOS_CATCH("") - - // Construct the constraint equation. - struct PretensionConstraintData { - std::vector::Pointer> dofs; - std::vector gradient; - double gap; - } constraint_data; - KRATOS_TRY - const std::array*,3> all_displacement_components { - &DISPLACEMENT_X, - &DISPLACEMENT_Y, - &DISPLACEMENT_Z - }; - - Condition::DofsVectorType dofs; - for (Condition& r_condition : r_model_part.Conditions()) { - r_condition.GetDofList(dofs, r_model_part.GetProcessInfo()); - for (Dof::Pointer p_negative_side_dof : dofs) { - const auto it_displacement_component = std::find_if( - all_displacement_components.begin(), - all_displacement_components.end(), - [p_negative_side_dof](const Variable* p_variable) { - return p_variable->Key() == p_negative_side_dof->GetVariable().Key(); - }); - if (it_displacement_component != all_displacement_components.end()) { - // Insert the negative side component. - const std::size_t i_displacement_component = std::distance(all_displacement_components.begin(), it_displacement_component); - - { - const double gradient_component = pretension_surface.normal[i_displacement_component]; - constraint_data.dofs.push_back(p_negative_side_dof); - constraint_data.gradient.push_back(gradient_component); - } - - // Find the duplicate node on the positive side. - { - const Node::IndexType negative_side_node_id = p_negative_side_dof->Id(); - const auto it_positive_side_node = duplicated_nodes.find(negative_side_node_id); - KRATOS_ERROR_IF(it_positive_side_node == duplicated_nodes.end()); - Node& r_positive_side_node = *it_positive_side_node->second; - const auto& r_positive_side_dofs = r_positive_side_node.GetDofs(); - const auto it_positive_side_dof = std::find_if( - r_positive_side_dofs.begin(), - r_positive_side_dofs.end(), - [it_displacement_component](const auto& rp_dof) { - return rp_dof->GetVariable().Key() == (*it_displacement_component)->Key(); - }); - KRATOS_ERROR_IF(it_positive_side_dof == r_positive_side_dofs.end()); - const double gradient_component = -pretension_surface.normal[i_displacement_component]; - constraint_data.dofs.push_back(it_positive_side_dof->get()); - constraint_data.gradient.push_back(gradient_component); - } - } // if p_dof in all_displacement_components - } // for p_negative_side_dof in dofs - } // for r_condition in r_model_part.Conditions - - constraint_data.gap = mpImpl->mPretensionValue; - KRATOS_CATCH("") - - // Insert the constraint. - const std::size_t constraint_id = r_model_part.MasterSlaveConstraints().empty() ? - 1ul : - r_model_part.MasterSlaveConstraints().back().Id() + 1; - - MasterSlaveConstraint::DofPointerVectorType masters, slaves; - MasterSlaveConstraint::MatrixType relation_matrix; - MasterSlaveConstraint::VectorType constants; - - slaves.push_back(constraint_data.dofs.back()); - masters.reserve(constraint_data.dofs.size() - 1); - std::copy( - constraint_data.dofs.begin(), - constraint_data.dofs.end() - 1, - std::back_inserter(masters)); - - relation_matrix.resize(1, constraint_data.dofs.size() - 1, false); - std::transform( - constraint_data.gradient.begin(), - constraint_data.gradient.end() - 1, - relation_matrix.begin2(), - [&constraint_data](double component){return -component / constraint_data.gradient.back();}); - - constants.resize(1, false); - constants[0] = -constraint_data.gap; - - KRATOS_ERROR_IF_NOT(KratosComponents::Has("LinearMasterSlaveConstraint")); - MasterSlaveConstraint::Pointer p_constraint; - - KRATOS_TRY - p_constraint = KratosComponents::Get("LinearMasterSlaveConstraint").Create( - constraint_id, - masters, - slaves, - relation_matrix, - constants); - r_model_part.AddMasterSlaveConstraint(p_constraint); - KRATOS_CATCH("") -} - - -const Parameters InsertPretensionOperation::GetDefaultParameters() const { - return Parameters(R"({ - "model_part_name" : "", - "pretension_value" : 0.0, - "imposition" : "lagrange_element" - })"); -} - - -} // namespace Kratos diff --git a/applications/StructuralMechanicsApplication/custom_operations/insert_pretension_operation.hpp b/applications/StructuralMechanicsApplication/custom_operations/insert_pretension_operation.hpp deleted file mode 100644 index 73ada50485fa..000000000000 --- a/applications/StructuralMechanicsApplication/custom_operations/insert_pretension_operation.hpp +++ /dev/null @@ -1,46 +0,0 @@ -// KRATOS ___| | | | -// \___ \ __| __| | | __| __| | | __| _` | | -// | | | | | ( | | | | ( | | -// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS -// -// License: BSD License -// license: StructuralMechanicsApplication/license.txt -// -// Main authors: Máté Kelemen -// - -#pragma once - -// --- Core Includes --- -#include "operations/operation.h" // Operation -#include "containers/model.h" // Model - - -namespace Kratos { - - -class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertPretensionOperation final : public Operation { -public: - KRATOS_CLASS_POINTER_DEFINITION(InsertPretensionOperation); - - InsertPretensionOperation(Model& rModel, Parameters Settings); - - InsertPretensionOperation(InsertPretensionOperation&&) noexcept; - - ~InsertPretensionOperation(); - - void Execute() override; - - const Parameters GetDefaultParameters() const override; - -private: - InsertPretensionOperation(const InsertPretensionOperation&) = delete; - - InsertPretensionOperation& operator=(const InsertPretensionOperation&) = delete; - - struct Impl; - std::unique_ptr mpImpl; -}; // class InsertPretensionOperation - - -} // namespace Kratos diff --git a/applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.cpp b/applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.cpp index 61551b8c3e4b..da7bedfe5c3e 100644 --- a/applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.cpp +++ b/applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.cpp @@ -1,15 +1,28 @@ // --- Structural Includes --- #include "custom_python/add_custom_operations_to_python.hpp" // AddCustomOperationsToPython -#include "custom_operations/insert_pretension_operation.hpp" // InsertPretensionOperation +#include "custom_operations/insert_pre_tension_operation.hpp" // InsertNeumannPreTensionOperation, InsertDirichletPreTensionOperation namespace Kratos::Python { void AddCustomOperationsToPython(pybind11::module& rModule) { - pybind11::class_(rModule, "InsertPretensionOperation") - .def(pybind11::init()) - ; + pybind11::class_< + InsertDirichletPreTensionOperation, + InsertDirichletPreTensionOperation::Pointer, + Operation>( + rModule, + "InsertDirichletPreTensionOperation") + .def(pybind11::init()) + ; + pybind11::class_< + InsertNeumannPreTensionOperation, + InsertNeumannPreTensionOperation::Pointer, + Operation>( + rModule, + "InsertNeumannPreTensionOperation") + .def(pybind11::init()) + ; } diff --git a/applications/StructuralMechanicsApplication/python_scripts/insert_pretension_process.py b/applications/StructuralMechanicsApplication/python_scripts/dirichlet_pre_tension_process.py similarity index 79% rename from applications/StructuralMechanicsApplication/python_scripts/insert_pretension_process.py rename to applications/StructuralMechanicsApplication/python_scripts/dirichlet_pre_tension_process.py index 6b9b27044397..a865e40cd920 100644 --- a/applications/StructuralMechanicsApplication/python_scripts/insert_pretension_process.py +++ b/applications/StructuralMechanicsApplication/python_scripts/dirichlet_pre_tension_process.py @@ -5,14 +5,14 @@ from KratosMultiphysics import StructuralMechanicsApplication as StructuralMechanics -class PretensionProcess(KratosMultiphysics.Process): +class DirichletPretensionProcess(KratosMultiphysics.Process): def __init__(self, model: KratosMultiphysics.Model, parameters: KratosMultiphysics.Parameters) -> None: + super().__init__() parameters.ValidateAndAssignDefaults(self.GetDefaultParameters()) - #self.__model_part = model.GetModelPart(parameters["model_part_name"].GetString()) - self.__insert_pretension_operation = StructuralMechanics.InsertPretensionOperation( + self.__insert_pretension_operation = StructuralMechanics.InsertDirichletPreTensionOperation( model, parameters) @@ -22,12 +22,13 @@ def ExecuteBeforeSolutionLoop(self) -> None: def GetDefaultParameters(self) -> KratosMultiphysics.Parameters: return KratosMultiphysics.Parameters(R"""{ "model_part_name" : "", - "pretension_value" : 0.0 + "magnitude" : 0.0, + "verbosity" : 1 }""") def Factory(parameters: KratosMultiphysics.Parameters, model: KratosMultiphysics.Model) -> KratosMultiphysics.Process: - return PretensionProcess( + return DirichletPretensionProcess( model, parameters["Parameters"]) diff --git a/applications/StructuralMechanicsApplication/python_scripts/neumann_pre_tension_process.py b/applications/StructuralMechanicsApplication/python_scripts/neumann_pre_tension_process.py new file mode 100644 index 000000000000..da2c16ba8f35 --- /dev/null +++ b/applications/StructuralMechanicsApplication/python_scripts/neumann_pre_tension_process.py @@ -0,0 +1,34 @@ +# --- Core Imports --- +import KratosMultiphysics + +# --- Structural Imports --- +from KratosMultiphysics import StructuralMechanicsApplication as StructuralMechanics + + +class NeumannPretensionProcess(KratosMultiphysics.Process): + + def __init__(self, + model: KratosMultiphysics.Model, + parameters: KratosMultiphysics.Parameters) -> None: + super().__init__() + parameters.ValidateAndAssignDefaults(self.GetDefaultParameters()) + self.__insert_pre_tension_operation = StructuralMechanics.InsertNeumannPreTensionOperation( + model, + parameters) + + def ExecuteBeforeSolutionLoop(self) -> None: + self.__insert_pre_tension_operation.Execute() + + def GetDefaultParameters(self) -> KratosMultiphysics.Parameters: + return KratosMultiphysics.Parameters(R"""{ + "model_part_name" : "", + "magnitude" : 0.0, + "verbosity" : 1 + }""") + + +def Factory(parameters: KratosMultiphysics.Parameters, + model: KratosMultiphysics.Model) -> KratosMultiphysics.Process: + return NeumannPretensionProcess( + model, + parameters["Parameters"]) diff --git a/applications/StructuralMechanicsApplication/tests/constraints/pre_tensioned_beam.json b/applications/StructuralMechanicsApplication/tests/constraints/pre_tensioned_beam.json new file mode 100644 index 000000000000..55e4afdb8c6d --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/constraints/pre_tensioned_beam.json @@ -0,0 +1,124 @@ +{ + "analysis_stage": "KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis", + "problem_data": { + "problem_name": "pre_tensioned_beam", + "parallel_type": "OpenMP", + "echo_level": 1, + "start_time": 0, + "end_time": 1 + }, + "solver_settings": { + "time_stepping": { + "time_step": 1.0 + }, + "solver_type": "Static", + "model_part_name": "root", + "domain_size": 2, + "echo_level": 0, + + "analysis_type": "non_linear", + "max_iteration" : 1e2, + "convergence_criterion": "displacement_criterion", + "displacement_relative_tolerance": 1e-4, + "displacement_absolute_tolerance": 1e-9, + + "rotation_dofs" : true, + "move_mesh_flag" : true, + + "model_import_settings": { + "input_type": "mdpa", + "input_filename" : "pre_tensioned_beam" + }, + + "material_import_settings": { + "materials_filename": "pre_tensioned_beam_materials.json" + }, + + "builder_and_solver_settings" : { + "type" : "p_multigrid", + "advanced_settings" : { + "constraint_imposition_settings" : { + "method" : "augmented_lagrange", + "max_iterations" : 20, + "penalty_factor" : "max", + "tolerance" : 1e-9, + "verbosity" : 1 + } + } + } + }, + "processes": { + "constraints_process_list": [ + { + "python_module": "assign_vector_variable_process", + "kratos_module": "KratosMultiphysics", + "process_name": "AssignVectorVariableProcess", + "Parameters": { + "model_part_name": "root.dirichlet", + "variable_name": "DISPLACEMENT", + "interval": [0, "End"], + "constrained": [true, true, false], + "value": [0, 0, 0] + } + }, + { + "python_module": "assign_vector_variable_process", + "kratos_module": "KratosMultiphysics", + "process_name": "AssignVectorVariableProcess", + "Parameters": { + "model_part_name": "root.dirichlet", + "variable_name": "ROTATION", + "interval": [0, "End"], + "constrained": [false, false, true], + "value": [null, null, 0] + } + }, + { + "python_module" : "neumann_pre_tension_process", + "kratos_module" : "KratosMultiphysics.StructuralMechanicsApplication", + "Parameters" : { + "model_part_name" : "root.pre_tension_surface", + "magnitude" : 2e7, + "verbosity" : 1 + } + } + ], + "loads_process_list": [ + { + "python_module": "assign_vector_by_direction_to_condition_process", + "kratos_module": "KratosMultiphysics", + "process_name": "AssignVectorByDirectionToConditionProcess", + "Parameters": { + "model_part_name": "root.neumann", + "variable_name": "POINT_LOAD", + "interval": [ + 0, "End" + ], + "modulus": 0, + "direction": [-2, 7, 0] + } + } + ] + } + //,"output_processes": { + //"vtk_output": [ + // { + // "Parameters": { + // "folder_name": "vtk", + // "model_part_name": "root", + // "nodal_data_value_variables": [], + // "nodal_solution_step_data_variables": [ + // "DISPLACEMENT", + // "ROTATION", + // "REACTION" + // ], + // "output_sub_model_parts": false, + // "save_output_files_in_folder": true + // }, + // "kratos_module": "KratosMultiphysics", + // "process_name": "VTKOutputProcess", + // "python_module": "vtk_output_process" + // } + //] + //} +} diff --git a/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.mdpa b/applications/StructuralMechanicsApplication/tests/constraints/pre_tensioned_beam.mdpa similarity index 64% rename from applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.mdpa rename to applications/StructuralMechanicsApplication/tests/constraints/pre_tensioned_beam.mdpa index 1c4188693421..6e2306a44233 100644 --- a/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.mdpa +++ b/applications/StructuralMechanicsApplication/tests/constraints/pre_tensioned_beam.mdpa @@ -2,15 +2,19 @@ Begin Properties 0 End Properties Begin Nodes - 1 0 0 0 - 2 1 0 0 - 3 2 0 0 - 4 0 2 0 - 5 1 2 0 - 6 2 2 0 - 7 1 1 0 + 1 0.0 2.0 0.0 + 2 1.0 1.0 0.0 + 3 2.0 0.0 0.0 + 4 0.5 4.0 0.0 + 5 3.0 4.0 0.0 + 6 4.5 4.0 0.0 + 7 2.0 2.5 0.0 End Nodes +Begin Geometries Point2D + 1 7 +End Geometries + Begin Elements CrBeamElement2D2N 1 0 1 2 2 0 2 3 @@ -21,7 +25,7 @@ Begin Elements CrBeamElement2D2N End Elements Begin Conditions PointLoadCondition2D1N - 1 0 7 + 1 0 5 End Conditions Begin SubModelPart neumann @@ -39,8 +43,8 @@ Begin SubModelPart dirichlet End SubModelPartNodes End SubModelPart -Begin SubModelPart pretension_surface - Begin SubModelPartConditions +Begin SubModelPart pre_tension_surface + Begin SubModelPartGeometries 1 - End SubModelPartConditions + End SubModelPartGeometries End SubModelPart diff --git a/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d_materials.json b/applications/StructuralMechanicsApplication/tests/constraints/pre_tensioned_beam_materials.json similarity index 59% rename from applications/StructuralMechanicsApplication/tests/constraints/pretension_1d_materials.json rename to applications/StructuralMechanicsApplication/tests/constraints/pre_tensioned_beam_materials.json index 1c61b4120912..d5f25791474b 100644 --- a/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d_materials.json +++ b/applications/StructuralMechanicsApplication/tests/constraints/pre_tensioned_beam_materials.json @@ -8,11 +8,11 @@ "name": "BeamConstitutiveLaw" }, "Variables": { - "DENSITY": 1e0, - "YOUNG_MODULUS": 1e0, - "POISSON_RATIO": 3e-1, - "CROSS_AREA": 1e-3, - "I33": 1e-5 + "DENSITY": 7.85e3, + "YOUNG_MODULUS": 2.069e11, + "POISSON_RATIO": 2.9e-1, + "CROSS_AREA": 7.68e-3, + "I33" : 7.763e-5 } } } diff --git a/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.json b/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.json deleted file mode 100644 index 21c74d374031..000000000000 --- a/applications/StructuralMechanicsApplication/tests/constraints/pretension_1d.json +++ /dev/null @@ -1,92 +0,0 @@ -{ - "analysis_stage": "KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis", - "problem_data": { - "problem_name": "test_pretension_1d", - "parallel_type": "OpenMP", - "echo_level": 1, - "start_time": 0, - "end_time": 1 - }, - "solver_settings": { - "time_stepping": { - "time_step": 0.5 - }, - "solver_type": "Static", - "model_part_name": "root", - "domain_size": 2, - "echo_level": 0, - "analysis_type": "linear", - "move_mesh_flag" : false, - "model_import_settings": { - "input_type": "mdpa", - "input_filename" : "pretension_1d" - }, - "material_import_settings": { - "materials_filename": "pretension_1d_materials.json" - }, - "convergence_criterion": "displacement_criterion", - "displacement_relative_tolerance": 1e-4, - "displacement_absolute_tolerance": 1e-9, - "rotation_dofs" : true - }, - "processes": { - "constraints_process_list": [ - { - "python_module": "assign_vector_variable_process", - "kratos_module": "KratosMultiphysics", - "process_name": "AssignVectorVariableProcess", - "Parameters": { - "model_part_name": "root.dirichlet", - "variable_name": "DISPLACEMENT", - "interval": [0, "End"], - "constrained": [true, true, false], - "value": [0, 0, 0] - } - }, - { - "python_module" : "insert_pretension_process", - "kratos_module" : "KratosMultiphysics.StructuralMechanicsApplication", - "Parameters" : { - "model_part_name" : "root.pretension_surface", - "pretension_value" : 5e-1 - } - } - ], - "loads_process_list": [ - { - "python_module": "assign_vector_by_direction_to_condition_process", - "kratos_module": "KratosMultiphysics", - "process_name": "AssignVectorByDirectionToConditionProcess", - "Parameters": { - "model_part_name": "root.neumann", - "variable_name": "POINT_LOAD", - "interval": [ - 0, "End" - ], - "modulus": 0.0, - "direction": [-2, 7, 0] - } - } - ] - }, - "output_processes": { - "vtk_output": [ - { - "Parameters": { - "folder_name": "vtk", - "model_part_name": "root", - "nodal_data_value_variables": [], - "nodal_solution_step_data_variables": [ - "DISPLACEMENT", - "REACTION" - ], - "output_sub_model_parts": false, - "save_output_files_in_folder": true - }, - "kratos_module": "KratosMultiphysics", - "process_name": "VTKOutputProcess", - "python_module": "vtk_output_process" - } - ] - } -} diff --git a/applications/StructuralMechanicsApplication/tests/test_pretension.py b/applications/StructuralMechanicsApplication/tests/test_pre_tension.py similarity index 86% rename from applications/StructuralMechanicsApplication/tests/test_pretension.py rename to applications/StructuralMechanicsApplication/tests/test_pre_tension.py index cef9e70f547d..31e7232e2db7 100644 --- a/applications/StructuralMechanicsApplication/tests/test_pretension.py +++ b/applications/StructuralMechanicsApplication/tests/test_pre_tension.py @@ -10,10 +10,10 @@ class TestPretension(KratosMultiphysics.KratosUnittest.TestCase): - def test_pretension_1d(self) -> None: + def test_pretensioned_beam(self) -> None: with WorkFolderScope("constraints", pathlib.Path(__file__).absolute()): # Load config - with open("pretension_1d.json", "r") as project_parameters_file: + with open("pre_tensioned_beam.json", "r") as project_parameters_file: parameters = KratosMultiphysics.Parameters(project_parameters_file.read()) model = KratosMultiphysics.Model() analysis = StructuralMechanicsAnalysis(model, parameters) From c0b42d78ff7c56ad8b6b50e7274e7e9cdcd48d54 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Thu, 23 Apr 2026 15:52:28 +0200 Subject: [PATCH 10/18] add missing STL includes --- .../custom_operations/insert_pre_tension_operation.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp index 218913c3c162..cf702183606c 100644 --- a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp +++ b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp @@ -26,6 +26,15 @@ // --- STL Includes --- #include // std::array #include // std::stringstream +#include // std::optional +#include // std::shared_ptr +#include // std::array +#include // std::vector +#include +#include +#include +#include +#include namespace Kratos { From 27166b2b2c345f36773f641c79babd73212555a4 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Wed, 29 Apr 2026 09:28:07 +0200 Subject: [PATCH 11/18] add condition for a single DoF --- .../point_load_condition_1d_1n.h | 114 ++++++++++++++++++ .../structural_mechanics_application.cpp | 4 +- .../structural_mechanics_application.h | 2 + 3 files changed, 119 insertions(+), 1 deletion(-) create mode 100644 applications/StructuralMechanicsApplication/custom_conditions/point_load_condition_1d_1n.h diff --git a/applications/StructuralMechanicsApplication/custom_conditions/point_load_condition_1d_1n.h b/applications/StructuralMechanicsApplication/custom_conditions/point_load_condition_1d_1n.h new file mode 100644 index 000000000000..4c69a5777a54 --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_conditions/point_load_condition_1d_1n.h @@ -0,0 +1,114 @@ +// KRATOS ___| | | | +// \___ \ __| __| | | __| __| | | __| _` | | +// | | | | | ( | | | | ( | | +// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS +// +// License: BSD License +// license: StructuralMechanicsApplication/license.txt +// +// Main authors: Máté Kelemen +// + +// --- Core Includes --- +#include "includes/condition.h" +#include "geometries/point_2d.h" +#include "includes/variables.h" + +// --- Structural Includes --- +#include "structural_mechanics_application_variables.h" // POINT_LOAD_X + + +namespace Kratos { + + +class PointLoadCondition1D1N : public Condition { +public: + KRATOS_CLASS_POINTER_DEFINITION(PointLoadCondition1D1N); + + PointLoadCondition1D1N() = default; + + PointLoadCondition1D1N( + Condition::IndexType Id, + Geometry::Pointer pGeometry, + Properties::Pointer pProperties) + : Condition(Id, pGeometry, pProperties) + {} + + Condition::Pointer Clone( + Condition::IndexType Id, + const Condition::NodesArrayType& rNodes) const override { + KRATOS_ERROR_IF_NOT(rNodes.size() == 1) + << "PointLoadCondition1D1N::Clone expects 1 node but got " << rNodes.size(); + return Condition::Pointer(new PointLoadCondition1D1N( + Id, + Geometry::Pointer(new Point2D(rNodes)), + this->pGetProperties())); + } + + void EquationIdVector( + Condition::EquationIdVectorType& rIndices, + const ProcessInfo&) const override { + rIndices.resize(1); + rIndices[0] = this->GetGeometry()[0].GetDofs()[0]->EquationId(); + } + + void GetDofList( + Condition::DofsVectorType& rDofs, + const ProcessInfo&) const override { + rDofs.resize(1); + rDofs[0] = this->GetGeometry()[0].GetDofs()[0].get(); + } + + void GetValuesVector( + Vector& rOutput, + int Step) const override { + rOutput.resize(1); + rOutput[0] = this->GetGeometry()[0].FastGetSolutionStepValue(DISPLACEMENT_X, Step); + } + + void GetFirstDerivativesVector( + Vector& rOutput, + int Step) const override { + rOutput.resize(1); + rOutput[0] = this->GetGeometry()[0].FastGetSolutionStepValue(VELOCITY_X, Step); + } + + void GetSecondDerivativesVector( + Vector& rOutput, + int Step) const override { + rOutput.resize(1); + rOutput[0] = this->GetGeometry()[0].FastGetSolutionStepValue(ACCELERATION_X, Step); + } + + void CalculateRightHandSide( + Condition::VectorType& rRhs, + const ProcessInfo&) override { + rRhs.resize(1); + rRhs[0] = this->GetValue(POINT_LOAD_X); + } + + void CalculateLocalSystem( + Condition::MatrixType& rLhs, + Condition::VectorType& rRhs, + const ProcessInfo& rProcessInfo) override { + rLhs.resize(1, 1, false); + rLhs(0, 0) = 0.0; + this->CalculateRightHandSide(rRhs, rProcessInfo); + } + + void CalculateMassMatrix( + Condition::MatrixType& rMatrix, + const ProcessInfo&) override { + rMatrix.resize(1, 1, false); + rMatrix(0, 0) = 0.0; + } + + void CalculateDampingMatrix( + Condition::MatrixType& rMatrix, + const ProcessInfo&) override { + rMatrix.resize(1, 1, false); + rMatrix(0, 0) = 0.0; + } +}; // class PointLoadCondition1D1N + +} // namespace Kratos diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp b/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp index 5abad476935f..564907bae53c 100644 --- a/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp +++ b/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp @@ -239,6 +239,7 @@ KratosStructuralMechanicsApplication::KratosStructuralMechanicsApplication() /* CONDITIONS */ // Adding point load conditions + mPointLoadCondition1D1N(0, Condition::GeometryType::Pointer(new Point2D(Condition::GeometryType::PointsArrayType(1))), Properties::Pointer(new Properties(0))), mPointLoadCondition2D1N(0, Condition::GeometryType::Pointer(new Point2D(Condition::GeometryType::PointsArrayType(1)))), mPointLoadCondition3D1N(0, Condition::GeometryType::Pointer(new Point3D(Condition::GeometryType::PointsArrayType(1)))), mPointContactCondition2D1N(0, Condition::GeometryType::Pointer(new Point2D(Condition::GeometryType::PointsArrayType(1)))), @@ -509,7 +510,7 @@ void KratosStructuralMechanicsApplication::Register() { KRATOS_REGISTER_3D_VARIABLE_WITH_COMPONENTS(MOVING_LOAD) KRATOS_REGISTER_VARIABLE(MOVING_LOAD_LOCAL_DISTANCE) KRATOS_REGISTER_VARIABLE(MOTION_TYPE) - + // Condition load variables KRATOS_REGISTER_VARIABLE(POINT_LOADS_VECTOR) KRATOS_REGISTER_VARIABLE(LINE_LOADS_VECTOR) @@ -756,6 +757,7 @@ void KratosStructuralMechanicsApplication::Register() { // Register the conditions // Point loads + KRATOS_REGISTER_CONDITION("PointLoadCondition1D1N", mPointLoadCondition1D1N) KRATOS_REGISTER_CONDITION("PointLoadCondition2D1N", mPointLoadCondition2D1N) KRATOS_REGISTER_CONDITION("PointLoadCondition3D1N", mPointLoadCondition3D1N) KRATOS_REGISTER_CONDITION("PointContactCondition2D1N", mPointContactCondition2D1N) diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application.h b/applications/StructuralMechanicsApplication/structural_mechanics_application.h index f2a7b67d0e3b..4350f4289029 100644 --- a/applications/StructuralMechanicsApplication/structural_mechanics_application.h +++ b/applications/StructuralMechanicsApplication/structural_mechanics_application.h @@ -110,6 +110,7 @@ #include "custom_conditions/point_moment_condition_3d.h" #include "custom_conditions/displacement_control_condition.h" #include "custom_conditions/moving_load_condition.h" +#include "custom_conditions/point_load_condition_1d_1n.h" /* Adding the displacement-based SBM condition */ #include "custom_conditions/displacement_shifted_boundary_condition.h" @@ -469,6 +470,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) KratosStructuralMechanicsAppl /* CONDITIONS*/ // Point load + const PointLoadCondition1D1N mPointLoadCondition1D1N; const PointLoadCondition mPointLoadCondition2D1N; const PointLoadCondition mPointLoadCondition3D1N; const PointContactCondition mPointContactCondition2D1N; From 19eb838a98085bfa2869782f65fe95a0e36601af Mon Sep 17 00:00:00 2001 From: matekelemen Date: Wed, 29 Apr 2026 09:29:40 +0200 Subject: [PATCH 12/18] add missing name for Quadrilateral2D4 --- kratos/geometries/quadrilateral_2d_4.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/kratos/geometries/quadrilateral_2d_4.h b/kratos/geometries/quadrilateral_2d_4.h index a9f6bd90b5b4..253e03151775 100644 --- a/kratos/geometries/quadrilateral_2d_4.h +++ b/kratos/geometries/quadrilateral_2d_4.h @@ -663,6 +663,10 @@ template class Quadrilateral2D4 return "2 dimensional quadrilateral with four nodes in 2D space"; } + std::string Name() const override { + return "Quadrilateral2D4"; + } + /** * Print information about this object. * @param rOStream Stream to print into it. From 27b911df776282683a1c85debfaaabee50a74123 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Wed, 29 Apr 2026 09:30:27 +0200 Subject: [PATCH 13/18] bugfixes to pre-tensioning in 2D and 3D --- .../insert_pre_tension_operation.cpp | 421 +++++++----------- 1 file changed, 160 insertions(+), 261 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp index cf702183606c..b71ff67a1ef6 100644 --- a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp +++ b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp @@ -11,6 +11,7 @@ // --- Structural Includes --- #include "custom_operations/insert_pre_tension_operation.hpp" // InsertPreTensionOperation +#include "custom_conditions/point_load_condition_1d_1n.h" // PointLoadCondition1D1N #include "structural_mechanics_application_variables.h" // POINT_LOAD_X // --- Core Includes --- @@ -21,7 +22,6 @@ #include "includes/element.h" // Element #include "input_output/logger.h" // KRATOS_INFO #include "solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.hpp" // LinearMultifreedomConstraint -#include "geometries/point_2d.h" // Point2D // --- STL Includes --- #include // std::array @@ -42,8 +42,15 @@ namespace Kratos { struct InsertPreTensionOperation::Impl { ModelPart* mpModelPart; - double mPretensionValue; + double mPreTensionValue; int mVerbosity; + + // Setting NEIGHBOUR_ELEMENTS on data value containers is fucked, + // so geometry-element adjacency has to be stored in another manner. + std::shared_ptr::IndexType, + std::vector> + > mpAdjacencyMap; }; // struct InsertPreTensionOperation::Impl @@ -54,7 +61,7 @@ InsertPreTensionOperation::InsertPreTensionOperation(Model& rModel, Parameters S KRATOS_TRY mpImpl->mpModelPart = &rModel.GetModelPart(Settings["model_part_name"].GetString()); - mpImpl->mPretensionValue = Settings["magnitude"].GetDouble(); + mpImpl->mPreTensionValue = Settings["magnitude"].GetDouble(); mpImpl->mVerbosity = Settings["verbosity"].GetInt(); KRATOS_CATCH("") @@ -64,6 +71,10 @@ InsertPreTensionOperation::InsertPreTensionOperation(Model& rModel, Parameters S KRATOS_ERROR_IF(mpImpl->mpModelPart->Geometries().empty()) << "InsertPreTensionOperation requires at least one Geometry in the provided ModelPart, " << "but there aren't any in '" << mpImpl->mpModelPart->Name() << "'."; + + mpImpl->mpAdjacencyMap = std::make_shared::IndexType, + std::vector>>(); } @@ -75,117 +86,79 @@ InsertPreTensionOperation::~InsertPreTensionOperation() = default; struct PreTensionSurfacePartition { array_1d normal; - std::vector positive_side, negative_side; + std::vector positive_side, negative_side; }; // struct SurfacePartition -void FindElementsAdjacentToGeometries(ModelPart& rModelPart) { - using NeighborElements = decltype(NEIGHBOUR_ELEMENTS)::Type; - - // Make sure every node on the pre-tension surface is in the model part. - KRATOS_TRY - ModelPart::NodesContainerType nodes; - for (Geometry& r_geometry : rModelPart.Geometries()) { - nodes.insert(r_geometry.begin(), r_geometry.end()); - } // for r_geometry in rModelPart.Conditions - rModelPart.Nodes().insert(nodes); - KRATOS_CATCH("") - - KRATOS_TRY - // Find elements containing each node in the model part. - FindGlobalNodalElementalNeighboursProcess(rModelPart.GetRootModelPart()).Execute(); - - block_for_each( - rModelPart.Geometries(), - [] (Geometry& r_geometry) -> void { - NeighborElements neighbor_elements; - - // Collect all neighbor elements from the geometry's nodes. - for (Node& r_node : r_geometry) { - NeighborElements& r_containing_elements = r_node.GetValue(NEIGHBOUR_ELEMENTS); - for (unsigned i_element=0u; i_element bool { - return rp_other_element->Id() == rp_element->Id(); - }); - if (it_element == neighbor_elements.ptr_end()) { - neighbor_elements.push_back(rp_element); - } - } - } +template +void FindElementsAdjacentToGeometries( + ModelPart& rModelPart, + std::unordered_map::IndexType,std::vector>& rAdjacencyMap) { + for (const auto& r_geometry : rModelPart.Geometries()) + rAdjacencyMap.emplace(r_geometry.Id(), std::vector {}); - // Filter elements that don't contain all nodes of the geometry. - neighbor_elements.erase(std::remove_if( - neighbor_elements.begin(), - neighbor_elements.end(), - [&r_geometry] (const Element& r_element) -> bool { - for (const Node& r_geometry_node : r_geometry) { - const auto it_node = std::find_if( - r_element.GetGeometry().begin(), - r_element.GetGeometry().end(), - [&r_geometry_node] (const Node& r_element_node) -> bool { - return r_element_node.Id() == r_geometry_node.Id(); + // Find elements containing geometries. + KRATOS_TRY + block_for_each( + rModelPart.Geometries(), + [&rAdjacencyMap] (Geometry& r_geometry) -> void { + std::vector adjacent_elements; + + // Collect all neighbor elements from the geometry's nodes. + for (Node& r_node : r_geometry) { + auto& r_containing_elements = r_node.GetValue(NEIGHBOUR_ELEMENTS); + for (unsigned i_element=0u; i_element bool { + return rp_other_element->Id() == rp_element->Id(); }); - if (it_node == r_element.GetGeometry().end()) - return true; - } - return false; - } - ), neighbor_elements.end()); - - // Set neighbor elements on the condition. - r_geometry.SetValue(NEIGHBOUR_ELEMENTS, neighbor_elements); - }); - KRATOS_CATCH("") + if (it_element == adjacent_elements.end()) { + adjacent_elements.emplace_back(Element::Pointer(rp_element.get())); + } + } // for i_element in r_containing_elements.size() + } // for node in geometry + + // Set neighbor elements on the geometry. + rAdjacencyMap[r_geometry.Id()] = std::move(adjacent_elements); + }); + KRATOS_CATCH("") } -template +template PreTensionSurfacePartition PartitionPreTensionSurface( const ModelPart::GeometryContainerType::iterator itGeometryBegin, - const ModelPart::GeometryContainerType::iterator itGeometryEnd) { + const ModelPart::GeometryContainerType::iterator itGeometryEnd, + const std::unordered_map< + Geometry::IndexType, + std::vector + >& rAdjacencyMap) { PreTensionSurfacePartition surface_partition; + std::fill( + surface_partition.normal.begin(), + surface_partition.normal.end(), + 0.0); // There are three valid possibilities here, separately for split and unsplit cases. // 1) All surfaces are 0-dimensional (corners) // => expecting exactly 1 geometry - // - unsplit case - // => expecting exactly 2 connected elements, both - // must have line geometries - // => the normal is not well-defined, so it is chosen - // as the average of the connected elements' directions - // - split case - // => expecting exactly 1 connected element that must have a line geometry - // => the normal is not well-defined, so it is chosen - // as the connected element's direction + // => all connected elements must have line geometries + // => the normal is not well-defined, so it is chosen + // as the average of the connected elements' directions // 2) All surfaces are 1-dimensional (edges) // => expecting at least 1 geometry - // - unsplit case - // => expecting at least 2 connected elements, all of which - // must be 2-dimensional - // => the normal is not well-defined, so it is chosen as the - // average of the rotated endpoints in 2D, in positive - // direction. This assumes coordinates in z-direction vanish. - // This assumption is checked. - // - split case - // => expecting at least 1 connected element, all of which - // must be 2-dimensional - // => the normal is not well-defined, so it is chosen as the - // average of the rotated endpoints in 2D, in positive - // direction. This assumes coordinates in z-direcion vanish. - // This assumption is checked. + // => all connected elements must be 2-dimensional + // => the normal is not well-defined, so it is chosen as the + // average of the rotated endpoints in 2D, in positive + // direction. This assumes coordinates in z-direction vanish. + // This assumption is checked. // 3) All surfaces are 2-dimensional (surfaces) // => expecting at least 1 geometry // => normals are well-defined - // - unsplit case - // => expecting at least 2 connected elements, all - // of which must be 3-dimensional - // - split case - // => expecting at least 1 connected element, all - // of which must be 3-dimensional + // => all connected elements must be 3-dimensional // // In all unsplit cases, at least one connected element is required on // either side of the average pre-tension plane. @@ -216,44 +189,15 @@ PreTensionSurfacePartition PartitionPreTensionSurface( << "but geometry " << r_geometry.Id() << " on geometry " << r_geometry.Id() << " (" << r_geometry.Name() << ") " << "is " << r_geometry.LocalSpaceDimension() << "-dimensional."; - if constexpr (IsSplit) { - KRATOS_ERROR_IF_NOT(r_geometry.Has(NEIGHBOUR_ELEMENTS)) - << "Expecting geometry " << r_geometry.Id() << " to have exactly " - << "1 neighboring element, but found none."; - } else { - KRATOS_ERROR_IF_NOT(r_geometry.Has(NEIGHBOUR_ELEMENTS)) - << "Expecting geometry " << r_geometry.Id() << " to have exactly " - << "2 neighboring elements, but found none."; - } - - auto& r_neighbor_elements = r_geometry.GetValue(NEIGHBOUR_ELEMENTS); - if (r_neighbor_elements.size() != (IsSplit ? 1 : 2)) { - std::stringstream message; - message << "Expecting geometry " << r_geometry.Id() << " to have exactly " - << (IsSplit ? '1' : '2') - << " neighboring elements, but found " << r_neighbor_elements.size() << ": ["; - for (const Element& r_element : r_neighbor_elements) message << r_element.Id() << " "; - message << "]."; - KRATOS_ERROR << message.str(); - } + const auto it_adjacency = rAdjacencyMap.find(r_geometry.Id()); + KRATOS_ERROR_IF(it_adjacency == rAdjacencyMap.end()) + << "no adjacency map is stored for geometry " << r_geometry.Id(); - // Check whether the 2 connected elements lie on opposite sides of the surface (only in the unsplit case). - const array_1d - first_connected_direction = r_neighbor_elements[0].GetGeometry().Center() - r_geometry.Center(), - second_connected_direction = r_neighbor_elements[IsSplit ? 0 : 1].GetGeometry().Center() - r_geometry.Center(); - const double direction_inner_product = std::inner_product( - first_connected_direction.begin(), - first_connected_direction.end(), - second_connected_direction.begin(), - 0.0); - KRATOS_ERROR_IF_NOT(IsSplit || direction_inner_product < 0.0) - << "Element " << r_neighbor_elements[0].Id() << " and " << r_neighbor_elements[1].Id() << " " - << "do not lie on opposite sides of condition " << r_geometry.Id() << " " - << "defining a pre-tension surface."; + const std::vector& r_adjacent_elements = it_adjacency->second; if constexpr (SurfaceDimension == 0u) { - for (unsigned i_element=0u; i_element @@ -300,21 +244,17 @@ PreTensionSurfacePartition PartitionPreTensionSurface( surface_partition.normal += surface_normal; // Sort connected elements. - for (unsigned i_element=0u; i_element element_direction = r_element.GetGeometry().Center() - surface_center; const double direction_inner_product = std::inner_product( surface_normal.begin(), surface_normal.end(), element_direction.begin(), 0.0); - if (0 < direction_inner_product) surface_partition.negative_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); - else surface_partition.positive_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); - } // for i_element in range(r_neighbor_elements.size()) - - KRATOS_ERROR_IF_NOT(IsSplit || surface_partition.positive_side.size() == surface_partition.negative_side.size()) - << "Elements " << r_neighbor_elements[0].Id() << " and " << r_neighbor_elements[1].Id() << " " - << "lie on the same side of geometry " << r_geometry.Id() << " that defines a pre-tension surface."; + if (0 < direction_inner_product) surface_partition.negative_side.push_back(*(r_adjacent_elements.begin() + i_element)); + else surface_partition.positive_side.push_back(*(r_adjacent_elements.begin() + i_element)); + } // for i_element in range(r_adjacent_elements.size()) } /*else if SurfaceDimension == 1*/ else if constexpr (SurfaceDimension == 2u) { // Define the surface's plane. const array_1d surface_center = r_geometry.Center(); @@ -322,21 +262,17 @@ PreTensionSurfacePartition PartitionPreTensionSurface( surface_partition.normal += surface_normal; // Sort connected elements. - for (unsigned i_element=0u; i_element element_direction = r_element.GetGeometry().Center() - surface_center; const double direction_inner_product = std::inner_product( surface_normal.begin(), surface_normal.end(), element_direction.begin(), 0.0); - if (0 < direction_inner_product) surface_partition.negative_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); - else surface_partition.positive_side.push_back(*(r_neighbor_elements.ptr_begin() + i_element)); - } // for i_element in range(r_neighbor_elements.size()) - - KRATOS_ERROR_IF_NOT(IsSplit || surface_partition.positive_side.size() == surface_partition.negative_side.size()) - << "Elements " << r_neighbor_elements[0].Id() << " and " << r_neighbor_elements[1].Id() << " " - << "lie on the same side of condition " << r_geometry.Id() << " that defines a pre-tension surface."; + if (0 < direction_inner_product) surface_partition.negative_side.push_back(*(r_adjacent_elements.begin() + i_element)); + else surface_partition.positive_side.push_back(*(r_adjacent_elements.begin() + i_element)); + } // for i_element in range(r_adjacent_elements.size()) } /*else if SurfaceDimension == 2*/ } // for r_condition in r_model_part.Conditions @@ -376,6 +312,7 @@ class PlaneDisplacementConstraint : public MultifreedomConstraint { DofPointerVectorType&& rDofs, const std::vector& rConstraintLabels, std::shared_ptr pSharedSurfaceNormal, + std::shared_ptr::IndexType,std::vector>> pAdjacencyMap, int Verbosity) : MultifreedomConstraint( Id, @@ -383,6 +320,7 @@ class PlaneDisplacementConstraint : public MultifreedomConstraint { rConstraintLabels), mVerbosity(Verbosity), mSurfaceGeometries(rSurfaceGeometries), + mpAdjacencyMap(pAdjacencyMap), mpSharedSurfaceNormal(pSharedSurfaceNormal) { KRATOS_ERROR_IF(rSurfaceGeometries.empty()); } @@ -445,21 +383,24 @@ class PlaneDisplacementConstraint : public MultifreedomConstraint { PreTensionSurfacePartition surface_partition; switch (mSurfaceGeometries.front().LocalSpaceDimension()) { case 0: { - surface_partition = PartitionPreTensionSurface<0,true>( + surface_partition = PartitionPreTensionSurface<0>( mSurfaceGeometries.begin(), - mSurfaceGeometries.end()); + mSurfaceGeometries.end(), + *mpAdjacencyMap); break; } case 1: { - surface_partition = PartitionPreTensionSurface<1,true>( + surface_partition = PartitionPreTensionSurface<1>( mSurfaceGeometries.begin(), - mSurfaceGeometries.end()); + mSurfaceGeometries.end(), + *mpAdjacencyMap); break; } case 2: { - surface_partition = PartitionPreTensionSurface<2,true>( + surface_partition = PartitionPreTensionSurface<2>( mSurfaceGeometries.begin(), - mSurfaceGeometries.end()); + mSurfaceGeometries.end(), + *mpAdjacencyMap); break; } default: KRATOS_ERROR << "Invalid pre-tension surface dimension: " << mSurfaceGeometries.front().LocalSpaceDimension() << "."; @@ -480,6 +421,11 @@ class PlaneDisplacementConstraint : public MultifreedomConstraint { ModelPart::GeometryContainerType mSurfaceGeometries; + std::shared_ptr::IndexType, + std::vector> + > mpAdjacencyMap; + /// @details The surface's normal is potentially shared by a set of constraints, /// since each constraint only restricts one pair of nodes. If the surface /// contains more than a single pair, computing its normal should not be @@ -509,6 +455,7 @@ class InPlaneRelativeDisplacementConstraint : public PlaneDisplacementConstraint DofPointerVectorType(this->GetDofs()), constraint_labels, mpSharedSurfaceNormal, + mpAdjacencyMap, mVerbosity)); } @@ -572,11 +519,12 @@ void InsertPreTensionOperation::Execute() { ModelPart& r_model_part = *mpImpl->mpModelPart; // Find elements connected to the pre-tension surface. - // => each node will have a NEIGHBOR_ELEMENTS variable + // => each node will have a NEIGHBOUR_ELEMENTS variable // in its non-historical storage that lists elements // containing them. KRATOS_TRY - FindElementsAdjacentToGeometries(r_model_part); + FindGlobalNodalElementalNeighboursProcess(r_model_part.GetRootModelPart()).Execute(); + FindElementsAdjacentToGeometries(r_model_part, *mpImpl->mpAdjacencyMap); KRATOS_CATCH("") // Sort connected elements into 2 categories, lying on either side of @@ -591,19 +539,22 @@ void InsertPreTensionOperation::Execute() { case 0u: { surface_partition = PartitionPreTensionSurface<0>( r_model_part.Geometries().begin(), - r_model_part.Geometries().end()); + r_model_part.Geometries().end(), + *mpImpl->mpAdjacencyMap); break; } case 1u: { surface_partition = PartitionPreTensionSurface<1>( r_model_part.Geometries().begin(), - r_model_part.Geometries().end()); + r_model_part.Geometries().end(), + *mpImpl->mpAdjacencyMap); break; } case 2u: { surface_partition = PartitionPreTensionSurface<2>( r_model_part.Geometries().begin(), - r_model_part.Geometries().end()); + r_model_part.Geometries().end(), + *mpImpl->mpAdjacencyMap); break; } default: KRATOS_ERROR << "Invalid pre-tension surface dimension: " << pre_tension_surface_dimension << "."; @@ -685,9 +636,32 @@ void InsertPreTensionOperation::Execute() { } // for rp_element in surface_partition.negative_side KRATOS_CATCH("") - // Reset NEIGHBOR_ELEMENTS after duplicating nodes. + // Prune elements that don't have at least one face on the pre-tension surface. + KRATOS_TRY + for (auto& [id_geometry, r_adjacent_elements] : *mpImpl->mpAdjacencyMap) { + const Geometry& r_geometry = r_model_part.GetGeometry(id_geometry); + r_adjacent_elements.erase(std::remove_if( + r_adjacent_elements.begin(), + r_adjacent_elements.end(), + [&r_geometry] (const Element::Pointer& rp_element) -> bool { + for (const Node& r_geometry_node : r_geometry) { + const auto it_node = std::find_if( + rp_element->GetGeometry().begin(), + rp_element->GetGeometry().end(), + [&r_geometry_node] (const Node& r_element_node) -> bool { + return r_element_node.Id() == r_geometry_node.Id(); + }); + if (it_node == rp_element->GetGeometry().end()) + return true; + } // for r_geometry_node in r_geometry + return false; + } + ), r_adjacent_elements.end()); + } // for id_geometry, r_adjacency_elements + KRATOS_CATCH("") + KRATOS_TRY - FindElementsAdjacentToGeometries(r_model_part); + FindGlobalNodalElementalNeighboursProcess(r_model_part.GetRootModelPart()).Execute(); KRATOS_CATCH("") // Collect DoFs from elements on the positive side of the pre-tension surface. @@ -752,6 +726,7 @@ void InsertPreTensionOperation::Execute() { std::move(dofs), constraint_labels, p_protected_surface_normal, + mpImpl->mpAdjacencyMap, mpImpl->mVerbosity)); KRATOS_INFO_IF(this->Info(), 3 <= mpImpl->mVerbosity) @@ -844,6 +819,7 @@ class OutOfPlaneDisplacementConstraint : public PlaneDisplacementConstraint { DofPointerVectorType(this->GetDofs()), constraint_labels, mpSharedSurfaceNormal, + mpAdjacencyMap, mVerbosity)); } @@ -854,18 +830,19 @@ class OutOfPlaneDisplacementConstraint : public PlaneDisplacementConstraint { // DoFs are assumed in a specific layout. // - DISPLACEMENT_X of node 0 - // - DISPLACEMENT_Y of node 0 - // ... - // - DISPLACEMENT_Z (depending on the plane's dimension) of node n // - DISPLACEMENT_X of node 0's pair (only if relative) + // - DISPLACEMENT_Y of node 0 // - DISPLACEMENT_Y of node 0's pair (only if relative) // ... + // - DISPLACEMENT_Z (depending on the plane's dimension) of node n // - DISPLACEMENT_Z (depending on the plane's dimension) of node n's pair (only if relative) // - DoF of the control node auto& r_dofs = this->GetDofs(); KRATOS_ERROR_IF(r_dofs.empty()); const std::size_t controlled_dof_count = r_dofs.size() - 1; + KRATOS_ERROR_IF(IsRelative && controlled_dof_count % 2); + const std::size_t positive_side_dof_count = IsRelative ? controlled_dof_count / 2 : controlled_dof_count; @@ -886,10 +863,11 @@ class OutOfPlaneDisplacementConstraint : public PlaneDisplacementConstraint { Matrix constraint_gradient(1, r_dofs.size()); Vector constraint_gap(1), displacements(r_dofs.size()); + constexpr std::size_t dof_stride = IsRelative ? 2ul : 1ul; for (std::size_t i_dof=0ul; i_dofmpAdjacencyMap, mpImpl->mVerbosity)); mpImpl->mpModelPart->AddMasterSlaveConstraint(p_constraint); // Set a dirichlet condition on the control node's DoF. - p_control_dof->GetSolutionStepValue() = (mpImpl->mPretensionValue - p_control_dof->GetSolutionStepValue()) * rDuplicateNodeMap.size(); + p_control_dof->GetSolutionStepValue() = (mpImpl->mPreTensionValue - p_control_dof->GetSolutionStepValue()) * rDuplicateNodeMap.size(); p_control_dof->FixDof(); KRATOS_CATCH("") } @@ -984,93 +963,6 @@ std::string InsertDirichletPreTensionOperation::Info() const { } -class PointLoadCondition1D1N : public Condition { -public: - KRATOS_CLASS_POINTER_DEFINITION(PointLoadCondition1D1N); - - PointLoadCondition1D1N( - Condition::IndexType Id, - Geometry::Pointer pGeometry) - : Condition(Id, pGeometry) - {} - - Condition::Pointer Clone( - Condition::IndexType Id, - const Condition::NodesArrayType& rNodes) const override { - KRATOS_ERROR_IF_NOT(rNodes.size() == 1) - << "PointLoadCondition1D1N::Clone expects 1 node but got " << rNodes.size(); - return Condition::Pointer(new PointLoadCondition1D1N( - Id, - Geometry::Pointer(new Point2D(rNodes)))); - } - - void EquationIdVector( - Condition::EquationIdVectorType& rIndices, - const ProcessInfo&) const override { - rIndices.resize(1); - rIndices[0] = this->GetGeometry()[0].GetDofs()[0]->EquationId(); - } - - void GetDofList( - Condition::DofsVectorType& rDofs, - const ProcessInfo&) const override { - rDofs.resize(1); - rDofs[0] = this->GetGeometry()[0].GetDofs()[0].get(); - } - - void GetValuesVector( - Vector& rOutput, - int Step) const override { - rOutput.resize(1); - rOutput[0] = this->GetGeometry()[0].FastGetSolutionStepValue(DISPLACEMENT_X, Step); - } - - void GetFirstDerivativesVector( - Vector& rOutput, - int Step) const override { - rOutput.resize(1); - rOutput[0] = this->GetGeometry()[0].FastGetSolutionStepValue(VELOCITY_X, Step); - } - - void GetSecondDerivativesVector( - Vector& rOutput, - int Step) const override { - rOutput.resize(1); - rOutput[0] = this->GetGeometry()[0].FastGetSolutionStepValue(ACCELERATION_X, Step); - } - - void CalculateRightHandSide( - Condition::VectorType& rRhs, - const ProcessInfo&) override { - rRhs.resize(1); - rRhs[0] = this->GetValue(POINT_LOAD_X); - } - - void CalculateLocalSystem( - Condition::MatrixType& rLhs, - Condition::VectorType& rRhs, - const ProcessInfo& rProcessInfo) override { - rLhs.resize(1, 1, false); - rLhs(0, 0) = 0.0; - this->CalculateRightHandSide(rRhs, rProcessInfo); - } - - void CalculateMassMatrix( - Condition::MatrixType& rMatrix, - const ProcessInfo&) override { - rMatrix.resize(1, 1, false); - rMatrix(0, 0) = 0.0; - } - - void CalculateDampingMatrix( - Condition::MatrixType& rMatrix, - const ProcessInfo&) override { - rMatrix.resize(1, 1, false); - rMatrix(0, 0) = 0.0; - } -}; // class PointLoadCondition1D1N - - void InsertNeumannPreTensionOperation::InsertControlNodeConstraints( ModelPart& rModelPart, array_1d SurfaceNormal, @@ -1124,6 +1016,9 @@ void InsertNeumannPreTensionOperation::InsertControlNodeConstraints( std::optional>, LockObject>>(); + Properties::Pointer p_properties = mpImpl->mpModelPart->PropertiesArray().empty() + ? Properties::Pointer(new Properties(1)) + : *mpImpl->mpModelPart->PropertiesArray().begin(); // Insert the positive side constraint and condition. { MasterSlaveConstraint::Pointer p_constraint(new OutOfPlaneDisplacementConstraint( @@ -1132,11 +1027,13 @@ void InsertNeumannPreTensionOperation::InsertControlNodeConstraints( std::move(positive_side_dofs), {constraint_id}, p_protected_surface_normal, + mpImpl->mpAdjacencyMap, mpImpl->mVerbosity)); Condition::Pointer p_condition(new PointLoadCondition1D1N( condition_id, - Geometry::Pointer(new Point2D(pPositiveSideControlNode)))); - p_condition->SetValue(POINT_LOAD_X, mpImpl->mPretensionValue); + Geometry::Pointer(new Point2D(pPositiveSideControlNode)), + p_properties)); + p_condition->SetValue(POINT_LOAD_X, mpImpl->mPreTensionValue); KRATOS_INFO_IF(this->Info(), 3 <= mpImpl->mVerbosity) << "insert constraint " << p_constraint->Id() <<' ' @@ -1163,11 +1060,13 @@ void InsertNeumannPreTensionOperation::InsertControlNodeConstraints( std::move(negative_side_dofs), {constraint_id}, p_protected_surface_normal, + mpImpl->mpAdjacencyMap, mpImpl->mVerbosity)); Condition::Pointer p_condition(new PointLoadCondition1D1N( condition_id, - Geometry::Pointer(new Point2D(p_negative_side_control_node)))); - p_condition->SetValue(POINT_LOAD_X, mpImpl->mPretensionValue); + Geometry::Pointer(new Point2D(p_negative_side_control_node)), + p_properties)); + p_condition->SetValue(POINT_LOAD_X, -mpImpl->mPreTensionValue); KRATOS_INFO_IF(this->Info(), 3 <= mpImpl->mVerbosity) << "insert constraint " << p_constraint->Id() <<' ' From 3e0a613175d0a1b392334de95f7c9a18d84c0c70 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Wed, 29 Apr 2026 14:02:40 +0200 Subject: [PATCH 14/18] extend pre-tensioning to support transient loading --- .../insert_pre_tension_operation.cpp | 50 +++++++---- .../insert_pre_tension_operation.hpp | 21 ++++- .../add_custom_operations_to_python.cpp | 2 + .../dirichlet_pre_tension_process.py | 26 +++--- .../neumann_pre_tension_process.py | 26 +++--- .../pre_tension_process_base.py | 86 +++++++++++++++++++ 6 files changed, 158 insertions(+), 53 deletions(-) create mode 100644 applications/StructuralMechanicsApplication/python_scripts/pre_tension_process_base.py diff --git a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp index b71ff67a1ef6..eddedb1aa3be 100644 --- a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp +++ b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp @@ -42,7 +42,6 @@ namespace Kratos { struct InsertPreTensionOperation::Impl { ModelPart* mpModelPart; - double mPreTensionValue; int mVerbosity; // Setting NEIGHBOUR_ELEMENTS on data value containers is fucked, @@ -61,7 +60,6 @@ InsertPreTensionOperation::InsertPreTensionOperation(Model& rModel, Parameters S KRATOS_TRY mpImpl->mpModelPart = &rModel.GetModelPart(Settings["model_part_name"].GetString()); - mpImpl->mPreTensionValue = Settings["magnitude"].GetDouble(); mpImpl->mVerbosity = Settings["verbosity"].GetInt(); KRATOS_CATCH("") @@ -791,7 +789,6 @@ void InsertPreTensionOperation::Execute() { const Parameters InsertPreTensionOperation::GetDefaultParameters() const { return Parameters(R"({ "model_part_name" : "", - "magnitude" : 0.0, "verbosity" : 1 })"); } @@ -906,7 +903,10 @@ void InsertDirichletPreTensionOperation::InsertControlNodeConstraints( array_1d SurfaceNormal, const std::unordered_map rDuplicateNodeMap, Node::Pointer pControlNode, - const std::unordered_set*> rPositiveSideDofs) const { + const std::unordered_set*> rPositiveSideDofs) { + mpControlNode = pControlNode; + mPreTensionSurfaceSize = rDuplicateNodeMap.size(); + KRATOS_TRY // Insert a constraint that ties the average out-of-plane relative displacement // to a prescribed value. @@ -931,7 +931,7 @@ void InsertDirichletPreTensionOperation::InsertControlNodeConstraints( } // for rp_positive_side_node, rp_negative_side_node // Add the control node's DoF to the constraint equation. - Dof* p_control_dof = pControlNode->GetDofs()[0].get(); + Dof* p_control_dof = mpControlNode->GetDofs()[0].get(); dofs.push_back(p_control_dof); // Construct and register the constraint equation. @@ -950,11 +950,18 @@ void InsertDirichletPreTensionOperation::InsertControlNodeConstraints( mpImpl->mpAdjacencyMap, mpImpl->mVerbosity)); mpImpl->mpModelPart->AddMasterSlaveConstraint(p_constraint); +KRATOS_CATCH("") +} - // Set a dirichlet condition on the control node's DoF. - p_control_dof->GetSolutionStepValue() = (mpImpl->mPreTensionValue - p_control_dof->GetSolutionStepValue()) * rDuplicateNodeMap.size(); - p_control_dof->FixDof(); - KRATOS_CATCH("") + +void InsertDirichletPreTensionOperation::Apply(double Magnitude) { + KRATOS_ERROR_IF_NOT(mpControlNode); + KRATOS_ERROR_IF_NOT(mpControlNode->GetDofs().size() == 1); + KRATOS_TRY + Dof& r_control_dof = *mpControlNode->GetDofs()[0].get(); + r_control_dof.GetSolutionStepValue() = (Magnitude - r_control_dof.GetSolutionStepValue()) * mPreTensionSurfaceSize; + r_control_dof.FixDof(); + KRATOS_CATCH("") } @@ -968,7 +975,7 @@ void InsertNeumannPreTensionOperation::InsertControlNodeConstraints( array_1d SurfaceNormal, const std::unordered_map rDuplicateNodeMap, Node::Pointer pPositiveSideControlNode, - const std::unordered_set*> rPositiveSideDofs) const { + const std::unordered_set*> rPositiveSideDofs) { KRATOS_TRY const std::array*,3> all_displacement_components { &DISPLACEMENT_X, @@ -1029,11 +1036,10 @@ void InsertNeumannPreTensionOperation::InsertControlNodeConstraints( p_protected_surface_normal, mpImpl->mpAdjacencyMap, mpImpl->mVerbosity)); - Condition::Pointer p_condition(new PointLoadCondition1D1N( + mpPositiveSideLoad = Condition::Pointer(new PointLoadCondition1D1N( condition_id, Geometry::Pointer(new Point2D(pPositiveSideControlNode)), p_properties)); - p_condition->SetValue(POINT_LOAD_X, mpImpl->mPreTensionValue); KRATOS_INFO_IF(this->Info(), 3 <= mpImpl->mVerbosity) << "insert constraint " << p_constraint->Id() <<' ' @@ -1043,10 +1049,10 @@ void InsertNeumannPreTensionOperation::InsertControlNodeConstraints( mpImpl->mpModelPart->AddMasterSlaveConstraint(p_constraint); KRATOS_INFO_IF(this->Info(), 3 <= mpImpl->mVerbosity) - << "insert condition " << p_condition->Id() <<' ' + << "insert condition " << mpPositiveSideLoad->Id() <<' ' << "loading " << pPositiveSideControlNode->GetDofs()[0]->GetVariable().Name() << ' ' << "of node " << pPositiveSideControlNode->Id() << "\n"; - mpImpl->mpModelPart->AddCondition(p_condition); + mpImpl->mpModelPart->AddCondition(mpPositiveSideLoad); ++constraint_id; ++condition_id; @@ -1062,11 +1068,10 @@ void InsertNeumannPreTensionOperation::InsertControlNodeConstraints( p_protected_surface_normal, mpImpl->mpAdjacencyMap, mpImpl->mVerbosity)); - Condition::Pointer p_condition(new PointLoadCondition1D1N( + mpNegativeSideLoad = Condition::Pointer(new PointLoadCondition1D1N( condition_id, Geometry::Pointer(new Point2D(p_negative_side_control_node)), p_properties)); - p_condition->SetValue(POINT_LOAD_X, -mpImpl->mPreTensionValue); KRATOS_INFO_IF(this->Info(), 3 <= mpImpl->mVerbosity) << "insert constraint " << p_constraint->Id() <<' ' @@ -1076,10 +1081,10 @@ void InsertNeumannPreTensionOperation::InsertControlNodeConstraints( mpImpl->mpModelPart->AddMasterSlaveConstraint(p_constraint); KRATOS_INFO_IF(this->Info(), 3 <= mpImpl->mVerbosity) - << "insert condition " << p_condition->Id() <<' ' + << "insert condition " << mpNegativeSideLoad->Id() <<' ' << "loading " << p_negative_side_control_node->GetDofs()[0]->GetVariable().Name() << ' ' << "of node " << p_negative_side_control_node->Id() << "\n"; - mpImpl->mpModelPart->AddCondition(p_condition); + mpImpl->mpModelPart->AddCondition(mpNegativeSideLoad); ++constraint_id; ++condition_id; @@ -1088,6 +1093,15 @@ void InsertNeumannPreTensionOperation::InsertControlNodeConstraints( } +void InsertNeumannPreTensionOperation::Apply(double Magnitude) { + KRATOS_ERROR_IF_NOT(mpPositiveSideLoad && mpNegativeSideLoad); + KRATOS_TRY + mpPositiveSideLoad->SetValue(POINT_LOAD_X, Magnitude); + mpNegativeSideLoad->SetValue(POINT_LOAD_X, -Magnitude); + KRATOS_CATCH("") +} + + std::string InsertNeumannPreTensionOperation::Info() const { return "InsertNeumannPreTensionOperation"; } diff --git a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp index ac4ecfa6f5b2..3af177c84d97 100644 --- a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp +++ b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp @@ -23,6 +23,7 @@ namespace Kratos { +/// @brief Operation class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertPreTensionOperation : public Operation { public: KRATOS_CLASS_POINTER_DEFINITION(InsertPreTensionOperation); @@ -35,6 +36,8 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertPreTensionOperation : p void Execute() override; + virtual void Apply(double Magnitude) = 0; + const Parameters GetDefaultParameters() const override; virtual std::string Info() const; @@ -45,7 +48,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertPreTensionOperation : p array_1d SurfaceNormal, const std::unordered_map rDuplicateNodeMap, Node::Pointer pControlNode, - const std::unordered_set*> rPositiveSideDofs) const = 0; + const std::unordered_set*> rPositiveSideDofs) = 0; struct Impl; std::unique_ptr mpImpl; @@ -63,6 +66,8 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertDirichletPreTensionOper using InsertPreTensionOperation::InsertPreTensionOperation; + void Apply(double Magnitude) override; + std::string Info() const override; protected: @@ -73,7 +78,12 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertDirichletPreTensionOper array_1d SurfaceNormal, const std::unordered_map rDuplicateNodeMap, Node::Pointer pControlNode, - const std::unordered_set*> rPositiveSideDofs) const override; + const std::unordered_set*> rPositiveSideDofs) override; + +private: + Node::Pointer mpControlNode; + + std::size_t mPreTensionSurfaceSize; }; // class InsertDirichletPreTensionOperation @@ -83,6 +93,8 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertNeumannPreTensionOperat using InsertPreTensionOperation::InsertPreTensionOperation; + void Apply(double Magnitude) override; + std::string Info() const override; protected: @@ -94,7 +106,10 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertNeumannPreTensionOperat array_1d SurfaceNormal, const std::unordered_map rDuplicateNodeMap, Node::Pointer pControlNode, - const std::unordered_set*> rPositiveSideDofs) const override; + const std::unordered_set*> rPositiveSideDofs) override; + +private: + Condition::Pointer mpPositiveSideLoad, mpNegativeSideLoad; }; // class InsertNeumannPreTensionOperation diff --git a/applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.cpp b/applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.cpp index da7bedfe5c3e..dd15d39a9377 100644 --- a/applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.cpp +++ b/applications/StructuralMechanicsApplication/custom_python/add_custom_operations_to_python.cpp @@ -14,6 +14,7 @@ void AddCustomOperationsToPython(pybind11::module& rModule) { rModule, "InsertDirichletPreTensionOperation") .def(pybind11::init()) + .def("Apply", &InsertDirichletPreTensionOperation::Apply, pybind11::arg("Magnitude")) ; pybind11::class_< InsertNeumannPreTensionOperation, @@ -22,6 +23,7 @@ void AddCustomOperationsToPython(pybind11::module& rModule) { rModule, "InsertNeumannPreTensionOperation") .def(pybind11::init()) + .def("Apply", &InsertNeumannPreTensionOperation::Apply, pybind11::arg("Magnitude")) ; } diff --git a/applications/StructuralMechanicsApplication/python_scripts/dirichlet_pre_tension_process.py b/applications/StructuralMechanicsApplication/python_scripts/dirichlet_pre_tension_process.py index a865e40cd920..28ade66f40e7 100644 --- a/applications/StructuralMechanicsApplication/python_scripts/dirichlet_pre_tension_process.py +++ b/applications/StructuralMechanicsApplication/python_scripts/dirichlet_pre_tension_process.py @@ -3,28 +3,22 @@ # --- Structural Imports --- from KratosMultiphysics import StructuralMechanicsApplication as StructuralMechanics +from KratosMultiphysics.StructuralMechanicsApplication.pre_tension_process_base import PreTensionProcessBase, InsertDirichletPreTensionOperation -class DirichletPretensionProcess(KratosMultiphysics.Process): +class DirichletPretensionProcess(PreTensionProcessBase): def __init__(self, model: KratosMultiphysics.Model, parameters: KratosMultiphysics.Parameters) -> None: - super().__init__() - parameters.ValidateAndAssignDefaults(self.GetDefaultParameters()) - self.__insert_pretension_operation = StructuralMechanics.InsertDirichletPreTensionOperation( - model, - parameters) - - def ExecuteBeforeSolutionLoop(self) -> None: - self.__insert_pretension_operation.Execute() - - def GetDefaultParameters(self) -> KratosMultiphysics.Parameters: - return KratosMultiphysics.Parameters(R"""{ - "model_part_name" : "", - "magnitude" : 0.0, - "verbosity" : 1 - }""") + super().__init__(model, parameters) + + @classmethod + def _MakeOperation( + cls, + model: KratosMultiphysics.Model, + parameters: KratosMultiphysics.Parameters) -> InsertDirichletPreTensionOperation: + return InsertDirichletPreTensionOperation(model, parameters) def Factory(parameters: KratosMultiphysics.Parameters, diff --git a/applications/StructuralMechanicsApplication/python_scripts/neumann_pre_tension_process.py b/applications/StructuralMechanicsApplication/python_scripts/neumann_pre_tension_process.py index da2c16ba8f35..4b46419680b0 100644 --- a/applications/StructuralMechanicsApplication/python_scripts/neumann_pre_tension_process.py +++ b/applications/StructuralMechanicsApplication/python_scripts/neumann_pre_tension_process.py @@ -3,32 +3,26 @@ # --- Structural Imports --- from KratosMultiphysics import StructuralMechanicsApplication as StructuralMechanics +from KratosMultiphysics.StructuralMechanicsApplication.pre_tension_process_base import PreTensionProcessBase, InsertNeumannPreTensionOperation -class NeumannPretensionProcess(KratosMultiphysics.Process): +class NeumannPreTensionProcess(PreTensionProcessBase): def __init__(self, model: KratosMultiphysics.Model, parameters: KratosMultiphysics.Parameters) -> None: - super().__init__() - parameters.ValidateAndAssignDefaults(self.GetDefaultParameters()) - self.__insert_pre_tension_operation = StructuralMechanics.InsertNeumannPreTensionOperation( - model, - parameters) + super().__init__(model, parameters) - def ExecuteBeforeSolutionLoop(self) -> None: - self.__insert_pre_tension_operation.Execute() - - def GetDefaultParameters(self) -> KratosMultiphysics.Parameters: - return KratosMultiphysics.Parameters(R"""{ - "model_part_name" : "", - "magnitude" : 0.0, - "verbosity" : 1 - }""") + @classmethod + def _MakeOperation( + cls, + model: KratosMultiphysics.Model, + parameters: KratosMultiphysics.Parameters) -> InsertNeumannPreTensionOperation: + return InsertNeumannPreTensionOperation(model, parameters) def Factory(parameters: KratosMultiphysics.Parameters, model: KratosMultiphysics.Model) -> KratosMultiphysics.Process: - return NeumannPretensionProcess( + return NeumannPreTensionProcess( model, parameters["Parameters"]) diff --git a/applications/StructuralMechanicsApplication/python_scripts/pre_tension_process_base.py b/applications/StructuralMechanicsApplication/python_scripts/pre_tension_process_base.py new file mode 100644 index 000000000000..e6850820d8af --- /dev/null +++ b/applications/StructuralMechanicsApplication/python_scripts/pre_tension_process_base.py @@ -0,0 +1,86 @@ +# --- Core Imports --- +import KratosMultiphysics + +# --- Structural Imports --- +from KratosMultiphysics.StructuralMechanicsApplication import InsertDirichletPreTensionOperation, InsertNeumannPreTensionOperation + +# --- STD Imports --- +from typing import Union +import abc + + +class __PreTensionProcessMeta(type(KratosMultiphysics.Process), type(abc.ABC)): + pass + + +class PreTensionProcessBase(KratosMultiphysics.Process, abc.ABC, metaclass = __PreTensionProcessMeta): + """ @brief Common base class of @ref NeumannPreTensionProcess and @ref DirichletPreTensionProcess. + @see @ref InsertNeumannPreTensionOperation + @see @ref InsertDirichletPreTensionProcess + """ + + def __init__( + self, + model: KratosMultiphysics.Model, + parameters: KratosMultiphysics.Parameters) -> None: + super().__init__() + parameters.AddMissingParameters(self.GetDefaultParameters()) + + # Strip the "magnitude" setting because it has to be handled separately, + # since it can either be a numeric value or a string. + magnitude_setting: KratosMultiphysics.Parameters = parameters["magnitude"].Clone() + parameters.RemoveValue("magnitude") + + self.__magnitude: Union[float,KratosMultiphysics.GenericFunctionUtility] + if magnitude_setting.IsDouble(): + self.__magnitude = magnitude_setting.GetDouble() + elif magnitude_setting.IsString(): + self.__magnitude = KratosMultiphysics.GenericFunctionUtility(magnitude_setting.GetString()) + if self.__magnitude.DependsOnSpace(): + raise ValueError(f"\"magnitude\" cannot depend on spatial variables, only (pseudo-) time") + else: + raise ValueError(f"expecting either a numeric value or a string defining a function for \"magnitude\", but got {magnitude_setting}") + + self.__model_part: KratosMultiphysics.ModelPart = model.GetModelPart(parameters["model_part_name"].GetString()) + self.__operation: Union[InsertDirichletPreTensionOperation,InsertNeumannPreTensionOperation] = self._MakeOperation(model, parameters) + + + def ExecuteBeforeSolutionLoop(self) -> None: + self.__operation.Execute() + + + def ExecuteInitializeSolutionStep(self) -> None: + self.__operation.Apply(self.__GetMagnitude()) + + + def GetDefaultParameters(self) -> KratosMultiphysics.Parameters: + return KratosMultiphysics.Parameters(R"""{ + "model_part_name" : "", + "magnitude" : 0.0, + "verbosity" : 1 + }""") + + + @classmethod + @abc.abstractmethod + def _MakeOperation( + cls, + model: KratosMultiphysics.Model, + parameters: KratosMultiphysics.Parameters) -> Union[InsertDirichletPreTensionOperation,InsertNeumannPreTensionOperation]: + pass + + + def __GetMagnitude(self) -> float: + if isinstance(self.__magnitude, float): + return self.__magnitude + elif isinstance(self.__magnitude, KratosMultiphysics.GenericFunctionUtility): + return self.__magnitude.CallFunction( + 0.0, + 0.0, + 0.0, + self.__model_part.ProcessInfo[KratosMultiphysics.TIME], + 0.0, + 0.0, + 0.0) + else: + raise RuntimeError(f"unhandled magnitude type {type(self.__magnitude)}") From c0c6de8d5e30a5610670aa803a3d567f61d36942 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Wed, 29 Apr 2026 17:11:32 +0200 Subject: [PATCH 15/18] add docs for pre-tensioning --- .../insert_pre_tension_operation.hpp | 123 +++++++++++++++++- .../dirichlet_pre_tension_process.py | 9 ++ .../neumann_pre_tension_process.py | 9 ++ .../pre_tension_process_base.py | 4 +- 4 files changed, 143 insertions(+), 2 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp index 3af177c84d97..7de4226fa857 100644 --- a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp +++ b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp @@ -23,7 +23,115 @@ namespace Kratos { -/// @brief Operation +/// @defgroup pre_tensioning Pre-Tensioning +/// @brief Apply pre-tensioning to 1, 2, or 3D structural parts. +/// @details Pre-tensioning approximately models connector parts (e.g.: bolts, screws, etc.) subject to +/// initial loading within an assembly. It is meant to be a simplified alternative to high-fidelity +/// analyses of such assemblies, with less focus around the connector part but similar effects on +/// the rest of the structure. +/// +/// @subsection pre_tensioning_surface Pre-Tensioning Surface +/// The pre-tensioning system acts on a surface defined by a set of @ref Geometry "geometries" +/// within a @ref ModelPart "model part". The following requirements apply to this set of geometries: +/// - they @b must be defined on element boundaries (faces of polyhedra, edges of polygons, vertices of +/// lines), +/// - they @b must partition the connector part's elements into 2 distinct sets, +/// - they @b should lie on the same plane as much as possible, and +/// - they @b must be stored in a separate sub model part of the structure that only contains +/// these geometries. +/// +/// A pre-tensioning plane is computed by averaging the pre-tensioning surface. +/// The normal of this plane will be the basis of displacement constraints introduced later in this +/// process. +/// +/// @subsection pre_tensioning_node_duplication Node Duplication +/// The pre-tensioning surface partitions adjacent elements into two distinct groups: +/// the @a positive side and the @a negative side. Nodes lying on the surface are duplicated, and +/// negative side elements' nodes are replaced with them, while elements on the positive side +/// are left unchanged. This effectively means that the connector part gets cut in half along +/// the pre-tensioning surface. Duplicated nodes are inserted into the sub model part +/// containing the geometries that define the surface. +/// +/// @subsection pre_tensioning_in_plane_constraints In-Plane Constraints +/// After duplicating nodes, some @ref MultifreedomConstraint "constraints" must be inserted +/// that restrict displacements (and rotations). First of all, if the nodes have DoFs for rotations, +/// the rotations of the duplicates must match their original counterparts. +/// @f[ +/// \phi_o^i - \phi_d^i = 0 \quad \forall \quad 1 \leq i \leq m +/// @f] +/// where +/// - @f$\phi_o^i@f$ is the rotation of the original node @f$i@f$, +/// - @f$\phi_d^i@f$ is the rotation of node @f$i@f$'s duplicate, and +/// - @f$m@f$ is the number of original nodes on the pre-tensioning surface. +/// +/// Also, the relative +/// in-plane displacement component of original-duplicate node pairs must vanish. +/// @f[ +/// u_o^i - u_d^i - \left< u_o^i - u_d^i, n \right> n = 0 \quad \forall \quad 1 \leq i \leq m +/// @f] +/// where +/// - @f$u_o^i@f$ is the displacement of the original node @f$i@f$, +/// - @f$u_d^i@f$ is the displacement of node @f$i@f$'s duplicate, +/// - @f$n@f$ is the unit normal of the pre-tensioning plane, and +/// - @f$\left< \cdot , \cdot \right>@f$ denotes an inner product. +/// +/// @subsection pre_tensioning_out_of_plane_constraints Out-of-Plane Constraints +/// To apply a single pre-tensioning value (load or prescribed displacement) to the surface, +/// the out-of-plane displacement components on each set of nodes (original or duplicated) are +/// averaged and tied to newly inserted virtual DoFs. +/// @f[ +/// \begin{align} +/// u_v^o - \frac{1}{m} \sum_{i=1}^m{ \left< u_o^i, n \right> } &= 0 \\ +/// u_v^d - \frac{1}{m} \sum_{i=1}^m{ \left< u_d^i, n \right> } &= 0 +/// \end{align} +/// @f] +/// where +/// - @f$u_v^o@f$ is the newly inserted virtual DoF belonging to the original set of nodes, and +/// - @f$u_v^d@f$ is the newly inserted virtual DoF belonging to the duplicate set of nodes. +/// +/// Once the virtual degrees-of-freedom @f$u_v^o@f$ and @f$u_v^d@f$ are defined, they can either +/// be loaded (Neumann-type) or fixed (Dirichlet-type). +/// +/// @subsection pre_tensioning_neumann_type Neumann-Type Pre-Tensioning +/// Given a pre-tensioning force @f$f@f$, fix the reaction of @f$u_v^o@f$ to @f$f@f$ and +/// the reaction of @f$u_v^d@f$ to @f$-f@f$. +/// @f[ +/// \begin{align} +/// r_v^o - f &= 0 \\ +/// r_v^d + f &= 0 +/// \end{align} +/// @f] +/// where +/// - @f$r_v^0@f$ is the reaction of @f$u_v^o@f$, and +/// - @f$r_v^d@f$ is the reaction of @f$u_v^d@f$. +/// +/// @see @ref NeumannPreTensionProcess +/// @see @ref InsertNeumannPreTensionOperation +/// +/// @subsection pre_tensioning_dirichlet_type Dirichlet-Type Pre-Tensioning +/// Fix the relative average out-of-plane displacement to a prescribed value @f$\alpha@f$. +/// @f[ +/// u_v^o - u_v^d - \alpha = 0 +/// @f] +/// +/// @see @ref DirichletPreTensionProcess +/// @see @ref InsertDirichletPreTensionOperation +/// +/// @subsection pre_tensioning_implementation_notes Implementation Notes +/// Since in-plane displacement constraints reuse DoFs and are dense, they cannot be imposed +/// via master-slave elimination unless DoFs are rotated to line up with the +/// pre-tensioning plane. Since there is currently no robust way of doing this in +/// Kratos, constraints must be imposed via the method of augmented lagrange multipliers, +/// which is implemented in @ref p_multigrid "PMultigridBuilderAndSolver". +/// @note Any analysis involving pre-tensioning must use @ref PMultigridBuilderAndSolver. + + + +/// @brief Base class for @ref InsertDirichletPreTensionOperation and @ref InsertNeumannPreTensionOperation. +/// @see @ref InsertDirichletPreTensionOperation +/// @see @ref InsertNeumannPreTensionOperation +/// @see @ref pre_tensioning "Pre-Tensioning" +/// @ingroup pre_tensioning class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertPreTensionOperation : public Operation { public: KRATOS_CLASS_POINTER_DEFINITION(InsertPreTensionOperation); @@ -43,6 +151,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertPreTensionOperation : p virtual std::string Info() const; protected: + /// @brief Handle the average out-of-plane displacement. virtual void InsertControlNodeConstraints( ModelPart& rModelPart, array_1d SurfaceNormal, @@ -60,6 +169,12 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertPreTensionOperation : p }; // class InsertPreTensionOperation +/// @brief Pre-tensioning defined by a surface and a prescribed displacement. +/// @details Cut the mesh along the provided surface and apply fix the average +/// out-of-plane displacement to the provided value while forbidding relative +/// in-plane displacements. +/// @see @ref pre_tensioning "Pre-Tensioning" +/// @ingroup pre_tensioning class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertDirichletPreTensionOperation final : public InsertPreTensionOperation { public: KRATOS_CLASS_POINTER_DEFINITION(InsertDirichletPreTensionOperation); @@ -87,6 +202,12 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertDirichletPreTensionOper }; // class InsertDirichletPreTensionOperation +/// @brief Pre-tensioning defined by a surface and a force. +/// @details Cut the mesh along the provided surface and apply a force to +/// out-of-plane displacement components while forbidding relative +/// in-plane displacements. +/// @see @ref pre_tensioning "Pre-Tensioning" +/// @ingroup pre_tensioning class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) InsertNeumannPreTensionOperation final : public InsertPreTensionOperation { public: KRATOS_CLASS_POINTER_DEFINITION(InsertNeumannPreTensionOperation); diff --git a/applications/StructuralMechanicsApplication/python_scripts/dirichlet_pre_tension_process.py b/applications/StructuralMechanicsApplication/python_scripts/dirichlet_pre_tension_process.py index 28ade66f40e7..7012770baca5 100644 --- a/applications/StructuralMechanicsApplication/python_scripts/dirichlet_pre_tension_process.py +++ b/applications/StructuralMechanicsApplication/python_scripts/dirichlet_pre_tension_process.py @@ -7,6 +7,15 @@ class DirichletPretensionProcess(PreTensionProcessBase): + """ + @brief Pre-tensioning defined by a surface and a prescribed displacement. + @details Cut the mesh along the provided surface and fix the average + out-of-plane displacement to the provided value while forbidding relative + in-plane displacements. + @see @ref pre_tensioning "Pre-Tensioning" + @classname DirichletPreTensionProcess + @ingroup pre_tensioning + """ def __init__(self, model: KratosMultiphysics.Model, diff --git a/applications/StructuralMechanicsApplication/python_scripts/neumann_pre_tension_process.py b/applications/StructuralMechanicsApplication/python_scripts/neumann_pre_tension_process.py index 4b46419680b0..f598bfae8d2c 100644 --- a/applications/StructuralMechanicsApplication/python_scripts/neumann_pre_tension_process.py +++ b/applications/StructuralMechanicsApplication/python_scripts/neumann_pre_tension_process.py @@ -7,6 +7,15 @@ class NeumannPreTensionProcess(PreTensionProcessBase): + """ + @brief Pre-tensioning defined by a surface and a force. + @details Cut the mesh along the provided surface and apply a force to + out-of-plane displacement components while forbidding relative + in-plane displacements. + @see @ref pre_tensioning "Pre-Tensioning" + @classname NeumannPreTensionProcess + @ingroup pre_tensioning + """ def __init__(self, model: KratosMultiphysics.Model, diff --git a/applications/StructuralMechanicsApplication/python_scripts/pre_tension_process_base.py b/applications/StructuralMechanicsApplication/python_scripts/pre_tension_process_base.py index e6850820d8af..d3a658bc1bc2 100644 --- a/applications/StructuralMechanicsApplication/python_scripts/pre_tension_process_base.py +++ b/applications/StructuralMechanicsApplication/python_scripts/pre_tension_process_base.py @@ -10,13 +10,15 @@ class __PreTensionProcessMeta(type(KratosMultiphysics.Process), type(abc.ABC)): - pass + pass class PreTensionProcessBase(KratosMultiphysics.Process, abc.ABC, metaclass = __PreTensionProcessMeta): """ @brief Common base class of @ref NeumannPreTensionProcess and @ref DirichletPreTensionProcess. @see @ref InsertNeumannPreTensionOperation @see @ref InsertDirichletPreTensionProcess + @see @ref pre_tensioning "Pre-Tensioning" + @ingroup pre_tensioning """ def __init__( From 0792664a03d96ffbd139671e7b07acd9f36aba23 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Wed, 29 Apr 2026 18:42:38 +0200 Subject: [PATCH 16/18] fix bugs in dirichlet pre-tensioning --- .../insert_pre_tension_operation.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp index eddedb1aa3be..fde0a5e1943e 100644 --- a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp +++ b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.cpp @@ -218,6 +218,7 @@ PreTensionSurfacePartition PartitionPreTensionSurface( surface_partition.normal = opposite_node_position - surface_node_position; } else { surface_partition.negative_side.push_back(*(r_adjacent_elements.begin() + i_element)); + surface_partition.normal = surface_node_position - opposite_node_position; } } // for i_element in range(r_adjacent_elements.size()) } /*if SurfaceDimension == 0*/ else if constexpr (SurfaceDimension == 1u) { @@ -848,9 +849,9 @@ class OutOfPlaneDisplacementConstraint : public PlaneDisplacementConstraint { std::size_t dimension_count = 0ul; { std::unordered_set keys; - for (std::size_t i_dof=0ul; i_dofGetVariable().Key()); - } // for i_dof in positive_side_dof_count + } // for i_dof in r_dofs.size()-1 dimension_count = keys.size(); KRATOS_ERROR_IF_NOT(0 < dimension_count && dimension_count < 4) << "unexpected surface dimension " << dimension_count; @@ -862,9 +863,9 @@ class OutOfPlaneDisplacementConstraint : public PlaneDisplacementConstraint { constexpr std::size_t dof_stride = IsRelative ? 2ul : 1ul; for (std::size_t i_dof=0ul; i_dofGetDofs().size() == 1); KRATOS_TRY Dof& r_control_dof = *mpControlNode->GetDofs()[0].get(); - r_control_dof.GetSolutionStepValue() = (Magnitude - r_control_dof.GetSolutionStepValue()) * mPreTensionSurfaceSize; + r_control_dof.GetSolutionStepValue() = 0.5 * mPreTensionSurfaceSize * Magnitude - r_control_dof.GetSolutionStepValue(); r_control_dof.FixDof(); KRATOS_CATCH("") } From 26b1d8dcbcc3760a5f17ba80aeeb19eb68370312 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Wed, 29 Apr 2026 18:43:00 +0200 Subject: [PATCH 17/18] update pre-tensioning test --- .../tests/constraints/pre_tensioned_beam.json | 54 +++++------- .../tests/test_pre_tension.py | 85 ++++++++++++++++++- 2 files changed, 104 insertions(+), 35 deletions(-) diff --git a/applications/StructuralMechanicsApplication/tests/constraints/pre_tensioned_beam.json b/applications/StructuralMechanicsApplication/tests/constraints/pre_tensioned_beam.json index 55e4afdb8c6d..c0eff9ef2288 100644 --- a/applications/StructuralMechanicsApplication/tests/constraints/pre_tensioned_beam.json +++ b/applications/StructuralMechanicsApplication/tests/constraints/pre_tensioned_beam.json @@ -23,7 +23,8 @@ "displacement_absolute_tolerance": 1e-9, "rotation_dofs" : true, - "move_mesh_flag" : true, + "move_mesh_flag" : false, + "compute_reactions" : true, "model_import_settings": { "input_type": "mdpa", @@ -72,15 +73,6 @@ "constrained": [false, false, true], "value": [null, null, 0] } - }, - { - "python_module" : "neumann_pre_tension_process", - "kratos_module" : "KratosMultiphysics.StructuralMechanicsApplication", - "Parameters" : { - "model_part_name" : "root.pre_tension_surface", - "magnitude" : 2e7, - "verbosity" : 1 - } } ], "loads_process_list": [ @@ -100,25 +92,25 @@ } ] } - //,"output_processes": { - //"vtk_output": [ - // { - // "Parameters": { - // "folder_name": "vtk", - // "model_part_name": "root", - // "nodal_data_value_variables": [], - // "nodal_solution_step_data_variables": [ - // "DISPLACEMENT", - // "ROTATION", - // "REACTION" - // ], - // "output_sub_model_parts": false, - // "save_output_files_in_folder": true - // }, - // "kratos_module": "KratosMultiphysics", - // "process_name": "VTKOutputProcess", - // "python_module": "vtk_output_process" - // } - //] - //} + ,"output_processes": { + "vtk_output": [ + { + "Parameters": { + "folder_name": "vtk", + "model_part_name": "root", + "nodal_data_value_variables": [], + "nodal_solution_step_data_variables": [ + "DISPLACEMENT", + "ROTATION", + "REACTION" + ], + "output_sub_model_parts": false, + "save_output_files_in_folder": true + }, + "kratos_module": "KratosMultiphysics", + "process_name": "VTKOutputProcess", + "python_module": "vtk_output_process" + } + ] + } } diff --git a/applications/StructuralMechanicsApplication/tests/test_pre_tension.py b/applications/StructuralMechanicsApplication/tests/test_pre_tension.py index 31e7232e2db7..de008c42780b 100644 --- a/applications/StructuralMechanicsApplication/tests/test_pre_tension.py +++ b/applications/StructuralMechanicsApplication/tests/test_pre_tension.py @@ -1,3 +1,6 @@ +# --- External Imports --- +import numpy + # --- Kratos Imports --- import KratosMultiphysics import KratosMultiphysics.KratosUnittest @@ -6,19 +9,93 @@ # --- STD Imports --- import pathlib +import json class TestPretension(KratosMultiphysics.KratosUnittest.TestCase): - def test_pretensioned_beam(self) -> None: + def test_NeumannPreTensionedBeam(self) -> None: + with WorkFolderScope("constraints", pathlib.Path(__file__).absolute()): + # Load the configuration. + with open("pre_tensioned_beam.json", "r") as file: + configuration = json.loads(file.read()) + + # Insert pre-tensioning. + pre_tensioning_magnitude: float = 2e7 + configuration["processes"]["constraints_process_list"].append({ + "python_module" : "neumann_pre_tension_process", + "kratos_module" : "KratosMultiphysics.StructuralMechanicsApplication", + "Parameters" : { + "model_part_name" : "root.pre_tension_surface", + "magnitude" : pre_tensioning_magnitude, + "verbosity" : 3 + }}) + + # Construct and run the analysis. + parameters = KratosMultiphysics.Parameters(json.dumps(configuration)) + model = KratosMultiphysics.Model() + analysis = StructuralMechanicsAnalysis(model, parameters) + analysis.Run() + + # Check the reaction on the virtual nodes. + model_part: KratosMultiphysics.ModelPart = model.GetModelPart("root").GetSubModelPart("pre_tension_surface") + positive_side_virtual_node: KratosMultiphysics.Node = model_part.GetNode(8) + negative_side_virtual_node: KratosMultiphysics.Node = model_part.GetNode(10) + + self.assertAlmostEqual( + positive_side_virtual_node.GetSolutionStepValue(KratosMultiphysics.REACTION_X), + -pre_tensioning_magnitude) + self.assertAlmostEqual( + negative_side_virtual_node.GetSolutionStepValue(KratosMultiphysics.REACTION_X), + pre_tensioning_magnitude) + + + def test_DirichletPreTensionedBeam(self) -> None: with WorkFolderScope("constraints", pathlib.Path(__file__).absolute()): - # Load config - with open("pre_tensioned_beam.json", "r") as project_parameters_file: - parameters = KratosMultiphysics.Parameters(project_parameters_file.read()) + # Load the configuration. + with open("pre_tensioned_beam.json", "r") as file: + configuration = json.loads(file.read()) + + # Insert pre-tensioning. + pre_tensioning_magnitude: float = 1e-1 + configuration["processes"]["constraints_process_list"].append({ + "python_module" : "dirichlet_pre_tension_process", + "kratos_module" : "KratosMultiphysics.StructuralMechanicsApplication", + "Parameters" : { + "model_part_name" : "root.pre_tension_surface", + "magnitude" : pre_tensioning_magnitude, + "verbosity" : 3 + }}) + + # Construct and run the analysis. + parameters = KratosMultiphysics.Parameters(json.dumps(configuration)) model = KratosMultiphysics.Model() analysis = StructuralMechanicsAnalysis(model, parameters) analysis.Run() + # Check the displacements on the duplicate and virtual nodes. + model_part: KratosMultiphysics.ModelPart = model.GetModelPart("root") + original_node: KratosMultiphysics.Node = model_part.GetNode(7) + duplicate_node: KratosMultiphysics.Node = model_part.GetNode(9) + + normal: numpy.ndarray = numpy.array([2.0, 3.0, 0.0]) + normal /= numpy.linalg.norm(normal) + + relative_displacement: numpy.ndarray = numpy.array( + duplicate_node.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT) + - original_node.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT)) + relative_out_of_plane_displacement = numpy.dot(relative_displacement, normal) * normal + relative_in_plane_displacement = relative_displacement - relative_out_of_plane_displacement + + #raise RuntimeError(relative_displacement) + + self.assertLessEqual( + numpy.linalg.norm(relative_in_plane_displacement), + 1e-6) + self.assertAlmostEqual( + numpy.dot(relative_out_of_plane_displacement, normal), + pre_tensioning_magnitude) + if __name__ == "__main__": KratosMultiphysics.KratosUnittest.main() From 370d1d6b8a6c01b95f22c5355157aa76902faf15 Mon Sep 17 00:00:00 2001 From: matekelemen Date: Wed, 29 Apr 2026 23:37:26 +0200 Subject: [PATCH 18/18] remove multiline doxygen equations --- .../insert_pre_tension_operation.hpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp index 7de4226fa857..3f19124422a0 100644 --- a/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp +++ b/applications/StructuralMechanicsApplication/custom_operations/insert_pre_tension_operation.hpp @@ -80,10 +80,10 @@ namespace Kratos { /// the out-of-plane displacement components on each set of nodes (original or duplicated) are /// averaged and tied to newly inserted virtual DoFs. /// @f[ -/// \begin{align} -/// u_v^o - \frac{1}{m} \sum_{i=1}^m{ \left< u_o^i, n \right> } &= 0 \\ -/// u_v^d - \frac{1}{m} \sum_{i=1}^m{ \left< u_d^i, n \right> } &= 0 -/// \end{align} +/// u_v^o - \frac{1}{m} \sum_{i=1}^m{ \left< u_o^i, n \right> } = 0 +/// @f] +/// @f[ +/// u_v^d - \frac{1}{m} \sum_{i=1}^m{ \left< u_d^i, n \right> } = 0 /// @f] /// where /// - @f$u_v^o@f$ is the newly inserted virtual DoF belonging to the original set of nodes, and @@ -96,10 +96,10 @@ namespace Kratos { /// Given a pre-tensioning force @f$f@f$, fix the reaction of @f$u_v^o@f$ to @f$f@f$ and /// the reaction of @f$u_v^d@f$ to @f$-f@f$. /// @f[ -/// \begin{align} -/// r_v^o - f &= 0 \\ -/// r_v^d + f &= 0 -/// \end{align} +/// r_v^o - f = 0 +/// @f] +/// @f[ +/// r_v^d + f = 0 /// @f] /// where /// - @f$r_v^0@f$ is the reaction of @f$u_v^o@f$, and