1 /* Copyright (c) 1997-2021
2 Ewgenij Gawrilow, Michael Joswig, and the polymake team
3 Technische Universität Berlin, Germany
4 https://polymake.org
5
6 This program is free software; you can redistribute it and/or modify it
7 under the terms of the GNU General Public License as published by the
8 Free Software Foundation; either version 2, or (at your option) any
9 later version: http://www.gnu.org/licenses/gpl.txt.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15 --------------------------------------------------------------------------------
16 */
17
18 #include "polymake/client.h"
19 #include "polymake/PowerSet.h"
20 #include "polymake/Array.h"
21 #include "polymake/Rational.h"
22 #include "polymake/Matrix.h"
23 #include "polymake/Vector.h"
24 #include "polymake/list"
25 #include "polymake/hash_set"
26 #include <sstream>
27
28 namespace polymake { namespace topaz {
29
stellar_subdivision(BigObject p_in,const Array<Set<Int>> & subd_faces,OptionSet options)30 BigObject stellar_subdivision(BigObject p_in, const Array<Set<Int>>& subd_faces, OptionSet options)
31 {
32 const bool is_PC= !p_in.isa("topaz::SimplicialComplex");
33
34 Array<Set<Int>> C_in = p_in.give(is_PC ? Str("TRIANGULATION.FACETS") : Str("FACETS"));
35 Int n_vert = p_in.give(is_PC ? Str("N_POINTS") : Str("N_VERTICES"));
36
37 // compute new complex
38 std::list<Set<Int>> C;
39
40 for (Int i = 0; i < subd_faces.size(); ++i) {
41 const Set<Int> F = subd_faces[i];
42 bool facet_found = false;
43 for (auto f=entire(C_in);!f.at_end(); ++f) {
44 Set<Int> facet = *f;
45
46 if (incl(F,facet) < 1) { // F contained in the facet
47 facet_found = true;
48 const Int size = facet.size();
49 facet -= F;
50 facet += n_vert+i;
51
52 // add new facets
53 for (auto s_it = entire(all_subsets_of_k(F, size-facet.size())); !s_it.at_end(); ++s_it)
54 C.push_back(facet + *s_it);
55
56 } else {
57 C.push_back(facet);
58 }
59 }
60
61 if (!facet_found) {
62 throw std::runtime_error("stellar_subdivision: Input does not specify a face (of the complex generated so far).");
63 }
64 C_in = Array<Set<Int>>(C.size(),C.begin());
65 C.clear();
66 }
67
68 BigObject p_out(p_in.type());
69 p_out.set_description()<<"Obtained from " << p_in.name() << " by barycentric subdivision of the faces\n" << subd_faces << ".\n";
70 p_out.take(is_PC ? Str("TRIANGULATION.FACETS") : Str("FACETS")) << as_array(C_in);
71
72 // compute new coordinates
73 if (is_PC) {
74 Matrix<Rational> Coord = p_in.give("POINTS");
75 Coord.resize(n_vert+subd_faces.size(),Coord.cols());
76 for (Int i = 0; i < subd_faces.size(); ++i)
77 Coord[i+n_vert] = average( rows(Coord.minor(subd_faces[i], All)) );
78
79 p_out.take("POINTS") << Coord;
80 }
81 else
82 if (options["geometric_realization"]) {
83 Matrix<Rational> Coord = p_in.give("COORDINATES");
84 Coord.resize(n_vert+subd_faces.size(),Coord.cols());
85 for (Int i = 0; i < subd_faces.size(); ++i)
86 Coord[i+n_vert] = average( rows(Coord.minor(subd_faces[i], All)) );
87
88 p_out.take("COORDINATES") << Coord;
89 }
90
91 // compute new label
92 if (!options["no_labels"]) {
93 Array<std::string> L = p_in.give(is_PC ? Str("LABELS") : Str("VERTEX_LABELS"));
94 hash_set<std::string> old_L(n_vert);
95 for (auto l=entire(L); !l.at_end(); ++l)
96 old_L.insert(*l);
97
98 L.resize(n_vert+subd_faces.size());
99 for (Int i = 0; i < subd_faces.size(); ++i) {
100 std::ostringstream label;
101 auto v=entire(subd_faces[i]);
102 label << "{" << L[*v]; ++v;
103 for ( ; !v.at_end(); ++v)
104 label << "," << L[*v];
105 label << "}";
106 std::string l=label.str(), ll=l;
107
108 // test if ll is unique
109 Int j = 0;
110 while (old_L.find(ll) != old_L.end()) {
111 ++j;
112 label.str("");
113 label << l << "_" << j;
114 ll=label.str();
115 }
116
117 L[n_vert+i] = ll;
118 }
119 p_out.take(is_PC ? Str("LABELS") : Str("VERTEX_LABELS")) << L;
120 }
121 return p_out;
122 }
123
124 UserFunction4perl("# @category Producing a new simplicial complex from others"
125 "# Computes the complex obtained by stellar subdivision of the given //faces// of the //complex//."
126 "# @param SimplicialComplex complex"
127 "# @param Array<Set<Int>> faces"
128 "# @option Bool no_labels Do not create [[VERTEX_LABELS]]. default: 0"
129 "# @option Bool geometric_realization default 0"
130 "# @return SimplicialComplex",
131 &stellar_subdivision, "stellar_subdivision($,Array<Set<Int> > { no_labels => 0, geometric_realization => 0})");
132
133 InsertEmbeddedRule("# @category Producing a new simplicial complex from others"
134 "# Computes the complex obtained by stellar subdivision of the given //face// of the //complex//."
135 "# @param SimplicialComplex complex"
136 "# @param Set<Int> face"
137 "# @option Bool no_labels Do not create [[VERTEX_LABELS]]. default: 0"
138 "# @option Bool geometric_realization default 0"
139 "# @return SimplicialComplex\n"
140 "user_function stellar_subdivision(SimplicialComplex, Set<Int> { no_labels => 0, geometric_realization => 0}) {\n"
141 " my $a=new Array<Set<Int> >(1);\n"
142 " my $p=shift;\n"
143 " $a->[0]=shift;\n"
144 "stellar_subdivision($p,$a,@_); }\n");
145 } }
146
147 // Local Variables:
148 // mode:C++
149 // c-basic-offset:3
150 // indent-tabs-mode:nil
151 // End:
152