From 0bc4a2bb789081f41c49e1bdfafeaad66831a95c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Mon, 18 Oct 2021 13:11:40 +0200 Subject: [PATCH 01/43] add ray locator utility --- .../find_intersected_objects_utility.cpp | 156 ++++++++++++ .../find_intersected_objects_utility.h | 229 ++++++++++++++++++ 2 files changed, 385 insertions(+) create mode 100644 applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.cpp create mode 100644 applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.h diff --git a/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.cpp b/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.cpp new file mode 100644 index 000000000000..e16bb1d0f2bc --- /dev/null +++ b/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.cpp @@ -0,0 +1,156 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Miguel Maso Sotomayor +// + + +// System includes + + +// External includes + + +// Project includes +#include "includes/model_part.h" +#include "utilities/parallel_utilities.h" +#include "find_intersected_objects_utility.h" + +namespace Kratos +{ + +/// Local Flags +KRATOS_CREATE_LOCAL_FLAG(FindIntersectedObjectsUtility, INTERSECT_ELEMENTS, 0); +KRATOS_CREATE_LOCAL_FLAG(FindIntersectedObjectsUtility, INTERSECT_CONDITIONS, 1); + +FindIntersectedObjectsUtility::FindIntersectedObjectsUtility( + ModelPart& rThisModelPart, + Parameters ThisParameters +) : mrModelPart(rThisModelPart) +{ + Parameters default_parameters(R"({ + "intersect_elements" : true, + "intersect_conditions" : false + })"); + + // Set the flags + const bool intersect_elements = ThisParameters["intersect_elements"].GetBool(); + const bool intersect_conditions = ThisParameters["intersect_conditions"].GetBool(); + + mOptions.Set(INTERSECT_ELEMENTS, intersect_elements); + mOptions.Set(INTERSECT_CONDITIONS, intersect_conditions); + + if (mrModelPart.NumberOfNodes() > 0) { + GenerateOctree(); + } +} + +void FindIntersectedObjectsUtility::UpdateSearchStructure() +{ + if (mrModelPart.NumberOfNodes() > 0) { + GenerateOctree(); + } +} + +void FindIntersectedObjectsUtility::FindIntersectedObjects( + const GeometryType& rGeometry, + PointerVector& rResults) +{ + OctreeCellVectorType leaves; + leaves.clear(); + auto p_geometrical_object = Kratos::make_intrusive(0, Kratos::make_shared(rGeometry)); + mpOctree->GetIntersectedLeaves(p_geometrical_object, leaves); + FindIntersectedObjects(rGeometry, leaves, rResults); +} + +void FindIntersectedObjectsUtility::FindIntersectedObjects( + const GeometryType& rIntersectionGeometry, + OctreeCellVectorType& rLeaves, + PointerVector& rResults) +{ + for (auto p_leaf : rLeaves) { + for (auto p_volume_object : *(p_leaf->pGetObjects())) { + if (p_volume_object->GetGeometry().HasIntersection(rIntersectionGeometry)) { + rResults.push_back(p_volume_object); + } + } + } +} + +void FindIntersectedObjectsUtility::GenerateOctree() +{ + SetOctreeBoundingBox(); + + // Adding the volume model part to the octree + for (auto it_node = mrModelPart.NodesBegin(); it_node != mrModelPart.NodesEnd(); it_node++) { +#ifdef KRATOS_USE_AMATRIX // This macro definition is for the migration period and to be removed afterward please do not use it + mpOctree->Insert(it_node->Coordinates().data()); +#else + mpOctree->Insert(it_node->Coordinates().data().data()); +#endif // ifdef KRATOS_USE_AMATRIX + } + + // Add elements + if (mOptions.Is(INTERSECT_ELEMENTS)) { + auto& r_elements_array = mrModelPart.Elements(); + const auto it_elements_begin = r_elements_array.begin(); + const auto number_of_elements = r_elements_array.size(); + + // Iterate over the elements + for (int i = 0; i < static_cast(number_of_elements); i++) { + auto it_elements = it_elements_begin + i; + mpOctree->Insert(*(it_elements).base()); + } + } + + // Add conditions + if (mOptions.Is(INTERSECT_CONDITIONS)) { + auto& r_conditions_array = mrModelPart.Conditions(); + const auto it_conditions_begin = r_conditions_array.begin(); + const auto number_of_conditions = r_conditions_array.size(); + + // Iterate over the conditions + for (int i = 0; i < static_cast(number_of_conditions); i++) { + auto it_conditions = it_conditions_begin + i; + mpOctree->Insert(*(it_conditions).base()); + } + } +} + +void FindIntersectedObjectsUtility::SetOctreeBoundingBox() +{ + PointType low(mrModelPart.NodesBegin()->Coordinates()); + PointType high(mrModelPart.NodesBegin()->Coordinates()); + + // Loop over all nodes in the volume model part + for (auto it_node = mrModelPart.NodesBegin(); it_node != mrModelPart.NodesEnd(); it_node++) { + const array_1d& r_coordinates = it_node->Coordinates(); + for (IndexType i = 0; i < 3; i++) { + low[i] = r_coordinates[i] < low[i] ? r_coordinates[i] : low[i]; + high[i] = r_coordinates[i] > high[i] ? r_coordinates[i] : high[i]; + } + } + + // Slightly increase the bounding box size to avoid problems with geometries in the borders + // Note that std::numeric_limits::double() is added for the 2D cases. Otherwise, the + // third component will be 0, breaking the octree behaviour. + for(IndexType i = 0 ; i < 3; i++) { + low[i] -= std::abs(high[i] - low[i])*1e-3 + std::numeric_limits::epsilon(); + high[i] += std::abs(high[i] - low[i])*1e-3 + std::numeric_limits::epsilon(); + } + + // TODO: Octree needs refactoring to work with BoundingBox. Pooyan. +#ifdef KRATOS_USE_AMATRIX // This macro definition is for the migration period and to be removed afterward please do not use it + mpOctree->SetBoundingBox(low.data(), high.data()); +#else + mpOctree->SetBoundingBox(low.data().data(), high.data().data()); +#endif // ifdef KRATOS_USE_AMATRIX +} + +} // namespace Kratos. diff --git a/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.h b/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.h new file mode 100644 index 000000000000..49ecf1930375 --- /dev/null +++ b/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.h @@ -0,0 +1,229 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Miguel Maso Sotomayor +// + + +#ifndef KRATOS_FIND_INTERSECTED_ELEMENTS_UTILITY_H_INCLUDED +#define KRATOS_FIND_INTERSECTED_ELEMENTS_UTILITY_H_INCLUDED + + +// System includes + + +// External includes + + +// Project includes +#include "includes/kratos_parameters.h" +#include "spatial_containers/octree_binary.h" +#include "processes/find_intersected_geometrical_objects_process.h" + +namespace Kratos +{ +///@name Kratos Globals +///@{ + +///@} +///@name Type Definitions +///@{ + +///@} +///@name Enum's +///@{ + +///@} +///@name Functions +///@{ + +///@} +///@name Kratos Classes +///@{ + +/** + * @brief Forward declaration of ModelPart + */ +class ModelPart; + +/** + * @brief Find the objects in a volume model part that are intersected ty the given geometry + * @author Miguel Maso Sotomayor + * @ingroup ShallowWaterApplication +*/ +class KRATOS_API(SHALLOW_WATER_APPLICATION) FindIntersectedObjectsUtility +{ +public: + ///@name Type Definitions + ///@{ + + /// Pointer definition of FindIntersectedObjectsUtility + KRATOS_CLASS_POINTER_DEFINITION(FindIntersectedObjectsUtility); + + /// Local Flags + KRATOS_DEFINE_LOCAL_FLAG(INTERSECT_ELEMENTS); + KRATOS_DEFINE_LOCAL_FLAG(INTERSECT_CONDITIONS); + + /// Definition of the point type + using PointType = Point; + + /// Definition of the node type + using NodeType = Node<3>; + + /// Definition of the geometry type + using GeometryType = Geometry; + + /// Octree definitions + using ConfigurationType = Internals::DistanceSpatialContainersConfigure; + using CellType = OctreeBinaryCell; + using OctreeType = OctreeBinary; + using OctreePointerType = unique_ptr; + using CellNodeDataType = typename ConfigurationType::cell_node_data_type; + using OctreeCellVectorType = std::vector; + + ///@} + ///@name Life Cycle + ///@{ + + /// Constructor. + FindIntersectedObjectsUtility(ModelPart& rThisModelPart, Parameters ThisParameters = Parameters()); + + /// Destructor. + ~FindIntersectedObjectsUtility(){} + + ///@} + ///@name Operators + ///@{ + + + ///@} + ///@name Operations + ///@{ + + void UpdateSearchStructure(); + + void FindIntersectedObjects(const GeometryType& rGeometry, PointerVector& rResults); + + ///@} + ///@name Access + ///@{ + + + ///@} + ///@name Inquiry + ///@{ + + + ///@} + ///@name Input and output + ///@{ + + /// Turn back information as a string. + std::string Info() const { + return "FindIntersectedObjectsUtility"; + } + + /// Print information about this object. + void PrintInfo(std::ostream& rOStream) const { + rOStream << Info(); + } + + /// Print object's data. + void PrintData(std::ostream& rOStream) const {} + + ///@} + ///@name Friends + ///@{ + + + ///@} + +private: + ///@name Static Member Variables + ///@{ + + + ///@} + ///@name Member Variables + ///@{ + + Flags mOptions; + ModelPart& mrModelPart; + OctreePointerType mpOctree; + + ///@} + ///@name Private Operators + ///@{ + + + ///@} + ///@name Private Operations + ///@{ + + void GenerateOctree(); + + void SetOctreeBoundingBox(); + + void FindIntersectedObjects( + const GeometryType& rIntersectionGeometry, + OctreeCellVectorType& rLeaves, + PointerVector& rResults); + + ///@} + ///@name Private Access + ///@{ + + + ///@} + ///@name Private Inquiry + ///@{ + + + ///@} + ///@name Un accessible methods + ///@{ + + /// Assignment operator. + FindIntersectedObjectsUtility& operator=(FindIntersectedObjectsUtility const& rOther); + + /// Copy constructor. + FindIntersectedObjectsUtility(FindIntersectedObjectsUtility const& rOther); + + ///@} + +}; // Class FindIntersectedObjectsUtility + +///@} + +///@name Type Definitions +///@{ + + +///@} +///@name Input and output +///@{ + + +/// input stream function +inline std::istream& operator >> (std::istream& rIStream, FindIntersectedObjectsUtility& rThis); + +/// output stream function +inline std::ostream& operator << (std::ostream& rOStream, const FindIntersectedObjectsUtility& rThis) +{ + rThis.PrintInfo(rOStream); + rOStream << std::endl; + rThis.PrintData(rOStream); + + return rOStream; +} +///@} + +} // namespace Kratos. + +#endif // KRATOS_FIND_INTERSECTED_ELEMENTS_UTILITY_H_INCLUDED defined From 5638eff4f2a16c4c1c8732e52a0776e6ed40b618 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Mon, 18 Oct 2021 13:11:52 +0200 Subject: [PATCH 02/43] add c++ process --- .../depth_integration_process.cpp | 123 ++++++++++ .../depth_integration_process.h | 224 ++++++++++++++++++ 2 files changed, 347 insertions(+) create mode 100644 applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp create mode 100644 applications/ShallowWaterApplication/custom_processes/depth_integration_process.h diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp new file mode 100644 index 000000000000..b700a83f63fe --- /dev/null +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -0,0 +1,123 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Miguel Maso Sotomayor +// + +// System includes + + +// External includes + + +// Project includes +#include "includes/model_part.h" +#include "geometries/line_3d_2.h" +#include "depth_integration_process.h" +#include "utilities/parallel_utilities.h" +#include "shallow_water_application_variables.h" +#include "custom_utilities/find_intersected_objects_utility.h" + +namespace Kratos +{ + +DepthIntegrationProcess::DepthIntegrationProcess( + Model& rModel, + Parameters ThisParameters + ) : Process(), + mrVolumeModelPart(rModel.GetModelPart(ThisParameters["volume_model_part_name"].GetString())), + mrInterfaceModelPart(rModel.GetModelPart(ThisParameters["interface_model_part_name"].GetString())) +{ + ThisParameters.ValidateAndAssignDefaults(this->GetDefaultParameters()); + mDirection = ThisParameters["direction_of_integration"].GetVector(); + mDirection /= norm_2(mDirection); +} + +void DepthIntegrationProcess::Execute() +{ + double bottom, top; + GetVolumePartBounds(bottom, top); + FindIntersectedObjectsUtility intersections(mrVolumeModelPart); + for (auto& node : mrInterfaceModelPart.Nodes()) { + auto integration_line = CreateIntegrationLine(node, bottom, top); + PointerVector intersected_objects; + intersections.FindIntersectedObjects(*integration_line, intersected_objects); + Integrate(intersected_objects, node); + } +} + +void DepthIntegrationProcess::Integrate(PointerVector& rObjects, NodeType& rNode) +{ + array_1d velocity = ZeroVector(3); + double min_elevation = 1e6; + double max_elevation = -1e6; + int num_nodes = 0; + for (auto& object : rObjects) { + array_1d obj_velocity = ZeroVector(3); + double obj_min_elevation = 1e6; + double obj_max_elevation = -1e6; + int obj_num_nodes = object.GetGeometry().size(); + for (auto& node : object.GetGeometry()) { + velocity += node.FastGetSolutionStepValue(VELOCITY); + obj_min_elevation = std::min(obj_min_elevation, inner_prod(mDirection, node)); + obj_max_elevation = std::max(obj_max_elevation, inner_prod(mDirection, node)); + } + velocity += obj_velocity; + min_elevation = std::min(min_elevation, obj_min_elevation); + max_elevation = std::min(max_elevation, obj_max_elevation); + num_nodes += obj_num_nodes; + } + velocity /= num_nodes; + double height = max_elevation - min_elevation; + rNode.FastGetSolutionStepValue(VELOCITY) = velocity; + rNode.FastGetSolutionStepValue(HEIGHT) = height; +} + +void DepthIntegrationProcess::GetVolumePartBounds(double& rMin, double& rMax) +{ + using MultipleReduction = CombinedReduction,MaxReduction>; + + std::tie(rMin, rMax) = block_for_each(mrInterfaceModelPart.Nodes(), [&](NodeType& node){ + const double distance = inner_prod(mDirection, node); + return std::make_tuple(distance, distance); + }); +} + +Geometry>::Pointer DepthIntegrationProcess::CreateIntegrationLine( + const NodeType& rNode, + const double Bottom, + const double Top) +{ + array_1d origin = rNode - mDirection * inner_prod(rNode, mDirection); + array_1d start = origin + Top * mDirection; + array_1d end = origin + Bottom * mDirection; + PointerVector> nodes; + nodes.push_back(NodeType::Pointer(new NodeType(0, start))); + nodes.push_back(NodeType::Pointer(new NodeType(0, end))); + return Kratos::make_shared>(nodes); +} + +int DepthIntegrationProcess::Check() +{ + KRATOS_ERROR_IF(mDirection.size() != 3) << "DepthIntegrationProcess: The direction of integration must have three coordinates." << std::endl; + return 0; +} + +const Parameters DepthIntegrationProcess::GetDefaultParameters() const +{ + auto default_parameters = Parameters(R"( + { + "volume_model_part_name" : "", + "interface_model_part_name" : "", + "direction_of_integration" : [0.0, 0.0, 1.0] + })"); + return default_parameters; +} + +} // namespace Kratos. diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h new file mode 100644 index 000000000000..94ff26e24e9c --- /dev/null +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h @@ -0,0 +1,224 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Miguel Maso Sotomayor +// + +#ifndef KRATOS_DEPTH_INTEGRATION_PROCESS_H_INCLUDED +#define KRATOS_DEPTH_INTEGRATION_PROCESS_H_INCLUDED + + +// System includes + + +// External includes + + +// Project includes +#include "containers/model.h" +#include "processes/process.h" +#include "includes/kratos_parameters.h" + +namespace Kratos +{ +///@name Kratos Globals +///@{ + +///@} +///@name Type Definitions +///@{ + +///@} +///@name Enum's +///@{ + +///@} +///@name Functions +///@{ + +///@} +///@name Kratos Classes +///@{ + +/** + * @brief Forward declaration of ModelPart + */ +class ModelPart; + +/** + * @class DepthIntegrationProcess + * @ingroup ShallowWaterApplication + * @brief Calculate the minimum distance from all the nodes to a boundary condition in 2D + * @details The boundary conditions are assumed to be contained in a line + * @author Miguel Maso Sotomayor + */ +class KRATOS_API(SHALLOW_WATER_APPLICATION) DepthIntegrationProcess : public Process +{ +public: + ///@name Type Definitions + ///@{ + + /// Pointer definition of DepthIntegrationProcess + KRATOS_CLASS_POINTER_DEFINITION(DepthIntegrationProcess); + + /// Definition of the node type + using NodeType = Node<3>; + + /// Definition of the geometry type + using GeometryType = Geometry; + + ///@} + ///@name Life Cycle + ///@{ + + /** + * @brief Default constructor. + * @details Removed + */ + DepthIntegrationProcess() = delete; + + /** + * Constructor with Model and Parameters + */ + DepthIntegrationProcess(Model& rModel, Parameters ThisParameters = Parameters()); + + /// Destructor. + ~DepthIntegrationProcess() override = default; + + ///@} + ///@name Operators + ///@{ + + + ///@} + ///@name Operations + ///@{ + + void Execute() override; + + int Check() override; + + const Parameters GetDefaultParameters() const override; + + ///@} + ///@name Access + ///@{ + + + ///@} + ///@name Inquiry + ///@{ + + + ///@} + ///@name Input and output + ///@{ + + /// Turn back information as a string. + std::string Info() const override { + std::stringstream buffer; + buffer << "DepthIntegrationProcess"; + return buffer.str(); + } + + /// Print information about this object. + void PrintInfo(std::ostream& rOStream) const override { + rOStream << Info(); + } + + /// Print object's data. + void PrintData(std::ostream& rOStream) const override {} + + ///@} + ///@name Friends + ///@{ + + + ///@} + +private: + ///@name Static Member Variables + ///@{ + + ModelPart& mrVolumeModelPart; + ModelPart& mrInterfaceModelPart; + array_1d mDirection; + + ///@} + ///@name Member Variables + ///@{ + + + ///@} + ///@name Private Operators + ///@{ + + + ///@} + ///@name Private Operations + ///@{ + + void Integrate(PointerVector& rObjects, NodeType& rNode); + + void GetVolumePartBounds(double& rMin, double& rMax); + + GeometryType::Pointer CreateIntegrationLine(const NodeType& rNode, const double Bottom, const double Top); + + ///@} + ///@name Private Access + ///@{ + + + ///@} + ///@name Private Inquiry + ///@{ + + + ///@} + ///@name Un accessible methods + ///@{ + + /// Assignment operator. + DepthIntegrationProcess& operator=(DepthIntegrationProcess const& rOther) = delete; + + /// Copy constructor. + DepthIntegrationProcess(DepthIntegrationProcess const& rOther) = delete; + + ///@} + +}; // Class DepthIntegrationProcess + +///@} + +///@name Type Definitions +///@{ + + +///@} +///@name Input and output +///@{ + +/// input stream function +inline std::istream& operator >> (std::istream& rIStream, DepthIntegrationProcess& rThis); + +/// output stream function +inline std::ostream& operator << (std::ostream& rOStream, const DepthIntegrationProcess& rThis) +{ + rThis.PrintInfo(rOStream); + rOStream << std::endl; + rThis.PrintData(rOStream); + + return rOStream; +} + +///@} + +} // namespace Kratos. + +#endif // KRATOS_DEPTH_INTEGRATION_PROCESS_H_INCLUDED defined From f81698bedce529b32f7543fe5006b28f953d99b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Fri, 22 Oct 2021 16:47:50 +0200 Subject: [PATCH 03/43] export to tpython --- .../custom_python/add_custom_processes_to_python.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/applications/ShallowWaterApplication/custom_python/add_custom_processes_to_python.cpp b/applications/ShallowWaterApplication/custom_python/add_custom_processes_to_python.cpp index 271c9860b575..2b9ff76326ac 100644 --- a/applications/ShallowWaterApplication/custom_python/add_custom_processes_to_python.cpp +++ b/applications/ShallowWaterApplication/custom_python/add_custom_processes_to_python.cpp @@ -24,6 +24,7 @@ #include "custom_processes/apply_perturbation_function_process.h" #include "custom_processes/apply_sinusoidal_function_process.h" #include "custom_processes/calculate_distance_to_boundary_process.h" +#include "custom_processes/depth_integration_process.h" namespace Kratos @@ -68,6 +69,11 @@ namespace Python .def(py::init()) ; + py::class_ + (m, "DepthIntegrationProcess") + .def(py::init()) + ; + } } // namespace Python. From 7b6aec8db0257ba7fd8c1771f82c906979e1b28b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Fri, 22 Oct 2021 16:48:28 +0200 Subject: [PATCH 04/43] fix pointer --- .../custom_processes/depth_integration_process.cpp | 8 ++++---- .../find_intersected_objects_utility.cpp | 10 ++++++---- .../find_intersected_objects_utility.h | 2 +- 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index b700a83f63fe..d7dcd4829361 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -45,9 +45,9 @@ void DepthIntegrationProcess::Execute() GetVolumePartBounds(bottom, top); FindIntersectedObjectsUtility intersections(mrVolumeModelPart); for (auto& node : mrInterfaceModelPart.Nodes()) { - auto integration_line = CreateIntegrationLine(node, bottom, top); + GeometryType::Pointer integration_line = CreateIntegrationLine(node, bottom, top); PointerVector intersected_objects; - intersections.FindIntersectedObjects(*integration_line, intersected_objects); + intersections.FindIntersectedObjects(integration_line, intersected_objects); Integrate(intersected_objects, node); } } @@ -70,7 +70,7 @@ void DepthIntegrationProcess::Integrate(PointerVector& rObjec } velocity += obj_velocity; min_elevation = std::min(min_elevation, obj_min_elevation); - max_elevation = std::min(max_elevation, obj_max_elevation); + max_elevation = std::max(max_elevation, obj_max_elevation); num_nodes += obj_num_nodes; } velocity /= num_nodes; @@ -83,7 +83,7 @@ void DepthIntegrationProcess::GetVolumePartBounds(double& rMin, double& rMax) { using MultipleReduction = CombinedReduction,MaxReduction>; - std::tie(rMin, rMax) = block_for_each(mrInterfaceModelPart.Nodes(), [&](NodeType& node){ + std::tie(rMin, rMax) = block_for_each(mrVolumeModelPart.Nodes(), [&](NodeType& node){ const double distance = inner_prod(mDirection, node); return std::make_tuple(distance, distance); }); diff --git a/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.cpp b/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.cpp index e16bb1d0f2bc..97b7e56da106 100644 --- a/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.cpp +++ b/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.cpp @@ -32,12 +32,14 @@ KRATOS_CREATE_LOCAL_FLAG(FindIntersectedObjectsUtility, INTERSECT_CONDITIONS, 1) FindIntersectedObjectsUtility::FindIntersectedObjectsUtility( ModelPart& rThisModelPart, Parameters ThisParameters -) : mrModelPart(rThisModelPart) +) : mrModelPart(rThisModelPart), + mpOctree(new OctreeType()) { Parameters default_parameters(R"({ "intersect_elements" : true, "intersect_conditions" : false })"); + ThisParameters.ValidateAndAssignDefaults(default_parameters); // Set the flags const bool intersect_elements = ThisParameters["intersect_elements"].GetBool(); @@ -59,14 +61,14 @@ void FindIntersectedObjectsUtility::UpdateSearchStructure() } void FindIntersectedObjectsUtility::FindIntersectedObjects( - const GeometryType& rGeometry, + GeometryType::Pointer pGeometry, PointerVector& rResults) { OctreeCellVectorType leaves; leaves.clear(); - auto p_geometrical_object = Kratos::make_intrusive(0, Kratos::make_shared(rGeometry)); + auto p_geometrical_object = Kratos::make_intrusive(0, pGeometry); mpOctree->GetIntersectedLeaves(p_geometrical_object, leaves); - FindIntersectedObjects(rGeometry, leaves, rResults); + FindIntersectedObjects(*pGeometry, leaves, rResults); } void FindIntersectedObjectsUtility::FindIntersectedObjects( diff --git a/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.h b/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.h index 49ecf1930375..e6382e4e6aab 100644 --- a/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.h +++ b/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.h @@ -108,7 +108,7 @@ class KRATOS_API(SHALLOW_WATER_APPLICATION) FindIntersectedObjectsUtility void UpdateSearchStructure(); - void FindIntersectedObjects(const GeometryType& rGeometry, PointerVector& rResults); + void FindIntersectedObjects(GeometryType::Pointer rGeometry, PointerVector& rResults); ///@} ///@name Access From 3a347ea88158b3509b944fc1dc2202716e3cb04c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Fri, 22 Oct 2021 17:04:29 +0200 Subject: [PATCH 05/43] Add new intersections type in tetrahedra --- kratos/geometries/tetrahedra_3d_4.h | 50 ++++++++++++++++++++--------- 1 file changed, 34 insertions(+), 16 deletions(-) diff --git a/kratos/geometries/tetrahedra_3d_4.h b/kratos/geometries/tetrahedra_3d_4.h index 09d51d3ba18f..fdf441f5b314 100644 --- a/kratos/geometries/tetrahedra_3d_4.h +++ b/kratos/geometries/tetrahedra_3d_4.h @@ -1447,29 +1447,47 @@ template class Tetrahedra3D4 : public Geometry } - /// detect if two tetrahedra are intersected - bool HasIntersection( const BaseType& rThisGeometry) override + /// Detect if the tetrahedra is intersected with another geometry + bool HasIntersection(const BaseType& rThisGeometry) override { + if (rThisGeometry.LocalSpaceDimension() < this->LocalSpaceDimension()) { + // Check the intersection of each face against the intersecting object + const auto faces = this->GenerateFaces(); + for (auto& face : faces) { + if (face.HasIntersection(rThisGeometry)){ + return true; + } + } - array_1d plane; - std::vector Intersection; + // Let check second geometry is inside. + // Considering that there are no intersection, if one point is inside all of it is inside. + array_1d local_point; + if (this->IsInside(rThisGeometry.GetPoint(0), local_point)){ + return true; + } - //const BaseType& geom_1 = *this; - const BaseType& geom_2 = rThisGeometry; + return false; + } + else { + array_1d plane; + std::vector Intersection; - GetPlanes(plane); - Intersection.push_back(geom_2); - for (unsigned int i = 0; i < 4; ++i) - { - std::vector inside; - for (unsigned int j = 0; j < Intersection.size(); ++j) + const BaseType& geom_2 = rThisGeometry; + + GetPlanes(plane); + Intersection.push_back(geom_2); + for (unsigned int i = 0; i < 4; ++i) { - SplitAndDecompose(Intersection[j], plane[i], inside); + std::vector inside; + for (unsigned int j = 0; j < Intersection.size(); ++j) + { + SplitAndDecompose(Intersection[j], plane[i], inside); + } + Intersection = inside; } - Intersection = inside; - } - return bool (Intersection.size() > 0); + return bool (Intersection.size() > 0); + } } From c8af23594174eed0fa71a985eee5741473c25243 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Tue, 26 Oct 2021 14:37:17 +0200 Subject: [PATCH 06/43] DRAFT migrating to find_intersected_process --- .../depth_integration_process.cpp | 72 ++++++++++++++++++- .../depth_integration_process.h | 9 ++- 2 files changed, 78 insertions(+), 3 deletions(-) diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index d7dcd4829361..0a603b4ffa45 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -20,6 +20,7 @@ #include "includes/model_part.h" #include "geometries/line_3d_2.h" #include "depth_integration_process.h" +#include "utilities/variable_utils.h" #include "utilities/parallel_utilities.h" #include "shallow_water_application_variables.h" #include "custom_utilities/find_intersected_objects_utility.h" @@ -37,12 +38,28 @@ DepthIntegrationProcess::DepthIntegrationProcess( ThisParameters.ValidateAndAssignDefaults(this->GetDefaultParameters()); mDirection = ThisParameters["direction_of_integration"].GetVector(); mDirection /= norm_2(mDirection); + + if (rModel.HasModelPart("integration_auxiliary_model_part")) { + mpIntegrationModelPart = &rModel.GetModelPart("integration_auxiliary_model_part"); + } else { + mpIntegrationModelPart = &rModel.CreateModelPart("integration_auxiliary_model_part"); + } } void DepthIntegrationProcess::Execute() { double bottom, top; - GetVolumePartBounds(bottom, top); + GetBoundingVolumeLimits(bottom, top); + // InitializeIntegrationModelPart(); + // FindIntersectedGeometricalObjectsProcess find_intersected_objects_process(*mpIntegrationModelPart, mrVolumeModelPart); + // find_intersected_objects_process.ExecuteInitialize(); + // for (auto& node : mrInterfaceModelPart.Nodes()) { + // InitializeIntegrationLine(); + // SetIntegrationLine(node, bottom, top); + // find_intersected_objects_process.FindIntersections(); + // auto intersected_objects = find_intersected_objects_process.GetIntersections(); + // Integrate(intersected_objects[0], node); + // } FindIntersectedObjectsUtility intersections(mrVolumeModelPart); for (auto& node : mrInterfaceModelPart.Nodes()) { GeometryType::Pointer integration_line = CreateIntegrationLine(node, bottom, top); @@ -79,7 +96,7 @@ void DepthIntegrationProcess::Integrate(PointerVector& rObjec rNode.FastGetSolutionStepValue(HEIGHT) = height; } -void DepthIntegrationProcess::GetVolumePartBounds(double& rMin, double& rMax) +void DepthIntegrationProcess::GetBoundingVolumeLimits(double& rMin, double& rMax) { using MultipleReduction = CombinedReduction,MaxReduction>; @@ -89,6 +106,57 @@ void DepthIntegrationProcess::GetVolumePartBounds(double& rMin, double& rMax) }); } +void DepthIntegrationProcess::InitializeIntegrationModelPart() +{ + // Empty the model part + VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Nodes()); + VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Elements()); + VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Conditions()); + mpIntegrationModelPart->RemoveNodesFromAllLevels(); + mpIntegrationModelPart->RemoveElementsFromAllLevels(); + mpIntegrationModelPart->RemoveConditionsFromAllLevels(); + + // Add a dummy node in order to allow the generation of the octree + if (mrVolumeModelPart.NumberOfNodes() > 0) { + auto i_node = mrVolumeModelPart.NodesBegin(); + mpIntegrationModelPart->AddNode(&*i_node); + } else { + KRATOS_ERROR << "DepthIntegrationProcess: The volume model part is empty." << std::endl; + } +} + +void DepthIntegrationProcess::InitializeIntegrationLine() +{ + // Remove the dummy node + VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Nodes()); + mpIntegrationModelPart->RemoveNodesFromAllLevels(); + + // Create the integration lines + ModelPart::PropertiesType::Pointer p_prop; + if (mpIntegrationModelPart->HasProperties(1)) { + p_prop = mpIntegrationModelPart->pGetProperties(1); + } else { + p_prop = mpIntegrationModelPart->CreateNewProperties(1); + } + mpIntegrationModelPart->CreateNewNode(1, 0.0, 0.0, 0.0); + mpIntegrationModelPart->CreateNewNode(2, 0.0, 0.0, 1.0); + mpIntegrationModelPart->CreateNewElement("Element3D2N", 1, {{1, 2}}, p_prop); +} + +void DepthIntegrationProcess::SetIntegrationLine( + const NodeType& rNode, + const double Bottom, + const double Top) +{ + // Set the integration limits to the lines + array_1d base_point = rNode - mDirection * inner_prod(rNode, mDirection); + array_1d start = base_point + Top * mDirection; + array_1d end = base_point + Bottom * mDirection; + auto& integration_line = mpIntegrationModelPart->GetElement(1); + integration_line.GetGeometry()[0].Coordinates() = start; + integration_line.GetGeometry()[1].Coordinates() = end; +} + Geometry>::Pointer DepthIntegrationProcess::CreateIntegrationLine( const NodeType& rNode, const double Bottom, diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h index 94ff26e24e9c..05a1fe2ff703 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h @@ -148,6 +148,7 @@ class KRATOS_API(SHALLOW_WATER_APPLICATION) DepthIntegrationProcess : public Pro ModelPart& mrVolumeModelPart; ModelPart& mrInterfaceModelPart; + ModelPart* mpIntegrationModelPart; array_1d mDirection; ///@} @@ -166,7 +167,13 @@ class KRATOS_API(SHALLOW_WATER_APPLICATION) DepthIntegrationProcess : public Pro void Integrate(PointerVector& rObjects, NodeType& rNode); - void GetVolumePartBounds(double& rMin, double& rMax); + void GetBoundingVolumeLimits(double& rMin, double& rMax); + + void InitializeIntegrationModelPart(); + + void InitializeIntegrationLine(); + + void SetIntegrationLine(const NodeType& rNode, const double Bottom, const double Top); GeometryType::Pointer CreateIntegrationLine(const NodeType& rNode, const double Bottom, const double Top); From 8e9c1fdaa883cff2935a45d121723988d056fe90 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Tue, 26 Oct 2021 16:41:48 +0200 Subject: [PATCH 07/43] allow nonhistorical --- .../custom_processes/depth_integration_process.cpp | 8 +++++--- .../custom_processes/depth_integration_process.h | 10 ++++++++++ 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index 0a603b4ffa45..85abcbcae521 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -36,6 +36,7 @@ DepthIntegrationProcess::DepthIntegrationProcess( mrInterfaceModelPart(rModel.GetModelPart(ThisParameters["interface_model_part_name"].GetString())) { ThisParameters.ValidateAndAssignDefaults(this->GetDefaultParameters()); + mStoreHistorical = ThisParameters["store_historical_database"].GetBool(); mDirection = ThisParameters["direction_of_integration"].GetVector(); mDirection /= norm_2(mDirection); @@ -92,8 +93,8 @@ void DepthIntegrationProcess::Integrate(PointerVector& rObjec } velocity /= num_nodes; double height = max_elevation - min_elevation; - rNode.FastGetSolutionStepValue(VELOCITY) = velocity; - rNode.FastGetSolutionStepValue(HEIGHT) = height; + SetValue(rNode, VELOCITY, velocity); + SetValue(rNode, HEIGHT, height); } void DepthIntegrationProcess::GetBoundingVolumeLimits(double& rMin, double& rMax) @@ -183,7 +184,8 @@ const Parameters DepthIntegrationProcess::GetDefaultParameters() const { "volume_model_part_name" : "", "interface_model_part_name" : "", - "direction_of_integration" : [0.0, 0.0, 1.0] + "direction_of_integration" : [0.0, 0.0, 1.0], + "store_historical_database" : false, })"); return default_parameters; } diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h index 05a1fe2ff703..4e8ba3a5753d 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h @@ -150,6 +150,7 @@ class KRATOS_API(SHALLOW_WATER_APPLICATION) DepthIntegrationProcess : public Pro ModelPart& mrInterfaceModelPart; ModelPart* mpIntegrationModelPart; array_1d mDirection; + bool mStoreHistorical; ///@} ///@name Member Variables @@ -177,6 +178,15 @@ class KRATOS_API(SHALLOW_WATER_APPLICATION) DepthIntegrationProcess : public Pro GeometryType::Pointer CreateIntegrationLine(const NodeType& rNode, const double Bottom, const double Top); + template> + void SetValue(NodeType& rNode, const TVarType& rVariable, TDataType& rValue) + { + if (mStoreHistorical) + rNode.FastGetSolutionStepValue(rVariable) = rValue; + else + rNode.GetValue(rVariable) = rValue; + } + ///@} ///@name Private Access ///@{ From 4a6de59a9a3d154470c254719e36f73666302e80 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Tue, 26 Oct 2021 16:42:04 +0200 Subject: [PATCH 08/43] first version python process --- .../depth_integration_output_process.py | 77 +++++++++++++++++++ 1 file changed, 77 insertions(+) create mode 100644 applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py new file mode 100644 index 000000000000..76bf25b3b833 --- /dev/null +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py @@ -0,0 +1,77 @@ +import KratosMultiphysics as KM +import KratosMultiphysics.ShallowWaterApplication as SW +from KratosMultiphysics.HDF5Application import single_mesh_temporal_output_process + +def Factory(settings, model): + if not isinstance(settings, KM.Parameters): + raise Exception("expected input shall be a Parameters object, encapsulating a json string") + return DepthIntegrationOutputProcess(model, settings["Parameters"]) + +class DepthIntegrationOutputProcess(KM.OutputProcess): + """DepthIntegrationOutputProcess + + This process performs a depth integration from a Navier-Stokes domain to a shallow water domain. + The depth integration values are stored in the nodes of the shallow water domain and + printed in HDF5 format. + """ + + @staticmethod + def GetDefaultParameters(): + return KM.Parameters("""{ + "volume_model_part_name" : "", + "interface_model_part_name" : "", + "store_historical_database" : false, + "direction_of_integration" : [0.0, 0.0, 1.0], + "interval" : [0.0,"End"], + "file_settings" : {}, + "output_time_settings" : {}, + }""") + + def __init__(self, model, settings): + """The constructor of the DepthIntegrationOutputProcess""" + + KM.OutputProcess.__init__(self) + + self.settings = settings + self.settings.ValidateAndAssignDefaults(self.GetDefaultParameters()) + + self.volume_model_part = model[self.settings["volume_model_part_name"].GetString()] + self.interface_model_part = model[self.settings["interface_model_part_name"].GetString()] + self.store_historical = settings["store_historical_database"].GetBool() + self.interval = KM.IntervalUtility(settings) + + integration_settings = KM.Parameters() + integration_settings.AddValue("volume_model_part_name", settings["volume_model_part_name"]) + integration_settings.AddValue("interface_model_part_name", settings["interface_model_part_name"]) + integration_settings.AddValue("direction_of_integration", settings["direction_of_integration"]) + integration_settings.AddValue("store_historical_database", settings["store_historical_database"]) + + self.integration_process = SW.DepthIntegrationProcess(model, integration_settings) + + hdf5_settings = KM.Parameters() + hdf5_settings.AddValue("model_part_name", settings["interface_model_part_name"]) + hdf5_settings.AddValue("file_settings", settings["file_settings"]) + hdf5_settings.AddValue("output_time_settings", settings["output_time_settings"]) + data_settings = KM.Parameters("""{"list_of_variables" : ["VELOCITY","HEIGHT"]}""") + if self.store_historical: + hdf5_settings.AddValue("nodal_solution_step_data_settings", data_settings) + else: + hdf5_settings.AddValue("nodal_data_value_settings", data_settings) + + self.hdf5_process = single_mesh_temporal_output_process.Factory(hdf5_settings, model) + + def Check(self): + self.integration_process.Check() + self.hdf5_process.Check() + + def ExecuteInitialize(self): + self.integration_process.ExecuteInitialize() + self.hdf5_process.ExecuteInitialize() + + def ExecuteBeforeSolutionLoop(self): + self.integration_process.ExecuteBeforeSolutionLoop() + self.hdf5_process.ExecuteBeforeSolutionLoop() + + def ExecuteInitializeSolutionStep(self): + self.integration_process.ExecuteInitialize() + self.hdf5_process.ExecuteInitialize() \ No newline at end of file From f57bac248204cc2ea0827fb526ea736c221738a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Tue, 26 Oct 2021 19:15:11 +0200 Subject: [PATCH 09/43] move default parmeters --- .../depth_integration_process.cpp | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index 85abcbcae521..5b9ef014ef00 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -28,6 +28,18 @@ namespace Kratos { +const Parameters DepthIntegrationProcess::GetDefaultParameters() const +{ + auto default_parameters = Parameters(R"( + { + "volume_model_part_name" : "", + "interface_model_part_name" : "", + "direction_of_integration" : [0.0, 0.0, 1.0], + "store_historical_database" : false + })"); + return default_parameters; +} + DepthIntegrationProcess::DepthIntegrationProcess( Model& rModel, Parameters ThisParameters @@ -178,16 +190,4 @@ int DepthIntegrationProcess::Check() return 0; } -const Parameters DepthIntegrationProcess::GetDefaultParameters() const -{ - auto default_parameters = Parameters(R"( - { - "volume_model_part_name" : "", - "interface_model_part_name" : "", - "direction_of_integration" : [0.0, 0.0, 1.0], - "store_historical_database" : false, - })"); - return default_parameters; -} - } // namespace Kratos. From 43b1da71bb6120a39672b06f0df783c5e28b30f5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Wed, 27 Oct 2021 16:14:24 +0200 Subject: [PATCH 10/43] minor fixes for the output --- .../depth_integration_process.cpp | 14 +++++++++-- .../depth_integration_process.h | 1 + .../depth_integration_output_process.py | 24 ++++++++++++++----- 3 files changed, 31 insertions(+), 8 deletions(-) diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index 5b9ef014ef00..6fa8a63095ef 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -18,6 +18,7 @@ // Project includes #include "includes/model_part.h" +#include "geometries/line_2d_2.h" #include "geometries/line_3d_2.h" #include "depth_integration_process.h" #include "utilities/variable_utils.h" @@ -57,6 +58,7 @@ DepthIntegrationProcess::DepthIntegrationProcess( } else { mpIntegrationModelPart = &rModel.CreateModelPart("integration_auxiliary_model_part"); } + mDimension = mrVolumeModelPart.GetProcessInfo()[DOMAIN_SIZE]; } void DepthIntegrationProcess::Execute() @@ -105,6 +107,8 @@ void DepthIntegrationProcess::Integrate(PointerVector& rObjec } velocity /= num_nodes; double height = max_elevation - min_elevation; + const array_1d momentum = height*velocity; + SetValue(rNode, MOMENTUM, momentum); SetValue(rNode, VELOCITY, velocity); SetValue(rNode, HEIGHT, height); } @@ -144,6 +148,9 @@ void DepthIntegrationProcess::InitializeIntegrationLine() VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Nodes()); mpIntegrationModelPart->RemoveNodesFromAllLevels(); + // Set the element name + std::string element_name = (mDimension == 2) ? "Element2D2N" : "Element3D2N"; + // Create the integration lines ModelPart::PropertiesType::Pointer p_prop; if (mpIntegrationModelPart->HasProperties(1)) { @@ -153,7 +160,7 @@ void DepthIntegrationProcess::InitializeIntegrationLine() } mpIntegrationModelPart->CreateNewNode(1, 0.0, 0.0, 0.0); mpIntegrationModelPart->CreateNewNode(2, 0.0, 0.0, 1.0); - mpIntegrationModelPart->CreateNewElement("Element3D2N", 1, {{1, 2}}, p_prop); + mpIntegrationModelPart->CreateNewElement(element_name, 1, {{1, 2}}, p_prop); } void DepthIntegrationProcess::SetIntegrationLine( @@ -181,7 +188,10 @@ Geometry>::Pointer DepthIntegrationProcess::CreateIntegrationLine( PointerVector> nodes; nodes.push_back(NodeType::Pointer(new NodeType(0, start))); nodes.push_back(NodeType::Pointer(new NodeType(0, end))); - return Kratos::make_shared>(nodes); + if (mDimension == 2) + return Kratos::make_shared>(nodes); + else + return Kratos::make_shared>(nodes); } int DepthIntegrationProcess::Check() diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h index 4e8ba3a5753d..d50189f1a8bb 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h @@ -151,6 +151,7 @@ class KRATOS_API(SHALLOW_WATER_APPLICATION) DepthIntegrationProcess : public Pro ModelPart* mpIntegrationModelPart; array_1d mDirection; bool mStoreHistorical; + int mDimension; ///@} ///@name Member Variables diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py index 76bf25b3b833..03b6681dddf9 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py @@ -24,7 +24,7 @@ def GetDefaultParameters(): "direction_of_integration" : [0.0, 0.0, 1.0], "interval" : [0.0,"End"], "file_settings" : {}, - "output_time_settings" : {}, + "output_time_settings" : {} }""") def __init__(self, model, settings): @@ -52,13 +52,15 @@ def __init__(self, model, settings): hdf5_settings.AddValue("model_part_name", settings["interface_model_part_name"]) hdf5_settings.AddValue("file_settings", settings["file_settings"]) hdf5_settings.AddValue("output_time_settings", settings["output_time_settings"]) - data_settings = KM.Parameters("""{"list_of_variables" : ["VELOCITY","HEIGHT"]}""") + data_settings = KM.Parameters("""{"list_of_variables" : ["MOMENTUM","VELOCITY","HEIGHT"]}""") if self.store_historical: hdf5_settings.AddValue("nodal_solution_step_data_settings", data_settings) else: hdf5_settings.AddValue("nodal_data_value_settings", data_settings) + hdf5_process_settings = KM.Parameters() + hdf5_process_settings.AddValue("Parameters", hdf5_settings) - self.hdf5_process = single_mesh_temporal_output_process.Factory(hdf5_settings, model) + self.hdf5_process = single_mesh_temporal_output_process.Factory(hdf5_process_settings, model) def Check(self): self.integration_process.Check() @@ -67,11 +69,21 @@ def Check(self): def ExecuteInitialize(self): self.integration_process.ExecuteInitialize() self.hdf5_process.ExecuteInitialize() + if not self.store_historical: + KM.VariableUtils().SetNonHistoricalVariableToZero(KM.VELOCITY, self.interface_model_part.Nodes) + KM.VariableUtils().SetNonHistoricalVariableToZero(KM.MOMENTUM, self.interface_model_part.Nodes) + KM.VariableUtils().SetNonHistoricalVariableToZero(SW.HEIGHT, self.interface_model_part.Nodes) def ExecuteBeforeSolutionLoop(self): self.integration_process.ExecuteBeforeSolutionLoop() self.hdf5_process.ExecuteBeforeSolutionLoop() - def ExecuteInitializeSolutionStep(self): - self.integration_process.ExecuteInitialize() - self.hdf5_process.ExecuteInitialize() \ No newline at end of file + def ExecuteBeforeOutputStep(self): + self.integration_process.Execute() + + def IsOutputStep(self): + # return self.hdf5_process.IsOutputstep() + return True + + def PrintOutput(self): + self.hdf5_process.ExecuteFinalizeSolutionStep() From 8fc1a63fe135aff5d8e17ffa1a123278297d646a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Wed, 27 Oct 2021 16:14:41 +0200 Subject: [PATCH 11/43] draft for input process --- .../depth_integration_input_process.py | 76 +++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py new file mode 100644 index 000000000000..5a44f3bf1204 --- /dev/null +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py @@ -0,0 +1,76 @@ +import KratosMultiphysics as KM +import KratosMultiphysics.ShallowWaterApplication as SW +from KratosMultiphysics.HDF5Application import single_mesh_temporal_input_process + +def Factory(settings, model): + if not isinstance(settings, KM.Parameters): + raise Exception("expected input shall be a Parameters object, encapsulating a json string") + return DepthIntegrationOutputProcess(model, settings["Parameters"]) + +class DepthIntegrationOutputProcess(KM.OutputProcess): + """DepthIntegrationOutputProcess + + This process performs a depth integration from a Navier-Stokes domain to a shallow water domain. + The depth integration values are stored in the nodes of the shallow water domain and + printed in HDF5 format. + """ + + @staticmethod + def GetDefaultParameters(): + return KM.Parameters("""{ + "model_part_name" : "", + "interface_model_part_name" : "interface_model_part", + "read_historical_database" : false, + "interval" : [0.0,"End"], + "invert_yz_axis" : false, + "ignore_vertical_component" : true, + "file_settings" : {} + }""") + + def __init__(self, model, settings): + """The constructor of the DepthIntegrationOutputProcess""" + + KM.OutputProcess.__init__(self) + + self.settings = settings + self.settings.ValidateAndAssignDefaults(self.GetDefaultParameters()) + + self.model_part = model[self.settings["volume_model_part_name"].GetString()] + self.read_historical = settings["read_historical_database"].GetBool() + self.interval = KM.IntervalUtility(settings) + interface_model_part_name = self.settings["volume_model_part_name"].GetString() + if model.HasModelPart(interface_model_part_name): + raise Exception("DepthIntegrationInputProcess: There is an existing interface model part with name '{}'".format(interface_model_part_name)) + self.interface_model_part = model.CreateModelPart(interface_model_part_name) + + hdf5_settings = KM.Parameters() + hdf5_settings.AddValue("model_part_name", settings["interface_model_part_name"]) + hdf5_settings.AddValue("file_settings", settings["file_settings"]) + data_settings = KM.Parameters("""{"list_of_variables" : ["MOMENTUM","VELOCITY","HEIGHT"]}""") + if self.read_historical: + hdf5_settings.AddValue("nodal_solution_step_data_settings", data_settings) + else: + hdf5_settings.AddValue("nodal_data_value_settings", data_settings) + hdf5_process_settings = KM.Parameters() + hdf5_process_settings.AddValue("Parameters", hdf5_settings) + + self.hdf5_process = single_mesh_temporal_input_process.Factory(hdf5_process_settings, model) + + def Check(self): + self.hdf5_process.Check() + + def ExecuteInitialize(self): + self.hdf5_process.ExecuteInitialize() + if not self.read_historical: + KM.VariableUtils().SetNonHistoricalVariableToZero(KM.VELOCITY, self.interface_model_part.Nodes) + KM.VariableUtils().SetNonHistoricalVariableToZero(KM.MOMENTUM, self.interface_model_part.Nodes) + KM.VariableUtils().SetNonHistoricalVariableToZero(SW.HEIGHT, self.interface_model_part.Nodes) + + def ExecuteBeforeSolutionLoop(self): + self.hdf5_process.ExecuteBeforeSolutionLoop() + + def ExecuteBeforeSolutionStep(self): + self.hdf5_process.ExecuteBeforeSolutionStep() + print(self.interface_model_part) + for node in self.interface_model_part.Nodes: + print(node.GetSolutionStepValue(KM.VELOCITY)) From 90c31de7838f1cda1eeb98db3137d57efb50ee8a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Thu, 28 Oct 2021 11:40:00 +0200 Subject: [PATCH 12/43] some docstrings --- .../depth_integration_output_process.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py index 03b6681dddf9..992d9eb646cc 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py @@ -31,12 +31,10 @@ def __init__(self, model, settings): """The constructor of the DepthIntegrationOutputProcess""" KM.OutputProcess.__init__(self) + settings.ValidateAndAssignDefaults(self.GetDefaultParameters()) - self.settings = settings - self.settings.ValidateAndAssignDefaults(self.GetDefaultParameters()) - - self.volume_model_part = model[self.settings["volume_model_part_name"].GetString()] - self.interface_model_part = model[self.settings["interface_model_part_name"].GetString()] + self.volume_model_part = model[settings["volume_model_part_name"].GetString()] + self.interface_model_part = model[settings["interface_model_part_name"].GetString()] self.store_historical = settings["store_historical_database"].GetBool() self.interval = KM.IntervalUtility(settings) @@ -63,27 +61,30 @@ def __init__(self, model, settings): self.hdf5_process = single_mesh_temporal_output_process.Factory(hdf5_process_settings, model) def Check(self): + '''Check the processes.''' self.integration_process.Check() self.hdf5_process.Check() def ExecuteInitialize(self): - self.integration_process.ExecuteInitialize() - self.hdf5_process.ExecuteInitialize() + '''Initialize the variables.''' if not self.store_historical: KM.VariableUtils().SetNonHistoricalVariableToZero(KM.VELOCITY, self.interface_model_part.Nodes) KM.VariableUtils().SetNonHistoricalVariableToZero(KM.MOMENTUM, self.interface_model_part.Nodes) KM.VariableUtils().SetNonHistoricalVariableToZero(SW.HEIGHT, self.interface_model_part.Nodes) def ExecuteBeforeSolutionLoop(self): - self.integration_process.ExecuteBeforeSolutionLoop() + '''Write the interface_model_part in HDF5 format.''' self.hdf5_process.ExecuteBeforeSolutionLoop() def ExecuteBeforeOutputStep(self): + '''Perform the depth integration over the interface_model_part.''' self.integration_process.Execute() def IsOutputStep(self): + '''IsOutputStep.''' # return self.hdf5_process.IsOutputstep() return True def PrintOutput(self): + '''Print the integrated variables from interface_model_part.''' self.hdf5_process.ExecuteFinalizeSolutionStep() From c5f3c54c99c40e82df9a4fa2958fe13b52cd1e68 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Thu, 28 Oct 2021 11:41:07 +0200 Subject: [PATCH 13/43] add new variables tool --- .../add_custom_utilities_to_python.cpp | 4 +++ .../shallow_water_utilities.h | 25 +++++++++++++++++++ 2 files changed, 29 insertions(+) diff --git a/applications/ShallowWaterApplication/custom_python/add_custom_utilities_to_python.cpp b/applications/ShallowWaterApplication/custom_python/add_custom_utilities_to_python.cpp index 67522810293c..8e682294101d 100644 --- a/applications/ShallowWaterApplication/custom_python/add_custom_utilities_to_python.cpp +++ b/applications/ShallowWaterApplication/custom_python/add_custom_utilities_to_python.cpp @@ -105,6 +105,10 @@ void AddCustomUtilitiesToPython(pybind11::module& m) .def("ComputeL2Norm", &ShallowWaterUtilities::ComputeL2NormAABB) .def("ComputeL2NormNonHistorical", &ShallowWaterUtilities::ComputeL2Norm) .def("ComputeL2NormNonHistorical", &ShallowWaterUtilities::ComputeL2NormAABB) + .def("SwapYZComponents", &ShallowWaterUtilities::SwapYZComponents) + .def("SwapYZComponentsNonHistorical", &ShallowWaterUtilities::SwapYZComponentsNonHistorical) + .def("SwapYZComponentsNonHistorical", &ShallowWaterUtilities::SwapYZComponentsNonHistorical) + .def("SwapYZComponentsNonHistorical", &ShallowWaterUtilities::SwapYZComponentsNonHistorical) .def("ComputeHydrostaticForces", ComputeHydrostaticForces1) .def("ComputeHydrostaticForces", ComputeHydrostaticForces2) .def("ComputeHydrostaticForces", ComputeHydrostaticForces1) diff --git a/applications/ShallowWaterApplication/custom_utilities/shallow_water_utilities.h b/applications/ShallowWaterApplication/custom_utilities/shallow_water_utilities.h index 7d9febfb93e7..fbab1247960d 100644 --- a/applications/ShallowWaterApplication/custom_utilities/shallow_water_utilities.h +++ b/applications/ShallowWaterApplication/custom_utilities/shallow_water_utilities.h @@ -65,6 +65,8 @@ class KRATOS_API(SHALLOW_WATER_APPLICATION) ShallowWaterUtilities typedef Geometry GeometryType; + typedef ModelPart::NodesContainerType NodesContainerType; + ///@} ///@name Pointer definition ///@{ @@ -140,6 +142,29 @@ class KRATOS_API(SHALLOW_WATER_APPLICATION) ShallowWaterUtilities */ void SetMeshZCoordinate(ModelPart& rModelPart, const Variable& rVariable); + /** + * @brief Swap the Y and Z components of a vector variable + */ + void SwapYZComponents(const Variable>& rVariable, NodesContainerType& rNodes) + { + block_for_each(rNodes, [&](NodeType& rNode){ + array_1d& r_value = rNode.FastGetSolutionStepValue(rVariable); + std::swap(r_value[1], r_value[2]); + }); + } + + /** + * @brief Swap the Y and Z components of a vector variable + */ + template + void SwapYZComponentsNonHistorical(const Variable>& rVariable, TContainerType& rContainer) + { + block_for_each(rContainer, [&](typename TContainerType::value_type& rEntity){ + array_1d& r_value = rEntity.GetValue(rVariable); + std::swap(r_value[1], r_value[2]); + }); + } + /** * @brief Compute the L-2 norm for the given double variable */ From 7a0b718008a37632d61d867592132ce5a71f0889 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Thu, 28 Oct 2021 13:12:42 +0200 Subject: [PATCH 14/43] wip --- .../depth_integration_input_process.py | 74 ++++++++++++------- 1 file changed, 49 insertions(+), 25 deletions(-) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py index 5a44f3bf1204..799f74317905 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py @@ -1,18 +1,19 @@ import KratosMultiphysics as KM import KratosMultiphysics.ShallowWaterApplication as SW +import KratosMultiphysics.PfemFluidDynamicsApplication as PFEM +import KratosMultiphysics.DelaunayMeshingApplication as Delaunay +from KratosMultiphysics.HDF5Application import read_model_part_from_hdf5_process from KratosMultiphysics.HDF5Application import single_mesh_temporal_input_process def Factory(settings, model): if not isinstance(settings, KM.Parameters): raise Exception("expected input shall be a Parameters object, encapsulating a json string") - return DepthIntegrationOutputProcess(model, settings["Parameters"]) + return DepthIntegrationInputProcess(model, settings["Parameters"]) -class DepthIntegrationOutputProcess(KM.OutputProcess): - """DepthIntegrationOutputProcess +class DepthIntegrationInputProcess(KM.OutputProcess): + """DepthIntegrationInputProcess - This process performs a depth integration from a Navier-Stokes domain to a shallow water domain. - The depth integration values are stored in the nodes of the shallow water domain and - printed in HDF5 format. + Read the depth integrated values from an HDF5 file and set them as boundary conditions. """ @staticmethod @@ -22,26 +23,26 @@ def GetDefaultParameters(): "interface_model_part_name" : "interface_model_part", "read_historical_database" : false, "interval" : [0.0,"End"], - "invert_yz_axis" : false, + "swap_yz_axis" : false, "ignore_vertical_component" : true, "file_settings" : {} }""") def __init__(self, model, settings): - """The constructor of the DepthIntegrationOutputProcess""" + """The constructor of the DepthIntegrationInputProcess""" KM.OutputProcess.__init__(self) + settings.ValidateAndAssignDefaults(self.GetDefaultParameters()) - self.settings = settings - self.settings.ValidateAndAssignDefaults(self.GetDefaultParameters()) - - self.model_part = model[self.settings["volume_model_part_name"].GetString()] + self.model_part = model[settings["model_part_name"].GetString()] self.read_historical = settings["read_historical_database"].GetBool() self.interval = KM.IntervalUtility(settings) - interface_model_part_name = self.settings["volume_model_part_name"].GetString() + interface_model_part_name = settings["interface_model_part_name"].GetString() if model.HasModelPart(interface_model_part_name): raise Exception("DepthIntegrationInputProcess: There is an existing interface model part with name '{}'".format(interface_model_part_name)) self.interface_model_part = model.CreateModelPart(interface_model_part_name) + self.swap_yz_axis = settings["swap_yz_axis"].GetBool() + self.ignore_vertical_component = settings["ignore_vertical_component"].GetBool() hdf5_settings = KM.Parameters() hdf5_settings.AddValue("model_part_name", settings["interface_model_part_name"]) @@ -53,24 +54,47 @@ def __init__(self, model, settings): hdf5_settings.AddValue("nodal_data_value_settings", data_settings) hdf5_process_settings = KM.Parameters() hdf5_process_settings.AddValue("Parameters", hdf5_settings) - + + self.hdf5_read = read_model_part_from_hdf5_process.Factory(hdf5_process_settings.Clone(), model) self.hdf5_process = single_mesh_temporal_input_process.Factory(hdf5_process_settings, model) def Check(self): + '''Check the processes.''' + self.hdf5_read.Check() self.hdf5_process.Check() def ExecuteInitialize(self): - self.hdf5_process.ExecuteInitialize() - if not self.read_historical: - KM.VariableUtils().SetNonHistoricalVariableToZero(KM.VELOCITY, self.interface_model_part.Nodes) - KM.VariableUtils().SetNonHistoricalVariableToZero(KM.MOMENTUM, self.interface_model_part.Nodes) - KM.VariableUtils().SetNonHistoricalVariableToZero(SW.HEIGHT, self.interface_model_part.Nodes) + '''Read the interface_model_part and set the variables.''' + self.hdf5_read.ExecuteInitialize() + self._CheckInputVariables() + self._MapToBoundaryCondition() - def ExecuteBeforeSolutionLoop(self): - self.hdf5_process.ExecuteBeforeSolutionLoop() + def ExecuteInitializeSolutionStep(self): + '''Set the variables in interface_model_part at the current time step.''' + self.hdf5_process.ExecuteInitializeSolutionStep() # look for the current time + self._CheckInputVariables() + self._MapToBoundaryCondition() - def ExecuteBeforeSolutionStep(self): - self.hdf5_process.ExecuteBeforeSolutionStep() - print(self.interface_model_part) + ### + print("ExecuteInitializeSolutionStep") for node in self.interface_model_part.Nodes: - print(node.GetSolutionStepValue(KM.VELOCITY)) + print(node.GetValue(KM.VELOCITY)) + + def _CheckInputVariables(self): + if self.swap_yz_axis: + if self.read_historical: + SW.ShallowWaterUtilities().SwapYZComponents(KM.MOMENTUM, self.interface_model_part.Nodes) + SW.ShallowWaterUtilities().SwapYZComponents(KM.VELOCITY, self.interface_model_part.Nodes) + else: + SW.ShallowWaterUtilities().SwapYZComponentsNonHistorical(KM.MOMENTUM, self.interface_model_part.Nodes) + SW.ShallowWaterUtilities().SwapYZComponentsNonHistorical(KM.VELOCITY, self.interface_model_part.Nodes) + if self.ignore_vertical_component: + if self.read_historical: + KM.VariableUtils().SetVariableToZero(KM.MOMENTUM_Z, self.interface_model_part.Nodes) + KM.VariableUtils().SetVariableToZero(KM.VELOCITY_Z, self.interface_model_part.Nodes) + else: + KM.VariableUtils().SetNonHistoricalVariableToZero(KM.MOMENTUM_Z, self.interface_model_part.Nodes) + KM.VariableUtils().SetNonHistoricalVariableToZero(KM.VELOCITY_Z, self.interface_model_part.Nodes) + + def _MapToBoundaryCondition(self): + pass From ef461a61774636b498e1cc9e72817c43dbf0f0d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Fri, 29 Oct 2021 12:04:04 +0200 Subject: [PATCH 15/43] read hdf5 file closest time --- .../depth_integration_input_process.py | 36 +++++++++++++++++-- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py index 799f74317905..c485288dec74 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py @@ -1,9 +1,12 @@ +from re import A import KratosMultiphysics as KM +from KratosMultiphysics.HDF5Application.core.controllers import DefaultController import KratosMultiphysics.ShallowWaterApplication as SW import KratosMultiphysics.PfemFluidDynamicsApplication as PFEM import KratosMultiphysics.DelaunayMeshingApplication as Delaunay from KratosMultiphysics.HDF5Application import read_model_part_from_hdf5_process from KratosMultiphysics.HDF5Application import single_mesh_temporal_input_process +from os import path, listdir def Factory(settings, model): if not isinstance(settings, KM.Parameters): @@ -56,7 +59,8 @@ def __init__(self, model, settings): hdf5_process_settings.AddValue("Parameters", hdf5_settings) self.hdf5_read = read_model_part_from_hdf5_process.Factory(hdf5_process_settings.Clone(), model) - self.hdf5_process = single_mesh_temporal_input_process.Factory(hdf5_process_settings, model) + self.hdf5_process = single_mesh_temporal_input_process.Factory(hdf5_process_settings.Clone(), model) + self._GetInputTimes(settings['file_settings']) def Check(self): '''Check the processes.''' @@ -70,8 +74,9 @@ def ExecuteInitialize(self): self._MapToBoundaryCondition() def ExecuteInitializeSolutionStep(self): - '''Set the variables in interface_model_part at the current time step.''' - self.hdf5_process.ExecuteInitializeSolutionStep() # look for the current time + '''Set the variables in the interface_model_part at the current time.''' + self._SetCurrentTime() + self.hdf5_process.ExecuteInitializeSolutionStep() self._CheckInputVariables() self._MapToBoundaryCondition() @@ -80,6 +85,31 @@ def ExecuteInitializeSolutionStep(self): for node in self.interface_model_part.Nodes: print(node.GetValue(KM.VELOCITY)) + def _GetInputTimes(self, file_settings): + # Get all the file names + file_name = file_settings["file_name"].GetString() + folder_name = path.dirname(file_name) + file_names = [f for f in listdir(folder_name) if path.isfile(path.join(folder_name, f))] + + # Find the common parts of the found names and the file name. The difference are the times + self.file_pattern = path.basename(file_name) + self.file_pattern = self.file_pattern.replace('', self.interface_model_part.Name) + prefix = path.commonprefix([self.file_pattern, file_names[0]]) + suffix = self.file_pattern.replace(prefix, '') + suffix = path.commonprefix([''.join(reversed(suffix)), ''.join(reversed(file_names[0]))]) + suffix = ''.join(reversed(suffix)) + self.times = [] + for f in file_names: + f = f.replace(prefix, '') + f = f.replace(suffix, '') + self.times.append(float(f)) + self.times.sort() + + def _SetCurrentTime(self): + current_time = self.model_part.ProcessInfo.GetValue(KM.TIME) + closest_time = next(filter(lambda x: x>current_time, self.times)) + self.interface_model_part.ProcessInfo.SetValue(KM.TIME, closest_time) + def _CheckInputVariables(self): if self.swap_yz_axis: if self.read_historical: From 6fba32bd8bc6608ef4d078344fb9bec3905074eb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Fri, 29 Oct 2021 12:26:15 +0200 Subject: [PATCH 16/43] map to boundary conditions --- .../depth_integration_input_process.py | 37 +++++++++++++++---- 1 file changed, 29 insertions(+), 8 deletions(-) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py index c485288dec74..2edb6641c23e 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py @@ -1,9 +1,10 @@ from re import A import KratosMultiphysics as KM -from KratosMultiphysics.HDF5Application.core.controllers import DefaultController import KratosMultiphysics.ShallowWaterApplication as SW +import KratosMultiphysics.MappingApplication as Mapping import KratosMultiphysics.PfemFluidDynamicsApplication as PFEM import KratosMultiphysics.DelaunayMeshingApplication as Delaunay +from KratosMultiphysics.kratos_utilities import GenerateVariableListFromInput from KratosMultiphysics.HDF5Application import read_model_part_from_hdf5_process from KratosMultiphysics.HDF5Application import single_mesh_temporal_input_process from os import path, listdir @@ -26,6 +27,7 @@ def GetDefaultParameters(): "interface_model_part_name" : "interface_model_part", "read_historical_database" : false, "interval" : [0.0,"End"], + "list_of_variables" : ["MOMENTUM"], "swap_yz_axis" : false, "ignore_vertical_component" : true, "file_settings" : {} @@ -46,6 +48,7 @@ def __init__(self, model, settings): self.interface_model_part = model.CreateModelPart(interface_model_part_name) self.swap_yz_axis = settings["swap_yz_axis"].GetBool() self.ignore_vertical_component = settings["ignore_vertical_component"].GetBool() + self.variables = GenerateVariableListFromInput(settings["list_of_variables"]) hdf5_settings = KM.Parameters() hdf5_settings.AddValue("model_part_name", settings["interface_model_part_name"]) @@ -71,19 +74,22 @@ def ExecuteInitialize(self): '''Read the interface_model_part and set the variables.''' self.hdf5_read.ExecuteInitialize() self._CheckInputVariables() + self._CreateMapper() self._MapToBoundaryCondition() def ExecuteInitializeSolutionStep(self): '''Set the variables in the interface_model_part at the current time.''' - self._SetCurrentTime() - self.hdf5_process.ExecuteInitializeSolutionStep() - self._CheckInputVariables() - self._MapToBoundaryCondition() + current_time = self.model_part.ProcessInfo.GetValue(KM.TIME) + if self.interval.IsInInterval(current_time): + self._SetCurrentTime() + self.hdf5_process.ExecuteInitializeSolutionStep() + self._CheckInputVariables() + self._MapToBoundaryCondition() - ### + ###################################### print("ExecuteInitializeSolutionStep") for node in self.interface_model_part.Nodes: - print(node.GetValue(KM.VELOCITY)) + print(node.GetValue(KM.MOMENTUM)) def _GetInputTimes(self, file_settings): # Get all the file names @@ -127,4 +133,19 @@ def _CheckInputVariables(self): KM.VariableUtils().SetNonHistoricalVariableToZero(KM.VELOCITY_Z, self.interface_model_part.Nodes) def _MapToBoundaryCondition(self): - pass + for variable in self.variables: + if self.read_historical: + self.mapper.Map(variable, variable) + else: + self.mapper.Map(variable, variable, KM.Mapper.FROM_NON_HISTORICAL) + + def _CreateMapper(self): + mapper_settings = KM.Parameters("""{ + "mapper_type": "nearest_neighbor", + "echo_level" : 0 + }""") + + self.mapper = KM.MapperFactory.CreateMapper( + self.interface_model_part, + self.model_part, + mapper_settings) From 1265bb1b1729bea32ba0987fb0f83147c4db137c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Fri, 29 Oct 2021 15:33:45 +0200 Subject: [PATCH 17/43] fixity and serach radius --- .../depth_integration_input_process.py | 20 +++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py index 2edb6641c23e..2aa06d5468f5 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py @@ -1,4 +1,3 @@ -from re import A import KratosMultiphysics as KM import KratosMultiphysics.ShallowWaterApplication as SW import KratosMultiphysics.MappingApplication as Mapping @@ -28,6 +27,7 @@ def GetDefaultParameters(): "read_historical_database" : false, "interval" : [0.0,"End"], "list_of_variables" : ["MOMENTUM"], + "list_of_variables_to_fix" : ["MOMENTUM_X","MOMENTUM_Y"], "swap_yz_axis" : false, "ignore_vertical_component" : true, "file_settings" : {} @@ -49,6 +49,7 @@ def __init__(self, model, settings): self.swap_yz_axis = settings["swap_yz_axis"].GetBool() self.ignore_vertical_component = settings["ignore_vertical_component"].GetBool() self.variables = GenerateVariableListFromInput(settings["list_of_variables"]) + self.variables_to_fix = GenerateVariableListFromInput(settings["list_of_variables_to_fix"]) hdf5_settings = KM.Parameters() hdf5_settings.AddValue("model_part_name", settings["interface_model_part_name"]) @@ -139,11 +140,26 @@ def _MapToBoundaryCondition(self): else: self.mapper.Map(variable, variable, KM.Mapper.FROM_NON_HISTORICAL) + for variable in self.variables_to_fix: + KM.VariableUtils().ApplyFixity(variable, True, self.model_part.Nodes) + def _CreateMapper(self): mapper_settings = KM.Parameters("""{ "mapper_type": "nearest_neighbor", - "echo_level" : 0 + "echo_level" : 0, + "search_settings" : { + "search_radius" : 0.0 + } }""") + min_point = KM.Point([ 1e6, 1e6, 1e6]) + max_point = KM.Point([-1e6, -1e6, -1e6]) + for node in self.model_part.Nodes: + for i in range(3): + point = KM.Point(node) + min_point[i] = min([min_point[i], point[i]]) + max_point[i] = max([max_point[i], point[i]]) + distance = 1.05 * (max_point - min_point).norm_2() + mapper_settings["search_settings"]["search_radius"].SetDouble(distance) self.mapper = KM.MapperFactory.CreateMapper( self.interface_model_part, From af675009cbb61233c227c4afc4df2aba1147276b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Fri, 29 Oct 2021 15:46:40 +0200 Subject: [PATCH 18/43] settings --- .../depth_integration_input_process.py | 58 +++++++++---------- 1 file changed, 26 insertions(+), 32 deletions(-) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py index 2aa06d5468f5..ba5751cba514 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py @@ -22,40 +22,39 @@ class DepthIntegrationInputProcess(KM.OutputProcess): @staticmethod def GetDefaultParameters(): return KM.Parameters("""{ - "model_part_name" : "", - "interface_model_part_name" : "interface_model_part", - "read_historical_database" : false, - "interval" : [0.0,"End"], - "list_of_variables" : ["MOMENTUM"], - "list_of_variables_to_fix" : ["MOMENTUM_X","MOMENTUM_Y"], - "swap_yz_axis" : false, - "ignore_vertical_component" : true, - "file_settings" : {} + "model_part_name" : "", + "interface_model_part_name" : "interface_model_part", + "read_historical_database" : false, + "interval" : [0.0,"End"], + "list_of_variables" : ["MOMENTUM"], + "list_of_variables_to_fix" : ["MOMENTUM_X","MOMENTUM_Y"], + "default_time_after_interval" : 0.0, + "swap_yz_axis" : false, + "ignore_vertical_component" : true, + "file_settings" : {} }""") def __init__(self, model, settings): """The constructor of the DepthIntegrationInputProcess""" KM.OutputProcess.__init__(self) - settings.ValidateAndAssignDefaults(self.GetDefaultParameters()) + self.settings = settings + self.settings.ValidateAndAssignDefaults(self.GetDefaultParameters()) - self.model_part = model[settings["model_part_name"].GetString()] - self.read_historical = settings["read_historical_database"].GetBool() - self.interval = KM.IntervalUtility(settings) - interface_model_part_name = settings["interface_model_part_name"].GetString() + self.model_part = model[self.settings["model_part_name"].GetString()] + self.interval = KM.IntervalUtility(self.settings) + interface_model_part_name = self.settings["interface_model_part_name"].GetString() if model.HasModelPart(interface_model_part_name): raise Exception("DepthIntegrationInputProcess: There is an existing interface model part with name '{}'".format(interface_model_part_name)) self.interface_model_part = model.CreateModelPart(interface_model_part_name) - self.swap_yz_axis = settings["swap_yz_axis"].GetBool() - self.ignore_vertical_component = settings["ignore_vertical_component"].GetBool() - self.variables = GenerateVariableListFromInput(settings["list_of_variables"]) - self.variables_to_fix = GenerateVariableListFromInput(settings["list_of_variables_to_fix"]) + self.variables = GenerateVariableListFromInput(self.settings["list_of_variables"]) + self.variables_to_fix = GenerateVariableListFromInput(self.settings["list_of_variables_to_fix"]) hdf5_settings = KM.Parameters() - hdf5_settings.AddValue("model_part_name", settings["interface_model_part_name"]) - hdf5_settings.AddValue("file_settings", settings["file_settings"]) + hdf5_settings.AddString("model_part_name", interface_model_part_name) + hdf5_settings.AddValue("file_settings", self.settings["file_settings"]) data_settings = KM.Parameters("""{"list_of_variables" : ["MOMENTUM","VELOCITY","HEIGHT"]}""") - if self.read_historical: + if self.settings["read_historical_database"].GetBool(): hdf5_settings.AddValue("nodal_solution_step_data_settings", data_settings) else: hdf5_settings.AddValue("nodal_data_value_settings", data_settings) @@ -64,7 +63,7 @@ def __init__(self, model, settings): self.hdf5_read = read_model_part_from_hdf5_process.Factory(hdf5_process_settings.Clone(), model) self.hdf5_process = single_mesh_temporal_input_process.Factory(hdf5_process_settings.Clone(), model) - self._GetInputTimes(settings['file_settings']) + self._GetInputTimes(self.settings['file_settings']) def Check(self): '''Check the processes.''' @@ -87,11 +86,6 @@ def ExecuteInitializeSolutionStep(self): self._CheckInputVariables() self._MapToBoundaryCondition() - ###################################### - print("ExecuteInitializeSolutionStep") - for node in self.interface_model_part.Nodes: - print(node.GetValue(KM.MOMENTUM)) - def _GetInputTimes(self, file_settings): # Get all the file names file_name = file_settings["file_name"].GetString() @@ -118,15 +112,15 @@ def _SetCurrentTime(self): self.interface_model_part.ProcessInfo.SetValue(KM.TIME, closest_time) def _CheckInputVariables(self): - if self.swap_yz_axis: - if self.read_historical: + if self.settings["swap_yz_axis"].GetBool(): + if self.settings["read_historical_database"].GetBool(): SW.ShallowWaterUtilities().SwapYZComponents(KM.MOMENTUM, self.interface_model_part.Nodes) SW.ShallowWaterUtilities().SwapYZComponents(KM.VELOCITY, self.interface_model_part.Nodes) else: SW.ShallowWaterUtilities().SwapYZComponentsNonHistorical(KM.MOMENTUM, self.interface_model_part.Nodes) SW.ShallowWaterUtilities().SwapYZComponentsNonHistorical(KM.VELOCITY, self.interface_model_part.Nodes) - if self.ignore_vertical_component: - if self.read_historical: + if self.settings["ignore_vertical_component"].GetBool(): + if self.settings["read_historical_database"].GetBool(): KM.VariableUtils().SetVariableToZero(KM.MOMENTUM_Z, self.interface_model_part.Nodes) KM.VariableUtils().SetVariableToZero(KM.VELOCITY_Z, self.interface_model_part.Nodes) else: @@ -135,7 +129,7 @@ def _CheckInputVariables(self): def _MapToBoundaryCondition(self): for variable in self.variables: - if self.read_historical: + if self.settings["read_historical_database"].GetBool(): self.mapper.Map(variable, variable) else: self.mapper.Map(variable, variable, KM.Mapper.FROM_NON_HISTORICAL) From 2c95d25d853a389af964985e63b2fc5a2a3d5918 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Fri, 29 Oct 2021 16:51:41 +0200 Subject: [PATCH 19/43] add smoothing function after interval --- .../depth_integration_input_process.py | 22 +++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py index ba5751cba514..243dba223e06 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py @@ -7,6 +7,7 @@ from KratosMultiphysics.HDF5Application import read_model_part_from_hdf5_process from KratosMultiphysics.HDF5Application import single_mesh_temporal_input_process from os import path, listdir +from math import exp def Factory(settings, model): if not isinstance(settings, KM.Parameters): @@ -29,6 +30,7 @@ def GetDefaultParameters(): "list_of_variables" : ["MOMENTUM"], "list_of_variables_to_fix" : ["MOMENTUM_X","MOMENTUM_Y"], "default_time_after_interval" : 0.0, + "semi_period_after_interval" : 1.0, "swap_yz_axis" : false, "ignore_vertical_component" : true, "file_settings" : {} @@ -85,6 +87,12 @@ def ExecuteInitializeSolutionStep(self): self.hdf5_process.ExecuteInitializeSolutionStep() self._CheckInputVariables() self._MapToBoundaryCondition() + else: + self._SetDefaultTime() + self.hdf5_process.ExecuteInitializeSolutionStep() + self._CheckInputVariables() + self._MapToBoundaryCondition() + self._SmoothDefaultValue() def _GetInputTimes(self, file_settings): # Get all the file names @@ -111,6 +119,10 @@ def _SetCurrentTime(self): closest_time = next(filter(lambda x: x>current_time, self.times)) self.interface_model_part.ProcessInfo.SetValue(KM.TIME, closest_time) + def _SetDefaultTime(self): + default_time = self.settings["default_time_after_interval"].GetDouble() + self.interface_model_part.ProcessInfo.SetValue(KM.TIME, default_time) + def _CheckInputVariables(self): if self.settings["swap_yz_axis"].GetBool(): if self.settings["read_historical_database"].GetBool(): @@ -159,3 +171,13 @@ def _CreateMapper(self): self.interface_model_part, self.model_part, mapper_settings) + + def _SmoothDefaultValue(self): + elapsed_time = self.model_part.ProcessInfo.GetValue(KM.DELTA_TIME) + semi_period = self.settings["semi_period_after_interval"].GetDouble() + smooth = exp(elapsed_time / semi_period) + for node in self.model_part.Nodes: + for variable in self.variables: + initial = node.GetSolutionStepValue(variable, 1) + increment = node.GetSolutionStepValue(variable) - initial + node.SetSolutionStepValue(variable, initial + increment * smooth) From 94de459c7cc2b123b3baa84c6bc88a34fe813064 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Mon, 8 Nov 2021 14:07:12 +0100 Subject: [PATCH 20/43] avoid dependencies between simulations --- .../depth_integration_output_process.py | 52 +++++++++++++++---- 1 file changed, 43 insertions(+), 9 deletions(-) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py index 992d9eb646cc..dabb7c83b896 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py @@ -20,6 +20,7 @@ def GetDefaultParameters(): return KM.Parameters("""{ "volume_model_part_name" : "", "interface_model_part_name" : "", + "output_model_part_name" : "", "store_historical_database" : false, "direction_of_integration" : [0.0, 0.0, 1.0], "interval" : [0.0,"End"], @@ -35,6 +36,7 @@ def __init__(self, model, settings): self.volume_model_part = model[settings["volume_model_part_name"].GetString()] self.interface_model_part = model[settings["interface_model_part_name"].GetString()] + self.output_model_part = model.CreateModelPart(settings["output_model_part_name"].GetString()) self.store_historical = settings["store_historical_database"].GetBool() self.interval = KM.IntervalUtility(settings) @@ -47,7 +49,7 @@ def __init__(self, model, settings): self.integration_process = SW.DepthIntegrationProcess(model, integration_settings) hdf5_settings = KM.Parameters() - hdf5_settings.AddValue("model_part_name", settings["interface_model_part_name"]) + hdf5_settings.AddValue("model_part_name", settings["output_model_part_name"]) hdf5_settings.AddValue("file_settings", settings["file_settings"]) hdf5_settings.AddValue("output_time_settings", settings["output_time_settings"]) data_settings = KM.Parameters("""{"list_of_variables" : ["MOMENTUM","VELOCITY","HEIGHT"]}""") @@ -57,34 +59,66 @@ def __init__(self, model, settings): hdf5_settings.AddValue("nodal_data_value_settings", data_settings) hdf5_process_settings = KM.Parameters() hdf5_process_settings.AddValue("Parameters", hdf5_settings) - + self.hdf5_process = single_mesh_temporal_output_process.Factory(hdf5_process_settings, model) + def Check(self): '''Check the processes.''' self.integration_process.Check() self.hdf5_process.Check() + def ExecuteInitialize(self): - '''Initialize the variables.''' + '''Initialize the output model part.''' + self._InitializeOutputModelPart() + self._SetOutputProcessInfo() if not self.store_historical: - KM.VariableUtils().SetNonHistoricalVariableToZero(KM.VELOCITY, self.interface_model_part.Nodes) - KM.VariableUtils().SetNonHistoricalVariableToZero(KM.MOMENTUM, self.interface_model_part.Nodes) - KM.VariableUtils().SetNonHistoricalVariableToZero(SW.HEIGHT, self.interface_model_part.Nodes) + variables = [KM.VELOCITY, KM.MOMENTUM, SW.HEIGHT] + for var in variables: + KM.VariableUtils().SetNonHistoricalVariableToZero(var, self.interface_model_part.Nodes) + KM.VariableUtils().SetNonHistoricalVariableToZero(var, self.output_model_part.Nodes) + def ExecuteBeforeSolutionLoop(self): - '''Write the interface_model_part in HDF5 format.''' + '''Write the interface model part in HDF5 format.''' + self._SetOutputProcessInfo() self.hdf5_process.ExecuteBeforeSolutionLoop() + + def ExecuteInitializeSolutionStep(self): + '''Synchronize the ProcessInfo of the output and interface model part.''' + self._SetOutputProcessInfo() + + def ExecuteBeforeOutputStep(self): - '''Perform the depth integration over the interface_model_part.''' + '''Perform the depth integration over the interface model part.''' self.integration_process.Execute() + def IsOutputStep(self): '''IsOutputStep.''' # return self.hdf5_process.IsOutputstep() return True + def PrintOutput(self): - '''Print the integrated variables from interface_model_part.''' + '''Print the integrated variables from interface model part.''' self.hdf5_process.ExecuteFinalizeSolutionStep() + + + def _InitializeOutputModelPart(self): + KM.MergeVariableListsUtility().Merge(self.interface_model_part, self.output_model_part) + domain_size = self.volume_model_part.ProcessInfo[KM.DOMAIN_SIZE] + element_name = "Element{}D2N".format(domain_size) + condition_name = "LineCondition{}D2N".format(domain_size) + KM.ConnectivityPreserveModeler().GenerateModelPart( + self.interface_model_part, self.output_model_part, element_name, condition_name) + self.output_model_part.ProcessInfo.Clear() + + + def _SetOutputProcessInfo(self): + time = self.interface_model_part.ProcessInfo[KM.TIME] + step = self.interface_model_part.ProcessInfo[KM.STEP] + self.output_model_part.ProcessInfo[KM.TIME] = time + self.output_model_part.ProcessInfo[KM.STEP] = step From 8dd4a8dffccadbed3cf83d3c035523fe97c90d48 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Tue, 9 Nov 2021 12:37:55 +0100 Subject: [PATCH 21/43] create copy of process info --- .../python_scripts/coupling/depth_integration_output_process.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py index dabb7c83b896..de5d90aa4dad 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py @@ -114,7 +114,7 @@ def _InitializeOutputModelPart(self): condition_name = "LineCondition{}D2N".format(domain_size) KM.ConnectivityPreserveModeler().GenerateModelPart( self.interface_model_part, self.output_model_part, element_name, condition_name) - self.output_model_part.ProcessInfo.Clear() + self.output_model_part.ProcessInfo = KM.ProcessInfo() def _SetOutputProcessInfo(self): From 8334e6188a8b881f3a5e3a784d5d38621e497bff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Tue, 9 Nov 2021 12:48:48 +0100 Subject: [PATCH 22/43] style --- .../coupling/depth_integration_input_process.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py index 243dba223e06..849a33c15c93 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py @@ -67,11 +67,13 @@ def __init__(self, model, settings): self.hdf5_process = single_mesh_temporal_input_process.Factory(hdf5_process_settings.Clone(), model) self._GetInputTimes(self.settings['file_settings']) + def Check(self): '''Check the processes.''' self.hdf5_read.Check() self.hdf5_process.Check() + def ExecuteInitialize(self): '''Read the interface_model_part and set the variables.''' self.hdf5_read.ExecuteInitialize() @@ -79,6 +81,7 @@ def ExecuteInitialize(self): self._CreateMapper() self._MapToBoundaryCondition() + def ExecuteInitializeSolutionStep(self): '''Set the variables in the interface_model_part at the current time.''' current_time = self.model_part.ProcessInfo.GetValue(KM.TIME) @@ -94,6 +97,7 @@ def ExecuteInitializeSolutionStep(self): self._MapToBoundaryCondition() self._SmoothDefaultValue() + def _GetInputTimes(self, file_settings): # Get all the file names file_name = file_settings["file_name"].GetString() @@ -114,15 +118,18 @@ def _GetInputTimes(self, file_settings): self.times.append(float(f)) self.times.sort() + def _SetCurrentTime(self): current_time = self.model_part.ProcessInfo.GetValue(KM.TIME) closest_time = next(filter(lambda x: x>current_time, self.times)) self.interface_model_part.ProcessInfo.SetValue(KM.TIME, closest_time) + def _SetDefaultTime(self): default_time = self.settings["default_time_after_interval"].GetDouble() self.interface_model_part.ProcessInfo.SetValue(KM.TIME, default_time) + def _CheckInputVariables(self): if self.settings["swap_yz_axis"].GetBool(): if self.settings["read_historical_database"].GetBool(): @@ -139,6 +146,7 @@ def _CheckInputVariables(self): KM.VariableUtils().SetNonHistoricalVariableToZero(KM.MOMENTUM_Z, self.interface_model_part.Nodes) KM.VariableUtils().SetNonHistoricalVariableToZero(KM.VELOCITY_Z, self.interface_model_part.Nodes) + def _MapToBoundaryCondition(self): for variable in self.variables: if self.settings["read_historical_database"].GetBool(): @@ -149,6 +157,7 @@ def _MapToBoundaryCondition(self): for variable in self.variables_to_fix: KM.VariableUtils().ApplyFixity(variable, True, self.model_part.Nodes) + def _CreateMapper(self): mapper_settings = KM.Parameters("""{ "mapper_type": "nearest_neighbor", @@ -172,6 +181,7 @@ def _CreateMapper(self): self.model_part, mapper_settings) + def _SmoothDefaultValue(self): elapsed_time = self.model_part.ProcessInfo.GetValue(KM.DELTA_TIME) semi_period = self.settings["semi_period_after_interval"].GetDouble() From dcf4d3645525838c1140ec96fb36ca0ac22d96d5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Tue, 9 Nov 2021 12:49:35 +0100 Subject: [PATCH 23/43] fix smooth function --- .../coupling/depth_integration_input_process.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py index 849a33c15c93..f4c3610d37a3 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py @@ -7,7 +7,7 @@ from KratosMultiphysics.HDF5Application import read_model_part_from_hdf5_process from KratosMultiphysics.HDF5Application import single_mesh_temporal_input_process from os import path, listdir -from math import exp +from numpy import expm1 def Factory(settings, model): if not isinstance(settings, KM.Parameters): @@ -185,7 +185,7 @@ def _CreateMapper(self): def _SmoothDefaultValue(self): elapsed_time = self.model_part.ProcessInfo.GetValue(KM.DELTA_TIME) semi_period = self.settings["semi_period_after_interval"].GetDouble() - smooth = exp(elapsed_time / semi_period) + smooth = -expm1(-elapsed_time / semi_period) for node in self.model_part.Nodes: for variable in self.variables: initial = node.GetSolutionStepValue(variable, 1) From 27fbd2a7422fb83c0108ccad2e25434e9afaac9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Wed, 10 Nov 2021 16:20:29 +0100 Subject: [PATCH 24/43] possible solution to gid_output_process --- kratos/python_scripts/gid_output_process.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/kratos/python_scripts/gid_output_process.py b/kratos/python_scripts/gid_output_process.py index 6db237f037a2..21cbf154b44a 100644 --- a/kratos/python_scripts/gid_output_process.py +++ b/kratos/python_scripts/gid_output_process.py @@ -29,6 +29,7 @@ class GiDOutputProcess(KM.Process): "MultiFileFlag": "SingleFile" }, "file_label": "time", + "time_label_format": "{:.12f}", "output_control_type": "step", "output_interval": 1.0, "flush_after_output": false, @@ -198,6 +199,8 @@ def ExecuteInitialize(self): msg = "{0} Error: Unknown value \"{1}\" read for parameter \"{2}\"".format(self.__class__.__name__,output_file_label,"file_label") raise Exception(msg) + self.time_label_format = result_file_configuration["time_label_format"].GetString() + output_control_type = result_file_configuration["output_control_type"].GetString() if output_control_type == "time": self.output_control_is_time = True @@ -361,7 +364,7 @@ def _InitializeGiDIO(self,gidpost_flags,param): WriteConditionsFlag.WriteConditionsOnly) # Cuts are conditions, so we always print conditions in the cut ModelPart def __get_pretty_time(self,time): - pretty_time = "{0:.12g}".format(time) + pretty_time = self.time_label_format.format(time) pretty_time = float(pretty_time) return pretty_time From cffaa3854e2a8676a4d58736e2f7b6fa90151868 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Thu, 11 Nov 2021 12:44:38 +0100 Subject: [PATCH 25/43] new intersection type --- kratos/geometries/triangle_2d_3.h | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/kratos/geometries/triangle_2d_3.h b/kratos/geometries/triangle_2d_3.h index 3b8d965727e1..60e4ffdc7b4c 100644 --- a/kratos/geometries/triangle_2d_3.h +++ b/kratos/geometries/triangle_2d_3.h @@ -26,6 +26,7 @@ #include "geometries/line_2d_2.h" #include "integration/triangle_gauss_legendre_integration_points.h" #include "integration/triangle_collocation_integration_points.h" +#include "utilities/intersection_utilities.h" namespace Kratos { @@ -497,9 +498,24 @@ template class Triangle2D3 /// detect if two triangle are intersected bool HasIntersection( const BaseType& rThisGeometry ) const override { - const BaseType& geom_1 = *this; - const BaseType& geom_2 = rThisGeometry; - return NoDivTriTriIsect(geom_1[0], geom_1[1], geom_1[2], geom_2[0], geom_2[1], geom_2[2]); + const auto geometry_type = rThisGeometry.GetGeometryType(); + + if (geometry_type == GeometryData::KratosGeometryType::Kratos_Line2D2) { + Point result; + for (auto& edge : this->GenerateEdges()) { + if (IntersectionUtilities::ComputeLineLineIntersection(rThisGeometry, edge[0], edge[1], result)) + return true; + } + return this->IsInside(rThisGeometry[0], result); + } + else if(geometry_type == GeometryData::KratosGeometryType::Kratos_Triangle2D3) { + const BaseType& geom_1 = *this; + const BaseType& geom_2 = rThisGeometry; + return NoDivTriTriIsect(geom_1[0], geom_1[1], geom_1[2], geom_2[0], geom_2[1], geom_2[2]); + } + else { + KRATOS_ERROR << "Triangle2D3::HasIntersection : Geometry cannot be identified, please, check the intersecting geometry type." << std::endl; + } } /** From 41b3fbc635cdceac38b3b6920267ac5ea872a6b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Thu, 11 Nov 2021 15:07:01 +0100 Subject: [PATCH 26/43] remove dependency between apps --- .../depth_integration_output_process.py | 36 +++++++++++++++---- 1 file changed, 29 insertions(+), 7 deletions(-) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py index de5d90aa4dad..11a45f1ac5d4 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py @@ -39,6 +39,7 @@ def __init__(self, model, settings): self.output_model_part = model.CreateModelPart(settings["output_model_part_name"].GetString()) self.store_historical = settings["store_historical_database"].GetBool() self.interval = KM.IntervalUtility(settings) + self.variables = [KM.VELOCITY, KM.MOMENTUM, SW.HEIGHT] integration_settings = KM.Parameters() integration_settings.AddValue("volume_model_part_name", settings["volume_model_part_name"]) @@ -74,15 +75,13 @@ def ExecuteInitialize(self): self._InitializeOutputModelPart() self._SetOutputProcessInfo() if not self.store_historical: - variables = [KM.VELOCITY, KM.MOMENTUM, SW.HEIGHT] - for var in variables: + for var in self.variables: KM.VariableUtils().SetNonHistoricalVariableToZero(var, self.interface_model_part.Nodes) KM.VariableUtils().SetNonHistoricalVariableToZero(var, self.output_model_part.Nodes) def ExecuteBeforeSolutionLoop(self): '''Write the interface model part in HDF5 format.''' - self._SetOutputProcessInfo() self.hdf5_process.ExecuteBeforeSolutionLoop() @@ -94,6 +93,7 @@ def ExecuteInitializeSolutionStep(self): def ExecuteBeforeOutputStep(self): '''Perform the depth integration over the interface model part.''' self.integration_process.Execute() + self._MapToOutputModelPart() def IsOutputStep(self): @@ -108,13 +108,35 @@ def PrintOutput(self): def _InitializeOutputModelPart(self): - KM.MergeVariableListsUtility().Merge(self.interface_model_part, self.output_model_part) + if self.store_historical: + self.output_model_part.AddNodalSolutionStepVariable(SW.HEIGHT) + self.output_model_part.AddNodalSolutionStepVariable(KM.MOMENTUM) + self.output_model_part.AddNodalSolutionStepVariable(KM.VELOCITY) domain_size = self.volume_model_part.ProcessInfo[KM.DOMAIN_SIZE] element_name = "Element{}D2N".format(domain_size) condition_name = "LineCondition{}D2N".format(domain_size) - KM.ConnectivityPreserveModeler().GenerateModelPart( - self.interface_model_part, self.output_model_part, element_name, condition_name) - self.output_model_part.ProcessInfo = KM.ProcessInfo() + KM.DuplicateMeshModeler(self.interface_model_part).GenerateMesh( + self.output_model_part, element_name, condition_name) + + + def _MapToOutputModelPart(self): + if self.store_historical: + for variable in self.variables: + KM.VariableUtils().CopyModelPartNodalVar( + variable, + self.interface_model_part, + self.output_model_part, + 0) + else: + for variable in self.variables: + for node_src, node_dest in zip(self.interface_model_part.Nodes, self.output_model_part.Nodes): + node_dest.SetValue(variable, node_src.GetValue(variable)) + #TODO: implement this function in VariableUtils + # KM.VariableUtils().CopyModelPartFlaggedNodalNonHistoricalVarToNonHistoricalVar( + # variable, + # self.interface_model_part, + # self.output_model_part, + # KM.Flags(), True) def _SetOutputProcessInfo(self): From 88cb7498cd9ee0f8a0ba43ee3e62273215afd7b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Thu, 11 Nov 2021 16:29:56 +0100 Subject: [PATCH 27/43] clean dependencies, rename var --- .../depth_integration_input_process.py | 95 +++++++++---------- 1 file changed, 47 insertions(+), 48 deletions(-) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py index f4c3610d37a3..52ad4229ef44 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py @@ -1,13 +1,11 @@ import KratosMultiphysics as KM import KratosMultiphysics.ShallowWaterApplication as SW import KratosMultiphysics.MappingApplication as Mapping -import KratosMultiphysics.PfemFluidDynamicsApplication as PFEM -import KratosMultiphysics.DelaunayMeshingApplication as Delaunay from KratosMultiphysics.kratos_utilities import GenerateVariableListFromInput -from KratosMultiphysics.HDF5Application import read_model_part_from_hdf5_process +from KratosMultiphysics.HDF5Application import import_model_part_from_hdf5_process from KratosMultiphysics.HDF5Application import single_mesh_temporal_input_process from os import path, listdir -from numpy import expm1 +from numpy import expm1 #TODO: implement the related function in shallow water utilities def Factory(settings, model): if not isinstance(settings, KM.Parameters): @@ -23,8 +21,8 @@ class DepthIntegrationInputProcess(KM.OutputProcess): @staticmethod def GetDefaultParameters(): return KM.Parameters("""{ - "model_part_name" : "", - "interface_model_part_name" : "interface_model_part", + "interface_model_part_name" : "", + "input_model_part_name" : "input_model_part", "read_historical_database" : false, "interval" : [0.0,"End"], "list_of_variables" : ["MOMENTUM"], @@ -37,54 +35,40 @@ def GetDefaultParameters(): }""") def __init__(self, model, settings): - """The constructor of the DepthIntegrationInputProcess""" + """The constructor of the DepthIntegrationInputProcess.""" KM.OutputProcess.__init__(self) self.settings = settings self.settings.ValidateAndAssignDefaults(self.GetDefaultParameters()) - self.model_part = model[self.settings["model_part_name"].GetString()] + self.interface_model_part = model[self.settings["interface_model_part_name"].GetString()] + self.input_model_part = model.CreateModelPart(self.settings["input_model_part_name"].GetString()) self.interval = KM.IntervalUtility(self.settings) - interface_model_part_name = self.settings["interface_model_part_name"].GetString() - if model.HasModelPart(interface_model_part_name): - raise Exception("DepthIntegrationInputProcess: There is an existing interface model part with name '{}'".format(interface_model_part_name)) - self.interface_model_part = model.CreateModelPart(interface_model_part_name) self.variables = GenerateVariableListFromInput(self.settings["list_of_variables"]) self.variables_to_fix = GenerateVariableListFromInput(self.settings["list_of_variables_to_fix"]) - hdf5_settings = KM.Parameters() - hdf5_settings.AddString("model_part_name", interface_model_part_name) - hdf5_settings.AddValue("file_settings", self.settings["file_settings"]) - data_settings = KM.Parameters("""{"list_of_variables" : ["MOMENTUM","VELOCITY","HEIGHT"]}""") - if self.settings["read_historical_database"].GetBool(): - hdf5_settings.AddValue("nodal_solution_step_data_settings", data_settings) - else: - hdf5_settings.AddValue("nodal_data_value_settings", data_settings) - hdf5_process_settings = KM.Parameters() - hdf5_process_settings.AddValue("Parameters", hdf5_settings) - - self.hdf5_read = read_model_part_from_hdf5_process.Factory(hdf5_process_settings.Clone(), model) - self.hdf5_process = single_mesh_temporal_input_process.Factory(hdf5_process_settings.Clone(), model) + self.hdf5_import = import_model_part_from_hdf5_process.Factory(self._CreateHDF5Parameters(), model) + self.hdf5_process = single_mesh_temporal_input_process.Factory(self._CreateHDF5Parameters(), model) self._GetInputTimes(self.settings['file_settings']) def Check(self): '''Check the processes.''' - self.hdf5_read.Check() + self.hdf5_import.Check() self.hdf5_process.Check() def ExecuteInitialize(self): - '''Read the interface_model_part and set the variables.''' - self.hdf5_read.ExecuteInitialize() + '''Read the input_model_part and set the variables.''' + self.hdf5_import.ExecuteInitialize() self._CheckInputVariables() self._CreateMapper() self._MapToBoundaryCondition() def ExecuteInitializeSolutionStep(self): - '''Set the variables in the interface_model_part at the current time.''' - current_time = self.model_part.ProcessInfo.GetValue(KM.TIME) + '''Set the variables in the input_model_part at the current time.''' + current_time = self.interface_model_part.ProcessInfo.GetValue(KM.TIME) if self.interval.IsInInterval(current_time): self._SetCurrentTime() self.hdf5_process.ExecuteInitializeSolutionStep() @@ -104,9 +88,10 @@ def _GetInputTimes(self, file_settings): folder_name = path.dirname(file_name) file_names = [f for f in listdir(folder_name) if path.isfile(path.join(folder_name, f))] - # Find the common parts of the found names and the file name. The difference are the times + # Find the common parts (prefix and suffix) of the found names and the file pattern + # The different part is the time, we need to store all the available times self.file_pattern = path.basename(file_name) - self.file_pattern = self.file_pattern.replace('', self.interface_model_part.Name) + self.file_pattern = self.file_pattern.replace('', self.input_model_part.Name) prefix = path.commonprefix([self.file_pattern, file_names[0]]) suffix = self.file_pattern.replace(prefix, '') suffix = path.commonprefix([''.join(reversed(suffix)), ''.join(reversed(file_names[0]))]) @@ -120,31 +105,31 @@ def _GetInputTimes(self, file_settings): def _SetCurrentTime(self): - current_time = self.model_part.ProcessInfo.GetValue(KM.TIME) + current_time = self.interface_model_part.ProcessInfo.GetValue(KM.TIME) closest_time = next(filter(lambda x: x>current_time, self.times)) - self.interface_model_part.ProcessInfo.SetValue(KM.TIME, closest_time) + self.input_model_part.ProcessInfo.SetValue(KM.TIME, closest_time) def _SetDefaultTime(self): default_time = self.settings["default_time_after_interval"].GetDouble() - self.interface_model_part.ProcessInfo.SetValue(KM.TIME, default_time) + self.input_model_part.ProcessInfo.SetValue(KM.TIME, default_time) def _CheckInputVariables(self): if self.settings["swap_yz_axis"].GetBool(): if self.settings["read_historical_database"].GetBool(): - SW.ShallowWaterUtilities().SwapYZComponents(KM.MOMENTUM, self.interface_model_part.Nodes) - SW.ShallowWaterUtilities().SwapYZComponents(KM.VELOCITY, self.interface_model_part.Nodes) + SW.ShallowWaterUtilities().SwapYZComponents(KM.MOMENTUM, self.input_model_part.Nodes) + SW.ShallowWaterUtilities().SwapYZComponents(KM.VELOCITY, self.input_model_part.Nodes) else: - SW.ShallowWaterUtilities().SwapYZComponentsNonHistorical(KM.MOMENTUM, self.interface_model_part.Nodes) - SW.ShallowWaterUtilities().SwapYZComponentsNonHistorical(KM.VELOCITY, self.interface_model_part.Nodes) + SW.ShallowWaterUtilities().SwapYZComponentsNonHistorical(KM.MOMENTUM, self.input_model_part.Nodes) + SW.ShallowWaterUtilities().SwapYZComponentsNonHistorical(KM.VELOCITY, self.input_model_part.Nodes) if self.settings["ignore_vertical_component"].GetBool(): if self.settings["read_historical_database"].GetBool(): - KM.VariableUtils().SetVariableToZero(KM.MOMENTUM_Z, self.interface_model_part.Nodes) - KM.VariableUtils().SetVariableToZero(KM.VELOCITY_Z, self.interface_model_part.Nodes) + KM.VariableUtils().SetVariableToZero(KM.MOMENTUM_Z, self.input_model_part.Nodes) + KM.VariableUtils().SetVariableToZero(KM.VELOCITY_Z, self.input_model_part.Nodes) else: - KM.VariableUtils().SetNonHistoricalVariableToZero(KM.MOMENTUM_Z, self.interface_model_part.Nodes) - KM.VariableUtils().SetNonHistoricalVariableToZero(KM.VELOCITY_Z, self.interface_model_part.Nodes) + KM.VariableUtils().SetNonHistoricalVariableToZero(KM.MOMENTUM_Z, self.input_model_part.Nodes) + KM.VariableUtils().SetNonHistoricalVariableToZero(KM.VELOCITY_Z, self.input_model_part.Nodes) def _MapToBoundaryCondition(self): @@ -155,7 +140,7 @@ def _MapToBoundaryCondition(self): self.mapper.Map(variable, variable, KM.Mapper.FROM_NON_HISTORICAL) for variable in self.variables_to_fix: - KM.VariableUtils().ApplyFixity(variable, True, self.model_part.Nodes) + KM.VariableUtils().ApplyFixity(variable, True, self.interface_model_part.Nodes) def _CreateMapper(self): @@ -168,7 +153,7 @@ def _CreateMapper(self): }""") min_point = KM.Point([ 1e6, 1e6, 1e6]) max_point = KM.Point([-1e6, -1e6, -1e6]) - for node in self.model_part.Nodes: + for node in self.interface_model_part.Nodes: for i in range(3): point = KM.Point(node) min_point[i] = min([min_point[i], point[i]]) @@ -177,16 +162,30 @@ def _CreateMapper(self): mapper_settings["search_settings"]["search_radius"].SetDouble(distance) self.mapper = KM.MapperFactory.CreateMapper( + self.input_model_part, self.interface_model_part, - self.model_part, mapper_settings) + def _CreateHDF5Parameters(self): + hdf5_settings = KM.Parameters() + hdf5_settings.AddValue("model_part_name", self.settings["input_model_part_name"]) + hdf5_settings.AddValue("file_settings", self.settings["file_settings"]) + data_settings = KM.Parameters("""{"list_of_variables" : ["MOMENTUM","VELOCITY","HEIGHT"]}""") + if self.settings["read_historical_database"].GetBool(): + hdf5_settings.AddValue("nodal_solution_step_data_settings", data_settings) + else: + hdf5_settings.AddValue("nodal_data_value_settings", data_settings) + hdf5_process_settings = KM.Parameters() + hdf5_process_settings.AddValue("Parameters", hdf5_settings) + return hdf5_process_settings + + def _SmoothDefaultValue(self): - elapsed_time = self.model_part.ProcessInfo.GetValue(KM.DELTA_TIME) + elapsed_time = self.interface_model_part.ProcessInfo.GetValue(KM.DELTA_TIME) semi_period = self.settings["semi_period_after_interval"].GetDouble() smooth = -expm1(-elapsed_time / semi_period) - for node in self.model_part.Nodes: + for node in self.interface_model_part.Nodes: for variable in self.variables: initial = node.GetSolutionStepValue(variable, 1) increment = node.GetSolutionStepValue(variable) - initial From 3cad360bae6bf376257b8beac5aebba3abde878b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Fri, 12 Nov 2021 14:45:34 +0100 Subject: [PATCH 28/43] less restrictive threshold --- .../custom_utilities/shallow_water_utilities.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/ShallowWaterApplication/custom_utilities/shallow_water_utilities.cpp b/applications/ShallowWaterApplication/custom_utilities/shallow_water_utilities.cpp index 1514d4a7fdf2..11894869899b 100644 --- a/applications/ShallowWaterApplication/custom_utilities/shallow_water_utilities.cpp +++ b/applications/ShallowWaterApplication/custom_utilities/shallow_water_utilities.cpp @@ -381,7 +381,7 @@ bool ShallowWaterUtilities::IsWet(const GeometryType& rGeometry, const double He bool ShallowWaterUtilities::IsWet(const double Height, const double DryHeight) { const double wet_fraction = PhaseFunction::WetFraction(Height, DryHeight); - const double threshold = 1.0 - 1e-16; + const double threshold = 1.0 - 1e-6; return (wet_fraction >= threshold); } From c6603e70749c6175b345ee861cbd39c8f5ec15b2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Fri, 12 Nov 2021 14:47:25 +0100 Subject: [PATCH 29/43] remove duplicated after merge --- .../custom_python/add_custom_utilities_to_python.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/applications/ShallowWaterApplication/custom_python/add_custom_utilities_to_python.cpp b/applications/ShallowWaterApplication/custom_python/add_custom_utilities_to_python.cpp index faf2a0d61299..5fa1f9cd2bc9 100644 --- a/applications/ShallowWaterApplication/custom_python/add_custom_utilities_to_python.cpp +++ b/applications/ShallowWaterApplication/custom_python/add_custom_utilities_to_python.cpp @@ -110,10 +110,6 @@ void AddCustomUtilitiesToPython(pybind11::module& m) .def("ComputeL2Norm", &ShallowWaterUtilities::ComputeL2NormAABB) .def("ComputeL2NormNonHistorical", &ShallowWaterUtilities::ComputeL2Norm) .def("ComputeL2NormNonHistorical", &ShallowWaterUtilities::ComputeL2NormAABB) - .def("SwapYZComponents", &ShallowWaterUtilities::SwapYZComponents) - .def("SwapYZComponentsNonHistorical", &ShallowWaterUtilities::SwapYZComponentsNonHistorical) - .def("SwapYZComponentsNonHistorical", &ShallowWaterUtilities::SwapYZComponentsNonHistorical) - .def("SwapYZComponentsNonHistorical", &ShallowWaterUtilities::SwapYZComponentsNonHistorical) .def("ComputeHydrostaticForces", ComputeHydrostaticForces1) .def("ComputeHydrostaticForces", ComputeHydrostaticForces2) .def("ComputeHydrostaticForces", ComputeHydrostaticForces1) From cb7dc8d858be6eb1e1bac4b8991f2f926a94c1e3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Fri, 12 Nov 2021 15:13:13 +0100 Subject: [PATCH 30/43] remove python import, update defaults --- .../depth_integration_input_process.py | 34 ++++++++++--------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py index 52ad4229ef44..b2783d6ebc90 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py @@ -5,7 +5,6 @@ from KratosMultiphysics.HDF5Application import import_model_part_from_hdf5_process from KratosMultiphysics.HDF5Application import single_mesh_temporal_input_process from os import path, listdir -from numpy import expm1 #TODO: implement the related function in shallow water utilities def Factory(settings, model): if not isinstance(settings, KM.Parameters): @@ -18,21 +17,24 @@ class DepthIntegrationInputProcess(KM.OutputProcess): Read the depth integrated values from an HDF5 file and set them as boundary conditions. """ - @staticmethod - def GetDefaultParameters(): - return KM.Parameters("""{ + def GetDefaultParameters(self): + default_parameters = KM.Parameters("""{ "interface_model_part_name" : "", "input_model_part_name" : "input_model_part", "read_historical_database" : false, "interval" : [0.0,"End"], "list_of_variables" : ["MOMENTUM"], "list_of_variables_to_fix" : ["MOMENTUM_X","MOMENTUM_Y"], - "default_time_after_interval" : 0.0, + "default_time_after_interval" : null, "semi_period_after_interval" : 1.0, "swap_yz_axis" : false, "ignore_vertical_component" : true, "file_settings" : {} }""") + if self.settings.Has("default_time_after_interval"): + if self.settings["default_time_after_interval"].IsDouble(): + default_parameters["default_time_after_interval"].SetDouble(0.0) + return default_parameters def __init__(self, model, settings): """The constructor of the DepthIntegrationInputProcess.""" @@ -75,11 +77,12 @@ def ExecuteInitializeSolutionStep(self): self._CheckInputVariables() self._MapToBoundaryCondition() else: - self._SetDefaultTime() - self.hdf5_process.ExecuteInitializeSolutionStep() - self._CheckInputVariables() - self._MapToBoundaryCondition() - self._SmoothDefaultValue() + if self.settings["default_time_after_interval"] is not None: + self._SetDefaultTime() + self.hdf5_process.ExecuteInitializeSolutionStep() + self._CheckInputVariables() + self._MapToBoundaryCondition() + self._SmoothDefaultValue() def _GetInputTimes(self, file_settings): @@ -184,9 +187,8 @@ def _CreateHDF5Parameters(self): def _SmoothDefaultValue(self): elapsed_time = self.interface_model_part.ProcessInfo.GetValue(KM.DELTA_TIME) semi_period = self.settings["semi_period_after_interval"].GetDouble() - smooth = -expm1(-elapsed_time / semi_period) - for node in self.interface_model_part.Nodes: - for variable in self.variables: - initial = node.GetSolutionStepValue(variable, 1) - increment = node.GetSolutionStepValue(variable) - initial - node.SetSolutionStepValue(variable, initial + increment * smooth) + for variable in self.variables: + SW.ShallowWaterUtilities().SmoothTemporalVariable( + self.interface_model_part, + variable, + semi_period) From 1fd9e7bc5109bf01c759319c83fbe5499f3420e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Mon, 15 Nov 2021 13:58:01 +0100 Subject: [PATCH 31/43] avoid recursive deletion --- .../custom_processes/depth_integration_process.cpp | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index 6fa8a63095ef..9af55cfdfb40 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -133,13 +133,10 @@ void DepthIntegrationProcess::InitializeIntegrationModelPart() mpIntegrationModelPart->RemoveElementsFromAllLevels(); mpIntegrationModelPart->RemoveConditionsFromAllLevels(); - // Add a dummy node in order to allow the generation of the octree - if (mrVolumeModelPart.NumberOfNodes() > 0) { - auto i_node = mrVolumeModelPart.NodesBegin(); - mpIntegrationModelPart->AddNode(&*i_node); - } else { - KRATOS_ERROR << "DepthIntegrationProcess: The volume model part is empty." << std::endl; - } + // Check the volume model part is non empty and there is almost one node in the integration model part. + // The generation of the octree needs non empty model parts. + mpIntegrationModelPart->CreateNewNode(1, 0.0, 0.0, 0.0); + KRATOS_ERROR_IF(mrVolumeModelPart.NumberOfNodes() == 0) << "DepthIntegrationProcess: The volume model part is empty." << std::endl; } void DepthIntegrationProcess::InitializeIntegrationLine() From e135657d56e8a1e08c9d1e203f798457ecc98580 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Mon, 15 Nov 2021 13:58:39 +0100 Subject: [PATCH 32/43] first steps to remove custom utility --- .../depth_integration_process.cpp | 34 ++++++++++--------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index 9af55cfdfb40..d549f2f9d1f3 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -24,7 +24,8 @@ #include "utilities/variable_utils.h" #include "utilities/parallel_utilities.h" #include "shallow_water_application_variables.h" -#include "custom_utilities/find_intersected_objects_utility.h" +#include "processes/find_intersected_geometrical_objects_process.h" +// #include "custom_utilities/find_intersected_objects_utility.h" namespace Kratos { @@ -65,23 +66,24 @@ void DepthIntegrationProcess::Execute() { double bottom, top; GetBoundingVolumeLimits(bottom, top); - // InitializeIntegrationModelPart(); - // FindIntersectedGeometricalObjectsProcess find_intersected_objects_process(*mpIntegrationModelPart, mrVolumeModelPart); - // find_intersected_objects_process.ExecuteInitialize(); - // for (auto& node : mrInterfaceModelPart.Nodes()) { - // InitializeIntegrationLine(); - // SetIntegrationLine(node, bottom, top); - // find_intersected_objects_process.FindIntersections(); - // auto intersected_objects = find_intersected_objects_process.GetIntersections(); - // Integrate(intersected_objects[0], node); - // } - FindIntersectedObjectsUtility intersections(mrVolumeModelPart); + InitializeIntegrationModelPart(); + FindIntersectedGeometricalObjectsProcess find_intersected_objects_process(*mpIntegrationModelPart, mrVolumeModelPart); + find_intersected_objects_process.ExecuteInitialize(); + InitializeIntegrationLine(); for (auto& node : mrInterfaceModelPart.Nodes()) { - GeometryType::Pointer integration_line = CreateIntegrationLine(node, bottom, top); - PointerVector intersected_objects; - intersections.FindIntersectedObjects(integration_line, intersected_objects); - Integrate(intersected_objects, node); + SetIntegrationLine(node, bottom, top); + find_intersected_objects_process.FindIntersections(); + auto intersected_objects = find_intersected_objects_process.GetIntersections(); + Integrate(intersected_objects[0], node); } + + // FindIntersectedObjectsUtility intersections(mrVolumeModelPart); + // for (auto& node : mrInterfaceModelPart.Nodes()) { + // GeometryType::Pointer integration_line = CreateIntegrationLine(node, bottom, top); + // PointerVector intersected_objects; + // intersections.FindIntersectedObjects(integration_line, intersected_objects); + // Integrate(intersected_objects, node); + // } } void DepthIntegrationProcess::Integrate(PointerVector& rObjects, NodeType& rNode) From f2397b1fe9c65a1680b99e5f02af5695a764d223 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Mon, 15 Nov 2021 14:00:21 +0100 Subject: [PATCH 33/43] tmp fix bad copy --- .../postprocess/visualization_mesh_process.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/applications/ShallowWaterApplication/python_scripts/postprocess/visualization_mesh_process.py b/applications/ShallowWaterApplication/python_scripts/postprocess/visualization_mesh_process.py index 7d3941915e1c..759918a5ddca 100644 --- a/applications/ShallowWaterApplication/python_scripts/postprocess/visualization_mesh_process.py +++ b/applications/ShallowWaterApplication/python_scripts/postprocess/visualization_mesh_process.py @@ -130,11 +130,14 @@ def _TransferVariables(self): self.topographic_model_part, 0) for variable in self.nonhistorical_variables: - KM.VariableUtils().CopyModelPartFlaggedNodalNonHistoricalVarToHistoricalVar( - variable, - self.computing_model_part, - self.topographic_model_part, - KM.Flags()) + for node_src, node_dest in zip(self.computing_model_part.Nodes, self.topographic_model_part.Nodes): + node_dest.SetValue(variable, node_src.GetValue(variable)) + #TODO: implement this function in VariableUtils + # KM.VariableUtils().CopyModelPartFlaggedNodalNonHistoricalVarToNonHistoricalVar( + # variable, + # self.computing_model_part, + # self.topographic_model_part, + # KM.Flags(), False) def _FlattenMeshCoordinates(self): From 7f8a9f78b289b6b25eb2dbe2e6d00b66c151c1f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Mon, 15 Nov 2021 15:49:41 +0100 Subject: [PATCH 34/43] move check to Check method --- .../custom_processes/depth_integration_process.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index d549f2f9d1f3..e6281189f6d1 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -135,10 +135,8 @@ void DepthIntegrationProcess::InitializeIntegrationModelPart() mpIntegrationModelPart->RemoveElementsFromAllLevels(); mpIntegrationModelPart->RemoveConditionsFromAllLevels(); - // Check the volume model part is non empty and there is almost one node in the integration model part. - // The generation of the octree needs non empty model parts. + // Ensure there is almost one node in the integration model part in order to construct the octree. mpIntegrationModelPart->CreateNewNode(1, 0.0, 0.0, 0.0); - KRATOS_ERROR_IF(mrVolumeModelPart.NumberOfNodes() == 0) << "DepthIntegrationProcess: The volume model part is empty." << std::endl; } void DepthIntegrationProcess::InitializeIntegrationLine() @@ -195,7 +193,10 @@ Geometry>::Pointer DepthIntegrationProcess::CreateIntegrationLine( int DepthIntegrationProcess::Check() { + // Check dimension of the given direction, otherwise, the product operation will be undefined. KRATOS_ERROR_IF(mDirection.size() != 3) << "DepthIntegrationProcess: The direction of integration must have three coordinates." << std::endl; + // If the volume model part is empty, then it is not possible to construct the octree. + KRATOS_ERROR_IF(mrVolumeModelPart.NumberOfNodes() == 0) << "DepthIntegrationProcess: The volume model part is empty." << std::endl; return 0; } From 156d123d49d769f9fb20b7c123daf7a54ab375e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Mon, 15 Nov 2021 22:13:27 +0100 Subject: [PATCH 35/43] wip: towards using process from the core --- .../depth_integration_process.cpp | 79 +++++++++++++++---- .../depth_integration_process.h | 6 +- 2 files changed, 69 insertions(+), 16 deletions(-) diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index e6281189f6d1..5662f2f4c9c2 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -54,12 +54,12 @@ DepthIntegrationProcess::DepthIntegrationProcess( mDirection = ThisParameters["direction_of_integration"].GetVector(); mDirection /= norm_2(mDirection); - if (rModel.HasModelPart("integration_auxiliary_model_part")) { + if (rModel.HasModelPart("integration_auxiliary_model_part")) { // TODO: remove the if and create without check mpIntegrationModelPart = &rModel.GetModelPart("integration_auxiliary_model_part"); } else { mpIntegrationModelPart = &rModel.CreateModelPart("integration_auxiliary_model_part"); } - mDimension = mrVolumeModelPart.GetProcessInfo()[DOMAIN_SIZE]; + mDimension = mrVolumeModelPart.GetProcessInfo()[DOMAIN_SIZE]; //TODO: remove and use local var } void DepthIntegrationProcess::Execute() @@ -67,15 +67,25 @@ void DepthIntegrationProcess::Execute() double bottom, top; GetBoundingVolumeLimits(bottom, top); InitializeIntegrationModelPart(); + CreateIntegrationLines(); FindIntersectedGeometricalObjectsProcess find_intersected_objects_process(*mpIntegrationModelPart, mrVolumeModelPart); find_intersected_objects_process.ExecuteInitialize(); - InitializeIntegrationLine(); - for (auto& node : mrInterfaceModelPart.Nodes()) { - SetIntegrationLine(node, bottom, top); - find_intersected_objects_process.FindIntersections(); - auto intersected_objects = find_intersected_objects_process.GetIntersections(); - Integrate(intersected_objects[0], node); - } + find_intersected_objects_process.FindIntersections(); + auto intersected_objects = find_intersected_objects_process.GetIntersections(); + Integrate(intersected_objects); + + // double bottom, top; + // GetBoundingVolumeLimits(bottom, top); + // InitializeIntegrationModelPart(); + // FindIntersectedGeometricalObjectsProcess find_intersected_objects_process(*mpIntegrationModelPart, mrVolumeModelPart); + // find_intersected_objects_process.ExecuteInitialize(); + // InitializeIntegrationLine(); + // for (auto& node : mrInterfaceModelPart.Nodes()) { + // SetIntegrationLine(node, bottom, top); + // find_intersected_objects_process.FindIntersections(); + // auto intersected_objects = find_intersected_objects_process.GetIntersections(); + // Integrate(intersected_objects[0], node); + // } // FindIntersectedObjectsUtility intersections(mrVolumeModelPart); // for (auto& node : mrInterfaceModelPart.Nodes()) { @@ -135,8 +145,47 @@ void DepthIntegrationProcess::InitializeIntegrationModelPart() mpIntegrationModelPart->RemoveElementsFromAllLevels(); mpIntegrationModelPart->RemoveConditionsFromAllLevels(); - // Ensure there is almost one node in the integration model part in order to construct the octree. - mpIntegrationModelPart->CreateNewNode(1, 0.0, 0.0, 0.0); + // // The construction of the octree needs almost one node. + // mpIntegrationModelPart->CreateNewNode(1, 0.0, 0.0, 0.0); // TODO: remove +} + +void DepthIntegrationProcess::CreateIntegrationLines(const double Low, const double High) +{ + // Set the element name + dimension = mrVolumeModelPart.GetProcessInfo()[DOMAIN_SIZE]; + std::string element_name = (dimension == 2) ? "Element2D2N" : "Element3D2N"; + + // Get some properties + ModelPart::PropertiesType::Pointer p_prop; + if (mpIntegrationModelPart->HasProperties(1)) { + p_prop = mpIntegrationModelPart->pGetProperties(1); + } else { + p_prop = mpIntegrationModelPart->CreateNewProperties(1); + } + + // Initialize the ids counter + std::size_t node_id = 0; + std::size_t elem_id = 0; + + // Create the integration lines + for (auto& node : mrInterfaceModelPart.Nodes()) { + const double distance = inner_prod(node, mDirection); + array_1d start = node + (Low + distance) * mDirection; + array_1d end = node + (High - distance) * mDirection; + mpIntegrationModelPart->CreateNewNode(++node_id, start); + mpIntegrationModelPart->CreateNewNode(++node_id, end); + mpIntegrationModelPart->CreateNewElement(element_name, ++elem_id, {node_id-1, node_id}, p_prop); + } +} + +void DepthIntegrationProcess::Integrate(std::vector>& rResults) +{ + KRATOS_ERROR_IF(rResults.size() != mrInterfaceModelPart.NumberOfElements()) << "DepthIntegrationProcess: the number of nodes in the interface and the number of integration lines mismatch." + IndexPartition(static_cast(rResults.size())).for_each([&](int i){ + auto& objects_in_line = rResults[i]; + auto i_node = mrInterfaceModelPart.NodesBegin() + i; + Integrate(objects_in_line, node); + }); } void DepthIntegrationProcess::InitializeIntegrationLine() @@ -193,10 +242,10 @@ Geometry>::Pointer DepthIntegrationProcess::CreateIntegrationLine( int DepthIntegrationProcess::Check() { - // Check dimension of the given direction, otherwise, the product operation will be undefined. - KRATOS_ERROR_IF(mDirection.size() != 3) << "DepthIntegrationProcess: The direction of integration must have three coordinates." << std::endl; - // If the volume model part is empty, then it is not possible to construct the octree. - KRATOS_ERROR_IF(mrVolumeModelPart.NumberOfNodes() == 0) << "DepthIntegrationProcess: The volume model part is empty." << std::endl; + dimension = mrVolumeModelPart.GetProcessInfo()[DOMAIN_SIZE]; + KRATOS_ERROR_IF(dimension != 2 && dimension != 3) << Info() << ": Wrong DOMAIN_SIZE equal to " << dimesion << "in model part " << mrVolumeModelPart.Name() << std::endl; + KRATOS_ERROR_IF(mDirection.size() != 3) << Info() << ": The direction of integration must have three coordinates." << std::endl; + KRATOS_ERROR_IF(mrVolumeModelPart.NumberOfNodes() == 0) << Info() << ": The volume model part is empty. Not possible to construct the octree." << std::endl; return 0; } diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h index d50189f1a8bb..8371abf847d8 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h @@ -171,7 +171,11 @@ class KRATOS_API(SHALLOW_WATER_APPLICATION) DepthIntegrationProcess : public Pro void GetBoundingVolumeLimits(double& rMin, double& rMax); - void InitializeIntegrationModelPart(); + void InitializeIntegrationModelPart(const double Low, const double High); + + void CreateIntegrationLines(); + + void Integrate(std::vector>& rResults); void InitializeIntegrationLine(); From 0cf095066672acbad8c7dd49ae19dcde3c3d94f7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Tue, 16 Nov 2021 11:48:59 +0100 Subject: [PATCH 36/43] minor c++ --- .../depth_integration_process.cpp | 26 +++++++++---------- .../depth_integration_process.h | 4 +-- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index 5662f2f4c9c2..5c13c4eab2ca 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -65,9 +65,9 @@ DepthIntegrationProcess::DepthIntegrationProcess( void DepthIntegrationProcess::Execute() { double bottom, top; - GetBoundingVolumeLimits(bottom, top); InitializeIntegrationModelPart(); - CreateIntegrationLines(); + GetBoundingVolumeLimits(bottom, top); + CreateIntegrationLines(bottom, top); FindIntersectedGeometricalObjectsProcess find_intersected_objects_process(*mpIntegrationModelPart, mrVolumeModelPart); find_intersected_objects_process.ExecuteInitialize(); find_intersected_objects_process.FindIntersections(); @@ -152,7 +152,7 @@ void DepthIntegrationProcess::InitializeIntegrationModelPart() void DepthIntegrationProcess::CreateIntegrationLines(const double Low, const double High) { // Set the element name - dimension = mrVolumeModelPart.GetProcessInfo()[DOMAIN_SIZE]; + const auto dimension = mrVolumeModelPart.GetProcessInfo()[DOMAIN_SIZE]; std::string element_name = (dimension == 2) ? "Element2D2N" : "Element3D2N"; // Get some properties @@ -170,21 +170,21 @@ void DepthIntegrationProcess::CreateIntegrationLines(const double Low, const dou // Create the integration lines for (auto& node : mrInterfaceModelPart.Nodes()) { const double distance = inner_prod(node, mDirection); - array_1d start = node + (Low + distance) * mDirection; - array_1d end = node + (High - distance) * mDirection; - mpIntegrationModelPart->CreateNewNode(++node_id, start); - mpIntegrationModelPart->CreateNewNode(++node_id, end); - mpIntegrationModelPart->CreateNewElement(element_name, ++elem_id, {node_id-1, node_id}, p_prop); + array_1d start = node + mDirection * (Low + distance); + array_1d end = node + mDirection * (High - distance); + mpIntegrationModelPart->CreateNewNode(++node_id, start[0], start[1], start[2]); + mpIntegrationModelPart->CreateNewNode(++node_id, end[0], end[1], end[2]); + mpIntegrationModelPart->CreateNewElement(element_name, ++elem_id, {{node_id-1, node_id}}, p_prop); } } void DepthIntegrationProcess::Integrate(std::vector>& rResults) { - KRATOS_ERROR_IF(rResults.size() != mrInterfaceModelPart.NumberOfElements()) << "DepthIntegrationProcess: the number of nodes in the interface and the number of integration lines mismatch." + KRATOS_ERROR_IF(rResults.size() != mrInterfaceModelPart.NumberOfElements()) << "DepthIntegrationProcess: the number of nodes in the interface and the number of integration lines mismatch."; IndexPartition(static_cast(rResults.size())).for_each([&](int i){ auto& objects_in_line = rResults[i]; auto i_node = mrInterfaceModelPart.NodesBegin() + i; - Integrate(objects_in_line, node); + Integrate(objects_in_line, *i_node); }); } @@ -242,9 +242,9 @@ Geometry>::Pointer DepthIntegrationProcess::CreateIntegrationLine( int DepthIntegrationProcess::Check() { - dimension = mrVolumeModelPart.GetProcessInfo()[DOMAIN_SIZE]; - KRATOS_ERROR_IF(dimension != 2 && dimension != 3) << Info() << ": Wrong DOMAIN_SIZE equal to " << dimesion << "in model part " << mrVolumeModelPart.Name() << std::endl; - KRATOS_ERROR_IF(mDirection.size() != 3) << Info() << ": The direction of integration must have three coordinates." << std::endl; + const auto dimension = mrVolumeModelPart.GetProcessInfo()[DOMAIN_SIZE]; + KRATOS_ERROR_IF(dimension != 2 && dimension != 3) << Info() << ": Wrong DOMAIN_SIZE equal to " << dimension << "in model part " << mrVolumeModelPart.Name() << std::endl; + KRATOS_ERROR_IF(mDirection.size() != 3) << Info() << ": The direction of integration must be given with three coordinates." << std::endl; KRATOS_ERROR_IF(mrVolumeModelPart.NumberOfNodes() == 0) << Info() << ": The volume model part is empty. Not possible to construct the octree." << std::endl; return 0; } diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h index 8371abf847d8..011431be5c70 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h @@ -171,9 +171,9 @@ class KRATOS_API(SHALLOW_WATER_APPLICATION) DepthIntegrationProcess : public Pro void GetBoundingVolumeLimits(double& rMin, double& rMax); - void InitializeIntegrationModelPart(const double Low, const double High); + void InitializeIntegrationModelPart(); - void CreateIntegrationLines(); + void CreateIntegrationLines(const double Low, const double High); void Integrate(std::vector>& rResults); From 2f1b3129037528c17312b92b0270049644721e57 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Tue, 16 Nov 2021 12:51:22 +0100 Subject: [PATCH 37/43] fix process --- .../depth_integration_process.cpp | 46 +++++++++++-------- 1 file changed, 26 insertions(+), 20 deletions(-) diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index 5c13c4eab2ca..c6575fc7b342 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -99,27 +99,33 @@ void DepthIntegrationProcess::Execute() void DepthIntegrationProcess::Integrate(PointerVector& rObjects, NodeType& rNode) { array_1d velocity = ZeroVector(3); - double min_elevation = 1e6; - double max_elevation = -1e6; - int num_nodes = 0; - for (auto& object : rObjects) { - array_1d obj_velocity = ZeroVector(3); - double obj_min_elevation = 1e6; - double obj_max_elevation = -1e6; - int obj_num_nodes = object.GetGeometry().size(); - for (auto& node : object.GetGeometry()) { - velocity += node.FastGetSolutionStepValue(VELOCITY); - obj_min_elevation = std::min(obj_min_elevation, inner_prod(mDirection, node)); - obj_max_elevation = std::max(obj_max_elevation, inner_prod(mDirection, node)); + array_1d momentum = ZeroVector(3); + double height = 0.0; + + if (rObjects.size() > 0) + { + double min_elevation = 1e6; + double max_elevation = -1e6; + int num_nodes = 0; + for (auto& object : rObjects) { + array_1d obj_velocity = ZeroVector(3); + double obj_min_elevation = 1e6; + double obj_max_elevation = -1e6; + int obj_num_nodes = object.GetGeometry().size(); + for (auto& node : object.GetGeometry()) { + velocity += node.FastGetSolutionStepValue(VELOCITY); + obj_min_elevation = std::min(obj_min_elevation, inner_prod(mDirection, node)); + obj_max_elevation = std::max(obj_max_elevation, inner_prod(mDirection, node)); + } + velocity += obj_velocity; + min_elevation = std::min(min_elevation, obj_min_elevation); + max_elevation = std::max(max_elevation, obj_max_elevation); + num_nodes += obj_num_nodes; } - velocity += obj_velocity; - min_elevation = std::min(min_elevation, obj_min_elevation); - max_elevation = std::max(max_elevation, obj_max_elevation); - num_nodes += obj_num_nodes; + velocity /= num_nodes; + height = max_elevation - min_elevation; + momentum = height*velocity; } - velocity /= num_nodes; - double height = max_elevation - min_elevation; - const array_1d momentum = height*velocity; SetValue(rNode, MOMENTUM, momentum); SetValue(rNode, VELOCITY, velocity); SetValue(rNode, HEIGHT, height); @@ -180,7 +186,7 @@ void DepthIntegrationProcess::CreateIntegrationLines(const double Low, const dou void DepthIntegrationProcess::Integrate(std::vector>& rResults) { - KRATOS_ERROR_IF(rResults.size() != mrInterfaceModelPart.NumberOfElements()) << "DepthIntegrationProcess: the number of nodes in the interface and the number of integration lines mismatch."; + KRATOS_ERROR_IF(rResults.size() != mrInterfaceModelPart.NumberOfNodes()) << "DepthIntegrationProcess: the number of nodes in the interface and the number of integration lines mismatch."; IndexPartition(static_cast(rResults.size())).for_each([&](int i){ auto& objects_in_line = rResults[i]; auto i_node = mrInterfaceModelPart.NodesBegin() + i; From 7ef1060c3f74d9e09d14415d697464591180d143 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Tue, 16 Nov 2021 15:03:29 +0100 Subject: [PATCH 38/43] add c++ test --- .../depth_integration_process.cpp | 2 +- .../test_depth_integration_process.cpp | 168 ++++++++++++++++++ 2 files changed, 169 insertions(+), 1 deletion(-) create mode 100644 applications/ShallowWaterApplication/tests/cpp_tests/test_depth_integration_process.cpp diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index c6575fc7b342..1fe2d0a2a854 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -187,7 +187,7 @@ void DepthIntegrationProcess::CreateIntegrationLines(const double Low, const dou void DepthIntegrationProcess::Integrate(std::vector>& rResults) { KRATOS_ERROR_IF(rResults.size() != mrInterfaceModelPart.NumberOfNodes()) << "DepthIntegrationProcess: the number of nodes in the interface and the number of integration lines mismatch."; - IndexPartition(static_cast(rResults.size())).for_each([&](int i){ + IndexPartition(static_cast(rResults.size())).for_each([&](int i) { auto& objects_in_line = rResults[i]; auto i_node = mrInterfaceModelPart.NodesBegin() + i; Integrate(objects_in_line, *i_node); diff --git a/applications/ShallowWaterApplication/tests/cpp_tests/test_depth_integration_process.cpp b/applications/ShallowWaterApplication/tests/cpp_tests/test_depth_integration_process.cpp new file mode 100644 index 000000000000..01d974d8d2a7 --- /dev/null +++ b/applications/ShallowWaterApplication/tests/cpp_tests/test_depth_integration_process.cpp @@ -0,0 +1,168 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Miguel Maso Sotomayor +// +// + +// System includes + +// External includes + +// Project includes +#include "testing/testing.h" +#include "containers/model.h" +#include "utilities/parallel_utilities.h" +#include "geometries/hexahedra_3d_8.h" +#include "geometries/quadrilateral_2d_4.h" +#include "processes/structured_mesh_generator_process.h" +#include "custom_processes/depth_integration_process.h" +#include "shallow_water_application_variables.h" + +namespace Kratos { + +namespace Testing { + +typedef ModelPart::NodeType NodeType; + +void FillModelParts2D(ModelPart& rVolumeModelPart, ModelPart& rInterfaceModelPart) +{ + Node<3>::Pointer p_point_1 = Kratos::make_intrusive>(1, 0.0, 0.0, 0.0); + Node<3>::Pointer p_point_2 = Kratos::make_intrusive>(2, 0.0, 1.0, 0.0); + Node<3>::Pointer p_point_3 = Kratos::make_intrusive>(3, 1.0, 1.0, 0.0); + Node<3>::Pointer p_point_4 = Kratos::make_intrusive>(4, 1.0, 0.0, 0.0); + + Quadrilateral2D4> geometry(p_point_1, p_point_2, p_point_3, p_point_4); + + Parameters mesher_parameters(R"( + { + "number_of_divisions" : 5, + "element_name" : "Element2D3N", + "create_skin_sub_model_part" : false + })"); + + auto& root_model_part = rVolumeModelPart.GetRootModelPart(); + root_model_part.AddNodalSolutionStepVariable(VELOCITY); + root_model_part.AddNodalSolutionStepVariable(MOMENTUM); + root_model_part.AddNodalSolutionStepVariable(HEIGHT); + root_model_part.GetProcessInfo().SetValue(DOMAIN_SIZE, 2); + StructuredMeshGeneratorProcess(geometry, rVolumeModelPart, mesher_parameters).Execute(); + + auto id = root_model_part.NumberOfNodes(); + rInterfaceModelPart.CreateNewNode(++id, 0.5, 0.0, 0.0); +} + +void FillModelParts3D(ModelPart& rVolumeModelPart, ModelPart& rInterfaceModelPart) +{ + Node<3>::Pointer p_point_1 = Kratos::make_intrusive>(1, 0.0, 0.0, 0.0); + Node<3>::Pointer p_point_2 = Kratos::make_intrusive>(2, 1.0, 0.0, 0.0); + Node<3>::Pointer p_point_3 = Kratos::make_intrusive>(3, 1.0, 1.0, 0.0); + Node<3>::Pointer p_point_4 = Kratos::make_intrusive>(4, 0.0, 1.0, 0.0); + Node<3>::Pointer p_point_5 = Kratos::make_intrusive>(5, 0.0, 0.0, 1.0); + Node<3>::Pointer p_point_6 = Kratos::make_intrusive>(6, 1.0, 0.0, 1.0); + Node<3>::Pointer p_point_7 = Kratos::make_intrusive>(7, 1.0, 1.0, 1.0); + Node<3>::Pointer p_point_8 = Kratos::make_intrusive>(8, 0.0, 1.0, 1.0); + + Hexahedra3D8> geometry( + p_point_1, p_point_2, p_point_3, p_point_4, p_point_5, p_point_6, p_point_7, p_point_8); + + Parameters mesher_parameters(R"( + { + "number_of_divisions" : 3, + "element_name" : "Element3D4N", + "create_skin_sub_model_part" : false + })"); + + auto& root_model_part = rVolumeModelPart.GetRootModelPart(); + root_model_part.AddNodalSolutionStepVariable(VELOCITY); + root_model_part.AddNodalSolutionStepVariable(MOMENTUM); + root_model_part.AddNodalSolutionStepVariable(HEIGHT); + root_model_part.GetProcessInfo().SetValue(DOMAIN_SIZE, 3); + StructuredMeshGeneratorProcess(geometry, rVolumeModelPart, mesher_parameters).Execute(); + + auto id = root_model_part.NumberOfNodes(); + rInterfaceModelPart.CreateNewNode(++id, 0.5, 0.00, 0.0); + rInterfaceModelPart.CreateNewNode(++id, 0.5, 0.25, 0.0); + rInterfaceModelPart.CreateNewNode(++id, 0.5, 0.50, 0.0); + rInterfaceModelPart.CreateNewNode(++id, 0.5, 0.75, 0.0); + rInterfaceModelPart.CreateNewNode(++id, 0.5, 1.00, 0.0); +} + +void ApplyVelocityField(ModelPart& rModelPart) +{ + block_for_each(rModelPart.Nodes(), [](NodeType& rNode){ + array_1d& velocity = rNode.FastGetSolutionStepValue(VELOCITY); + const array_1d& coords = rNode.Coordinates(); + velocity[0] = coords[1] + coords[2] * coords[2]; + velocity[1] = 0.0; + velocity[2] = 0.0; + }); +} + +KRATOS_TEST_CASE_IN_SUITE(DepthIntegrationProcess2D, ShallowWaterApplicationFastSuite) +{ + Model model; + auto& r_model_part = model.CreateModelPart("model_part"); + auto& r_volume_model_part = r_model_part.CreateSubModelPart("volume"); + auto& r_interface_model_part = r_model_part.CreateSubModelPart("interface"); + FillModelParts2D(r_volume_model_part, r_interface_model_part); + ApplyVelocityField(r_volume_model_part); + + Parameters process_parameters(R"( + { + "volume_model_part_name" : "model_part.volume", + "interface_model_part_name" : "model_part.interface", + "direction_of_integration" : [0.0, 1.0, 0.0], + "store_historical_database" : false + })"); + DepthIntegrationProcess(model, process_parameters).Execute(); + + std::vector> reference; + reference.push_back({0.482051, 0.0, 0.0}); + + for(std::size_t i = 0; i < r_interface_model_part.NumberOfNodes(); ++i) { + auto i_node = r_interface_model_part.NodesBegin() + i; + KRATOS_CHECK_VECTOR_NEAR(i_node->GetValue(VELOCITY), reference[i], 1e-6); + } +} + +KRATOS_TEST_CASE_IN_SUITE(DepthIntegrationProcess3D, ShallowWaterApplicationFastSuite) +{ + Model model; + auto& r_model_part = model.CreateModelPart("model_part"); + auto& r_volume_model_part = r_model_part.CreateSubModelPart("volume"); + auto& r_interface_model_part = r_model_part.CreateSubModelPart("interface"); + FillModelParts3D(r_volume_model_part, r_interface_model_part); + ApplyVelocityField(r_volume_model_part); + + Parameters process_parameters(R"( + { + "volume_model_part_name" : "model_part.volume", + "interface_model_part_name" : "model_part.interface", + "direction_of_integration" : [0.0, 0.0, 1.0], + "store_historical_database" : false + })"); + DepthIntegrationProcess(model, process_parameters).Execute(); + + std::vector> reference; + reference.push_back({0.518519, 0.0, 0.0}); + reference.push_back({0.759259, 0.0, 0.0}); + reference.push_back({0.851852, 0.0, 0.0}); + reference.push_back({1.092592, 0.0, 0.0}); + reference.push_back({1.240740, 0.0, 0.0}); + + for(std::size_t i = 0; i < r_interface_model_part.NumberOfNodes(); ++i) { + auto i_node = r_interface_model_part.NodesBegin() + i; + KRATOS_CHECK_VECTOR_NEAR(i_node->GetValue(VELOCITY), reference[i], 1e-6); + } +} + +} // namespace Testing + +} // namespace Kratos From 17b83c1d7c13823c66c59e330bce047a95c4660c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Tue, 16 Nov 2021 16:31:49 +0100 Subject: [PATCH 39/43] towards removing c++ utility --- .../depth_integration_process.cpp | 162 +++++------------- .../depth_integration_process.h | 19 +- 2 files changed, 47 insertions(+), 134 deletions(-) diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index 1fe2d0a2a854..7be0e21f25ec 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -53,13 +53,7 @@ DepthIntegrationProcess::DepthIntegrationProcess( mStoreHistorical = ThisParameters["store_historical_database"].GetBool(); mDirection = ThisParameters["direction_of_integration"].GetVector(); mDirection /= norm_2(mDirection); - - if (rModel.HasModelPart("integration_auxiliary_model_part")) { // TODO: remove the if and create without check - mpIntegrationModelPart = &rModel.GetModelPart("integration_auxiliary_model_part"); - } else { - mpIntegrationModelPart = &rModel.CreateModelPart("integration_auxiliary_model_part"); - } - mDimension = mrVolumeModelPart.GetProcessInfo()[DOMAIN_SIZE]; //TODO: remove and use local var + mpIntegrationModelPart = &rModel.CreateModelPart("integration_auxiliary_model_part"); } void DepthIntegrationProcess::Execute() @@ -73,62 +67,17 @@ void DepthIntegrationProcess::Execute() find_intersected_objects_process.FindIntersections(); auto intersected_objects = find_intersected_objects_process.GetIntersections(); Integrate(intersected_objects); - - // double bottom, top; - // GetBoundingVolumeLimits(bottom, top); - // InitializeIntegrationModelPart(); - // FindIntersectedGeometricalObjectsProcess find_intersected_objects_process(*mpIntegrationModelPart, mrVolumeModelPart); - // find_intersected_objects_process.ExecuteInitialize(); - // InitializeIntegrationLine(); - // for (auto& node : mrInterfaceModelPart.Nodes()) { - // SetIntegrationLine(node, bottom, top); - // find_intersected_objects_process.FindIntersections(); - // auto intersected_objects = find_intersected_objects_process.GetIntersections(); - // Integrate(intersected_objects[0], node); - // } - - // FindIntersectedObjectsUtility intersections(mrVolumeModelPart); - // for (auto& node : mrInterfaceModelPart.Nodes()) { - // GeometryType::Pointer integration_line = CreateIntegrationLine(node, bottom, top); - // PointerVector intersected_objects; - // intersections.FindIntersectedObjects(integration_line, intersected_objects); - // Integrate(intersected_objects, node); - // } } -void DepthIntegrationProcess::Integrate(PointerVector& rObjects, NodeType& rNode) +void DepthIntegrationProcess::InitializeIntegrationModelPart() { - array_1d velocity = ZeroVector(3); - array_1d momentum = ZeroVector(3); - double height = 0.0; - - if (rObjects.size() > 0) - { - double min_elevation = 1e6; - double max_elevation = -1e6; - int num_nodes = 0; - for (auto& object : rObjects) { - array_1d obj_velocity = ZeroVector(3); - double obj_min_elevation = 1e6; - double obj_max_elevation = -1e6; - int obj_num_nodes = object.GetGeometry().size(); - for (auto& node : object.GetGeometry()) { - velocity += node.FastGetSolutionStepValue(VELOCITY); - obj_min_elevation = std::min(obj_min_elevation, inner_prod(mDirection, node)); - obj_max_elevation = std::max(obj_max_elevation, inner_prod(mDirection, node)); - } - velocity += obj_velocity; - min_elevation = std::min(min_elevation, obj_min_elevation); - max_elevation = std::max(max_elevation, obj_max_elevation); - num_nodes += obj_num_nodes; - } - velocity /= num_nodes; - height = max_elevation - min_elevation; - momentum = height*velocity; - } - SetValue(rNode, MOMENTUM, momentum); - SetValue(rNode, VELOCITY, velocity); - SetValue(rNode, HEIGHT, height); + // Empty the model part + VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Nodes()); + VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Elements()); + VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Conditions()); + mpIntegrationModelPart->RemoveNodesFromAllLevels(); + mpIntegrationModelPart->RemoveElementsFromAllLevels(); + mpIntegrationModelPart->RemoveConditionsFromAllLevels(); } void DepthIntegrationProcess::GetBoundingVolumeLimits(double& rMin, double& rMax) @@ -141,20 +90,6 @@ void DepthIntegrationProcess::GetBoundingVolumeLimits(double& rMin, double& rMax }); } -void DepthIntegrationProcess::InitializeIntegrationModelPart() -{ - // Empty the model part - VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Nodes()); - VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Elements()); - VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Conditions()); - mpIntegrationModelPart->RemoveNodesFromAllLevels(); - mpIntegrationModelPart->RemoveElementsFromAllLevels(); - mpIntegrationModelPart->RemoveConditionsFromAllLevels(); - - // // The construction of the octree needs almost one node. - // mpIntegrationModelPart->CreateNewNode(1, 0.0, 0.0, 0.0); // TODO: remove -} - void DepthIntegrationProcess::CreateIntegrationLines(const double Low, const double High) { // Set the element name @@ -194,56 +129,39 @@ void DepthIntegrationProcess::Integrate(std::vector& rObjects, NodeType& rNode) { - // Remove the dummy node - VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Nodes()); - mpIntegrationModelPart->RemoveNodesFromAllLevels(); - - // Set the element name - std::string element_name = (mDimension == 2) ? "Element2D2N" : "Element3D2N"; - - // Create the integration lines - ModelPart::PropertiesType::Pointer p_prop; - if (mpIntegrationModelPart->HasProperties(1)) { - p_prop = mpIntegrationModelPart->pGetProperties(1); - } else { - p_prop = mpIntegrationModelPart->CreateNewProperties(1); + array_1d velocity = ZeroVector(3); + array_1d momentum = ZeroVector(3); + double height = 0.0; + + if (rObjects.size() > 0) + { + double min_elevation = 1e6; + double max_elevation = -1e6; + int num_nodes = 0; + for (auto& object : rObjects) { + array_1d obj_velocity = ZeroVector(3); + double obj_min_elevation = 1e6; + double obj_max_elevation = -1e6; + int obj_num_nodes = object.GetGeometry().size(); + for (auto& node : object.GetGeometry()) { + velocity += node.FastGetSolutionStepValue(VELOCITY); + obj_min_elevation = std::min(obj_min_elevation, inner_prod(mDirection, node)); + obj_max_elevation = std::max(obj_max_elevation, inner_prod(mDirection, node)); + } + velocity += obj_velocity; + min_elevation = std::min(min_elevation, obj_min_elevation); + max_elevation = std::max(max_elevation, obj_max_elevation); + num_nodes += obj_num_nodes; + } + velocity /= num_nodes; + height = max_elevation - min_elevation; + momentum = height*velocity; } - mpIntegrationModelPart->CreateNewNode(1, 0.0, 0.0, 0.0); - mpIntegrationModelPart->CreateNewNode(2, 0.0, 0.0, 1.0); - mpIntegrationModelPart->CreateNewElement(element_name, 1, {{1, 2}}, p_prop); -} - -void DepthIntegrationProcess::SetIntegrationLine( - const NodeType& rNode, - const double Bottom, - const double Top) -{ - // Set the integration limits to the lines - array_1d base_point = rNode - mDirection * inner_prod(rNode, mDirection); - array_1d start = base_point + Top * mDirection; - array_1d end = base_point + Bottom * mDirection; - auto& integration_line = mpIntegrationModelPart->GetElement(1); - integration_line.GetGeometry()[0].Coordinates() = start; - integration_line.GetGeometry()[1].Coordinates() = end; -} - -Geometry>::Pointer DepthIntegrationProcess::CreateIntegrationLine( - const NodeType& rNode, - const double Bottom, - const double Top) -{ - array_1d origin = rNode - mDirection * inner_prod(rNode, mDirection); - array_1d start = origin + Top * mDirection; - array_1d end = origin + Bottom * mDirection; - PointerVector> nodes; - nodes.push_back(NodeType::Pointer(new NodeType(0, start))); - nodes.push_back(NodeType::Pointer(new NodeType(0, end))); - if (mDimension == 2) - return Kratos::make_shared>(nodes); - else - return Kratos::make_shared>(nodes); + SetValue(rNode, MOMENTUM, momentum); + SetValue(rNode, VELOCITY, velocity); + SetValue(rNode, HEIGHT, height); } int DepthIntegrationProcess::Check() diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h index 011431be5c70..4876c7b66015 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h @@ -78,17 +78,19 @@ class KRATOS_API(SHALLOW_WATER_APPLICATION) DepthIntegrationProcess : public Pro ///@{ /** - * @brief Default constructor. + * @brief Default constructor * @details Removed */ DepthIntegrationProcess() = delete; /** - * Constructor with Model and Parameters + * @brief Constructor with Model and Parameters */ DepthIntegrationProcess(Model& rModel, Parameters ThisParameters = Parameters()); - /// Destructor. + /** + * @brief Destructor + */ ~DepthIntegrationProcess() override = default; ///@} @@ -151,7 +153,6 @@ class KRATOS_API(SHALLOW_WATER_APPLICATION) DepthIntegrationProcess : public Pro ModelPart* mpIntegrationModelPart; array_1d mDirection; bool mStoreHistorical; - int mDimension; ///@} ///@name Member Variables @@ -167,21 +168,15 @@ class KRATOS_API(SHALLOW_WATER_APPLICATION) DepthIntegrationProcess : public Pro ///@name Private Operations ///@{ - void Integrate(PointerVector& rObjects, NodeType& rNode); + void InitializeIntegrationModelPart(); void GetBoundingVolumeLimits(double& rMin, double& rMax); - void InitializeIntegrationModelPart(); - void CreateIntegrationLines(const double Low, const double High); void Integrate(std::vector>& rResults); - void InitializeIntegrationLine(); - - void SetIntegrationLine(const NodeType& rNode, const double Bottom, const double Top); - - GeometryType::Pointer CreateIntegrationLine(const NodeType& rNode, const double Bottom, const double Top); + void Integrate(PointerVector& rObjects, NodeType& rNode); template> void SetValue(NodeType& rNode, const TVarType& rVariable, TDataType& rValue) From c2255d3875de81757ea5b94bf2e5948ed54c8175 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Tue, 16 Nov 2021 16:34:49 +0100 Subject: [PATCH 40/43] remove file --- .../depth_integration_process.cpp | 1 - .../find_intersected_objects_utility.cpp | 158 ------------ .../find_intersected_objects_utility.h | 229 ------------------ 3 files changed, 388 deletions(-) delete mode 100644 applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.cpp delete mode 100644 applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.h diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index 7be0e21f25ec..0e83cacd7b04 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -25,7 +25,6 @@ #include "utilities/parallel_utilities.h" #include "shallow_water_application_variables.h" #include "processes/find_intersected_geometrical_objects_process.h" -// #include "custom_utilities/find_intersected_objects_utility.h" namespace Kratos { diff --git a/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.cpp b/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.cpp deleted file mode 100644 index 97b7e56da106..000000000000 --- a/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.cpp +++ /dev/null @@ -1,158 +0,0 @@ -// | / | -// ' / __| _` | __| _ \ __| -// . \ | ( | | ( |\__ ` -// _|\_\_| \__,_|\__|\___/ ____/ -// Multi-Physics -// -// License: BSD License -// Kratos default license: kratos/license.txt -// -// Main authors: Miguel Maso Sotomayor -// - - -// System includes - - -// External includes - - -// Project includes -#include "includes/model_part.h" -#include "utilities/parallel_utilities.h" -#include "find_intersected_objects_utility.h" - -namespace Kratos -{ - -/// Local Flags -KRATOS_CREATE_LOCAL_FLAG(FindIntersectedObjectsUtility, INTERSECT_ELEMENTS, 0); -KRATOS_CREATE_LOCAL_FLAG(FindIntersectedObjectsUtility, INTERSECT_CONDITIONS, 1); - -FindIntersectedObjectsUtility::FindIntersectedObjectsUtility( - ModelPart& rThisModelPart, - Parameters ThisParameters -) : mrModelPart(rThisModelPart), - mpOctree(new OctreeType()) -{ - Parameters default_parameters(R"({ - "intersect_elements" : true, - "intersect_conditions" : false - })"); - ThisParameters.ValidateAndAssignDefaults(default_parameters); - - // Set the flags - const bool intersect_elements = ThisParameters["intersect_elements"].GetBool(); - const bool intersect_conditions = ThisParameters["intersect_conditions"].GetBool(); - - mOptions.Set(INTERSECT_ELEMENTS, intersect_elements); - mOptions.Set(INTERSECT_CONDITIONS, intersect_conditions); - - if (mrModelPart.NumberOfNodes() > 0) { - GenerateOctree(); - } -} - -void FindIntersectedObjectsUtility::UpdateSearchStructure() -{ - if (mrModelPart.NumberOfNodes() > 0) { - GenerateOctree(); - } -} - -void FindIntersectedObjectsUtility::FindIntersectedObjects( - GeometryType::Pointer pGeometry, - PointerVector& rResults) -{ - OctreeCellVectorType leaves; - leaves.clear(); - auto p_geometrical_object = Kratos::make_intrusive(0, pGeometry); - mpOctree->GetIntersectedLeaves(p_geometrical_object, leaves); - FindIntersectedObjects(*pGeometry, leaves, rResults); -} - -void FindIntersectedObjectsUtility::FindIntersectedObjects( - const GeometryType& rIntersectionGeometry, - OctreeCellVectorType& rLeaves, - PointerVector& rResults) -{ - for (auto p_leaf : rLeaves) { - for (auto p_volume_object : *(p_leaf->pGetObjects())) { - if (p_volume_object->GetGeometry().HasIntersection(rIntersectionGeometry)) { - rResults.push_back(p_volume_object); - } - } - } -} - -void FindIntersectedObjectsUtility::GenerateOctree() -{ - SetOctreeBoundingBox(); - - // Adding the volume model part to the octree - for (auto it_node = mrModelPart.NodesBegin(); it_node != mrModelPart.NodesEnd(); it_node++) { -#ifdef KRATOS_USE_AMATRIX // This macro definition is for the migration period and to be removed afterward please do not use it - mpOctree->Insert(it_node->Coordinates().data()); -#else - mpOctree->Insert(it_node->Coordinates().data().data()); -#endif // ifdef KRATOS_USE_AMATRIX - } - - // Add elements - if (mOptions.Is(INTERSECT_ELEMENTS)) { - auto& r_elements_array = mrModelPart.Elements(); - const auto it_elements_begin = r_elements_array.begin(); - const auto number_of_elements = r_elements_array.size(); - - // Iterate over the elements - for (int i = 0; i < static_cast(number_of_elements); i++) { - auto it_elements = it_elements_begin + i; - mpOctree->Insert(*(it_elements).base()); - } - } - - // Add conditions - if (mOptions.Is(INTERSECT_CONDITIONS)) { - auto& r_conditions_array = mrModelPart.Conditions(); - const auto it_conditions_begin = r_conditions_array.begin(); - const auto number_of_conditions = r_conditions_array.size(); - - // Iterate over the conditions - for (int i = 0; i < static_cast(number_of_conditions); i++) { - auto it_conditions = it_conditions_begin + i; - mpOctree->Insert(*(it_conditions).base()); - } - } -} - -void FindIntersectedObjectsUtility::SetOctreeBoundingBox() -{ - PointType low(mrModelPart.NodesBegin()->Coordinates()); - PointType high(mrModelPart.NodesBegin()->Coordinates()); - - // Loop over all nodes in the volume model part - for (auto it_node = mrModelPart.NodesBegin(); it_node != mrModelPart.NodesEnd(); it_node++) { - const array_1d& r_coordinates = it_node->Coordinates(); - for (IndexType i = 0; i < 3; i++) { - low[i] = r_coordinates[i] < low[i] ? r_coordinates[i] : low[i]; - high[i] = r_coordinates[i] > high[i] ? r_coordinates[i] : high[i]; - } - } - - // Slightly increase the bounding box size to avoid problems with geometries in the borders - // Note that std::numeric_limits::double() is added for the 2D cases. Otherwise, the - // third component will be 0, breaking the octree behaviour. - for(IndexType i = 0 ; i < 3; i++) { - low[i] -= std::abs(high[i] - low[i])*1e-3 + std::numeric_limits::epsilon(); - high[i] += std::abs(high[i] - low[i])*1e-3 + std::numeric_limits::epsilon(); - } - - // TODO: Octree needs refactoring to work with BoundingBox. Pooyan. -#ifdef KRATOS_USE_AMATRIX // This macro definition is for the migration period and to be removed afterward please do not use it - mpOctree->SetBoundingBox(low.data(), high.data()); -#else - mpOctree->SetBoundingBox(low.data().data(), high.data().data()); -#endif // ifdef KRATOS_USE_AMATRIX -} - -} // namespace Kratos. diff --git a/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.h b/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.h deleted file mode 100644 index e6382e4e6aab..000000000000 --- a/applications/ShallowWaterApplication/custom_utilities/find_intersected_objects_utility.h +++ /dev/null @@ -1,229 +0,0 @@ -// | / | -// ' / __| _` | __| _ \ __| -// . \ | ( | | ( |\__ ` -// _|\_\_| \__,_|\__|\___/ ____/ -// Multi-Physics -// -// License: BSD License -// Kratos default license: kratos/license.txt -// -// Main authors: Miguel Maso Sotomayor -// - - -#ifndef KRATOS_FIND_INTERSECTED_ELEMENTS_UTILITY_H_INCLUDED -#define KRATOS_FIND_INTERSECTED_ELEMENTS_UTILITY_H_INCLUDED - - -// System includes - - -// External includes - - -// Project includes -#include "includes/kratos_parameters.h" -#include "spatial_containers/octree_binary.h" -#include "processes/find_intersected_geometrical_objects_process.h" - -namespace Kratos -{ -///@name Kratos Globals -///@{ - -///@} -///@name Type Definitions -///@{ - -///@} -///@name Enum's -///@{ - -///@} -///@name Functions -///@{ - -///@} -///@name Kratos Classes -///@{ - -/** - * @brief Forward declaration of ModelPart - */ -class ModelPart; - -/** - * @brief Find the objects in a volume model part that are intersected ty the given geometry - * @author Miguel Maso Sotomayor - * @ingroup ShallowWaterApplication -*/ -class KRATOS_API(SHALLOW_WATER_APPLICATION) FindIntersectedObjectsUtility -{ -public: - ///@name Type Definitions - ///@{ - - /// Pointer definition of FindIntersectedObjectsUtility - KRATOS_CLASS_POINTER_DEFINITION(FindIntersectedObjectsUtility); - - /// Local Flags - KRATOS_DEFINE_LOCAL_FLAG(INTERSECT_ELEMENTS); - KRATOS_DEFINE_LOCAL_FLAG(INTERSECT_CONDITIONS); - - /// Definition of the point type - using PointType = Point; - - /// Definition of the node type - using NodeType = Node<3>; - - /// Definition of the geometry type - using GeometryType = Geometry; - - /// Octree definitions - using ConfigurationType = Internals::DistanceSpatialContainersConfigure; - using CellType = OctreeBinaryCell; - using OctreeType = OctreeBinary; - using OctreePointerType = unique_ptr; - using CellNodeDataType = typename ConfigurationType::cell_node_data_type; - using OctreeCellVectorType = std::vector; - - ///@} - ///@name Life Cycle - ///@{ - - /// Constructor. - FindIntersectedObjectsUtility(ModelPart& rThisModelPart, Parameters ThisParameters = Parameters()); - - /// Destructor. - ~FindIntersectedObjectsUtility(){} - - ///@} - ///@name Operators - ///@{ - - - ///@} - ///@name Operations - ///@{ - - void UpdateSearchStructure(); - - void FindIntersectedObjects(GeometryType::Pointer rGeometry, PointerVector& rResults); - - ///@} - ///@name Access - ///@{ - - - ///@} - ///@name Inquiry - ///@{ - - - ///@} - ///@name Input and output - ///@{ - - /// Turn back information as a string. - std::string Info() const { - return "FindIntersectedObjectsUtility"; - } - - /// Print information about this object. - void PrintInfo(std::ostream& rOStream) const { - rOStream << Info(); - } - - /// Print object's data. - void PrintData(std::ostream& rOStream) const {} - - ///@} - ///@name Friends - ///@{ - - - ///@} - -private: - ///@name Static Member Variables - ///@{ - - - ///@} - ///@name Member Variables - ///@{ - - Flags mOptions; - ModelPart& mrModelPart; - OctreePointerType mpOctree; - - ///@} - ///@name Private Operators - ///@{ - - - ///@} - ///@name Private Operations - ///@{ - - void GenerateOctree(); - - void SetOctreeBoundingBox(); - - void FindIntersectedObjects( - const GeometryType& rIntersectionGeometry, - OctreeCellVectorType& rLeaves, - PointerVector& rResults); - - ///@} - ///@name Private Access - ///@{ - - - ///@} - ///@name Private Inquiry - ///@{ - - - ///@} - ///@name Un accessible methods - ///@{ - - /// Assignment operator. - FindIntersectedObjectsUtility& operator=(FindIntersectedObjectsUtility const& rOther); - - /// Copy constructor. - FindIntersectedObjectsUtility(FindIntersectedObjectsUtility const& rOther); - - ///@} - -}; // Class FindIntersectedObjectsUtility - -///@} - -///@name Type Definitions -///@{ - - -///@} -///@name Input and output -///@{ - - -/// input stream function -inline std::istream& operator >> (std::istream& rIStream, FindIntersectedObjectsUtility& rThis); - -/// output stream function -inline std::ostream& operator << (std::ostream& rOStream, const FindIntersectedObjectsUtility& rThis) -{ - rThis.PrintInfo(rOStream); - rOStream << std::endl; - rThis.PrintData(rOStream); - - return rOStream; -} -///@} - -} // namespace Kratos. - -#endif // KRATOS_FIND_INTERSECTED_ELEMENTS_UTILITY_H_INCLUDED defined From 4b73f548000cb808e91192f99996795a21dfb189 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Tue, 16 Nov 2021 16:35:01 +0100 Subject: [PATCH 41/43] update python interface --- .../coupling/depth_integration_input_process.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py index b2783d6ebc90..8f1abbcc576e 100644 --- a/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py @@ -188,7 +188,8 @@ def _SmoothDefaultValue(self): elapsed_time = self.interface_model_part.ProcessInfo.GetValue(KM.DELTA_TIME) semi_period = self.settings["semi_period_after_interval"].GetDouble() for variable in self.variables: - SW.ShallowWaterUtilities().SmoothTemporalVariable( - self.interface_model_part, + SW.ShallowWaterUtilities().SmoothHistoricalVariable( variable, + self.interface_model_part.Nodes, + elapsed_time, semi_period) From 2fbbee3a6640023de42dbc2483cbc951f6abf55c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Wed, 17 Nov 2021 10:27:51 +0100 Subject: [PATCH 42/43] fix runtime error if multiple instances --- .../custom_processes/depth_integration_process.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp index 0e83cacd7b04..462201301c4a 100644 --- a/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -52,7 +52,11 @@ DepthIntegrationProcess::DepthIntegrationProcess( mStoreHistorical = ThisParameters["store_historical_database"].GetBool(); mDirection = ThisParameters["direction_of_integration"].GetVector(); mDirection /= norm_2(mDirection); - mpIntegrationModelPart = &rModel.CreateModelPart("integration_auxiliary_model_part"); + if (rModel.HasModelPart("integration_auxiliary_model_part")) { // This is to allow multiple instances of this process + mpIntegrationModelPart = &rModel.GetModelPart("integration_auxiliary_model_part"); + } else { + mpIntegrationModelPart = &rModel.CreateModelPart("integration_auxiliary_model_part"); + } } void DepthIntegrationProcess::Execute() @@ -70,7 +74,6 @@ void DepthIntegrationProcess::Execute() void DepthIntegrationProcess::InitializeIntegrationModelPart() { - // Empty the model part VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Nodes()); VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Elements()); VariableUtils().SetFlag(TO_ERASE, true, mpIntegrationModelPart->Conditions()); From a6aa9d9c735c6ec89bc892e51ff6fbf14a6bcae9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20Mas=C3=B3=20Sotomayor?= Date: Mon, 22 Nov 2021 12:39:32 +0100 Subject: [PATCH 43/43] update test --- .../tests/cpp_tests/test_depth_integration_process.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/applications/ShallowWaterApplication/tests/cpp_tests/test_depth_integration_process.cpp b/applications/ShallowWaterApplication/tests/cpp_tests/test_depth_integration_process.cpp index 01d974d8d2a7..0fb92582e5cb 100644 --- a/applications/ShallowWaterApplication/tests/cpp_tests/test_depth_integration_process.cpp +++ b/applications/ShallowWaterApplication/tests/cpp_tests/test_depth_integration_process.cpp @@ -151,13 +151,13 @@ KRATOS_TEST_CASE_IN_SUITE(DepthIntegrationProcess3D, ShallowWaterApplicationFast DepthIntegrationProcess(model, process_parameters).Execute(); std::vector> reference; - reference.push_back({0.518519, 0.0, 0.0}); - reference.push_back({0.759259, 0.0, 0.0}); + reference.push_back({0.462963, 0.0, 0.0}); + reference.push_back({0.574074, 0.0, 0.0}); reference.push_back({0.851852, 0.0, 0.0}); - reference.push_back({1.092592, 0.0, 0.0}); - reference.push_back({1.240740, 0.0, 0.0}); + reference.push_back({1.129630, 0.0, 0.0}); + reference.push_back({1.240741, 0.0, 0.0}); - for(std::size_t i = 0; i < r_interface_model_part.NumberOfNodes(); ++i) { + for (std::size_t i = 0; i < r_interface_model_part.NumberOfNodes(); ++i) { auto i_node = r_interface_model_part.NodesBegin() + i; KRATOS_CHECK_VECTOR_NEAR(i_node->GetValue(VELOCITY), reference[i], 1e-6); }