How to use the Spatial Containers

From KratosWiki
(Difference between revisions)
Jump to: navigation, search
(Using the search algorithm based in objects)
m (Calling the Constructor)
 
(7 intermediate revisions by 2 users not shown)
Line 7: Line 7:
 
comprising the problem.  
 
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
This would be very CPU intensive, and would limit the application that use and need to evaluate
+
the contact forces ,comprising a very small number (a few thousand) of separates bodies.  
+
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.
 
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.
With this purpose was created in Kratos Program  a contact search algorithm based on spatial decomposition,  
+
so it can be used for any application and is unique in that it was implemented in a generic way.
+
 
Therefore it requires a contact search algorithm having the following characteristics:
 
Therefore it requires a contact search algorithm having the following characteristics:
  
Line 23: Line 17:
 
*RAM efficient
 
*RAM efficient
  
The space decomposition implemented in Kratos are:
+
The space decomposition algorithms implemented in Kratos are:
  
* Trees : Kd-Trees, Octrees   
+
* '''Trees''' : Kd-Trees, Octrees   
* Bins  
+
* '''Bins'''  : Dynamic Bins, Static Bins
  
 
There are two  kind of searches algorithms in Kratos:
 
There are two  kind of searches algorithms in Kratos:
  
* Points Searches: Used commonly for searching near points in a radius  and can be used for
+
* '''Points Searches''': Used commonly for searching near points in a radius  and can be used for searches contacts for spheric and circle geometries. Use the Kd_tree,Octree and bins for spatial decomposition.     
searches contacts for spheric and circle geometries. Use the Kd_tree,Octree and bins for spatial
+
decomposition.     
+
  
* Object Searches: The searches is based in the concept of objects. Is a generic class, thats mean, we can use either objects type even points. For now, use the bins for spatial decomposition. It could be interesting proof hybrid methods would be a
+
* '''Object Searches''': The searches is based in the concept of objects. Is a generic class, thats mean, we can use either objects type even points. For now, use the bins for spatial decomposition.
and future investigation lines.
+
  
== Spatial decomposition based in bins. Main methods  ==
+
== The User Configure File ==
  
Exist implemented in Kratos  two versions of  domain decomposition based in bins:
+
The configuration file is an user-made file where the specific characteristics of the problem.
dynamic and static bins. The difference between them is the dynamic allocation of memory.
+
In simple terms, it can be considered as a small template class that the user needs to fill with
However both define the methods.  The most common methods used are presented and described below:
+
the types, operations and containers necessary for the strategy to search the contacts.
  
* Bins based in points :
+
This file is divided of four parts:
    *SearchNearestPoint : Calculates the nearest point of a point.
+
    *SearchInRadius    : Get the list of point in a certain radius. 
+
  
* Bins based in Objects : 
+
*Definition of the types and their containers.
    * SearchObjects      : Get the list of contacts of a certain objects even itself.
+
*Function that calculates bounding box
    * SearchObjectsInner : The same as above, just without return itself.
+
*Function of intersection between objects
    * SearchConcts      : Get the list of all pair contacts.  
+
*Function of intercession between the object and cells.
       
+
  
== Using the search algorithm based in objects ==
+
== An Example Configure File ==
  
Before start using this algorithm, is necessary to look some aspects:
+
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"*.
  
* The bins dynamic or static is a generic class. That means, a template
+
=== Creating the Skeleton ===
parameter is used. This template parameter is a user configure file. 
+
  
* The user need to create and extra file ('''user configure file'''). This fail describes the types
+
The first thing we have to do is changing the default names and directives of the file.  
and containers used by the bins. It has also the intersection method for contacts.
+
This normal implies changing two tokens that appear repeatedly across the file:
As There are many geometric figures, so the intersection algorithm differs each other.
+
Most of intersection algorithm can be found in computational geometric books and
+
more of them is written c++ or c. The file commonly is called the configure file.     
+
  
* Pointers and iterator: the bins works by pointers and references. So, in the configure file,
+
  KRATOS_FILENAME_H_INCLUDED
   we need to know what is the class, a pointer to the class, and iterator of pointers class
+
    
   containers.
+
   ClassName
  
 +
that in our case would become:
  
With this in hand, below describe how to create a configuration file and how to use the contact search algorithm.
+
  KRATOS_EXAMPLE_CONFIGURE_INCLUDED
 +
 
 +
  ExampleConfigure
  
== The User Configure File ==
+
=== Definition of the types and containers ===
  
It say that is the main file created by the user,
+
Once our file *"example_configure_file.h"* is prepared we need to fill it with the
which defines the types, and operations containers
+
types and operations.  
necessary to make the spatial decomposition and the search for contacts.
+
Is then the parameter that defines Bins.
+
   
+
This consists 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.
+
In our file called "discrete_configure_file.h" need to create the class,
+
types and methods before mentioned. Is useful to create this file using
+
the template of Kratos found in kratos/templates directory. Only we need to do is
+
changes the default names and directives by the name of our file. Normally this file is stored
+
in the custom_utilities directory of the application.
+
  
 +
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
  
* Change the name of directives 
+
In our example, a possible definition of these types can be:
  #if !defined(KRATOS_DISCRETE_PARTICLE_CONFIGURE_INCLUDED)
+
  #define  KRATOS_DISCRETE_PARTICLE_CONFIGURE_INCLUDED'''
+
  
*Change the name of class   
+
<!--THIS IS THE ORIGINAL CODE:
 
   ///@}
 
   ///@}
   ///@name Kratos Classes
+
   ///@name Type Definitions
 
   ///@{
 
   ///@{
   template <std::size_t TDimension>
+
   typedef  Point<Dimension, double>                        PointType;
   class DiscreteParticleConfigure{'''
+
  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;-->
 +
  <span style='color:#008000; '>///@}</span>
 +
  <span style='color:#008000; '>///@name Type Definitions</span>
 +
  <span style='color:#008000; '>///@{</span>
 
    
 
    
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span>  Point<span style='color:#0000ff; '>&lt;</span>Dimension<span style='color:#0000ff; '>,</span> <span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>></span>                        PointType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span>  std<span style='color:#0000ff; '>::</span>vector<span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>></span><span style='color:#0000ff; '>::</span><span style='color:#0000ff; font-weight:bold; '>iterator</span>                  DistanceIteratorType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span>  ModelPart<span style='color:#0000ff; '>::</span>ElementsContainerType<span style='color:#0000ff; '>::</span>ContainerType ContainerType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span>  ModelPart<span style='color:#0000ff; '>::</span>ElementsContainerType<span style='color:#0000ff; '>::</span>ContainerType ResultContainerType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span>  ContainerType<span style='color:#0000ff; '>::</span>value_type                      PointerType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span>  ContainerType<span style='color:#0000ff; '>::</span><span style='color:#0000ff; font-weight:bold; '>iterator</span>                        IteratorType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span>  ResultContainerType<span style='color:#0000ff; '>::</span>value_type                ResultPointerType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span>  ResultContainerType<span style='color:#0000ff; '>::</span><span style='color:#0000ff; font-weight:bold; '>iterator</span>                  ResultIteratorType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span>  ContactPair<span style='color:#0000ff; '>&lt;</span>PointerType<span style='color:#0000ff; '>></span>                        ContactPairType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span>  std<span style='color:#0000ff; '>::</span>vector<span style='color:#0000ff; '>&lt;</span>ContactPairType<span style='color:#0000ff; '>></span>                    ContainerContactType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span>  ContainerContactType<span style='color:#0000ff; '>::</span><span style='color:#0000ff; font-weight:bold; '>iterator</span>                  IteratorContactType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span>  ContainerContactType<span style='color:#0000ff; '>::</span>value_type                PointerContactType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span>  std<span style='color:#0000ff; '>::</span>vector<span style='color:#0000ff; '>&lt;</span>PointerType<span style='color:#0000ff; '>></span><span style='color:#0000ff; '>::</span><span style='color:#0000ff; font-weight:bold; '>iterator</span>              PointerTypeIterator<span style='color:#0000ff; '>;</span>
  
* Typing the typenames, objects and container of the particle
+
{{Note|
 +
Please, note that these are only suggestions and any other types and containers different to ModelPart can be used.
 +
}}
  
      ///@}
+
=== Definition of the operations ===
      ///@name Type Definitions
+
      ///@{
+
      typedef  Point<Dimension, double>                        PointType;
+
      typedef  std::vector<double>::iterator                  DistanceIteratorType;
+
      typedef  ModelPart::ElementsContainerType::ContainerType ContainerType;
+
      typedef  ContainerType::value_type                      PointerType;
+
      typedef  ContainerType::iterator                        IteratorType;
+
      typedef  ModelPart::ElementsContainerType::ContainerType ResultContainerType;
+
      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;
+
  
 +
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.
  
{{Note|
+
The list of functions to be defined are:
Another types and container different to the model part  can be used.
+
}}
+
  
 +
* Computing the bounding box.
 +
* The Intersection function between objects.
 +
* Intersection cell and objects.
 +
* Distance function between two objects.
  
* Computing the bounding box
+
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.
  
    '''static inline void CalculateBoundingBox'''(const PointerType& rObject, PointType& rLowPoint, PointType& rHighPoint)
+
==== BoundingBox ====
        {
+
        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;
+
            }
+
        }
+
  
{{Note|
+
*'''static inline void CalculateBoundingBox'''(const PointerType& rObject, PointType& rLowPoint, PointType& rHighPoint):
It is necessary the name  ''' static inline void CalculateBoundingBox ''' is well type, beacause
+
This method is used to calculate the boundingbox of the element (rObject). Once calculated the Low and High points of the
the function is called by the bins.
+
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'''
  
  
* The Intersection function between objects.
+
As an example, for our spheric elements:
 +
<!--THIS IS THE ORIGINAL CODE:
 +
  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 bool Intersection(const PointerType& rObj_1, const PointerType& rObj_2){
+
  static inline void CalculateBoundingBox(const PointerType& rObject, PointType& rLowPoint, PointType& rHighPoint, 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);
+
    rHighPoint = rLowPoint  = rObject->GetGeometry().GetPoint(0);
        const double& radius_1 = rObj_1->GetGeometry()(0)->GetSolutionStepValue(RADIUS);
+
   
        const double& radius_2 = rObj_2->GetGeometry()(0)->GetSolutionStepValue(RADIUS);
+
    for(std::size_t i = 0; i < 3; i++){
        double radius_sum      = radius_1 + radius_2;
+
      rLowPoint[i]  += -Radius;
        bool intersect        = (distance_2 - radius_sum * radius_sum) <= 0;
+
      rHighPoint[i] += Radius;
        return intersect;
+
    }
        }
+
  }-->
  
 +
<span style='color:#0000ff; font-weight:bold; '>static</span> <span style='color:#0000ff; font-weight:bold; '>inline</span> <span style='color:#0000ff; font-weight:bold; '>void</span> CalculateBoundingBox<span style='color:#0000ff; '>(</span><span style='color:#0000ff; font-weight:bold; '>const</span> PointerType<span style='color:#0000ff; '>&amp;</span> rObject<span style='color:#0000ff; '>,</span> PointType<span style='color:#0000ff; '>&amp;</span> rLowPoint<span style='color:#0000ff; '>,</span> PointType<span style='color:#0000ff; '>&amp;</span> rHighPoint<span style='color:#0000ff; '>)</span> <span style='color:#0000ff; '>{</span>
 +
 
 +
    rHighPoint <span style='color:#0000ff; '>=</span> rLowPoint  <span style='color:#0000ff; '>=</span> rObject<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>.</span>GetPoint<span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
    <span style='color:#0000ff; font-weight:bold; '>double</span> radius <span style='color:#0000ff; '>=</span> rObject<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetSolutionStepValue<span style='color:#0000ff; '>(</span>RADIUS<span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
 
 +
    <span style='color:#0000ff; font-weight:bold; '>for</span><span style='color:#0000ff; '>(</span>std<span style='color:#0000ff; '>::</span><span style='color:#0000ff; font-weight:bold; '>size_t</span> i <span style='color:#0000ff; '>=</span> <span style='color:#800000; '>0</span><span style='color:#0000ff; '>;</span> i <span style='color:#0000ff; '>&lt;</span> <span style='color:#800000; '>3</span><span style='color:#0000ff; '>;</span> i<span style='color:#0000ff; '>+</span><span style='color:#0000ff; '>+</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>{</span>
 +
      rLowPoint<span style='color:#0000ff; '>[</span>i<span style='color:#0000ff; '>]</span>  <span style='color:#0000ff; '>+</span><span style='color:#0000ff; '>=</span> <span style='color:#0000ff; '>-</span>radius<span style='color:#0000ff; '>;</span>
 +
      rHighPoint<span style='color:#0000ff; '>[</span>i<span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>+</span><span style='color:#0000ff; '>=</span> radius<span style='color:#0000ff; '>;</span>
 +
    <span style='color:#0000ff; '>}</span>
 +
<span style='color:#0000ff; '>}</span>
  
* Intersection cell and objects
+
<span style='color:#0000ff; font-weight:bold; '>static</span> <span style='color:#0000ff; font-weight:bold; '>inline</span> <span style='color:#0000ff; font-weight:bold; '>void</span> CalculateBoundingBox<span style='color:#0000ff; '>(</span><span style='color:#0000ff; font-weight:bold; '>const</span> PointerType<span style='color:#0000ff; '>&amp;</span> rObject<span style='color:#0000ff; '>,</span> PointType<span style='color:#0000ff; '>&amp;</span> rLowPoint<span style='color:#0000ff; '>,</span> PointType<span style='color:#0000ff; '>&amp;</span> rHighPoint<span style='color:#0000ff; '>,</span> <span style='color:#0000ff; font-weight:bold; '>const</span> <span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>&amp;</span> Radius<span style='color:#0000ff; '>)</span> <span style='color:#0000ff; '>{</span>
 +
   
 +
    rHighPoint <span style='color:#0000ff; '>=</span> rLowPoint  <span style='color:#0000ff; '>=</span> rObject<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>.</span>GetPoint<span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
   
 +
    <span style='color:#0000ff; font-weight:bold; '>for</span><span style='color:#0000ff; '>(</span>std<span style='color:#0000ff; '>::</span><span style='color:#0000ff; font-weight:bold; '>size_t</span> i <span style='color:#0000ff; '>=</span> <span style='color:#800000; '>0</span><span style='color:#0000ff; '>;</span> i <span style='color:#0000ff; '>&lt;</span> <span style='color:#800000; '>3</span><span style='color:#0000ff; '>;</span> i<span style='color:#0000ff; '>+</span><span style='color:#0000ff; '>+</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>{</span>
 +
      rLowPoint<span style='color:#0000ff; '>[</span>i<span style='color:#0000ff; '>]</span>  <span style='color:#0000ff; '>+</span><span style='color:#0000ff; '>=</span> <span style='color:#0000ff; '>-</span>Radius<span style='color:#0000ff; '>;</span>
 +
      rHighPoint<span style='color:#0000ff; '>[</span>i<span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>+</span><span style='color:#0000ff; '>=</span> Radius<span style='color:#0000ff; '>;</span>
 +
    <span style='color:#0000ff; '>}</span>
 +
<span style='color:#0000ff; '>}</span>
  
    static inline bool  IntersectionBox(const PointerType& rObject,  const PointType& rLowPoint, const PointType& rHighPoint){
+
==== Object-Object Intersection ====
        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 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'''.
  
  
{{Note|
+
As an example:
Intersection cell and objects: If we do not have this function a simple return true is enougth.  
+
<!--THIS IS THE ORIGINAL CODE:
}}
+
  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;
 +
  }
  
== Calling the Constructor ==
+
  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;
 +
  }-->
  
Once the configure file is done, we are able to use Contact Search Algorithm,
+
  <span style='color:#0000ff; font-weight:bold; '>static</span> <span style='color:#0000ff; font-weight:bold; '>inline</span> <span style='color:#0000ff; font-weight:bold; '>bool</span> Intersection<span style='color:#0000ff; '>(</span><span style='color:#0000ff; font-weight:bold; '>const</span> PointerType<span style='color:#0000ff; '>&amp;</span> rObj_1<span style='color:#0000ff; '>,</span> <span style='color:#0000ff; font-weight:bold; '>const</span> PointerType<span style='color:#0000ff; '>&amp;</span> rObj_2<span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>{</span>
based in  a domain decomposition. In the fail where it will be called , we need to add
+
   
in the directives the corresponding files, those are  the configure file and the
+
    array_1d<span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>,</span> <span style='color:#800000; '>3</span><span style='color:#0000ff; '>></span> rObj_2_to_rObj_1 <span style='color:#0000ff; '>=</span> rObj_1<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>.</span>GetPoint<span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span> <span style='color:#0000ff; '>-</span> rObj_2<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>.</span>GetPoint<span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
bins. Maybe is required some redefinition of the typedef defined above in the configure file.
+
    <span style='color:#0000ff; font-weight:bold; '>double</span> distance_2 <span style='color:#0000ff; '>=</span> inner_prod<span style='color:#0000ff; '>(</span>rObj_2_to_rObj_1<span style='color:#0000ff; '>,</span> rObj_2_to_rObj_1<span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
   
 +
    <span style='color:#0000ff; font-weight:bold; '>const</span> <span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>&amp;</span> radius_1 <span style='color:#0000ff; '>=</span> rObj_1<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetSolutionStepValue<span style='color:#0000ff; '>(</span>RADIUS<span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
    <span style='color:#0000ff; font-weight:bold; '>const</span> <span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>&amp;</span> radius_2 <span style='color:#0000ff; '>=</span> rObj_2<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetSolutionStepValue<span style='color:#0000ff; '>(</span>RADIUS<span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
   
 +
    <span style='color:#0000ff; font-weight:bold; '>double</span> radius_sum      <span style='color:#0000ff; '>=</span> radius_1 <span style='color:#0000ff; '>+</span> radius_2<span style='color:#0000ff; '>;</span>
 +
    <span style='color:#0000ff; font-weight:bold; '>bool</span> intersect        <span style='color:#0000ff; '>=</span> <span style='color:#0000ff; '>(</span>distance_2 <span style='color:#0000ff; '>-</span> radius_sum <span style='color:#0000ff; '>*</span> radius_sum<span style='color:#0000ff; '>)</span> <span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; '>=</span> <span style='color:#800000; '>0</span><span style='color:#0000ff; '>;</span>
 +
    <span style='color:#0000ff; font-weight:bold; '>return</span> intersect<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; '>}</span>
  
      static const std::size_t space_dim                = 2;
+
  <span style='color:#0000ff; font-weight:bold; '>static</span> <span style='color:#0000ff; font-weight:bold; '>inline</span> <span style='color:#0000ff; font-weight:bold; '>bool</span> Intersection<span style='color:#0000ff; '>(</span><span style='color:#0000ff; font-weight:bold; '>const</span> PointerType<span style='color:#0000ff; '>&amp;</span> rObj_1<span style='color:#0000ff; '>,</span> <span style='color:#0000ff; font-weight:bold; '>const</span> PointerType<span style='color:#0000ff; '>&amp;</span> rObj_2<span style='color:#0000ff; '>,</span> <span style='color:#0000ff; font-weight:bold; '>const</span> <span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>&amp;</span> Radius<span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>{</span>
      typedef DiscreteConfigure<space_dim>             Configure;
+
   
      typedef Configure::PointType                      PointType;
+
    array_1d<span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>,</span> <span style='color:#800000; '>3</span><span style='color:#0000ff; '>></span> rObj_2_to_rObj_1 <span style='color:#0000ff; '>=</span> rObj_1<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>.</span>GetPoint<span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span> <span style='color:#0000ff; '>-</span> rObj_2<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>.</span>GetPoint<span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
      typedef PointType::CoordinatesArrayType          CoordinatesArrayType;
+
    <span style='color:#0000ff; font-weight:bold; '>double</span> distance_2 <span style='color:#0000ff; '>=</span> inner_prod<span style='color:#0000ff; '>(</span>rObj_2_to_rObj_1<span style='color:#0000ff; '>,</span> rObj_2_to_rObj_1<span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
      typedef Configure::ContainerType                  ContainerType;
+
   
      typedef Configure::PointerType                    PointerType;
+
    <span style='color:#0000ff; font-weight:bold; '>const</span> <span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>&amp;</span> radius_2 <span style='color:#0000ff; '>=</span> rObj_2<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetSolutionStepValue<span style='color:#0000ff; '>(</span>RADIUS<span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
      typedef Configure::IteratorType                  IteratorType;
+
   
      typedef Configure::ResultContainerType            ResultContainerType;
+
    <span style='color:#0000ff; font-weight:bold; '>double</span> radius_sum      <span style='color:#0000ff; '>=</span> Radius <span style='color:#0000ff; '>+</span> radius_2<span style='color:#0000ff; '>;</span>
      typedef Configure::ResultPointerType              ResultPointerType;
+
    <span style='color:#0000ff; font-weight:bold; '>bool</span> intersect        <span style='color:#0000ff; '>=</span> <span style='color:#0000ff; '>(</span>distance_2 <span style='color:#0000ff; '>-</span> radius_sum <span style='color:#0000ff; '>*</span> radius_sum<span style='color:#0000ff; '>)</span> <span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; '>=</span> <span style='color:#800000; '>0</span><span style='color:#0000ff; '>;</span>
      typedef Configure::ResultIteratorType            ResultIteratorType;
+
    <span style='color:#0000ff; font-weight:bold; '>return</span> intersect<span style='color:#0000ff; '>;</span>
      typedef Configure::ContactPairType                ContactPairType;
+
  <span style='color:#0000ff; '>}</span>
      typedef Configure::ContainerContactType          ContainerContactType;
+
      typedef Configure::IteratorContactType            IteratorContactType;
+
      typedef Configure::PointerContactType            PointerContactType;
+
      typedef Configure::PointerTypeIterator            PointerTypeIterator;
+
      typedef ContainerContactType                      ContainerContactPair;
+
      typedef IteratorContactType                      IteratorContainerContactPair;
+
      typedef PointerContactType                        PointerContainerContactPair;
+
  
 +
==== Object-Cell Intersection ====
  
After that we call the constructor as follow:
+
*'''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'''.
  
* Using Dynamic Bins
+
*'''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).
    IteratorType it_begin    =  mBoundaryElements.begin();
+
If the intersection is not defined the default return value must be '''True'''.
    IteratorType it_end      =  mBoundaryElements.end()
+
    BinsObjectDynamic<Configure>  rBinsObjectDynamic(it_begin, it_end)
+
  
*Using Static Bins
 
  
     BinsObjectStatic<Configure> Bins(it_begin, it_end);
+
As an example:
     const std::size_t MaxNumberOfResults  = 1000;
+
<!--THIS IS THE ORIGINAL CODE:
     std::size_t NumberOfResults          = 0;
+
  static inline bool IntersectionBox(const PointerType& rObject,  const PointType& rLowPoint, const PointType& rHighPoint){
     ResultIteratorType begin;
+
      
 +
    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;
 +
  }-->
  
 +
  <span style='color:#0000ff; font-weight:bold; '>static</span> <span style='color:#0000ff; font-weight:bold; '>inline</span> <span style='color:#0000ff; font-weight:bold; '>bool</span> IntersectionBox<span style='color:#0000ff; '>(</span><span style='color:#0000ff; font-weight:bold; '>const</span> PointerType<span style='color:#0000ff; '>&amp;</span> rObject<span style='color:#0000ff; '>,</span>  <span style='color:#0000ff; font-weight:bold; '>const</span> PointType<span style='color:#0000ff; '>&amp;</span> rLowPoint<span style='color:#0000ff; '>,</span> <span style='color:#0000ff; font-weight:bold; '>const</span> PointType<span style='color:#0000ff; '>&amp;</span> rHighPoint<span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>{</span>
 +
   
 +
    array_1d<span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>,</span> <span style='color:#800000; '>3</span><span style='color:#0000ff; '>></span> center_of_particle <span style='color:#0000ff; '>=</span> rObject<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>.</span>GetPoint<span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
    <span style='color:#0000ff; font-weight:bold; '>const</span> <span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>&amp;</span> radius <span style='color:#0000ff; '>=</span> rObject<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetSolutionStepValue<span style='color:#0000ff; '>(</span>RADIUS<span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
   
 +
    <span style='color:#0000ff; font-weight:bold; '>bool</span> intersect <span style='color:#0000ff; '>=</span> <span style='color:#0000ff; '>(</span> rLowPoint<span style='color:#0000ff; '>[</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>]</span>  <span style='color:#0000ff; '>-</span> radius <span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; '>=</span> center_of_particle<span style='color:#0000ff; '>[</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>&amp;</span><span style='color:#0000ff; '>&amp;</span>
 +
                      rLowPoint<span style='color:#0000ff; '>[</span><span style='color:#800000; '>1</span><span style='color:#0000ff; '>]</span>  <span style='color:#0000ff; '>-</span> radius <span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; '>=</span> center_of_particle<span style='color:#0000ff; '>[</span><span style='color:#800000; '>1</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>&amp;</span><span style='color:#0000ff; '>&amp;</span>
 +
                      rLowPoint<span style='color:#0000ff; '>[</span><span style='color:#800000; '>2</span><span style='color:#0000ff; '>]</span>  <span style='color:#0000ff; '>-</span> radius <span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; '>=</span> center_of_particle<span style='color:#0000ff; '>[</span><span style='color:#800000; '>2</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>&amp;</span><span style='color:#0000ff; '>&amp;</span>
 +
                      rHighPoint<span style='color:#0000ff; '>[</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>+</span> radius <span style='color:#0000ff; '>></span><span style='color:#0000ff; '>=</span> center_of_particle<span style='color:#0000ff; '>[</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>&amp;</span><span style='color:#0000ff; '>&amp;</span>
 +
                      rHighPoint<span style='color:#0000ff; '>[</span><span style='color:#800000; '>1</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>+</span> radius <span style='color:#0000ff; '>></span><span style='color:#0000ff; '>=</span> center_of_particle<span style='color:#0000ff; '>[</span><span style='color:#800000; '>1</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>&amp;</span><span style='color:#0000ff; '>&amp;</span>
 +
                      rHighPoint<span style='color:#0000ff; '>[</span><span style='color:#800000; '>2</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>+</span> radius <span style='color:#0000ff; '>></span><span style='color:#0000ff; '>=</span> center_of_particle<span style='color:#0000ff; '>[</span><span style='color:#800000; '>2</span><span style='color:#0000ff; '>]</span>
 +
                      <span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
   
 +
    <span style='color:#0000ff; font-weight:bold; '>return</span>  intersect<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; '>}</span>
  
{{Note|
+
  <span style='color:#0000ff; font-weight:bold; '>static</span> <span style='color:#0000ff; font-weight:bold; '>inline</span> <span style='color:#0000ff; font-weight:bold; '>bool</span> IntersectionBox<span style='color:#0000ff; '>(</span><span style='color:#0000ff; font-weight:bold; '>const</span> PointerType<span style='color:#0000ff; '>&amp;</span> rObject<span style='color:#0000ff; '>,</span>  <span style='color:#0000ff; font-weight:bold; '>const</span> PointType<span style='color:#0000ff; '>&amp;</span> rLowPoint<span style='color:#0000ff; '>,</span> <span style='color:#0000ff; font-weight:bold; '>const</span> PointType<span style='color:#0000ff; '>&amp;</span> rHighPoint<span style='color:#0000ff; '>,</span> <span style='color:#0000ff; font-weight:bold; '>const</span> <span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>&amp;</span> Radius<span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>{</span>
     Note that we can use the dynamic or static bins. The difference between them
+
      
     is the dynamic allocation  of memory. Naturally, the static bins is faster
+
    array_1d<span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>,</span> <span style='color:#800000; '>3</span><span style='color:#0000ff; '>></span> center_of_particle <span style='color:#0000ff; '>=</span> rObject<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>.</span>GetPoint<span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
    because the we know size of containers.  
+
   
}}
+
    <span style='color:#0000ff; font-weight:bold; '>bool</span> intersect <span style='color:#0000ff; '>=</span> <span style='color:#0000ff; '>(</span> rLowPoint<span style='color:#0000ff; '>[</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>]</span>  <span style='color:#0000ff; '>-</span> Radius <span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; '>=</span> center_of_particle<span style='color:#0000ff; '>[</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>&amp;</span><span style='color:#0000ff; '>&amp;</span>
 +
                      rLowPoint<span style='color:#0000ff; '>[</span><span style='color:#800000; '>1</span><span style='color:#0000ff; '>]</span>  <span style='color:#0000ff; '>-</span> Radius <span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; '>=</span> center_of_particle<span style='color:#0000ff; '>[</span><span style='color:#800000; '>1</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>&amp;</span><span style='color:#0000ff; '>&amp;</span>
 +
                      rLowPoint<span style='color:#0000ff; '>[</span><span style='color:#800000; '>2</span><span style='color:#0000ff; '>]</span>  <span style='color:#0000ff; '>-</span> Radius <span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; '>=</span> center_of_particle<span style='color:#0000ff; '>[</span><span style='color:#800000; '>2</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>&amp;</span><span style='color:#0000ff; '>&amp;</span>
 +
                      rHighPoint<span style='color:#0000ff; '>[</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>+</span> Radius <span style='color:#0000ff; '>></span><span style='color:#0000ff; '>=</span> center_of_particle<span style='color:#0000ff; '>[</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>&amp;</span><span style='color:#0000ff; '>&amp;</span>
 +
                      rHighPoint<span style='color:#0000ff; '>[</span><span style='color:#800000; '>1</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>+</span> Radius <span style='color:#0000ff; '>></span><span style='color:#0000ff; '>=</span> center_of_particle<span style='color:#0000ff; '>[</span><span style='color:#800000; '>1</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>&amp;</span><span style='color:#0000ff; '>&amp;</span>
 +
                      rHighPoint<span style='color:#0000ff; '>[</span><span style='color:#800000; '>2</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>+</span> Radius <span style='color:#0000ff; '>></span><span style='color:#0000ff; '>=</span> center_of_particle<span style='color:#0000ff; '>[</span><span style='color:#800000; '>2</span><span style='color:#0000ff; '>]</span>
 +
                      <span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
   
 +
    <span style='color:#0000ff; font-weight:bold; '>return</span>  intersect<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; '>}</span>
 +
 
 +
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:
 +
<!--THIS IS THE ORIGINAL CODE:
 +
  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]) );
 +
  }-->
 +
 
 +
  <span style='color:#0000ff; font-weight:bold; '>static</span> <span style='color:#0000ff; font-weight:bold; '>inline</span> <span style='color:#0000ff; font-weight:bold; '>void</span> Distance<span style='color:#0000ff; '>(</span><span style='color:#0000ff; font-weight:bold; '>const</span> PointerType<span style='color:#0000ff; '>&amp;</span> rObj_1<span style='color:#0000ff; '>,</span> <span style='color:#0000ff; font-weight:bold; '>const</span> PointerType<span style='color:#0000ff; '>&amp;</span> rObj_2<span style='color:#0000ff; '>,</span> <span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>&amp;</span> distance<span style='color:#0000ff; '>)</span> <span style='color:#0000ff; '>{</span>
 +
   
 +
    array_1d<span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>,</span> <span style='color:#800000; '>3</span><span style='color:#0000ff; '>></span> center_of_particle1 <span style='color:#0000ff; '>=</span> rObj_1<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>.</span>GetPoint<span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
    array_1d<span style='color:#0000ff; '>&lt;</span><span style='color:#0000ff; font-weight:bold; '>double</span><span style='color:#0000ff; '>,</span> <span style='color:#800000; '>3</span><span style='color:#0000ff; '>></span> center_of_particle2 <span style='color:#0000ff; '>=</span> rObj_2<span style='color:#0000ff; '>-</span><span style='color:#0000ff; '>></span>GetGeometry<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>.</span>GetPoint<span style='color:#0000ff; '>(</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
   
 +
    distance <span style='color:#0000ff; '>=</span> <span style='color:#0000ff; font-weight:bold; '>sqrt</span><span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>(</span>center_of_particle1<span style='color:#0000ff; '>[</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>-</span> center_of_particle2<span style='color:#0000ff; '>[</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>]</span><span style='color:#0000ff; '>)</span> <span style='color:#0000ff; '>*</span> <span style='color:#0000ff; '>(</span>center_of_particle1<span style='color:#0000ff; '>[</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>-</span> center_of_particle2<span style='color:#0000ff; '>[</span><span style='color:#800000; '>0</span><span style='color:#0000ff; '>]</span><span style='color:#0000ff; '>)</span> <span style='color:#0000ff; '>+</span>
 +
                    <span style='color:#0000ff; '>(</span>center_of_particle1<span style='color:#0000ff; '>[</span><span style='color:#800000; '>1</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>-</span> center_of_particle2<span style='color:#0000ff; '>[</span><span style='color:#800000; '>1</span><span style='color:#0000ff; '>]</span><span style='color:#0000ff; '>)</span> <span style='color:#0000ff; '>*</span> <span style='color:#0000ff; '>(</span>center_of_particle1<span style='color:#0000ff; '>[</span><span style='color:#800000; '>1</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>-</span> center_of_particle2<span style='color:#0000ff; '>[</span><span style='color:#800000; '>1</span><span style='color:#0000ff; '>]</span><span style='color:#0000ff; '>)</span> <span style='color:#0000ff; '>+</span>
 +
                    <span style='color:#0000ff; '>(</span>center_of_particle1<span style='color:#0000ff; '>[</span><span style='color:#800000; '>2</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>-</span> center_of_particle2<span style='color:#0000ff; '>[</span><span style='color:#800000; '>2</span><span style='color:#0000ff; '>]</span><span style='color:#0000ff; '>)</span> <span style='color:#0000ff; '>*</span> <span style='color:#0000ff; '>(</span>center_of_particle1<span style='color:#0000ff; '>[</span><span style='color:#800000; '>2</span><span style='color:#0000ff; '>]</span> <span style='color:#0000ff; '>-</span> center_of_particle2<span style='color:#0000ff; '>[</span><span style='color:#800000; '>2</span><span style='color:#0000ff; '>]</span><span style='color:#0000ff; '>)</span> <span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; '>}</span>
 +
 
 +
== Calling the Constructor ==
 +
 
 +
Once the configure file is finished we are able to use Contact Search Algorithm.
 +
 
 +
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:
 +
<!--THIS IS THE ORIGINAL CODE:
 +
  #include "spatial_containers/bins_dinamyc_objects.h"
  
*Lastly, we can get the list of pairs of contacts
+
And define:
 +
<!--THIS IS THE ORIGINAL CODE:
 +
  static const std::size_t space_dim = 2;
 
    
 
    
           rBins.SearchContact(mPairContacts);
+
  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;-->
  
*However, is faster to do a loop and get the bodies that are in contact with a particular element.
+
  <span style='color:#0000ff; font-weight:bold; '>static</span> <span style='color:#0000ff; font-weight:bold; '>const</span> std<span style='color:#0000ff; '>::</span><span style='color:#0000ff; font-weight:bold; '>size_t</span> space_dim <span style='color:#0000ff; '>=</span> <span style='color:#800000; '>2</span><span style='color:#0000ff; '>;</span>
     
+
 
    For(IteratorType it =it_begin; it!=it_end; it++)
+
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> ExampleConfigure<span style='color:#0000ff; '>&lt;</span>space_dim<span style='color:#0000ff; '>></span>              Configure<span style='color:#0000ff; '>;</span> 
        rBinsObjectDynamic.SearchObjectsInner(*it, Result);   
+
 
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>PointType                      PointType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> PointType<span style='color:#0000ff; '>::</span>CoordinatesArrayType          CoordinatesArrayType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>ContainerType                  ContainerType<span style='color:#0000ff; '>;</span> 
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>PointerType                    PointerType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>IteratorType                   IteratorType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>ResultContainerType            ResultContainerType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>ResultPointerType              ResultPointerType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>ResultIteratorType            ResultIteratorType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>ContactPairType                ContactPairType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>ContainerContactType          ContainerContactType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>IteratorContactType            IteratorContactType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>PointerContactType            PointerContactType<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>PointerTypeIterator            PointerTypeIterator<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>ContainerContactType          ContainerContactPair<span style='color:#0000ff; '>;</span> 
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>IteratorContactType            IteratorContainerContactPair<span style='color:#0000ff; '>;</span>  
 +
  <span style='color:#0000ff; font-weight:bold; '>typedef</span> Configure<span style='color:#0000ff; '>::</span>PointerContactType            PointerContainerContactPair<span style='color:#0000ff; '>;</span>
  
where ''Result'' is the list of contact for a particular element.
+
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:
 +
<!--THIS IS THE ORIGINAL CODE:   
 +
  IteratorType it_begin =  mElements.begin();
 +
  IteratorType it_end =  mElements.end();
 +
 
 +
  BinsObjectDynamic<Configure>  rBinsObjectDynamic(it_begin, it_end);-->
 +
  IteratorType it_begin <span style='color:#0000ff; '>=</span> mElements<span style='color:#0000ff; '>.</span>begin<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
  IteratorType it_end <span style='color:#0000ff; '>=</span> mElements<span style='color:#0000ff; '>.</span>end<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
 +
   
 +
  BinsObjectDynamic<span style='color:#0000ff; '>&lt;</span>Configure<span style='color:#0000ff; '>></span>  rBinsObjectDynamic<span style='color:#0000ff; '>(</span>it_begin<span style='color:#0000ff; '>,</span> it_end<span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>
  
==References==
+
Finally we call the search. For example if we want to search the neighbours of our elements in a Radius N:
*  Munjiza. Ante, ''The Combined Finite-Discrete Element Method'', John Wiley & Sons, Ltd. 2004
+
<!--THIS IS THE ORIGINAL CODE:
 +
  ContainerType SearchElements;
 +
  double * Radius;
 +
  ResultContainerType Results;
 +
  double * ResultDistance;
 +
 
 +
  int MaxNumberOfElements = MAX_RES;
 +
 
 +
  ...
 +
 
 +
rBins.SearchObjectsInRadiusExclusive(SearElementPointerToGeometricalObjecPointerTemporalVector,Radius[i],Results.begin(),ResultsDistances.begin(),MaxNumberOfElements);-->
 +
 
 +
  ContainerType SearchElements<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>double</span> <span style='color:#0000ff; '>*</span> Radius<span style='color:#0000ff; '>;</span>
 +
  ResultContainerType Results<span style='color:#0000ff; '>;</span>
 +
  <span style='color:#0000ff; font-weight:bold; '>double</span> <span style='color:#0000ff; '>*</span> ResultDistance<span style='color:#0000ff; '>;</span>
 +
 
 +
  <span style='color:#0000ff; font-weight:bold; '>int</span> MaxNumberOfElements <span style='color:#0000ff; '>=</span> MAX_RES<span style='color:#0000ff; '>;</span>
 +
 
 +
  <span style='color:#0000ff; '>.</span><span style='color:#0000ff; '>.</span><span style='color:#0000ff; '>.</span>
 +
 
 +
  rBins<span style='color:#0000ff; '>.</span>SearchObjectsInRadiusExclusive<span style='color:#0000ff; '>(</span>SearElementPointerToGeometricalObjecPointerTemporalVector<span style='color:#0000ff; '>,</span>Radius<span style='color:#0000ff; '>[</span>i<span style='color:#0000ff; '>]</span><span style='color:#0000ff; '>,</span>Results<span style='color:#0000ff; '>.</span>begin<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>,</span>ResultsDistances<span style='color:#0000ff; '>.</span>begin<span style='color:#0000ff; '>(</span><span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>,</span>MaxNumberOfElements<span style='color:#0000ff; '>)</span><span style='color:#0000ff; '>;</span>

Latest revision as of 15:23, 17 July 2015

Contents

Using 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

The space decomposition algorithms implemented in Kratos are:

  • Trees : Kd-Trees, Octrees
  • Bins  : Dynamic Bins, Static Bins

There are two kind of searches algorithms in Kratos:

  • Points Searches: Used commonly for searching near points in a radius and can be used for searches contacts for spheric and circle geometries. Use the Kd_tree,Octree and bins for spatial decomposition.
  • Object Searches: The searches is based in the concept of objects. Is a generic class, thats mean, we can use either objects type even points. For now, use the bins for spatial decomposition.

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;

Warn_icon.gif Note:

Please, note that these are only suggestions and any other types and containers different to ModelPart can be used.

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.

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:

 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);
Personal tools
Categories