How to use the Spatial Search Algorithm
Contents |
The Spatial Search Algorithm
Large-scale in Finite Element Methods, Discrete Element Method and Combined Finite-Discrete Element Method simulations involve contact of a large number of separate bodies. Processing contact interaction for all possible contacts would involve a total number of operations proportional to N^2, where N is the total number of separates bodies comprising the problem.
This would be very CPU intensive, and would limit the overall performance of the application. To reduce CPU requirements of processing contact interaction, it is necessary to eliminate couples of discrete elements that are far from each other and are not in contact. A set procedures designed to detect bodies that are close to each other is usually called a contact detection algorithm, or sometimes a contact search algorithm.
With this purpose a set contact search algorithms based on spatial decomposition are provided in Kratos, so they can be used for any application. Therefore it requires a contact search algorithm having the following characteristics:
- Robust,
- CPU efficient,
- RAM efficient
Usage
Kratos search algorithm is presented in several variants but all of them follow a common structure that needs to be followed in order to make a correct use of the algorithms. All the search algorithms are composed of two parts, the *user configuration file* and the *strategy*
The User Configure File
The configuration file is an user-made file where the specific characteristics of the problem. In simple terms, it can be considered as a small template class that the user needs to fill with the types, operations and containers necessary for the strategy to search the contacts.
This file is divided of four parts:
- Definition of the types and their containers.
- Function that calculates bounding box
- Function of intersection between objects
- Function of intercession between the object and cells.
An Example Configure File
Suppose that we have a set of discrete elements like a spheres. The first step is to create a configuration file. The easiest way to do that is copying the template found in "KRATOS_ROOT/kratos/templates/header_template.h". In this example we are going to name the file *"example_configure_file.h"*.
Creating the Skeleton
The first thing we have to do is changing the default names and directives of the file. This normal implies changing two tokens that appear repeatedly across the file:
KRATOS_FILENAME_H_INCLUDED ClassName
that in our case would become:
KRATOS_EXAMPLE_CONFIGURE_INCLUDED ExampleConfigure
Definition of the types and containers
Once our file *"example_configure_file.h"* is prepared we need to fill it with the types and operations.
The proper place to define the files is under the "Type Definitions" section and the list of types that need to be defined is the following:
- PointType: Type of the points that form our element.
- DistanceIteratorType: Type of the iterator for the distances between contacts.
- ContainerType: Type of the elements container.
- ResultContainerType: Type of the elements container for the results. This is normally the same as ContainerType but it can vary in case our results are different elements. For example if we want to search all the tetrahedrons in a X distance from a set of triangles.
- PointerType: Type of the pointer to our elements
- IteratorType: Type of the iterator to our elements
- ResultPointerType: Type of the pointer to our result elements
- ResultIteratorType: Type of the iterator to our result elements
- ContactPairType: Type of the contact pair
- ContainerContactType:
- IteratorContactType:
- PointerContactType:
- PointerTypeIterator: Type of the pointer iterator
In our example, a possible definition of these types can be:
///@}
///@name Type Definitions
///@{
typedef Point<Dimension, double> PointType;
typedef std::vector<double>::iterator DistanceIteratorType;
typedef ModelPart::ElementsContainerType::ContainerType ContainerType;
typedef ModelPart::ElementsContainerType::ContainerType ResultContainerType;
typedef ContainerType::value_type PointerType;
typedef ContainerType::iterator IteratorType;
typedef ResultContainerType::value_type ResultPointerType;
typedef ResultContainerType::iterator ResultIteratorType;
typedef ContactPair<PointerType> ContactPairType;
typedef std::vector<ContactPairType> ContainerContactType;
typedef ContainerContactType::iterator IteratorContactType;
typedef ContainerContactType::value_type PointerContactType;
typedef std::vector<PointerType>::iterator PointerTypeIterator;
|
Please, note that these are only suggestions and any other types and containers different to ModelPart can be used. |
Note:
Definition of the operations
Along with the definition of the types of the previous section one also needs to define a set of basic operations that will be used by the strategy. This operations involve mainly interaction between objects which have to be defined apart as they are bound to the characteristics of the elements.
The list of functions to be defined are:
- Computing the bounding box.
- The Intersection function between objects.
- Intersection cell and objects.
- Distance function between two objects.
Please note that the signature of all functions is mandatory as the strategy will only be able to use them if they appear exactly as described.
BoundingBox
- static inline void CalculateBoundingBox(const PointerType& rObject, PointType& rLowPoint, PointType& rHighPoint):
This method is used to calculate the boundingbox of the element (rObject). Once calculated the Low and High points of the boundingbox need to be stored in rLowPoint and rHighPoint
- static inline void CalculateBoundingBox(const PointerType& rObject, PointType& rLowPoint, PointType& rHighPoint, const double& Radius):
This method is used to calculate the boundingbox of the element (rObject) given an imposed radius. Once calculated the Low and High points of the boundingbox need to be stored in rLowPoint and rHighPoint
As an example, for our spheric elements:
static inline void CalculateBoundingBox(const PointerType& rObject, PointType& rLowPoint, PointType& rHighPoint) {
rHighPoint = rLowPoint = rObject->GetGeometry().GetPoint(0);
double radius = rObject->GetGeometry()(0)->GetSolutionStepValue(RADIUS);
for(std::size_t i = 0; i < 3; i++){
rLowPoint[i] += -radius;
rHighPoint[i] += radius;
}
}
static inline void CalculateBoundingBox(const PointerType& rObject, PointType& rLowPoint, PointType& rHighPoint, const double& Radius) {
rHighPoint = rLowPoint = rObject->GetGeometry().GetPoint(0);
for(std::size_t i = 0; i < 3; i++){
rLowPoint[i] += -Radius;
rHighPoint[i] += Radius;
}
}
Object-Object Intersection
- static inline bool Intersection(const PointerType& rObj_1, const PointerType& rObj_2):
This method is used to test if two given elements (rObj_1, rObj_2) intersect with each other. If the intersection is not defined the default return value must be True.
- static inline bool Intersection(const PointerType& rObj_1, const PointerType& rObj_2, const double& Radius):
This method is used to test if elements rObj_2 intersect with rObj_1 in a radius (Radius) . If the intersection is not defined the default return value must be True.
As an example:
static inline bool Intersection(const PointerType& rObj_1, const PointerType& rObj_2){
array_1d<double, 3> rObj_2_to_rObj_1 = rObj_1->GetGeometry().GetPoint(0) - rObj_2->GetGeometry().GetPoint(0);
double distance_2 = inner_prod(rObj_2_to_rObj_1, rObj_2_to_rObj_1);
const double& radius_1 = rObj_1->GetGeometry()(0)->GetSolutionStepValue(RADIUS);
const double& radius_2 = rObj_2->GetGeometry()(0)->GetSolutionStepValue(RADIUS);
double radius_sum = radius_1 + radius_2;
bool intersect = (distance_2 - radius_sum * radius_sum) <= 0;
return intersect;
}
static inline bool Intersection(const PointerType& rObj_1, const PointerType& rObj_2, const double& Radius){
array_1d<double, 3> rObj_2_to_rObj_1 = rObj_1->GetGeometry().GetPoint(0) - rObj_2->GetGeometry().GetPoint(0);
double distance_2 = inner_prod(rObj_2_to_rObj_1, rObj_2_to_rObj_1);
const double& radius_2 = rObj_2->GetGeometry()(0)->GetSolutionStepValue(RADIUS);
double radius_sum = Radius + radius_2;
bool intersect = (distance_2 - radius_sum * radius_sum) <= 0;
return intersect;
}
Object-Cell Intersection
- static inline bool IntersectionBox(const PointerType& rObject, const PointType& rLowPoint, const PointType& rHighPoint):
This method is used to test if a given element (rObject) intersect with a cell defined by rLowPoint and rHighPoint. If the intersection is not defined the default return value must be True.
- static inline bool IntersectionBox(const PointerType& rObject, const PointType& rLowPoint, const PointType& rHighPoint, const double& Radius):
This method is used to test if a given element (rObject) intersect with a cell defined by rLowPoint and rHighPoint in a radius (Radius). If the intersection is not defined the default return value must be True.
As an example:
static inline bool IntersectionBox(const PointerType& rObject, const PointType& rLowPoint, const PointType& rHighPoint){
array_1d<double, 3> center_of_particle = rObject->GetGeometry().GetPoint(0);
const double& radius = rObject->GetGeometry()(0)->GetSolutionStepValue(RADIUS);
bool intersect = ( rLowPoint[0] - radius <= center_of_particle[0] &&
rLowPoint[1] - radius <= center_of_particle[1] &&
rLowPoint[2] - radius <= center_of_particle[2] &&
rHighPoint[0] + radius >= center_of_particle[0] &&
rHighPoint[1] + radius >= center_of_particle[1] &&
rHighPoint[2] + radius >= center_of_particle[2]
);
return intersect;
}
static inline bool IntersectionBox(const PointerType& rObject, const PointType& rLowPoint, const PointType& rHighPoint, const double& Radius){
array_1d<double, 3> center_of_particle = rObject->GetGeometry().GetPoint(0);
bool intersect = ( rLowPoint[0] - Radius <= center_of_particle[0] &&
rLowPoint[1] - Radius <= center_of_particle[1] &&
rLowPoint[2] - Radius <= center_of_particle[2] &&
rHighPoint[0] + Radius >= center_of_particle[0] &&
rHighPoint[1] + Radius >= center_of_particle[1] &&
rHighPoint[2] + Radius >= center_of_particle[2]
);
return intersect;
}
Please notice that if the object-cell function returns True the strategy will behave as a brute froce search.
Distance function
- static inline void Distance(const PointerType& rObj_1, const PointerType& rObj_2, double& distance):
This function returns the distance between to given elements (rObj_1, rObj_2). Result must be stored in "distance"
For instance:
static inline void Distance(const PointerType& rObj_1, const PointerType& rObj_2, double& distance) {
array_1d<double, 3> center_of_particle1 = rObj_1->GetGeometry().GetPoint(0);
array_1d<double, 3> center_of_particle2 = rObj_2->GetGeometry().GetPoint(0);
distance = sqrt((center_of_particle1[0] - center_of_particle2[0]) * (center_of_particle1[0] - center_of_particle2[0]) +
(center_of_particle1[1] - center_of_particle2[1]) * (center_of_particle1[1] - center_of_particle2[1]) +
(center_of_particle1[2] - center_of_particle2[2]) * (center_of_particle1[2] - center_of_particle2[2]) );
}
Calling the Constructor
Once the configure file is finished, we are able to use Contact Search Algorithm, based in a domain decomposition.
In the file where our search is going to be called we will need to include our configuration file along with the desired strategy. The available strategies can be found in "KRATOS_ROOT/kratos/spatial_containers" and basically are:
- Bins_static
- Bins_static_objects
- Bins_dynamic
- Bins_dynamic_objects
Typically dynamic strategies are slower but are not bound to a specific amount of memory. Static strategies are faster but we need to anticipate the memory requirements.
For our example we are going to select "bins_dinamyc_objects.h". Hence our file will have to include:
#include "spatial_containers/bins_dinamyc_objects.h"
And define:
static const std::size_t space_dim = 2; typedef ExampleConfigure<space_dim> Configure; typedef Configure::PointType PointType; typedef PointType::CoordinatesArrayType CoordinatesArrayType; typedef Configure::ContainerType ContainerType; typedef Configure::PointerType PointerType; typedef Configure::IteratorType IteratorType; typedef Configure::ResultContainerType ResultContainerType; typedef Configure::ResultPointerType ResultPointerType; typedef Configure::ResultIteratorType ResultIteratorType; typedef Configure::ContactPairType ContactPairType; typedef Configure::ContainerContactType ContainerContactType; typedef Configure::IteratorContactType IteratorContactType; typedef Configure::PointerContactType PointerContactType; typedef Configure::PointerTypeIterator PointerTypeIterator; typedef Configure::ContainerContactType ContainerContactPair; typedef Configure::IteratorContactType IteratorContainerContactPair; typedef Configure::PointerContactType PointerContainerContactPair;
After that we can create the search strategy itself. In order to create our search strategy we will need to pass all the elements in our domain:
IteratorType it_begin = mElements.begin(); IteratorType it_end = mElements.end(); BinsObjectDynamic<Configure> rBinsObjectDynamic(it_begin, it_end);
Finally we call the search. For example if we want to search the neighbours of our elements in a Radius N:
ContainerType SearchElements; double * Radius; ResultContainerType Results; double * ResultDistance; int MaxNumberOfElements = MAX_RES; ... rBins.SearchObjectsInRadiusExclusive(SearElementPointerToGeometricalObjecPointerTemporalVector,Radius[i],Results.begin(),ResultsDistances.begin(),MaxNumberOfElements);
