transfer_utility.h

Go to the documentation of this file.

/*
==============================================================================
KratosIncompressibleFluidApplication
A library based on:
Kratos
A General Purpose Software for Multi-Physics Finite Element Analysis
Version 1.0 (Released on march 05, 2007).
 
Copyright 2007
Pooyan Dadvand, Riccardo Rossi
pooyan@cimne.upc.edu
rrossi@cimne.upc.edu
- CIMNE (International Center for Numerical Methods in Engineering),
Gran Capita' s/n, 08034 Barcelona, Spain
 
 
Permission is hereby granted, free  of charge, to any person obtaining
a  copy  of this  software  and  associated  documentation files  (the
"Software"), to  deal in  the Software without  restriction, including
without limitation  the rights to  use, copy, modify,  merge, publish,
distribute,  sublicense and/or  sell copies  of the  Software,  and to
permit persons to whom the Software  is furnished to do so, subject to
the following condition:
 
Distribution of this code for  any  commercial purpose  is permissible
ONLY BY DIRECT ARRANGEMENT WITH THE COPYRIGHT OWNERS.
 
The  above  copyright  notice  and  this permission  notice  shall  be
included in all copies or substantial portions of the Software.
 
THE  SOFTWARE IS  PROVIDED  "AS  IS", WITHOUT  WARRANTY  OF ANY  KIND,
EXPRESS OR  IMPLIED, INCLUDING  BUT NOT LIMITED  TO THE  WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT  SHALL THE AUTHORS OR COPYRIGHT HOLDERS  BE LIABLE FOR ANY
CLAIM, DAMAGES OR  OTHER LIABILITY, WHETHER IN AN  ACTION OF CONTRACT,
TORT  OR OTHERWISE, ARISING  FROM, OUT  OF OR  IN CONNECTION  WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 
==============================================================================
*/
 
//
//   Project Name:        Kratos
//   Last Modified by:    $Author: pbecker $
//   Date:                $Date: 2011-09-21 12:30:32 $
//   Revision:            $Revision: 1.0 $
//
//
 
 
#if !defined(KRATOS_MOVE_PART_UTILITY_FLUID_ONLY_DIFF2_INCLUDED )
#define  KRATOS_MOVE_PART_UTILITY_FLUID_ONLY_DIFF2_INCLUDED
 
 
 
// System includes
#include <string>
#include <iostream>
#include <algorithm>
 
// External includes
 
 
// Project includes
#include "includes/define.h"
#include "includes/node.h"
 
#include "includes/dof.h"
#include "includes/variables.h"
#include "containers/array_1d.h"
#include "containers/data_value_container.h"
#include "includes/mesh.h"
#include "utilities/math_utils.h"
 
#include "utilities/geometry_utilities.h"
 
#include "includes/model_part.h"
 
 
#include "spatial_containers/spatial_containers.h"
#include "spatial_containers/cell.h"
#include "spatial_containers/bins_dynamic_objects.h"
 
#include "utilities/spatial_containers_configure.h"
 
#include "geometries/line_2d_2.h"
#include "geometries/triangle_2d_3.h"
#include "geometries/triangle_3d_3.h"
#include "geometries/point.h"
 
#include "pfem_2_application_variables.h"
#include "utilities/openmp_utils.h"
 
#include "time.h"
 
//#include "processes/process.h"
 
namespace Kratos
{
    //this class is to be modified by the user to customize the interpolation process
    template< unsigned int TDim>
    class MoveParticleUtilityDiffFluidOnly
    {
    public:
 
        typedef SpatialContainersConfigure<TDim>     Configure;
        typedef typename Configure::PointType                      PointType;
        //typedef PointType::CoordinatesArrayType           CoordinatesArrayType;
        typedef typename Configure::ContainerType                  ContainerType;
        //typedef Configure::PointerType                    PointerType;
        typedef typename Configure::IteratorType                   IteratorType;
        typedef typename Configure::ResultContainerType            ResultContainerType;
        //typedef Configure::ResultPointerType              ResultPointerType;
        typedef typename Configure::ResultIteratorType             ResultIteratorType;
        //typedef Configure::ContactPairType                ContactPairType;
        //typedef Configure::ContainerContactType           ContainerContactType;
        //typedef Configure::IteratorContactType            IteratorContactType;
        //typedef Configure::PointerContactType             PointerContactType;
        //typedef Configure::PointerTypeIterator            PointerTypeIterator;
 
        KRATOS_CLASS_POINTER_DEFINITION(TransferUtility);
 
        //template<unsigned int TDim>
        TransferUtility(ModelPart& calculation_model_part, ModelPart& topographic_model_part)
            : mcalculation_model_part(calculation_model_part) , mtopographic_model_part(topographic_model_part)
        {
            KRATOS_WATCH("initializing transfer utility")
 
            ProcessInfo& CurrentProcessInfo = mcalculation_model_part.GetProcessInfo();
 
 
 
            //loop in elements to change their ID to their position in the array. Easier to get information later.
            //DO NOT PARALELIZE THIS! IT MUST BE SERIAL!!!!!!!!!!!!!!!!!!!!!!
            /*
            ModelPart::ElementsContainerType::iterator ielembegin = mcalculation_model_part.ElementsBegin();
            for(unsigned int ii=0; ii<mr_model_part.Elements().size(); ii++)
            {
                ModelPart::ElementsContainerType::iterator ielem = ielembegin+ii;
                ielem->SetId(ii+1);
            }
            mlast_elem_id= (mr_model_part.ElementsEnd()-1)->Id();
            */
 
 
            //CONSTRUCTING BIN STRUCTURE
            ContainerType& rElements           =  mtopographic_model_part.ElementsArray();
            IteratorType it_begin              =  rElements.begin();
            IteratorType it_end                =  rElements.end();
            //const int number_of_elem         =   rElements.size();
 
            typename BinsObjectDynamic<Configure>::Pointer paux = typename BinsObjectDynamic<Configure>::Pointer(new BinsObjectDynamic<Configure>(it_begin, it_end  ) );
            paux.swap(mpBinsObjectDynamic);
 
        }
 
 
        ~TransferUtility()
        {}
 
 
        void GatherInformationFromTopographicDomain()
        {
            KRATOS_TRY
            KRATOS_WATCH("Gathering Information From Topographic Domain ")
            ProcessInfo& CurrentProcessInfo = mcalculation_model_part.GetProcessInfo();
            double delta_t = CurrentProcessInfo[DELTA_TIME];
            array_1d<double,3> & gravity= CurrentProcessInfo[GRAVITY];
 
            const unsigned int max_results = 1000;
 
            //array_1d<double,TDim+1> N;
            //const int max_results = 1000;
 
            ModelPart::NodesContainerType::iterator inodebegin = mcalculation_model_part.NodesBegin();
 
 
            vector<unsigned int> node_partition;
            #ifdef _OPENMP
                int number_of_threads = omp_get_max_threads();
            #else
                int number_of_threads = 1;
            #endif
            OpenMPUtils::CreatePartition(number_of_threads, mcalculation_model_part.Nodes().size(), node_partition);
 
            //before doing anything we must reset the vector of nodes contained by each element (particles that are inside each element.
            #pragma omp parallel for
            for(int kkk=0; kkk<number_of_threads; kkk++)
            {
                array_1d<double,TDim+1> N;
                ResultContainerType results(max_results);
                ResultIteratorType result_begin = results.begin();
 
                for(unsigned int ii=node_partition[kkk]; ii<node_partition[kkk+1]; ii++)
                {
                    if ( (results.size()) !=max_results)
                        results.resize(max_results);
 
                    //const int & elem_id = ielem->Id();
                    ModelPart::NodesContainerType::iterator inode = inodebegin+ii;
                    Element::Pointer pelement(*ielem.base());
                    Geometry<Node >& geom = ielem->GetGeometry();
 
                    ParticlePointerVector&  element_particle_pointers =  (ielem->GetValue(FLUID_PARTICLE_POINTERS));
                    int & number_of_particles_in_elem=ielem->GetValue(NUMBER_OF_PARTICLES);
                    //std::cout << "elem " << ii << " with " << (unsigned int)number_of_particles_in_elem << " particles" << std::endl;
 
 
                    is_found = FindNodeOnMesh(position, N ,pelement,result_begin,MaxNumberOfResults); //good, now we know where this point is:
 
 
                    }
 
 
 
 
                }
            }
            KRATOS_CATCH("")
        }
 
 
 
 
    protected:
 
 
    private:
 
 
 
 
 
 
    bool FindNodeOnMesh( //int last_element,
                         array_1d<double,3>& position,
                         array_1d<double,TDim+1>& N,
                         //Element::Pointer pelement,
                         Element::Pointer & pelement,
                         ResultIteratorType result_begin,
                         const unsigned int MaxNumberOfResults)
    {
        typedef std::size_t SizeType;
 
        const array_1d<double,3>& coords = position;
         array_1d<double,TDim+1> aux_N;
        //before using the bin to search for possible elements we check first the last element in which the particle was.
 
        //ModelPart::ElementsContainerType::iterator i = mr_model_part.ElementsBegin()+last_element;
        Geometry<Node >& geom_default = pelement->GetGeometry(); //(*(i))->GetGeometry();
        bool is_found_1 = CalculatePosition(geom_default,coords[0],coords[1],coords[2],N);
        if(is_found_1 == true)
        {
            //pelement = (*(i));
            return true;
        }
 
        //KRATOS_WATCH("will look in another element")
        //KRATOS_WATCH(TDim)
 
        //to begin with we check the neighbour elements:
        GlobalPointersVector< Element >& neighb_elems = pelement->GetValue(NEIGHBOUR_ELEMENTS);
        //the first we check is the one that has negative shape function, because it means it went outside in this direction:
        /*
        unsigned int checked_element=0;
        for (unsigned int i=0;i!=(TDim+1);i++)
        {
            if (N[i]<0.0)
            {
                checked_element=i;
                Geometry<Node >& geom = neighb_elems[i].GetGeometry();
                bool is_found_2 = CalculatePosition(geom,coords[0],coords[1],coords[2],aux_N);
                if (is_found_2)
                {
                    pelement=Element::Pointer(((neighb_elems(i))));
                    N=aux_N;
                    return true;
                }
                break;
            }
        }
        */
 
        for (unsigned int i=0;i!=(neighb_elems.size());i++)
        {
            if(neighb_elems(i).get()!=nullptr)
            {
                Geometry<Node >& geom = neighb_elems[i].GetGeometry();
                bool is_found_2 = CalculatePosition(geom,coords[0],coords[1],coords[2],N);
                if (is_found_2)
                {
                    pelement=Element::Pointer(((neighb_elems(i))));
                    return true;
                }
            }
        }
 
 
        //ask to the container for the list of candidate elements
        SizeType results_found = mpBinsObjectDynamic->SearchObjectsInCell(coords, result_begin, MaxNumberOfResults );
        //KRATOS_WATCH(results_found)
 
        if(results_found>0){
            //loop over the candidate elements and check if the particle falls within
            for(SizeType i = 0; i< results_found; i++)
            {
                //std::cout<< "KIIIIIIIIIIIIII" << std::endl;
                //KRATOS_WATCH((*(result_begin+i))->Id());
                Geometry<Node >& geom = (*(result_begin+i))->GetGeometry();
 
 
                //find local position
                bool is_found = CalculatePosition(geom,coords[0],coords[1],coords[2],N);
 
                //KRATOS_WATCH("ln243");
                //KRATOS_WATCH(N);
 
                if(is_found == true)
                {
                    //pelement.clear();
                    //pelement.push_back( Element::WeakPointer((*(result_begin+i).base())));
                    pelement=Element::Pointer((*(result_begin+i).base()));
                    return true;
                }
            }
        }
 
        //not found case
        return false;
    }
 
 
 
    //***************************************
        //***************************************
 
        inline bool CalculatePosition(Geometry<Node >&geom,
                const double xc, const double yc, const double zc,
                array_1d<double, 3 > & N
                )
        {
            double x0 = geom[0].X();
            double y0 = geom[0].Y();
            double x1 = geom[1].X();
            double y1 = geom[1].Y();
            double x2 = geom[2].X();
            double y2 = geom[2].Y();
 
            double area = CalculateVol(x0, y0, x1, y1, x2, y2);
            double inv_area = 0.0;
            if (area == 0.0)
            {
                KRATOS_THROW_ERROR(std::logic_error, "element with zero area found", "");
            } else
            {
                inv_area = 1.0 / area;
            }
 
 
            N[0] = CalculateVol(x1, y1, x2, y2, xc, yc) * inv_area;
            N[1] = CalculateVol(x2, y2, x0, y0, xc, yc) * inv_area;
            N[2] = CalculateVol(x0, y0, x1, y1, xc, yc) * inv_area;
            //KRATOS_WATCH(N);
 
            if (N[0] >= 0.0 && N[1] >= 0.0 && N[2] >= 0.0 && N[0] <= 1.0 && N[1] <= 1.0 && N[2] <= 1.0) //if the xc yc is inside the triangle return true
                return true;
 
            return false;
        }
 
            //***************************************
        //***************************************
 
        inline bool CalculatePosition(Geometry<Node >&geom,
                const double xc, const double yc, const double zc,
                array_1d<double, 4 > & N
                )
        {
 
            double x0 = geom[0].X();
            double y0 = geom[0].Y();
            double z0 = geom[0].Z();
            double x1 = geom[1].X();
            double y1 = geom[1].Y();
            double z1 = geom[1].Z();
            double x2 = geom[2].X();
            double y2 = geom[2].Y();
            double z2 = geom[2].Z();
            double x3 = geom[3].X();
            double y3 = geom[3].Y();
            double z3 = geom[3].Z();
 
            double vol = CalculateVol(x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3);
 
            double inv_vol = 0.0;
            if (vol < 0.0000000000001)
            {
                KRATOS_THROW_ERROR(std::logic_error, "element with zero vol found", "");
            } else
            {
                inv_vol = 1.0 / vol;
            }
 
            N[0] = CalculateVol(x1, y1, z1, x3, y3, z3, x2, y2, z2, xc, yc, zc) * inv_vol;
            N[1] = CalculateVol(x0, y0, z0, x1, y1, z1, x2, y2, z2, xc, yc, zc) * inv_vol;
            N[2] = CalculateVol(x3, y3, z3, x1, y1, z1, x0, y0, z0, xc, yc, zc) * inv_vol;
            N[3] = CalculateVol(x3, y3, z3, x0, y0, z0, x2, y2, z2, xc, yc, zc) * inv_vol;
 
 
            if (N[0] >= 0.0 && N[1] >= 0.0 && N[2] >= 0.0 && N[3] >= 0.0 &&
                    N[0] <= 1.0 && N[1] <= 1.0 && N[2] <= 1.0 && N[3] <= 1.0)
                //if the xc yc zc is inside the tetrahedron return true
                return true;
 
            return false;
        }
 
        inline double CalculateVol(const double x0, const double y0,
                const double x1, const double y1,
                const double x2, const double y2
                )
        {
            return 0.5 * ((x1 - x0)*(y2 - y0)- (y1 - y0)*(x2 - x0));
        }
        //***************************************
        //***************************************
 
        inline double CalculateVol(const double x0, const double y0, const double z0,
                const double x1, const double y1, const double z1,
                const double x2, const double y2, const double z2,
                const double x3, const double y3, const double z3
                )
        {
            double x10 = x1 - x0;
            double y10 = y1 - y0;
            double z10 = z1 - z0;
 
            double x20 = x2 - x0;
            double y20 = y2 - y0;
            double z20 = z2 - z0;
 
            double x30 = x3 - x0;
            double y30 = y3 - y0;
            double z30 = z3 - z0;
 
            double detJ = x10 * y20 * z30 - x10 * y30 * z20 + y10 * z20 * x30 - y10 * x20 * z30 + z10 * x20 * y30 - z10 * y20 * x30;
            return detJ * 0.1666666666666666666667;
        }
 
 
 
    ModelPart& mcalculation_model_part;
    ModelPart& mtopographic_model_part;
 
    typename BinsObjectDynamic<Configure>::Pointer  mpBinsObjectDynamic;
 
    };
 
}  // namespace Kratos.
 
#endif // KRATOS_MOVE_PART_UTILITY_DIFF2_INCLUDED  defined