Merge branch 'dev' of https://github.com/gimli-org/gimli into dev

src/integration.cpp

@@ -44,34 +44,103 @@ IntegrationRules::~IntegrationRules(){
 // std::cout << " delete IntegrationRules" << std::endl;
 }
 
+const R3Vector & IntegrationRules::gauAbscissa(Index order) const { 
+    ASSERT_SIZE(gauAbscissa_, order)
+    return gauAbscissa_[order]; 
+}
+
+const RVector & IntegrationRules::gauWeights(Index order) const { 
+    ASSERT_SIZE(gauWeights_, order)
+    return gauWeights_[order]; 
+}
+
+const R3Vector & IntegrationRules::triGLAbscissa(Index order) const { 
+    ASSERT_SIZE(triGLAbscissa_, order)
+    return triGLAbscissa_[order]; 
+}
+
+const RVector & IntegrationRules::triGLWeights(Index order) const { 
+    ASSERT_SIZE(triGLWeights_, order)
+    return triGLWeights_[order]; 
+}
+
+const R3Vector & IntegrationRules::edgAbscissa(Index order) const { 
+    ASSERT_SIZE(edgAbscissa_, order)
+    return edgAbscissa_[order]; 
+}
+const RVector & IntegrationRules::edgWeights(Index order) const { 
+    ASSERT_SIZE(edgWeights_, order)
+    return edgWeights_[order]; 
+}
+const R3Vector & IntegrationRules::triAbscissa(Index order) const { 
+    ASSERT_SIZE(triAbscissa_, order)
+    return triAbscissa_[order]; 
+}
+const RVector & IntegrationRules::triWeights(Index order) const { 
+    ASSERT_SIZE(triWeights_, order)
+    return triWeights_[order]; 
+}
+const R3Vector & IntegrationRules::tetAbscissa(Index order) const { 
+    ASSERT_SIZE(tetAbscissa_, order)
+    return tetAbscissa_[order]; 
+}
+const RVector & IntegrationRules::tetWeights(Index order) const { 
+    ASSERT_SIZE(tetWeights_, order)
+    return tetWeights_[order]; 
+}
+const R3Vector & IntegrationRules::quaAbscissa(Index order) const { 
+    ASSERT_SIZE(quaAbscissa_, order)
+    return quaAbscissa_[order]; 
+}
+const RVector & IntegrationRules::quaWeights(Index order) const { 
+    ASSERT_SIZE(quaWeights_, order)
+    return quaWeights_[order]; 
+}
+const R3Vector & IntegrationRules::hexAbscissa(Index order) const { 
+    ASSERT_SIZE(hexAbscissa_, order)
+    return hexAbscissa_[order]; 
+}
+const RVector & IntegrationRules::hexWeights(Index order) const { 
+    ASSERT_SIZE(hexWeights_, order)
+    return hexWeights_[order]; 
+}
+const R3Vector & IntegrationRules::priAbscissa(Index order) const { 
+    ASSERT_SIZE(priAbscissa_, order)
+    return priAbscissa_[order]; 
+}
+const RVector & IntegrationRules::priWeights(Index order) const { 
+    ASSERT_SIZE(priWeights_, order)
+    return priWeights_[order]; 
+}
+
 const R3Vector & IntegrationRules::abscissa(const Shape & shape, uint order) const {
     switch(shape.rtti()){
-        case MESH_SHAPE_EDGE_RTTI: return edgAbscissa_[order];
+        case MESH_SHAPE_EDGE_RTTI: return edgAbscissa(order);
         case MESH_SHAPE_TRIANGLE_RTTI:
-             if (triUseGaussLegendre_) return triGLAbscissa_[order];
+             if (triUseGaussLegendre_) return triGLAbscissa(order);
              else return triAbscissa_[order];
-        case MESH_SHAPE_QUADRANGLE_RTTI: return quaAbscissa_[order];
-        case MESH_SHAPE_TETRAHEDRON_RTTI: return tetAbscissa_[order];
-        case MESH_SHAPE_HEXAHEDRON_RTTI: return hexAbscissa_[order];
+        case MESH_SHAPE_QUADRANGLE_RTTI: return quaAbscissa(order);
+        case MESH_SHAPE_TETRAHEDRON_RTTI: return tetAbscissa(order);
+        case MESH_SHAPE_HEXAHEDRON_RTTI: return hexAbscissa(order);
         default:
             break;
     }
-    return gauAbscissa_[order];
+    return gauAbscissa(order);
 }
 
 const RVector & IntegrationRules::weights(const Shape & shape, uint order) const {
     switch(shape.rtti()){
-        case MESH_SHAPE_EDGE_RTTI: return edgWeights_[order];
+        case MESH_SHAPE_EDGE_RTTI: return edgWeights(order);
         case MESH_SHAPE_TRIANGLE_RTTI:
-            if (triUseGaussLegendre_) return triGLWeights_[order];
+            if (triUseGaussLegendre_) return triGLWeights(order);
             else return triWeights_[order];
-        case MESH_SHAPE_QUADRANGLE_RTTI: return quaWeights_[order];
-        case MESH_SHAPE_TETRAHEDRON_RTTI: return tetWeights_[order];
-        case MESH_SHAPE_HEXAHEDRON_RTTI: return hexWeights_[order];
+        case MESH_SHAPE_QUADRANGLE_RTTI: return quaWeights(order);
+        case MESH_SHAPE_TETRAHEDRON_RTTI: return tetWeights(order);
+        case MESH_SHAPE_HEXAHEDRON_RTTI: return hexWeights(order);
         default:
             break;
     }
-    return gauWeights_[order];
+    return gauWeights(order);
 }
 
 void IntegrationRules::initGau_(){

src/integration.h

@@ -29,43 +29,44 @@ class DLLEXPORT IntegrationRules : public Singleton< IntegrationRules >{
     friend class Singleton< IntegrationRules >;
 
     /*! Return Gauss-Legendre quadrature point upto order <10. */
-    inline const R3Vector & gauAbscissa(uint order) const { return gauAbscissa_[order]; }
-
+    const R3Vector & gauAbscissa(Index order) const;
+    
     /*! Return Gauss-Legendre quadrature weights upto order <10. */
-    inline const RVector & gauWeights(uint order) const { return gauWeights_[order]; }
+    const RVector & gauWeights(Index order) const;
 
     /*!
      * Generic quadrature positions for a triangle based on Gauss-Legendre quadrature
      H. T. RATHOD1*, K. V. NAGARAJA2, B. VENKATESUDU3 AND N. L. RAMESH4.
      Gauss Legendre quadrature over a triangle.
      J. Indian Inst. Sci., Sept.-Oct. 2004, 84, 183-188
     */
-    inline const R3Vector & triGLAbscissa(uint order) const { return triGLAbscissa_[order]; }
+    const R3Vector & triGLAbscissa(Index order) const;
+
     /*!
      * Generic quadrature weights for a triangle based on Gauss-Legendre quadrature
      H.T. RATHOD, K. V. NAGARAJA, B. VENKATESUDU AND N. L. RAMESH.
      Gauss Legendre quadrature over a triangle.
      J. Indian Inst. Sci., Sept.-Oct. 2004, 84, 183-188
     */
-    inline const RVector & triGLWeights(uint order) const { return triGLWeights_[order]; }
+    const RVector & triGLWeights(Index order) const;
 
-    inline const R3Vector & edgAbscissa(uint order) const { return edgAbscissa_[order]; }
-    inline const RVector & edgWeights(uint order) const { return edgWeights_[order]; }
+    const R3Vector & edgAbscissa(Index order) const;
+    const RVector & edgWeights(Index order) const;
 
-    inline const R3Vector & triAbscissa(uint order) const { return triAbscissa_[order]; }
-    inline const RVector & triWeights(uint order) const { return triWeights_[order]; }
+    const R3Vector & triAbscissa(Index order) const;
+    const RVector & triWeights(Index order) const;
 
-    inline const R3Vector & tetAbscissa(uint order) const { return tetAbscissa_[order]; }
-    inline const RVector & tetWeights(uint order) const { return tetWeights_[order]; }
+    const R3Vector & tetAbscissa(Index order) const;
+    const RVector & tetWeights(Index order) const;
 
-    inline const R3Vector & quaAbscissa(uint order) const { return quaAbscissa_[order]; }
-    inline const RVector & quaWeights(uint order) const { return quaWeights_[order]; }
+    const R3Vector & quaAbscissa(Index order) const;
+    const RVector & quaWeights(Index order) const;
 
-    inline const R3Vector & hexAbscissa(uint order) const { return hexAbscissa_[order]; }
-    inline const RVector & hexWeights(uint order) const { return hexWeights_[order]; }
+    const R3Vector & hexAbscissa(Index order) const;
+    const RVector & hexWeights(Index order) const;
 
-    inline const R3Vector & priAbscissa(uint order) const { return priAbscissa_[order]; }
-    inline const RVector & priWeights(uint order) const { return priWeights_[order]; }
+    const R3Vector & priAbscissa(Index order) const;
+    const RVector & priWeights(Index order) const;
     /*!
      * Return Gauss-Legendre quadrature positions for a given shape of the \ref MeshEntity upto order 10
      */

src/mesh.cpp

@@ -177,17 +177,6 @@ Node * Mesh::createNode(const RVector3 & pos, int marker){
     return createNode_(pos, marker, -1);
 }
 
-Node * Mesh::createSecondaryNode_(const RVector3 & pos){
-    Index id = this->secondaryNodeCount();
-    secNodeVector_.push_back(new Node(pos));
-    secNodeVector_.back()->setId(this->nodeCount() + id);
-    return secNodeVector_.back();
-}
-
-Node * Mesh::createSecondaryNode(const RVector3 & pos){
-    return createSecondaryNode_(pos);
-}
-
 Node & Mesh::secondaryNode(Index i) {
     ASSERT_RANGE(i, 0, this->secondaryNodeCount())
     return *secNodeVector_[i];
@@ -198,6 +187,31 @@ Node & Mesh::secondaryNode(Index i) const {
     return *secNodeVector_[i];
 }
 
+Node * Mesh::createSecondaryNode_(const RVector3 & pos){
+    Index id = this->secondaryNodeCount();
+    secNodeVector_.push_back(new Node(pos));
+    secNodeVector_.back()->setId(this->nodeCount() + id);
+    return secNodeVector_.back();
+}
+
+Node * Mesh::createSecondaryNode(const RVector3 & pos, double tol){
+    bool useTree = false;
+    if (tol > 0.0){
+        fillKDTree_();
+        useTree = true;
+
+        Node * refNode = tree_->nearest(pos);
+        if (refNode){
+            if (pos.distance(refNode->pos()) < tol) {
+                return refNode;
+            }
+        }
+    }
+    Node *newNode = createSecondaryNode_(pos);
+    if (useTree) tree_->insert(newNode);
+    return newNode;
+}
+
 Node * Mesh::createNodeWithCheck(const RVector3 & pos, double tol, bool warn, bool edgeCheck){
     bool useTree = false;
     if (tol > -1.0){

src/mesh.h

@@ -190,7 +190,10 @@ class DLLEXPORT Mesh {
     Node * createNode(const Node & node);
     Node * createNode(const RVector3 & pos, int marker=0);
 
-    Node * createSecondaryNode(const RVector3 & pos);
+    /*! Create a secondary node, which is stored in an aditional list for additional use. 
+    If tolerance tol set to a value > 0, then it will be checked if there is already a node 
+    at this position and return a ptr to the existing node instead of creating a new. */
+    Node * createSecondaryNode(const RVector3 & pos, double tol=-1);
 
     /*! Create new Node with duplication checks. Returns the already existing node when its within a tolerance distance to pos.
     If edgeCheck is set, any 2d (p1) boundary edges will be checked for any intersection with pos and splitted if necessary.*/

src/shape.cpp

@@ -643,7 +643,7 @@ double HexahedronShape::volume() const {
     TetrahedronShape tet;
     for (uint i = 0; i < 5; i ++){
         tet.setNodes(nodeVector_[HexahedronSplit5TetID[i][0]], nodeVector_[HexahedronSplit5TetID[i][1]],
-                      nodeVector_[HexahedronSplit5TetID[i][2]], nodeVector_[HexahedronSplit5TetID[i][3]]);
+                     nodeVector_[HexahedronSplit5TetID[i][2]], nodeVector_[HexahedronSplit5TetID[i][3]]);
         sum += tet.volume();
     }