3232#include "libmesh/point.h"
3333
3434// C++ includes
35-
35+ #include <numeric> // gcd
3636
3737using namespace libMesh ;
3838
@@ -45,6 +45,7 @@ void usage_error(const char * progname)
4545 << " --elem type elem type (e.g. TET14)\n"
4646 << " --childnum num child number\n"
4747 << " --denominator num denominator to use\n"
48+ << " --diff only output diffs from existing\n"
4849 << std ::endl ;
4950
5051 exit (1 );
@@ -83,7 +84,9 @@ int main(int argc, char ** argv)
8384 assert_argument (cl , "--childnum" , argv [0 ], 0 );
8485
8586 const int denominator =
86- assert_argument (cl , "--denominator" , argv [0 ], 0 );
87+ assert_argument (cl , "--denominator" , argv [0 ], 1 );
88+
89+ const bool diff = cl .search ("--diff" );
8790
8891 const ElemType elem_type =
8992 Utility ::string_to_enum < ElemType > (elem_type_string );
@@ -131,26 +134,97 @@ int main(int argc, char ** argv)
131134 elem -> set_node (n ) = & nodes [n ];
132135 }
133136
134- for (auto i : elem -> node_index_range ())
137+ const unsigned int denomdigits = std ::ceil (std ::log10 (denominator ));
138+ const unsigned int spacing = denomdigits * 2 + 3 ;
139+
140+ std ::cout .precision (17 );
141+
142+ std ::cout << " // embedding matrix for child " << childnum << '\n' ;
143+ std ::cout << " {\n" ;
144+ std ::cout << " //" ;
145+ for (auto i : make_range (n_nodes ))
146+ {
147+ const unsigned int indexdigits =
148+ std ::ceil (std ::log10 (i ));
149+ const int padding = spacing - indexdigits - 2 * (i == 0 );
150+ for (int k = 0 ; k < padding ; ++ k )
151+ std ::cout << ' ' ;
152+ std ::cout << i ;
153+ if (i + 1 == n_nodes )
154+ std ::cout << '\n' ;
155+ else
156+ std ::cout << ',' ;
157+ }
158+
159+ for (auto i : make_range (n_nodes ))
135160 {
136161 const Point & pt = elem -> point (i );
137- std ::cout << '{' ;
162+ std ::cout << " {" ;
138163 for (auto j : make_range (n_nodes ))
139164 {
140165 Real shape = FEInterface ::shape (dim , fe_type , elem_type , j , pt );
166+
167+ // Don't print -0 or 1e-17; we don't tolerate FP error at 0
168+ if (std ::abs (shape ) < TOLERANCE * std ::sqrt (TOLERANCE ))
169+ shape = 0 ;
170+
141171 const Real embed = ref .embedding_matrix (childnum , i , j );
142- if (std ::abs (shape - embed ) < TOLERANCE * TOLERANCE )
143- std ::cout << "++++++," ;
172+ if (diff &&
173+ std ::abs (shape - embed ) < TOLERANCE * std ::sqrt (TOLERANCE ))
174+ {
175+ for (unsigned int k = 0 ; k != spacing ; ++ k )
176+ std ::cout << '+' ;
177+ if (j + 1 != n_nodes )
178+ std ::cout << ',' ;
179+ }
144180 else
145181 {
146- std ::cout << (shape * denominator );
147- if (shape != 0.0 )
148- std ::cout << "/r" << denominator ;
149- std ::cout << ", " ;
182+ int oldnumerator = std ::round (shape * denominator );
183+ int newnumerator = oldnumerator ;
184+ int newdenominator = denominator ;
185+
186+ if (std ::abs (oldnumerator - shape * denominator ) < TOLERANCE * sqrt (TOLERANCE ))
187+ {
188+ int the_gcd = std ::gcd (newnumerator , newdenominator );
189+ newnumerator /= the_gcd ;
190+ newdenominator /= the_gcd ;
191+ }
192+
193+ const unsigned int newdenomdigits =
194+ std ::ceil (std ::log10 (newdenominator ));
195+ std ::ostringstream ostr ;
196+ ostr << (shape * newdenominator );
197+ const int padding =
198+ (shape != 0.0 && newdenominator != 1 ) ?
199+ int (spacing )- newdenomdigits - 2 - ostr .str ().size () :
200+ int (spacing )- ostr .str ().size ();
201+ for (int k = 0 ; k < padding ; ++ k )
202+ std ::cout << ' ' ;
203+ std ::cout << ostr .str ();
204+ if (shape != 0.0 && newdenominator != 1 )
205+ {
206+ if (1 << (int )std ::round (std ::log2 (newdenominator )) ==
207+ newdenominator )
208+ std ::cout << "/" << newdenominator << '.' ;
209+ // If we don't have a power of 2 we need to make
210+ // sure we're dividing at might-exceed-double Real
211+ // precision
212+ else
213+ std ::cout << "/r" << newdenominator ;
214+ }
215+ if (j + 1 != n_nodes )
216+ std ::cout << ',' ;
150217 }
151218 }
152- std ::cout << '}' << std ::endl ;
219+
220+ if (i + 1 == n_nodes )
221+ std ::cout << "} " ;
222+ else
223+ std ::cout << "}, " ;
224+
225+ std ::cout << "// " << i << '\n' ;
153226 }
227+ std ::cout << " },\n" ;
154228
155229 return 0 ;
156230}
0 commit comments