build_mesh_boundary_for_fluids_process.hpp

Go to the documentation of this file.

//
//   Project Name:        KratosPfemApplication     $
//   Created by:          $Author:          AFranci $
//   Last modified by:    $Co-Author:               $
//   Date:                $Date:     September 2018 $
//   Revision:            $Revision:            0.0 $
//
//
 
#if !defined(KRATOS_BUILD_MESH_BOUNDARY_FOR_FLUIDS_PROCESS_H_INCLUDED)
#define KRATOS_BUILD_MESH_BOUNDARY_FOR_FLUIDS_PROCESS_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "custom_processes/build_model_part_boundary_process.hpp"
// #include <boost/timer.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
{
 
 
  typedef ModelPart::NodesContainerType NodesContainerType;
  typedef ModelPart::ElementsContainerType ElementsContainerType;
  typedef ModelPart::ConditionsContainerType ConditionsContainerType;
 
  typedef GlobalPointersVector<Node> NodeWeakPtrVectorType;
  typedef GlobalPointersVector<Element> ElementWeakPtrVectorType;
 
 
 
 
  class BuildMeshBoundaryForFluidsProcess
      : public BuildModelPartBoundaryProcess
  {
  public:
 
    KRATOS_CLASS_POINTER_DEFINITION(BuildMeshBoundaryForFluidsProcess);
 
 
    BuildMeshBoundaryForFluidsProcess(ModelPart &rModelPart,
                                      MesherUtilities::MeshingParameters &rRemeshingParameters,
                                      int EchoLevel)
        : BuildModelPartBoundaryProcess(rModelPart, rModelPart.Name(), EchoLevel),
          mrRemesh(rRemeshingParameters)
    {
    }
 
    virtual ~BuildMeshBoundaryForFluidsProcess()
    {
    }
 
 
    void operator()()
    {
      Execute();
    }
 
 
    void Execute() override
    {
      KRATOS_TRY
 
      bool success = false;
 
      if (mEchoLevel > 0)
        std::cout << " [ Build Boundary on ModelPart [" << mrModelPart.Name() << "] ]" << std::endl;
 
      success = UniqueSkinSearch(mrModelPart);
 
      if (!success)
      {
        std::cout << "  ERROR:  BOUNDARY BUILD FAILED ModelPart : [" << mrModelPart << "] " << std::endl;
      }
 
      KRATOS_CATCH(" ")
    }
 
 
 
 
    std::string Info() const override
    {
      return "BuildMeshBoundaryForFluidsProcess";
    }
 
    void PrintInfo(std::ostream &rOStream) const override
    {
      rOStream << "BuildMeshBoundaryForFluidsProcess";
    }
 
    void PrintData(std::ostream &rOStream) const override
    {
    }
 
 
 
  protected:
 
 
 
    //**************************************************************************
    //**************************************************************************
 
    bool UniqueSkinSearch(ModelPart &rModelPart)
    {
 
      KRATOS_TRY
 
      if (mEchoLevel > 0)
      {
        std::cout << " [ Initial Conditions : " << rModelPart.Conditions().size() << std::endl;
      }
 
      if (!rModelPart.Elements().size() || (rModelPart.Is(ACTIVE)))
      {
        if (mEchoLevel > 0)
        {
          std::cout << " [ Final Conditions   : " << rModelPart.Conditions().size() << std::endl;
        }
        return true;
      }
      const unsigned int dimension = mrModelPart.ElementsBegin()->GetGeometry().WorkingSpaceDimension();
 
      // reset the boundary flag in all nodes and check if a remesh process has been performed
      bool any_node_to_erase = false;
      for (ModelPart::NodesContainerType::const_iterator in = rModelPart.NodesBegin(); in != rModelPart.NodesEnd(); in++)
      {
 
        NodeWeakPtrVectorType &nNodes = in->GetValue(NEIGHBOUR_NODES);
        unsigned int neighNodes = nNodes.size();
        if (neighNodes < (dimension + 1))
        {
          in->FastGetSolutionStepValue(ISOLATED_NODE) = 1;
        }
        else
        {
          in->FastGetSolutionStepValue(ISOLATED_NODE) = 0;
        }
 
        if (in->Is(ISOLATED)) // to record the information of the previous step
        {
          in->Set(PFEMFlags::PREVIOUS_ISOLATED, true);
        }
        else
        {
          in->Reset(PFEMFlags::PREVIOUS_ISOLATED);
        }
 
        if (in->Is(FREE_SURFACE)) // to record the information of the previous step
        {
          in->Set(PFEMFlags::PREVIOUS_FREESURFACE, true);
        }
        else
        {
          in->Reset(PFEMFlags::PREVIOUS_FREESURFACE);
        }
 
        in->Reset(BOUNDARY);
        in->Reset(FREE_SURFACE);
 
        if (any_node_to_erase == false)
          if (in->Is(TO_ERASE))
            any_node_to_erase = true;
      }
      SetBoundaryAndFreeSurface(rModelPart);
 
      return true;
 
      KRATOS_CATCH("")
    }
 
    // bool SetBoundaryAndFreeSurface( ModelPart& rModelPart, ModelPart::ConditionsContainerType& rTemporaryConditions, std::vector<int>& rPreservedConditions, unsigned int& rConditionId )
    bool SetBoundaryAndFreeSurface(ModelPart &rModelPart)
    {
 
      KRATOS_TRY
 
      // properties to be used in the generation
      int number_properties = rModelPart.GetParentModelPart().NumberOfProperties();
      Properties::Pointer properties = rModelPart.GetParentModelPart().pGetProperties(number_properties - 1);
 
      ModelPart::ElementsContainerType::iterator elements_begin = mrModelPart.ElementsBegin();
      ModelPart::ElementsContainerType::iterator elements_end = mrModelPart.ElementsEnd();
 
      for (ModelPart::ElementsContainerType::iterator ie = elements_begin; ie != elements_end; ie++)
      {
 
        Geometry<Node> &rElementGeometry = ie->GetGeometry();
 
        if (rElementGeometry.FacesNumber() >= 3)
        { // 3 or 4
 
          /*each face is opposite to the corresponding node number so in 2D triangle
        0 ----- 1 2
        1 ----- 2 0
        2 ----- 0 1
      */
 
          /*each face is opposite to the corresponding node number so in 3D tetrahedron
        0 ----- 1 2 3
        1 ----- 2 0 3
        2 ----- 0 1 3
        3 ----- 0 2 1
      */
 
          // finding boundaries and creating the "skin"
          //
          //********************************************************************
 
          boost::numeric::ublas::matrix<unsigned int> lpofa; // connectivities of points defining faces
          boost::numeric::ublas::vector<unsigned int> lnofa; // number of points defining faces
 
          ElementWeakPtrVectorType &rE = ie->GetValue(NEIGHBOUR_ELEMENTS);
 
          // get matrix nodes in faces
          rElementGeometry.NodesInFaces(lpofa);
          rElementGeometry.NumberNodesInFaces(lnofa);
 
          // loop on neighbour elements of an element
          unsigned int iface = 0;
          for (ElementWeakPtrVectorType::iterator ne = rE.begin(); ne != rE.end(); ne++)
          {
 
            unsigned int NumberNodesInFace = lnofa[iface];
 
            if ((ne)->Id() == ie->Id())
            {
 
              // if no neighbour is present => the face is free surface
              bool freeSurfaceFace = false;
              for (unsigned int j = 1; j <= NumberNodesInFace; j++)
              {
                rElementGeometry[lpofa(j, iface)].Set(BOUNDARY);
                if (rElementGeometry[lpofa(j, iface)].IsNot(RIGID))
                {
                  freeSurfaceFace = true;
                }
              }
              if (freeSurfaceFace == true)
              {
                for (unsigned int j = 1; j <= NumberNodesInFace; j++)
                {
                  rElementGeometry[lpofa(j, iface)].Set(FREE_SURFACE);
                }
              }
 
            } // end face condition
 
            iface += 1;
 
          } // end loop neighbours
        }
      }
 
      return true;
 
      KRATOS_CATCH("")
    }
 
    //**************************************************************************
    //**************************************************************************
 
    bool CheckAcceptedCondition(ModelPart &rModelPart, Condition &rCondition) override
    {
      KRATOS_TRY
 
      bool node_not_preserved = false;
      bool condition_not_preserved = false;
 
      Geometry<Node> &rConditionGeometry = rCondition.GetGeometry();
 
      for (unsigned int j = 0; j < rConditionGeometry.size(); j++)
      {
        if (rConditionGeometry[j].Is(TO_ERASE) || rConditionGeometry[j].Is(TO_REFINE))
          node_not_preserved = true;
 
        if (rConditionGeometry[j].Is(ISOLATED) || rConditionGeometry[j].IsNot(BOUNDARY))
          condition_not_preserved = true;
      }
 
      if (rCondition.Is(TO_ERASE))
        condition_not_preserved = true;
 
      if (rCondition.Is(BOUNDARY)) // flag for composite condition
        condition_not_preserved = true;
 
      if (node_not_preserved == true || condition_not_preserved == true)
        return false;
      else
        return true;
 
      KRATOS_CATCH("")
    }
 
    //**************************************************************************
    //**************************************************************************
 
    void AddConditionToModelPart(ModelPart &rModelPart, Condition::Pointer pCondition) override
    {
      KRATOS_TRY
 
      rModelPart.Conditions().push_back(pCondition);
 
      KRATOS_CATCH("")
    }
 
 
 
 
 
 
  private:
 
 
    MesherUtilities::MeshingParameters &mrRemesh;
 
 
 
    //**************************************************************************
    //**************************************************************************
 
    bool FindNodeInCondition(Geometry<Node> &rConditionGeometry, Geometry<Node> &rElementGeometry, boost::numeric::ublas::matrix<unsigned int> &lpofa, boost::numeric::ublas::vector<unsigned int> &lnofa, unsigned int &iface)
    {
      KRATOS_TRY
 
      // not equivalent geometry sizes for boundary conditions:
      if (rConditionGeometry.size() != lnofa[iface])
        return false;
 
      // line boundary condition:
      if (lnofa[iface] == 2)
      {
        if (rConditionGeometry[0].Id() == rElementGeometry[lpofa(1, iface)].Id() ||
            rConditionGeometry[1].Id() == rElementGeometry[lpofa(2, iface)].Id() ||
            rConditionGeometry[0].Id() == rElementGeometry[lpofa(2, iface)].Id() ||
            rConditionGeometry[1].Id() == rElementGeometry[lpofa(1, iface)].Id())
        {
          return true;
        }
        else
        {
          return false;
        }
      }
 
      // 3D faces:
      if (lnofa[iface] == 3)
      {
        if (rConditionGeometry[0].Id() == rElementGeometry[lpofa(1, iface)].Id() ||
            rConditionGeometry[1].Id() == rElementGeometry[lpofa(2, iface)].Id() ||
            rConditionGeometry[2].Id() == rElementGeometry[lpofa(3, iface)].Id() ||
            rConditionGeometry[0].Id() == rElementGeometry[lpofa(3, iface)].Id() ||
            rConditionGeometry[1].Id() == rElementGeometry[lpofa(1, iface)].Id() ||
            rConditionGeometry[2].Id() == rElementGeometry[lpofa(2, iface)].Id() ||
            rConditionGeometry[0].Id() == rElementGeometry[lpofa(2, iface)].Id() ||
            rConditionGeometry[1].Id() == rElementGeometry[lpofa(3, iface)].Id() ||
            rConditionGeometry[2].Id() == rElementGeometry[lpofa(1, iface)].Id())
        {
          return true;
        }
        else
        {
          return false;
        }
      }
 
      if (lnofa[iface] > 3)
      {
        KRATOS_THROW_ERROR(std::logic_error, "Wrong Condition Number of Face Nodes", *this);
      }
 
      return false;
 
      KRATOS_CATCH(" ")
    }
 
 
 
 
    BuildMeshBoundaryForFluidsProcess &operator=(BuildMeshBoundaryForFluidsProcess const &rOther);
 
 
  }; // Class BuildMeshBoundaryForFluidsProcess
 
 
 
 
  inline std::istream &operator>>(std::istream &rIStream,
                                  BuildMeshBoundaryForFluidsProcess &rThis);
 
  inline std::ostream &operator<<(std::ostream &rOStream,
                                  const BuildMeshBoundaryForFluidsProcess &rThis)
  {
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
  }
 
} // namespace Kratos.
 
#endif // KRATOS_BUILD_MESH_BOUNDARY_FOR_FLUIDS_PROCESS_H_INCLUDED  defined