subdivision/subUtil/Hexahedron.cpp

/*
XLiFE++ is an extended library of finite elements written in C++
    Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*!
  \file Hexahedron.cpp
  \author Y. Lafranche
  \since 3 Fev 2014
  \date 19 Mar 2014

  \brief Implementation of xlifepp::subdivision::Hexahedron class members and related functions
*/

#include "Hexahedron.hpp"

#include <algorithm>
using namespace std;

namespace xlifepp {
namespace subdivision {

//-------------------------------------------------------------------------------
//  Static variables
//-------------------------------------------------------------------------------
// Index array defining the edges according to the following numering convention:
// edge #1=(1,2), edge #2 =(3,4), edge #3 =(5,6), edge #4 =(7,8),
// edge #5=(1,3), edge #6 =(2,4), edge #7 =(5,7), edge #8 =(6,8),
// edge #9=(1,5), edge #10=(2,6), edge #11=(3,7), edge #12=(4,8),
short Hexahedron::rkEdge[/* nb_edges_ */][2] = {{0,1}, {2,3}, {4,5}, {6,7},
                                                {0,2}, {1,3}, {4,6}, {5,7},
                                                {0,4}, {1,5}, {2,6}, {3,7}};

// Index array defining the faces according to the following numering convention:
// face #1=(1,3,4,2), face #2=(5,7,8,6) : orthogonal to X
// face #3=(1,2,6,5), face #4=(3,4,8,7) : orthogonal to Y
// face #5=(1,5,7,3), face #6=(2,6,8,4) : orthogonal to Z
short Hexahedron::rkFace[/* nb_faces_ */][4] = {{0,2,3,1}, {4,6,7,5}, {0,1,5,4},
                                                {2,3,7,6}, {0,4,6,2}, {1,5,7,3}};

// Array containing the edge numbers (>=1) defining a face
// face #1=(5,2,6,1),   face #2=(7,4,8,3),  face #3=(1,10,3,9),
// face #4=(2,12,4,11), face #5=(9,7,11,5), face #6=(10,8,12,6)
short Hexahedron::nuEdge[/* nb_faces_ */][4] = {{5,2,6,1},   {7,4,8,3},  {1,10,3,9},
                                                {2,12,4,11}, {9,7,11,5}, {10,8,12,6}};
//-------------------------------------------------------------------------------
//  Constructors, Destructor
//-------------------------------------------------------------------------------
/*!
 Conversion constructor. Constructor by vertices, given by their rank in the list
 */
Hexahedron::Hexahedron(const number_t num,
                       const number_t rV1, const number_t rV2, const number_t rV3, const number_t rV4,
                       const number_t rV5, const number_t rV6, const number_t rV7, const number_t rV8,
                       const number_t bdSideNum)
:CartesianFig<Hexahedron>(num,bdSideNum){
   vertices_[0] = rV1;
   vertices_[1] = rV2;
   vertices_[2] = rV3;
   vertices_[3] = rV4;
   vertices_[4] = rV5;
   vertices_[5] = rV6;
   vertices_[6] = rV7;
   vertices_[7] = rV8;
}

//-------------------------------------------------------------------------------
//  Other public member functions
//-------------------------------------------------------------------------------
/*!
 returns the numbers (>= 1) of the edges defining the face number i
 */
vector<short> Hexahedron::numEdgesOfFace(const number_t i) {
   return vector<short> (nuEdge[i-1], nuEdge[i-1]+4);
}
/*!
 returns the ranks (>= 0) of the vertices defining the edges of the hexahedron
 */
vector<pair<short,short> > Hexahedron::rkEdgeVertices() {
   vector<pair<short,short> > V(nb_edges_);
   for (number_t i=0; i<nb_edges_; i++) { V[i] = make_pair(rkEdge[i][0],rkEdge[i][1]); }
   return V;
}
/*!
 This function gives the local position of the vertices of any hexahedron of
 the mesh correspondig to the numbering convention explained on top of the
 Hexahedron class header file (Hexahedron.hpp).

 The numbering convention of high order vertices is deduced from a triple tensor
 product of the standard numbering over the segment defined by the integers [0,1,...k].
 Thus, the entry BC[i] in the returned vector BC contains a vector of 3 integers whose
 components BC[i][j] are such that 0<= BC[i][j] <= k, where k is the order of the mesh.
 Each component is the result of a permutation of a 3D tensor product of the uniform
 subdivision of the interval [0,1,...k]. This gives an intrinsic way of numbering and
 localization of the points if these components are sorted in alphanumeric order.
 Nota: this does not match with the numbering of the vertices XLiFE++ uses.

 Moreover, the returned vector gives the location of the vertices according
 to the way they are encountered in the numbering process implemented by the
 functions GeomFigureMesh<T_>::createHOV, GeomFigureMesh<T_>::createHOeV, HexahedronMesh::createHOfV
 and HexahedronMesh::createHOiV.
 */
vector< vector<number_t> > Hexahedron::numberingOfVertices(const number_t Order) {
   number_t order = Order;
   if (order < 1) {order = 1;}

   number_t NbVert((order+1)*(order+1)*(order+1));
   vector<number_t> V(3,0);
   vector< vector<number_t> > BC(NbVert,V);

//    main vertices (order 1)
//    #1 -> (0, 0, 0), #2 -> (0, 0, k), #3 -> (0, k, 0), #4 -> (0, k, k)
//    #5 -> (k, 0, 0), #6 -> (k, 0, k), #7 -> (k, k, 0), #8 -> (k, k, k)
//    int rkVert[] = {3,7,2,6, 0,4,1,5}; // Vertices 4,8,3,7, 1,5,2,6 in Melina++
   const number_t ivx(0), ivy(1), ivz(2);
   number_t indBC(0);
   for (number_t ix=0; ix<2; ix++) {
      V[ivx] = ix*order;
      for (number_t iy=0; iy<2; iy++) {
         V[ivy] = iy*order;
         for (number_t iz=0; iz<2; iz++) {
            V[ivz] = iz*order;
//             BC[rkVert[indBC++]] = V;
            BC[indBC++] = V;
            }
         }
      }
// Order 2 : rkVert[] = {4, 12, 8, 14, 24, 17, 3, 18, 7, 13, 25, 19, 26, 27, 23, 16, 22, 10, 1, 15, 5, 20, 21, 11, 2, 9, 6}

/*
      edge vertices (cf. GeomFigureMesh<T_>::createHOeV)
      vertex numbering is decreasing along the edges
        1. 4 edges along Z : edge #1 =(1,2) -> (0, 0, .)
                             edge #2 =(3,4) -> (0, k, .)
                             edge #3 =(5,6) -> (k, 0, .)
                             edge #4 =(7,8) -> (k, k, .)
        2. 4 edges along Y : edge #5 =(1,3) -> (0, ., 0)
                             edge #6 =(2,4) -> (0, ., k)
                             edge #7 =(5,7) -> (k, ., 0)
                             edge #8 =(6,8) -> (k, ., k)
        3. 4 edges along X : edge #9 =(1,5) -> (., 0, 0)
                             edge #10=(2,6) -> (., 0, k)
                             edge #11=(3,7) -> (., k, 0)
                             edge #12=(4,8) -> (., k, k)
*/
   number_t tiv[]={2,1,0}; // Z, Y, X
   for (number_t indEdge=0; indEdge<nb_edges_; indEdge++) {
      V = BC[rkEdge[indEdge][0]];
      number_t iv = tiv[indEdge/4];
      for (number_t i1=1; i1<order; i1++) {
         V[iv] = order-i1;
         BC[indBC++] = V;
      }
   }
/*
      face vertices (cf. HexahedronMesh::createHOfV)
        1. 2 faces orthogonal to X : face #1=(1,3,4,2) -> (0, ., .)
                                     face #2=(5,7,8,6) -> (k, ., .)
                                                   increasing order along Y first then Z
        2. 2 faces orthogonal to Y : face #3=(1,2,6,5) -> (., 0, .)
                                     face #4=(3,4,8,7) -> (., k, .)
                                                   increasing order along Z first then X
        3. 2 faces orthogonal to Z : face #5=(1,5,7,3) -> (., ., 0)
                                     face #6=(2,6,8,4) -> (., ., k)
                                                   increasing order along X first then Y
*/
   number_t tiv1[]={2,0,1}, tiv2[]={1,2,0}; // couples (iv1,iv2) : (Z,Y), (X,Z) and (Y,X)
   for (number_t indFace=0; indFace<nb_faces_; indFace++) {
      V = BC[rkFace[indFace][0]];
      number_t iv1 = tiv1[indFace/2], iv2 = tiv2[indFace/2];
      for (number_t i1=1; i1<order; i1++) {
         V[iv1] = i1;
         for (number_t i2=1; i2<order; i2++) {
            V[iv2] = i2;
            BC[indBC++] = V;
         }
      }
   }
//    internal vertices (cf. HexahedronMesh::createHOiV)
//    increasing order along X first then Y then Z
   for (number_t iz=1; iz<order; iz++) {
      V[ivz] = iz;
      for (number_t iy=1; iy<order; iy++) {
         V[ivy] = iy;
         for (number_t ix=1; ix<order; ix++) {
            V[ivx] = ix;
            BC[indBC++] = V;
            }
         }
      }
   return BC;
}

} // end of namespace subdivision
} // end of namespace xlifepp