polyhedron.cpp

Go to the documentation of this file.

/*****************************************************************************
    
    polyhedron.cpp -- See declarations in polyhedron.hpp.

    This file is a part of the Arageli library.

    Copyright (C) Sergey S. Lyalin, 2006
    University of Nizhni Novgorod, Russia, 2005--2006

*****************************************************************************/

#include "config.hpp"

#if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || \
    defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_polyhedron)

// REFERENCE ADDITIONAL HEADERS HERE

#include "frwrddecl.hpp"
#include "polyhedron.hpp"
#include "vector.hpp"
#include "matrix.hpp"


namespace Arageli
{


template <typename T, typename M, typename C, typename CFG>
template <typename D>
polyhedron<T, M, C, CFG>::polyhedron (const D& dim, const fromspace_t&)
{
    cone_m.inequaton_matrix().assign_eye(1, dim + 1);   // one inequation x[0] >= 0
    cone_m.equation_matrix().assign(0, dim + 1);
    set_normal_implicit();
}


template <typename T, typename M, typename C, typename CFG>
template <typename M1>
polyhedron<T, M, C, CFG>::polyhedron (const M1& ineqmat, const fromineq_t&)
:   cone_m(ineqmat, fromineq)
{
    ARAGELI_ASSERT_0(cone_m.inequation_matrix().ncols() > 0);

    // Add one additional inequation to cut all points with x[0] < 0.
    // May be it's superfluous, but the determination of superfluity is expensive.
    
    cone_m.inequation_matrix().insert_row
    (
        cone_m.inequation_matrix().nrows(),
        Arageli::null<inequation_element_type>()
    );
    
    cone_m.inequation_matrix(cone_m.inequation_matrix().nrows(), 0) =
        unit<inequation_element_type>();
}

    
template <typename T, typename M, typename C, typename CFG>
template <typename M1>
polyhedron<T, M, C, CFG>::polyhedron (const M1& vert, const fromivert_t&)
{
    cone_m.generatrix_matrix() = vert;
    cone_m.generatrix_matrix().insert_col
        (0, vector<T, false>(vert.nrows(), unit<T>(), fromval));
    cone_m.basis_matrix().resize(0, cone_m.generatrix_matrix().ncols());
    output_aligned(std::cout << "\ncone_m.generatrix_matrix() = \n", cone_m.generatrix_matrix());
}


template <typename T, typename M, typename C, typename CFG>
template <typename M1>
polyhedron<T, M, C, CFG>::polyhedron (const M1& vert, const fromvert_t&)
{
    // WARNING! This is too expensive!
    // TODO: Make it faster!

    vector<T, false> b;
    b.reserve(vert.nrows());
    cone_m.generatrix_matrix().resize(0, vert.ncols());

    for(size_type i = 0; i < vert.nrows(); ++i)
    {
        T lcmval = unit<T>();
        for(size_type j = 0; j < vert.ncols(); ++j)
            lcmval = lcm(lcmval, vert(i, j).denominator());

        cone_m.generatrix_matrix().insert_row
            (cone_m.generatrix_matrix().nrows(), lcmval*vert.copy_row(i));
        b.push_back(lcmval);
    }

    cone_m.generatrix_matrix().insert_col(0, b);
    cone_m.basis_matrix().resize(0, cone_m.generatrix_matrix().ncols());
}


template <typename T, typename M, typename C, typename CFG>
sideset polyhedron<T, M, C, CFG>::sides () const
{
    matrix<T> a = cone_m.inequation_matrix(), f, q, e;
    skeleton(a, f, q, e);

    //output_aligned(std::cout << "\nq = \n", q);

    // we are intrested only in q matrix

    vector<vector<std::size_t> > res;

    for(std::size_t i = 0; i < q.ncols(); ++i)
    {
        vector<std::size_t> curres;

        for(std::size_t j = 0; j < q.nrows(); ++j)
            if(is_null(q(j, i)))curres.push_back(j);

        ARAGELI_ASSERT_1(!curres.is_empty());
        
        res.push_back(curres);
    }

    return sideset(res);
}




template
<
    typename Out, 
    typename T,
    typename R,
    typename M,
    typename CFG
>
void output_vrml (Out& out, const polyhedron<T, R, M, CFG>& p)
{
    ARAGELI_ASSERT_0(p.space_dim() == 3);

    typedef polyhedron<T, R, M, CFG> Polyhedra;
    typedef sideset Side_set;
    typedef typename Side_set::vertex_indices_set Vertex_set;
    
    // VRML header
    out <<
        "#VRML V2.0 utf8\n"
        "#This file is automatically generated by Arageli library.\n"
        "\n"
        "Shape { appearance Appearance { material Material { diffuseColor 1 1 1 } }"
        "geometry IndexedFaceSet { solid TRUE coord Coordinate { point [\n";

    output_aligned(out, p.template vertices<matrix<double> >(), "", ",", " ");

    out << "] } coordIndex [\n";

    Side_set sides = p.sides();

    for(typename Side_set::const_iterator i = sides.begin(); i != sides.end(); ++i)
    {
        Vertex_set vi = *i;
        for
        (
            typename Vertex_set::const_iterator j = vi.begin();
            j != vi.end(); ++j
        )
            out << *j << ", ";
        out << "-1,\n";
    }

    out << "] } }\n";
}


//template
//<
//  typename Out,
//  typename P,
//  typename X1, typename Y1, typename X2, typename Y2,
//  typename Viewdir,
//  typename Colormap
//>
//void output_polytope_pstricks_3d
//(
//  Out& out,
//  const P& p,
//  double x1,
//  double y1,
//  double x2,
//  double y2,
//  double linewidth,
//  const Viewdir& viewdir,
//  const Colormap& colormap
//)
//{
//  
//
//}


}


#else   // #if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || ...


namespace Arageli
{

// PLACE ALL NOT TEMPLATE NOT INLINE IMPLEMENTATIONS HERE

}


#endif  // #if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || ...

---