How to use the Spatial Search Algorithm

Revision as of 15:28, 4 May 2012

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 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.

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:

• Robust,
• CPU efficient,
• RAM efficient

The User Configure File

It say that is the main file created by the user, which defines the types, and operations contenerdores 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 Exemple 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.

• Change the name of directives

 #if !defined(KRATOS_DISCRETE_PARTICLE_CONFIGURE_INCLUDED)
 #define  KRATOS_DISCRETE_PARTICLE_CONFIGURE_INCLUDED

• Change the name of class

 ///@}
 ///@name Kratos Classes
 ///@{
 template <std::size_t TDimension>
 class DiscreteParticleConfigure{
 

• Typing the typenames, objects and container of the particle

     ///@}
     ///@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;

• Computing the bounding box

    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;
           }
       }

• The Intersection function between objects.

    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;
       }

• Intersection cell and objects

    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;
       }

Calling the Constructor

Once the configure file is done, we are able to use Contact Search Algorithm, 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 bins. Maybe is required some redefinition of the typedef defined above in the configure file.

static const std::size_t space_dim = 2; typedef DiscreteConfigure<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 ContainerContactType ContainerContactPair; typedef IteratorContactType IteratorContainerContactPair; typedef PointerContactType PointerContainerContactPair;

After that we call the constructor as follow:

• Using Dynamic Bins

   IteratorType it_begin     =   mBoundaryElements.begin();
   IteratorType it_end       =   mBoundaryElements.end();  
   BinsObjectDynamic<Configure>  rBinsObjectDynamic(it_begin, it_end);  

• Using Static Bins

   BinsObjectStatic<Configure>  Bins(it_begin, it_end);
   const std::size_t MaxNumberOfResults  = 1000;
   std::size_t  NumberOfResults          = 0;
   ResultIteratorType  begin;

   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 more fast
   because the we know size of containers.   

• Lastly, we can get the list of pairs of contacts

          rBins.SearchContact(mPairContacts);

• However, is more fast to do a loop and get the bodies that are in contact with a particular element.

    For(IteratorType it =it_begin; it!=it_end; it++) 
       rBinsObjectDynamic.SearchObjectsInner(*it, Result);  

where Result is the list of contact for a particular element.

References

•  Munjiza. Ante, The Combined Finite-Discrete Element Method, John Wiley & Sons, Ltd. 2004