Merge branch 'feature/CFI_Mesh_IO_fixed' into 'master'

feature/cfi-mesh-io

See merge request !864

CHANGELOG.md

@@ -27,6 +27,9 @@ v5.0.0
   function definitions for the Absorption Forcing (!769)
 - Improve performance of DisContField2D::v_ExtractTracePhys (!824)
 - Fix small bug in Jacobian Energy (!857)
+- Adds CFI CAD engine back-end (!864)
+- Adds CFI Mesh IO support (!864)
+- Cleanup of CAD system data structures (!864)
 
 **NekMesh**:
 - Add feature to read basic 2D geo files as CAD (!731)

CMakeLists.txt

@@ -247,6 +247,7 @@ INCLUDE (ThirdPartyVTK)
 INCLUDE (ThirdPartyOCE)
 INCLUDE (ThirdPartyTetGen)
 INCLUDE (ThirdPartyCCM)
+INCLUDE (FindCFI)
 
 INCLUDE (Doxygen)
 

cmake/FindCFI.cmake

+########################################################################
+#
+# ThirdParty configuration for Nektar++
+#
+# OpenCascade
+#
+########################################################################
+
+IF(NEKTAR_USE_MESHGEN)
+
+    OPTION(NEKTAR_USE_CFI "Use CFI as cad engine." OFF)
+
+    IF (NEKTAR_USE_CFI)
+
+        SET(TEST_ENV $ENV{CFI_DIR})
+        IF(NOT DEFINED TEST_ENV)
+            MESSAGE(FATAL_ERROR "Cannot build with CFI without environment variable CFI_DIR set which points to cadfix1100fcs folder in the CFI installation")
+        ENDIF()
+
+        FIND_LIBRARY(CFI_LIBRARY_API NAMES cadfixapi PATHS $ENV{CFI_DIR}/lib64)
+
+        IF(CFI_LIBRARY_API)
+            FIND_PATH (CFI_INCLUDE_DIR_HXX cadfixapi.hxx PATHS $ENV{CFI_DIR}/oocfi/cxx/cadfixapi)
+            FIND_PATH (CFI_INCLUDE_DIR cfiStandardFun.h PATHS $ENV{CFI_DIR}/include)
+
+            IF(CFI_INCLUDE_DIR)
+                SET(CFI_LIBRARIES_TMP cadfixapi extra)
+                FOREACH(CFI_LIBRARIES_TMP ${CFI_LIBRARIES_TMP})
+                    LIST(APPEND CFI_LIBRARIES $ENV{CFI_DIR}/lib64/lib${CFI_LIBRARIES_TMP}.so)
+                ENDFOREACH()
+
+                MESSAGE(STATUS "cfi libraries: ${CFI_LIBRARIES}")
+
+                INCLUDE_DIRECTORIES(NekMeshUtils ${CFI_INCLUDE_DIR_HXX})
+                INCLUDE_DIRECTORIES(NekMeshUtils ${CFI_INCLUDE_DIR})
+
+            ELSE()
+                MESSAGE(FATAL_ERROR "Cannot find cadfixapi headers")
+            ENDIF()
+        ELSE()
+            MESSAGE(FATAL_ERROR "Cannot find cadfixapi libraries")
+        ENDIF()
+    ENDIF()
+ENDIF()

library/NekMeshUtils/2DGenerator/2DGenerator.cpp

@@ -75,7 +75,7 @@ void Generator2D::Process()
     // Check that cad is 2D
     Array<OneD, NekDouble> bndBox = m_mesh->m_cad->GetBoundingBox();
     ASSERTL0(fabs(bndBox[5] - bndBox[4]) < 1.0e-7, "CAD isn't 2D");
-    
+
     if (m_mesh->m_verbose)
     {
         cout << endl << "2D meshing" << endl;
@@ -385,15 +385,14 @@ void Generator2D::MakeBL(int faceid)
         // constructed by computing a normal but found on the adjacent curve
         if (it->second.size() == 1)
         {
-            vector<pair<int, CADCurveSharedPtr> > curves =
-                it->first->GetCADCurves();
+            vector<CADCurveSharedPtr> curves = it->first->GetCADCurves();
 
             vector<EdgeSharedPtr> edges =
-                m_curvemeshes[curves[0].first]->GetMeshEdges();
+                m_curvemeshes[curves[0]->GetId()]->GetMeshEdges();
             vector<EdgeSharedPtr>::iterator ie =
                 find(edges.begin(), edges.end(), it->second[0]);
             int rightCurve =
-                (ie == edges.end()) ? curves[0].first : curves[1].first;
+                (ie == edges.end()) ? curves[0]->GetId() : curves[1]->GetId();
 
             vector<NodeSharedPtr> nodes =
                 m_curvemeshes[rightCurve]->GetMeshPoints();
@@ -430,7 +429,7 @@ void Generator2D::MakeBL(int faceid)
             new Node(m_mesh->m_numNodes++, n[0], n[1], 0.0));
         CADSurfSharedPtr s = m_mesh->m_cad->GetSurf(faceid);
         Array<OneD, NekDouble> uv = s->locuv(n);
-        nn->SetCADSurf(faceid, s, uv);
+        nn->SetCADSurf(s, uv);
         nodeNormals[it->first] = nn;
     }
 

library/NekMeshUtils/CADSystem/CADCurve.cpp

+////////////////////////////////////////////////////////////////////////////////
+//
+//  File: CADSystem.cpp
+//
+//  For more information, please see: http://www.nektar.info/
+//
+//  The MIT License
+//
+//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
+//  Department of Aeronautics, Imperial College London (UK), and Scientific
+//  Computing and Imaging Institute, University of Utah (USA).
+//
+//  License for the specific language governing rights and limitations under
+//  Permission is hereby granted, free of charge, to any person obtaining a
+//  copy of this software and associated documentation files (the "Software"),
+//  to deal in the Software without restriction, including without limitation
+//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
+//  and/or sell copies of the Software, and to permit persons to whom the
+//  Software is furnished to do so, subject to the following conditions:
+//
+//  The above copyright notice and this permission notice shall be included
+//  in all copies or substantial portions of the Software.
+//
+//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+//  DEALINGS IN THE SOFTWARE.
+//
+//  Description: cad object methods.
+//
+////////////////////////////////////////////////////////////////////////////////
+#include "CADCurve.h"
+#include "CADSurf.h"
+
+using namespace std;
+
+namespace Nektar
+{
+namespace NekMeshUtils
+{
+
+Array<OneD, NekDouble> CADCurve::NormalWRT(NekDouble t, int surf)
+{
+    Array<OneD, NekDouble> p = P(t);
+    pair<CADSurfSharedPtr, CADOrientation::Orientation> surface;
+    ASSERTL0(m_adjSurfs.size() == 1,
+             "This will only work in 2D for one surface at the moment");
+    surface = m_adjSurfs[0];
+
+    Array<OneD, NekDouble> uv = surface.first->locuv(p);
+    Array<OneD, NekDouble> d1 = surface.first->D1(uv);
+
+    NekDouble t1 = t - 1e-8;
+    NekDouble t2 = t + 1e-8;
+
+    if (surface.second == CADOrientation::eBackwards)
+    {
+        swap(t1, t2);
+    }
+
+    Array<OneD, NekDouble> uv1 = surface.first->locuv(P(t1));
+    Array<OneD, NekDouble> uv2 = surface.first->locuv(P(t2));
+
+    NekDouble du = uv2[1] - uv1[1];
+    NekDouble dv = -1.0 * (uv2[0] - uv1[0]);
+
+    Array<OneD, NekDouble> N(3, 0.0);
+    N[0] = (d1[3] * du + d1[6] * dv) / 2.0;
+    N[1] = (d1[4] * du + d1[7] * dv) / 2.0;
+    N[2] = (d1[5] * du + d1[8] * dv) / 2.0;
+
+    NekDouble mag = sqrt(N[0] * N[0] + N[1] * N[1] + N[2] * N[2]);
+    N[0] /= mag;
+    N[1] /= mag;
+    N[2] /= mag;
+
+    return N;
+}
+
+CADOrientation::Orientation CADCurve::GetOrienationWRT(int surf)
+{
+    for (int i = 0; i < m_adjSurfs.size(); i++)
+    {
+        if (m_adjSurfs[i].first->GetId() == surf)
+        {
+            return m_adjSurfs[i].second;
+        }
+    }
+
+    ASSERTL0(false, "surf not in adjecency list");
+    return CADOrientation::eUnknown;
+}
+}
+}

library/NekMeshUtils/CADSystem/CADCurve.h

@@ -36,15 +36,22 @@
 #ifndef NEKMESHUTILS_CADSYSTEM_CADCURVE
 #define NEKMESHUTILS_CADSYSTEM_CADCURVE
 
+#include <LibUtilities/BasicUtils/SharedArray.hpp>
+#include <LibUtilities/BasicUtils/NekFactory.hpp>
+
 #include <NekMeshUtils/CADSystem/CADObject.h>
-#include <NekMeshUtils/CADSystem/CADVert.h>
-#include <NekMeshUtils/CADSystem/CADSurf.h>
 
 namespace Nektar
 {
 namespace NekMeshUtils
 {
 
+// forward declarations
+class CADVert;
+typedef std::shared_ptr<CADVert> CADVertSharedPtr;
+class CADSurf;
+typedef std::shared_ptr<CADSurf> CADSurfSharedPtr;
+
 /**
  * @brief base class for CAD curves.
  *
@@ -97,6 +104,9 @@ public:
      */
     virtual Array<OneD, NekDouble> D2(NekDouble t) = 0;
 
+    /**
+     * @brief Calculates the radius of curvature of the curve at point t
+     */
     virtual NekDouble Curvature(NekDouble t) = 0;
 
     /**
@@ -126,9 +136,11 @@ public:
     }
 
     /*
-     * @brief returns the ids of neigbouring surfaces
+     * @brief returns the ids of surfaces bound by this curve as well as their
+     *        Orientation with respect to the loop of curves
      */
-    std::vector<std::pair<CADSurfSharedPtr, CADOrientation::Orientation> > GetAdjSurf()
+    std::vector<std::pair<CADSurfSharedPtr, CADOrientation::Orientation> >
+        GetAdjSurf()
     {
         return m_adjSurfs;
     }
@@ -159,36 +171,36 @@ public:
     }
 
     /*
-     * @brief locates a point in the parametric space
+     * @brief locates a point in the parametric space. returns the
+     * distance to the point and passes t by reference and updates it
      */
-    virtual NekDouble loct(Array<OneD, NekDouble> xyz) = 0;
-    
-    virtual NekDouble DistanceTo(Array<OneD, NekDouble> xyz) = 0;
+    virtual NekDouble loct(Array<OneD, NekDouble> xyz,
+                           NekDouble &t) = 0;
 
-    CADOrientation::Orientation GetOrienationWRT(int surf)
-    {
-        for(int i = 0; i < m_adjSurfs.size(); i++)
-        {
-            if(m_adjSurfs[i].first->GetId() == surf)
-            {
-                return m_adjSurfs[i].second;
-            }
-        }
-
-        ASSERTL0(false,"surf not in adjecency list");
-        return CADOrientation::eUnknown;
-    }
+    /**
+     * @brief Returns the orientation of the curve with respect to a given
+     * surface by id surf
+     */
+    CADOrientation::Orientation GetOrienationWRT(int surf);
 
-    virtual Array<OneD, NekDouble> NormalWRT(NekDouble t, int surf)=0;
-    virtual Array<OneD, NekDouble> N(NekDouble t)=0;
+    /**
+     * @brief Returns the normal to the curve which is orientate with respect
+     * to the surface surf
+     */
+    Array<OneD, NekDouble> NormalWRT(NekDouble t, int surf);
 
+    /**
+     * @brief Returns the normal to a curve, it will always point in the concave
+     * direction
+     */
+    virtual Array<OneD, NekDouble> N(NekDouble t) = 0;
 
 protected:
-
     /// Length of edge
     NekDouble m_length;
     /// List of surfaces which this curve belongs to.
-    std::vector<std::pair<CADSurfSharedPtr, CADOrientation::Orientation> > m_adjSurfs;
+    std::vector<std::pair<CADSurfSharedPtr, CADOrientation::Orientation> >
+        m_adjSurfs;
     /// list of end vertices
     std::vector<CADVertSharedPtr> m_mainVerts;
 };
@@ -197,8 +209,7 @@ typedef std::shared_ptr<CADCurve> CADCurveSharedPtr;
 
 typedef LibUtilities::NekFactory<std::string, CADCurve> CADCurveFactory;
 
-CADCurveFactory& GetCADCurveFactory();
-
+CADCurveFactory &GetCADCurveFactory();
 }
 }
 

library/NekMeshUtils/CADSystem/CADObject.h

@@ -80,24 +80,46 @@ public:
     }
 
     /**
-     * @brief Return ID of the vertex
+     * @brief Return ID of the CAD object
      */
     int GetId()
     {
         return m_id;
     }
 
+    /**
+     * @brief Get the type of the CAD object
+     */
     CADType::cadType GetType()
     {
         return m_type;
     }
 
+    /**
+     * @brief Get the Orientation of the CAD object
+     */
     virtual CADOrientation::Orientation Orientation()
     {
         ASSERTL0(false,"must be implemented at the cad object level");
         return CADOrientation::eUnknown;
     }
 
+    /**
+     * @brief Give the CAD object a string name
+     */
+    void SetName(std::string i)
+    {
+        m_name = i;
+    }
+
+    /**
+     * @brief Get the name of a CAD object
+     */
+    std::string GetName()
+    {
+        return m_name;
+    }
+
 protected:
     /// ID of the vert.
     int m_id;
@@ -105,6 +127,8 @@ protected:
     CADType::cadType m_type;
     /// orientation of the CADObject
     CADOrientation::Orientation m_orientation;
+    /// string name of the cad
+    std::string m_name;
 };
 
 typedef std::shared_ptr<CADObject> CADObjectSharedPtr;

library/NekMeshUtils/CADSystem/CADSurf.cpp

@@ -50,16 +50,26 @@ namespace Nektar
 namespace NekMeshUtils
 {
 
+Array<OneD, NekDouble> CADSurf::locuv(Array<OneD, NekDouble> p)
+{
+    NekDouble dist;
+    Array<OneD, NekDouble> uv = locuv(p,dist);
+
+    WARNINGL1(dist < 1e-3, "large locuv distance");
+
+    return uv;
+}
+
 void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
-                             vector<CADSystem::EdgeLoopSharedPtr> &ein)
+                             vector<EdgeLoopSharedPtr> &ein)
 {
     // this piece of code orientates the surface,
     // it used to be face mesh but its easier to have it here
     int np = 20;
-    vector<vector<Array<OneD, NekDouble> > > loopt;
+    vector<vector<Array<OneD, NekDouble>>> loopt;
     for (int i = 0; i < ein.size(); i++)
     {
-        vector<Array<OneD, NekDouble> > loop;
+        vector<Array<OneD, NekDouble>> loop;
         for (int j = 0; j < ein[i]->edges.size(); j++)
         {
             Array<OneD, NekDouble> bnds = ein[i]->edges[j]->GetBounds();
@@ -69,7 +79,7 @@ void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
                 for (int k = 0; k < np - 1; k++)
                 {
                     NekDouble t = bnds[0] + dt * k;
-                    Array<OneD, NekDouble> l  = ein[i]->edges[j]->P(t);
+                    Array<OneD, NekDouble> l = ein[i]->edges[j]->P(t);
                     Array<OneD, NekDouble> uv = surf->locuv(l);
                     loop.push_back(uv);
                 }
@@ -79,7 +89,7 @@ void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
                 for (int k = np - 1; k > 0; k--)
                 {
                     NekDouble t = bnds[0] + dt * k;
-                    Array<OneD, NekDouble> l  = ein[i]->edges[j]->P(t);
+                    Array<OneD, NekDouble> l = ein[i]->edges[j]->P(t);
                     Array<OneD, NekDouble> uv = surf->locuv(l);
                     loop.push_back(uv);
                 }
@@ -88,7 +98,7 @@ void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
         loopt.push_back(loop);
     }
 
-    vector<bg::model::polygon<point_xy, false, true> > polygons;
+    vector<bg::model::polygon<point_xy, false, true>> polygons;
 
     for (int i = 0; i < loopt.size(); i++)
     {
@@ -107,7 +117,7 @@ void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
 
         ein[i]->area = area;
 
-        point_xy cen;
+        point_xy cen(0.0, 0.0);
         bg::centroid(polygon, cen);
 
         ein[i]->center    = Array<OneD, NekDouble>(2);
@@ -161,7 +171,8 @@ void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
 
             p = point_xy(P[0], P[1]);
 
-            ASSERTL0(boost::geometry::within(p, polygons[i]), "point is not side loop");
+            ASSERTL0(boost::geometry::within(p, polygons[i]),
+                     "point is not side loop");
         }
     }
 }

library/NekMeshUtils/CADSystem/CADSurf.h

@@ -36,9 +36,10 @@
 #ifndef NekMeshUtils_CADSYSTEM_CADSURF
 #define NekMeshUtils_CADSYSTEM_CADSURF
 
+#include <LibUtilities/BasicUtils/NekFactory.hpp>
+#include <LibUtilities/BasicUtils/SharedArray.hpp>
+
 #include <NekMeshUtils/CADSystem/CADObject.h>
-#include <NekMeshUtils/CADSystem/CADSystem.h>
-#include <NekMeshUtils/CADSystem/CADVert.h>
 
 namespace Nektar
 {
@@ -47,6 +48,22 @@ namespace NekMeshUtils
 
 class CADCurve;
 typedef std::shared_ptr<CADCurve> CADCurveSharedPtr;
+class CADSurf;
+typedef std::shared_ptr<CADSurf> CADSurfSharedPtr;
+
+/**
+ * @brief struct which descibes a collection of cad edges which are a
+ *        loop on the cad surface
+ */
+struct EdgeLoop
+{
+    std::vector<CADCurveSharedPtr> edges;
+    std::vector<CADOrientation::Orientation> edgeo;
+    Array<OneD, NekDouble> center;
+    NekDouble area;
+};
+
+typedef std::shared_ptr<EdgeLoop> EdgeLoopSharedPtr;
 
 /**
  * @brief base class for a cad surface
@@ -70,18 +87,24 @@ public:
     {
     }
 
-    static void OrientateEdges(
-        CADSurfSharedPtr surf, std::vector<CADSystem::EdgeLoopSharedPtr> &ein);
+    /**
+     * @brief Static function which orientates the edge loop on a surface
+     */
+    static void OrientateEdges(CADSurfSharedPtr surf,
+                               std::vector<EdgeLoopSharedPtr> &ein);
 
     /**
      * @brief Get the loop structures which bound the cad surface
      */
-    std::vector<CADSystem::EdgeLoopSharedPtr> GetEdges()
+    std::vector<EdgeLoopSharedPtr> GetEdges()
     {
         return m_edges;
     }
 
-    void SetEdges(std::vector<CADSystem::EdgeLoopSharedPtr> ein)
+    /**
+     * @brief Set the edge loop
+     */
+    void SetEdges(std::vector<EdgeLoopSharedPtr> ein)
     {
         m_edges = ein;
     }
@@ -121,40 +144,42 @@ public:
      * @brief Get the x,y,z at parametric point u,v.
      *
      * @param uv Array of u and v parametric coords.
-     * @return Array of xyz location.
+     * @return Array of x,y,z location.
      */
     virtual Array<OneD, NekDouble> P(Array<OneD, NekDouble> uv) = 0;
 
     /**
      * @brief Performs a reverse look up to find u,v and x,y,z.
+     * if xyz is off the surface it will return the closest point
      *
      * @param p Array of xyz location
      * @return The parametric location of xyz on this surface
      */
-    virtual Array<OneD, NekDouble> locuv(Array<OneD, NekDouble> p) = 0;
+    virtual Array<OneD, NekDouble> locuv(Array<OneD, NekDouble> p,
+                                         NekDouble &dist) = 0;
 
     /**
-     * @brief does unconstrained locuv to project point from anywhere
-     * and calculate the distance between the orthonormal projection to the
-     * surface
-     * and the point
+     * @brief overload function of locuv ommiting the dist parameter
+     * in this case if the projection is greater than a certain distance
+     * it will produce a warning. To do large distance projection use the other
+     * locuv method
      */
-    virtual NekDouble DistanceTo(Array<OneD, NekDouble> p) = 0;
+    Array<OneD, NekDouble> locuv(Array<OneD, NekDouble> p);
 
     /**
-     * @brief takes a point from anywhere find the nearest surface point and its
-     * uv
+     * @brief Returns the bounding box of the surface
      */
-    virtual NekDouble ProjectTo(Array<OneD, NekDouble> &tp,
-                           Array<OneD, NekDouble> &uv) = 0;
-                           
     virtual Array<OneD, NekDouble> BoundingBox() = 0;
 
     /**
      * @brief returns curvature at point uv
      */
     virtual NekDouble Curvature(Array<OneD, NekDouble> uv) = 0;
-    
+
+    /**
+     * @brief Is the surface defined by a planar surface (i.e not nurbs and is
+     * flat)
+     */
     virtual bool IsPlanar() = 0;
 
     /**
@@ -164,25 +189,13 @@ public:
     {
         return m_orientation;
     }
-    
-    void SetName(std::string i)
-    {
-        m_name = i;
-    }
-    
-    std::string GetName()
-    {
-        return m_name;
-    }
 
 protected:
     /// List of bounding edges in loops with orientation.
-    std::vector<CADSystem::EdgeLoopSharedPtr> m_edges;
+    std::vector<EdgeLoopSharedPtr> m_edges;
 
     /// Function which tests the the value of uv used is within the surface
     virtual void Test(Array<OneD, NekDouble> uv) = 0;
-    
-    std::string m_name;
 };
 
 typedef std::shared_ptr<CADSurf> CADSurfSharedPtr;
@@ -190,7 +203,7 @@ typedef std::shared_ptr<CADSurf> CADSurfSharedPtr;
 typedef LibUtilities::NekFactory<std::string, CADSurf> CADSurfFactory;
 
 CADSurfFactory &GetCADSurfFactory();
-}
-}
+} // namespace NekMeshUtils
+} // namespace Nektar
 
 #endif

library/NekMeshUtils/CADSystem/CADSystem.cpp

@@ -102,11 +102,5 @@ Array<OneD, NekDouble> CADSystem::GetPeriodicTranslationVector(int first,
     return ret;
 }