diff --git a/applications/SystemIdentificationApplication/custom_python/add_custom_sensors_to_python.cpp b/applications/SystemIdentificationApplication/custom_python/add_custom_sensors_to_python.cpp index 3035f5373183..4105cfb54128 100644 --- a/applications/SystemIdentificationApplication/custom_python/add_custom_sensors_to_python.cpp +++ b/applications/SystemIdentificationApplication/custom_python/add_custom_sensors_to_python.cpp @@ -71,7 +71,7 @@ void AddCustomSensorsToPython(pybind11::module& m) py::arg("weight"), py::arg("error_threshold") = Sensor::DefaultErrorThreshold) .def_static("GetDefaultParameters", &DisplacementSensor::GetDefaultParameters) - .def_static("Create", &DisplacementSensor::Create, py::arg("domain_model_part"), py::arg("sensor_model_part"), py::arg("sensor_id"), py::arg("sensor_parameters")) + .def_static("Create", &DisplacementSensor::Create, py::arg("domain_model_part"), py::arg("sensor_model_part"), py::arg("sensor_id"), py::arg("sensor_parameters"), py::arg("domain_bins")) ; auto strain_sensor = py::class_(sensor_module, "StrainSensor"); @@ -93,7 +93,7 @@ void AddCustomSensorsToPython(pybind11::module& m) py::arg("weight"), py::arg("error_threshold") = Sensor::DefaultErrorThreshold) .def_static("GetDefaultParameters", &StrainSensor::GetDefaultParameters) - .def_static("Create", &StrainSensor::Create, py::arg("domain_model_part"), py::arg("sensor_model_part"), py::arg("sensor_id"), py::arg("sensor_parameters")) + .def_static("Create", &StrainSensor::Create, py::arg("domain_model_part"), py::arg("sensor_model_part"), py::arg("sensor_id"), py::arg("sensor_parameters"), py::arg("domain_bins")) ; } diff --git a/applications/SystemIdentificationApplication/custom_sensors/displacement_sensor.cpp b/applications/SystemIdentificationApplication/custom_sensors/displacement_sensor.cpp index 9de915670db2..ee42653f0279 100644 --- a/applications/SystemIdentificationApplication/custom_sensors/displacement_sensor.cpp +++ b/applications/SystemIdentificationApplication/custom_sensors/displacement_sensor.cpp @@ -16,7 +16,6 @@ // External includes // Project includes -#include "utilities/brute_force_point_locator.h" // Application includes #include "custom_utilities/sensor_utils.h" @@ -48,8 +47,10 @@ DisplacementSensor::DisplacementSensor( const auto& r_geometry = rElement.GetGeometry(); const auto& current_sensor_location = *(this->GetNode()); + // same local coordinate tolerance as the domain bins used to find rElement, otherwise a + // point lying on a face shared by two elements may be rejected by the element the bins found. Point local_point; - if (r_geometry.IsInside(current_sensor_location, local_point)) { + if (r_geometry.IsInside(current_sensor_location, local_point, 1e-6)) { // point is within the geometry. Use shape function evaluations // from the element geometry to get the shape function values. r_geometry.ShapeFunctionsValues(mNs, local_point); @@ -104,7 +105,8 @@ Sensor::Pointer DisplacementSensor::Create( ModelPart& rDomainModelPart, ModelPart& rSensorModelPart, const IndexType Id, - Parameters SensorParameters) + Parameters SensorParameters, + GeometricalObjectsBins& rDomainBins) { KRATOS_TRY @@ -120,12 +122,13 @@ Sensor::Pointer DisplacementSensor::Create( << "Location of the sensor \"" << SensorParameters["name"].GetString() << "\" should have 3 components. [ location = " << location << " ].\n"; - Point loc(location[0], location[1], location[2]); - - Vector dummy_shape_functions; + const auto search_result = rDomainBins.SearchIsInside(Point(location[0], location[1], location[2])); + KRATOS_ERROR_IF_NOT(search_result.GetIsObjectFound()) + << "Location of the sensor \"" << SensorParameters["name"].GetString() + << "\" is not inside any element of " << rDomainModelPart.FullName() + << ". [ location = " << location << " ].\n"; - const auto element_id = BruteForcePointLocator(rDomainModelPart).FindElement(loc, dummy_shape_functions); - const auto& r_element = rDomainModelPart.GetElement(element_id); + const auto& r_element = rDomainModelPart.GetElement(search_result.Get()->Id()); auto p_node = rSensorModelPart.CreateNewNode(Id, location[0], location[1], location[2]); diff --git a/applications/SystemIdentificationApplication/custom_sensors/displacement_sensor.h b/applications/SystemIdentificationApplication/custom_sensors/displacement_sensor.h index 1e7d84e18410..efa29c58f395 100644 --- a/applications/SystemIdentificationApplication/custom_sensors/displacement_sensor.h +++ b/applications/SystemIdentificationApplication/custom_sensors/displacement_sensor.h @@ -19,6 +19,7 @@ // Project includes #include "includes/ublas_interface.h" #include "includes/element.h" +#include "spatial_containers/geometrical_objects_bins.h" // Application includes #include "sensor.h" @@ -64,11 +65,17 @@ class KRATOS_API(SYSTEM_IDENTIFICATION_APPLICATION) DisplacementSensor : public ///@name Static operations ///@{ + /** + * @brief Creates the sensor at the location given in SensorParameters. + * @details rDomainBins must be built from the elements of rDomainModelPart. It is + * shared by all sensors, so each sensor is located without a linear search. + */ static Sensor::Pointer Create( ModelPart& rDomainModelPart, ModelPart& rSensorModelPart, const IndexType Id, - Parameters SensorParameters); + Parameters SensorParameters, + GeometricalObjectsBins& rDomainBins); static Parameters GetDefaultParameters(); diff --git a/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.cpp b/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.cpp index 2a802a13fcbb..949a43e60212 100644 --- a/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.cpp +++ b/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.cpp @@ -16,7 +16,6 @@ // Project includes #include "includes/kratos_components.h" -#include "utilities/brute_force_point_locator.h" // Application includes #include "custom_utilities/sensor_utils.h" @@ -41,7 +40,8 @@ StrainSensor::StrainSensor( mStrainType(rStrainType), mrStrainVariable(rStrainVariable) { - KRATOS_ERROR_IF_NOT(rElement.GetGeometry().IsInside(*(this->GetNode()), mLocalPoint)) + // same local coordinate tolerance as the domain bins used to find rElement + KRATOS_ERROR_IF_NOT(rElement.GetGeometry().IsInside(*(this->GetNode()), mLocalPoint, 1e-6)) << "The point " << this->GetNode()->Coordinates() << " is not inside or on the boundary of the geometry of element with id " << mElementId << "."; @@ -75,7 +75,8 @@ Sensor::Pointer StrainSensor::Create( ModelPart& rDomainModelPart, ModelPart& rSensorModelPart, const IndexType Id, - Parameters SensorParameters) + Parameters SensorParameters, + GeometricalObjectsBins& rDomainBins) { KRATOS_TRY @@ -86,12 +87,13 @@ Sensor::Pointer StrainSensor::Create( << "Location of the sensor \"" << SensorParameters["name"].GetString() << "\" should have 3 components. [ location = " << location << " ].\n"; - Point loc(location[0], location[1], location[2]); - - Vector dummy_shape_functions; + const auto search_result = rDomainBins.SearchIsInside(Point(location[0], location[1], location[2])); + KRATOS_ERROR_IF_NOT(search_result.GetIsObjectFound()) + << "Location of the sensor \"" << SensorParameters["name"].GetString() + << "\" is not inside any element of " << rDomainModelPart.FullName() + << ". [ location = " << location << " ].\n"; - const auto element_id = BruteForcePointLocator(rDomainModelPart).FindElement(loc, dummy_shape_functions); - const auto& r_element = rDomainModelPart.GetElement(element_id); + const auto& r_element = rDomainModelPart.GetElement(search_result.Get()->Id()); auto p_node = rSensorModelPart.CreateNewNode(Id, location[0], location[1], location[2]); diff --git a/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.h b/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.h index 8b310cc80d52..4f7255e518c1 100644 --- a/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.h +++ b/applications/SystemIdentificationApplication/custom_sensors/strain_sensor.h @@ -18,6 +18,7 @@ // Project includes #include "includes/element.h" +#include "spatial_containers/geometrical_objects_bins.h" // Application includes #include "sensor.h" @@ -78,11 +79,17 @@ class KRATOS_API(SYSTEM_IDENTIFICATION_APPLICATION) StrainSensor : public Sensor ///@name Static operations ///@{ + /** + * @brief Creates the sensor at the location given in SensorParameters. + * @details rDomainBins must be built from the elements of rDomainModelPart. It is + * shared by all sensors, so each sensor is located without a linear search. + */ static Sensor::Pointer Create( ModelPart& rDomainModelPart, ModelPart& rSensorModelPart, const IndexType Id, - Parameters SensorParameters); + Parameters SensorParameters, + GeometricalObjectsBins& rDomainBins); static Parameters GetDefaultParameters(); diff --git a/applications/SystemIdentificationApplication/python_scripts/sensor_generator_analysis.py b/applications/SystemIdentificationApplication/python_scripts/sensor_generator_analysis.py index addff98d8224..4847102dfeb2 100644 --- a/applications/SystemIdentificationApplication/python_scripts/sensor_generator_analysis.py +++ b/applications/SystemIdentificationApplication/python_scripts/sensor_generator_analysis.py @@ -63,7 +63,8 @@ def __GenerationBoundingSurfaceBased(self, sensors_dict, sensor_group_params: Kr }""") sensor_group_params.ValidateAndAssignDefaults(defaults) - point_locator = Kratos.BruteForcePointLocator(self.model_part) + # spatial search built once for all generated points (mesh is in its initial configuration here) + bins = Kratos.GeometricalObjectsBins(self.model_part.Elements, 1e-8) corner_1 = sensor_group_params["bounding_surface_corner_1"].GetVector() corner_2 = sensor_group_params["bounding_surface_corner_2"].GetVector() @@ -83,9 +84,8 @@ def __GenerationBoundingSurfaceBased(self, sensors_dict, sensor_group_params: Kr for i_z in range(number_of_sensors[2]): z_coord = corner_1[2] + distance[2] * (i_z + 1) / (number_of_sensors[2] + 1) loc = Kratos.Point(x_coord, y_coord, z_coord) - shape_funcs = Kratos.Vector() - elem_id = point_locator.FindElement(loc, shape_funcs, Kratos.Configuration.Initial, 1e-8) - if elem_id != -1 and KratosSI.SensorUtils.IsPointInGeometry(loc, self.model_part.GetElement(elem_id).GetGeometry()): + result = bins.SearchIsInside(loc) + if result.IsObjectFound() and KratosSI.SensorUtils.IsPointInGeometry(loc, self.model_part.GetElement(result.Get().Id).GetGeometry()): current_params = sensor_params.Clone() current_params["location"].SetVector(loc) json_dict = json.loads(current_params.WriteJsonString().replace("", str(index))) diff --git a/applications/SystemIdentificationApplication/python_scripts/utilities/sensor_utils.py b/applications/SystemIdentificationApplication/python_scripts/utilities/sensor_utils.py index 8f998e332ba1..07084a44cf67 100644 --- a/applications/SystemIdentificationApplication/python_scripts/utilities/sensor_utils.py +++ b/applications/SystemIdentificationApplication/python_scripts/utilities/sensor_utils.py @@ -38,6 +38,10 @@ def CreateSensors(sensor_model_part: Kratos.ModelPart, domain_model_part: Kratos except: pass + # One spatial search over the domain elements, shared by all sensors, so each sensor's Create + # locates its element(s) without a linear search. Built lazily since empty bins are not searchable. + domain_bins: 'typing.Optional[Kratos.GeometricalObjectsBins]' = None + list_of_sensors: 'list[KratosSI.Sensors.Sensor]' = [] for parameters in list_of_parameters: if not parameters.Has("type"): @@ -47,7 +51,13 @@ def CreateSensors(sensor_model_part: Kratos.ModelPart, domain_model_part: Kratos if not sensor_type_name in dict_of_sensor_types.keys(): raise RuntimeError(f"Unsupported sensor type = \"{sensor_type_name}\" requested. Followings are supported:\n\t" + "\n\t".join(dict_of_sensor_types.keys())) - sensor: KratosSI.Sensors.Sensor = dict_of_sensor_types[sensor_type_name].Create(domain_model_part, sensor_model_part, len(list_of_sensors) + 1, parameters) + if domain_bins is None: + if domain_model_part.NumberOfElements() == 0: + raise RuntimeError(f"The domain model part \"{domain_model_part.FullName()}\" has no elements to locate the sensors in.") + # same local coordinate tolerance as the former brute force point locator + domain_bins = Kratos.GeometricalObjectsBins(domain_model_part.Elements, 1e-6) + + sensor: KratosSI.Sensors.Sensor = dict_of_sensor_types[sensor_type_name].Create(domain_model_part, sensor_model_part, len(list_of_sensors) + 1, parameters, domain_bins) list_of_sensors.append(sensor) return list_of_sensors