diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp index c720228a455c28f3b1d6b08c55d7e8533e0320a5..8361c2f29cab7253f10c76c639016cad54d2fd73 100644 --- a/src/mesh/MeshNodeBoundary.hpp +++ b/src/mesh/MeshNodeBoundary.hpp @@ -211,6 +211,122 @@ _getNormal(const MeshType& mesh) return normal; } +template <> +template <typename MeshType> +inline TinyVector<3,double> +MeshFlatNodeBoundary<3>:: +_getNormal(const MeshType& mesh) +{ + static_assert(MeshType::dimension == 3); + typedef TinyVector<3,double> R3; + + + R3 xmin(std::numeric_limits<double>::max(), + std::numeric_limits<double>::max(), + std::numeric_limits<double>::max()); + R3 ymin = xmin; + R3 zmin = xmin;; + + R3 xmax =-xmin; + R3 ymax = xmax; + R3 zmax = xmax; + + const Kokkos::View<const R3*> xr = mesh.xr(); + + for (size_t r=0; r<m_node_list.extent(0); ++r) { + const R3& x = xr[m_node_list[r]]; + if (x[0]<xmin[0]) { + xmin = x; + } + if (x[1]<ymin[1]) { + ymin = x; + } + if (x[2]<zmin[2]) { + zmin = x; + } + if (x[0]>xmax[0]) { + xmax = x; + } + if (x[1]>ymax[1]) { + ymax = x; + } + if (x[2]>zmax[2]) { + zmax = x; + } + } + + const R3 u = xmax-xmin; + const R3 v = ymax-ymin; + const R3 w = zmax-zmin; + + const R3 uv = crossProduct(u,v); + const double uv_l2 = (uv, uv); + + R3 normal = uv; + double normal_l2 = uv_l2; + + const R3 uw = crossProduct(u,w); + const double uw_l2 = (uw,uw); + + if (uw_l2 > uv_l2) { + normal = uw; + normal_l2 = uw_l2; + } + + const R3 vw = crossProduct(v,w); + const double vw_l2 = (vw,vw); + + if (vw_l2 > normal_l2) { + normal = vw; + normal_l2 = vw_l2; + } + + + if (normal_l2 ==0) { + std::cerr << "not able to compute normal!\n"; + std::exit(1); + } + + normal *= 1./sqrt(normal_l2); + +#warning Add flatness test + // this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh); + + return normal; +} + +template <> +template <typename MeshType> +inline TinyVector<1,double> +MeshFlatNodeBoundary<1>:: +_getOutgoingNormal(const MeshType& mesh) +{ + static_assert(MeshType::dimension == 1); + typedef TinyVector<1,double> R; + + const R normal = this->_getNormal(mesh); + + const Kokkos::View<const R*>& xr = mesh.xr(); + const Kokkos::View<const unsigned int**>& node_cells = mesh.connectivity().nodeCells(); + const Kokkos::View<const unsigned int**>& cell_nodes = mesh.connectivity().cellNodes(); + const Kokkos::View<const unsigned short*>& cell_nb_nodes = mesh.connectivity().cellNbNodes(); + + const size_t r0 = m_node_list[0]; + const size_t j0 = node_cells(r0,0); + double max_height = 0; + for (int r=0; r<cell_nb_nodes(j0); ++r) { + const double height = (xr(cell_nodes(j0, r))-xr(r0), normal); + if (std::abs(height) > std::abs(max_height)) { + max_height = height; + } + } + if (max_height > 0) { + return -normal; + } else { + return normal; + } +} + template <> template <typename MeshType> inline TinyVector<2,double> @@ -245,16 +361,16 @@ _getOutgoingNormal(const MeshType& mesh) template <> template <typename MeshType> -inline TinyVector<1,double> -MeshFlatNodeBoundary<1>:: +inline TinyVector<3,double> +MeshFlatNodeBoundary<3>:: _getOutgoingNormal(const MeshType& mesh) { - static_assert(MeshType::dimension == 1); - typedef TinyVector<1,double> R; + static_assert(MeshType::dimension == 3); + typedef TinyVector<3,double> R3; - const R normal = this->_getNormal(mesh); + const R3 normal = this->_getNormal(mesh); - const Kokkos::View<const R*>& xr = mesh.xr(); + const Kokkos::View<const R3*>& xr = mesh.xr(); const Kokkos::View<const unsigned int**>& node_cells = mesh.connectivity().nodeCells(); const Kokkos::View<const unsigned int**>& cell_nodes = mesh.connectivity().cellNodes(); const Kokkos::View<const unsigned short*>& cell_nb_nodes = mesh.connectivity().cellNbNodes();