mesh.hpp

Go to the documentation of this file.

// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.

#ifndef MFEM_MESH
#define MFEM_MESH

#include "../config/config.hpp"
#include "../general/stable3d.hpp"
#include "triangle.hpp"
#include "tetrahedron.hpp"
#include "vertex.hpp"
#include "ncmesh.hpp"
#include "../fem/eltrans.hpp"
#include "../fem/coefficient.hpp"
#include <iostream>
#include <fstream>
#include <limits>

namespace mfem
{

// Data type mesh

class KnotVector;
class NURBSExtension;
class FiniteElementSpace;
class GridFunction;
struct Refinement;
class named_ifstream;

#ifdef MFEM_USE_MPI
class ParMesh;
class ParNCMesh;
#endif


class Mesh
{
#ifdef MFEM_USE_MPI
   friend class ParMesh;
   friend class ParNCMesh;
#endif
   friend class NURBSExtension;

protected:
   int Dim;
   int spaceDim;

   int NumOfVertices, NumOfElements, NumOfBdrElements;
   int NumOfEdges, NumOfFaces;

   int BaseGeom, BaseBdrGeom; // element base geometries, -1 if not all the same

   int meshgen; // see MeshGenerator()

   // Counter for Mesh transformations: refinement, derefinement, rebalancing.
   // Used for checking during Update operations on objects depending on the
   // Mesh, such as FiniteElementSpace, GridFunction, etc.
   long sequence;

   Array<Element *> elements;
   // Vertices are only at the corners of elements, where you would expect them
   // in the lowest-order mesh.
   Array<Vertex> vertices;
   Array<Element *> boundary;
   Array<Element *> faces;

   struct FaceInfo
   {
      // Inf = 64 * LocalFaceIndex + FaceOrientation
      int Elem1No, Elem2No, Elem1Inf, Elem2Inf;
      int NCFace; /* -1 if this is a regular conforming/boundary face;
                     index into 'nc_faces_info' if >= 0. */
   };
   // NOTE: in NC meshes, master faces have Elem2No == -1. Slave faces on the
   // other hand have Elem2No and Elem2Inf set to the master face's element and
   // its local face number.

   struct NCFaceInfo
   {
      bool Slave; // true if this is a slave face, false if master face
      int MasterFace; // if Slave, this is the index of the master face
      const DenseMatrix* PointMatrix; // if Slave, position within master face
      // (NOTE: PointMatrix points to a matrix owned by NCMesh.)

      NCFaceInfo(bool slave, int master, const DenseMatrix* pm)
         : Slave(slave), MasterFace(master), PointMatrix(pm) {}
   };

   Array<FaceInfo> faces_info;
   Array<NCFaceInfo> nc_faces_info;

   Table *el_to_edge;
   Table *el_to_face;
   Table *el_to_el;
   Array<int> be_to_edge;  // for 2D
   Table *bel_to_edge;     // for 3D
   Array<int> be_to_face;
   mutable Table *face_edge;
   mutable Table *edge_vertex;

   IsoparametricTransformation Transformation, Transformation2;
   IsoparametricTransformation FaceTransformation, EdgeTransformation;
   FaceElementTransformations FaceElemTr;

   // refinement embeddings for forward compatibility with NCMesh
   CoarseFineTransformations CoarseFineTr;

   // Nodes are only active for higher order meshes, and share locations with
   // the vertices, plus all the higher- order control points within the
   // element and along the edges and on the faces.
   GridFunction *Nodes;
   int own_nodes;

   static const int vtk_quadratic_tet[10];
   static const int vtk_quadratic_hex[27];

#ifdef MFEM_USE_MEMALLOC
   friend class Tetrahedron;
   MemAlloc <Tetrahedron, 1024> TetMemory;
#endif

public:
   typedef Geometry::Constants<Geometry::SEGMENT>     seg_t;
   typedef Geometry::Constants<Geometry::TRIANGLE>    tri_t;
   typedef Geometry::Constants<Geometry::SQUARE>      quad_t;
   typedef Geometry::Constants<Geometry::TETRAHEDRON> tet_t;
   typedef Geometry::Constants<Geometry::CUBE>        hex_t;

   enum Operation { NONE, REFINE, DEREFINE, REBALANCE };

   Array<int> attributes;
   Array<int> bdr_attributes;

   NURBSExtension *NURBSext;
   NCMesh *ncmesh;

protected:
   Operation last_operation;

   void Init();

   void InitTables();

   void DeleteTables();

   Element *ReadElementWithoutAttr(std::istream &);
   static void PrintElementWithoutAttr(const Element *, std::ostream &);

   Element *ReadElement(std::istream &);
   static void PrintElement(const Element *, std::ostream &);

   // Readers for different mesh formats, used in the Load() method.
   // The implementations of these methods are in mesh_readers.cpp.
   void ReadMFEMMesh(std::istream &input, bool mfem_v11, int &curved);
   void ReadLineMesh(std::istream &input);
   void ReadNetgen2DMesh(std::istream &input, int &curved);
   void ReadNetgen3DMesh(std::istream &input);
   void ReadTrueGridMesh(std::istream &input);
   void ReadVTKMesh(std::istream &input, int &curved, int &read_gf);
   void ReadNURBSMesh(std::istream &input, int &curved, int &read_gf);
   void ReadInlineMesh(std::istream &input, int generate_edges = 0);
   void ReadGmshMesh(std::istream &input);
   /* Note NetCDF (optional library) is used for reading cubit files */
#ifdef MFEM_USE_NETCDF
   void ReadCubit(named_ifstream &input, int &curved, int &read_gf);
#endif

   static void skip_comment_lines(std::istream &is, const char comment_char)
   {
      while (1)
      {
         is >> std::ws;
         if (is.peek() != comment_char) { break; }
         is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
      }
   }
   // Check for, and remove, a trailing '\r'.
   static void filter_dos(std::string &line)
   {
      if (!line.empty() && *line.rbegin() == '\r')
      { line.resize(line.size()-1); }
   }

   void SetMeshGen(); // set 'meshgen'

   double GetLength(int i, int j) const;

   void GetElementJacobian(int i, DenseMatrix &J);

   void MarkForRefinement();
   void MarkTriMeshForRefinement();
   void GetEdgeOrdering(DSTable &v_to_v, Array<int> &order);
   void MarkTetMeshForRefinement();

   void PrepareNodeReorder(DSTable **old_v_to_v, Table **old_elem_vert);
   void DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert);

   STable3D *GetFacesTable();
   STable3D *GetElementToFaceTable(int ret_ftbl = 0);

   void RedRefinement(int i, const DSTable &v_to_v,
                      int *edge1, int *edge2, int *middle)
   { UniformRefinement(i, v_to_v, edge1, edge2, middle); }

   void GreenRefinement(int i, const DSTable &v_to_v,
                        int *edge1, int *edge2, int *middle)
   { Bisection(i, v_to_v, edge1, edge2, middle); }

   void Bisection(int i, const DSTable &, int *, int *, int *);

   void Bisection(int i, const DSTable &, int *);

   void UniformRefinement(int i, const DSTable &, int *, int *, int *);

   void AverageVertices (int * indexes, int n, int result);

   void InitRefinementTransforms();
   int FindCoarseElement(int i);

   void UpdateNodes();

   virtual void QuadUniformRefinement();

   virtual void HexUniformRefinement();

   virtual void NURBSUniformRefinement();

   virtual void LocalRefinement(const Array<int> &marked_el, int type = 3);

   virtual void NonconformingRefinement(const Array<Refinement> &refinements,
                                        int nc_limit = 0);

   virtual bool NonconformingDerefinement(Array<double> &elem_error,
                                          double threshold, int nc_limit = 0,
                                          int op = 1);

   void DerefineMesh(const Array<int> &derefinements);

   void LoadPatchTopo(std::istream &input, Array<int> &edge_to_knot);

   void UpdateNURBS();

   void PrintTopo(std::ostream &out, const Array<int> &e_to_k) const;

   void GetLocalPtToSegTransformation(IsoparametricTransformation &, int);
   void GetLocalSegToTriTransformation (IsoparametricTransformation &loc,
                                        int i);
   void GetLocalSegToQuadTransformation (IsoparametricTransformation &loc,
                                         int i);
   void GetLocalTriToTetTransformation (IsoparametricTransformation &loc,
                                        int i);
   void GetLocalQuadToHexTransformation (IsoparametricTransformation &loc,
                                         int i);
   void GetLocalFaceTransformation(int face_type, int elem_type,
                                   IsoparametricTransformation &Transf,
                                   int inf);
   void ApplyLocalSlaveTransformation(IsoparametricTransformation &transf,
                                      const FaceInfo &fi);
   bool IsSlaveFace(const FaceInfo &fi) const;

   static int GetTriOrientation (const int * base, const int * test);
   static int GetQuadOrientation (const int * base, const int * test);

   static void GetElementArrayEdgeTable(const Array<Element*> &elem_array,
                                        const DSTable &v_to_v,
                                        Table &el_to_edge);

   void GetVertexToVertexTable(DSTable &) const;

   int GetElementToEdgeTable(Table &, Array<int> &);

   void AddPointFaceElement(int lf, int gf, int el);

   void AddSegmentFaceElement (int lf, int gf, int el, int v0, int v1);

   void AddTriangleFaceElement (int lf, int gf, int el,
                                int v0, int v1, int v2);

   void AddQuadFaceElement (int lf, int gf, int el,
                            int v0, int v1, int v2, int v3);
   bool FaceIsTrueInterior(int FaceNo) const
   {
      return FaceIsInterior(FaceNo) || (faces_info[FaceNo].Elem2Inf >= 0);
   }

   // shift cyclically 3 integers left-to-right
   inline static void ShiftL2R(int &, int &, int &);
   // shift cyclically 3 integers so that the smallest is first
   inline static void Rotate3(int &, int &, int &);

   void FreeElement (Element *E);

   void GenerateFaces();
   void GenerateNCFaceInfo();

   void InitMesh(int _Dim, int _spaceDim, int NVert, int NElem, int NBdrElem);

   void InitBaseGeom();

   void Make3D(int nx, int ny, int nz, Element::Type type, int generate_edges,
               double sx, double sy, double sz);

   void Make2D(int nx, int ny, Element::Type type, int generate_edges,
               double sx, double sy);

   void Make1D(int n, double sx = 1.0);

   void InitFromNCMesh(const NCMesh &ncmesh);

   Mesh(const NCMesh &ncmesh);

   void Swap(Mesh& other, bool non_geometry = false);

public:

   Mesh() { Init(); InitTables(); meshgen = 0; Dim = 0; }

   explicit Mesh(const Mesh &mesh, bool copy_nodes = true);

   Mesh(int _Dim, int NVert, int NElem, int NBdrElem = 0, int _spaceDim= -1)
   {
      if (_spaceDim == -1)
      {
         _spaceDim = _Dim;
      }
      InitMesh(_Dim, _spaceDim, NVert, NElem, NBdrElem);
   }

   Element *NewElement(int geom);

   void AddVertex(const double *);
   void AddTri(const int *vi, int attr = 1);
   void AddTriangle(const int *vi, int attr = 1);
   void AddQuad(const int *vi, int attr = 1);
   void AddTet(const int *vi, int attr = 1);
   void AddHex(const int *vi, int attr = 1);
   void AddHexAsTets(const int *vi, int attr = 1);
   // 'elem' should be allocated using the NewElement method
   void AddElement(Element *elem)     { elements[NumOfElements++] = elem; }
   void AddBdrElement(Element *elem)  { boundary[NumOfBdrElements++] = elem; }
   void AddBdrSegment(const int *vi, int attr = 1);
   void AddBdrTriangle(const int *vi, int attr = 1);
   void AddBdrQuad(const int *vi, int attr = 1);
   void AddBdrQuadAsTriangles(const int *vi, int attr = 1);
   void GenerateBoundaryElements();
   void FinalizeTriMesh(int generate_edges = 0, int refine = 0,
                        bool fix_orientation = true);
   void FinalizeQuadMesh(int generate_edges = 0, int refine = 0,
                         bool fix_orientation = true);
   void FinalizeTetMesh(int generate_edges = 0, int refine = 0,
                        bool fix_orientation = true);
   void FinalizeHexMesh(int generate_edges = 0, int refine = 0,
                        bool fix_orientation = true);

   void SetAttributes();

#ifdef MFEM_USE_GECKO

   void GetGeckoElementReordering(Array<int> &ordering);
#endif

   void ReorderElements(const Array<int> &ordering, bool reorder_vertices = true);

   Mesh(int nx, int ny, int nz, Element::Type type, int generate_edges = 0,
        double sx = 1.0, double sy = 1.0, double sz = 1.0)
   {
      Make3D(nx, ny, nz, type, generate_edges, sx, sy, sz);
   }

   Mesh(int nx, int ny, Element::Type type, int generate_edges = 0,
        double sx = 1.0, double sy = 1.0)
   {
      Make2D(nx, ny, type, generate_edges, sx, sy);
   }

   explicit Mesh(int n, double sx = 1.0)
   {
      Make1D(n, sx);
   }

   Mesh(const char *filename, int generate_edges = 0, int refine = 1,
        bool fix_orientation = true);

   Mesh(std::istream &input, int generate_edges = 0, int refine = 1,
        bool fix_orientation = true);

   Mesh(Mesh *mesh_array[], int num_pieces);

   /* This is similar to the above mesh constructor, but here the current
      mesh is destroyed and another one created based on the data stream
      again given in MFEM, netgen, or VTK format. If generate_edges = 0
      (default) edges are not generated, if 1 edges are generated. */
   void Load(std::istream &input, int generate_edges = 0, int refine = 1,
             bool fix_orientation = true);

   inline int MeshGenerator() { return meshgen; }

   inline int GetNV() const { return NumOfVertices; }

   inline int GetNE() const { return NumOfElements; }

   inline int GetNBE() const { return NumOfBdrElements; }

   inline int GetNEdges() const { return NumOfEdges; }

   inline int GetNFaces() const { return NumOfFaces; }

   int GetNumFaces() const;

   virtual long ReduceInt(int value) const { return value; }

   long GetGlobalNE() const { return ReduceInt(NumOfElements); }

   inline int EulerNumber() const
   { return NumOfVertices - NumOfEdges + NumOfFaces - NumOfElements; }
   inline int EulerNumber2D() const
   { return NumOfVertices - NumOfEdges + NumOfElements; }

   int Dimension() const { return Dim; }
   int SpaceDimension() const { return spaceDim; }

   const double *GetVertex(int i) const { return vertices[i](); }

   double *GetVertex(int i) { return vertices[i](); }

   const Element* const *GetElementsArray() const
   { return elements.GetData(); }

   const Element *GetElement(int i) const { return elements[i]; }

   Element *GetElement(int i) { return elements[i]; }

   const Element *GetBdrElement(int i) const { return boundary[i]; }

   Element *GetBdrElement(int i) { return boundary[i]; }

   const Element *GetFace(int i) const { return faces[i]; }

   int GetFaceBaseGeometry(int i) const;

   int GetElementBaseGeometry(int i = 0) const
   { return i < GetNE() ? elements[i]->GetGeometryType() : BaseGeom; }

   int GetBdrElementBaseGeometry(int i = 0) const
   { return i < GetNBE() ? boundary[i]->GetGeometryType() : BaseBdrGeom; }

   void GetElementVertices(int i, Array<int> &dofs) const
   { elements[i]->GetVertices(dofs); }

   void GetBdrElementVertices(int i, Array<int> &dofs) const
   { boundary[i]->GetVertices(dofs); }

   void GetElementEdges(int i, Array<int> &edges, Array<int> &cor) const;

   void GetBdrElementEdges(int i, Array<int> &edges, Array<int> &cor) const;

   void GetFaceEdges(int i, Array<int> &, Array<int> &) const;

   void GetFaceVertices(int i, Array<int> &vert) const
   {
      if (Dim == 1)
      {
         vert.SetSize(1); vert[0] = i;
      }
      else
      {
         faces[i]->GetVertices(vert);
      }
   }

   void GetEdgeVertices(int i, Array<int> &vert) const;

   Table *GetFaceEdgeTable() const;

   Table *GetEdgeVertexTable() const;

   void GetElementFaces(int i, Array<int> &, Array<int> &) const;

   void GetBdrElementFace(int i, int *, int *) const;

   int GetBdrElementEdgeIndex(int i) const;

   void GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const;

   int GetElementType(int i) const;

   int GetBdrElementType(int i) const;

   /* Return point matrix of element i of dimension Dim X #dofs, where for
      every degree of freedom we give its coordinates in space of dimension
      Dim. */
   void GetPointMatrix(int i, DenseMatrix &pointmat) const;

   /* Return point matrix of boundary element i of dimension Dim X #dofs,
      where for every degree of freedom we give its coordinates in space
      of dimension Dim. */
   void GetBdrPointMatrix(int i, DenseMatrix &pointmat) const;

   static FiniteElement *GetTransformationFEforElementType(int);

   void GetElementTransformation(int i, IsoparametricTransformation *ElTr);

   ElementTransformation *GetElementTransformation(int i);

   void GetElementTransformation(int i, const Vector &nodes,
                                 IsoparametricTransformation *ElTr);

   ElementTransformation * GetBdrElementTransformation(int i);
   void GetBdrElementTransformation(int i, IsoparametricTransformation *ElTr);

   void GetFaceTransformation(int i, IsoparametricTransformation *FTr);

   ElementTransformation *GetFaceTransformation(int FaceNo);

   void GetEdgeTransformation(int i, IsoparametricTransformation *EdTr);

   ElementTransformation *GetEdgeTransformation(int EdgeNo);

   FaceElementTransformations *GetFaceElementTransformations(int FaceNo,
                                                             int mask = 31);

   FaceElementTransformations *GetInteriorFaceTransformations (int FaceNo)
   {
      if (faces_info[FaceNo].Elem2No < 0) { return NULL; }
      return GetFaceElementTransformations (FaceNo);
   }

   FaceElementTransformations *GetBdrFaceTransformations (int BdrElemNo);

   bool FaceIsInterior(int FaceNo) const
   {
      return (faces_info[FaceNo].Elem2No >= 0);
   }
   void GetFaceElements (int Face, int *Elem1, int *Elem2);
   void GetFaceInfos (int Face, int *Inf1, int *Inf2);

   int GetFaceGeometryType(int Face) const;
   int GetFaceElementType(int Face) const;

   void CheckElementOrientation(bool fix_it = true);
   void CheckBdrElementOrientation(bool fix_it = true);

   int GetAttribute(int i) const { return elements[i]->GetAttribute();}

   int GetBdrAttribute(int i) const { return boundary[i]->GetAttribute(); }

   const Table &ElementToElementTable();

   const Table &ElementToFaceTable() const;

   const Table &ElementToEdgeTable() const;

   Table *GetVertexToElementTable();

   Table *GetFaceToElementTable() const;

   virtual void ReorientTetMesh();

   int *CartesianPartitioning(int nxyz[]);
   int *GeneratePartitioning(int nparts, int part_method = 1);
   void CheckPartitioning(int *partitioning);

   void CheckDisplacements(const Vector &displacements, double &tmax);

   // Vertices are only at the corners of elements, where you would expect them
   // in the lowest-order mesh.
   void MoveVertices(const Vector &displacements);
   void GetVertices(Vector &vert_coord) const;
   void SetVertices(const Vector &vert_coord);

   // Nodes are only active for higher order meshes, and share locations with
   // the vertices, plus all the higher- order control points within the element
   // and along the edges and on the faces.
   void GetNode(int i, double *coord);
   void SetNode(int i, const double *coord);

   // Node operations for curved mesh.
   // They call the corresponding '...Vertices' method if the
   // mesh is not curved (i.e. Nodes == NULL).
   void MoveNodes(const Vector &displacements);
   void GetNodes(Vector &node_coord) const;
   void SetNodes(const Vector &node_coord);

   GridFunction *GetNodes() { return Nodes; }
   const GridFunction *GetNodes() const { return Nodes; }
   void NewNodes(GridFunction &nodes, bool make_owner = false);
   void SwapNodes(GridFunction *&nodes, int &own_nodes_);

   void GetNodes(GridFunction &nodes) const;
   void SetNodalFESpace(FiniteElementSpace *nfes);
   void SetNodalGridFunction(GridFunction *nodes, bool make_owner = false);
   const FiniteElementSpace *GetNodalFESpace() const;

   void SetCurvature(int order, bool discont = false, int space_dim = -1,
                     int ordering = 1);

   void UniformRefinement();

   void GeneralRefinement(const Array<Refinement> &refinements,
                          int nonconforming = -1, int nc_limit = 0);

   void GeneralRefinement(const Array<int> &el_to_refine,
                          int nonconforming = -1, int nc_limit = 0);

   void RandomRefinement(double prob, bool aniso = false,
                         int nonconforming = -1, int nc_limit = 0);

   void RefineAtVertex(const Vertex& vert,
                       double eps = 0.0, int nonconforming = -1);

   bool RefineByError(const Array<double> &elem_error, double threshold,
                      int nonconforming = -1, int nc_limit = 0);

   bool RefineByError(const Vector &elem_error, double threshold,
                      int nonconforming = -1, int nc_limit = 0);

   bool DerefineByError(Array<double> &elem_error, double threshold,
                        int nc_limit = 0, int op = 1);

   bool DerefineByError(const Vector &elem_error, double threshold,
                        int nc_limit = 0, int op = 1);

   // NURBS mesh refinement methods
   void KnotInsert(Array<KnotVector *> &kv);
   void DegreeElevate(int t);

   void EnsureNCMesh(bool triangles_nonconforming = false);

   bool Conforming() const { return ncmesh == NULL; }
   bool Nonconforming() const { return ncmesh != NULL; }

   const CoarseFineTransformations &GetRefinementTransforms();

   Operation GetLastOperation() const { return last_operation; }

   long GetSequence() const { return sequence; }

   virtual void PrintXG(std::ostream &out = std::cout) const;

   virtual void Print(std::ostream &out = std::cout) const;

   void PrintVTK(std::ostream &out);

   void PrintVTK(std::ostream &out, int ref, int field_data=0);

   void GetElementColoring(Array<int> &colors, int el0 = 0);

   void PrintWithPartitioning (int *partitioning,
                               std::ostream &out, int elem_attr = 0) const;

   void PrintElementsWithPartitioning (int *partitioning,
                                       std::ostream &out,
                                       int interior_faces = 0);


   void PrintSurfaces(const Table &Aface_face, std::ostream &out) const;

   void ScaleSubdomains (double sf);
   void ScaleElements (double sf);

   void Transform(void (*f)(const Vector&, Vector&));
   void Transform(VectorCoefficient &deformation);

   void RemoveUnusedVertices();

   void RemoveInternalBoundaries();

   double GetElementSize(int i, int type = 0);

   double GetElementSize(int i, const Vector &dir);

   double GetElementVolume(int i);

   void GetBoundingBox(Vector &min, Vector &max, int ref = 2);

   void PrintCharacteristics(Vector *Vh = NULL, Vector *Vk = NULL,
                             std::ostream &out = std::cout);

   virtual void PrintInfo(std::ostream &out = std::cout)
   {
      PrintCharacteristics(NULL, NULL, out);
   }

   void MesquiteSmooth(const int mesquite_option = 0);

   virtual ~Mesh();
};

std::ostream &operator<<(std::ostream &out, const Mesh &mesh);


class NodeExtrudeCoefficient : public VectorCoefficient
{
private:
   int n, layer;
   double p[2], s;
   Vector tip;
public:
   NodeExtrudeCoefficient(const int dim, const int _n, const double _s);
   void SetLayer(const int l) { layer = l; }
   using VectorCoefficient::Eval;
   virtual void Eval(Vector &V, ElementTransformation &T,
                     const IntegrationPoint &ip);
   virtual ~NodeExtrudeCoefficient() { }
};


Mesh *Extrude1D(Mesh *mesh, const int ny, const double sy,
                const bool closed = false);


class named_ifstream : public std::ifstream
{
public:
   const char *filename;
   named_ifstream(const char *mesh_name) :
      std::ifstream(mesh_name), filename(mesh_name) {}
};


// inline functions
inline void Mesh::ShiftL2R(int &a, int &b, int &c)
{
   int t = a;
   a = c;  c = b;  b = t;
}

inline void Mesh::Rotate3(int &a, int &b, int &c)
{
   if (a < b)
   {
      if (a > c)
      {
         ShiftL2R(a, b, c);
      }
   }
   else
   {
      if (b < c)
      {
         ShiftL2R(c, b, a);
      }
      else
      {
         ShiftL2R(a, b, c);
      }
   }
}

}

#endif