Add Myles' fix to WA interface

Author: EduardoMolina

Common/include/interface_interpolation/CSlidingMesh.hpp

@@ -1,4 +1,4 @@
-/*!
+/*!
  * \file CSlidingMesh.hpp
  * \brief Sliding mesh interpolation.
  * \author H. Kline
@@ -67,6 +67,21 @@ class CSlidingMesh final : public CInterpolator {
                                       const unsigned long* nNeighbor, const su2double *coord,
                                       unsigned long centralNode, su2double** element);
 
+  /*!
+   * \brief For 3-Dimensional grids, build the dual surface element
+   * \param[in] map         - array containing the index of the boundary points connected to the node
+   * \param[in] startIndex  - for each vertex specifies the corresponding index in the global
+   *                          array containing the indexes of all its neighbouring vertexes
+   * \param[in] nNeighbour  - for each vertex specifies the number of its neighbouring vertexes (on the boundary)
+   * \param[in] coord       - array containing the coordinates of all the boundary vertexes
+   * \param[in] centralNode - label of the vertex around which the dual surface element is built
+   * \param[out] element  - double array where element node coordinates will be stored
+   * \return Number of points included in the element.
+   */
+  static int Build_DualElement(const unsigned long *map, const unsigned long *startIndex,
+                                 const unsigned long* nNeighbor, const su2double *coord,
+                                 unsigned long centralNode, su2double** element);
+
   /*!
    * \brief For 2-Dimensional grids, compute intersection length of two segments projected along a given direction
    * \param[in] nDim - Number of dimensions

Common/src/interface_interpolation/CSlidingMesh.cpp

@@ -469,14 +469,15 @@ void CSlidingMesh::SetTransferCoeff(const CConfig* const* config) {
 
         nEdges_target = Target_nLinkedNodes[target_iPoint];
 
-        nNode_target = 2*(nEdges_target + 1);
-
-        target_element = new su2double*[nNode_target];
-        for (ii = 0; ii < nNode_target; ii++)
+        target_element = new su2double*[ 2*nEdges_target + 2 ];
+        for (ii = 0; ii < 2*nEdges_target + 2; ii++)
           target_element[ii] = new su2double[nDim];
 
-        nNode_target = Build_3D_surface_element(Target_LinkedNodes, Target_StartLinkedNodes, Target_nLinkedNodes,
-                                                TargetPoint_Coord, target_iPoint, target_element);
+//      nNode_target = Build_3D_surface_element(Target_LinkedNodes, Target_StartLinkedNodes, Target_nLinkedNodes,
+//                                                TargetPoint_Coord, target_iPoint, target_element);
+
+        nNode_target = Build_DualElement(Target_LinkedNodes, Target_StartLinkedNodes, Target_nLinkedNodes,
+                                           TargetPoint_Coord, target_iPoint, target_element);
 
         /*--- Brute force to find the closest donor_node ---*/
 
@@ -508,8 +509,11 @@ void CSlidingMesh::SetTransferCoeff(const CConfig* const* config) {
         for (ii = 0; ii < 2*nEdges_donor + 2; ii++)
           donor_element[ii] = new su2double[nDim];
 
-        nNode_donor = Build_3D_surface_element(Donor_LinkedNodes, Donor_StartLinkedNodes, Donor_nLinkedNodes,
-                                               DonorPoint_Coord, donor_iPoint, donor_element);
+//      nNode_donor = Build_3D_surface_element(Donor_LinkedNodes, Donor_StartLinkedNodes, Donor_nLinkedNodes,
+//                                               DonorPoint_Coord, donor_iPoint, donor_element);
+
+        nNode_donor = Build_DualElement(Donor_LinkedNodes, Donor_StartLinkedNodes, Donor_nLinkedNodes,
+                                          DonorPoint_Coord, donor_iPoint, donor_element);
 
         Area = 0;
         for (ii = 1; ii < nNode_target-1; ii++){
@@ -609,8 +613,11 @@ void CSlidingMesh::SetTransferCoeff(const CConfig* const* config) {
                 for (ii = 0; ii < 2*nEdges_donor + 2; ii++)
                   donor_element[ii] = new su2double[nDim];
 
-                nNode_donor = Build_3D_surface_element(Donor_LinkedNodes, Donor_StartLinkedNodes, Donor_nLinkedNodes,
-                                                       DonorPoint_Coord, donor_iPoint, donor_element);
+//              nNode_donor = Build_3D_surface_element(Donor_LinkedNodes, Donor_StartLinkedNodes, Donor_nLinkedNodes,
+//                                                       DonorPoint_Coord, donor_iPoint, donor_element);
+
+                nNode_donor = Build_DualElement(Donor_LinkedNodes, Donor_StartLinkedNodes, Donor_nLinkedNodes,
+                                                  DonorPoint_Coord, donor_iPoint, donor_element);
 
                 tmp_Area = 0;
                 for (ii = 1; ii < nNode_target-1; ii++)
@@ -839,6 +846,127 @@ int CSlidingMesh::Build_3D_surface_element(const unsigned long *map, const unsig
 
 }
 
+int CSlidingMesh::Build_DualElement(const unsigned long *map, const unsigned long *startIndex,
+                                    const unsigned long* nNeighbor, const su2double *coord,
+                                    unsigned long centralNode, su2double** element){
+
+    /*--- Given a node "centralNode", this routines reconstruct the vertex centered surface quadrilateral element around the node and stores it into "element" ---*/
+    /*--- Returns the number of points included in the element ---*/
+
+    unsigned long iNode, jNode, kNode, ElementIndex, iPoint, jPoint, kPoint, nOuterNodes, nTriNodes;
+    unsigned short nDim = 3, iDim, nTmp;
+    const unsigned long *OuterNodes, *ptr;
+    int NodeIndex;
+    unsigned long **DualElement;
+
+    nTriNodes = 3;
+
+    /* --- Store central node as element first point --- */
+    for (iDim = 0; iDim < nDim; iDim++) {
+        element[0][iDim] = coord[centralNode * nDim + iDim];
+    }
+
+    /* --- Get the number of nodes directly connected to the central node  --- */
+    nOuterNodes = nNeighbor[centralNode];
+
+    /* --- Array containing the points connected to the central node --- */
+    OuterNodes = &map[startIndex[centralNode]];
+
+    /* --- Allocate auxiliary structure, vectors are longer than needed but this avoid further re-allocations due to length variation --- */
+    DualElement = new unsigned long*[nOuterNodes];
+    for ( int i = 0; i < nOuterNodes; i++ ) {
+        DualElement[i] = new unsigned long[nTriNodes];
+    }
+
+    /*--- Initialise dual elements ---*/
+    for ( int i = 0; i < nOuterNodes; i++ ){
+        for (int j = 0; j < nTriNodes; j++) {
+            DualElement[i][j] = -1;
+        }
+    }
+
+    NodeIndex = 0;
+
+    /* --- Build and order tri segments of a hexa dual element --- */
+    while ((nOuterNodes - NodeIndex) > 0) {
+
+        for (iNode = 0; iNode < nOuterNodes && ((nOuterNodes - NodeIndex) > 0); iNode++) {
+            iPoint = OuterNodes[iNode];
+            ptr = &map[startIndex[iPoint]];
+            nTmp = nNeighbor[iPoint];
+
+            for (jNode = 0; jNode < nTmp && ((nOuterNodes - NodeIndex) > 0); jNode++) {
+                jPoint = ptr[jNode];
+
+                for (kNode = 0; kNode < nOuterNodes && ((nOuterNodes - NodeIndex) > 0); kNode++) {
+                    kPoint = OuterNodes[kNode];
+
+                    /*--- Find the shared outer nodes in order ---*/
+                    if (jPoint == kPoint) {
+
+                        if (NodeIndex == 0) {
+                            DualElement[NodeIndex][0] = centralNode;
+                            DualElement[NodeIndex][1] = iPoint;
+                            DualElement[NodeIndex][2] = jPoint;
+                            NodeIndex++;
+                            break;
+                        }
+
+                        if (iPoint == DualElement[NodeIndex - 1][2] && jPoint != DualElement[NodeIndex - 1][1]) {
+                            DualElement[NodeIndex][0] = centralNode;
+                            DualElement[NodeIndex][1] = iPoint;
+                            DualElement[NodeIndex][2] = jPoint;
+                            NodeIndex++;
+                            break;
+                        }
+
+                    }
+                }
+            }
+        }
+    }
+
+    NodeIndex = 0;
+    ElementIndex = 1;
+
+    /* --- Build array containing the quad dual element by finding the mid point of each edge and the baricenter of each face.
+     * Each quad is split through its diagonal connecting the central node, baricenter and 4th node.   --- */
+    while ((nOuterNodes - NodeIndex) > 0) {
+        for (iDim = 0; iDim < nDim; iDim++) {
+            element[ElementIndex][iDim] = (element[0][iDim] + coord[DualElement[NodeIndex][1] * nDim + iDim]) / 2;
+        }
+        ElementIndex++;
+
+        for (iDim = 0; iDim < nDim; iDim++) {
+            element[ElementIndex][iDim] = (element[0][iDim] + coord[DualElement[NodeIndex][1] * nDim + iDim] + coord[DualElement[NodeIndex][2] * nDim + iDim]) / 3;
+        }
+        ElementIndex++;
+        NodeIndex++;
+    }
+
+    // This is a closed element, so add again element 1 to the end of the structure, useful later
+    if(DualElement[NodeIndex - 1][2] == DualElement[0][1]){
+
+        for (iDim = 0; iDim < nDim; iDim++) {
+            element[ElementIndex][iDim] = element[1][iDim];
+        }
+        ElementIndex++;
+    }
+    else {
+        for (iDim = 0; iDim < nDim; iDim++) {
+            element[ElementIndex][iDim] = (element[0][iDim] + coord[DualElement[NodeIndex - 1][2] * nDim + iDim]) / 2;
+        }
+        ElementIndex++;
+    }
+
+    for ( int i = 0; i < nOuterNodes; i++ ) {
+        delete[] DualElement[i];
+    }
+    delete [] DualElement;
+
+    return (int)ElementIndex;
+}
+
 su2double CSlidingMesh::ComputeLineIntersectionLength(unsigned short nDim, const su2double* A1, const su2double* A2,
                                                       const su2double* B1, const su2double* B2, const su2double* Direction) {