ref_element.cc

Go to the documentation of this file.

/*!
 *
 * Copyright (C) 2015 Technical University of Liberec.  All rights reserved.
 * 
 * This program is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License version 3 as published by the
 * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
 * 
 * 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.
 *
 * 
 * @file    ref_element.cc
 * @brief   Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.
 * @author  Jan Stebel
 */
 
#include "system/global_defs.h"
#include "system/system.hh"
#include "mesh/ref_element.hh"
#include "fem/arena_vec.hh"
 
 
 
using namespace arma;
using namespace std;
 
constexpr IdxVector<3> _AuxInteract::S3_node_sides [2][4];
constexpr IdxVector<2> _AuxInteract::S3_line_sides [2][6];
constexpr IdxVector<3> _AuxInteract::S3_node_lines [2][4];
 
 
 
template<std::size_t n>
std::vector< std::vector<unsigned int> > _array_to_vec( const std::vector<IdxVector<n>> array_vec) {
    std::vector< std::vector<unsigned int> > vec(array_vec.size());
    for(unsigned int i=0; i<array_vec.size(); i++)
        for(unsigned int j=0;j<n; j++)
            vec[i].push_back(array_vec[i][j]);
    return vec;
}
 
 
 
template<> const std::vector<IdxVector<2>> RefElement<1>::line_nodes_
    = { {0,1} };
 
template<> const std::vector<IdxVector<2>> RefElement<2>::line_nodes_
    = { {0,1},
        {0,2},
        {1,2}};
 
template<> const std::vector<IdxVector<2>> RefElement<3>::line_nodes_
    = { {0,1},  //0
        {0,2},  //1
        {0,3},  //2 <-3 (fixed order)
        {1,2},  //3 <-2
        {1,3},  //4
        {2,3}}; //5
 
 
template<> const std::vector<IdxVector<1>> RefElement<1>::node_lines_
    = { {0},
        {0}};
 
template<> const std::vector<IdxVector<2>> RefElement<2>::node_lines_
    = { {1,0},
        {0,2},
        {2,1}};
 
 
 
 
 
// Lexicographic order.
template<> const std::vector<IdxVector<3>> RefElement<3>::side_nodes_
    = { { 0, 1, 2 },
        { 0, 1, 3 },
        { 0, 2, 3 },
        { 1, 2, 3 }};
 
//// Order clockwise, faces opposite to the lines from node_lines.
//// !! dependes on S3 inversion
//template<> const std::vector<IdxVector<3>> RefElement<3>::node_sides_
//    = { { { 2, 1, 0 },
//          { 3, 0, 1 },
//          { 3, 2, 0 },
//          { 3, 1, 2 }},
//        { { 2, 0, 1 },
//          { 3, 1, 0 },
//          { 3, 0, 2 },
//          { 3, 2, 1 }};
 
 
 
// lexicographic order
template<> const std::vector<IdxVector<3>> RefElement<3>::side_lines_
    = { {0,1,3},
        {0,2,4},
        {1,2,5},
        {3,4,5}};
 
 
template<> const IdxVector<1> RefElement<0>::topology_zeros_[] = {
   {(1 << 1)}  //node 0
};
 
 
template<> const IdxVector<2> RefElement<1>::topology_zeros_[] = {
   {(1 << 1),  //node 0
    (1 << 0)}, //node 1
   {0,         //the element
    0}
};
 
 
template<> const IdxVector<3> RefElement<2>::topology_zeros_[] = {
   {(1 << 1) | (1 << 2),  //node 0
    (1 << 0) | (1 << 2),  //node 1
    (1 << 0) | (1 << 1)}, //node 2
   {(1 << 2),  //line 0
    (1 << 1),  //line 1
    (1 << 0)}, //line 2
   {0,         //the element
    0,
    0}
};
 
template<> const IdxVector<6> RefElement<3>::topology_zeros_[] = {
   {(1 << 1) | (1 << 2) | (1 << 3),  //node 0
    (1 << 0) | (1 << 2) | (1 << 3),  //node 1
    (1 << 0) | (1 << 1) | (1 << 3),  //node 2
    (1 << 0) | (1 << 1) | (1 << 2),  //node 3
    0,
    0},
   {(1 << 2) | (1 << 3),    //line 0
    (1 << 1) | (1 << 3),    //line 1
    (1 << 1) | (1 << 2),    //line 2
    (1 << 0) | (1 << 3),    //line 3
    (1 << 0) | (1 << 2),    //line 4
    (1 << 0) | (1 << 1)},   //line 5
   {1 << 3,    //side 0
    1 << 2,    //side 1
    1 << 1,    //side 2
    1 << 0,    //side 3
    0,
    0},
   {0,    //the element
    0,
    0,
    0,
    0,
    0}
};
 
    
 
 
// 0: nodes of nodes
// 1: nodes of lines
// 2: nodes of sides
// 3: nodes of tetrahedron
template<> const std::vector< std::vector< std::vector<unsigned int> > > RefElement<0>::nodes_of_subelements = {
        { {0} }
};
 
template<> const std::vector< std::vector< std::vector<unsigned int> > > RefElement<1>::nodes_of_subelements = {
        { {0}, {1} },
        _array_to_vec(line_nodes_)
};
 
template<> const std::vector< std::vector< std::vector<unsigned int> > > RefElement<2>::nodes_of_subelements = {
        { {0}, {1}, {2} },
        _array_to_vec(line_nodes_),
        { {0,1,2} }
};
 
template<> const std::vector< std::vector< std::vector<unsigned int> > > RefElement<3>::nodes_of_subelements = {
        { {0}, {1}, {2}, {3} },
        _array_to_vec(line_nodes_),
        _array_to_vec(side_nodes_),
        { {0,1,2,3} }
};
 
 
template<unsigned int dim>
auto RefElement<dim>::local_to_bary(const LocalPoint& lp) -> BaryPoint
{
    ASSERT_EQ(lp.n_rows, dim);
    BaryPoint bp;
    bp.rows(1, dim ) = lp;
    bp( 0 ) = 1.0 - arma::sum(lp);
    return bp;
 
    // new armadillo
    // return arma::join_col( arma::vec::fixed<1>( { 1.0 - arma::sum( local )} ), local);
 
}
 
template<unsigned int dim>
auto RefElement<dim>::bary_to_local(const BaryPoint& bp) -> LocalPoint
{
    ASSERT_EQ(bp.n_rows, dim+1);
    LocalPoint lp = bp.rows(1, dim);
    return lp;
}
 
template<unsigned int dim>
inline unsigned int RefElement<dim>::oposite_node(unsigned int sid)
{
    return n_sides - sid - 1;
}
 
 
template<unsigned int dim>
unsigned int RefElement<dim>::normal_orientation(unsigned int sid)
{
    ASSERT_LT(sid, n_sides).error("Side number is out of range!");
 
    return sid % 2;
}
 
 
template<>
vec::fixed<1> RefElement<1>::normal_vector(unsigned int sid)
{
    ASSERT_LT(sid, n_sides).error("Side number is out of range!");
 
    return node_coords(sid) - node_coords(1-sid);
}
 
template<>
vec::fixed<2> RefElement<2>::normal_vector(unsigned int sid)
{
    ASSERT_LT(sid, n_sides).error("Side number is out of range!");
    vec::fixed<2> barycenter, bar_side, n, t;
 
    // tangent vector along line
    t = node_coords(line_nodes_[sid][1]) - node_coords(line_nodes_[sid][0]);
    // barycenter coordinates
    barycenter.fill(1./3);
    // vector from barycenter to the side
    bar_side = node_coords(line_nodes_[sid][0]) - barycenter;
    // normal vector to side (modulo sign)
    n(0) = -t(1);
    n(1) = t(0);
    n /= norm(n,2);
    // check sign of normal vector
    if (dot(n,bar_side) < 0) n *= -1;
 
    return n;
}
 
template<>
vec::fixed<3> RefElement<3>::normal_vector(unsigned int sid)
{
    ASSERT_LT(sid, n_sides).error("Side number is out of range!");
    vec::fixed<3> barycenter, bar_side, n, t1, t2;
 
    // tangent vectors of side
    t1 = node_coords(side_nodes_[sid][1]) - node_coords(side_nodes_[sid][0]);
    t2 = node_coords(side_nodes_[sid][2]) - node_coords(side_nodes_[sid][0]);
    // baryucenter coordinates
    barycenter.fill(0.25);
    // vector from barycenter to the side
    bar_side = node_coords(side_nodes_[sid][0]) - barycenter;
    // normal vector (modulo sign)
    n = cross(t1,t2);
    n /= norm(n,2);
    // check sign of normal vector
    if (dot(n,bar_side) < 0) n = -n;
 
    return n;
}
 
 
 
template<unsigned int dim>
auto RefElement<dim>::normal_vector_array(ArenaVec<uint> loc_side_idx_vec) -> Eigen::Matrix<ArenaVec<double>, dim, 1>
{
    Eigen::Matrix<ArenaVec<double>, dim, 1> normal_vec_array;
    for (uint i=0; i<dim; ++i)
        normal_vec_array(i) = ArenaVec<double>( loc_side_idx_vec.data_size(), loc_side_idx_vec.arena() );
    for (uint sid=0; sid<loc_side_idx_vec.data_size(); ++sid) {
        arma::vec::fixed<dim> n_vec = RefElement<dim>::normal_vector( loc_side_idx_vec(sid) );
        for (uint i=0; i<dim; ++i) {
            normal_vec_array(i)(sid) = n_vec(i);
        }
    }
    return normal_vec_array;
}
 
 
 
template<>
auto RefElement<0>::normal_vector_array(FMT_UNUSED ArenaVec<uint> loc_side_idx_array) -> Eigen::Matrix<ArenaVec<double>, 0, 1>
{
    Eigen::Matrix<ArenaVec<double>, 0, 1> normal_vec_array;
    return normal_vec_array;
}
 
 
 
template<unsigned int dim>
auto RefElement<dim>::barycentric_on_face(const BaryPoint &barycentric, unsigned int i_face) -> FaceBaryPoint
{
    ASSERT_EQ(barycentric.n_rows, dim+1);
    FaceBaryPoint face_barycentric;
    for(unsigned int i=0; i < dim; i++) {
//        unsigned int i_sub_node = (i+1)%dim;
//         unsigned int i_bary = (dim + side_nodes_[i_face][i_sub_node])%(dim+1);
        unsigned int i_bary = interact_<0,dim-1>(i_face)[i];
        face_barycentric[i] = barycentric[ i_bary ];
    }
    return face_barycentric;
}
 
 
template<unsigned int dim>
std::pair<unsigned int, unsigned int> RefElement<dim>::zeros_positions(const BaryPoint &barycentric,
                                                                       double tolerance)
{
    unsigned int zeros = 0;
    unsigned int n_zeros = 0;
    for(unsigned int i=0; i < dim+1; i++){
        if(std::fabs(barycentric[i]) < tolerance)
        {
            zeros = zeros | (1 << i);
            n_zeros++;
        }
    }
    
    return std::make_pair(n_zeros, zeros);
}
 
 
template<>
auto RefElement<0>::clip(const BaryPoint &barycentric) -> BaryPoint
{
    return barycentric;
}
 
template<unsigned int dim>
auto  RefElement<dim>::make_bary_unit_vec()->BarycentricUnitVec
{
    std::vector<arma::vec::fixed<dim+1> > bary_unit_vec(dim+1, arma::zeros(dim+1));
    for(unsigned int i=0; i<dim; i++) {
        bary_unit_vec[i][i] = 1.0;
        bary_unit_vec[i][dim] = -1.0;
        bary_unit_vec[dim][i] = -1.0 / dim;
    }
    bary_unit_vec[dim][dim] = 1.0;
    return bary_unit_vec;
}
 
 
template<unsigned int dim>
auto RefElement<dim>::clip(const BaryPoint &barycentric) -> BaryPoint
{
    static BarycentricUnitVec bary_unit_vec = make_bary_unit_vec();
    ASSERT_EQ(barycentric.n_rows, dim+1);
    for(unsigned int i_bary=0; i_bary < dim +1; i_bary ++) {
        if (barycentric[i_bary] < 0.0) {
            // index of barycentric coord that is constant on the face i_side
            // as we use barycentric coords starting with local coordinates:
            // TODO: rather work only with local coords and/or with canonical barycentric coords
            unsigned int i_side = (dim - i_bary);
            // project to face
            arma::vec projection_to_face(dim+1);
            //barycentric.print(cout, "input");
            //cout << "is: " << i_side << endl;
            //cout << "ibary: " << i_bary << endl;
            //bary_unit_vec[i_bary].print(cout, "normal");
            //barycentric.subvec(0, dim-1).print(cout, "bary sub");
            projection_to_face = barycentric - barycentric[i_bary]*bary_unit_vec[i_bary];
            //projection_to_face(dim) = 1.0 - arma::sum(projection_to_face.subvec(0, dim-1));
            //projection_to_face.print(cout, "projection");
            auto bary_on_face = barycentric_on_face(projection_to_face, i_side);
            //bary_on_face.print(cout, "b on f");
            auto sub_clip = RefElement<dim-1>::clip(bary_on_face);
            //sub_clip.print(cout, "sub clip");
            return interpolate<dim-1>(sub_clip, i_side);
        }
    }
    return barycentric;
 
}
 
 
 
template<unsigned int dim>
auto RefElement<dim>::centers_of_subelements(unsigned int sub_dim)->CentersList
{
    static std::vector< std::vector<LocalPoint> > list;
    if (list.size() == 0) {
        list.resize(dim+1);
        for(unsigned int sdim=0; sdim < dim+1; sdim++) {
            // Temporary solution until we have unified interface to
            // the internal indexing.
            // We use the fact that numbering of subelements goes as numbering of
            // k combinations over nodes.
//             std::vector<unsigned int> subel_comb(sdim+2);
            for(auto &sub_el_nodes : nodes_of_subelements[sdim]) {
                ASSERT_EQ(sub_el_nodes.size(), sdim+1);
                LocalPoint center = arma::zeros(dim);
                for( unsigned int i_node : sub_el_nodes)
                    center+=node_coords( i_node );
                center/=(sdim+1);
                list[sdim].push_back(center);
            }
        }
    }
 
    ASSERT_LE(sub_dim, dim);
    return list[sub_dim];
}
 
 
template<>
double RefElement<1>::side_measure(unsigned int sid)
{
    ASSERT_LT(sid, n_sides).error("Side number is out of range!");
 
    return 1;
}
 
 
template<>
double RefElement<2>::side_measure(unsigned int sid)
{
    ASSERT_LT(sid, n_sides).error("Side number is out of range!");
 
    return norm(node_coords(line_nodes_[sid][1]) - node_coords(line_nodes_[sid][0]),2);
}
 
 
template<>
double RefElement<3>::side_measure(unsigned int sid)
{
    ASSERT_LT(sid, n_sides).error("Side number is out of range!");
 
    return 0.5*norm(cross(node_coords(side_nodes_[sid][1]) - node_coords(side_nodes_[sid][0]),
            node_coords(side_nodes_[sid][2]) - node_coords(side_nodes_[sid][0])),2);
}
 
 
 
template <>
unsigned int RefElement<3>::line_between_faces(unsigned int f1, unsigned int f2) {
    unsigned int i,j;
    i=j=0;
    while (side_lines_[f1][i] != side_lines_[f2][j])
        if (side_lines_[f1][i] < side_lines_[f2][j]) i++;
        else j++;
    return side_lines_[f1][i];
}
 
 
/**
     * Basic line interpolation
     */
template<unsigned int dim>
arma::vec::fixed<dim+1> RefElement<dim>::line_barycentric_interpolation(
                                                       arma::vec::fixed<dim+1> first_coords, 
                                                       arma::vec::fixed<dim+1> second_coords, 
                                                       double first_theta, double second_theta, double theta){
 
    arma::vec::fixed<dim+1> bary_interpolated_coords;
    bary_interpolated_coords = ((theta - first_theta) * second_coords + (second_theta - theta) * first_coords)
                               /(second_theta - first_theta);
    return bary_interpolated_coords;
}
 
template class RefElement<0>;
template class RefElement<1>;
template class RefElement<2>;
template class RefElement<3>;