gmCellGeometry.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 _GEMA_CELL_GEOMETRY_H_
#define _GEMA_CELL_GEOMETRY_H_
 
#include "gmCoreConfig.h"
#include "gmMatrix.h"
#include "gmVector.h"
 
#include <QAtomicPointer>
#include <QMutex>
 
#include "gmCellType.h"
#include "gmCellGeometryInfo.h" 
 
class GmIntegrationRule;
 
class GMC_API_EXPORT GmCellGeometry
{
public:
  GmCellGeometry(GmCellType type, int P = 0, int Q = 0) 
    : _info(_infoRegistry[type]), _P(P), _Q(Q)
  { 
    assert(type < GM_NUM_CELL_TYPES);  
    assert(_info); 
  }
 
  GmCellType type() const { return _info->_type; }
 
  const char* typeName() const { return _info->_name; }
 
  GmCellFamilyType family() const { return _info->_family; }
 
  bool isInterface() const { return _info->_interface; }
 
  bool isHierarchical() const { return _info->_hierarchical; }
 
  int order() const { return _info->_order; }
 
  GmCellType linearElement() const { return _info->_eqLinearType; }
 
  int numNodes() const { return _info->_nnodes; }
 
  int numVertices() const { return _info->_nvertices; }
 
  int numExtraDofNodes() const { return _info->_nextraDofNodes; }
 
  int numCoord() const { return _info->_ncoord; }
 
  int numEdges() const { return _info->_nedges; }
 
  int numFaces() const { return _info->_nfaces; }
 
  bool isLine() const { return numFaces() == 0; }
 
  bool isSurface() const { return numFaces() == 1; }
 
  bool isSolid() const { return numFaces() > 1; }
 
  int numEdgeNodes(int edgeIndex) const
  {
    assert(edgeIndex >= 0 && edgeIndex < numEdges());
    return _info->_edgeNodeInfo[edgeIndex]._edgeNodes.size();
  }
 
  int edgeNode(int edgeIndex, int i) const
  {
    assert(i >= 0 && i < numEdgeNodes(edgeIndex));
    return _info->_edgeNodeInfo[edgeIndex]._edgeNodes[i];
  }
 
  int firstEdgeNode(int edgeIndex) const
  {
    assert(edgeIndex >= 0 && edgeIndex < numEdges());
    return _info->_edgeNodeInfo[edgeIndex]._edgeNodes.first();
  }
 
  int lastEdgeNode(int edgeIndex) const
  {
    assert(edgeIndex >= 0 && edgeIndex < numEdges());
    return _info->_edgeNodeInfo[edgeIndex]._edgeNodes.last();
  }
 
  bool nodeOnEdge(int nodeIndex, int edgeIndex) const
  {
    assert(nodeIndex >= 0 && nodeIndex < numNodes());
    assert(edgeIndex >= 0 && edgeIndex < numEdges());
    return _info->_edgeNodeInfo[edgeIndex]._edgeNodes.indexOf(nodeIndex) >= 0;
  }
 
  int edgeFromNodes(int n1, int n2) const
  {
    assert(n1 >= 0 && n1 < numNodes() && n2 >= 0 && n2 < numNodes());
    for(int e = 0; e < numEdges(); e++)
    {
      int en1 = _info->_edgeNodeInfo[e]._edgeNodes.first();
      int en2 = _info->_edgeNodeInfo[e]._edgeNodes.last();
      if((n1 == en1 && n2 == en2) || (n1 == en2 && n2 == en1))
        return e;
    }
    return -1;
  }
 
  int numFaceEdges(int faceIndex) const
  {
    assert(faceIndex >= 0 && faceIndex < numFaces());
    return _info->_faceEdgeInfo[faceIndex].size();
  }
 
  int faceEdge(int faceIndex, int i) const
  {
    assert(i >= 0 && i < numFaceEdges(faceIndex));
    return _info->_faceEdgeInfo[faceIndex][i];
  }
 
  int numEdgeFaces(int edgeIndex) const
  {
      assert(edgeIndex >= 0 && edgeIndex < numEdges());
      return _info->_edgeFaceInfo[edgeIndex].size();
  }
 
  int edgeFace(int edgeIndex, int i) const
  {
      assert(i >= 0 && i < numEdgeFaces(edgeIndex));
      return _info->_edgeFaceInfo[edgeIndex][i];
  }
 
  int numFaceBoundaryNodes(int faceIndex) const
  {
    assert(faceIndex >= 0 && faceIndex < numFaces());
    return _info->_faceNodeInfo[faceIndex]._nboundaryNodes;
  }
 
  int numFaceNodes(int faceIndex) const
  {
    assert(faceIndex >= 0 && faceIndex < numFaces());
    return _info->_faceNodeInfo[faceIndex]._faceNodes.size();
  }
 
  int faceNode(int faceIndex, int i) const
  {
    assert(i >= 0 && i < numFaceNodes(faceIndex));
    return _info->_faceNodeInfo[faceIndex]._faceNodes[i];
  }
 
  int numIncidenceNodes(int nodeIndex) const
  {
    assert(nodeIndex >= 0 && nodeIndex < numNodes());
    return _info->_nodeIncidenceInfo[nodeIndex].size();
  }
 
  int incidenceNode(int nodeIndex, int i) const
  {
    assert(i >= 0 && i < numIncidenceNodes(nodeIndex));
    return _info->_nodeIncidenceInfo[nodeIndex][i];
  }
 
  int numVolumeInternalNodes() const { return _info->_volumeNodeInfo.size(); }
 
  int volumeInternalNode(int i) const
  {
    assert(i >= 0 && i < numVolumeInternalNodes());
    return _info->_volumeNodeInfo[i];
  }
 
  int numFaceTypes() const { return _info->_nfaceTypes; }
 
  GmCellFaceType faceType(int faceIndex) const
  {
    int nedges = numFaceEdges(faceIndex);
    assert(nedges == 3 || nedges == 4);
    return (nedges == 3 ? GM_TRI_FACE : GM_QUAD_FACE);
  }
 
  GmCellType edgeElement(int edgeIndex, bool mesh = false) const
  {
    assert(edgeIndex >= 0 && edgeIndex < numEdges());
    GmCellType barType = mesh ? _info->_edgeNodeInfo[edgeIndex]._meshEdgeType : 
                                _info->_edgeNodeInfo[edgeIndex]._eqEdgeType;
    assert(barType == GM_INV_CELL_TYPE ||
           (GmCellGeometry(barType).isLine() &&
            GmCellGeometry(barType).numCoord() == numCoord() &&
            (mesh || GmCellGeometry(barType).numNodes() == numEdgeNodes(edgeIndex))));
    return barType;
  }
 
  GmCellType faceElement(int faceIndex, bool mesh = false) const
  {
    assert(faceIndex >= 0 && faceIndex < numFaces());
    GmCellType faceType = mesh ? _info->_faceNodeInfo[faceIndex]._meshFaceType :
                                 _info->_faceNodeInfo[faceIndex]._eqFaceType;
    assert(faceType == GM_INV_CELL_TYPE ||
           (GmCellGeometry(faceType).isSurface() &&
            GmCellGeometry(faceType).numCoord() == numCoord() &&
            (mesh || GmCellGeometry(faceType).numNodes() == numFaceNodes(faceIndex))));
      return faceType;
  }
 
  int edgeElementNode(int edgeIndex, int i, bool mesh = false) const
  {
    Q_UNUSED(mesh);
    assert(edgeElement(edgeIndex, mesh) != GM_INV_CELL_TYPE);
    assert(i >= 0 && i < GmCellGeometry(edgeElement(edgeIndex, mesh)).numNodes());
    if(i == 0)
      return firstEdgeNode(edgeIndex);
    else if(i == 1)
      return lastEdgeNode(edgeIndex);
    else
    { 
      assert(i == 2 && (mesh || numEdgeNodes(edgeIndex) == 3));
      return numEdgeNodes(edgeIndex) == 3 ? edgeNode(edgeIndex, 1) : -1;
    }
  }
 
  int faceElementNode(int faceIndex, int i, bool mesh = false) const
  {
    assert(faceElement(faceIndex, mesh) != GM_INV_CELL_TYPE);
    assert(i >= 0 && i < GmCellGeometry(faceElement(faceIndex, mesh)).numNodes());
    return mesh ? _info->_faceNodeInfo[faceIndex]._meshTypeNodes[i] :
                  _info->_faceNodeInfo[faceIndex]._faceTypeNodes[i];
  }
 
  GmCellType edgeLinearElement(int edgeIndex, bool mesh = false) const { return GmCellGeometry(linearElement()).edgeElement(edgeIndex, mesh); }
 
  GmCellType faceLinearElement(int faceIndex, bool mesh = false) const { return GmCellGeometry(linearElement()).faceElement(faceIndex, mesh); }
 
  int numInterfaceContactSides() const { return _info->_interfaceInfo._contactSides.size(); }
 
  int interfaceContactSide(int i) const
  {
    assert(i >= 0 && i < numInterfaceContactSides());
    return _info->_interfaceInfo._contactSides[i];
  }
 
  int numInterfaceApertureSides() const { return _info->_interfaceInfo._apertureSides.size(); }
 
  int interfaceApertureSide(int i) const
  {
    assert(i >= 0 && i < numInterfaceApertureSides());
    return _info->_interfaceInfo._apertureSides[i];
  }
 
  int numInterfaceContactEdges(int faceIndex) const
  { 
    assert(faceIndex >= 0 && faceIndex < numFaces());
    return _info->_interfaceInfo._contactEdges[faceIndex].size();
  }
 
  int interfaceContactEdge(int faceIndex, int i) const
  {
    assert(faceIndex >= 0 && faceIndex < numFaces());
    assert(i >= 0 && i < numInterfaceContactEdges(faceIndex));
    return _info->_interfaceInfo._contactEdges[faceIndex][i];
  }
 
  int numInterfaceApertureEdges(int faceIndex) const
  {
    assert(faceIndex >= 0 && faceIndex < numFaces());
    return _info->_interfaceInfo._apertureEdges[faceIndex].size();
  }
 
  int interfaceApertureEdge(int faceIndex, int i) const
  {
    assert(faceIndex >= 0 && faceIndex < numFaces());
    assert(i >= 0 && i < numInterfaceApertureEdges(faceIndex));
    return _info->_interfaceInfo._apertureEdges[faceIndex][i];
  }
 
  static int maxNumNodes() { return _maxNumNodes; }
 
  static const char* familyTypeToStr(GmCellFamilyType family);
  static int         strToFamilyType(QString str);
 
  static const char* faceTypeToStr(GmCellFaceType type);
  static int         strToFaceType(QString str);
 
  static int name2Type(const QString name) {
    static QHash<QString, int> nameToType;
    static QMutex mtx;
    //build entries
    if(nameToType.size() != GM_NUM_CELL_TYPES) {
      mtx.lock();
      //In case it was being built while we entered the previous if-clause
      if(nameToType.isEmpty()) {
        for(int i = 0; i < GM_NUM_CELL_TYPES; ++i) {
          nameToType[_infoRegistry[i]->_name] = i;
        }
      }
      mtx.unlock();
    }
    assert(nameToType.size() == GM_NUM_CELL_TYPES);
    auto it = nameToType.find(name);
    if(it == nameToType.end()) {
      return -1;
    }
    return it.value();
  }
 
  static const char* type2Name(GmCellType t) {
    return _infoRegistry[t]->_name;
  }
 
  double length(const GmMatrix& X) const
  {
    checkX(X);
    if(isSolid() || (isSurface() && !isInterface()))
      return 0.0;
    return _info->dimension(X);
  }
 
  double area(const GmMatrix& X) const
  {
    checkX(X);
    if(isLine() || (isSurface() && isInterface()) || (isSolid() && !isInterface()))
      return 0.0;
    return _info->dimension(X);
  }
 
  double volume(const GmMatrix& X) const
  {
    checkX(X);
    if(!isSolid() || isInterface())
      return 0.0;
    return _info->dimension(X);
  }
 
  double characteristicLength(const GmMatrix& X) const
  {
    checkX(X);
    if(isLine() || (isSurface() && isInterface()))
      return _info->dimension(X);
    else if(isSurface() || (isSolid() && isInterface()))
      return std::sqrt(_info->dimension(X));
    else
      return std::pow(_info->dimension(X), 1.0/3.0);
  }
 
  double characteristicDimension(const GmMatrix& X) const { checkX(X); return _info->dimension(X); }
 
  bool localAxis(const GmMatrix& X, GmMatrix& R) const { checkX(X); return _info->localAxis(X, R); }
 
  void centroidCartesian(const GmMatrix& X, GmVector& coord) const { checkX(X); return _info->centroidCartesian(X, coord); }
 
  // Capability related methods.
  //------------------------------------------------------------------------------------
 
  bool hasCapability(GmCellGeometryCapabilities capability) const 
  { 
    assert(capability >= 0 && capability < GM_CELL_GEOMETRY_CAPABILITY_COUNT); 
    return _info->_capabilities[capability]; 
  }
 
  int isValid(const GmMatrix& X, double tol) const { checkX(X); return _info->isValid(X, tol); }
 
  double quality(const GmMatrix& J, double tol) const { return _info->quality(J, tol); }
 
  bool contains(const GmMatrix& X, const GmVector& coord) const { checkX(X); return _info->contains(X, coord); }
 
  //------------------------------------------------------------------------------------
  // CellGeometryInfo registering and query functions. Not for general use
  //------------------------------------------------------------------------------------
 
  static void registerCellType(const GmCellGeometryInfo* cellInfo);
 
  static const GmCellGeometryInfo* geometryInfo(GmCellType type) 
  {
    assert(type >= 0 && type < GM_NUM_CELL_TYPES);
    return _infoRegistry[type];
  }
 
  //------------------------------------------------------------------------------------
  // Shape & Integration rule factory functions.
  // Marked as private and accessed by selected friend classes since they should NOT be 
  // used by regular code. Shape functions and integration rules should be obtained through
  // the Shape and IntegrationRule factories.
  //------------------------------------------------------------------------------------
 
private:
  friend class GmShape;
  friend class GmIntegrationRule;
  friend class GmBorderIntegrationRule;
 
  GmShape* shapeInstance() const 
  { 
    assert((!_info->_hierarchical && _P == 0 && _Q == 0) || (_info->_hierarchical && _P > 0 && _Q > 0));
    return _info->shapeInstance(_P, _Q); 
  }
 
  GmIntegrationRule* integrationRule(GmIntegrationRuleType irType, int rule1 = -1, int rule2 = -1, int rule3 = -1) const
  {
    assert((!_info->_hierarchical && _P == 0 && _Q == 0) || (_info->_hierarchical && _P > 0 && _Q > 0));
    return _info->integrationRule(irType, rule1, rule2, rule3, _P, _Q);
  }
 
  GmBorderIntegrationRule* edgeIntegrationRule(GmIntegrationRuleType irType, int rule1 = -1) const
  {
    assert((!_info->_hierarchical && _P == 0 && _Q == 0) || (_info->_hierarchical && _P > 0 && _Q > 0));
    return _info->edgeIntegrationRule(irType, rule1, _P, _Q);
  }
 
  GmBorderIntegrationRule* faceIntegrationRule(int faceType, GmIntegrationRuleType irType, int rule1 = -1, int rule2 = -1) const
  {
    assert((!_info->_hierarchical && _P == 0 && _Q == 0) || (_info->_hierarchical && _P > 0 && _Q > 0));
    return _info->faceIntegrationRule(faceType, irType, rule1, rule2, _P, _Q);
  }
 
  GmBorderIntegrationRule* borderIntegrationRule(int faceType, GmIntegrationRuleType irType, int rule1 = -1, int rule2 = -1) const
  {
    if(isSurface())
      return edgeIntegrationRule(irType, rule1);
    else if(isSolid())
      return faceIntegrationRule(faceType, irType, rule1, rule2);
    else  // line
      return NULL;
  }
 
private:
  inline void checkX(const GmMatrix& X) const { Q_UNUSED(X); assert(X.n_cols == numNodes() && X.n_rows == numCoord()); }
 
#ifndef NDEBUG
  static void checkMetadata(const GmCellGeometryMetadata* data, bool registryFullyInitialised);
  static void checkGeometryIEFCompatibility(GmCellType type);
  static void checkElementRegistry();
#endif
  
  const GmCellGeometryInfo* _info;   
  int                       _P, _Q;  
 
  static const GmCellGeometryInfo* _infoRegistry[GM_NUM_CELL_TYPES]; 
  static int                       _maxNumNodes;                     
};
 
 
 
#endif