uibhmTopology.h

Go to the documentation of this file.

/************************************************************************
**
** Copyright (C) 2014 by Carlos Augusto Teixera Mendes
** All rights reserved.
**
** This file is part of the "GeMA" software. It's use should respect
** the terms in the license agreement that can be found together
** with this source code.
** It is provided AS IS, with NO WARRANTY OF ANY KIND,
** INCLUDING THE WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR
** A PARTICULAR PURPOSE.
**
************************************************************************/
 
#ifndef _UIBHM_TOPOLOGY_H_
#define _UIBHM_TOPOLOGY_H_
 
#include <assert.h>
 
#include "gmpGemaMeshConfig.h"
#include "uibhmQueue.h"
 
#include <gmLog.h>
 
class GmCellMesh;
class GmPCR;
template <template <class> class Vector> class GmpNodeMeshData;
template <template <class> class Vector> class GmpCellMeshData;
 
template <class T> class UibhmQuery;
template <class T> class UibhmSurfaceQuery;
template <class T> class UibhmSolidQuery;
 
#define GMP_UIBHM_SIDE_BITS 4
 
#define GM_UIBHM_SIDE_MASK (~(~0u >> GMP_UIBHM_SIDE_BITS))
 
#define GMP_UIBHM_MAX_SIDE_VALUE ((1 << GMP_UIBHM_SIDE_BITS) - 1)
 
#define GMP_UIBHM_BORDER_SIDE  GMP_UIBHM_MAX_SIDE_VALUE
 
#define GMP_UIBHM_MID_SIDE (GMP_UIBHM_MAX_SIDE_VALUE - 1)
 
#define GMP_UIBHM_MAX_NUM_SIDES GMP_UIBHM_MID_SIDE
 
#define GMP_UIBHM_CELL_BITS (32 - GMP_UIBHM_SIDE_BITS)
 
#define GM_UIBHM_CELL_MASK (~GM_UIBHM_SIDE_MASK)
 
#define GMP_UIBHM_MAX_NUM_CELLS (1 << GMP_UIBHM_CELL_BITS)
 
#define GMP_UIBHM_INVALID_VALUE (~0u)
 
#define GMP_UIBHM_INTERFACE_BORDER (GMP_UIBHM_INVALID_VALUE << 1)
 
 
// We don't want the topology structure to be able to store less cells than the mesh 
// representation
static_assert(GMP_UIBHM_CELL_BITS >= GMP_GEMAMESH_CELL_ID_BITS, "Unexpected number of cell id bits");
 
 
class GMP_GEMAMESH_API_EXPORT UibhmTopologyBase
{
public:
 
  enum CheckOptions
  {
    BasicCheck,       
    SimpleCheck,      
    FullCheck,        
    //-----------------------------------------
    // IMPORTANT: The code expects that all Warn options are placed AFTER the error options (above). Please RESPECT that.
    //            It also expects that FirstWarnOption is equal to the first warn option!!!!
    //-----------------------------------------
    FirstWarnOption,
    BasicCheckWarn = FirstWarnOption, 
    SimpleCheckWarn,                  
    FullCheckWarn,                    
  };
 
  virtual ~UibhmTopologyBase() {}
 
  virtual bool init(bool buildEdgeMap, CheckOptions extraChecks) = 0;
  virtual bool addNewCells() = 0;
  virtual bool cellIsIndexed(int cellId) = 0;
};
 
 
template <class Geometry> 
class GMP_GEMAMESH_API_EXPORT UibhmTopology : public UibhmTopologyBase, public Geometry
{
public:
 
  using HalfSide = quint32;
 
  using HSVT = typename Geometry::template VectorType<HalfSide>;
 
  using IVT = typename Geometry::template VectorType<int>;
 
  UibhmTopology(typename Geometry::NodeData* nodeData, typename Geometry::CellData* cellData, 
                const GmPCR* pcr, const GmLogCategory& logger)
    : Geometry(nodeData, cellData, pcr), _nhalf(0), _lastAddedCells(0), _hasOffset(false), _needsEH(false),
      _extraChecks(BasicCheck), _checking(0), _checkLogger(logger) {}
 
  virtual bool init(bool buildEdgeMap, CheckOptions extraChecks);
  virtual bool addNewCells();
 
  virtual bool cellIsIndexed(int cellId) { return cellId < _lastAddedCells; }
 
  void printParameters(const GmLogCategory& logger, GmLogLevel level) const;
  void print          (const GmLogCategory& logger, GmLogLevel level) const;
 
  size_t printMemoryStatistics(const GmLogCategory& logger, GmLogLevel level) const;
  size_t usedMemory() const;
 
  static HalfSide hsMake(int cellId, int side)
  {
    assert(cellId >= 0 && cellId < GMP_UIBHM_MAX_NUM_CELLS);
    assert(side   >= 0 && side   < GMP_UIBHM_MAX_NUM_SIDES);
    return (side << GMP_UIBHM_CELL_BITS) | cellId;
  }
 
  static HalfSide hsMakeSpecial(int id, int side)
  {
    assert(id >= 0 && id < GMP_UIBHM_MAX_NUM_CELLS);
    assert(side == GMP_UIBHM_BORDER_SIDE || side == GMP_UIBHM_MID_SIDE);
    return (side << GMP_UIBHM_CELL_BITS) | id;
  }
 
  static int hsCell(HalfSide hs) { return hs & GM_UIBHM_CELL_MASK; }
 
  static int hsSide(HalfSide hs) { return hs >> GMP_UIBHM_CELL_BITS; }
 
  void hsBreak(HalfSide hs, int* cellId, int* side) const
  {
    *cellId = hsCell(hs);
    *side   = hsSide(hs);
    assert(*cellId >= 0 && *cellId < this->numCells());
    assert(*side   >= 0 && *side   < this->numSides(*cellId));
  }
 
  static bool hsIsValid(HalfSide hs) { return hsSide(hs) < GMP_UIBHM_MAX_NUM_SIDES; }
 
  static bool hsIsMidElement(HalfSide hs) { return hsSide(hs) == GMP_UIBHM_MID_SIDE; }
 
  int linearHs(int cellId, int side) const
  {
    int lhs;
    if constexpr(Geometry::IsHybrid)
      lhs = (_hasOffset ? _Of[cellId] : this->co(cellId)) + side;
    else
      lhs = cellId * this->_typeTemplate->numSides() + side;
    assert(lhs >= 0 && lhs < _nhalf);
    return lhs;
  }
 
  int linearHs(HalfSide hs) const
  {
    int cellId, side;
    hsBreak(hs, &cellId, &side);
    return linearHs(cellId, side);
  }
 
  HalfSide oppositeHs(int cellId, int side) const { return _O[linearHs(cellId, side)]; }
 
  HalfSide oppositeHs(HalfSide hs) const { return _O[linearHs(hs)]; }
 
  bool isBorderHs(int cellId, int side) const { return hsSide(oppositeHs(cellId, side)) == GMP_UIBHM_BORDER_SIDE; }
 
  bool isBorderHs(HalfSide hs) const { assert(hsIsValid(hs));  return hsSide(oppositeHs(hs)) == GMP_UIBHM_BORDER_SIDE; }
 
  bool isNonInterfaceBorderHs(int cellId, int side) const
  {
    HalfSide o = oppositeHs(cellId, side);
    return (hsSide(o) == GMP_UIBHM_BORDER_SIDE) && (o != GMP_UIBHM_INTERFACE_BORDER);
  }
 
  bool isNonInterfaceBorderHs(HalfSide hs) const
  {
    assert(hsIsValid(hs));
    HalfSide o = oppositeHs(hs);
    return (hsSide(o) == GMP_UIBHM_BORDER_SIDE) && (o != GMP_UIBHM_INTERFACE_BORDER);
  }
 
  bool isInterfaceBorderHs(HalfSide hs) const
  {
    assert(hsIsValid(hs));
    return oppositeHs(hs) == GMP_UIBHM_INTERFACE_BORDER;
  }
 
  void fillBorderKind(int cellId, int side, bool* isExtBorder, bool* isInterfaceBorder)
  {
    HalfSide o = oppositeHs(cellId, side);
    *isExtBorder       = (hsSide(o) == GMP_UIBHM_BORDER_SIDE) && (o != GMP_UIBHM_INTERFACE_BORDER);
    *isInterfaceBorder = (o == GMP_UIBHM_INTERFACE_BORDER);
  }
 
  HalfSide firstHs(int cellId) const
  {
    assert(cellId >= 0 && cellId < this->numCells());
    return hsMake(cellId, 0);
  }
 
  HalfSide nextHs(HalfSide hs) const
  {
    int cellId, side;
    hsBreak(hs, &cellId, &side);
    return hsMake(cellId, (side + 1) % this->numSides(cellId)); 
  }
 
  HalfSide prevHs(HalfSide hs) const
  {
    int cellId, side;
    hsBreak(hs, &cellId, &side);
    int ncs = this->numSides(cellId);
    return hsMake(cellId, (side + ncs - 1) % ncs);
  }
 
  HalfSide adjacentHf(HalfSide hf, int faceEdge) const
  {
    assert(this->solid());
    int cellId, face;
    hsBreak(hf, &cellId, &face);
    assert(faceEdge >= 0 && faceEdge < this->tt(cellId)->numFaceVertices(face));
    return hsMake(cellId, this->tt(cellId)->adjacentFaces(face)[faceEdge]);
  }
 
  HalfSide radialHf(HalfSide hf, int faceEdge) const { return oppositeHs(adjacentHf(hf, faceEdge)); }
 
 protected:
  bool doInit();
  bool doAddNewCells();
 
  int  buildOffsetTable();
  int  updateOffsetTable(int numAddedCells);
 
  bool buildOppositeTable();
  bool updateOppositeTable(int numAddedCells);
  template <class Key, class MakeKey> bool oppositeTableTraverseCells(MakeKey makeKey, int first, int last, QMap <Key, HalfSide>& adjMap);
 
  bool buildVertexHalfSideTable();
  bool updateVertexHalfSideTable(int numAddedCells, QVector<HalfSide>& borderToIntBorderList);
  void vertexHalfSideTableTraverseCells(int first, int last, const QSet<int>& oldCells = QSet<int>());
 
  bool buildEdgeHalfFaceMap();
  bool updateEdgeHalfFaceMap(int numAddedCells);
  bool edgeHalfFaceMapTraverseCells(int first, int last);
 
  bool adjustSolidInterfaceElementBorders();
  bool updateSolidInterfaceElementBorders(int numAddedCells, QVector<HalfSide>& borderToIntBorderList);
  bool solidInterfaceElementBordersTraverseCells(int first, int last, QVector<HalfSide>& borderToIntBorderList);
 
  bool buildSurfaceBordersTable();
  bool updateSurfaceBordersTable(int numAddedCells);
 
  bool buildSolidBordersTable();
  bool updateSolidBordersTable(int numAddedCells);
 
  struct BSIterator
  {
    UibhmQueue<HalfSide> _queue;   
    QVector<bool>        _visited; 
  };
 
  BSIterator* borderSurfaceIterator(int borderIndex);
  HalfSide    nextBorderSurface(BSIterator* context);
 
  bool numNodesDifferentFromNumEdges() const;
 
  bool checkLimits(int ncells);
  bool checkOppositeSides(HalfSide s1, HalfSide s2);
  template <class Key> bool checkOppositeMap(const QMap<Key, int>& countMap);
  bool fullNodeCheck();
 
  bool logError(QString msg);
  bool logFindFaceEdgeError(int cellId, int face, int v1, int v2, bool bidir = false);
  bool logFindNodeIndexError(int cellId, int v);
 
  // The auxiliar map used to build the opposite table must be indexed by a key
  // that uniquely identifies a side. Since we want to support large meshes we 
  // can NOT use a simple coding method like just doing v1*(nv^2) + v2*nv + v3
  // for a face identified by v1, v2 & v3 (with a maximum of nv vertices), since 
  // that will result on an overflow, even if storing the result in a 64 bits 
  // integer, with around 2.6 million nodes.
  // - For surface meshes we encode the edge with a size_t as proposed by 
  // Geometry::encodeEdge
  // - For solid meshes without interface elements, we use a vector with 3 numbers
  // filled with the smallest vertices from the face. Using only 3 vertices is possible
  // since we do NOT support non-manifold meshes and so if two elements share a face,
  // both faces should have the exact same vertices (although in different orders) and
  // there can not be a "partial" intersection between faces (other than at an edge or 
  // a vertex).
  // - For solid meshes with interface elements, there is an additional complication 
  // since when "closing the fracture" at the mesh interior, we can have interface elements 
  // sharing nodes at both its sides. On that scenario, using only the 3 smallest vertices 
  // might not really describe a face. Look at the example below for two Hex elements 
  //        H__________G......M__________L
  //        /|        /|      /|        /|
  //       / |       / |     / |       / |
  //     E/__|_____F/......F/__|______/K |
  //      |  |______|__|....|..|______|__|
  //      |  /D     |  /C   |  /C     |  /J
  //      | /       | /     | /       | /
  //      |/________|/......|/________|/ 
  //      A         B       B         I
  // separated by an interface element (the dotted box). Notice that nodes B, C and F are common 
  // to the two hexs, but M and G are distinct. The interface element between them is given by 
  // nodes B,B,C,C,F,F,M,G. Since the only difference between the right face from the first 
  // hex and the left face from the second hex is the G x M node, they can not be distinguished 
  // by 3 nodes only. At present, we assume that no face has more than 4 vertices and handle that 
   // case by the V4 class. It also stores an ordered set of nodes. If the face has only 3 nodes, the 
  // extra node is completed with -1.
  struct V3
  {
    int v[3];
    bool operator<(const struct V3& other) const
    {
      return (v[0] < other.v[0]) ||
             (v[0] == other.v[0] && v[1] < other.v[1]) ||
             (v[0] == other.v[0] && v[1] == other.v[1] && v[2] < other.v[2]);
    }
    bool operator==(const struct V3& other) const
    {
      return v[0] == other.v[0] && v[1] == other.v[1] && v[2] == other.v[2];
    }
  };
 
  struct V4
  {
    int v[4];
    bool operator<(const struct V4& other) const
    {
      return (v[0] < other.v[0]) ||
             (v[0] == other.v[0] && v[1] < other.v[1]) ||
             (v[0] == other.v[0] && v[1] == other.v[1] && v[2] < other.v[2]) ||
             (v[0] == other.v[0] && v[1] == other.v[1] && v[2] == other.v[2] && v[3] < other.v[3]);
    }
    bool operator==(const struct V4& other) const
    {
      return v[0] == other.v[0] && v[1] == other.v[1] && v[2] == other.v[2] && v[3] == other.v[3];
    }
  };
 
  static size_t makeSurfaceKey(QVarLengthArray<int, 10>& vertexList)
  {
    assert(vertexList.size() == 2);
    return UibhmTopology<Geometry>::encodeEdge(vertexList[0], vertexList[1]);
  };
 
  static V3 makeSolidKey3(QVarLengthArray<int, 10>& vertexList)
  {
    assert(vertexList.size() > 2);
    qSort(vertexList);
    return V3{ vertexList[0], vertexList[1], vertexList[2] };
  };
 
  static V4 makeSolidKey4(QVarLengthArray<int, 10>& vertexList)
  {
    assert(vertexList.size() > 2 && vertexList.size() <= 4);
    qSort(vertexList);
    return V4{ vertexList[0], vertexList[1], vertexList[2], vertexList.size() == 4 ? vertexList[3] : -1};
  };
 
  // Let the query classes access our internal structure
  friend class UibhmQuery<Geometry>;
  friend class UibhmSurfaceQuery<Geometry>;
  friend class UibhmSolidQuery<Geometry>;
 
  int  _nhalf;          
  int  _lastAddedCells; 
  bool _hasOffset;      
  bool _needsEH;        
 
  HSVT _O;
 
  IVT _Of; 
 
  HSVT _VH;
 
  QVector<HalfSide> _B;
 
  QMap<size_t, HalfSide> _EH;
 
  CheckOptions         _extraChecks;   
  int                  _checking;      
  const GmLogCategory& _checkLogger;   
};
 
 
#endif