generate_new_elements_mesher_process.hpp

Go to the documentation of this file.

//
//   Project Name:        KratosDelaunayMeshingApplication $
//   Created by:          $Author:             JMCarbonell $
//   Last modified by:    $Co-Author:                      $
//   Date:                $Date:                April 2018 $
//   Revision:            $Revision:                   0.0 $
//
//
 
#if !defined(KRATOS_GENERATE_NEW_ELEMENTS_MESHER_PROCESS_H_INCLUDED)
#define KRATOS_GENERATE_NEW_ELEMENTS_MESHER_PROCESS_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "custom_meshers/laplacian_smoothing.hpp"
#include "custom_processes/mesher_process.hpp"
 
//Data:     MASTER_ELEMENTS(set), MASTER_NODES(set)
//StepData:
//Flags:    (checked) TO_ERASE, TO_REFINE, CONTACT, NEW_ENTITY
//          (set)     BOUNDARY(set),  [TO_REFINE(nodes), TO_ERASE(condition)]->locally to not preserve condition
//          (modified)
//          (reset)
// (set):=(set in this process)
 
namespace Kratos
{
 
 
 
 
 
 
 
  class GenerateNewElementsMesherProcess
      : public MesherProcess
  {
  public:
 
    KRATOS_CLASS_POINTER_DEFINITION(GenerateNewElementsMesherProcess);
 
    typedef ModelPart::NodeType NodeType;
    typedef ModelPart::ConditionType ConditionType;
    typedef ModelPart::PropertiesType PropertiesType;
    typedef ConditionType::GeometryType GeometryType;
 
    typedef ModelPart::ElementsContainerType ElementsContainerType;
    typedef GlobalPointersVector<Element> ElementWeakPtrVectorType;
 
 
    GenerateNewElementsMesherProcess(ModelPart &rModelPart,
                                     MesherUtilities::MeshingParameters &rRemeshingParameters,
                                     int EchoLevel)
        : mrModelPart(rModelPart),
          mrRemesh(rRemeshingParameters)
    {
      mEchoLevel = EchoLevel;
    }
 
    virtual ~GenerateNewElementsMesherProcess()
    {
    }
 
 
    void operator()()
    {
      Execute();
    }
 
 
    void Execute() override
    {
      KRATOS_TRY
 
      if (mEchoLevel > 0)
        std::cout << " [ GENERATE NEW ELEMENTS: " << std::endl;
 
      if (mrModelPart.Name() != mrRemesh.SubModelPartName)
        std::cout << " ModelPart Supplied do not corresponds to the Meshing Domain: (" << mrModelPart.Name() << " != " << mrRemesh.SubModelPartName << ")" << std::endl;
 
      //*******************************************************************
      //selecting elements
 
      if (!mrRemesh.MeshElementsSelectedFlag) //Select Mesh Elements not performed  ... is needed to be done before building new elements
      {
        std::cout << " ERROR : no selection of elements performed before building the elements " << std::endl;
        SelectElementsMesherProcess SelectElements(mrModelPart, mrRemesh, mEchoLevel);
        SelectElements.Execute();
      }
 
      //*******************************************************************
      //setting new elements
      //(mrModelPart.Elements()).reserve(mrRemesh.Info->NumberOfElements);
 
      //*******************************************************************
      // mine 2016 TIP
      // //All nodes in boundary element change
      // if(mrRemesh.AvoidTipElementsFlag){ //is not working correctly some dispositions not considered
      //   if( mEchoLevel > 0 )
      //    std::cout<<"[   AVOID TIP ELEMENTS START ]"<<std::endl;
 
      //   ChangeTipElementsUtilities TipElements;
      //   //TipElements.SwapDiagonals(mrModelPart,out,mrRemesh.PreservedElements);
 
      //   if( mEchoLevel > 0 )
      //    std::cout<<"[   AVOID TIP ELEMENTS END ]"<<std::endl;
      // }
      //*******************************************************************
 
      //properties to be used in the generation
      ModelPart::ElementsContainerType::iterator element_begin = mrModelPart.ElementsBegin();
      ModelPart::NodesContainerType::iterator nodes_begin = mrModelPart.NodesBegin();
 
      // int number_properties = mrModelPart.GetParentModelPart().NumberOfProperties();
      // if(number_properties<0)
      //   KRATOS_ERROR<<" number of properties is "<<number_properties<<std::endl;
      // Properties::Pointer properties = mrModelPart.GetParentModelPart().GetMesh().pGetProperties(number_properties-1);
      // properties->PrintData(std::cout);
      // std::cout<<std::endl;
 
      MeshDataTransferUtilities DataTransferUtilities;
 
      const Element &rReferenceElement = mrRemesh.GetReferenceElement();
 
      const unsigned int nds = element_begin->GetGeometry().size();
 
      int &OutNumberOfElements = mrRemesh.OutMesh.GetNumberOfElements();
      int *OutElementList = mrRemesh.OutMesh.GetElementList();
 
      std::vector<NodeType::Pointer> list_of_element_centers;
      std::vector<Geometry<NodeType>> list_of_element_vertices; //is this list needed?
      //find the center and "radius" of the element
      double xc = 0;
      double yc = 0;
      double zc = 0;
      double radius = 0;
 
      //generate kratos elements (conditions are not touched)
      int id = 0;
      // std::vector<std::vector<int> > EmptyNeighList;
      // mrRemesh.NeighbourList.swap(EmptyNeighList);
      // mrRemesh.NeighbourList.clear(); //destroy all elements
 
      for (int el = 0; el < OutNumberOfElements; ++el)
      {
        if (mrRemesh.PreservedElements[el])
        {
          Geometry<NodeType> vertices;
          std::vector<int> neighbours(nds);
 
          for (unsigned int i = 0; i < nds; ++i)
          {
            //note that OutElementList, starts from node 1, not from node 0, it can be directly assigned to mrRemesh.NodalPreIds.
            vertices.push_back(*(nodes_begin + OutElementList[el * nds + i] - 1).base());
            //vertices.push_back(mrModelPart.pGetNode(OutElementList[el*3+pn]));
 
            if (vertices.back().Is(TO_ERASE))
              std::cout << " WARNING:: mesh vertex RELEASED " << vertices.back().Id() << std::endl;
          }
 
          id += 1;
 
          mrRemesh.PreservedElements[el] = id;
 
          //*******************************************************************
          //1) Store Preserved elements in an array of vertices (Geometry<NodeType > vertices;)
 
          if (vertices.size() == 3)
          {
            DataTransferUtilities.CalculateCenterAndSearchRadius(vertices[0].X(), vertices[0].Y(),
                                                                 vertices[1].X(), vertices[1].Y(),
                                                                 vertices[2].X(), vertices[2].Y(),
                                                                 xc, yc, radius);
          }
          else if (vertices.size() == 4)
          {
            DataTransferUtilities.CalculateCenterAndSearchRadius(vertices[0].X(), vertices[0].Y(), vertices[0].Z(),
                                                                 vertices[1].X(), vertices[1].Y(), vertices[1].Z(),
                                                                 vertices[2].X(), vertices[2].Y(), vertices[2].Z(),
                                                                 vertices[3].X(), vertices[3].Y(), vertices[3].Z(),
                                                                 xc, yc, zc, radius);
          }
 
          //std::cout<<" XC ["<<id<<"]: ("<<xc<<" "<<yc<<") "<<std::endl;
          //std::cout<<" vertices "<<vertices[0].X()<<" "<<vertices[2].X()<<std::endl;
          //*******************************************************************
 
          NodeType::Pointer p_center = Kratos::make_intrusive<NodeType>(id, xc, yc, zc);
 
          //*******************************************************************
          //2) Create list_of_centers
 
          list_of_element_centers.push_back(p_center);
          list_of_element_vertices.push_back(vertices);
 
          //*******************************************************************
 
          // std::cout<<" list of centers "<<list_of_element_centers.back()->X()<<" "<<list_of_element_centers.back()->Y()<<std::endl;
          // std::cout<<" list of vertices ";
          // std::cout.flush();
          // std::cout<<" vertices "<<list_of_element_vertices.back()[0].X()<<" "<<list_of_element_vertices.back()[2].X()<<std::endl;
          // std::cout.flush();
        }
        else
        {
          mrRemesh.PreservedElements[el] = -1;
        }
 
        // if(mrRemesh.PreservedElements[el])
        //   {
        //  std::cout<<" Neighbours ["<<id-1<<"] :";
        //  for(int pn=0; pn<3; ++pn)
        //    {
        //      std::cout<<" neighborlist ("<<pn<<") "<<mrRemesh.NeighbourList[id-1][pn]<<std::endl;
        //    }
        //   }
 
        //std::cout<<" NodalPreIds ["<<el<<"] :"<<NodalPreIds[el]<<std::endl;
      }
 
      //*******************************************************************
      //5) Laplacian Smoothing
      // check geometrical smoothing is better for solid cutting problems commented July 2018 for fluid testing
      //if( mrRemesh.Options.Is(MesherUtilities::MESH_SMOOTHING) && mrRemesh.Info->GeometricalSmoothingRequired ){
      if (mrRemesh.Options.Is(MesherUtilities::MESH_SMOOTHING))
      {
        //Check Mesh Info to perform smoothing:
        if (mrRemesh.Options.Is(MesherUtilities::REFINE))
          mrRemesh.Info->CheckGeometricalSmoothing();
        else
          mrRemesh.Info->GeometricalSmoothingRequired = true;
 
        int &OutNumberOfPoints = mrRemesh.OutMesh.GetNumberOfPoints();
        LaplacianSmoothing MeshGeometricSmoothing(mrModelPart);
        MeshGeometricSmoothing.SetEchoLevel(mEchoLevel);
        MeshGeometricSmoothing.ApplyMeshSmoothing(mrModelPart, mrRemesh.PreservedElements, OutElementList, OutNumberOfPoints);
      }
      //*******************************************************************
 
      //*******************************************************************
      //6) Pass  rReferenceElement and transfer variables
      DataTransferUtilities.TransferData(mrModelPart, rReferenceElement, list_of_element_centers, list_of_element_vertices, MeshDataTransferUtilities::ELEMENT_TO_ELEMENT);
      //*******************************************************************
 
      //*******************************************************************w
      //std::cout<<" Number of Nodes "<<mrModelPart.Nodes().size()<<" Number Of Ids "<<mrRemesh.NodalPreIds.size()<<std::endl;
 
      //7) Restore global ID's
      for (ModelPart::NodesContainerType::iterator in = mrModelPart.NodesBegin(); in != mrModelPart.NodesEnd(); ++in)
      {
        //std::cout<<" node (local:"<<in->Id()<<", global:"<<mrRemesh.NodalPreIds[ in->Id() ]<<")"<<std::endl;
        in->SetId(mrRemesh.NodalPreIds[in->Id()]);
      }
      //*******************************************************************
 
      //*******************************************************************
 
      //8) Filling the neighbour list
      SetElementNeighbours(mrModelPart);
 
      //*******************************************************************
 
      if (mEchoLevel > 0)
        std::cout << "   GENERATE NEW ELEMENTS ]; " << std::endl;
 
      KRATOS_CATCH(" ")
    }
 
 
 
 
    std::string Info() const override
    {
      return "GenerateNewElementsMesherProcess";
    }
 
    void PrintInfo(std::ostream &rOStream) const override
    {
      rOStream << "GenerateNewElementsMesherProcess";
    }
 
    void PrintData(std::ostream &rOStream) const override
    {
    }
 
 
 
  protected:
 
 
 
 
 
 
 
 
  private:
 
    ModelPart &mrModelPart;
 
    MesherUtilities::MeshingParameters &mrRemesh;
 
    int mEchoLevel;
 
 
    //**************************************************************************
    //**************************************************************************
 
    void SetElementNeighbours(ModelPart &rModelPart)
    {
      KRATOS_TRY
 
      if (mEchoLevel > 0)
      {
        std::cout << " [ SET ELEMENT NEIGHBOURS : " << std::endl;
        std::cout << "   Initial Faces : " << rModelPart.Conditions().size() << std::endl;
      }
 
      ModelPart::ElementsContainerType::iterator element_begin = rModelPart.ElementsBegin();
 
      const unsigned int nds = element_begin->GetGeometry().size();
 
      int *OutElementNeighbourList = mrRemesh.OutMesh.GetElementNeighbourList();
 
      int facecounter = 0;
      int Id = 0;
      for (auto i_elem(rModelPart.ElementsBegin()); i_elem != rModelPart.ElementsEnd(); ++i_elem)
      {
 
        for (unsigned int i = 0; i < mrRemesh.PreservedElements.size(); ++i)
        {
          if (mrRemesh.PreservedElements[Id] == -1)
            Id++;
          else
            break;
        }
 
        //Id = i_elem->Id()-1;
 
        unsigned int number_of_faces = i_elem->GetGeometry().FacesNumber(); //defined for triangles and tetrahedra
 
        ElementWeakPtrVectorType &nElements = i_elem->GetValue(NEIGHBOUR_ELEMENTS);
        nElements.resize(number_of_faces);
 
        int index = 0;
        for (unsigned int iface = 0; iface < number_of_faces; ++iface)
        {
 
          index = OutElementNeighbourList[Id * nds + iface];
 
          if (index > 0)
          {
            //std::cout<<" Element "<<Id<<" size "<<mrRemesh.PreservedElements.size()<<std::endl;
            //std::cout<<" Index pre "<<index<<" size "<<mrRemesh.PreservedElements.size()<<std::endl;
            index = mrRemesh.PreservedElements[index - 1];
            //std::cout<<" Index post "<<index<<std::endl;
          }
 
          if (index > 0)
          {
            nElements(iface) = *(element_begin + index - 1).base();
          }
          else
          {
            nElements(iface) = *i_elem.base();
            facecounter++;
          }
        }
 
        Id++;
      }
 
      if (mEchoLevel > 0)
      {
        std::cout << "   Final Faces : " << facecounter << std::endl;
        std::cout << "   SET ELEMENT NEIGHBORS ]; " << std::endl;
      }
 
      KRATOS_CATCH("")
    }
 
 
 
 
 
    GenerateNewElementsMesherProcess &operator=(GenerateNewElementsMesherProcess const &rOther);
 
    //GenerateNewElementsMesherProcess(GenerateNewElementsMesherProcess const& rOther);
 
 
  }; // Class GenerateNewElementsMesherProcess
 
 
 
 
  inline std::istream &operator>>(std::istream &rIStream,
                                  GenerateNewElementsMesherProcess &rThis);
 
  inline std::ostream &operator<<(std::ostream &rOStream,
                                  const GenerateNewElementsMesherProcess &rThis)
  {
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
  }
 
} // namespace Kratos.
 
#endif // KRATOS_GENERATE_NEW_ELEMENTS_MESHER_PROCESS_H_INCLUDED  defined