diff --git a/applications/CoSimulationApplication/custom_processes/data_transfer_3D_1D_process.cpp b/applications/CoSimulationApplication/custom_processes/data_transfer_3D_1D_process.cpp new file mode 100644 index 000000000000..0088ebc0d8ad --- /dev/null +++ b/applications/CoSimulationApplication/custom_processes/data_transfer_3D_1D_process.cpp @@ -0,0 +1,353 @@ +// KRATOS / ___|___/ ___|(_)_ __ ___ _ _| | __ _| |_(_) ___ _ ___ +// | | / _ \___ \| | '_ ` _ \| | | | |/ _` | __| |/ _ \| '_ | +// | |__| (_) |__) | | | | | | | |_| | | (_| | |_| | (_) | | | | +// \____\___/____/|_|_| |_| |_|\__,_|_|\__,_|\__|_|\___/|_| |_| +// +// License: BSD License +// license: CoSimulationApplication/license.txt +// +// Main authors: Vicente Mataix Ferrandiz +// + +// System includes + +// External includes + +// Project includes +#include "custom_processes/data_transfer_3D_1D_process.h" +#include "utilities/parallel_utilities.h" +#include "utilities/reduction_utilities.h" +#include "utilities/atomic_utilities.h" +#include "utilities/variable_utils.h" +#include "spatial_containers/spatial_containers.h" // kd-tree + +namespace Kratos +{ + +#define KRATOS_KDTREE_POINTELEMENT_DEFINITION \ +using PointType = PointElement; \ +using PointTypePointer = PointType::Pointer; \ +using PointVector = std::vector; \ +using PointIterator = PointVector::iterator; \ +using DistanceVector = std::vector; \ +using DistanceIterator = DistanceVector::iterator; \ +using BucketType = Bucket<3ul, PointType, PointVector, PointTypePointer, PointIterator, DistanceIterator>; \ +using KDTree = Tree>; + +/** + * @brief This method determines the 1D model part + * @param rFirstModelPart The first ModelPart + * @param rSecondModelPart The second ModelPart + * @return The 1D model part + */ +ModelPart& Determine1DModelPart(ModelPart& rFirstModelPart, ModelPart& rSecondModelPart) +{ + // In MPI could be a bit troublesome + if (rFirstModelPart.IsDistributed()) { + KRATOS_ERROR << "Not compatible with MPI" << std::endl; + } else { // In serial we just check model part has 1D entities (local space) + // Determine if the modelparts have entities + + // First model part + if (rFirstModelPart.NumberOfElements() > 0 || rFirstModelPart.NumberOfConditions() > 0 ) { + if (rFirstModelPart.NumberOfConditions() > 0) { + if (rFirstModelPart.ConditionsBegin()->GetGeometry().LocalSpaceDimension() == 1) { + return rFirstModelPart; + } + } else { + if (rFirstModelPart.ElementsBegin()->GetGeometry().LocalSpaceDimension() == 1) { + return rFirstModelPart; + } + } + } + + // Second model part + if (rSecondModelPart.NumberOfElements() > 0 || rSecondModelPart.NumberOfConditions() > 0 ) { + if (rSecondModelPart.NumberOfConditions() > 0) { + if (rSecondModelPart.ConditionsBegin()->GetGeometry().LocalSpaceDimension() == 1) { + return rSecondModelPart; + } + } else { + if (rSecondModelPart.ElementsBegin()->GetGeometry().LocalSpaceDimension() == 1) { + return rSecondModelPart; + } + } + } + } + + KRATOS_ERROR << "Impossible to detect 1D model part" << std::endl; + return rFirstModelPart; +} + +/** + * @brief This method determines the 3D model part + * @details The counter part of the previous method + * @param rFirstModelPart The first ModelPart + * @param rSecondModelPart The second ModelPart + * @return The 3D model part + */ +ModelPart& Determine3DModelPart(ModelPart& rFirstModelPart, ModelPart& rSecondModelPart) +{ + // It is the counter part of the previous method + ModelPart* p_1d_model_part = &Determine1DModelPart(rFirstModelPart, rSecondModelPart); + if (p_1d_model_part == &rFirstModelPart) { + return rSecondModelPart; + } else { + return rFirstModelPart; + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +DataTransfer3D1DProcess::DataTransfer3D1DProcess( + ModelPart& rFirstModelPart, + ModelPart& rSecondModelPart, + Parameters ThisParameters + ) : mr3DModelPart(Determine3DModelPart(rFirstModelPart, rSecondModelPart)), + mr1DModelPart(Determine1DModelPart(rFirstModelPart, rSecondModelPart)), + mThisParameters(ThisParameters) +{ + // Validate deafult parameters + mThisParameters.ValidateAndAssignDefaults(GetDefaultParameters()); + + // Throwing and error if MPI is used + KRATOS_ERROR_IF(mr3DModelPart.IsDistributed() || mr1DModelPart.IsDistributed()) << "This implementation is not available for MPI model parts" << std::endl; + + // Getting mThisParameters + GetVariablesList(mOriginListVariables, mDestinationListVariables); + + // We generate auxiliary 3D model part + const bool has_intersected_model_part = mr3DModelPart.HasSubModelPart("IntersectedElements3D"); + auto& r_aux_model_part_3D = has_intersected_model_part ? mr3DModelPart.GetSubModelPart("IntersectedElements3D") : mr3DModelPart.CreateSubModelPart("IntersectedElements3D"); + if (has_intersected_model_part) { + VariableUtils().SetFlag(TO_ERASE, true, r_aux_model_part_3D.Nodes()); + VariableUtils().SetFlag(TO_ERASE, true, r_aux_model_part_3D.Elements()); + r_aux_model_part_3D.RemoveNodes(TO_ERASE); + r_aux_model_part_3D.RemoveElements(TO_ERASE); + } + + + // Max length of the elements considered + const double max_length_1d = GetMaxLength(mr1DModelPart); + const double max_length_3d = GetMaxLength(mr3DModelPart); + const double max_length = std::max(max_length_1d, max_length_3d); + + // The elements array + auto& r_elements_array = mr1DModelPart.Elements(); + const auto it_elem_begin = r_elements_array.begin(); + + /// KDtree definitions + KRATOS_KDTREE_POINTELEMENT_DEFINITION + + // Some auxiliary values + const IndexType allocation_size = ThisParameters["search_parameters"]["allocation_size"].GetInt(); // Allocation size for the vectors and max number of potential results + const double search_factor = ThisParameters["search_parameters"]["search_factor"].GetDouble(); // The search factor to be considered + const double search_increment_factor = ThisParameters["search_parameters"]["search_increment_factor"].GetDouble(); // The search increment factor to be considered + IndexType bucket_size = ThisParameters["search_parameters"]["bucket_size"].GetInt(); // Bucket size for kd-tree + + PointVector points_vector; + points_vector.reserve(r_elements_array.size()); + for (std::size_t i = 0; i < r_elements_array.size(); ++i) { + auto it_elem = it_elem_begin + i; + points_vector.push_back(PointTypePointer(new PointType(*(it_elem.base())))); + } + KDTree tree_points(points_vector.begin(), points_vector.end(), bucket_size); + + // Iterate over the elements (first assign NODAL_VOLUME) + const std::size_t aux_elem_number = block_for_each>(mr3DModelPart.Elements(), [&](Element& rElement) { + double search_radius = search_factor * max_length; + + // Initialize values + PointVector points_found(allocation_size); + IndexType number_points_found = 0; + auto& r_geometry_tetra = rElement.GetGeometry(); + + // Iterate doing the search + const Point center = r_geometry_tetra.Center(); + while (number_points_found == 0) { + search_radius *= search_increment_factor; + const PointElement point(center.X(), center.Y(), center.Z()); + number_points_found = tree_points.SearchInRadius(point, search_radius, points_found.begin(), allocation_size); + } + + // Getting the intersection points + bool intersected = false; + for (IndexType i = 0; i < number_points_found; ++i) { + auto p_point = points_found[i]; + auto p_geometry = p_point->pGetElement()->pGetGeometry(); + if (r_geometry_tetra.HasIntersection(*p_geometry)) { + intersected = true; + break; + } + } + // 1 point at least is required + if (intersected) { + rElement.Set(VISITED); + return 1; + } else { + return 0; + } + }); + + // We add the elements to the auxiliary model part + std::vector elements_id; + elements_id.reserve(aux_elem_number); + for (auto& r_elem : mr3DModelPart.Elements()) { + if (r_elem.Is(VISITED)) { + elements_id.push_back(r_elem.Id()); + } + } + r_aux_model_part_3D.AddElements(elements_id); + + // We add the nodes to the model part + for (auto& r_elem : r_aux_model_part_3D.Elements()) { + auto& r_geometry = r_elem.GetGeometry(); + for (auto& r_node : r_geometry) { + if (!r_aux_model_part_3D.HasNode(r_node.Id())) { + r_aux_model_part_3D.AddNode(&r_node); + } + } + } + + // Clear VISITED flag + block_for_each(r_aux_model_part_3D.Elements(), [](Element& rElement){ + rElement.Reset(VISITED); + }); + + // Check mapper exists + KRATOS_ERROR_IF_NOT(MapperFactoryType::HasMapper("nearest_element")) << "Mapper \"nearest_element\" not available. Please import MappingApplication" << std::endl; + + // Generate the mapper + Parameters mapping_variables = mThisParameters["mapping_variables"]; + mpMapper = MapperFactoryType::CreateMapper(r_aux_model_part_3D, mr1DModelPart, mapping_variables); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void DataTransfer3D1DProcess::Execute() +{ + KRATOS_ERROR_IF(mpMapper == nullptr) << "The mapper is not initialized" << std::endl; + + // From 3D to 1D + if (this->Is(MapperFlags::USE_TRANSPOSE)) { + InterpolateFrom3Dto1D(); + } else { // From 1D to 3D + InterpolateFrom1Dto3D(); + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void DataTransfer3D1DProcess::GetVariablesList( + std::vector*>& rOriginListVariables, + std::vector*>& rDestinationListVariables + ) +{ + // The components as an array of strings + const std::array direction_string({"X", "Y", "Z"}); + + // Getting variables + const std::vector origin_variables_names = mThisParameters["origin_variables"].GetStringArray(); + const std::vector destination_variables_names = mThisParameters["destination_variables"].GetStringArray(); + KRATOS_ERROR_IF(origin_variables_names.size() == 0) << "No variables defined" << std::endl; + KRATOS_ERROR_IF(origin_variables_names.size() != destination_variables_names.size()) << "Origin and destination variables do not coincide in size" << std::endl; + for (IndexType i_var = 0; i_var < origin_variables_names.size(); ++i_var) { + if (KratosComponents>::Has(origin_variables_names[i_var])) { + KRATOS_ERROR_IF_NOT(KratosComponents>::Has(destination_variables_names[i_var])) << "Destination variable is not a scalar" << std::endl; + rOriginListVariables.push_back(&KratosComponents>::Get(origin_variables_names[i_var])); + rDestinationListVariables.push_back(&KratosComponents>::Get(destination_variables_names[i_var])); + } else if (KratosComponents>>::Has(origin_variables_names[i_var])) { + const bool check_variable = KratosComponents>>::Get(destination_variables_names[i_var]); + KRATOS_ERROR_IF_NOT(check_variable) << "Destination variable is not an array of 3 components" << std::endl; + const auto& r_origin_var_name = origin_variables_names[i_var]; + const auto& r_destination_var_name = destination_variables_names[i_var]; + for (unsigned int i_comp = 0; i_comp < 3; ++i_comp) { + rOriginListVariables.push_back(&KratosComponents>::Get(r_origin_var_name + "_" + direction_string[i_comp])); + rDestinationListVariables.push_back(&KratosComponents>::Get(r_destination_var_name + "_" + direction_string[i_comp])); + } + } else { + KRATOS_ERROR << "Variable " << origin_variables_names[i_var] << " not available in KratosComponents as double or 3 components array" << std::endl; + } + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void DataTransfer3D1DProcess::InterpolateFrom1Dto3D() +{ + // Interpolate + Kratos::Flags mapper_flags = Kratos::Flags(); + const bool swap_sign = mThisParameters["swap_sign"].GetBool(); + if (swap_sign) mapper_flags.Set(MapperFlags::SWAP_SIGN); + for (std::size_t i_var = 0; i_var < mOriginListVariables.size(); ++i_var) { + mpMapper->InverseMap(*mDestinationListVariables[i_var], *mOriginListVariables[i_var], mapper_flags); + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void DataTransfer3D1DProcess::InterpolateFrom3Dto1D() +{ + // Interpolate + Kratos::Flags mapper_flags = Kratos::Flags(); + const bool swap_sign = mThisParameters["swap_sign"].GetBool(); + if (swap_sign) mapper_flags.Set(MapperFlags::SWAP_SIGN); + for (std::size_t i_var = 0; i_var < mOriginListVariables.size(); ++i_var) { + mpMapper->Map(*mOriginListVariables[i_var], *mDestinationListVariables[i_var], mapper_flags); + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +double DataTransfer3D1DProcess::GetMaxLength(const ModelPart& rModelPart) +{ + KRATOS_ERROR_IF(rModelPart.NumberOfElements() == 0 && rModelPart.NumberOfConditions() == 0) << "Empty model part" << std::endl; + if (rModelPart.NumberOfConditions() > 0) { + auto& r_conditions_array = rModelPart.Conditions(); + return block_for_each>(r_conditions_array, [&](const auto& rCondition) { + return rCondition.GetGeometry().Length(); + }); + } else { + auto& r_elements_array = rModelPart.Elements(); + return block_for_each>(r_elements_array, [&](const auto& rElement) { + return rElement.GetGeometry().Length(); + }); + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +const Parameters DataTransfer3D1DProcess::GetDefaultParameters() const +{ + Parameters default_parameters = Parameters(R"( + { + "origin_variables" : [], + "destination_variables" : [], + "swap_sign" : false, + "mapping_variables" : { + "mapper_type" : "nearest_element", + "echo_level" : 0, + "search_settings" : { + "max_num_search_iterations" : 8, + "echo_level" : 0 + } + }, + "search_parameters" : { + "allocation_size" : 100, + "bucket_size" : 4, + "search_factor" : 2.0, + "search_increment_factor" : 1.5 + } + })" ); + + return default_parameters; +} + +} // namespace Kratos \ No newline at end of file diff --git a/applications/CoSimulationApplication/custom_processes/data_transfer_3D_1D_process.h b/applications/CoSimulationApplication/custom_processes/data_transfer_3D_1D_process.h new file mode 100644 index 000000000000..9e9f2cb3e4df --- /dev/null +++ b/applications/CoSimulationApplication/custom_processes/data_transfer_3D_1D_process.h @@ -0,0 +1,232 @@ +// KRATOS / ___|___/ ___|(_)_ __ ___ _ _| | __ _| |_(_) ___ _ ___ +// | | / _ \___ \| | '_ ` _ \| | | | |/ _` | __| |/ _ \| '_ | +// | |__| (_) |__) | | | | | | | |_| | | (_| | |_| | (_) | | | | +// \____\___/____/|_|_| |_| |_|\__,_|_|\__,_|\__|_|\___/|_| |_| +// +// License: BSD License +// license: CoSimulationApplication/license.txt +// +// Main authors: Vicente Mataix Ferrandiz +// + +#pragma once + +// System includes + +// External includes + +// Project includes +#include "includes/define.h" +#include "includes/model_part.h" +#include "includes/kratos_parameters.h" +#include "processes/process.h" + +/* The mappers includes */ +#include "spaces/ublas_space.h" +#include "mappers/mapper_flags.h" +#include "factories/mapper_factory.h" + +namespace Kratos +{ +///@addtogroup CoSimulationApplication +///@{ + +///@} +///@name Functions +///@{ + +///@} +///@name Kratos Classes +///@{ + +/** + * @class PointElement + * @ingroup CoSimulationApplication + * @brief Custom Point container to be used by the search + * @details It stores the pointer of a certain element + * @author Vicente Mataix Ferrandiz + */ +class PointElement + : public Point +{ +public: + + ///@name Type Definitions + ///@{ + + /// Base class definition + typedef Point BaseType; + + /// Counted pointer of PointElement + KRATOS_CLASS_POINTER_DEFINITION( PointElement ); + + ///@} + ///@name Life Cycle + ///@{ + + /// Default constructors + PointElement(): + BaseType() + {} + + PointElement(const double X, const double Y, const double Z) + : BaseType(X, Y, Z) + {} + + PointElement(Element::Pointer pElement): + mpElement(pElement) + { + UpdatePoint(); + } + + + ///@} + ///@name Operators + ///@{ + + ///@} + ///@name Operations + ///@{ + + /** + * @brief Returns the element associated to the point + * @return mpElement The reference to the element associated to the point + */ + Element::Pointer pGetElement() + { + return mpElement; + } + + /** + * @brief This function updates the database, using as base for the coordinates the condition center + */ + void UpdatePoint() + { + noalias(this->Coordinates()) = mpElement->GetGeometry().Center().Coordinates(); + } + +private: + ///@} + ///@name Member Variables + ///@{ + + Element::Pointer mpElement = nullptr; // The element instance + + ///@} + +}; // Class PointElement + +#define DEFINE_MAPPER_FACTORY_SERIAL \ +using SparseSpace = UblasSpace, boost::numeric::ublas::vector>; \ +using DenseSpace = UblasSpace, DenseVector>; \ +using MapperFactoryType = MapperFactory; \ +using MapperType = Mapper; + +/** + * @class DataTransfer3D1DProcess + * @ingroup CoSimulationApplication + * @brief This utility includes auxiliary methods to transfer from 3D domains to 1D domains and viceversa + * @author Vicente Mataix Ferrandiz + */ +class KRATOS_API(CO_SIMULATION_APPLICATION) DataTransfer3D1DProcess + : public Process +{ +public: + ///@name Type Definitions + ///@{ + + /// Pointer definition of DataTransfer3D1DProcess + KRATOS_CLASS_POINTER_DEFINITION(DataTransfer3D1DProcess); + + /// Geometry definition + using GeometryType = Geometry; + + // Define mapper factory + DEFINE_MAPPER_FACTORY_SERIAL + + ///@} + ///@name Life Cycle + ///@{ + + /// Default constructor. + DataTransfer3D1DProcess( + ModelPart& rFirstModelPart, + ModelPart& rSecondModelPart, + Parameters ThisParameters = Parameters(R"({})") + ); + + + ///@} + ///@name Operators + ///@{ + + void operator()() + { + Execute(); + } + + ///@} + ///@name Operations + ///@{ + + /** + * @brief This method executes the process + */ + void Execute() override; + + /** + * @brief This method provides the defaults parameters to avoid conflicts between the different constructors + */ + const Parameters GetDefaultParameters() const override; + + ///@} +private: + ///@name Private member Variables + ///@{ + + ModelPart& mr3DModelPart; /// The 3D model part + ModelPart& mr1DModelPart; /// The 1D model part + Parameters mThisParameters; /// The parameters containing the configuration + MapperType::Pointer mpMapper = nullptr; /// The mapper pointer + std::vector*> mOriginListVariables; /// The list of variables to be transferred (origin) + std::vector*> mDestinationListVariables; /// The list of variables to be transferred (destination) + + ///@} + ///@name Private Operations + ///@{ + + /** + * @brief This method interpolates from the 3D to the 1D + */ + void InterpolateFrom3Dto1D(); + + /** + * @brief This method interpolates from the 1D to the 3D + */ + void InterpolateFrom1Dto3D(); + + /** + * @brief This method computes the list of variables to interpolate + * @param rOriginListVariables The list of origin variables to interpolate + * @param rDestinationListVariables The list of destination variables to interpolate + */ + void GetVariablesList( + std::vector*>& rOriginListVariables, + std::vector*>& rDestinationListVariables + ); + + /** + * @brief This method computes maximum length of the elements + * @param rModelPart The model part to compute + * @return The maximum length + */ + static double GetMaxLength(const ModelPart& rModelPart); + + ///@} +}; // Class DataTransfer3D1DProcess + +///@} + +///@} addtogroup block + +} // namespace Kratos. \ No newline at end of file diff --git a/applications/CoSimulationApplication/custom_python/add_custom_processes_to_python.cpp b/applications/CoSimulationApplication/custom_python/add_custom_processes_to_python.cpp new file mode 100644 index 000000000000..45b4d1077261 --- /dev/null +++ b/applications/CoSimulationApplication/custom_python/add_custom_processes_to_python.cpp @@ -0,0 +1,33 @@ +// KRATOS / ___|___/ ___|(_)_ __ ___ _ _| | __ _| |_(_) ___ _ ___ +// | | / _ \___ \| | '_ ` _ \| | | | |/ _` | __| |/ _ \| '_ | +// | |__| (_) |__) | | | | | | | |_| | | (_| | |_| | (_) | | | | +// \____\___/____/|_|_| |_| |_|\__,_|_|\__,_|\__|_|\___/|_| |_| +// +// License: BSD License +// license: CoSimulationApplication/license.txt +// +// Main authors: +// + +// System includes + +// External includes + +// Project includes +#include "includes/define.h" +#include "custom_python/add_custom_processes_to_python.h" +#include "custom_processes/data_transfer_3D_1D_process.h" + +namespace Kratos::Python{ + + void AddCustomProcessesToPython(pybind11::module& m) + { + namespace py = pybind11; + + py::class_< DataTransfer3D1DProcess, DataTransfer3D1DProcess::Pointer, Process>(m, "DataTransfer3D1DProcess") + .def(py::init()) + .def(py::init()) + ; + } + +} // namespace Kratos::Python. diff --git a/applications/CoSimulationApplication/custom_python/add_custom_processes_to_python.h b/applications/CoSimulationApplication/custom_python/add_custom_processes_to_python.h new file mode 100644 index 000000000000..9674290a87f5 --- /dev/null +++ b/applications/CoSimulationApplication/custom_python/add_custom_processes_to_python.h @@ -0,0 +1,26 @@ +// KRATOS / ___|___/ ___|(_)_ __ ___ _ _| | __ _| |_(_) ___ _ ___ +// | | / _ \___ \| | '_ ` _ \| | | | |/ _` | __| |/ _ \| '_ | +// | |__| (_) |__) | | | | | | | |_| | | (_| | |_| | (_) | | | | +// \____\___/____/|_|_| |_| |_|\__,_|_|\__,_|\__|_|\___/|_| |_| +// +// License: BSD License +// license: CoSimulationApplication/license.txt +// +// Main authors: +// + +#pragma once + +// System includes +#include + +// External includes + +// Project includes +#include "includes/define.h" + +namespace Kratos::Python { + +void AddCustomProcessesToPython(pybind11::module& m); + +} // namespace Kratos::Python. diff --git a/applications/CoSimulationApplication/custom_python/co_simulation_python_application.cpp b/applications/CoSimulationApplication/custom_python/co_simulation_python_application.cpp index cc150245d6f3..1ef1887e5137 100644 --- a/applications/CoSimulationApplication/custom_python/co_simulation_python_application.cpp +++ b/applications/CoSimulationApplication/custom_python/co_simulation_python_application.cpp @@ -23,6 +23,7 @@ #include "custom_python/add_co_sim_io_to_python.h" #include "custom_python/add_custom_io_to_python.h" #include "custom_python/add_custom_utilities_to_python.h" +#include "custom_python/add_custom_processes_to_python.h" namespace Kratos::Python { @@ -39,6 +40,7 @@ PYBIND11_MODULE(KratosCoSimulationApplication,m) AddCoSimIOToPython(m); AddCustomIOToPython(m); AddCustomUtilitiesToPython(m); + AddCustomProcessesToPython(m); KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, SCALAR_DISPLACEMENT ); KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, SCALAR_ROOT_POINT_DISPLACEMENT ); diff --git a/applications/CoSimulationApplication/python_scripts/data_transfer_operators/3d_1d_data_transfer.py b/applications/CoSimulationApplication/python_scripts/data_transfer_operators/3d_1d_data_transfer.py new file mode 100644 index 000000000000..d6c1d0b864fd --- /dev/null +++ b/applications/CoSimulationApplication/python_scripts/data_transfer_operators/3d_1d_data_transfer.py @@ -0,0 +1,288 @@ +# Importing the base class +from KratosMultiphysics.CoSimulationApplication.base_classes.co_simulation_data_transfer_operator import CoSimulationDataTransferOperator + +# Importing the Kratos Library +import KratosMultiphysics as KM +import KratosMultiphysics.MappingApplication as KratosMapping + +# CoSimulation imports +import KratosMultiphysics.CoSimulationApplication as KratosCoSim +import KratosMultiphysics.CoSimulationApplication.co_simulation_tools as cs_tools +import KratosMultiphysics.CoSimulationApplication.colors as colors + +def Create(*args): + return Kratos3D1DDataTransferOperator(*args) + +class Kratos3D1DDataTransferOperator(CoSimulationDataTransferOperator): + """DataTransferOperator that transfers values from one 3D interface (ModelPart) to another 1D interface (ModelPart). + The data_transfer_3d_1ds of the Kratos-MappingApplication are used + """ + + def __init__(self, settings, parent_coupled_solver_data_communicator): + """ + Initializes the data transfer object with the given settings and sets up + internal structures for managing data transfer processes between solvers. + + Parameters: + ---------- + settings : KratosMultiphysics.Parameters + A Parameters object containing the configuration for data transfer. + Must include `"3d_1d_data_transfer_settings"` to specify relevant + transfer settings. + + parent_coupled_solver_data_communicator : CoupledSolverDataCommunicator + The communicator object that facilitates data sharing between solvers + in the coupled analysis. + + Raises: + ------ + Exception + If the function is not called within an MPI distributed environment. + Exception + If the required `"3d_1d_data_transfer_settings"` is missing in `settings`. + + Notes: + ------ + - This constructor is intended for usage in an MPI environment only. It + checks for distributed mode (`KM.IsDistributedRun()`), ensuring the + function is not misused outside of MPI. + - The `__data_transfer_process` dictionary is initialized as an empty + structure for storing and managing the data transfer processes. + - The `origin_is_3d` attribute is set to `None` initially and will be determined + based on the geometry of the model part in subsequent processes. + """ + if KM.IsDistributedRun(): + raise Exception("This function can only be called when Kratos is running in MPI!") + + if not settings.Has("3d_1d_data_transfer_settings"): + raise Exception('No "3d_1d_data_transfer_settings" provided!') + + # Initialize the parent class with the provided settings and communicator + super().__init__(settings, parent_coupled_solver_data_communicator) + + # Initialize attributes for tracking geometry dimension and data transfer processes + self.origin_is_3d = None + self.__data_transfer_process = {} + + def _ExecuteTransferData(self, from_solver_data, to_solver_data, transfer_options): + """ + Executes the data transfer between two solvers, either by retrieving and using an existing + data transfer process or creating a new one if none exists. The function also allows for + specific transfer options such as toggling the sign of the transfer data. + + Parameters: + ---------- + from_solver_data : SolverData + Data object containing information about the origin solver's model part and variable to transfer. + + to_solver_data : SolverData + Data object containing information about the destination solver's model part and variable to receive data. + + transfer_options : TransferOptions + Additional options for data transfer, such as swapping the sign of transferred values. + + Notes: + ------ + The function uses internal dictionaries to manage transfer processes, utilizing unique identifiers + based on model part names and variables to ensure consistent data mapping between solvers. + + Steps: + ------ + 1. **Identify or Create Transfer Process**: + - Identifiers for the origin and destination model parts are generated, combining solver names, + model part names, and variable names. + - If a transfer process already exists between the identified model parts, it is reused. + - If no transfer process exists, a new one is created with the provided settings and options. + + 2. **Set Transfer Parameters**: + - Custom settings are created and added to manage specific transfer configurations, + such as swapping signs or adjusting for 3D/1D transfers. + + 3. **Execute Transfer**: + - Depending on whether the origin is in 3D, specific transpose settings are applied, and the + transfer is executed. + + Raises: + ------ + ValueError + If the `transfer_options` provided are invalid or unsupported. + """ + model_part_origin_name = from_solver_data.model_part_name + variable_origin = from_solver_data.variable + identifier_origin = from_solver_data.solver_name + "." + model_part_origin_name + "." + variable_origin.Name() + + model_part_destination_name = to_solver_data.model_part_name + variable_destination = to_solver_data.variable + identifier_destination = to_solver_data.solver_name + "." + model_part_destination_name + "." + variable_destination.Name() + + identifier_tuple = (identifier_origin, identifier_destination) + inverse_identifier_tuple = (identifier_destination, identifier_origin) + + if identifier_tuple in self.__data_transfer_process: + if self.origin_is_3d == True: + self.__data_transfer_process[identifier_tuple].Set(KM.Mapper.USE_TRANSPOSE) + self.__data_transfer_process[identifier_tuple].Execute() + if self.origin_is_3d == True: + self.__data_transfer_process[identifier_tuple].Reset(KM.Mapper.USE_TRANSPOSE) + elif inverse_identifier_tuple in self.__data_transfer_process: + if self.origin_is_3d == False: + self.__data_transfer_process[inverse_identifier_tuple].Set(KM.Mapper.USE_TRANSPOSE) + self.__data_transfer_process[inverse_identifier_tuple].Execute() + if self.origin_is_3d == False: + self.__data_transfer_process[inverse_identifier_tuple].Reset(KM.Mapper.USE_TRANSPOSE) + else: + model_part_origin = from_solver_data.GetModelPart() + model_part_destination = to_solver_data.GetModelPart() + + # Logging information if echo level is set + if self.echo_level > 0: + info_msg = ( + f"Creating Mapper:\n Origin: ModelPart '{model_part_origin_name}' " + f"of solver '{from_solver_data.solver_name}'\n" + f" Destination: ModelPart '{model_part_destination_name}' " + f"of solver '{to_solver_data.solver_name}'" + ) + cs_tools.cs_print_info(colors.bold(self._ClassName()), info_msg) + + # Creating parameters for the data transfer + # TODO: I should check if the parameters change between different data transfers + parameters = KM.Parameters(self.settings["3d_1d_data_transfer_settings"].WriteJsonString()) + if not parameters.Has("origin_variables"): + parameters.AddEmptyArray("origin_variables") + parameters["origin_variables"].Append(variable_origin.Name()) + if not parameters.Has("destination_variables"): + parameters.AddEmptyArray("destination_variables") + parameters["destination_variables"].Append(variable_destination.Name()) + for transfer_option in transfer_options.GetStringArray(): + if transfer_option == "swap_sign": + parameters["swap_sign"].SetBool(True) + + # Check if the origin model part is 3D or 1D + self.origin_is_3d = self.__check_model_part_3D(model_part_origin) + # Clone is necessary because the settings are validated and defaults assigned, which could influence the creation of other data transfers + self.__data_transfer_process[identifier_tuple] = KratosCoSim.DataTransfer3D1DProcess(model_part_origin, model_part_destination, parameters.Clone()) + self.__data_transfer_process[inverse_identifier_tuple] = KratosCoSim.DataTransfer3D1DProcess(model_part_origin, model_part_destination, parameters.Clone()) + + # Execute the data transfer + if self.origin_is_3d: + self.__data_transfer_process[identifier_tuple].Set(KM.Mapper.USE_TRANSPOSE) + self.__data_transfer_process[identifier_tuple].Execute() + if self.origin_is_3d: + self.__data_transfer_process[identifier_tuple].Reset(KM.Mapper.USE_TRANSPOSE) + + def _Check(self, from_solver_data, to_solver_data): + """ + Validates that the `from_solver_data` and `to_solver_data` contain only + nodal data for transfer operations. Raises an exception if the data + location is not at the nodes. + + Parameters: + ---------- + from_solver_data : SolverData + The source data to check, containing information about the solver and + data location in the model part. + + to_solver_data : SolverData + The destination data to check, containing information about the solver + and data location in the model part. + + Raises: + ------ + Exception + If the location in `from_solver_data` or `to_solver_data` is not at + the nodes, an exception is raised specifying the issue, the model part, + and the solver name. + + Notes: + ------ + This function assumes `from_solver_data` and `to_solver_data` are instances + containing a `location` attribute to specify where the data is located, as + well as `model_part_name` and `solver_name` attributes for error reporting. + """ + def CheckData(data_to_check): + if "node" not in data_to_check.location: + raise Exception( + f'Transfer only supports nodal values! "{self._ClassName()}"\n' + f'Checking ModelPart "{data_to_check.model_part_name}" of solver "{data_to_check.solver_name}"' + ) + + # Perform the nodal data location check on both the source and destination data + CheckData(from_solver_data) + CheckData(to_solver_data) + + def __check_model_part_3D(self, model_part): + """ + Checks whether the given `model_part` consists entirely of 3D elements and + verifies that all elements have the same geometry type. + + Parameters: + ---------- + model_part : ModelPart + The model part to check, which should contain elements representing a 3D structure. + Returns: + ------- + bool + True if all elements are 3D and have the same geometry type, False if any + element is 1D. + Raises: + ------ + ValueError + If elements in `model_part` have different geometry types, an exception is raised. + """ + first_geometry_type = None + + # First check the conditions + for cond in model_part.Conditions: + geom = cond.GetGeometry() + # Check for consistent geometry type + if first_geometry_type is None: + first_geometry_type = geom.LocalSpaceDimension() + elif geom.LocalSpaceDimension() != first_geometry_type: + raise ValueError("Inconsistent geometry types found in model_part.") + + # Now checking elements + if first_geometry_type is None: + for elem in model_part.Elements: + geom = elem.GetGeometry() + # Check for consistent geometry type + if first_geometry_type is None: + first_geometry_type = geom.LocalSpaceDimension() + elif geom.LocalSpaceDimension() != first_geometry_type: + raise ValueError("Inconsistent geometry types found in model_part.") + + # Raise error if no geometries + if first_geometry_type is None: + raise ValueError("Neither elements or conditions defined.") + + # Check if the element is 1D + return first_geometry_type != 1 + + @classmethod + def _GetDefaultParameters(cls): + this_defaults = KM.Parameters("""{ + "3d_1d_data_transfer_settings" : { + "origin_variables" : [], + "destination_variables" : [], + "swap_sign" : false, + "interpolate_parameters" : { + "data_transfer_3d_1d_type" : "nearest_element", + "echo_level" : 0, + "search_settings" : { + "max_num_search_iterations" : 8, + "echo_level" : 0 + } + }, + "search_parameters" : { + "allocation_size" : 100, + "bucket_size" : 4, + "search_factor" : 2.0, + "search_increment_factor" : 1.5 + } + } + }""") + this_defaults.AddMissingParameters(super()._GetDefaultParameters()) + return this_defaults + + @classmethod + def _GetListAvailableTransferOptions(cls): + return ["swap_sign"] diff --git a/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_mdpa/block.mdpa b/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_mdpa/block.mdpa new file mode 100644 index 000000000000..aaae8e037523 --- /dev/null +++ b/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_mdpa/block.mdpa @@ -0,0 +1,166 @@ +Begin Properties 0 +End Properties + +Begin Nodes + 1 5.0000000000e-01 7.4019237890e-01 5.0000000000e+00 + 2 2.5980762110e-01 5.0000000000e-01 5.0000000000e+00 + 3 3.4281892070e-01 4.1828680480e-01 4.7441145114e+00 + 4 5.6443695940e-01 5.7364978990e-01 4.6327669521e+00 + 5 5.9938576100e-01 4.8603785390e-01 4.8753320697e+00 + 6 5.0000000000e-01 2.5980762110e-01 5.0000000000e+00 + 7 7.4019237890e-01 5.0000000000e-01 5.0000000000e+00 + 8 4.4003055460e-01 4.5291315690e-01 4.4212714469e+00 + 9 5.0835079890e-01 6.8726318170e-01 4.2911963632e+00 + 10 3.7200237220e-01 4.5478541020e-01 4.3394586377e+00 + 11 6.1208638080e-01 6.6161285160e-01 4.3784836326e+00 + 12 5.5135720070e-01 2.5046482880e-01 4.8020602923e+00 + 13 4.1520677510e-01 2.2360679770e-01 4.5649314626e+00 + 14 7.0868671920e-01 4.1898814340e-01 4.4606845428e+00 + 15 3.9082250900e-01 5.5226420310e-01 4.0843225002e+00 + 16 6.1027293620e-01 4.4820298840e-01 4.2068619868e+00 + 17 5.7722290800e-01 7.3480283290e-01 4.0288610617e+00 + 18 4.0998842460e-01 2.8315208640e-01 4.1301895383e+00 + 19 7.7639320230e-01 6.6666666670e-01 4.1176470588e+00 + 20 5.7181411580e-01 4.0285112420e-01 3.9447709968e+00 + 21 3.3254720760e-01 4.9054913050e-01 3.8037787124e+00 + 22 5.7040135210e-01 6.0815454050e-01 3.8450816703e+00 + 23 4.6946784150e-01 5.9085366100e-01 3.6290076285e+00 + 24 4.6828132290e-01 3.2180228850e-01 3.6746049744e+00 + 25 7.1882727580e-01 5.1415189520e-01 3.6074842591e+00 + 26 3.3627944000e-01 5.1366730810e-01 3.3549744710e+00 + 27 5.0806806330e-01 3.2597022250e-01 3.4562997362e+00 + 28 6.0495306620e-01 5.6578248810e-01 3.3662547425e+00 + 29 4.5852415320e-01 6.4795236530e-01 3.1500500928e+00 + 30 5.0857858630e-01 3.8145404180e-01 3.1881231920e+00 + 31 5.9314095450e-01 7.7639320230e-01 3.0882352941e+00 + 32 3.5104778870e-01 7.8603846070e-01 2.9499073328e+00 + 33 7.7639320230e-01 4.0685904550e-01 3.0882352941e+00 + 34 5.0041761040e-01 4.4381291700e-01 2.9449284529e+00 + 35 6.7339901000e-01 5.3280185780e-01 2.8701266480e+00 + 36 3.3940215500e-01 4.9924503800e-01 2.7986033692e+00 + 37 4.1813630890e-01 2.2360679770e-01 2.7955905052e+00 + 38 5.4918201330e-01 3.5907848490e-01 2.7090064436e+00 + 39 4.9955044090e-01 5.5407116860e-01 2.6504288930e+00 + 40 3.4961282290e-01 3.7527708350e-01 2.5199722258e+00 + 41 6.6093036990e-01 2.1955882950e-01 2.6199498219e+00 + 42 6.6546140670e-01 5.4313693130e-01 2.4248122610e+00 + 43 4.4218446050e-01 6.6834201630e-01 2.3249399338e+00 + 44 4.9188520700e-01 4.2581410360e-01 2.3250365404e+00 + 45 6.1806236640e-01 6.1093696420e-01 2.1859281881e+00 + 46 3.3355134810e-01 7.7639320230e-01 2.0622910796e+00 + 47 4.9356598460e-01 3.9738209530e-01 2.0993851876e+00 + 48 5.7208948270e-01 5.8953548620e-01 1.9251963759e+00 + 49 5.7148434340e-01 3.2725635420e-01 1.8476097830e+00 + 50 3.8564387860e-01 4.4053408670e-01 1.7928026981e+00 + 51 4.0650607010e-01 7.0926222380e-01 1.7467250078e+00 + 52 6.2847509890e-01 5.2515517220e-01 1.6217611273e+00 + 53 4.7602604430e-01 3.0467495140e-01 1.5957714247e+00 + 54 3.5448682630e-01 5.2651637820e-01 1.5383051287e+00 + 55 1.6155591140e-01 3.1203915160e-01 1.5158200880e+00 + 56 5.9145297970e-01 5.3722735940e-01 1.4093874710e+00 + 57 3.4827223890e-01 3.7631316350e-01 1.3471348899e+00 + 58 4.3326439360e-01 5.9628718520e-01 1.2601791012e+00 + 59 1.9840518100e-01 4.1283099640e-01 1.1128282095e+00 + 60 7.1427636360e-01 5.1190773800e-01 1.1272310392e+00 + 61 5.2466539710e-01 3.8136616280e-01 1.1142290547e+00 + 62 4.0806450050e-01 5.9596043410e-01 1.0134037105e+00 + 63 1.6155591140e-01 3.1077229470e-01 9.3008177740e-01 + 64 3.9980022880e-01 2.9355734770e-01 8.7587741100e-01 + 65 5.9714137200e-01 4.7200176320e-01 8.6831493500e-01 + 66 3.5612183670e-01 5.6361992650e-01 7.4619923280e-01 + 67 3.0183624550e-01 4.4538042710e-01 6.9455637370e-01 + 68 5.6392535720e-01 3.1100278840e-01 6.5604399090e-01 + 69 5.7989863740e-01 5.7132310000e-01 6.0844736300e-01 + 70 3.2912448310e-01 2.2360679770e-01 5.9410829900e-01 + 71 3.6000041600e-01 5.3030603290e-01 4.5859261560e-01 + 72 5.6730310670e-01 6.8440108380e-01 3.7149917810e-01 + 73 6.0620616480e-01 4.1520222180e-01 3.9102231700e-01 + 74 7.4468079850e-01 5.5938095630e-01 2.0463118700e-01 + 75 4.7562166220e-01 5.1857178180e-01 1.7408453330e-01 + 76 6.5273608740e-01 2.6357219280e-01 1.2089350140e-01 + 77 2.5980762110e-01 5.0000000000e-01 0.0000000000e+00 + 78 5.0000000000e-01 7.4019237890e-01 0.0000000000e+00 + 79 5.0000000000e-01 2.5980762110e-01 0.0000000000e+00 + 80 7.4019237890e-01 5.0000000000e-01 0.0000000000e+00 +End Nodes + +Begin Elements Element3D4N + 1 0 5 6 7 2 + 2 0 14 4 11 8 + 3 0 14 4 8 13 + 4 0 12 3 5 4 + 5 0 11 14 8 16 + 6 0 18 16 20 15 + 7 0 16 20 15 19 + 8 0 18 16 15 10 + 9 0 10 8 16 11 + 10 0 16 15 10 9 + 11 0 10 16 9 11 + 12 0 12 13 3 4 + 13 0 5 7 1 2 + 14 0 45 48 47 46 + 15 0 52 48 51 50 + 16 0 42 43 45 47 + 17 0 47 48 49 50 + 18 0 48 49 50 52 + 19 0 50 52 53 54 + 20 0 54 53 55 56 + 21 0 56 52 54 53 + 22 0 58 56 54 55 + 23 0 51 52 50 54 + 24 0 62 58 61 65 + 25 0 61 62 65 59 + 26 0 60 56 58 61 + 27 0 56 58 61 57 + 28 0 56 58 57 55 + 29 0 58 60 61 65 + 30 0 62 66 65 63 + 31 0 69 72 73 71 + 32 0 73 69 71 70 + 33 0 73 72 74 75 + 34 0 73 72 75 71 + 35 0 74 73 75 76 + 36 0 75 74 76 80 + 37 0 75 76 79 80 + 38 0 68 70 67 69 + 39 0 67 68 69 64 + 40 0 65 66 69 67 + 41 0 65 66 67 63 + 42 0 69 65 67 64 + 43 0 38 35 39 37 + 44 0 39 38 37 40 + 45 0 39 38 40 41 + 46 0 3 5 4 1 + 47 0 39 35 36 37 + 48 0 34 37 36 35 + 49 0 36 34 35 32 + 50 0 44 42 47 43 + 51 0 44 42 43 40 + 52 0 62 59 63 65 + 53 0 69 67 71 70 + 54 0 21 24 23 20 + 55 0 23 26 28 24 + 56 0 31 29 32 34 + 57 0 28 26 29 30 + 58 0 29 28 30 33 + 59 0 28 26 30 27 + 60 0 28 26 27 24 + 61 0 23 28 25 24 + 62 0 21 15 20 17 + 63 0 31 29 34 33 + 64 0 21 23 22 20 + 65 0 22 21 20 17 + 66 0 47 48 50 46 + 67 0 25 23 24 20 + 68 0 15 20 17 19 + 69 0 13 8 3 4 + 70 0 29 30 34 33 + 71 0 43 45 47 46 + 72 0 43 42 39 40 + 73 0 42 39 40 41 + 74 0 75 77 78 80 + 75 0 35 34 31 32 + 76 0 75 79 77 80 + 77 0 5 3 2 1 +End Elements + diff --git a/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_mdpa/torus3d.mdpa b/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_mdpa/torus3d.mdpa new file mode 100644 index 000000000000..3ce2c1cb4410 --- /dev/null +++ b/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_mdpa/torus3d.mdpa @@ -0,0 +1,154 @@ +Begin Properties 0 +End Properties + +Begin Nodes + 1 2.9808074840e-01 -1.1021175054e+00 -7.0205989000e-03 + 2 -0.0000000000e+00 -1.2174811220e+00 2.0664452950e-01 + 3 -0.0000000000e+00 -1.1568605599e+00 -2.5572400110e-01 + 4 4.8516237780e-01 -1.0708726508e+00 -2.4320238390e-01 + 5 1.0637993110e-01 -9.0650575860e-01 5.8233869500e-02 + 6 5.1255932080e-01 -9.2964041530e-01 8.4175679600e-02 + 7 -1.9946137300e-01 -1.0790266848e+00 -4.0864338300e-02 + 8 1.8982901160e-01 -8.0989941880e-01 -2.4844539150e-01 + 9 -0.0000000000e+00 -8.5400868260e-01 -2.6208116160e-01 + 10 2.6787840270e-01 -6.4671567270e-01 0.0000000000e+00 + 11 -0.0000000000e+00 -7.0000000000e-01 0.0000000000e+00 + 12 -1.8982901160e-01 -8.0989941880e-01 -2.4844539150e-01 + 13 5.9451875670e-01 -7.3099586720e-01 -2.9438622350e-01 + 14 -2.1424131440e-01 -7.5845619390e-01 2.1239769560e-01 + 15 8.1947977940e-01 -7.0316920800e-01 5.7577770000e-04 + 16 4.9497474680e-01 -4.9497474680e-01 0.0000000000e+00 + 17 -2.6787840270e-01 -6.4671567270e-01 0.0000000000e+00 + 18 -5.4568911840e-01 -9.1850197500e-01 4.6675922300e-02 + 19 -5.9844766760e-01 -7.2778286530e-01 -2.9438622350e-01 + 20 7.4360019570e-01 -3.7286843680e-01 -2.4844539150e-01 + 21 6.4671567270e-01 -2.6787840270e-01 0.0000000000e+00 + 22 -4.9497474680e-01 -4.9497474680e-01 0.0000000000e+00 + 23 9.7705908600e-01 -4.2033683210e-01 5.0830205500e-02 + 24 -8.1947984150e-01 -7.0316915120e-01 5.7579440000e-04 + 25 8.0566980400e-01 -1.4037124450e-01 2.3833931350e-01 + 26 8.1944637150e-01 -1.2013846130e-01 -2.4594086730e-01 + 27 7.0000000000e-01 -0.0000000000e+00 0.0000000000e+00 + 28 -6.4671567270e-01 -2.6787840270e-01 0.0000000000e+00 + 29 -7.4360019580e-01 -3.7286843680e-01 -2.4844539150e-01 + 30 1.0512719803e+00 -1.0547482640e-01 -8.9361513000e-03 + 31 -8.8506975170e-01 -4.4556477460e-01 2.9986185390e-01 + 32 7.8730877470e-01 1.0277656680e-01 2.1808109360e-01 + 33 -9.7705908600e-01 -4.2033683210e-01 5.0830205500e-02 + 34 8.8324301090e-01 1.2240540840e-01 -2.7976375550e-01 + 35 -7.8730877470e-01 -1.0277656680e-01 2.1808109360e-01 + 36 6.4671567270e-01 2.6787840270e-01 0.0000000000e+00 + 37 -7.0000000000e-01 0.0000000000e+00 0.0000000000e+00 + 38 1.0512719804e+00 1.0547482640e-01 -8.9361513000e-03 + 39 -8.8324301090e-01 -1.2240540840e-01 -2.7976375550e-01 + 40 7.4360019580e-01 3.7286843680e-01 -2.4844539150e-01 + 41 -1.0512719804e+00 -1.0547482640e-01 -8.9361513000e-03 + 42 -8.1944637150e-01 1.2013846130e-01 -2.4594086730e-01 + 43 -8.0566980400e-01 1.4037124450e-01 2.3833931350e-01 + 44 4.9497474680e-01 4.9497474680e-01 0.0000000000e+00 + 45 -6.4671567270e-01 2.6787840270e-01 0.0000000000e+00 + 46 9.7705908600e-01 4.2033683210e-01 5.0830205500e-02 + 47 8.8506975170e-01 4.4556477460e-01 2.9986185390e-01 + 48 -1.0512719803e+00 1.0547482640e-01 -8.9361513000e-03 + 49 2.6787840270e-01 6.4671567270e-01 0.0000000000e+00 + 50 -4.9497474680e-01 4.9497474680e-01 0.0000000000e+00 + 51 -7.4360019570e-01 3.7286843680e-01 -2.4844539150e-01 + 52 0.0000000000e+00 7.0000000000e-01 0.0000000000e+00 + 53 -2.6787840270e-01 6.4671567270e-01 0.0000000000e+00 + 54 2.1424131440e-01 7.5845619390e-01 2.1239769560e-01 + 55 5.9844766760e-01 7.2778286530e-01 -2.9438622350e-01 + 56 8.1947984150e-01 7.0316915120e-01 5.7579440000e-04 + 57 -9.7705908600e-01 4.2033683210e-01 5.0830205500e-02 + 58 1.8982901160e-01 8.0989941880e-01 -2.4844539150e-01 + 59 -1.8982901160e-01 8.0989941880e-01 -2.4844539150e-01 + 60 0.0000000000e+00 8.5400868260e-01 -2.6208116160e-01 + 61 -5.9451875670e-01 7.3099586720e-01 -2.9438622350e-01 + 62 -1.0637993110e-01 9.0650575860e-01 5.8233869500e-02 + 63 5.4568911840e-01 9.1850197500e-01 4.6675922300e-02 + 64 -8.1947977940e-01 7.0316920800e-01 5.7577770000e-04 + 65 -5.1255932080e-01 9.2964041530e-01 8.4175679600e-02 + 66 1.9946137300e-01 1.0790266848e+00 -4.0864338300e-02 + 67 -2.9808074840e-01 1.1021175054e+00 -7.0205989000e-03 + 68 0.0000000000e+00 1.1568605599e+00 -2.5572400110e-01 + 69 -4.8516237780e-01 1.0708726508e+00 -2.4320238390e-01 + 70 0.0000000000e+00 1.2174811220e+00 2.0664452950e-01 +End Nodes + +Begin Elements Element3D4N + 1 0 43 37 35 48 + 2 0 25 27 32 30 + 3 0 11 14 17 7 + 4 0 52 54 49 66 + 5 0 37 28 35 41 + 6 0 27 36 32 38 + 7 0 45 37 43 48 + 8 0 21 27 25 30 + 9 0 45 51 42 57 + 10 0 21 20 26 23 + 11 0 60 59 62 67 + 12 0 59 62 67 53 + 13 0 9 8 5 1 + 14 0 8 5 1 10 + 15 0 62 67 53 65 + 16 0 5 1 10 6 + 17 0 45 50 51 57 + 18 0 21 16 20 23 + 19 0 22 28 29 33 + 20 0 44 36 40 46 + 21 0 52 60 62 66 + 22 0 11 9 5 7 + 23 0 62 67 70 68 + 24 0 5 1 2 3 + 25 0 67 65 69 59 + 26 0 1 6 4 8 + 27 0 48 41 37 35 + 28 0 30 38 27 32 + 29 0 62 60 67 68 + 30 0 5 9 1 3 + 31 0 14 17 7 18 + 32 0 54 49 66 63 + 33 0 29 28 39 33 + 34 0 40 36 34 46 + 35 0 14 11 5 7 + 36 0 54 52 62 66 + 37 0 17 12 7 18 + 38 0 49 58 66 63 + 39 0 67 59 53 65 + 40 0 1 8 10 6 + 41 0 62 60 68 66 + 42 0 5 9 3 7 + 43 0 57 43 48 45 + 44 0 23 25 30 21 + 45 0 31 33 41 28 + 46 0 47 46 38 36 + 47 0 69 65 61 59 + 48 0 4 6 13 8 + 49 0 28 35 41 31 + 50 0 36 32 38 47 + 51 0 45 42 48 57 + 52 0 21 26 30 23 + 53 0 39 41 33 28 + 54 0 34 38 46 36 + 55 0 53 61 50 65 + 56 0 10 13 16 6 + 57 0 2 7 5 3 + 58 0 70 66 62 68 + 59 0 19 17 22 18 + 60 0 55 49 44 63 + 61 0 19 24 18 22 + 62 0 55 56 63 44 + 63 0 65 64 61 50 + 64 0 6 15 13 16 + 65 0 33 24 29 22 + 66 0 46 56 40 44 + 67 0 64 57 51 50 + 68 0 15 23 20 16 + 69 0 19 24 22 29 + 70 0 55 56 44 40 + 71 0 64 61 50 51 + 72 0 15 13 16 20 + 73 0 53 59 61 65 + 74 0 10 8 13 6 + 75 0 17 12 18 19 + 76 0 49 58 63 55 +End Elements \ No newline at end of file diff --git a/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_solutions/block_from_1d_to_3d_data_transfer.json b/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_solutions/block_from_1d_to_3d_data_transfer.json new file mode 100644 index 000000000000..d87d64087edd --- /dev/null +++ b/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_solutions/block_from_1d_to_3d_data_transfer.json @@ -0,0 +1,405 @@ +{ + "TIME": [ + 0.0 + ], + "NODE_1": { + "TEMPERATURE": [ + -5.000000000000016 + ] + }, + "NODE_2": { + "TEMPERATURE": [ + -5.000000000000016 + ] + }, + "NODE_3": { + "TEMPERATURE": [ + -4.744114511400073 + ] + }, + "NODE_4": { + "TEMPERATURE": [ + -4.632766952099969 + ] + }, + "NODE_5": { + "TEMPERATURE": [ + -4.8753320696999936 + ] + }, + "NODE_6": { + "TEMPERATURE": [ + -5.000000000000016 + ] + }, + "NODE_7": { + "TEMPERATURE": [ + -5.000000000000016 + ] + }, + "NODE_8": { + "TEMPERATURE": [ + -4.421271446899998 + ] + }, + "NODE_9": { + "TEMPERATURE": [ + -4.291196363199978 + ] + }, + "NODE_10": { + "TEMPERATURE": [ + -4.339458637699915 + ] + }, + "NODE_11": { + "TEMPERATURE": [ + -4.378483632599995 + ] + }, + "NODE_12": { + "TEMPERATURE": [ + -4.802060292300001 + ] + }, + "NODE_13": { + "TEMPERATURE": [ + -4.564931462600029 + ] + }, + "NODE_14": { + "TEMPERATURE": [ + -4.460684542800019 + ] + }, + "NODE_15": { + "TEMPERATURE": [ + -4.084322500200041 + ] + }, + "NODE_16": { + "TEMPERATURE": [ + -4.2068619867999315 + ] + }, + "NODE_17": { + "TEMPERATURE": [ + -4.028861061700022 + ] + }, + "NODE_18": { + "TEMPERATURE": [ + -4.130189538300072 + ] + }, + "NODE_19": { + "TEMPERATURE": [ + -4.1176470587999905 + ] + }, + "NODE_20": { + "TEMPERATURE": [ + -3.9447709967999094 + ] + }, + "NODE_21": { + "TEMPERATURE": [ + -3.803778712400012 + ] + }, + "NODE_22": { + "TEMPERATURE": [ + -3.8450816702999906 + ] + }, + "NODE_23": { + "TEMPERATURE": [ + -3.629007628499963 + ] + }, + "NODE_24": { + "TEMPERATURE": [ + -3.674604974400008 + ] + }, + "NODE_25": { + "TEMPERATURE": [ + -3.6074842590999663 + ] + }, + "NODE_26": { + "TEMPERATURE": [ + -3.354974470999953 + ] + }, + "NODE_27": { + "TEMPERATURE": [ + -3.4562997361999694 + ] + }, + "NODE_28": { + "TEMPERATURE": [ + -3.366254742499978 + ] + }, + "NODE_29": { + "TEMPERATURE": [ + -3.1500500928000235 + ] + }, + "NODE_30": { + "TEMPERATURE": [ + -3.188123192000007 + ] + }, + "NODE_31": { + "TEMPERATURE": [ + -3.0882352940999827 + ] + }, + "NODE_32": { + "TEMPERATURE": [ + -2.949907332800022 + ] + }, + "NODE_33": { + "TEMPERATURE": [ + -3.0882352940999827 + ] + }, + "NODE_34": { + "TEMPERATURE": [ + -2.9449284529000144 + ] + }, + "NODE_35": { + "TEMPERATURE": [ + -2.8701266479999834 + ] + }, + "NODE_36": { + "TEMPERATURE": [ + -2.798603369199957 + ] + }, + "NODE_37": { + "TEMPERATURE": [ + -2.7955905052000016 + ] + }, + "NODE_38": { + "TEMPERATURE": [ + -2.7090064435999954 + ] + }, + "NODE_39": { + "TEMPERATURE": [ + -2.650428893000013 + ] + }, + "NODE_40": { + "TEMPERATURE": [ + -2.519972225799978 + ] + }, + "NODE_41": { + "TEMPERATURE": [ + -2.619949821899999 + ] + }, + "NODE_42": { + "TEMPERATURE": [ + -2.4248122609999996 + ] + }, + "NODE_43": { + "TEMPERATURE": [ + -2.324939933800011 + ] + }, + "NODE_44": { + "TEMPERATURE": [ + -2.325036540399993 + ] + }, + "NODE_45": { + "TEMPERATURE": [ + -2.185928188099998 + ] + }, + "NODE_46": { + "TEMPERATURE": [ + -2.062291079600019 + ] + }, + "NODE_47": { + "TEMPERATURE": [ + -2.0993851876000025 + ] + }, + "NODE_48": { + "TEMPERATURE": [ + -1.9251963758999755 + ] + }, + "NODE_49": { + "TEMPERATURE": [ + -1.8476097829999927 + ] + }, + "NODE_50": { + "TEMPERATURE": [ + -1.7928026980999978 + ] + }, + "NODE_51": { + "TEMPERATURE": [ + -1.746725007799981 + ] + }, + "NODE_52": { + "TEMPERATURE": [ + -1.6217611272999952 + ] + }, + "NODE_53": { + "TEMPERATURE": [ + -1.5957714246999872 + ] + }, + "NODE_54": { + "TEMPERATURE": [ + -1.5383051286999883 + ] + }, + "NODE_55": { + "TEMPERATURE": [ + -1.515820087999995 + ] + }, + "NODE_56": { + "TEMPERATURE": [ + -1.409387471 + ] + }, + "NODE_57": { + "TEMPERATURE": [ + -1.347134889900002 + ] + }, + "NODE_58": { + "TEMPERATURE": [ + -1.2601791011999965 + ] + }, + "NODE_59": { + "TEMPERATURE": [ + -1.1128282094999995 + ] + }, + "NODE_60": { + "TEMPERATURE": [ + -1.1272310391999967 + ] + }, + "NODE_61": { + "TEMPERATURE": [ + -1.1142290547000029 + ] + }, + "NODE_62": { + "TEMPERATURE": [ + -1.013403710499997 + ] + }, + "NODE_63": { + "TEMPERATURE": [ + -0.9300817773999981 + ] + }, + "NODE_64": { + "TEMPERATURE": [ + -0.8758774109999918 + ] + }, + "NODE_65": { + "TEMPERATURE": [ + -0.868314934999993 + ] + }, + "NODE_66": { + "TEMPERATURE": [ + -0.7461992327999923 + ] + }, + "NODE_67": { + "TEMPERATURE": [ + -0.6945563736999933 + ] + }, + "NODE_68": { + "TEMPERATURE": [ + -0.6560439909000009 + ] + }, + "NODE_69": { + "TEMPERATURE": [ + -0.6084473629999976 + ] + }, + "NODE_70": { + "TEMPERATURE": [ + -0.5941082989999897 + ] + }, + "NODE_71": { + "TEMPERATURE": [ + -0.45859261559999875 + ] + }, + "NODE_72": { + "TEMPERATURE": [ + -0.37149917809999305 + ] + }, + "NODE_73": { + "TEMPERATURE": [ + -0.3910223169999885 + ] + }, + "NODE_74": { + "TEMPERATURE": [ + -0.20463118699999838 + ] + }, + "NODE_75": { + "TEMPERATURE": [ + -0.17408453329999468 + ] + }, + "NODE_76": { + "TEMPERATURE": [ + -0.12089350139999472 + ] + }, + "NODE_77": { + "TEMPERATURE": [ + 0.0 + ] + }, + "NODE_78": { + "TEMPERATURE": [ + 0.0 + ] + }, + "NODE_79": { + "TEMPERATURE": [ + 0.0 + ] + }, + "NODE_80": { + "TEMPERATURE": [ + 0.0 + ] + } +} \ No newline at end of file diff --git a/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_solutions/block_from_3d_to_1d_data_transfer.json b/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_solutions/block_from_3d_to_1d_data_transfer.json new file mode 100644 index 000000000000..a373704da08a --- /dev/null +++ b/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_solutions/block_from_3d_to_1d_data_transfer.json @@ -0,0 +1,510 @@ +{ + "TIME": [ + 0.0 + ], + "NODE_1": { + "TEMPERATURE": [ + 1.1596359426047086e-16 + ] + }, + "NODE_2": { + "TEMPERATURE": [ + 0.05000000000000011 + ] + }, + "NODE_3": { + "TEMPERATURE": [ + 0.10000000000000012 + ] + }, + "NODE_4": { + "TEMPERATURE": [ + 0.14999999999999997 + ] + }, + "NODE_5": { + "TEMPERATURE": [ + 0.19999999999999982 + ] + }, + "NODE_6": { + "TEMPERATURE": [ + 0.24999999999999972 + ] + }, + "NODE_7": { + "TEMPERATURE": [ + 0.2999999999999997 + ] + }, + "NODE_8": { + "TEMPERATURE": [ + 0.34999999999999964 + ] + }, + "NODE_9": { + "TEMPERATURE": [ + 0.3999999999999997 + ] + }, + "NODE_10": { + "TEMPERATURE": [ + 0.4500000000000002 + ] + }, + "NODE_11": { + "TEMPERATURE": [ + 0.49999999999999983 + ] + }, + "NODE_12": { + "TEMPERATURE": [ + 0.5499999999999998 + ] + }, + "NODE_13": { + "TEMPERATURE": [ + 0.6000000000000002 + ] + }, + "NODE_14": { + "TEMPERATURE": [ + 0.6499999999999997 + ] + }, + "NODE_15": { + "TEMPERATURE": [ + 0.7000000000000004 + ] + }, + "NODE_16": { + "TEMPERATURE": [ + 0.75 + ] + }, + "NODE_17": { + "TEMPERATURE": [ + 0.8 + ] + }, + "NODE_18": { + "TEMPERATURE": [ + 0.85 + ] + }, + "NODE_19": { + "TEMPERATURE": [ + 0.9000000000000001 + ] + }, + "NODE_20": { + "TEMPERATURE": [ + 0.9500000000000008 + ] + }, + "NODE_21": { + "TEMPERATURE": [ + 0.9999999999999993 + ] + }, + "NODE_22": { + "TEMPERATURE": [ + 1.0499999999999994 + ] + }, + "NODE_23": { + "TEMPERATURE": [ + 1.0999999999999994 + ] + }, + "NODE_24": { + "TEMPERATURE": [ + 1.1500000000000004 + ] + }, + "NODE_25": { + "TEMPERATURE": [ + 1.2000000000000013 + ] + }, + "NODE_26": { + "TEMPERATURE": [ + 1.2499999999999996 + ] + }, + "NODE_27": { + "TEMPERATURE": [ + 1.2999999999999996 + ] + }, + "NODE_28": { + "TEMPERATURE": [ + 1.3499999999999999 + ] + }, + "NODE_29": { + "TEMPERATURE": [ + 1.3999999999999997 + ] + }, + "NODE_30": { + "TEMPERATURE": [ + 1.4499999999999997 + ] + }, + "NODE_31": { + "TEMPERATURE": [ + 1.4999999999999982 + ] + }, + "NODE_32": { + "TEMPERATURE": [ + 1.5499999999999985 + ] + }, + "NODE_33": { + "TEMPERATURE": [ + 1.5999999999999988 + ] + }, + "NODE_34": { + "TEMPERATURE": [ + 1.6499999999999984 + ] + }, + "NODE_35": { + "TEMPERATURE": [ + 1.7000000000000006 + ] + }, + "NODE_36": { + "TEMPERATURE": [ + 1.7500000000000018 + ] + }, + "NODE_37": { + "TEMPERATURE": [ + 1.7999999999999996 + ] + }, + "NODE_38": { + "TEMPERATURE": [ + 1.8499999999999994 + ] + }, + "NODE_39": { + "TEMPERATURE": [ + 1.9000000000000006 + ] + }, + "NODE_40": { + "TEMPERATURE": [ + 1.95 + ] + }, + "NODE_41": { + "TEMPERATURE": [ + 2.0 + ] + }, + "NODE_42": { + "TEMPERATURE": [ + 2.0499999999999994 + ] + }, + "NODE_43": { + "TEMPERATURE": [ + 2.099999999999999 + ] + }, + "NODE_44": { + "TEMPERATURE": [ + 2.1499999999999995 + ] + }, + "NODE_45": { + "TEMPERATURE": [ + 2.2 + ] + }, + "NODE_46": { + "TEMPERATURE": [ + 2.249999999999999 + ] + }, + "NODE_47": { + "TEMPERATURE": [ + 2.2999999999999985 + ] + }, + "NODE_48": { + "TEMPERATURE": [ + 2.3500000000000005 + ] + }, + "NODE_49": { + "TEMPERATURE": [ + 2.4 + ] + }, + "NODE_50": { + "TEMPERATURE": [ + 2.4499999999999993 + ] + }, + "NODE_51": { + "TEMPERATURE": [ + 2.4999999999999987 + ] + }, + "NODE_52": { + "TEMPERATURE": [ + 2.549999999999999 + ] + }, + "NODE_53": { + "TEMPERATURE": [ + 2.5999999999999996 + ] + }, + "NODE_54": { + "TEMPERATURE": [ + 2.65 + ] + }, + "NODE_55": { + "TEMPERATURE": [ + 2.6999999999999997 + ] + }, + "NODE_56": { + "TEMPERATURE": [ + 2.75 + ] + }, + "NODE_57": { + "TEMPERATURE": [ + 2.8 + ] + }, + "NODE_58": { + "TEMPERATURE": [ + 2.849999999999998 + ] + }, + "NODE_59": { + "TEMPERATURE": [ + 2.9000000000000017 + ] + }, + "NODE_60": { + "TEMPERATURE": [ + 2.9500000000000033 + ] + }, + "NODE_61": { + "TEMPERATURE": [ + 3.0000000000000004 + ] + }, + "NODE_62": { + "TEMPERATURE": [ + 3.050000000000001 + ] + }, + "NODE_63": { + "TEMPERATURE": [ + 3.1000000000000005 + ] + }, + "NODE_64": { + "TEMPERATURE": [ + 3.1500000000000012 + ] + }, + "NODE_65": { + "TEMPERATURE": [ + 3.2 + ] + }, + "NODE_66": { + "TEMPERATURE": [ + 3.25 + ] + }, + "NODE_67": { + "TEMPERATURE": [ + 3.3 + ] + }, + "NODE_68": { + "TEMPERATURE": [ + 3.350000000000002 + ] + }, + "NODE_69": { + "TEMPERATURE": [ + 3.4000000000000017 + ] + }, + "NODE_70": { + "TEMPERATURE": [ + 3.4500000000000006 + ] + }, + "NODE_71": { + "TEMPERATURE": [ + 3.5000000000000004 + ] + }, + "NODE_72": { + "TEMPERATURE": [ + 3.5500000000000007 + ] + }, + "NODE_73": { + "TEMPERATURE": [ + 3.6000000000000005 + ] + }, + "NODE_74": { + "TEMPERATURE": [ + 3.6499999999999995 + ] + }, + "NODE_75": { + "TEMPERATURE": [ + 3.6999999999999993 + ] + }, + "NODE_76": { + "TEMPERATURE": [ + 3.7499999999999973 + ] + }, + "NODE_77": { + "TEMPERATURE": [ + 3.7999999999999976 + ] + }, + "NODE_78": { + "TEMPERATURE": [ + 3.8499999999999983 + ] + }, + "NODE_79": { + "TEMPERATURE": [ + 3.8999999999999986 + ] + }, + "NODE_80": { + "TEMPERATURE": [ + 3.9499999999999997 + ] + }, + "NODE_81": { + "TEMPERATURE": [ + 4.000000000000001 + ] + }, + "NODE_82": { + "TEMPERATURE": [ + 4.05 + ] + }, + "NODE_83": { + "TEMPERATURE": [ + 4.1000000000000005 + ] + }, + "NODE_84": { + "TEMPERATURE": [ + 4.150000000000002 + ] + }, + "NODE_85": { + "TEMPERATURE": [ + 4.200000000000002 + ] + }, + "NODE_86": { + "TEMPERATURE": [ + 4.250000000000002 + ] + }, + "NODE_87": { + "TEMPERATURE": [ + 4.3000000000000025 + ] + }, + "NODE_88": { + "TEMPERATURE": [ + 4.349999999999999 + ] + }, + "NODE_89": { + "TEMPERATURE": [ + 4.4 + ] + }, + "NODE_90": { + "TEMPERATURE": [ + 4.45 + ] + }, + "NODE_91": { + "TEMPERATURE": [ + 4.500000000000001 + ] + }, + "NODE_92": { + "TEMPERATURE": [ + 4.549999999999997 + ] + }, + "NODE_93": { + "TEMPERATURE": [ + 4.599999999999996 + ] + }, + "NODE_94": { + "TEMPERATURE": [ + 4.650000000000001 + ] + }, + "NODE_95": { + "TEMPERATURE": [ + 4.699999999999999 + ] + }, + "NODE_96": { + "TEMPERATURE": [ + 4.750000000000003 + ] + }, + "NODE_97": { + "TEMPERATURE": [ + 4.8000000000000025 + ] + }, + "NODE_98": { + "TEMPERATURE": [ + 4.850000000000001 + ] + }, + "NODE_99": { + "TEMPERATURE": [ + 4.900000000000002 + ] + }, + "NODE_100": { + "TEMPERATURE": [ + 4.949999999999999 + ] + }, + "NODE_101": { + "TEMPERATURE": [ + 5.0 + ] + } + } \ No newline at end of file diff --git a/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_solutions/torus_from_1d_to_3d_data_transfer.json b/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_solutions/torus_from_1d_to_3d_data_transfer.json new file mode 100644 index 000000000000..f753c8fc05f0 --- /dev/null +++ b/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_solutions/torus_from_1d_to_3d_data_transfer.json @@ -0,0 +1,355 @@ +{ + "TIME": [ + 0.0 + ], + "NODE_1": { + "TEMPERATURE": [ + 1.2242051915147831 + ] + }, + "NODE_2": { + "TEMPERATURE": [ + 1.0070424509768041 + ] + }, + "NODE_3": { + "TEMPERATURE": [ + 0.9952300862169075 + ] + }, + "NODE_4": { + "TEMPERATURE": [ + 1.324647709945493 + ] + }, + "NODE_5": { + "TEMPERATURE": [ + 1.1077008105716535 + ] + }, + "NODE_6": { + "TEMPERATURE": [ + 1.397961686300839 + ] + }, + "NODE_7": { + "TEMPERATURE": [ + 0.7985694731629032 + ] + }, + "NODE_8": { + "TEMPERATURE": [ + 1.1999749511432667 + ] + }, + "NODE_9": { + "TEMPERATURE": [ + 1.004439395074985 + ] + }, + "NODE_10": { + "TEMPERATURE": [ + 1.3085295567077497 + ] + }, + "NODE_11": { + "TEMPERATURE": [ + 0.9902854313348358 + ] + }, + "NODE_12": { + "TEMPERATURE": [ + 0.7471507135297519 + ] + }, + "NODE_13": { + "TEMPERATURE": [ + 1.406221240670888 + ] + }, + "NODE_14": { + "TEMPERATURE": [ + 0.6882384603612067 + ] + }, + "NODE_15": { + "TEMPERATURE": [ + 1.4094409639610816 + ] + }, + "NODE_16": { + "TEMPERATURE": [ + 1.4135157333501 + ] + }, + "NODE_17": { + "TEMPERATURE": [ + 0.5348073607089424 + ] + }, + "NODE_18": { + "TEMPERATURE": [ + 0.34859297543875467 + ] + }, + "NODE_19": { + "TEMPERATURE": [ + 0.13955398293077756 + ] + }, + "NODE_20": { + "TEMPERATURE": [ + 1.3419912864298753 + ] + }, + "NODE_21": { + "TEMPERATURE": [ + 1.3085295567077546 + ] + }, + "NODE_22": { + "TEMPERATURE": [ + -7.854827899222983e-15 + ] + }, + "NODE_23": { + "TEMPERATURE": [ + 1.3130714543045237 + ] + }, + "NODE_24": { + "TEMPERATURE": [ + -0.10918824437398825 + ] + }, + "NODE_25": { + "TEMPERATURE": [ + 1.1540391285180764 + ] + }, + "NODE_26": { + "TEMPERATURE": [ + 1.1356361535543054 + ] + }, + "NODE_27": { + "TEMPERATURE": [ + 1.0091225871936453 + ] + }, + "NODE_28": { + "TEMPERATURE": [ + -0.534807360708959 + ] + }, + "NODE_29": { + "TEMPERATURE": [ + -0.4440085111696036 + ] + }, + "NODE_30": { + "TEMPERATURE": [ + 1.0946073831013212 + ] + }, + "NODE_31": { + "TEMPERATURE": [ + -0.44327269438953026 + ] + }, + "NODE_32": { + "TEMPERATURE": [ + 0.8556151165694356 + ] + }, + "NODE_33": { + "TEMPERATURE": [ + -0.5233358368565266 + ] + }, + "NODE_34": { + "TEMPERATURE": [ + 0.8506041366239743 + ] + }, + "NODE_35": { + "TEMPERATURE": [ + -0.8556151165694376 + ] + }, + "NODE_36": { + "TEMPERATURE": [ + 0.5348073607089591 + ] + }, + "NODE_37": { + "TEMPERATURE": [ + -0.9902854313348459 + ] + }, + "NODE_38": { + "TEMPERATURE": [ + 0.8943909520281689 + ] + }, + "NODE_39": { + "TEMPERATURE": [ + -0.8506041366239765 + ] + }, + "NODE_40": { + "TEMPERATURE": [ + 0.4440085111695957 + ] + }, + "NODE_41": { + "TEMPERATURE": [ + -0.894390952028169 + ] + }, + "NODE_42": { + "TEMPERATURE": [ + -1.1356361535543038 + ] + }, + "NODE_43": { + "TEMPERATURE": [ + -1.1540391285180762 + ] + }, + "NODE_44": { + "TEMPERATURE": [ + 7.185224637495935e-15 + ] + }, + "NODE_45": { + "TEMPERATURE": [ + -1.3085295567077562 + ] + }, + "NODE_46": { + "TEMPERATURE": [ + 0.5233358368565268 + ] + }, + "NODE_47": { + "TEMPERATURE": [ + 0.4432726943895222 + ] + }, + "NODE_48": { + "TEMPERATURE": [ + -1.0946073831013226 + ] + }, + "NODE_49": { + "TEMPERATURE": [ + -0.5348073607089437 + ] + }, + "NODE_50": { + "TEMPERATURE": [ + -1.4135157333501 + ] + }, + "NODE_51": { + "TEMPERATURE": [ + -1.3419912864298753 + ] + }, + "NODE_52": { + "TEMPERATURE": [ + -0.990285431334835 + ] + }, + "NODE_53": { + "TEMPERATURE": [ + -1.30852955670775 + ] + }, + "NODE_54": { + "TEMPERATURE": [ + -0.6882384603612071 + ] + }, + "NODE_55": { + "TEMPERATURE": [ + -0.13955398293077223 + ] + }, + "NODE_56": { + "TEMPERATURE": [ + 0.10918824437398837 + ] + }, + "NODE_57": { + "TEMPERATURE": [ + -1.3130714543045245 + ] + }, + "NODE_58": { + "TEMPERATURE": [ + -0.7471507135297535 + ] + }, + "NODE_59": { + "TEMPERATURE": [ + -1.1999749511432671 + ] + }, + "NODE_60": { + "TEMPERATURE": [ + -0.9952725244086595 + ] + }, + "NODE_61": { + "TEMPERATURE": [ + -1.4062212406708874 + ] + }, + "NODE_62": { + "TEMPERATURE": [ + -1.107700810571653 + ] + }, + "NODE_63": { + "TEMPERATURE": [ + -0.34859297543875395 + ] + }, + "NODE_64": { + "TEMPERATURE": [ + -1.4094409639610816 + ] + }, + "NODE_65": { + "TEMPERATURE": [ + -1.39796168630084 + ] + }, + "NODE_66": { + "TEMPERATURE": [ + -0.7985694731629069 + ] + }, + "NODE_67": { + "TEMPERATURE": [ + -1.2242051915147831 + ] + }, + "NODE_68": { + "TEMPERATURE": [ + -1.0050794422666673 + ] + }, + "NODE_69": { + "TEMPERATURE": [ + -1.3246477099454892 + ] + }, + "NODE_70": { + "TEMPERATURE": [ + -1.0070424509768032 + ] + } +} \ No newline at end of file diff --git a/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_solutions/torus_from_3d_to_1d_data_transfer.json b/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_solutions/torus_from_3d_to_1d_data_transfer.json new file mode 100644 index 000000000000..6975a158e642 --- /dev/null +++ b/applications/CoSimulationApplication/tests/3d_1d_data_transfer/3d_1d_data_transfer_solutions/torus_from_3d_to_1d_data_transfer.json @@ -0,0 +1,505 @@ +{ + "TIME": [ + 0.0 + ], + "NODE_1": { + "TEMPERATURE": [ + -1.0000000000000002 + ] + }, + "NODE_2": { + "TEMPERATURE": [ + -0.9352362088989583 + ] + }, + "NODE_3": { + "TEMPERATURE": [ + -0.8667814677501736 + ] + }, + "NODE_4": { + "TEMPERATURE": [ + -0.7949059361429641 + ] + }, + "NODE_5": { + "TEMPERATURE": [ + -0.7198932739637763 + ] + }, + "NODE_6": { + "TEMPERATURE": [ + -0.6420395219202062 + ] + }, + "NODE_7": { + "TEMPERATURE": [ + -0.5616519332035734 + ] + }, + "NODE_8": { + "TEMPERATURE": [ + -0.47904776090094686 + ] + }, + "NODE_9": { + "TEMPERATURE": [ + -0.3945530059421484 + ] + }, + "NODE_10": { + "TEMPERATURE": [ + -0.30850113052301836 + ] + }, + "NODE_11": { + "TEMPERATURE": [ + -0.22123174208247426 + ] + }, + "NODE_12": { + "TEMPERATURE": [ + -0.13308925302709967 + ] + }, + "NODE_13": { + "TEMPERATURE": [ + -0.04442152149272273 + ] + }, + "NODE_14": { + "TEMPERATURE": [ + 0.04442152149272308 + ] + }, + "NODE_15": { + "TEMPERATURE": [ + 0.13308925302709967 + ] + }, + "NODE_16": { + "TEMPERATURE": [ + 0.2212317420824742 + ] + }, + "NODE_17": { + "TEMPERATURE": [ + 0.30850113052301875 + ] + }, + "NODE_18": { + "TEMPERATURE": [ + 0.3945530059421486 + ] + }, + "NODE_19": { + "TEMPERATURE": [ + 0.4790477609009469 + ] + }, + "NODE_20": { + "TEMPERATURE": [ + 0.5616519332035732 + ] + }, + "NODE_21": { + "TEMPERATURE": [ + 0.6420395219202063 + ] + }, + "NODE_22": { + "TEMPERATURE": [ + 0.7198932739637764 + ] + }, + "NODE_23": { + "TEMPERATURE": [ + 0.794905936142964 + ] + }, + "NODE_24": { + "TEMPERATURE": [ + 0.8667814677501733 + ] + }, + "NODE_25": { + "TEMPERATURE": [ + 0.9352362088989581 + ] + }, + "NODE_26": { + "TEMPERATURE": [ + 0.9999999999999999 + ] + }, + "NODE_27": { + "TEMPERATURE": [ + 1.0608172479575848 + ] + }, + "NODE_28": { + "TEMPERATURE": [ + 1.117447934878782 + ] + }, + "NODE_29": { + "TEMPERATURE": [ + 1.169668565314413 + ] + }, + "NODE_30": { + "TEMPERATURE": [ + 1.2172730482934861 + ] + }, + "NODE_31": { + "TEMPERATURE": [ + 1.2600735106701006 + ] + }, + "NODE_32": { + "TEMPERATURE": [ + 1.297901038572929 + ] + }, + "NODE_33": { + "TEMPERATURE": [ + 1.3306063440310922 + ] + }, + "NODE_34": { + "TEMPERATURE": [ + 1.3580603541455791 + ] + }, + "NODE_35": { + "TEMPERATURE": [ + 1.3801547204810118 + ] + }, + "NODE_36": { + "TEMPERATURE": [ + 1.3968022466674208 + ] + }, + "NODE_37": { + "TEMPERATURE": [ + 1.4079372325244792 + ] + }, + "NODE_38": { + "TEMPERATURE": [ + 1.4135157333501005 + ] + }, + "NODE_39": { + "TEMPERATURE": [ + 1.4135157333501003 + ] + }, + "NODE_40": { + "TEMPERATURE": [ + 1.407937232524479 + ] + }, + "NODE_41": { + "TEMPERATURE": [ + 1.3968022466674204 + ] + }, + "NODE_42": { + "TEMPERATURE": [ + 1.3801547204810116 + ] + }, + "NODE_43": { + "TEMPERATURE": [ + 1.3580603541455787 + ] + }, + "NODE_44": { + "TEMPERATURE": [ + 1.330606344031092 + ] + }, + "NODE_45": { + "TEMPERATURE": [ + 1.2979010385729297 + ] + }, + "NODE_46": { + "TEMPERATURE": [ + 1.260073510670101 + ] + }, + "NODE_47": { + "TEMPERATURE": [ + 1.2172730482934864 + ] + }, + "NODE_48": { + "TEMPERATURE": [ + 1.1696685653144134 + ] + }, + "NODE_49": { + "TEMPERATURE": [ + 1.1174479348787822 + ] + }, + "NODE_50": { + "TEMPERATURE": [ + 1.0608172479575853 + ] + }, + "NODE_51": { + "TEMPERATURE": [ + 1.0000000000000004 + ] + }, + "NODE_52": { + "TEMPERATURE": [ + 0.9352362088989583 + ] + }, + "NODE_53": { + "TEMPERATURE": [ + 0.8667814677501736 + ] + }, + "NODE_54": { + "TEMPERATURE": [ + 0.7949059361429641 + ] + }, + "NODE_55": { + "TEMPERATURE": [ + 0.7198932739637767 + ] + }, + "NODE_56": { + "TEMPERATURE": [ + 0.6420395219202065 + ] + }, + "NODE_57": { + "TEMPERATURE": [ + 0.5616519332035734 + ] + }, + "NODE_58": { + "TEMPERATURE": [ + 0.47904776090094753 + ] + }, + "NODE_59": { + "TEMPERATURE": [ + 0.3945530059421483 + ] + }, + "NODE_60": { + "TEMPERATURE": [ + 0.30850113052301886 + ] + }, + "NODE_61": { + "TEMPERATURE": [ + 0.22123174208247512 + ] + }, + "NODE_62": { + "TEMPERATURE": [ + 0.13308925302709967 + ] + }, + "NODE_63": { + "TEMPERATURE": [ + 0.044421521492723375 + ] + }, + "NODE_64": { + "TEMPERATURE": [ + -0.044421521492723076 + ] + }, + "NODE_65": { + "TEMPERATURE": [ + -0.13308925302709995 + ] + }, + "NODE_66": { + "TEMPERATURE": [ + -0.22123174208247423 + ] + }, + "NODE_67": { + "TEMPERATURE": [ + -0.3085011305230191 + ] + }, + "NODE_68": { + "TEMPERATURE": [ + -0.39455300594214854 + ] + }, + "NODE_69": { + "TEMPERATURE": [ + -0.47904776090094736 + ] + }, + "NODE_70": { + "TEMPERATURE": [ + -0.5616519332035737 + ] + }, + "NODE_71": { + "TEMPERATURE": [ + -0.6420395219202061 + ] + }, + "NODE_72": { + "TEMPERATURE": [ + -0.719893273963777 + ] + }, + "NODE_73": { + "TEMPERATURE": [ + -0.7949059361429641 + ] + }, + "NODE_74": { + "TEMPERATURE": [ + -0.8667814677501731 + ] + }, + "NODE_75": { + "TEMPERATURE": [ + -0.9352362088989584 + ] + }, + "NODE_76": { + "TEMPERATURE": [ + -0.9999999999999998 + ] + }, + "NODE_77": { + "TEMPERATURE": [ + -1.0608172479575844 + ] + }, + "NODE_78": { + "TEMPERATURE": [ + -1.1174479348787822 + ] + }, + "NODE_79": { + "TEMPERATURE": [ + -1.1696685653144125 + ] + }, + "NODE_80": { + "TEMPERATURE": [ + -1.2172730482934861 + ] + }, + "NODE_81": { + "TEMPERATURE": [ + -1.260073510670101 + ] + }, + "NODE_82": { + "TEMPERATURE": [ + -1.297901038572929 + ] + }, + "NODE_83": { + "TEMPERATURE": [ + -1.3306063440310916 + ] + }, + "NODE_84": { + "TEMPERATURE": [ + -1.3580603541455796 + ] + }, + "NODE_85": { + "TEMPERATURE": [ + -1.3801547204810116 + ] + }, + "NODE_86": { + "TEMPERATURE": [ + -1.3968022466674204 + ] + }, + "NODE_87": { + "TEMPERATURE": [ + -1.407937232524479 + ] + }, + "NODE_88": { + "TEMPERATURE": [ + -1.4135157333501 + ] + }, + "NODE_89": { + "TEMPERATURE": [ + -1.4135157333501 + ] + }, + "NODE_90": { + "TEMPERATURE": [ + -1.4079372325244788 + ] + }, + "NODE_91": { + "TEMPERATURE": [ + -1.3968022466674204 + ] + }, + "NODE_92": { + "TEMPERATURE": [ + -1.3801547204810118 + ] + }, + "NODE_93": { + "TEMPERATURE": [ + -1.3580603541455791 + ] + }, + "NODE_94": { + "TEMPERATURE": [ + -1.330606344031092 + ] + }, + "NODE_95": { + "TEMPERATURE": [ + -1.2979010385729293 + ] + }, + "NODE_96": { + "TEMPERATURE": [ + -1.2600735106701013 + ] + }, + "NODE_97": { + "TEMPERATURE": [ + -1.2172730482934864 + ] + }, + "NODE_98": { + "TEMPERATURE": [ + -1.1696685653144134 + ] + }, + "NODE_99": { + "TEMPERATURE": [ + -1.1174479348787825 + ] + }, + "NODE_100": { + "TEMPERATURE": [ + -1.060817247957585 + ] + } + } \ No newline at end of file diff --git a/applications/CoSimulationApplication/tests/test_3d_1d_data_transfer_process.py b/applications/CoSimulationApplication/tests/test_3d_1d_data_transfer_process.py new file mode 100644 index 000000000000..ded91e5e6a30 --- /dev/null +++ b/applications/CoSimulationApplication/tests/test_3d_1d_data_transfer_process.py @@ -0,0 +1,279 @@ +# Import the necessary Kratos modules for multiphysics simulations and testing. +import KratosMultiphysics as KM +import KratosMultiphysics.KratosUnittest as KratosUnittest + +# Import specific applications for mapping and co-simulation. +import KratosMultiphysics.MappingApplication +import KratosMultiphysics.CoSimulationApplication as CoSimulationApplication + +# Import utilities for processing JSON results. +from KratosMultiphysics.from_json_check_result_process import FromJsonCheckResultProcess +from KratosMultiphysics.json_output_process import JsonOutputProcess + +# Import basic dependencies for math and file path handling. +import math +from pathlib import Path + +# Function to get the full file path for a given file name. +def GetFilePath(fileName): + return Path(__file__).parent / fileName + +class Test3D1DDataTransferProcessBlock(KratosUnittest.TestCase): + """Test class for the 3D to 1D data transfer process.""" + + def setUp(self): + """Setup function to initialize model parts for testing.""" + self.current_model = KM.Model() + + # Initialize a 3D model part called "Block". + self.model_part_block = self.current_model.CreateModelPart("Block") + self.model_part_block.ProcessInfo[KM.STEP] = 0 + self.model_part_block.ProcessInfo.SetValue(KM.DOMAIN_SIZE, 3) + self.model_part_block.ProcessInfo.SetValue(KM.TIME, 0.0) + self.model_part_block.ProcessInfo.SetValue(KM.DELTA_TIME, 1.0) + self.model_part_block.AddNodalSolutionStepVariable(KM.TEMPERATURE) + input_mdpa = GetFilePath("3d_1d_data_transfer/3d_1d_data_transfer_mdpa/block") + model_part_io_block = KM.ModelPartIO(input_mdpa) + model_part_io_block.ReadModelPart(self.model_part_block) + + # Initialize a 1D model part called "Line". + self.model_part_line = self.current_model.CreateModelPart("Line") + self.model_part_line.ProcessInfo[KM.STEP] = 0 + self.model_part_line.ProcessInfo.SetValue(KM.DOMAIN_SIZE, 3) + self.model_part_line.ProcessInfo.SetValue(KM.TIME, 0.0) + self.model_part_line.ProcessInfo.SetValue(KM.DELTA_TIME, 1.0) + self.model_part_line.AddNodalSolutionStepVariable(KM.TEMPERATURE) + + # Create nodes and elements in the 1D model part. + number_of_line_elements = 100 + initial_z_coordinate = 0.0 + end_z_coordinate = 5.0 + self.model_part_line.CreateNewNode(1, 0.5, 0.5, initial_z_coordinate) + prop = self.model_part_line.GetProperties()[1] + for i in range(number_of_line_elements): + z = initial_z_coordinate + (i + 1) * (end_z_coordinate - initial_z_coordinate) / number_of_line_elements + self.model_part_line.CreateNewNode(i + 2, 0.5, 0.5, z) + self.model_part_line.CreateNewElement("Element3D2N", i + 1, [i + 1, i + 2], prop) + + # Set temperature values for nodes in the block and line model parts. + for node in self.model_part_block.Nodes: + node.SetSolutionStepValue(KM.TEMPERATURE, node.Z) # Set temperature based on Z coordinate. + del node + + for node in self.model_part_line.Nodes: + node.SetSolutionStepValue(KM.TEMPERATURE, -node.Z) # Set negative temperature based on Z coordinate. + del node + + def test_block_from_1d_to_3d(self): + """Test data transfer from the 1D line to the 3D block.""" + # Define parameters for the data transfer process. + parameters = KM.Parameters("""{ + "origin_variables" : ["TEMPERATURE"], + "destination_variables" : ["TEMPERATURE"] + }""") + + # Execute the data transfer process. + process = CoSimulationApplication.DataTransfer3D1DProcess(self.model_part_block, self.model_part_line, parameters) + process.Execute() + + # Define result file path for the output of the transfer. + result_file = GetFilePath("3d_1d_data_transfer/3d_1d_data_transfer_solutions/block_from_1d_to_3d") + + # # Generate output results in JSON format. + # generate_result(self.current_model, result_file, "Block") + + # Check the results against expected outcomes. + check_results(self.current_model, result_file, "Block") + + # Debug output for visualization. + # debug_vtk(self.current_model, ["Block", "Line"]) + + def test_block_from_3d_to_1d(self): + """Test data transfer from the 3D block to the 1D line.""" + # Define parameters for the data transfer process. + parameters = KM.Parameters("""{ + "origin_variables" : ["TEMPERATURE"], + "destination_variables" : ["TEMPERATURE"] + }""") + + # Execute the data transfer process with transpose setting. + process = CoSimulationApplication.DataTransfer3D1DProcess(self.model_part_block, self.model_part_line, parameters) + process.Set(KratosMultiphysics.Mapper.USE_TRANSPOSE) + process.Execute() + + # Define result file path for the output of the transfer. + result_file = GetFilePath("3d_1d_data_transfer/3d_1d_data_transfer_solutions/block_from_3d_to_1d") + + # # Generate output results in JSON format. + # generate_result(self.current_model, result_file, "Line") + + # Check the results against expected outcomes. + check_results(self.current_model, result_file, "Line") + + # Debug output for visualization. + # debug_vtk(self.current_model, ["Block", "Line"]) + +# Test class for the toroidal data transfer process. +class Test3D1DDataTransferProcessTorus(KratosUnittest.TestCase): + def setUp(self): + """Setup function to initialize model parts for testing with a torus.""" + self.current_model = KM.Model() + + # Initialize a 3D model part called "Torus". + self.model_part_torus = self.current_model.CreateModelPart("Torus") + self.model_part_torus.ProcessInfo[KM.STEP] = 0 + self.model_part_torus.ProcessInfo.SetValue(KM.DOMAIN_SIZE, 3) + self.model_part_torus.ProcessInfo.SetValue(KM.TIME, 0.0) + self.model_part_torus.ProcessInfo.SetValue(KM.DELTA_TIME, 1.0) + self.model_part_torus.AddNodalSolutionStepVariable(KM.TEMPERATURE) + input_mdpa = GetFilePath("3d_1d_data_transfer/3d_1d_data_transfer_mdpa/torus3d") + model_part_io_torus = KM.ModelPartIO(input_mdpa) + model_part_io_torus.ReadModelPart(self.model_part_torus) + + # Initialize a 1D model part called "Circle". + self.model_part_circle = self.current_model.CreateModelPart("Circle") + self.model_part_circle.ProcessInfo[KM.STEP] = 0 + self.model_part_circle.ProcessInfo.SetValue(KM.DOMAIN_SIZE, 3) + self.model_part_circle.ProcessInfo.SetValue(KM.TIME, 0.0) + self.model_part_circle.ProcessInfo.SetValue(KM.DELTA_TIME, 1.0) + self.model_part_circle.AddNodalSolutionStepVariable(KM.TEMPERATURE) + + # Create nodes and elements in the circular model part. + number_of_circle_elements = 100 + radius = 1.0 + self.model_part_circle.CreateNewNode(1, radius, 0.0, 0) + prop = self.model_part_circle.GetProperties()[1] + for i in range(number_of_circle_elements - 1): + angle = 2 * math.pi * (i + 1) / number_of_circle_elements + self.model_part_circle.CreateNewNode(i + 2, radius * math.cos(angle), radius * math.sin(angle), 0) + self.model_part_circle.CreateNewElement("Element3D2N", i + 1, [i + 1, i + 2], prop) + self.model_part_circle.CreateNewElement("Element3D2N", number_of_circle_elements, [i + 2, 1], prop) + + # Set temperature values for nodes in the torus and circle model parts. + for node in self.model_part_torus.Nodes: + node.SetSolutionStepValue(KM.TEMPERATURE, -node.X + node.Y) # Example temperature calculation. + del node + + for node in self.model_part_circle.Nodes: + node.SetSolutionStepValue(KM.TEMPERATURE, node.X - node.Y) # Example temperature calculation. + del node + + def test_torus_from_1d_to_3d(self): + """Test data transfer from the 1D circle to the 3D torus.""" + # Define parameters for the data transfer process. + parameters = KM.Parameters("""{ + "origin_variables" : ["TEMPERATURE"], + "destination_variables" : ["TEMPERATURE"] + }""") + + # Execute the data transfer process. + process = CoSimulationApplication.DataTransfer3D1DProcess(self.model_part_torus, self.model_part_circle, parameters) + process.Execute() + + # Define result file path for the output of the transfer. + result_file = GetFilePath("3d_1d_data_transfer/3d_1d_data_transfer_solutions/torus_from_1d_to_3d") + + # # Generate output results in JSON format. + # generate_result(self.current_model, result_file, "Torus") + + # Check the results against expected outcomes. + check_results(self.current_model, result_file, "Torus") + + # Debug output for visualization. + # debug_vtk(self.current_model, ["Torus", "Circle"]) + + def test_torus_from_3d_to_1d(self): + """Test data transfer from the 3D torus to the 1D circle.""" + # Define parameters for the data transfer process. + parameters = KM.Parameters("""{ + "origin_variables" : ["TEMPERATURE"], + "destination_variables" : ["TEMPERATURE"] + }""") + + # Execute the data transfer process with transpose setting. + process = CoSimulationApplication.DataTransfer3D1DProcess(self.model_part_torus, self.model_part_circle, parameters) + process.Set(KratosMultiphysics.Mapper.USE_TRANSPOSE) + process.Execute() + + # Define result file path for the output of the transfer. + result_file = GetFilePath("3d_1d_data_transfer/3d_1d_data_transfer_solutions/torus_from_3d_to_1d") + + # # Generate output results in JSON format. + # generate_result(self.current_model, result_file, "Circle") + + # Check the results against expected outcomes. + check_results(self.current_model, result_file, "Circle") + + # Debug output for visualization. + # debug_vtk(self.current_model, ["Torus", "Circle"]) + +def debug_vtk(model, list_names): + """Generate VTK output for visualization of the model parts.""" + import KratosMultiphysics.vtk_output_process as vtk_output_process + + for name in list_names: + # Define parameters for VTK output. + vtk_output_parameters = KratosMultiphysics.Parameters("""{ + "Parameters" : { + "model_part_name" : "", + "output_precision" : 8, + "output_interval" : 2, + "output_sub_model_parts" : true, + "output_path" : "", + "nodal_solution_step_data_variables" : ["TEMPERATURE"], + "element_flags" : [] + } + }""") + vtk_output_parameters["Parameters"]["model_part_name"].SetString(name) + vtk_output_parameters["Parameters"]["output_path"].SetString("vtk_output_" + name) + + # Initialize and execute VTK output process. + vtk_output_process_torus = vtk_output_process.Factory(vtk_output_parameters, model) + vtk_output_process_torus.ExecuteInitializeSolutionStep() + vtk_output_process_torus.ExecuteFinalizeSolutionStep() + vtk_output_process_torus.PrintOutput() + +def check_results(model, input_filename, domain): + """Check the results against expected output using JSON files.""" + # Define parameters for result checking. + check_parameters = KM.Parameters("""{ + "check_variables" : ["TEMPERATURE"], + "input_file_name" : "", + "model_part_name" : "", + "historical_value" : true, + "time_frequency" : 0.0 + }""") + + check_parameters["model_part_name"].SetString(domain) + check_parameters["input_file_name"].SetString(str(input_filename) + "_data_transfer.json") + + # Execute the checking process. + check = FromJsonCheckResultProcess(model, check_parameters) + check.ExecuteInitialize() + check.ExecuteBeforeSolutionLoop() + check.ExecuteFinalizeSolutionStep() + +def generate_result(model, output_filename, domain): + """Generate output results in JSON format.""" + # Define parameters for output generation. + out_parameters = KM.Parameters("""{ + "output_variables" : ["TEMPERATURE"], + "output_file_name" : "", + "model_part_name" : "", + "historical_value" : true, + "time_frequency" : 0.0 + }""") + + out_parameters["model_part_name"].SetString(domain) + out_parameters["output_file_name"].SetString(f"{output_filename}_data_transfer.json") + + # Execute the output generation process. + out = JsonOutputProcess(model, out_parameters) + out.ExecuteInitialize() + out.ExecuteBeforeSolutionLoop() + out.ExecuteFinalizeSolutionStep() + +# Main execution block to run the unit tests. +if __name__ == '__main__': + KM.Logger.GetDefaultOutput().SetSeverity(KM.Logger.Severity.WARNING) # Set the logging severity level. + KratosUnittest.main() # Execute the unit tests. \ No newline at end of file diff --git a/applications/CoSimulationApplication/tests/test_CoSimulationApplication.py b/applications/CoSimulationApplication/tests/test_CoSimulationApplication.py index ff9c17975150..d44643c7aed3 100644 --- a/applications/CoSimulationApplication/tests/test_CoSimulationApplication.py +++ b/applications/CoSimulationApplication/tests/test_CoSimulationApplication.py @@ -23,6 +23,8 @@ from test_co_simulation_coupled_solver import TestCoupledSolverCouplingInterfaceDataAccess from test_model_part_utilties import TestModelPartUtiliites from test_thermal_rom_co_sim import TestThermalRomCoSim +from test_3d_1d_data_transfer_process import Test3D1DDataTransferProcessBlock +from test_3d_1d_data_transfer_process import Test3D1DDataTransferProcessTorus from test_cosim_EMPIRE_API import TestCoSim_EMPIRE_API from test_co_sim_io_py_exposure import TestCoSimIOPyExposure @@ -65,6 +67,8 @@ def AssembleTestSuites(): smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TestConvergenceAcceleratorWrapper])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TestTinyFetiCoSimulationCases])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TestThermalRomCoSim])) + smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([Test3D1DDataTransferProcessBlock])) + smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([Test3D1DDataTransferProcessTorus])) ################################################################################ diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/small_strains/fatigue/generic_small_strain_high_cycle_fatigue_law.h b/applications/ConstitutiveLawsApplication/custom_constitutive/small_strains/fatigue/generic_small_strain_high_cycle_fatigue_law.h index 23fb370a8ddc..b2526fbe580d 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/small_strains/fatigue/generic_small_strain_high_cycle_fatigue_law.h +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/small_strains/fatigue/generic_small_strain_high_cycle_fatigue_law.h @@ -451,7 +451,7 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) GenericSmallStrainHighCycleFatig void save(Serializer &rSerializer) const override { - KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, ConstitutiveLaw) + KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, BaseType) rSerializer.save("FatigueReductionFactor", mFatigueReductionFactor); rSerializer.save("PreviousStresses", mPreviousStresses); rSerializer.save("MaxStress", mMaxStress); @@ -476,7 +476,7 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) GenericSmallStrainHighCycleFatig void load(Serializer &rSerializer) override { - KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, ConstitutiveLaw) + KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, BaseType) rSerializer.load("FatigueReductionFactor", mFatigueReductionFactor); rSerializer.load("PreviousStresses", mPreviousStresses); rSerializer.load("MaxStress", mMaxStress); diff --git a/applications/ConstitutiveLawsApplication/custom_processes/advance_in_time_high_cycle_fatigue_process.h b/applications/ConstitutiveLawsApplication/custom_processes/advance_in_time_high_cycle_fatigue_process.h index 7699ea420dc6..e0a51f63e624 100644 --- a/applications/ConstitutiveLawsApplication/custom_processes/advance_in_time_high_cycle_fatigue_process.h +++ b/applications/ConstitutiveLawsApplication/custom_processes/advance_in_time_high_cycle_fatigue_process.h @@ -52,8 +52,6 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION)AdvanceInTimeHighCycleFatigueProc ///@name Enum's ///@{ -protected: - public: static constexpr double tolerance = std::numeric_limits::epsilon(); @@ -98,10 +96,26 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION)AdvanceInTimeHighCycleFatigueProc */ void TimeAndCyclesUpdate(const double Increment); -protected: +private: // Member Variables - ModelPart& mrModelPart; // The model part to compute - Parameters mThisParameters; + ModelPart& mrModelPart; // The model part to compute + Parameters mThisParameters; // The project parameters + + friend class Serializer; + + void save(Serializer &rSerializer) const override + { + KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, Process) + rSerializer.save("ModelPart", mrModelPart); + rSerializer.save("ThisParameters", mThisParameters); + } + + void load(Serializer &rSerializer) override + { + KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, Process) + rSerializer.load("ModelPart", mrModelPart); + rSerializer.load("ThisParameters", mThisParameters); + } }; // Class AdvanceInTimeHighCycleFatigueProcess diff --git a/applications/DEMApplication/DEM_application.cpp b/applications/DEMApplication/DEM_application.cpp index 6cb1ef5dafba..6a05ec219e23 100644 --- a/applications/DEMApplication/DEM_application.cpp +++ b/applications/DEMApplication/DEM_application.cpp @@ -101,6 +101,7 @@ KRATOS_CREATE_VARIABLE(int, ROTATION_OPTION) KRATOS_CREATE_VARIABLE(int, VIRTUAL_MASS_OPTION) KRATOS_CREATE_VARIABLE(int, SEARCH_CONTROL) KRATOS_CREATE_VARIABLE(bool, IS_TIME_TO_PRINT) +KRATOS_CREATE_VARIABLE(bool, IS_TIME_TO_UPDATE_CONTACT_ELEMENT) KRATOS_CREATE_VARIABLE(double, COORDINATION_NUMBER) KRATOS_CREATE_VARIABLE(double, CONTINUUM_SEARCH_RADIUS_AMPLIFICATION_FACTOR) KRATOS_CREATE_VARIABLE(double, MAX_AMPLIFICATION_RATIO_OF_THE_SEARCH_RADIUS) @@ -153,6 +154,7 @@ KRATOS_CREATE_VARIABLE(bool, DELTA_OPTION) KRATOS_CREATE_VARIABLE(double, SKIN_SPHERE) KRATOS_CREATE_VARIABLE(int, PROPERTIES_ID) KRATOS_CREATE_VARIABLE(int, CONTACT_MESH_OPTION) +KRATOS_CREATE_VARIABLE(int, BOUNDING_BOX_SERVO_LOADING_OPTION) KRATOS_CREATE_VARIABLE(double, MAX_NUMBER_OF_INTACT_BONDS_TO_CONSIDER_A_SPHERE_BROKEN) //KRATOS_CREATE_VARIABLE(int, FAILURE_CRITERION_OPTION) KRATOS_CREATE_VARIABLE(int, CONCRETE_TEST_OPTION) @@ -612,6 +614,7 @@ void KratosDEMApplication::Register() { KRATOS_REGISTER_VARIABLE(VIRTUAL_MASS_OPTION) KRATOS_REGISTER_VARIABLE(SEARCH_CONTROL) KRATOS_REGISTER_VARIABLE(IS_TIME_TO_PRINT) + KRATOS_REGISTER_VARIABLE(IS_TIME_TO_UPDATE_CONTACT_ELEMENT) KRATOS_REGISTER_VARIABLE(COORDINATION_NUMBER) KRATOS_REGISTER_VARIABLE(CONTINUUM_SEARCH_RADIUS_AMPLIFICATION_FACTOR) KRATOS_REGISTER_VARIABLE(MAX_AMPLIFICATION_RATIO_OF_THE_SEARCH_RADIUS) @@ -660,6 +663,7 @@ void KratosDEMApplication::Register() { KRATOS_REGISTER_VARIABLE(SKIN_SPHERE) KRATOS_REGISTER_VARIABLE(PROPERTIES_ID) KRATOS_REGISTER_VARIABLE(CONTACT_MESH_OPTION) + KRATOS_REGISTER_VARIABLE(BOUNDING_BOX_SERVO_LOADING_OPTION) KRATOS_REGISTER_VARIABLE(MAX_NUMBER_OF_INTACT_BONDS_TO_CONSIDER_A_SPHERE_BROKEN) //KRATOS_REGISTER_VARIABLE(FAILURE_CRITERION_OPTION) KRATOS_REGISTER_VARIABLE(CONCRETE_TEST_OPTION) diff --git a/applications/DEMApplication/DEM_application_variables.h b/applications/DEMApplication/DEM_application_variables.h index 2cf219a098cb..a45f1d4753a1 100644 --- a/applications/DEMApplication/DEM_application_variables.h +++ b/applications/DEMApplication/DEM_application_variables.h @@ -68,6 +68,7 @@ namespace Kratos KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, VIRTUAL_MASS_OPTION) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, SEARCH_CONTROL) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, bool, IS_TIME_TO_PRINT) + KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, bool, IS_TIME_TO_UPDATE_CONTACT_ELEMENT) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, COORDINATION_NUMBER) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, CONTINUUM_SEARCH_RADIUS_AMPLIFICATION_FACTOR) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, MAX_AMPLIFICATION_RATIO_OF_THE_SEARCH_RADIUS) @@ -125,6 +126,7 @@ namespace Kratos KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, PROPERTIES_ID) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, CONTACT_MESH_OPTION) + KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, BOUNDING_BOX_SERVO_LOADING_OPTION) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, MAX_NUMBER_OF_INTACT_BONDS_TO_CONSIDER_A_SPHERE_BROKEN) //KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, FAILURE_CRITERION_OPTION) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, CONCRETE_TEST_OPTION) diff --git a/applications/DEMApplication/custom_elements/spheric_particle.cpp b/applications/DEMApplication/custom_elements/spheric_particle.cpp index cbb3004337dd..ce7366a02521 100644 --- a/applications/DEMApplication/custom_elements/spheric_particle.cpp +++ b/applications/DEMApplication/custom_elements/spheric_particle.cpp @@ -894,6 +894,11 @@ void SphericParticle::ComputeBallToBallContactForceAndMoment(SphericParticle::Pa if ((i < (int)mNeighbourElements.size()) && this->Id() < neighbour_iterator_id) { CalculateOnContactElements(i, LocalContactForce, GlobalContactForce); } + } else if (r_process_info[IS_TIME_TO_UPDATE_CONTACT_ELEMENT] && r_process_info[CONTACT_MESH_OPTION] == 1) { + unsigned int neighbour_iterator_id = data_buffer.mpOtherParticle->Id(); + if ((i < (int)mNeighbourElements.size()) && this->Id() < neighbour_iterator_id) { + CalculateOnContactElements(i, LocalContactForce, GlobalContactForce); + } } if (r_process_info[ENERGY_CALCULATION_OPTION]){ diff --git a/applications/DEMApplication/custom_python/DEM_python_application.cpp b/applications/DEMApplication/custom_python/DEM_python_application.cpp index 44bbc81b2c0f..0c3fa1f75fcc 100644 --- a/applications/DEMApplication/custom_python/DEM_python_application.cpp +++ b/applications/DEMApplication/custom_python/DEM_python_application.cpp @@ -71,6 +71,7 @@ PYBIND11_MODULE(KratosDEMApplication,m) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, AUTOMATIC_SKIN_COMPUTATION) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, SKIN_FACTOR_RADIUS) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, IS_TIME_TO_PRINT) + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, IS_TIME_TO_UPDATE_CONTACT_ELEMENT) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, COORDINATION_NUMBER) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, CONTINUUM_SEARCH_RADIUS_AMPLIFICATION_FACTOR) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, MAX_AMPLIFICATION_RATIO_OF_THE_SEARCH_RADIUS) @@ -125,6 +126,7 @@ PYBIND11_MODULE(KratosDEMApplication,m) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, COHESIVE_GROUP) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PROPERTIES_ID) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, CONTACT_MESH_OPTION) + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, BOUNDING_BOX_SERVO_LOADING_OPTION) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, MAX_NUMBER_OF_INTACT_BONDS_TO_CONSIDER_A_SPHERE_BROKEN) //KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, FAILURE_CRITERION_OPTION) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, CONCRETE_TEST_OPTION) diff --git a/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.cpp b/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.cpp index 74c115e07a96..154bd71697be 100644 --- a/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.cpp +++ b/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.cpp @@ -459,6 +459,9 @@ namespace Kratos { if (is_time_to_print_results && r_process_info[CONTACT_MESH_OPTION] == 1) { CreateContactElements(); InitializeContactElements(); + }else if (r_process_info[IS_TIME_TO_UPDATE_CONTACT_ELEMENT] && r_process_info[CONTACT_MESH_OPTION] == 1) { + CreateContactElements(); + InitializeContactElements(); } //RebuildPropertiesProxyPointers(mListOfSphericParticles); diff --git a/applications/DEMApplication/python_scripts/DEM_analysis_stage.py b/applications/DEMApplication/python_scripts/DEM_analysis_stage.py index 94f647c3c0e6..8614ea26f5e7 100644 --- a/applications/DEMApplication/python_scripts/DEM_analysis_stage.py +++ b/applications/DEMApplication/python_scripts/DEM_analysis_stage.py @@ -245,6 +245,7 @@ def FillAnalyticSubModelPartsWithNewParticles(self): def Initialize(self): self.time = 0.0 self.time_old_print = 0.0 + self.time_old_update_contact_element = 0.0 self.ReadModelParts() @@ -309,6 +310,13 @@ def Initialize(self): self.BoundingBoxMaxY_update = self.DEM_parameters["BoundingBoxMaxY"].GetDouble() self.BoundingBoxMaxZ_update = self.DEM_parameters["BoundingBoxMaxZ"].GetDouble() + self.bounding_box_servo_loading_option = False + if "BoundingBoxServoLoadingOption" in self.DEM_parameters.keys(): + if self.DEM_parameters["BoundingBoxServoLoadingOption"].GetBool(): + self.bounding_box_servo_loading_option = True + + self.bounding_box_move_velocity = [0.0, 0.0, 0.0] + def SetMaterials(self): self.ReadMaterialsFile() @@ -498,6 +506,10 @@ def RunAnalytics(self, time): def IsTimeToPrintPostProcess(self): return self.do_print_results_option and self.DEM_parameters["OutputTimeStep"].GetDouble() - (self.time - self.time_old_print) < 1e-2 * self._GetSolver().dt + def IsTimeToUpdateContactElementForServo(self): + return self.bounding_box_servo_loading_option and self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingFrequency"].GetInt() * \ + self._GetSolver().dt - (self.time - self.time_old_update_contact_element) < 1e-2 * self._GetSolver().dt + def PrintResults(self): #### GiD IO ########################################## if self.IsTimeToPrintPostProcess(): @@ -523,6 +535,8 @@ def InitializeSolutionStep(self): self.FillAnalyticSubModelPartsWithNewParticles() if self.DEM_parameters["ContactMeshOption"].GetBool(): self.UpdateIsTimeToPrintInModelParts(self.IsTimeToPrintPostProcess()) + if self.bounding_box_servo_loading_option: + self.UpdateIsTimeToUpdateContactElementForServo(self.IsTimeToUpdateContactElementForServo()) if self.DEM_parameters["Dimension"].GetInt() == 2: self.spheres_model_part.ProcessInfo[IMPOSED_Z_STRAIN_OPTION] = self.DEM_parameters["ImposeZStrainIn2DOption"].GetBool() @@ -530,32 +544,98 @@ def InitializeSolutionStep(self): if self.spheres_model_part.ProcessInfo[IMPOSED_Z_STRAIN_OPTION]: self.spheres_model_part.ProcessInfo.SetValue(IMPOSED_Z_STRAIN_VALUE, eval(self.DEM_parameters["ZStrainValue"].GetString())) + bounding_box_servo_loading_option = False + if "BoundingBoxServoLoadingOption" in self.DEM_parameters.keys(): + if self.DEM_parameters["BoundingBoxServoLoadingOption"].GetBool(): + bounding_box_servo_loading_option = True + if "BoundingBoxMoveOption" in self.DEM_parameters.keys(): if self.DEM_parameters["BoundingBoxMoveOption"].GetBool(): time_step = self.spheres_model_part.ProcessInfo[TIME_STEPS] - NStepSearch = self.DEM_parameters["NeighbourSearchFrequency"].GetInt() - if (time_step + 1) % NStepSearch == 0 and (time_step > 0): - self.UpdateSearchStartegyAndCPlusPlusStrategy() - self.procedures.UpdateBoundingBox(self.spheres_model_part, self.creator_destructor) + if bounding_box_servo_loading_option: + NStepSearch = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingFrequency"].GetInt() + if (time_step + 1) % NStepSearch == 0 and (time_step > 0): + measured_global_stress = self.MeasureSphereForGettingGlobalStressTensor() + self.CalculateBoundingBoxMoveVelocity(measured_global_stress) + self.UpdateSearchStartegyAndCPlusPlusStrategy(self.bounding_box_move_velocity) + self.procedures.UpdateBoundingBox(self.spheres_model_part, self.creator_destructor, self.bounding_box_move_velocity) + else: + NStepSearch = self.DEM_parameters["NeighbourSearchFrequency"].GetInt() + if (time_step + 1) % NStepSearch == 0 and (time_step > 0): + bounding_box_move_velocity = self.DEM_parameters["BoundingBoxMoveVelocity"].GetVector() + self.UpdateSearchStartegyAndCPlusPlusStrategy(bounding_box_move_velocity) + self.procedures.UpdateBoundingBox(self.spheres_model_part, self.creator_destructor, bounding_box_move_velocity) + + def CalculateBoundingBoxMoveVelocity(self, measured_global_stress): + + # note: servo_loading_type = "isotropic" or "anisotropic" + servo_loading_type = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingType"].GetString() + + if servo_loading_type == "isotropic": + self.CalculateBoundingBoxMoveVelocityIsotropic(measured_global_stress) + elif servo_loading_type == "anisotropic": + self.CalculateBoundingBoxMoveVelocityAnisotropic(measured_global_stress) + + def CalculateBoundingBoxMoveVelocityIsotropic(self, measured_global_stress): + + servo_loading_stress = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingStress"].GetVector() + servo_loading_factor = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingFactor"].GetDouble() + d50 = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["MeanParticleDiameterD50"].GetDouble() + Youngs_modulus = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["ParticleYoungsModulus"].GetDouble() + + servo_loading_stress_mean = (servo_loading_stress[0] + servo_loading_stress[1] + servo_loading_stress[2]) / 3.0 + measured_global_stress_mean = (measured_global_stress[0][0] + measured_global_stress[1][1] + measured_global_stress[2][2]) / 3.0 + + delta_time = self.spheres_model_part.ProcessInfo.GetValue(DELTA_TIME) + servo_loading_coefficient = servo_loading_factor * d50 / (delta_time * Youngs_modulus) + self.bounding_box_move_velocity[0] = (servo_loading_stress_mean - measured_global_stress_mean) * servo_loading_coefficient + + bonuding_box_serva_loading_velocity_max = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingVelocityMax"].GetDouble() + if abs(self.bounding_box_move_velocity[0]) > bonuding_box_serva_loading_velocity_max: + self.bounding_box_move_velocity[0] = bonuding_box_serva_loading_velocity_max * np.sign(self.bounding_box_move_velocity[0]) + + self.bounding_box_move_velocity[1] = self.bounding_box_move_velocity[0] + self.bounding_box_move_velocity[2] = self.bounding_box_move_velocity[0] - def UpdateSearchStartegyAndCPlusPlusStrategy(self): + def CalculateBoundingBoxMoveVelocityAnisotropic(self, measured_global_stress): + + servo_loading_stress = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingStress"].GetVector() + servo_loading_factor = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingFactor"].GetDouble() + d50 = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["MeanParticleDiameterD50"].GetDouble() + Youngs_modulus = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["ParticleYoungsModulus"].GetDouble() + + delta_time = self.spheres_model_part.ProcessInfo.GetValue(DELTA_TIME) + servo_loading_coefficient = servo_loading_factor * d50 / (delta_time * Youngs_modulus) + self.bounding_box_move_velocity[0] = (servo_loading_stress[0] - measured_global_stress[0][0]) * servo_loading_coefficient + self.bounding_box_move_velocity[1] = (servo_loading_stress[1] - measured_global_stress[1][1]) * servo_loading_coefficient + self.bounding_box_move_velocity[2] = (servo_loading_stress[2] - measured_global_stress[2][2]) * servo_loading_coefficient + + bonuding_box_move_velocity_max = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingVelocityMax"].GetDouble() + if abs(self.bounding_box_move_velocity[0]) > bonuding_box_move_velocity_max: + self.bounding_box_move_velocity[0] = bonuding_box_move_velocity_max * np.sign(self.bounding_box_move_velocity[0]) + if abs(self.bounding_box_move_velocity[1]) > bonuding_box_move_velocity_max: + self.bounding_box_move_velocity[1] = bonuding_box_move_velocity_max * np.sign(self.bounding_box_move_velocity[1]) + if abs(self.bounding_box_move_velocity[2]) > bonuding_box_move_velocity_max: + self.bounding_box_move_velocity[2] = bonuding_box_move_velocity_max * np.sign(self.bounding_box_move_velocity[2]) + + def UpdateSearchStartegyAndCPlusPlusStrategy(self, move_velocity): delta_time = self.spheres_model_part.ProcessInfo.GetValue(DELTA_TIME) - move_velocity = self.DEM_parameters["BoundingBoxMoveVelocity"].GetDouble() + # note: control_bool_vector = [MoveMinX, MoveMinY, MoveMinZ, MoveMaxX, MoveMaxY, MoveMaxZ] control_bool_vector = self.DEM_parameters["BoundingBoxMoveOptionDetail"].GetVector() if control_bool_vector[0]: - self.BoundingBoxMinX_update += delta_time * move_velocity + self.BoundingBoxMinX_update += delta_time * move_velocity[0] if control_bool_vector[1]: - self.BoundingBoxMinY_update += delta_time * move_velocity + self.BoundingBoxMinY_update += delta_time * move_velocity[1] if control_bool_vector[2]: - self.BoundingBoxMinZ_update += delta_time * move_velocity + self.BoundingBoxMinZ_update += delta_time * move_velocity[2] if control_bool_vector[3]: - self.BoundingBoxMaxX_update -= delta_time * move_velocity + self.BoundingBoxMaxX_update -= delta_time * move_velocity[0] if control_bool_vector[4]: - self.BoundingBoxMaxY_update -= delta_time * move_velocity + self.BoundingBoxMaxY_update -= delta_time * move_velocity[1] if control_bool_vector[5]: - self.BoundingBoxMaxZ_update -= delta_time * move_velocity + self.BoundingBoxMaxZ_update -= delta_time * move_velocity[2] self._GetSolver().search_strategy = OMP_DEMSearch(self.BoundingBoxMinX_update, self.BoundingBoxMinY_update, @@ -571,6 +651,13 @@ def UpdateIsTimeToPrintInModelParts(self, is_time_to_print): self.UpdateIsTimeToPrintInOneModelPart(self.dem_inlet_model_part, is_time_to_print) self.UpdateIsTimeToPrintInOneModelPart(self.rigid_face_model_part, is_time_to_print) + def UpdateIsTimeToUpdateContactElementForServo(self, is_time_to_update_contact_element): + self.spheres_model_part.ProcessInfo[IS_TIME_TO_UPDATE_CONTACT_ELEMENT] = is_time_to_update_contact_element + self.cluster_model_part.ProcessInfo[IS_TIME_TO_UPDATE_CONTACT_ELEMENT] = is_time_to_update_contact_element + self.dem_inlet_model_part.ProcessInfo[IS_TIME_TO_UPDATE_CONTACT_ELEMENT] = is_time_to_update_contact_element + self.rigid_face_model_part.ProcessInfo[IS_TIME_TO_UPDATE_CONTACT_ELEMENT] = is_time_to_update_contact_element + self.time_old_update_contact_element = self.time + def UpdateIsTimeToPrintInOneModelPart(self, model_part, is_time_to_print): model_part.ProcessInfo[IS_TIME_TO_PRINT] = is_time_to_print @@ -742,7 +829,7 @@ def OutputSingleTimeLoop(self): self.PrintResultsForGid(self.time) self.time_old_print = self.time - def MeasureSphereForGettingPackingProperties(self, radius, center_x, center_y, center_z, type): + def MeasureSphereForGettingPackingProperties(self, radius, center_x, center_y, center_z, type, domain_size=[1,1,1]): ''' This is a function to establish a sphere to measure local packing properties The type could be "porosity", "averaged_coordination_number", "fabric_tensor", "stress_tensor" or "strain" @@ -781,6 +868,103 @@ def MeasureSphereForGettingPackingProperties(self, radius, center_x, center_y, c return measured_porosity + if type == "porosity_periodic": + + measure_sphere_volume = 4.0 / 3.0 * math.pi * radius * radius * radius + sphere_volume_inside_range = 0.0 + measured_porosity = 0.0 + + center_list = [] + center_list.append([center_x, center_y, center_z]) + if center_x + 0.5 * domain_size[0] < 2 * radius: + center_list.append([center_x + domain_size[0], center_y, center_z]) + if center_y + 0.5 * domain_size[1] < 2 * radius: + center_list.append([center_x + domain_size[0], center_y + domain_size[1], center_z]) + if center_y - 0.5 * domain_size[1] > -2 * radius: + center_list.append([center_x + domain_size[0], center_y - domain_size[1], center_z]) + if center_y + 0.5 * domain_size[2] < 2 * radius: + center_list.append([center_x + domain_size[0], center_y, center_z + domain_size[2]]) + if center_y - 0.5 * domain_size[2] > -2 * radius: + center_list.append([center_x + domain_size[0], center_y, center_z - domain_size[2]]) + if center_x - 0.5 * domain_size[0] > -2 * radius: + center_list.append([center_x - domain_size[0], center_y, center_z]) + if center_y + 0.5 * domain_size[1] < 2 * radius: + center_list.append([center_x - domain_size[0], center_y + domain_size[1], center_z]) + if center_y - 0.5 * domain_size[1] > -2 * radius: + center_list.append([center_x - domain_size[0], center_y - domain_size[1], center_z]) + if center_y + 0.5 * domain_size[2] < 2 * radius: + center_list.append([center_x - domain_size[0], center_y, center_z + domain_size[2]]) + if center_y - 0.5 * domain_size[2] > -2 * radius: + center_list.append([center_x - domain_size[0], center_y, center_z - domain_size[2]]) + if center_y + 0.5 * domain_size[1] < 2 * radius: + center_list.append([center_x, center_y + domain_size[1], center_z]) + if center_x + 0.5 * domain_size[0] < 2 * radius: + center_list.append([center_x + domain_size[0], center_y + domain_size[1], center_z]) + if center_x - 0.5 * domain_size[0] > -2 * radius: + center_list.append([center_x - domain_size[0], center_y + domain_size[1], center_z]) + if center_z + 0.5 * domain_size[2] < 2 * radius: + center_list.append([center_x, center_y + domain_size[1], center_z + domain_size[2]]) + if center_z - 0.5 * domain_size[2] > -2 * radius: + center_list.append([center_x, center_y + domain_size[1], center_z - domain_size[2]]) + if center_y - 0.5 * domain_size[1] > -2 * radius: + center_list.append([center_x, center_y - domain_size[1], center_z]) + if center_x + 0.5 * domain_size[0] < 2 * radius: + center_list.append([center_x + domain_size[0], center_y - domain_size[1], center_z]) + if center_x - 0.5 * domain_size[0] > -2 * radius: + center_list.append([center_x - domain_size[0], center_y - domain_size[1], center_z]) + if center_z + 0.5 * domain_size[2] < 2 * radius: + center_list.append([center_x, center_y - domain_size[1], center_z + domain_size[2]]) + if center_z - 0.5 * domain_size[2] > -2 * radius: + center_list.append([center_x, center_y - domain_size[1], center_z - domain_size[2]]) + if center_z + 0.5 * domain_size[2] < 2 * radius: + center_list.append([center_x, center_y, center_z + domain_size[2]]) + if center_x + 0.5 * domain_size[0] < 2 * radius: + center_list.append([center_x + domain_size[0], center_y, center_z + domain_size[2]]) + if center_x - 0.5 * domain_size[0] > -2 * radius: + center_list.append([center_x - domain_size[0], center_y, center_z + domain_size[2]]) + if center_y + 0.5 * domain_size[1] < 2 * radius: + center_list.append([center_x, center_y + domain_size[1], center_z + domain_size[2]]) + if center_y - 0.5 * domain_size[1] > -2 * radius: + center_list.append([center_x, center_y - domain_size[1], center_z + domain_size[2]]) + if center_z - 0.5 * domain_size[2] > -2 * radius: + center_list.append([center_x, center_y, center_z - domain_size[2]]) + if center_x + 0.5 * domain_size[0] < 2 * radius: + center_list.append([center_x + domain_size[0], center_y, center_z - domain_size[2]]) + if center_x - 0.5 * domain_size[0] > -2 * radius: + center_list.append([center_x - domain_size[0], center_y, center_z - domain_size[2]]) + if center_y + 0.5 * domain_size[1] < 2 * radius: + center_list.append([center_x, center_y + domain_size[1], center_z - domain_size[2]]) + if center_y - 0.5 * domain_size[1] > -2 * radius: + center_list.append([center_x, center_y - domain_size[1], center_z - domain_size[2]]) + + for center in center_list: + for node in self.spheres_model_part.Nodes: + + r = node.GetSolutionStepValue(RADIUS) + x = node.X + y = node.Y + z = node.Z + + center_to_sphere_distance = ((x - center[0])**2 + (y - center[1])**2 + (z - center[2])**2)**0.5 + + if center_to_sphere_distance < (radius - r): + + sphere_volume_inside_range += 4/3 * math.pi * r * r * r + + elif center_to_sphere_distance <= (radius + r): + + other_part_d = radius - (radius * radius + center_to_sphere_distance * center_to_sphere_distance - r * r) / (center_to_sphere_distance * 2) + + my_part_d = r - (r * r + center_to_sphere_distance * center_to_sphere_distance - radius * radius) / (center_to_sphere_distance * 2) + + cross_volume = math.pi * other_part_d * other_part_d * (radius - 1/3 * other_part_d) + math.pi * my_part_d * my_part_d * (r - 1/3 * my_part_d) + + sphere_volume_inside_range += cross_volume + + measured_porosity = 1.0 - (sphere_volume_inside_range / measure_sphere_volume) + + return measured_porosity + if type == "porosity_kdtree": measure_sphere_volume = 4.0 / 3.0 * math.pi * radius * radius * radius @@ -1107,6 +1291,47 @@ def MeasureSphereForGettingPackingProperties(self, radius, center_x, center_y, c if type == "strain": pass + def MeasureSphereForGettingGlobalStressTensor(self): + + if self.DEM_parameters["PostStressStrainOption"].GetBool() and self.DEM_parameters["ContactMeshOption"].GetBool(): + + bounding_box_volume = (self.BoundingBoxMaxX_update - self.BoundingBoxMinX_update) * \ + (self.BoundingBoxMaxY_update - self.BoundingBoxMinY_update) * \ + (self.BoundingBoxMaxZ_update - self.BoundingBoxMinZ_update) + + total_tensor = np.empty((3, 3)) + total_tensor[:] = 0.0 + + for element in self.contact_model_part.Elements: + + x_0 = element.GetNode(0).X + x_1 = element.GetNode(1).X + y_0 = element.GetNode(0).Y + y_1 = element.GetNode(1).Y + z_0 = element.GetNode(0).Z + z_1 = element.GetNode(1).Z + + if abs(x_0 - x_1) > 0.5 * (self.BoundingBoxMaxX_update - self.BoundingBoxMinX_update) or \ + abs(y_0 - y_1) > 0.5 * (self.BoundingBoxMaxY_update - self.BoundingBoxMinY_update) or \ + abs(z_0 - z_1) > 0.5 * (self.BoundingBoxMaxZ_update - self.BoundingBoxMinZ_update): + continue + + local_contact_force_X = element.GetValue(GLOBAL_CONTACT_FORCE)[0] + local_contact_force_Y = element.GetValue(GLOBAL_CONTACT_FORCE)[1] + local_contact_force_Z = element.GetValue(GLOBAL_CONTACT_FORCE)[2] + contact_force_vector = np.array([local_contact_force_X , local_contact_force_Y, local_contact_force_Z]) + vector_l = np.array([x_0 - x_1 , y_0 - y_1, z_0 - z_1]) + tensor = np.outer(contact_force_vector, vector_l) + total_tensor += tensor + + total_tensor = total_tensor / bounding_box_volume + + return total_tensor + + else: + + raise Exception('The \"PostStressStrainOption\" and \"ContactMeshOption\" in the [ProjectParametersDEM.json] should be [True].') + def MeasureSphereForGettingRadialDistributionFunction(self, radius, center_x, center_y, center_z, delta_r, d_mean): min_reference_particle_to_center_distance = 1e10 @@ -1643,9 +1868,6 @@ def MeasureCubicForGettingPackingProperties(self, side_length, center_x, center_ if total_contact_number: averaged_contact_force_modulus_square = total_contact_force_tensor_modulus_square / total_contact_number - print(averaged_total_particle_force_tensor_modulus_square) - print(averaged_contact_force_modulus_square) - if averaged_contact_force_modulus_square: return (averaged_total_particle_force_tensor_modulus_square / averaged_contact_force_modulus_square)**0.5 else: diff --git a/applications/DEMApplication/python_scripts/DEM_procedures.py b/applications/DEMApplication/python_scripts/DEM_procedures.py index ba5e3df9dec8..5201c4266098 100644 --- a/applications/DEMApplication/python_scripts/DEM_procedures.py +++ b/applications/DEMApplication/python_scripts/DEM_procedures.py @@ -826,28 +826,25 @@ def SetBoundingBox(self, spheres_model_part, clusters_model_part, rigid_faces_mo creator_destructor.SetHighNode(b_box_high) creator_destructor.CalculateSurroundingBoundingBox(spheres_model_part, clusters_model_part, rigid_faces_model_part, dem_inlet_model_part, self.bounding_box_enlargement_factor, self.automatic_bounding_box_OPTION) - def UpdateBoundingBox(self, spheres_model_part, creator_destructor): + def UpdateBoundingBox(self, spheres_model_part, creator_destructor, move_velocity): delta_time = spheres_model_part.ProcessInfo.GetValue(DELTA_TIME) - #NStepSearch = self.DEM_parameters["NeighbourSearchFrequency"].GetInt() - #delta_time_real = delta_time * NStepSearch - move_velocity = self.DEM_parameters["BoundingBoxMoveVelocity"].GetDouble() b_box_low = Array3() b_box_high = Array3() control_bool_vector = self.DEM_parameters["BoundingBoxMoveOptionDetail"].GetVector() if control_bool_vector[0]: - self.b_box_minX += delta_time * move_velocity + self.b_box_minX += delta_time * move_velocity[0] if control_bool_vector[1]: - self.b_box_minY += delta_time * move_velocity + self.b_box_minY += delta_time * move_velocity[1] if control_bool_vector[2]: - self.b_box_minZ += delta_time * move_velocity + self.b_box_minZ += delta_time * move_velocity[2] if control_bool_vector[3]: - self.b_box_maxX -= delta_time * move_velocity + self.b_box_maxX -= delta_time * move_velocity[0] if control_bool_vector[4]: - self.b_box_maxY -= delta_time * move_velocity + self.b_box_maxY -= delta_time * move_velocity[1] if control_bool_vector[5]: - self.b_box_maxZ -= delta_time * move_velocity + self.b_box_maxZ -= delta_time * move_velocity[2] b_box_low[0] = self.b_box_minX b_box_low[1] = self.b_box_minY b_box_low[2] = self.b_box_minZ diff --git a/applications/DEMApplication/python_scripts/dem_default_input_parameters.py b/applications/DEMApplication/python_scripts/dem_default_input_parameters.py index 52476b1f7d73..e35bae680caa 100644 --- a/applications/DEMApplication/python_scripts/dem_default_input_parameters.py +++ b/applications/DEMApplication/python_scripts/dem_default_input_parameters.py @@ -20,7 +20,17 @@ def GetDefaultInputParameters(): "BoundingBoxMinZ" : -10.0, "BoundingBoxMoveOption" : false, "BoundingBoxMoveOptionDetail" : [1, 1, 1, 1, 1, 1], - "BoundingBoxMoveVelocity" : 0.001, + "BoundingBoxMoveVelocity" : [0.01, 0.01, 0.01], + "BoundingBoxServoLoadingOption" : false, + "BoundingBoxServoLoadingSettings":{ + "BoundingBoxServoLoadingType" : "isotropic", + "BoundingBoxServoLoadingStress" : [0.0, 0.0, 0.0], + "BoundingBoxServoLoadingFactor" : 0.8, + "BoundingBoxServoLoadingFrequency" : 50, + "BoundingBoxServoLoadingVelocityMax": 0.05, + "MeanParticleDiameterD50" : 0.1, + "ParticleYoungsModulus" : 1e9 + }, "dem_inlet_option" : true, "dem_inlets_settings" : {}, "seed" : 42, diff --git a/applications/DEMApplication/python_scripts/sphere_strategy.py b/applications/DEMApplication/python_scripts/sphere_strategy.py index d9ccadeebd4e..3d8599bcbfe5 100644 --- a/applications/DEMApplication/python_scripts/sphere_strategy.py +++ b/applications/DEMApplication/python_scripts/sphere_strategy.py @@ -81,6 +81,10 @@ def __init__(self, all_model_parts, creator_destructor, dem_fem_search, DEM_para self.contact_mesh_option = DEM_parameters["ContactMeshOption"].GetBool() self.automatic_bounding_box_option = DEM_parameters["AutomaticBoundingBoxOption"].GetBool() + self.bounding_box_servo_loading_option = 0 + if "BoundingBoxServoLoadingOption" in DEM_parameters.keys(): + self.bounding_box_servo_loading_option = DEM_parameters["BoundingBoxServoLoadingOption"].GetBool() + self.delta_option = DEM_parameters["DeltaOption"].GetString() #TODO: this is not an option (bool) let's change the name to something including 'type' self.search_increment = 0.0 @@ -322,6 +326,11 @@ def SetVariablesAndOptions(self): else: self.spheres_model_part.ProcessInfo.SetValue(CONTACT_MESH_OPTION, 0) + if self.bounding_box_servo_loading_option: + self.spheres_model_part.ProcessInfo.SetValue(BOUNDING_BOX_SERVO_LOADING_OPTION, 1) + else: + self.spheres_model_part.ProcessInfo.SetValue(BOUNDING_BOX_SERVO_LOADING_OPTION, 0) + # PRINTING VARIABLES self.spheres_model_part.ProcessInfo.SetValue(PRINT_EXPORT_ID, self.print_export_id) diff --git a/applications/DEMApplication/tests/moving_periodic_boundary_tests_files/ProjectParametersDEM.json b/applications/DEMApplication/tests/moving_periodic_boundary_tests_files/ProjectParametersDEM.json index f3762b36437a..7f94c9e90c2e 100644 --- a/applications/DEMApplication/tests/moving_periodic_boundary_tests_files/ProjectParametersDEM.json +++ b/applications/DEMApplication/tests/moving_periodic_boundary_tests_files/ProjectParametersDEM.json @@ -13,7 +13,7 @@ "BoundingBoxMinY" : 0.0, "BoundingBoxMinZ" : -0.0015, "BoundingBoxMoveOption" : true, - "BoundingBoxMoveVelocity" : 0.001, + "BoundingBoxMoveVelocity" : [0.001, 0.001, 0.001], "dem_inlet_option" : false, "GravityX" : 0.0, "GravityY" : 0.0, diff --git a/applications/DEMApplication/tests/servo_control_tests_files/MaterialsDEM.json b/applications/DEMApplication/tests/servo_control_tests_files/MaterialsDEM.json new file mode 100644 index 000000000000..2573997bcfd9 --- /dev/null +++ b/applications/DEMApplication/tests/servo_control_tests_files/MaterialsDEM.json @@ -0,0 +1,29 @@ +{ + "materials" : [{ + "material_name" : "DEM-DefaultMaterial", + "material_id" : 1, + "Variables" : { + "PARTICLE_DENSITY" : 2650.0, + "YOUNG_MODULUS" : 3.5e7, + "POISSON_RATIO" : 0.2, + "PARTICLE_SPHERICITY" : 1.0 + } + }], + "material_relations" : [{ + "material_names_list" : ["DEM-DefaultMaterial","DEM-DefaultMaterial"], + "material_ids_list" : [1,1], + "Variables" : { + "DEM_DISCONTINUUM_CONSTITUTIVE_LAW_NAME" : "DEM_D_Hertz_viscous_Coulomb", + "PARTICLE_COHESION" : 0.0, + "STATIC_FRICTION" : 0.6, + "DYNAMIC_FRICTION" : 0.6, + "FRICTION_DECAY" : 500, + "COEFFICIENT_OF_RESTITUTION" : 0.01, + "ROLLING_FRICTION" : 0.01, + "ROLLING_FRICTION_WITH_WALLS" : 0.01, + "DEM_ROLLING_FRICTION_MODEL_NAME" : "DEMRollingFrictionModelConstantTorque" + } + }], + "material_assignation_table" : [["SpheresPart.DEMParts_Body","DEM-DefaultMaterial"]] +} + diff --git a/applications/DEMApplication/tests/servo_control_tests_files/ProjectParametersDEM.json b/applications/DEMApplication/tests/servo_control_tests_files/ProjectParametersDEM.json new file mode 100644 index 000000000000..b6a68a9500be --- /dev/null +++ b/applications/DEMApplication/tests/servo_control_tests_files/ProjectParametersDEM.json @@ -0,0 +1,87 @@ +{ + "Dimension" : 3, + "PeriodicDomainOption" : true, + "BoundingBoxOption" : true, + "AutomaticBoundingBoxOption" : false, + "BoundingBoxEnlargementFactor" : 1.0, + "BoundingBoxStartTime" : 0.0, + "BoundingBoxStopTime" : 1000.0, + "BoundingBoxMaxX" : 0.2, + "BoundingBoxMaxY" : 0.241421356237, + "BoundingBoxMaxZ" : 0.2, + "BoundingBoxMinX" : -0.2, + "BoundingBoxMinY" : -0.241421356237, + "BoundingBoxMinZ" : -0.2, + "BoundingBoxMoveOption" : true, + "BoundingBoxMoveVelocity" : [1.0, 1.0, 1.0], + "BoundingBoxMoveOptionDetail" : [1, 1, 1, 1, 1, 1], + "BoundingBoxServoLoadingOption" : true, + "BoundingBoxServoLoadingSettings":{ + "BoundingBoxServoLoadingType" : "isotropic", + "BoundingBoxServoLoadingStress" : [5e3, 5e3, 5e3], + "BoundingBoxServoLoadingFactor" : 5, + "BoundingBoxServoLoadingFrequency" : 10, + "BoundingBoxServoLoadingVelocityMax": 1.0, + "MeanParticleDiameterD50" : 2.1e-4, + "ParticleYoungsModulus" : 3.5e10 + }, + "dem_inlet_option" : false, + "dem_inlets_settings" : {}, + "GravityX" : 0.0, + "GravityY" : 0.0, + "GravityZ" : 0.0, + "RotationOption" : true, + "CleanIndentationsOption" : false, + "RadiusExpansionOption" : false, + "RadiusExpansionRate" : 1000.0, + "RadiusMultiplierMax" : 2.0, + "RadiusExpansionRateChangeOption": false, + "RadiusExpansionAcceleration" : -4e5, + "RadiusExpansionRateMin" : 100.0, + "EnergyCalculationOption" : true, + "ComputeStressTensorOption" : true, + "solver_settings" : { + "RemoveBallsInitiallyTouchingWalls" : false, + "strategy" : "sphere_strategy", + "material_import_settings" : { + "materials_filename" : "MaterialsDEM.json" + } + }, + "VirtualMassCoefficient" : 1.0, + "RollingFrictionOption" : true, + "GlobalDamping" : 0.0, + "GlobalViscousDamping" : 0.0, + "ContactMeshOption" : true, + "OutputFileType" : "Ascii", + "Multifile" : "multiple_files", + "ElementType" : "SphericPartDEMElement3D", + "TranslationalIntegrationScheme" : "Symplectic_Euler", + "RotationalIntegrationScheme" : "Direct_Integration", + "MaxTimeStep" : 2e-8, + "FinalTime" : 0.00006, + "NeighbourSearchFrequency" : 50, + "SearchTolerance" : 0.0000002, + "GraphExportFreq" : 0.00004, + "VelTrapGraphExportFreq" : 0.00004, + "OutputTimeStep" : 0.00004, + "PostBoundingBox" : true, + "PostLocalContactForce" : false, + "PostDisplacement" : true, + "PostRadius" : true, + "PostVelocity" : true, + "PostAngularVelocity" : false, + "PostElasticForces" : false, + "PostContactForces" : false, + "PostRigidElementForces" : false, + "PostStressStrainOption" : true, + "PostTangentialElasticForces" : false, + "PostTotalForces" : false, + "PostPressure" : false, + "PostShearStress" : false, + "PostSkinSphere" : false, + "PostNonDimensionalVolumeWear" : false, + "PostParticleMoment" : false, + "PostEulerAngles" : false, + "PostRollingResistanceMoment" : false, + "problem_name" : "testSC" +} diff --git a/applications/DEMApplication/tests/servo_control_tests_files/testSCDEM.mdpa b/applications/DEMApplication/tests/servo_control_tests_files/testSCDEM.mdpa new file mode 100644 index 000000000000..84ae1326d776 --- /dev/null +++ b/applications/DEMApplication/tests/servo_control_tests_files/testSCDEM.mdpa @@ -0,0 +1,55 @@ +Begin ModelPartData +// VARIABLE_NAME value +End ModelPartData + +Begin Properties 0 +End Properties + +Begin Nodes +1 -0.1 0.0 0.1 +2 0.1 0.0 0.1 +3 -0.1 0.0 -0.1 +4 0.1 0.0 -0.1 +5 0.0 0.141421356237 0.0 +6 0.0 -0.141421356237 0.0 +End Nodes + + +Begin Elements SphericParticle3D// GUI group identifier: Body +1 0 1 +2 0 2 +3 0 3 +4 0 4 +5 0 5 +6 0 6 +End Elements + +Begin NodalData RADIUS // GUI group identifier: Body + 1 0 0.1 + 2 0 0.1 + 3 0 0.1 + 4 0 0.1 + 5 0 0.1 + 6 0 0.1 +End NodalData + +Begin SubModelPart DEMParts_Body // Group Body // Subtree DEMParts +Begin SubModelPartNodes +1 +2 +3 +4 +5 +6 +End SubModelPartNodes +Begin SubModelPartElements +1 +2 +3 +4 +5 +6 +End SubModelPartElements +Begin SubModelPartConditions +End SubModelPartConditions +End SubModelPart diff --git a/applications/DEMApplication/tests/servo_control_tests_files/testSCDEM_FEM_boundary.mdpa b/applications/DEMApplication/tests/servo_control_tests_files/testSCDEM_FEM_boundary.mdpa new file mode 100644 index 000000000000..c3ac404e2a68 --- /dev/null +++ b/applications/DEMApplication/tests/servo_control_tests_files/testSCDEM_FEM_boundary.mdpa @@ -0,0 +1,6 @@ +Begin ModelPartData + // VARIABLE_NAME value +End ModelPartData + +Begin Properties 0 +End Properties diff --git a/applications/DEMApplication/tests/servo_control_tests_files/testSC_Graphs/EnergyPlot.grf b/applications/DEMApplication/tests/servo_control_tests_files/testSC_Graphs/EnergyPlot.grf new file mode 100644 index 000000000000..355c1f541d85 --- /dev/null +++ b/applications/DEMApplication/tests/servo_control_tests_files/testSC_Graphs/EnergyPlot.grf @@ -0,0 +1,2 @@ + Time Trans kinematic energy Rot kinematic energy Kinematic energy Gravitational energy Elastic energy Frictional energy Viscodamping energy Rolling resistance energy Total energy +4.002e-05 6.17221e-19 1.64829e-28 6.17221e-19 0 9.58515e-12 0 1.26916e-21 2.09691e-25 9.58515e-12 diff --git a/applications/DEMApplication/tests/test_DEMApplication.py b/applications/DEMApplication/tests/test_DEMApplication.py index 987099fd52ef..579f6a9077b6 100644 --- a/applications/DEMApplication/tests/test_DEMApplication.py +++ b/applications/DEMApplication/tests/test_DEMApplication.py @@ -35,6 +35,7 @@ import test_dem_3d_parallel_bond_model import test_dem_3d_smooth_joint_model import test_moving_periodic_boundary +import test_servo_control import sys sys.path.append('DEM3D_chung_ooi_tests/test1_data') sys.path.append('DEM3D_chung_ooi_tests/test2_data') @@ -107,6 +108,7 @@ def AssembleTestSuites(): smallSuite.addTest(test_dem_3d_parallel_bond_model.TestParallelBondModel("test_ParallelBondModel_1")) smallSuite.addTest(test_dem_3d_smooth_joint_model.TestSmoothJointModel("test_SmoothJointModel_1")) smallSuite.addTest(test_moving_periodic_boundary.TestMovingPeriodicBoundary("test_MovingPeriodicBoundary")) + smallSuite.addTest(test_servo_control.TestServoControl("test_ServoControl")) # Create a test suit with the selected tests plus all small tests nightSuite = suites['nightly'] diff --git a/applications/DEMApplication/tests/test_servo_control.py b/applications/DEMApplication/tests/test_servo_control.py new file mode 100644 index 000000000000..92156c684b40 --- /dev/null +++ b/applications/DEMApplication/tests/test_servo_control.py @@ -0,0 +1,52 @@ +import os +import KratosMultiphysics as Kratos +from KratosMultiphysics import Logger +Logger.GetDefaultOutput().SetSeverity(Logger.Severity.WARNING) +import KratosMultiphysics.KratosUnittest as KratosUnittest +import KratosMultiphysics.DEMApplication.DEM_analysis_stage + +import auxiliary_functions_for_tests + +this_working_dir_backup = os.getcwd() + +def GetFilePath(fileName): + return os.path.join(os.path.dirname(os.path.realpath(__file__)), fileName) + +class ServoConrolTestSolution(KratosMultiphysics.DEMApplication.DEM_analysis_stage.DEMAnalysisStage, KratosUnittest.TestCase): + + @classmethod + def GetMainPath(self): + return os.path.join(os.path.dirname(os.path.realpath(__file__)), "servo_control_tests_files") + + def GetProblemNameWithPath(self): + return os.path.join(self.main_path, self.DEM_parameters["problem_name"].GetString()) + + def FinalizeSolutionStep(self): + super().FinalizeSolutionStep() + tolerance = 1e-8 + stress_tensor = self.MeasureSphereForGettingGlobalStressTensor() + mean_stress = (stress_tensor[0][0]+stress_tensor[1][1]+stress_tensor[2][2])/3 + + if self.time >= 0.000048 and self.time < 0.000049: + expected_value = 5.007421025614752e-08 + self.assertAlmostEqual(mean_stress, expected_value, delta=tolerance) + + def Finalize(self): + self.procedures.RemoveFoldersWithResults(str(self.main_path), str(self.problem_name), '') + super().Finalize() + +class TestServoControl(KratosUnittest.TestCase): + + def setUp(self): + pass + + @classmethod + def test_ServoControl(self): + path = os.path.join(os.path.dirname(os.path.realpath(__file__)), "servo_control_tests_files") + parameters_file_name = os.path.join(path, "ProjectParametersDEM.json") + model = Kratos.Model() + auxiliary_functions_for_tests.CreateAndRunStageInSelectedNumberOfOpenMPThreads(ServoConrolTestSolution, model, parameters_file_name, 1) + +if __name__ == "__main__": + Kratos.Logger.GetDefaultOutput().SetSeverity(Logger.Severity.WARNING) + KratosUnittest.main() diff --git a/applications/GeoMechanicsApplication/custom_conditions/general_U_Pw_diff_order_condition.hpp b/applications/GeoMechanicsApplication/custom_conditions/general_U_Pw_diff_order_condition.hpp index 4f5d534532a9..969d4d019b77 100644 --- a/applications/GeoMechanicsApplication/custom_conditions/general_U_Pw_diff_order_condition.hpp +++ b/applications/GeoMechanicsApplication/custom_conditions/general_U_Pw_diff_order_condition.hpp @@ -42,7 +42,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeneralUPwDiffOrderCondition : publi KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(GeneralUPwDiffOrderCondition); - GeneralUPwDiffOrderCondition() : GeneralUPwDiffOrderCondition(0, nullptr, nullptr){}; + GeneralUPwDiffOrderCondition() : GeneralUPwDiffOrderCondition(0, nullptr, nullptr) {}; GeneralUPwDiffOrderCondition(IndexType NewId, GeometryType::Pointer pGeometry) : GeneralUPwDiffOrderCondition(NewId, pGeometry, nullptr) diff --git a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.hpp b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.hpp index ce994beee9d4..af23d882729c 100644 --- a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.hpp @@ -121,8 +121,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SteadyStatePwElement : public Transi void load(Serializer& rSerializer) override{KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, Element)} // Assignment operator. - SteadyStatePwElement& - operator=(SteadyStatePwElement const& rOther); + SteadyStatePwElement& operator=(SteadyStatePwElement const& rOther); // Copy constructor. SteadyStatePwElement(SteadyStatePwElement const& rOther); diff --git a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.cpp b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.cpp index eb6f8f7658ca..47a46b5506d6 100644 --- a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.cpp @@ -13,12 +13,13 @@ // Application includes #include "custom_elements/steady_state_Pw_piping_element.hpp" #include "custom_utilities/element_utilities.hpp" +#include "custom_utilities/transport_equation_utilities.hpp" #include "utilities/math_utils.h" + #include namespace Kratos { - template Element::Pointer SteadyStatePwPipingElement::Create(IndexType NewId, NodesArrayType const& ThisNodes, @@ -273,23 +274,6 @@ double SteadyStatePwPipingElement<3, 8>::CalculateHeadGradient(const PropertiesT << std::endl; } -/// -/// Calculate the particle diameter for the particles in the pipe. The particle diameter equals d70, -/// when the unmodified sellmeijer piping rule is used. -/// -/// -/// -/// -template -double SteadyStatePwPipingElement::CalculateParticleDiameter(const PropertiesType& Prop) -{ - double diameter; - - if (Prop[PIPE_MODIFIED_D]) diameter = 2.08e-4 * pow((Prop[PIPE_D_70] / 2.08e-4), 0.4); - else diameter = Prop[PIPE_D_70]; - return diameter; -} - /// /// Calculates the equilibrium pipe height of a piping element according to Sellmeijers rule /// @@ -311,7 +295,7 @@ double SteadyStatePwPipingElement::CalculateEquilibriumPipeHeig double dhdx = CalculateHeadGradient(Prop, Geom, pipe_length); // calculate particle diameter - double particle_d = CalculateParticleDiameter(Prop); + double particle_d = GeoTransportEquationUtilities::CalculateParticleDiameter(Prop); // todo calculate slope of pipe, currently pipe is assumed to be horizontal const double pipeSlope = 0; @@ -341,5 +325,4 @@ bool SteadyStatePwPipingElement::InEquilibrium(const Properties template class SteadyStatePwPipingElement<2, 4>; template class SteadyStatePwPipingElement<3, 6>; template class SteadyStatePwPipingElement<3, 8>; - } // Namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.hpp b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.hpp index 2533c6c801f0..b25819675176 100644 --- a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.hpp @@ -119,8 +119,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SteadyStatePwPipingElement std::vector& rValues, const ProcessInfo& rCurrentProcessInfo) override; - double CalculateParticleDiameter(const PropertiesType& Prop); - double pipe_initialised = false; private: diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.hpp b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.hpp index baddcafbf21f..602d9a7cc88e 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.hpp @@ -166,8 +166,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwElement : public UPwSmall void load(Serializer& rSerializer) override{KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, Element)} // Assignment operator. - TransientPwElement& - operator=(TransientPwElement const& rOther); + TransientPwElement& operator=(TransientPwElement const& rOther); // Copy constructor. TransientPwElement(TransientPwElement const& rOther); diff --git a/applications/GeoMechanicsApplication/custom_geometries/line_interface_geometry.h b/applications/GeoMechanicsApplication/custom_geometries/line_interface_geometry.h index f545efec4416..9c9de98e2d5e 100644 --- a/applications/GeoMechanicsApplication/custom_geometries/line_interface_geometry.h +++ b/applications/GeoMechanicsApplication/custom_geometries/line_interface_geometry.h @@ -261,19 +261,22 @@ class LineInterfaceGeometry : public Geometry { const auto points = this->Points(); const auto number_of_midline_nodes = std::size_t{points.size() / 2}; + auto result = PointerVector{number_of_midline_nodes}; auto is_null = [](const auto& rNodePtr) { return rNodePtr == nullptr; }; if (std::any_of(points.ptr_begin(), points.ptr_end(), is_null)) { // At least one point is not defined, so the points of the mid-line can't be computed. // As a result, all the mid-line points will be undefined. - return PointerVector{number_of_midline_nodes}; + return result; } - auto result = PointerVector{}; - for (std::size_t i = 0; i < number_of_midline_nodes; ++i) { - const auto mid_point = (points[i] + points[i + number_of_midline_nodes]) / 2; - result.push_back(make_intrusive(points[i].Id(), mid_point)); - } + auto begin_of_first_side = points.ptr_begin(); + auto begin_of_second_side = begin_of_first_side + number_of_midline_nodes; + auto make_mid_point = [](auto pPoint1, auto pPoint2) { + return make_intrusive(pPoint1->Id(), Point{(*pPoint1 + *pPoint2) / 2}); + }; + std::transform(begin_of_first_side, begin_of_second_side, begin_of_second_side, + result.ptr_begin(), make_mid_point); return result; } diff --git a/applications/GeoMechanicsApplication/custom_processes/README.md b/applications/GeoMechanicsApplication/custom_processes/README.md index 808729de488b..c7d3b681dd0a 100644 --- a/applications/GeoMechanicsApplication/custom_processes/README.md +++ b/applications/GeoMechanicsApplication/custom_processes/README.md @@ -89,11 +89,19 @@ When the overconsolidation ratio ("OCR") and optionally the unloading-reloading $$K_0 = OCR.K_0^{nc} + \frac{\nu_{ur}}{1 - \nu_{ur}} ( OCR - 1 )$$ -![Initial effective stress tensor](initial_effective_stress_tensor.png) +```math +\sigma^{'}_{initial} = \begin{bmatrix} {K_0 \sigma^{'}_{zz}} & 0 & 0 \\ + 0 & {K_0 \sigma^{'}_{zz}} & 0 \\ + 0 & 0 & {\sigma^{'}_{zz}} \end{bmatrix} +``` Alternatively, when the pre-overburden pressure "POP" is specified, the initial stress tensor becomes: -![Initial effective stress tensor with POP](initial_effective_stress_tensor_with_POP.png) +```math +\sigma^{'}_{initial} = \begin{bmatrix} {K_0^{nc} (\sigma^{'}_{zz} + POP )} - \frac{\nu_{ur}}{1 - \nu_{ur}} POP & 0 & 0 \\ + 0 & {K_0^{nc} (\sigma^{'}_{zz} + POP)} - \frac{\nu_{ur}}{1 - \nu_{ur}} POP & 0 \\ + 0 & 0 & {\sigma^{'}_{zz}} \end{bmatrix} +``` When the optional unloading-reloading Poisson's ratio is omitted, a default value $\nu_{ur} = 0$ is used such that the correction term drops to 0. diff --git a/applications/GeoMechanicsApplication/custom_processes/initial_effective_stress_tensor.png b/applications/GeoMechanicsApplication/custom_processes/initial_effective_stress_tensor.png deleted file mode 100644 index 76ef4f9bafcf..000000000000 Binary files a/applications/GeoMechanicsApplication/custom_processes/initial_effective_stress_tensor.png and /dev/null differ diff --git a/applications/GeoMechanicsApplication/custom_processes/initial_effective_stress_tensor_with_POP.png b/applications/GeoMechanicsApplication/custom_processes/initial_effective_stress_tensor_with_POP.png deleted file mode 100644 index 6f9650b372da..000000000000 Binary files a/applications/GeoMechanicsApplication/custom_processes/initial_effective_stress_tensor_with_POP.png and /dev/null differ diff --git a/applications/GeoMechanicsApplication/custom_python/add_custom_strategies_to_python.cpp b/applications/GeoMechanicsApplication/custom_python/add_custom_strategies_to_python.cpp index a985ff46abcf..2e1e5183b4b4 100644 --- a/applications/GeoMechanicsApplication/custom_python/add_custom_strategies_to_python.cpp +++ b/applications/GeoMechanicsApplication/custom_python/add_custom_strategies_to_python.cpp @@ -110,18 +110,18 @@ void AddCustomStrategiesToPython(pybind11::module& m) py::class_( m, "GeoMechanicsNewtonRaphsonStrategy") - .def(py::init()); py::class_( m, "GeoMechanicsNewtonRaphsonErosionProcessStrategy") - .def(py::init()); py::class_( m, "GeoMechanicsRammArcLengthStrategy") - .def(py::init()) .def("UpdateLoads", &GeoMechanicsRammArcLengthStrategyType::UpdateLoads); diff --git a/applications/GeoMechanicsApplication/custom_strategies/strategies/README.md b/applications/GeoMechanicsApplication/custom_strategies/strategies/README.md index b819a0ebaa4b..52febd5a1e65 100644 --- a/applications/GeoMechanicsApplication/custom_strategies/strategies/README.md +++ b/applications/GeoMechanicsApplication/custom_strategies/strategies/README.md @@ -54,6 +54,12 @@ This function performs an iterative process towards certain equilibrium values o the modelling is performed. Then a pipe height is corrected by a given increment for each open piping element. The maximum number of iterations is the user input, "max_piping_iterations". +The typical case is the following. Initially, the equilibrium pipe height is small and quickly the pipe height is +increased above the equilibrium height. At this moment PIPE_EROSION is set to true. Then due to the flow development, +the water pressure gradient decreases, which leads to a higher value of the equilibrium height. As a result the +equilibrium pipe height grows faster than the pipe height. When both of the heights become equal, the iterative process +is stopped. + ![check_pipe_equilibrium.svg](check_pipe_equilibrium.svg) ### Recalculate function diff --git a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp index 2431333c34f5..5e3aa66fc3bd 100644 --- a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp +++ b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp @@ -22,6 +22,7 @@ #include "boost/range/adaptor/filtered.hpp" #include "custom_elements/steady_state_Pw_piping_element.hpp" #include "custom_strategies/strategies/geo_mechanics_newton_raphson_strategy.hpp" +#include "custom_utilities/transport_equation_utilities.hpp" #include "geo_mechanics_application_variables.h" #include @@ -31,7 +32,6 @@ namespace Kratos { - template class GeoMechanicsNewtonRaphsonErosionProcessStrategy : public GeoMechanicsNewtonRaphsonStrategy @@ -46,15 +46,12 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy using DofsArrayType = typename BaseType::DofsArrayType; using TSystemMatrixType = typename BaseType::TSystemMatrixType; using TSystemVectorType = typename BaseType::TSystemVectorType; - using filtered_elements = - typename boost::range_detail::filtered_range, std::vector>; - using PropertiesType = Properties; - using NodeType = Node; - using GeometryType = Geometry; + using PropertiesType = Properties; + using NodeType = Node; + using GeometryType = Geometry; GeoMechanicsNewtonRaphsonErosionProcessStrategy(ModelPart& model_part, typename TSchemeType::Pointer pScheme, - typename TLinearSolver::Pointer pNewLinearSolver, typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria, typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver, Parameters& rParameters, @@ -63,16 +60,7 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy bool ReformDofSetAtEachStep = false, bool MoveMeshFlag = false) : GeoMechanicsNewtonRaphsonStrategy( - model_part, - pScheme, - pNewLinearSolver, - pNewConvergenceCriteria, - pNewBuilderAndSolver, - rParameters, - MaxIterations, - CalculateReactions, - ReformDofSetAtEachStep, - MoveMeshFlag) + model_part, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, rParameters, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag) { rank = model_part.GetCommunicator().MyPID(); mPipingIterations = rParameters["max_piping_iterations"].GetInt(); @@ -87,7 +75,17 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy // get piping elements std::vector PipeElements = GetPipingElements(); - unsigned int n_el = PipeElements.size(); // number of piping elements + + auto piping_interface_elements = std::vector*>{}; + std::transform( + PipeElements.begin(), PipeElements.end(), std::back_inserter(piping_interface_elements), + [](auto p_element) { return dynamic_cast*>(p_element); }); + KRATOS_DEBUG_ERROR_IF(std::any_of(piping_interface_elements.begin(), + piping_interface_elements.end(), [](auto p_element) { + return p_element == nullptr; + })) << "Not all open piping elements could be downcast to SteadyStatePwPipingElement<2, 4>*\n"; + + unsigned int n_el = PipeElements.size(); // number of piping elements // get initially open pipe elements unsigned int openPipeElements = this->InitialiseNumActivePipeElements(PipeElements); @@ -95,7 +93,7 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy if (PipeElements.empty()) { KRATOS_INFO_IF("PipingLoop", this->GetEchoLevel() > 0 && rank == 0) << "No Pipe Elements -> Finalizing Solution " << std::endl; - GeoMechanicsNewtonRaphsonStrategy::FinalizeSolutionStep(); + this->BaseClassFinalizeSolutionStep(); return; } // calculate max pipe height and pipe increment @@ -115,7 +113,7 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy std::function filter = [](Element* i) { return i->Has(PIPE_ACTIVE) && i->GetValue(PIPE_ACTIVE); }; - filtered_elements OpenPipeElements = PipeElements | boost::adaptors::filtered(filter); + auto OpenPipeElements = piping_interface_elements | boost::adaptors::filtered(filter); KRATOS_INFO_IF("PipingLoop", this->GetEchoLevel() > 0 && rank == 0) << "Number of Open Pipe Elements: " << boost::size(OpenPipeElements) << std::endl; @@ -134,13 +132,13 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy save_or_reset_pipe_heights(OpenPipeElements, grow); } // recalculate groundwater flow - bool converged = Recalculate(); + bool converged = this->Recalculate(); // error check KRATOS_ERROR_IF_NOT(converged) << "Groundwater flow calculation failed to converge." << std::endl; } - GeoMechanicsNewtonRaphsonStrategy::FinalizeSolutionStep(); + this->BaseClassFinalizeSolutionStep(); } /** @@ -150,7 +148,6 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy int Check() override { KRATOS_TRY - BaseType::Check(); this->GetBuilderAndSolver()->Check(BaseType::GetModelPart()); this->GetScheme()->Check(BaseType::GetModelPart()); @@ -254,21 +251,6 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy return nOpenElements; } - /// - /// Calculates the pipe particle diameter according to the modified Sellmeijer rule if selected. - /// Else the pipe particle diameter is equal to the d70. - /// - /// - /// - double CalculateParticleDiameter(const PropertiesType& Prop) - { - double diameter; - - if (Prop[PIPE_MODIFIED_D]) diameter = 2.08e-4 * pow((Prop[PIPE_D_70] / 2.08e-4), 0.4); - else diameter = Prop[PIPE_D_70]; - return diameter; - } - /// /// Calculates the maximum pipe height which the algorithm will allow /// @@ -282,8 +264,8 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy // loop over all elements for (Element* pipe_element : pipe_elements) { // calculate pipe particle diameter of pipe element - PropertiesType prop = pipe_element->GetProperties(); - double particle_diameter = this->CalculateParticleDiameter(prop); + PropertiesType prop = pipe_element->GetProperties(); + double particle_diameter = GeoTransportEquationUtilities::CalculateParticleDiameter(prop); // get maximum pipe particle diameter of all pipe elements if (particle_diameter > max_diameter) { @@ -306,26 +288,24 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy return max_pipe_height / (n_steps - 1); } - bool Recalculate() + virtual bool Recalculate() { KRATOS_INFO_IF("ResidualBasedNewtonRaphsonStrategy", this->GetEchoLevel() > 0 && rank == 0) << "Recalculating" << std::endl; - // KRATOS_INFO_IF("PipingLoop") << "Recalculating" << std::endl; - // ModelPart& CurrentModelPart = this->GetModelPart(); - // this->Clear(); - - // Reset displacements to the initial (Assumes Water Pressure is the convergence criteria) - /* block_for_each(CurrentModelPart.Nodes(), [&](Node& rNode) { - auto dold = rNode.GetSolutionStepValue(WATER_PRESSURE, 1); - rNode.GetSolutionStepValue(WATER_PRESSURE, 0) = dold; - });*/ GeoMechanicsNewtonRaphsonStrategy::InitializeSolutionStep(); GeoMechanicsNewtonRaphsonStrategy::Predict(); return GeoMechanicsNewtonRaphsonStrategy::SolveSolutionStep(); } - bool check_pipe_equilibrium(filtered_elements open_pipe_elements, double amax, unsigned int mPipingIterations) + virtual void BaseClassFinalizeSolutionStep() + { + // to override in a unit test + GeoMechanicsNewtonRaphsonStrategy::FinalizeSolutionStep(); + } + + template + bool check_pipe_equilibrium(const FilteredElementType& open_pipe_elements, double amax, unsigned int mPipingIterations) { bool equilibrium = false; bool converged = true; @@ -341,7 +321,7 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy equilibrium = true; // perform a flow calculation and stop growing if the calculation doesn't converge - converged = Recalculate(); + converged = this->Recalculate(); // todo: JDN (20220817) : grow not used. // if (!converged) @@ -353,14 +333,12 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy // Update depth of open piping Elements equilibrium = true; for (auto OpenPipeElement : open_pipe_elements) { - auto pElement = static_cast*>(OpenPipeElement); - // get open pipe element geometry and properties auto& Geom = OpenPipeElement->GetGeometry(); auto& prop = OpenPipeElement->GetProperties(); // calculate equilibrium pipe height and get current pipe height - double eq_height = pElement->CalculateEquilibriumPipeHeight( + double eq_height = OpenPipeElement->CalculateEquilibriumPipeHeight( prop, Geom, OpenPipeElement->GetValue(PIPE_ELEMENT_LENGTH)); double current_height = OpenPipeElement->GetValue(PIPE_HEIGHT); @@ -383,7 +361,6 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy equilibrium = true; OpenPipeElement->SetValue(PIPE_HEIGHT, small_pipe_height); } - // calculate difference between equilibrium height and current pipe height OpenPipeElement->SetValue(DIFF_PIPE_HEIGHT, eq_height - current_height); } @@ -437,17 +414,16 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy /// open pipe elements /// boolean to check if pipe grows /// - void save_or_reset_pipe_heights(filtered_elements open_pipe_elements, bool grow) + template + void save_or_reset_pipe_heights(const FilteredElementType& open_pipe_elements, bool grow) { - for (Element* OpenPipeElement : open_pipe_elements) { + for (auto p_element : open_pipe_elements) { if (grow) { - OpenPipeElement->SetValue(PREV_PIPE_HEIGHT, OpenPipeElement->GetValue(PIPE_HEIGHT)); + p_element->SetValue(PREV_PIPE_HEIGHT, p_element->GetValue(PIPE_HEIGHT)); } else { - OpenPipeElement->SetValue(PIPE_HEIGHT, OpenPipeElement->GetValue(PREV_PIPE_HEIGHT)); + p_element->SetValue(PIPE_HEIGHT, p_element->GetValue(PREV_PIPE_HEIGHT)); } } } - -}; // Class GeoMechanicsNewtonRaphsonStrategy - -} // namespace Kratos \ No newline at end of file +}; // Class GeoMechanicsNewtonRaphsonErosionProcessStrategy +} // namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_strategy.hpp b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_strategy.hpp index f8360def0513..c9017fe3965d 100644 --- a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_strategy.hpp +++ b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_strategy.hpp @@ -52,16 +52,14 @@ class GeoMechanicsNewtonRaphsonStrategy * @brief Default constructor * @param rModelPart The model part of the problem * @param pScheme The integration scheme - * @param pNewLinearSolver The linear solver employed * @param pNewConvergenceCriteria The convergence criteria employed * @param MaxIterations The maximum number of iterations * @param CalculateReactions The flag for the reaction calculation * @param ReformDofSetAtEachStep The flag that allows to compute the modification of the DOF * @param MoveMeshFlag The flag that allows to move the mesh */ - GeoMechanicsNewtonRaphsonStrategy(ModelPart& model_part, - typename TSchemeType::Pointer pScheme, - typename TLinearSolver::Pointer pNewLinearSolver, + GeoMechanicsNewtonRaphsonStrategy(ModelPart& model_part, + typename TSchemeType::Pointer pScheme, typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria, typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver, Parameters& rParameters, @@ -70,15 +68,7 @@ class GeoMechanicsNewtonRaphsonStrategy bool ReformDofSetAtEachStep = false, bool MoveMeshFlag = false) : ResidualBasedNewtonRaphsonStrategy( - model_part, - pScheme, - /*pNewLinearSolver,*/ - pNewConvergenceCriteria, - pNewBuilderAndSolver, - MaxIterations, - CalculateReactions, - ReformDofSetAtEachStep, - MoveMeshFlag) + model_part, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag) { // only include validation with c++11 since raw_literals do not exist in c++03 Parameters default_parameters(R"( @@ -100,62 +90,6 @@ class GeoMechanicsNewtonRaphsonStrategy // Validate against defaults -- this also ensures no type mismatch rParameters.ValidateAndAssignDefaults(default_parameters); - - // Set Load SubModelParts and Variable names - if (rParameters["loads_sub_model_part_list"].size() > 0) { - mSubModelPartList.resize(rParameters["loads_sub_model_part_list"].size()); - mVariableNames.resize(rParameters["loads_variable_list"].size()); - - if (mSubModelPartList.size() != mVariableNames.size()) - KRATOS_ERROR << "For each SubModelPart there must be a corresponding nodal Variable" - << std::endl; - - for (unsigned int i = 0; i < mVariableNames.size(); ++i) { - mSubModelPartList[i] = - &(model_part.GetSubModelPart(rParameters["loads_sub_model_part_list"][i].GetString())); - mVariableNames[i] = rParameters["loads_variable_list"][i].GetString(); - } - } - } - -protected: - /// Member Variables - std::vector mSubModelPartList; /// List of every SubModelPart associated to an external load - std::vector mVariableNames; /// Name of the nodal variable associated to every SubModelPart - - int Check() override - { - KRATOS_TRY - - return MotherType::Check(); - - KRATOS_CATCH("") - } - - double CalculateReferenceDofsNorm(DofsArrayType& rDofSet) - { - double ReferenceDofsNorm = 0.0; - - int NumThreads = ParallelUtilities::GetNumThreads(); - OpenMPUtils::PartitionVector DofSetPartition; - OpenMPUtils::DivideInPartitions(rDofSet.size(), NumThreads, DofSetPartition); - -#pragma omp parallel reduction(+ : ReferenceDofsNorm) - { - int k = OpenMPUtils::ThisThread(); - - typename DofsArrayType::iterator DofsBegin = rDofSet.begin() + DofSetPartition[k]; - typename DofsArrayType::iterator DofsEnd = rDofSet.begin() + DofSetPartition[k + 1]; - - for (typename DofsArrayType::iterator itDof = DofsBegin; itDof != DofsEnd; ++itDof) { - if (itDof->IsFree()) { - const double& temp = itDof->GetSolutionStepValue(); - ReferenceDofsNorm += temp * temp; - } - } - } - - return sqrt(ReferenceDofsNorm); } }; // Class GeoMechanicsNewtonRaphsonStrategy diff --git a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_ramm_arc_length_strategy.hpp b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_ramm_arc_length_strategy.hpp index 2ecd12520b26..e13287c0bcaf 100644 --- a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_ramm_arc_length_strategy.hpp +++ b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_ramm_arc_length_strategy.hpp @@ -19,6 +19,9 @@ // Application includes #include "geo_mechanics_application_variables.h" +#include +#include + namespace Kratos { @@ -50,12 +53,9 @@ class GeoMechanicsRammArcLengthStrategy using GrandMotherType::mpDx; // Delta x of iteration i using GrandMotherType::mpScheme; using GrandMotherType::mReformDofSetAtEachStep; - using MotherType::mSubModelPartList; - using MotherType::mVariableNames; - GeoMechanicsRammArcLengthStrategy(ModelPart& model_part, - typename TSchemeType::Pointer pScheme, - typename TLinearSolver::Pointer pNewLinearSolver, + GeoMechanicsRammArcLengthStrategy(ModelPart& model_part, + typename TSchemeType::Pointer pScheme, typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria, typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver, Parameters& rParameters, @@ -64,17 +64,24 @@ class GeoMechanicsRammArcLengthStrategy bool ReformDofSetAtEachStep = false, bool MoveMeshFlag = false) : GeoMechanicsNewtonRaphsonStrategy( - model_part, - pScheme, - pNewLinearSolver, - pNewConvergenceCriteria, - pNewBuilderAndSolver, - rParameters, - MaxIterations, - CalculateReactions, - ReformDofSetAtEachStep, - MoveMeshFlag) + model_part, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, rParameters, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag) { + // Set Load SubModelParts and Variable names + if (rParameters["loads_sub_model_part_list"].size() > 0) { + mSubModelPartList.resize(rParameters["loads_sub_model_part_list"].size()); + mVariableNames.resize(rParameters["loads_variable_list"].size()); + + if (mSubModelPartList.size() != mVariableNames.size()) + KRATOS_ERROR << "For each SubModelPart there must be a corresponding nodal Variable" + << std::endl; + + for (unsigned int i = 0; i < mVariableNames.size(); ++i) { + mSubModelPartList[i] = + &(model_part.GetSubModelPart(rParameters["loads_sub_model_part_list"][i].GetString())); + mVariableNames[i] = rParameters["loads_variable_list"][i].GetString(); + } + } + mDesiredIterations = rParameters["desired_iterations"].GetInt(); mMaxRadiusFactor = rParameters["max_radius_factor"].GetDouble(); mMinRadiusFactor = rParameters["min_radius_factor"].GetDouble(); @@ -306,7 +313,7 @@ class GeoMechanicsRammArcLengthStrategy else if (mRadius < mMinRadiusFactor * mRadius_0) mRadius = mMinRadiusFactor * mRadius_0; // Update Norm of x - mNormxEquilibrium = this->CalculateReferenceDofsNorm(rDofSet); + mNormxEquilibrium = CalculateReferenceDofsNorm(rDofSet); } else { std::cout << "************ NO CONVERGENCE: restoring equilibrium path ************" << std::endl; @@ -534,6 +541,23 @@ class GeoMechanicsRammArcLengthStrategy // Save the applied Lambda factor mLambda_old = mLambda; } + +private: + static double CalculateReferenceDofsNorm(DofsArrayType& rDofSet) + { + auto is_free_dof = [](const auto& rDof) { return rDof.IsFree(); }; + auto free_dofs = rDofSet | boost::adaptors::filtered(is_free_dof); + + auto free_dof_values = Vector{boost::size(free_dofs)}; + auto get_dof_value = [](const auto& rDof) { return rDof.GetSolutionStepValue(); }; + std::transform(std::begin(free_dofs), std::end(free_dofs), free_dof_values.begin(), get_dof_value); + + return norm_2(free_dof_values); + } + + std::vector mSubModelPartList; // List of every SubModelPart associated to an external load + std::vector mVariableNames; // Name of the nodal variable associated to every SubModelPart + }; // Class GeoMechanicsRammArcLengthStrategy } // namespace Kratos \ No newline at end of file diff --git a/applications/GeoMechanicsApplication/custom_utilities/solving_strategy_factory.hpp b/applications/GeoMechanicsApplication/custom_utilities/solving_strategy_factory.hpp index 03296529d279..6e21ce30c7c4 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/solving_strategy_factory.hpp +++ b/applications/GeoMechanicsApplication/custom_utilities/solving_strategy_factory.hpp @@ -69,7 +69,7 @@ class SolvingStrategyFactory ParametersUtilities::CopyOptionalParameters(rSolverSettings, strategy_entries); auto result = std::make_unique>( - rModelPart, scheme, solver, criteria, builder_and_solver, strategy_parameters, + rModelPart, scheme, criteria, builder_and_solver, strategy_parameters, max_iterations, calculate_reactions, reform_dof_set_at_each_step, move_mesh_flag); result->SetEchoLevel(echo_level); return result; diff --git a/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.cpp b/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.cpp new file mode 100644 index 000000000000..a6857bc732b0 --- /dev/null +++ b/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.cpp @@ -0,0 +1,25 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Anne van de Graaf +// + +#include "transport_equation_utilities.hpp" + +namespace Kratos +{ + +double GeoTransportEquationUtilities::CalculateParticleDiameter(const Properties& rProperties) +{ + return rProperties.Has(PIPE_MODIFIED_D) && rProperties[PIPE_MODIFIED_D] + ? 2.08e-4 * std::pow((rProperties[PIPE_D_70] / 2.08e-4), 0.4) + : rProperties[PIPE_D_70]; +} + +} // namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp b/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp index ed2d630fb5d0..069cb4e949f4 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp +++ b/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp @@ -18,12 +18,13 @@ #include "custom_retention/retention_law.h" #include "custom_utilities/stress_strain_utilities.h" #include "geo_mechanics_application_variables.h" +#include "includes/kratos_export_api.h" #include "includes/variables.h" namespace Kratos { -class GeoTransportEquationUtilities +class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoTransportEquationUtilities { public: template @@ -156,6 +157,8 @@ class GeoTransportEquationUtilities return result; } + [[nodiscard]] static double CalculateParticleDiameter(const Properties& rProperties); + private: [[nodiscard]] static double CalculateBiotCoefficient(const Matrix& rConstitutiveMatrix, const Properties& rProperties) diff --git a/applications/GeoMechanicsApplication/custom_workflows/dgeoflow.cpp b/applications/GeoMechanicsApplication/custom_workflows/dgeoflow.cpp index dbc7c4350839..e0abc0639d35 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/dgeoflow.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/dgeoflow.cpp @@ -138,8 +138,8 @@ KratosExecute::GeoMechanicsNewtonRaphsonErosionProcessStrategyType::Pointer Krat auto pSolvingStrategy = Kratos::make_unique>( - rModelPart, p_scheme, p_solver, p_criteria, p_builder_and_solver, p_parameters, - MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag); + rModelPart, p_scheme, p_criteria, p_builder_and_solver, p_parameters, MaxIterations, + CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag); pSolvingStrategy->Check(); return pSolvingStrategy; diff --git a/applications/GeoMechanicsApplication/python_scripts/geomechanics_solver.py b/applications/GeoMechanicsApplication/python_scripts/geomechanics_solver.py index c5580cb9f167..4a0d43e24856 100644 --- a/applications/GeoMechanicsApplication/python_scripts/geomechanics_solver.py +++ b/applications/GeoMechanicsApplication/python_scripts/geomechanics_solver.py @@ -489,7 +489,6 @@ def _ConstructSolver(self, builder_and_solver, strategy_type): self.strategy_params.AddValue("loads_variable_list",self.settings["loads_variable_list"]) solving_strategy = GeoMechanicsApplication.GeoMechanicsNewtonRaphsonStrategy(self.computing_model_part, self.scheme, - self.linear_solver, self.convergence_criterion, builder_and_solver, self.strategy_params, @@ -504,7 +503,6 @@ def _ConstructSolver(self, builder_and_solver, strategy_type): self.strategy_params.AddValue("max_piping_iterations", self.settings["max_piping_iterations"]) solving_strategy = GeoMechanicsApplication.GeoMechanicsNewtonRaphsonErosionProcessStrategy(self.computing_model_part, self.scheme, - self.linear_solver, self.convergence_criterion, builder_and_solver, self.strategy_params, @@ -545,7 +543,6 @@ def _ConstructSolver(self, builder_and_solver, strategy_type): self.strategy_params.AddValue("loads_variable_list",self.settings["loads_variable_list"]) solving_strategy = GeoMechanicsApplication.GeoMechanicsRammArcLengthStrategy(self.computing_model_part, self.scheme, - self.linear_solver, self.convergence_criterion, builder_and_solver, self.strategy_params, diff --git a/applications/GeoMechanicsApplication/tests/test_parameter_field/parameter_field_python_umat_parameters/custom_field_umat_parameters.py b/applications/GeoMechanicsApplication/python_scripts/user_defined_scripts/custom_field_umat_parameters.py similarity index 100% rename from applications/GeoMechanicsApplication/tests/test_parameter_field/parameter_field_python_umat_parameters/custom_field_umat_parameters.py rename to applications/GeoMechanicsApplication/python_scripts/user_defined_scripts/custom_field_umat_parameters.py diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_dgeosettlement.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_dgeosettlement.cpp index cc94499c9f60..5e1c9a5c74e1 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/test_dgeosettlement.cpp +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_dgeosettlement.cpp @@ -43,10 +43,10 @@ const std::string parameter_json_settings = R"( int RunStage(KratosGeoSettlement& rSettlement) { - return rSettlement.RunStage( - "", "", [](const char*) { /* kept empty as a stub method */ }, - [](const double) { /* kept empty as a stub method */ }, - [](const char*) { /* kept empty as a stub method */ }, []() { return false; }); + return rSettlement.RunStage("", "", [](const char*) { /* kept empty as a stub method */ }, + [](const double) { /* kept empty as a stub method */ }, + [](const char*) { /* kept empty as a stub method */ }, + []() { return false; }); } void ExpectNumberOfReadCallsIsEqualToOne(const KratosGeoSettlement& rSettlement) diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_geo_mechanics_newton_rapson_erosion_processs_strategy.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_geo_mechanics_newton_rapson_erosion_processs_strategy.cpp new file mode 100644 index 000000000000..bbe6e65e0697 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_geo_mechanics_newton_rapson_erosion_processs_strategy.cpp @@ -0,0 +1,216 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Gennady Markelov +// + +// Project includes +#include "custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp" +#include "geo_mechanics_fast_suite.h" +#include "test_utilities.h" + +#include +#include +#include + +namespace Kratos +{ +template +class MockGeoMechanicsNewtonRaphsonErosionProcessStrategy + : public GeoMechanicsNewtonRaphsonErosionProcessStrategy +{ +public: + using SolvingStrategyType = ImplicitSolvingStrategy; + using TConvergenceCriteriaType = ConvergenceCriteria; + using TBuilderAndSolverType = typename SolvingStrategyType::TBuilderAndSolverType; + using TSchemeType = typename SolvingStrategyType::TSchemeType; + + MockGeoMechanicsNewtonRaphsonErosionProcessStrategy(ModelPart& rModelPart, + typename TSchemeType::Pointer pScheme, + typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria, + typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver, + Parameters& rParameters) + : GeoMechanicsNewtonRaphsonErosionProcessStrategy( + rModelPart, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, rParameters), + mrModelPart(rModelPart) + { + } + +private: + bool Recalculate() override + { + // The function provides a pressure change over iterations that leads to a non-linear growth of an equilibrium pipe height. + // As a result the height curve crosses twice a pipe height growth curve in the search algorithm implemented + // in GeoMechanicsNewtonRaphsonErosionProcessStrategy. The search stops at the second point. + for (auto& node : mrModelPart.Nodes()) { + node.FastGetSolutionStepValue(WATER_PRESSURE) = + std::pow(node.FastGetSolutionStepValue(WATER_PRESSURE), 0.9725); + } + return true; + } + + void BaseClassFinalizeSolutionStep() override + { + // to avoid calling the BaseClassFinalizeSolutionStep in tested class. + } + + ModelPart& mrModelPart; +}; + +Node::Pointer CreateNodeWithSolutionStepValues(ModelPart& rModelPart, int Id, double CoordinateX, double CoordinateY, double WaterPressure) +{ + constexpr double coordinate_z = 0.0; + auto p_node = rModelPart.CreateNewNode(Id, CoordinateX, CoordinateY, coordinate_z); + + p_node->FastGetSolutionStepValue(WATER_PRESSURE) = WaterPressure; + const array_1d GravitationalAcceleration{0, -9.81, 0}; + p_node->FastGetSolutionStepValue(VOLUME_ACCELERATION) = GravitationalAcceleration; + + return p_node; +} + +Geometry::Pointer CreateQuadrilateral2D4N(ModelPart& rModelPart, + const std::vector& rNodeIds, + double Xmin, + double Xmax, + double WaterPressureLeft, + double WaterPressureRight) +{ + constexpr double y_min = -0.1; + constexpr double y_max = 0.1; + + return make_shared>( + CreateNodeWithSolutionStepValues(rModelPart, rNodeIds[0], Xmin, y_min, WaterPressureLeft), + CreateNodeWithSolutionStepValues(rModelPart, rNodeIds[1], Xmax, y_min, WaterPressureRight), + CreateNodeWithSolutionStepValues(rModelPart, rNodeIds[2], Xmax, y_max, WaterPressureRight), + CreateNodeWithSolutionStepValues(rModelPart, rNodeIds[3], Xmin, y_max, WaterPressureLeft)); +} + +auto SetupPipingStrategy(Model& rModel) +{ + using SparseSpaceType = UblasSpace; + using LocalSpaceType = UblasSpace; + using LinearSolverType = SkylineLUFactorizationSolver; + using ConvergenceCriteriaType = ConvergenceCriteria; + using MixedGenericCriteriaType = MixedGenericCriteria; + using ConvergenceVariableListType = MixedGenericCriteriaType::ConvergenceVariableListType; + + auto& r_model_part = rModel.CreateModelPart("ModelPart", 1); + r_model_part.AddNodalSolutionStepVariable(WATER_PRESSURE); + r_model_part.AddNodalSolutionStepVariable(VOLUME_ACCELERATION); + const auto& r_process_info = r_model_part.GetProcessInfo(); + + auto p_element_props = r_model_part.CreateNewProperties(0); + p_element_props->SetValue(CONSTITUTIVE_LAW, LinearElastic2DInterfaceLaw().Clone()); + + auto p_element = make_intrusive>( + 0, CreateQuadrilateral2D4N(r_model_part, std::vector{13, 14, 15, 16}, 3.0, 4.0, 2000.0, 2000.0), + p_element_props, std::make_unique()); + p_element->Initialize(r_process_info); + r_model_part.AddElement(p_element); + + // Set the pipe element properties + auto p_piping_element_props = r_model_part.CreateNewProperties(1); + p_piping_element_props->SetValue(PIPE_D_70, 2e-4); + p_piping_element_props->SetValue(PIPE_ETA, 0.25); + p_piping_element_props->SetValue(PIPE_THETA, 30); + p_piping_element_props->SetValue(DENSITY_SOLID, 2650); + p_piping_element_props->SetValue(DENSITY_WATER, 1000); + p_piping_element_props->SetValue(PIPE_ELEMENT_LENGTH, 1.0); + p_piping_element_props->SetValue(PIPE_MODIFIED_D, false); + p_piping_element_props->SetValue(PIPE_MODEL_FACTOR, 1); + p_piping_element_props->SetValue(PIPE_START_ELEMENT, 1); + p_piping_element_props->SetValue(CONSTITUTIVE_LAW, LinearElastic2DInterfaceLaw().Clone()); + + // Create the start piping element + auto p_piping_element = make_intrusive>( + 1, CreateQuadrilateral2D4N(r_model_part, std::vector{1, 2, 3, 4}, 0.0, 1.0, 0.0, 500.0), + p_piping_element_props, std::make_unique()); + p_piping_element->Initialize(r_process_info); + r_model_part.AddElement(p_piping_element); + + // Create other piping elements + p_piping_element = make_intrusive>( + 2, CreateQuadrilateral2D4N(r_model_part, std::vector{5, 6, 7, 8}, 2.0, 3.0, 500.0, 1000.0), + p_piping_element_props, std::make_unique()); + p_piping_element->Initialize(r_process_info); + r_model_part.AddElement(p_piping_element); + + p_piping_element = make_intrusive>( + 3, CreateQuadrilateral2D4N(r_model_part, std::vector{9, 10, 11, 12}, 1.0, 2.0, 1000.0, 1500.0), + p_piping_element_props, std::make_unique()); + p_piping_element->Initialize(r_process_info); + r_model_part.AddElement(p_piping_element); + + Parameters p_parameters(R"( + { + "max_piping_iterations": 25 + } )"); + + const auto p_builder_and_solver = + Kratos::make_shared>(nullptr); + const ConvergenceVariableListType convergence_settings{}; + const ConvergenceCriteriaType::Pointer p_criteria = + std::make_shared(convergence_settings); + + return std::make_unique>( + r_model_part, nullptr, p_criteria, p_builder_and_solver, p_parameters); +} +} // namespace Kratos + +namespace Kratos::Testing +{ +KRATOS_TEST_CASE_IN_SUITE(CheckGetPipingElements, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + // Arrange + Model model; + auto p_solving_strategy = SetupPipingStrategy(model); + + // Act + const auto piping_elements = p_solving_strategy->GetPipingElements(); + + // Assert + const auto elements = model.GetModelPart("ModelPart").Elements(); + constexpr int number_of_elements = 4; + KRATOS_EXPECT_EQ(elements.size(), number_of_elements); + + constexpr int number_of_piping_elements = 3; + KRATOS_EXPECT_EQ(piping_elements.size(), number_of_piping_elements); + KRATOS_EXPECT_EQ(piping_elements[0]->GetId(), 1); + KRATOS_EXPECT_EQ(piping_elements[1]->GetId(), 3); + KRATOS_EXPECT_EQ(piping_elements[2]->GetId(), 2); +} + +KRATOS_TEST_CASE_IN_SUITE(CheckFinalizeSolutionStep, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + // Arrange + Model model; + auto p_solving_strategy = SetupPipingStrategy(model); + + // Act + p_solving_strategy->FinalizeSolutionStep(); + + // Assert + auto elements = model.GetModelPart("ModelPart").Elements(); + + constexpr int active = 1; + KRATOS_EXPECT_EQ(elements[1].GetValue(PIPE_ACTIVE), active); + + constexpr double expected_active_pipe_height = 0.0183333; + KRATOS_EXPECT_NEAR(elements[1].GetValue(PIPE_HEIGHT), expected_active_pipe_height, Defaults::relative_tolerance); + + constexpr int inactive = 0; + KRATOS_EXPECT_EQ(elements[2].GetValue(PIPE_ACTIVE), inactive); + KRATOS_EXPECT_EQ(elements[3].GetValue(PIPE_ACTIVE), inactive); + + constexpr double expected_inactive_pipe_height = 0.0; + KRATOS_EXPECT_NEAR(elements[2].GetValue(PIPE_HEIGHT), expected_inactive_pipe_height, Defaults::relative_tolerance); + KRATOS_EXPECT_NEAR(elements[3].GetValue(PIPE_HEIGHT), expected_inactive_pipe_height, Defaults::relative_tolerance); +} +} // namespace Kratos::Testing \ No newline at end of file diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_time_stepping.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_time_stepping.cpp index 4eeb1fd3a039..a4eaf1a63d46 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/test_time_stepping.cpp +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_time_stepping.cpp @@ -74,22 +74,22 @@ class DummyStrategyWrapper : public StrategyWrapper // intentionally empty } - void CloneTimeStep() override{ + void CloneTimeStep() override { // intentionally empty }; - void RestorePositionsAndDOFVectorToStartOfStep() override{ + void RestorePositionsAndDOFVectorToStartOfStep() override { // intentionally empty }; - void SaveTotalDisplacementFieldAtStartOfTimeLoop() override{ + void SaveTotalDisplacementFieldAtStartOfTimeLoop() override { // intentionally empty }; - void AccumulateTotalDisplacementField() override{ + void AccumulateTotalDisplacementField() override { // intentionally empty }; - void ComputeIncrementalDisplacementField() override{ + void ComputeIncrementalDisplacementField() override { // intentionally empty }; - void OutputProcess() override{ + void OutputProcess() override { // intentionally empty }; diff --git a/applications/GeoMechanicsApplication/tests/test_parameter_field.py b/applications/GeoMechanicsApplication/tests/test_parameter_field.py index 09f3036e3bc7..06ada66ff958 100644 --- a/applications/GeoMechanicsApplication/tests/test_parameter_field.py +++ b/applications/GeoMechanicsApplication/tests/test_parameter_field.py @@ -61,40 +61,8 @@ def test_parameter_field_with_python_umat_parameters(self): test_name = os.path.join("test_parameter_field", "parameter_field_python_umat_parameters") file_path = test_helper.get_file_path(test_name) - custom_script_name = "custom_field_umat_parameters.py" - custom_python_file = os.path.join(file_path, custom_script_name) - if not os.path.isfile(custom_python_file): - raise RuntimeError(f"Source file does not exist. {custom_python_file}") - - # copy user defined python script to installation folder - new_custom_script_path = os.path.join(os.path.dirname(KratosGeo.__file__), "user_defined_scripts") - print(f"Source file path: {custom_python_file}") - print(f"Destination file path: {new_custom_script_path}") - try: - shutil.copy(custom_python_file, new_custom_script_path) - except shutil.SameFileError as e: - raise RuntimeError(f"Source and destination represents the same file.") from e - except PermissionError as e: - raise RuntimeError(f"Permission denied.") from e - except IOError as e: - print(f"Try to change the destination file permission {e}") - os.chmod(new_custom_script_path, stat.S_IRWXU) - shutil.copy(custom_python_file, new_custom_script_path) - except Exception as e: - raise RuntimeError(f"Error occurred while copying the file.") from e - except BaseException as e: - raise RuntimeError(f"An unknown error occurred while copying the file.") from e - - custom_python_file_new_location = os.path.join(new_custom_script_path, custom_script_name) - - if not os.path.isfile(custom_python_file_new_location): - raise RuntimeError(f"File {custom_python_file_new_location} does not exist after copying the source.") - simulation = test_helper.run_kratos(file_path) - if not os.path.isfile(custom_python_file_new_location): - raise RuntimeError(f"File {custom_python_file_new_location} does not exist after run_kratos.") - # get element centers elements = simulation._list_of_output_processes[0].model_part.Elements center_coords = [element.GetGeometry().Center() for element in elements] @@ -107,12 +75,6 @@ def test_parameter_field_with_python_umat_parameters(self): expected_res = 20000 * center_coord[0] + 30000 * center_coord[1] self.assertAlmostEqual(expected_res, res[0][0]) - if not os.path.isfile(custom_python_file_new_location): - raise RuntimeError(f"File {custom_python_file_new_location} does not exist before removal.") - - # remove user defined python script from installation folder - os.remove(os.path.join(new_custom_script_path, custom_script_name)) - def test_parameter_field_with_json_umat_parameters(self): """ Test to check if values from a parameter field stored in a json file are correctly added to diff --git a/applications/OptimizationApplication/python_scripts/optimization_analysis.py b/applications/OptimizationApplication/python_scripts/optimization_analysis.py index 43d703f4c8fc..e772cc297313 100644 --- a/applications/OptimizationApplication/python_scripts/optimization_analysis.py +++ b/applications/OptimizationApplication/python_scripts/optimization_analysis.py @@ -78,7 +78,7 @@ def _CreateModelPartControllers(self): "type": "mdpa_model_part_controller", "module": "KratosMultiphysics.OptimizationApplication.model_part_controllers" }""") - for model_part_controller_settings in self.project_parameters["model_parts"]: + for model_part_controller_settings in self.project_parameters["model_parts"].values(): model_part_controller_settings.AddMissingParameters(default_settings) model_part_controller: ModelPartController = OptimizationComponentFactory(self.model, model_part_controller_settings, self.optimization_problem) self.__list_of_model_part_controllers.append(model_part_controller) @@ -87,7 +87,7 @@ def _CreateAnalyses(self): default_settings = Kratos.Parameters("""{ "module": "KratosMultiphysics.OptimizationApplication.execution_policies" }""") - for analyses_settings in self.project_parameters["analyses"]: + for analyses_settings in self.project_parameters["analyses"].values(): analyses_settings.AddMissingParameters(default_settings) execution_policy = OptimizationComponentFactory(self.model, analyses_settings, self.optimization_problem) self.optimization_problem.AddComponent(execution_policy) @@ -96,7 +96,7 @@ def _CreateResponses(self): default_settings = Kratos.Parameters("""{ "module": "KratosMultiphysics.OptimizationApplication.responses" }""") - for response_settings in self.project_parameters["responses"]: + for response_settings in self.project_parameters["responses"].values(): response_settings.AddMissingParameters(default_settings) response_function: ResponseFunction = OptimizationComponentFactory(self.model, response_settings, self.optimization_problem) self.optimization_problem.AddComponent(response_function) @@ -105,7 +105,7 @@ def _CreateControls(self): default_settings = Kratos.Parameters("""{ "module" : "KratosMultiphysics.OptimizationApplication.controls" }""") - for control_settings in self.project_parameters["controls"]: + for control_settings in self.project_parameters["controls"].values(): control_settings.AddMissingParameters(default_settings) control = OptimizationComponentFactory(self.model, control_settings, self.optimization_problem) self.optimization_problem.AddComponent(control) @@ -129,7 +129,7 @@ def _CreateProcesses(self): for process in factory.ConstructListOfProcesses(kratos_processes[process_type]): self.optimization_problem.AddProcess(process_type, process) if optimization_data_processes.Has(process_type): - for process_settings in optimization_data_processes[process_type]: + for process_settings in optimization_data_processes[process_type].values(): process_settings.AddMissingParameters(optimization_data_process_default_settings) process = OptimizationComponentFactory(self.model, process_settings, self.optimization_problem) self.optimization_problem.AddProcess(process_type, process) diff --git a/applications/SystemIdentificationApplication/tests/test_system_identification.py b/applications/SystemIdentificationApplication/tests/test_system_identification.py index a8a16e956794..1ea3f43da52c 100644 --- a/applications/SystemIdentificationApplication/tests/test_system_identification.py +++ b/applications/SystemIdentificationApplication/tests/test_system_identification.py @@ -33,8 +33,8 @@ def test_SystemIdentification(self): "output_file_name" : "auxiliary_files/summary.csv", "remove_output_file" : true, "comparison_type" : "csv_file", - "tolerance" : 1e-5, - "relative_tolerance" : 1e-6, + "tolerance" : 1e-2, + "relative_tolerance" : 1e-3, "dimension" : 3 }""") CompareTwoFilesCheckProcess(params).Execute() @@ -52,8 +52,8 @@ def test_SystemIdentificationPNorm(self): "output_file_name" : "auxiliary_files/summary_p_norm.csv", "remove_output_file" : true, "comparison_type" : "csv_file", - "tolerance" : 1e-5, - "relative_tolerance" : 1e-6, + "tolerance" : 1e-2, + "relative_tolerance" : 1e-3, "dimension" : 3 }""") CompareTwoFilesCheckProcess(params).Execute() diff --git a/kratos/python/add_model_part_to_python.cpp b/kratos/python/add_model_part_to_python.cpp index c0af5931f51d..2d838a5488d3 100644 --- a/kratos/python/add_model_part_to_python.cpp +++ b/kratos/python/add_model_part_to_python.cpp @@ -744,7 +744,7 @@ void AddModelPartToPython(pybind11::module& m) MapInterface().CreateInterface(m,"GeometriesMapType"); PointerVectorSetPythonInterface().CreateInterface(m,"MasterSlaveConstraintsArray"); - py::class_, DataValueContainer, Flags >(m,"ModelPart") + py::class_, DataValueContainer, Flags>(m, "ModelPart") .def_property("Name", GetModelPartName, SetModelPartName) .def("FullName", &ModelPart::FullName) // .def_property("ProcessInfo", GetProcessInfo, SetProcessInfo) @@ -797,18 +797,30 @@ void AddModelPartToPython(pybind11::module& m) .def("NumberOfTables", &ModelPart::NumberOfTables) .def("AddTable", &ModelPart::AddTable) .def("GetTable", &ModelPart::pGetTable) - .def("HasProperties", [](ModelPart& rSelf, int Id) {return rSelf.HasProperties(Id);}) - .def("HasProperties", [](ModelPart& rSelf, const std::string& rAddress) {return rSelf.HasProperties(rAddress);}) - .def("RecursivelyHasProperties", [](ModelPart& rSelf, int Id) {return rSelf.RecursivelyHasProperties(Id);}) - .def("CreateNewProperties", [](ModelPart& rSelf, int Id) {return rSelf.CreateNewProperties(Id);}) - .def("GetProperties", [](ModelPart& rSelf, int Id) {return rSelf.pGetProperties(Id);}) - .def("GetProperties", [](ModelPart& rSelf, const std::string& rAddress) {return rSelf.pGetProperties(rAddress);}) - .def("GetProperties", [](ModelPart& rSelf) {return rSelf.pProperties();}) - .def("AddProperties", [](ModelPart& rSelf, Properties::Pointer pProperties) {rSelf.AddProperties(pProperties);}) - .def("RemoveProperties", [](ModelPart& rSelf, int Id) {return rSelf.RemoveProperties(Id);}) - .def("RemoveProperties", [](ModelPart& rSelf, Properties::Pointer pProperties) {return rSelf.RemoveProperties(pProperties);}) - .def("RemovePropertiesFromAllLevels", [](ModelPart& rSelf, int Id) {return rSelf.RemovePropertiesFromAllLevels(Id);}) - .def("RemovePropertiesFromAllLevels", [](ModelPart& rSelf, Properties::Pointer pProperties) {return rSelf.RemovePropertiesFromAllLevels(pProperties);}) + .def("HasProperties", [](ModelPart &rSelf, int Id) + { return rSelf.HasProperties(Id); }) + .def("HasProperties", [](ModelPart &rSelf, const std::string &rAddress) + { return rSelf.HasProperties(rAddress); }) + .def("RecursivelyHasProperties", [](ModelPart &rSelf, int Id) + { return rSelf.RecursivelyHasProperties(Id); }) + .def("CreateNewProperties", [](ModelPart &rSelf, int Id) + { return rSelf.CreateNewProperties(Id); }) + .def("GetProperties", [](ModelPart &rSelf, int Id) + { return rSelf.pGetProperties(Id); }) + .def("GetProperties", [](ModelPart &rSelf, const std::string &rAddress) + { return rSelf.pGetProperties(rAddress); }) + .def("GetProperties", [](ModelPart &rSelf) + { return rSelf.pProperties(); }) + .def("AddProperties", [](ModelPart &rSelf, Properties::Pointer pProperties) + { rSelf.AddProperties(pProperties); }) + .def("RemoveProperties", [](ModelPart &rSelf, int Id) + { return rSelf.RemoveProperties(Id); }) + .def("RemoveProperties", [](ModelPart &rSelf, Properties::Pointer pProperties) + { return rSelf.RemoveProperties(pProperties); }) + .def("RemovePropertiesFromAllLevels", [](ModelPart &rSelf, int Id) + { return rSelf.RemovePropertiesFromAllLevels(Id); }) + .def("RemovePropertiesFromAllLevels", [](ModelPart &rSelf, Properties::Pointer pProperties) + { return rSelf.RemovePropertiesFromAllLevels(pProperties); }) .def_property("Properties", ModelPartGetPropertiesContainer, ModelPartSetPropertiesContainer) .def("SetProperties", ModelPartSetPropertiesContainer) .def("PropertiesArray", &ModelPart::PropertiesArray, py::return_value_policy::reference_internal) @@ -864,12 +876,12 @@ void AddModelPartToPython(pybind11::module& m) .def("RemoveGeometry", ModelPartRemoveGeometry2) .def("RemoveGeometryFromAllLevels", ModelPartRemoveGeometryFromAllLevels1) .def("RemoveGeometryFromAllLevels", ModelPartRemoveGeometryFromAllLevels2) - .def_property("Geometries", [](ModelPart& self) { return self.Geometries(); }, - [](ModelPart& self, ModelPart::GeometriesMapType& geometries) { - KRATOS_ERROR << "Setting geometries is not allowed! Trying to set value of ModelPart::Geometries."; }) + .def_property("Geometries", [](ModelPart &self) + { return self.Geometries(); }, [](ModelPart &self, ModelPart::GeometriesMapType &geometries) + { KRATOS_ERROR << "Setting geometries is not allowed! Trying to set value of ModelPart::Geometries."; }) .def("CreateSubModelPart", &ModelPart::CreateSubModelPart, py::return_value_policy::reference_internal) .def("NumberOfSubModelParts", &ModelPart::NumberOfSubModelParts) - .def("GetSubModelPart", py::overload_cast(&ModelPart::GetSubModelPart), py::return_value_policy::reference_internal) // non-const version + .def("GetSubModelPart", py::overload_cast(&ModelPart::GetSubModelPart), py::return_value_policy::reference_internal) // non-const version .def("RemoveSubModelPart", RemoveSubModelPart1) .def("RemoveSubModelPart", RemoveSubModelPart2) .def("HasSubModelPart", &ModelPart::HasSubModelPart) @@ -878,17 +890,17 @@ void AddModelPartToPython(pybind11::module& m) .def("AddNodalSolutionStepVariable", AddNodalSolutionStepVariable) .def("AddNodalSolutionStepVariable", AddNodalSolutionStepVariable) .def("AddNodalSolutionStepVariable", AddNodalSolutionStepVariable) - .def("AddNodalSolutionStepVariable", AddNodalSolutionStepVariable >) + .def("AddNodalSolutionStepVariable", AddNodalSolutionStepVariable>) .def("AddNodalSolutionStepVariable", AddNodalSolutionStepVariable) .def("AddNodalSolutionStepVariable", AddNodalSolutionStepVariable) - .def("AddNodalSolutionStepVariable", AddNodalSolutionStepVariable >) + .def("AddNodalSolutionStepVariable", AddNodalSolutionStepVariable>) .def("HasNodalSolutionStepVariable", HasNodalSolutionStepVariable) .def("HasNodalSolutionStepVariable", HasNodalSolutionStepVariable) .def("HasNodalSolutionStepVariable", HasNodalSolutionStepVariable) - .def("HasNodalSolutionStepVariable", HasNodalSolutionStepVariable >) + .def("HasNodalSolutionStepVariable", HasNodalSolutionStepVariable>) .def("HasNodalSolutionStepVariable", HasNodalSolutionStepVariable) .def("HasNodalSolutionStepVariable", HasNodalSolutionStepVariable) - .def("HasNodalSolutionStepVariable", HasNodalSolutionStepVariable >) + .def("HasNodalSolutionStepVariable", HasNodalSolutionStepVariable>) .def("GetNodalSolutionStepDataSize", &ModelPart::GetNodalSolutionStepDataSize) .def("GetNodalSolutionStepTotalDataSize", &ModelPart::GetNodalSolutionStepTotalDataSize) .def("OverwriteSolutionStepData", &ModelPart::OverwriteSolutionStepData) @@ -900,45 +912,44 @@ void AddModelPartToPython(pybind11::module& m) .def("CreateNewGeometry", ModelPartCreateNewGeometry5) .def("CreateNewGeometry", ModelPartCreateNewGeometry6) .def("CreateNewElement", ModelPartCreateNewElement) - .def("CreateNewElement", [](ModelPart& rModelPart, const std::string& ElementName, ModelPart::IndexType Id, - ModelPart::GeometryType::Pointer pGeometry, ModelPart::PropertiesType::Pointer pProperties) - {return rModelPart.CreateNewElement(ElementName, Id, pGeometry, pProperties);}) + .def("CreateNewElement", [](ModelPart &rModelPart, const std::string &ElementName, ModelPart::IndexType Id, ModelPart::GeometryType::Pointer pGeometry, ModelPart::PropertiesType::Pointer pProperties) + { return rModelPart.CreateNewElement(ElementName, Id, pGeometry, pProperties); }) .def("CreateNewCondition", ModelPartCreateNewCondition) - .def("CreateNewCondition", [](ModelPart& rModelPart, const std::string& ConditionName, ModelPart::IndexType Id, - ModelPart::GeometryType::Pointer pGeometry, ModelPart::PropertiesType::Pointer pProperties) - {return rModelPart.CreateNewCondition(ConditionName, Id, pGeometry, pProperties);}) + .def("CreateNewCondition", [](ModelPart &rModelPart, const std::string &ConditionName, ModelPart::IndexType Id, ModelPart::GeometryType::Pointer pGeometry, ModelPart::PropertiesType::Pointer pProperties) + { return rModelPart.CreateNewCondition(ConditionName, Id, pGeometry, pProperties); }) .def("GetCommunicator", ModelPartGetCommunicator, py::return_value_policy::reference_internal) .def("Check", &ModelPart::Check) .def("IsSubModelPart", &ModelPart::IsSubModelPart) .def("IsDistributed", &ModelPart::IsDistributed) - .def("AddNodes",AddNodesByIds) - .def("AddConditions",AddConditionsByIds) - .def("AddElements",AddElementsByIds) - .def("GetParentModelPart", [](ModelPart& self) -> ModelPart& {return self.GetParentModelPart();}, py::return_value_policy::reference_internal) - .def("GetRootModelPart", [](ModelPart& self) -> ModelPart& {return self.GetRootModelPart();}, py::return_value_policy::reference_internal) - .def("GetModel", [](ModelPart& self) -> Model& {return self.GetModel();}, py::return_value_policy::reference_internal) - .def_property("SubModelParts", [](ModelPart& self){ return self.SubModelParts(); }, - [](ModelPart& self, ModelPart::SubModelPartsContainerType& subs){ KRATOS_ERROR << "setting submodelparts is not allowed"; }) + .def("AddNodes", AddNodesByIds) + .def("AddConditions", AddConditionsByIds) + .def("AddElements", AddElementsByIds) + .def("AddGeometries", [](ModelPart &rModelPart, std::vector& rGeometriesIds) {rModelPart.AddGeometries(rGeometriesIds);}) + .def("GetParentModelPart", [](ModelPart &self) -> ModelPart & + { return self.GetParentModelPart(); }, py::return_value_policy::reference_internal) + .def("GetRootModelPart", [](ModelPart &self) -> ModelPart & + { return self.GetRootModelPart(); }, py::return_value_policy::reference_internal) + .def("GetModel", [](ModelPart &self) -> Model & + { return self.GetModel(); }, py::return_value_policy::reference_internal) + .def_property("SubModelParts", [](ModelPart &self) + { return self.SubModelParts(); }, [](ModelPart &self, ModelPart::SubModelPartsContainerType &subs) + { KRATOS_ERROR << "setting submodelparts is not allowed"; }) .def_property_readonly("MasterSlaveConstraints", ModelPartGetMasterSlaveConstraints1) - .def("GetHistoricalVariablesNames", [](ModelPart& rModelPart) -> std::unordered_set { + .def("GetHistoricalVariablesNames", [](ModelPart &rModelPart) -> std::unordered_set + { std::unordered_set variable_names; for(auto & variable: rModelPart.GetNodalSolutionStepVariablesList()) { variable_names.insert(variable.Name()); } - return variable_names; - }) - .def("GetNonHistoricalVariablesNames", [](ModelPart& rModelPart, ModelPart::NodesContainerType& rContainer, bool doFullSearch=false) -> std::unordered_set { - return GetNonHistoricalVariablesNames(rModelPart, rContainer, doFullSearch); - }) - .def("GetNonHistoricalVariablesNames", [](ModelPart& rModelPart, ModelPart::ElementsContainerType& rContainer, bool doFullSearch=false) -> std::unordered_set { - return GetNonHistoricalVariablesNames(rModelPart, rContainer, doFullSearch); - }) - .def("GetNonHistoricalVariablesNames", [](ModelPart& rModelPart, ModelPart::ConditionsContainerType& rContainer, bool doFullSearch=false) -> std::unordered_set { - return GetNonHistoricalVariablesNames(rModelPart, rContainer, doFullSearch); - }) - .def("HasMasterSlaveConstraint", [](ModelPart& rModelPart, ModelPart::IndexType MasterSlaveConstraintId) -> bool { - return rModelPart.HasMasterSlaveConstraint(MasterSlaveConstraintId); - }) + return variable_names; }) + .def("GetNonHistoricalVariablesNames", [](ModelPart &rModelPart, ModelPart::NodesContainerType &rContainer, bool doFullSearch = false) -> std::unordered_set + { return GetNonHistoricalVariablesNames(rModelPart, rContainer, doFullSearch); }) + .def("GetNonHistoricalVariablesNames", [](ModelPart &rModelPart, ModelPart::ElementsContainerType &rContainer, bool doFullSearch = false) -> std::unordered_set + { return GetNonHistoricalVariablesNames(rModelPart, rContainer, doFullSearch); }) + .def("GetNonHistoricalVariablesNames", [](ModelPart &rModelPart, ModelPart::ConditionsContainerType &rContainer, bool doFullSearch = false) -> std::unordered_set + { return GetNonHistoricalVariablesNames(rModelPart, rContainer, doFullSearch); }) + .def("HasMasterSlaveConstraint", [](ModelPart &rModelPart, ModelPart::IndexType MasterSlaveConstraintId) -> bool + { return rModelPart.HasMasterSlaveConstraint(MasterSlaveConstraintId); }) .def("GetMasterSlaveConstraint", ModelPartGetMasterSlaveConstraint1) .def("GetMasterSlaveConstraints", ModelPartGetMasterSlaveConstraints1) .def("RemoveMasterSlaveConstraint", ModelPartRemoveMasterSlaveConstraint1) @@ -949,8 +960,8 @@ void AddModelPartToPython(pybind11::module& m) .def("RemoveMasterSlaveConstraintsFromAllLevels", &ModelPart::RemoveMasterSlaveConstraintsFromAllLevels) .def("AddMasterSlaveConstraint", ModelPartAddMasterSlaveConstraint) .def("AddMasterSlaveConstraints", AddMasterSlaveConstraintsByIds) - .def("CreateNewMasterSlaveConstraint",CreateNewMasterSlaveConstraint1, py::return_value_policy::reference_internal) - .def("CreateNewMasterSlaveConstraint",CreateNewMasterSlaveConstraint2, py::return_value_policy::reference_internal) + .def("CreateNewMasterSlaveConstraint", CreateNewMasterSlaveConstraint1, py::return_value_policy::reference_internal) + .def("CreateNewMasterSlaveConstraint", CreateNewMasterSlaveConstraint2, py::return_value_policy::reference_internal) .def("__str__", PrintObject) ; }