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) const
  {
    assert(edgeIndex >= 0 && edgeIndex < numEdges());
    GmCellType barType = _info->_edgeNodeInfo[edgeIndex]._eqEdgeType;
    assert(barType == GM_INV_CELL_TYPE ||
           (GmCellGeometry(barType).isLine() &&
            GmCellGeometry(barType).numNodes() == numEdgeNodes(edgeIndex) &&
            GmCellGeometry(barType).numCoord() == numCoord()));
    return barType;
  }

  GmCellType faceElement(int faceIndex) const
  {
    assert(faceIndex >= 0 && faceIndex < numFaces());
    GmCellType faceType = _info->_faceNodeInfo[faceIndex]._eqFaceType;
    assert(faceType == GM_INV_CELL_TYPE ||
           (GmCellGeometry(faceType).isSurface() &&
            GmCellGeometry(faceType).numNodes() == numFaceNodes(faceIndex) &&
            GmCellGeometry(faceType).numCoord() == numCoord()));
    return faceType;
  }

  int edgeElementNode(int edgeIndex, int i) const
  {
    assert(edgeElement(edgeIndex) != GM_INV_CELL_TYPE);
    assert(i >= 0 && i < GmCellGeometry(edgeElement(edgeIndex)).numNodes());
    static int bar3NodeTranslation[3] = { 0, 2, 1 }; // Edge ordering to Bar3 ordering
    if(numEdgeNodes(edgeIndex) == 2)
      return edgeNode(edgeIndex, i);
    else
      return edgeNode(edgeIndex, bar3NodeTranslation[i]);
  }

  int faceElementNode(int faceIndex, int i) const
  {
    assert(faceElement(faceIndex) != GM_INV_CELL_TYPE);
    assert(i >= 0 && i < GmCellGeometry(faceElement(faceIndex)).numNodes());
    return _info->_faceNodeInfo[faceIndex]._faceTypeNodes[i];
  }

  GmCellType edgeLinearElement(int edgeIndex) const { return GmCellGeometry(linearElement()).edgeElement(edgeIndex); }

  GmCellType faceLinearElement(int faceIndex) const { return GmCellGeometry(linearElement()).faceElement(faceIndex); }

  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);

  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); }

  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]; 
  }

  bool 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