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..462201301c4a --- /dev/null +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.cpp @@ -0,0 +1,178 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// 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_2d_2.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 "processes/find_intersected_geometrical_objects_process.h" + +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 + ) : Process(), + mrVolumeModelPart(rModel.GetModelPart(ThisParameters["volume_model_part_name"].GetString())), + 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); + 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() +{ + double bottom, top; + InitializeIntegrationModelPart(); + GetBoundingVolumeLimits(bottom, top); + CreateIntegrationLines(bottom, top); + FindIntersectedGeometricalObjectsProcess find_intersected_objects_process(*mpIntegrationModelPart, mrVolumeModelPart); + find_intersected_objects_process.ExecuteInitialize(); + find_intersected_objects_process.FindIntersections(); + auto intersected_objects = find_intersected_objects_process.GetIntersections(); + Integrate(intersected_objects); +} + +void DepthIntegrationProcess::InitializeIntegrationModelPart() +{ + 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) +{ + using MultipleReduction = CombinedReduction,MaxReduction>; + + std::tie(rMin, rMax) = block_for_each(mrVolumeModelPart.Nodes(), [&](NodeType& node){ + const double distance = inner_prod(mDirection, node); + return std::make_tuple(distance, distance); + }); +} + +void DepthIntegrationProcess::CreateIntegrationLines(const double Low, const double High) +{ + // Set the element name + const auto 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 + 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.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; + Integrate(objects_in_line, *i_node); + }); +} + +void DepthIntegrationProcess::Integrate(PointerVector& rObjects, NodeType& rNode) +{ + 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); +} + +int DepthIntegrationProcess::Check() +{ + 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; +} + +} // 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..4876c7b66015 --- /dev/null +++ b/applications/ShallowWaterApplication/custom_processes/depth_integration_process.h @@ -0,0 +1,241 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// 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; + + /** + * @brief Constructor with Model and Parameters + */ + DepthIntegrationProcess(Model& rModel, Parameters ThisParameters = Parameters()); + + /** + * @brief 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; + ModelPart* mpIntegrationModelPart; + array_1d mDirection; + bool mStoreHistorical; + + ///@} + ///@name Member Variables + ///@{ + + + ///@} + ///@name Private Operators + ///@{ + + + ///@} + ///@name Private Operations + ///@{ + + void InitializeIntegrationModelPart(); + + void GetBoundingVolumeLimits(double& rMin, double& rMax); + + void CreateIntegrationLines(const double Low, const double High); + + void Integrate(std::vector>& rResults); + + void Integrate(PointerVector& rObjects, NodeType& rNode); + + template> + void SetValue(NodeType& rNode, const TVarType& rVariable, TDataType& rValue) + { + if (mStoreHistorical) + rNode.FastGetSolutionStepValue(rVariable) = rValue; + else + rNode.GetValue(rVariable) = rValue; + } + + ///@} + ///@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 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. 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..8f1abbcc576e --- /dev/null +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_input_process.py @@ -0,0 +1,195 @@ +import KratosMultiphysics as KM +import KratosMultiphysics.ShallowWaterApplication as SW +import KratosMultiphysics.MappingApplication as Mapping +from KratosMultiphysics.kratos_utilities import GenerateVariableListFromInput +from KratosMultiphysics.HDF5Application import import_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): + raise Exception("expected input shall be a Parameters object, encapsulating a json string") + return DepthIntegrationInputProcess(model, settings["Parameters"]) + +class DepthIntegrationInputProcess(KM.OutputProcess): + """DepthIntegrationInputProcess + + Read the depth integrated values from an HDF5 file and set them as boundary conditions. + """ + + 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" : 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.""" + + KM.OutputProcess.__init__(self) + self.settings = settings + self.settings.ValidateAndAssignDefaults(self.GetDefaultParameters()) + + 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) + self.variables = GenerateVariableListFromInput(self.settings["list_of_variables"]) + self.variables_to_fix = GenerateVariableListFromInput(self.settings["list_of_variables_to_fix"]) + + 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_import.Check() + self.hdf5_process.Check() + + + def ExecuteInitialize(self): + '''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 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() + self._CheckInputVariables() + self._MapToBoundaryCondition() + else: + 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): + # 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 (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.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]))]) + 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.interface_model_part.ProcessInfo.GetValue(KM.TIME) + closest_time = next(filter(lambda x: x>current_time, self.times)) + self.input_model_part.ProcessInfo.SetValue(KM.TIME, closest_time) + + + def _SetDefaultTime(self): + default_time = self.settings["default_time_after_interval"].GetDouble() + 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.input_model_part.Nodes) + SW.ShallowWaterUtilities().SwapYZComponents(KM.VELOCITY, self.input_model_part.Nodes) + else: + 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.input_model_part.Nodes) + KM.VariableUtils().SetVariableToZero(KM.VELOCITY_Z, self.input_model_part.Nodes) + else: + 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): + for variable in self.variables: + if self.settings["read_historical_database"].GetBool(): + self.mapper.Map(variable, variable) + else: + self.mapper.Map(variable, variable, KM.Mapper.FROM_NON_HISTORICAL) + + for variable in self.variables_to_fix: + KM.VariableUtils().ApplyFixity(variable, True, self.interface_model_part.Nodes) + + + def _CreateMapper(self): + mapper_settings = KM.Parameters("""{ + "mapper_type": "nearest_neighbor", + "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.interface_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.input_model_part, + self.interface_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.interface_model_part.ProcessInfo.GetValue(KM.DELTA_TIME) + semi_period = self.settings["semi_period_after_interval"].GetDouble() + for variable in self.variables: + SW.ShallowWaterUtilities().SmoothHistoricalVariable( + variable, + self.interface_model_part.Nodes, + elapsed_time, + semi_period) 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..11a45f1ac5d4 --- /dev/null +++ b/applications/ShallowWaterApplication/python_scripts/coupling/depth_integration_output_process.py @@ -0,0 +1,146 @@ +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" : "", + "output_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) + settings.ValidateAndAssignDefaults(self.GetDefaultParameters()) + + 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) + self.variables = [KM.VELOCITY, KM.MOMENTUM, SW.HEIGHT] + + 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["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"]}""") + 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_process_settings, model) + + + def Check(self): + '''Check the processes.''' + self.integration_process.Check() + self.hdf5_process.Check() + + + def ExecuteInitialize(self): + '''Initialize the output model part.''' + self._InitializeOutputModelPart() + self._SetOutputProcessInfo() + if not self.store_historical: + 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.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.''' + self.integration_process.Execute() + self._MapToOutputModelPart() + + + 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() + + + def _InitializeOutputModelPart(self): + 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.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): + 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 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..0fb92582e5cb --- /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.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.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) { + 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