26#include "rheolef/iorheo.h"
35basis_rep<T>::put (ostream& os, reference_element hat_K)
const
40 default:
error_macro (
"unimplemented gnuplot visualization of "<<valued()<<
"-valued polynomials");
54 bool verbose = iorheo::getverbose(os);
55 bool clean = iorheo::getclean(os);
56 bool execute = iorheo::getexecute(os);
57 size_type nsub = iorheo::getsubdivide(os);
58 if (nsub <= 1) nsub = 20;
61 string basename =
"basis-" + name() +
"-" + hat_K.
name();
66 string plot_name = basename +
".plot";
67 string gdat_name = basename +
".gdat";
68 ofstream plot (plot_name.c_str());
69 cerr <<
"! file \"" << plot_name <<
"\" created" << endl;
70 filelist +=
" \"" + plot_name +
"\"";
73 plot <<
"gdat = \"" << gdat_name <<
"\"" << endl
74 <<
"set colors classic" << endl
75 <<
"set size square" << endl
83 for (
size_type loc_isid = 0, loc_nsid = hat_K.
n_subgeo(1); loc_isid < loc_nsid; ++loc_isid) {
88 plot <<
"set arrow from ";
92 if (mu+1 != d2) plot <<
", ";
93 xmin[mu] = min(xmin[mu],
T(a[mu]));
94 xmax[mu] = max(xmax[mu],
T(a[mu]));
99 if (mu+1 != d2) plot <<
", ";
100 xmin[mu] = min(xmin[mu],
T(b[mu]));
101 xmax[mu] = max(xmax[mu],
T(b[mu]));
103 plot <<
" nohead lc 0 lw 1" << endl;
107 <<
"set xrange ["<<xmin[0]-
delta<<
":"<<xmax[0]+
delta<<
"]"<<endl;
109 plot <<
"set yrange ["<<xmin[1]-
delta<<
":"<<xmax[1]+
delta<<
"]"<<endl;
111 plot <<
"set yrange [-1:1]"<<endl;
114 plot <<
"set zrange ["<<xmin[2]-
delta<<
":"<<xmax[2]+
delta<<
"]"<<endl
115 <<
"set xyplane at "<< xmin[2]-
delta << endl;
118 plot <<
"plot \\" << endl
120 for (
size_t loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
121 plot <<
" gdat u 1:"<<loc_idof+2 <<
" t \"L" << loc_idof+1 <<
"\" w l";
122 if (loc_idof+1 != loc_ndof) { plot <<
",\\"; }
126 plot <<
"pause_duration = 0.5" << endl
128 <<
"scale= 1./" << nsub << endl
130 <<
"idxmax = " << loc_ndof << endl
131 <<
"plot gdat i idx u 1:2:(scale*$3):(scale*$4) t sprintf(\"P%d\",idx+1) w vec lc 0 lw 2" << endl
132 <<
"load \"" << basename <<
".loop\"" << endl
134 string loop_name = basename +
".loop";
135 ofstream loop (loop_name.c_str());
136 cerr <<
"! file \"" << loop_name <<
"\" created" << endl;
137 filelist +=
" \"" + loop_name +
"\"";
138 loop <<
"idx = idx+1" << endl
139 <<
"if (idx < idxmax) \\" << endl
140 <<
" pause -1 \"<return>\"; \\" << endl
141 <<
" replot; \\" << endl
145 plot <<
"pause -1 \"<return>\"" << endl;
151 ofstream gdat (gdat_name.c_str());
152 cerr <<
"! file \"" << gdat_name <<
"\" created" << endl;
153 filelist +=
" \"" + gdat_name +
"\"";
154 gdat << setprecision(std::numeric_limits<T>::digits10)
155 <<
"# basis " << name() << endl
156 <<
"# element " << hat_K.
name() << endl
157 <<
"# degree " << degree() << endl
165 gdat <<
"# ndof " << loc_ndof << endl
167 for (
size_t loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
168 gdat <<
" L" << loc_idof+1;
173 eval (hat_K, hat_x, value);
175 for (
size_t loc_idof = 0, loc_ndof = value.size(); loc_idof < loc_ndof; ++loc_idof) {
176 for (
size_t mu = 0; mu <
d; ++mu) {
177 gdat <<
" " << value[loc_idof][mu];
187 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> hat_node;
188 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,Eigen::Dynamic> value;
191 size_t loc_nnod = value.rows();
192 size_t loc_ndof = value.cols();
193 for (
size_t loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
194 for (
size_t loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
195 for (
size_t mu = 0; mu <
d; ++mu) {
196 gdat << hat_node[loc_inod][mu] <<
" ";
198 for (
size_t mu = 0; mu <
d; ++mu) {
199 gdat << value(loc_inod,loc_idof)[mu] <<
" ";
203 gdat << endl << endl;
210 for (
size_type j = 0; j+k <= nsub; j++) {
211 for (
size_type i = 0; i+j+k <= nsub; i++) {
213 eval_lagrange (hat_x, hat_K, degree, sopt, inv_vdm, value);
214 Lambda = max(Lambda,
norm(value,1));
221 for (
size_type k = 0; k <= degree; k++) {
222 for (
size_type j = 0; j <= degree; j++) {
223 for (
size_type i = 0; i+j <= degree; i++) {
225 eval_lagrange (hat_x, hat_K, degree, sopt, inv_vdm, value);
226 Lambda = max(Lambda,
norm(value,1));
233 for (
size_type k = 0; k <= degree; k++) {
234 for (
size_type j = 0; j <= degree; j++) {
235 for (
size_type i = 0; i <= degree; i++) {
237 eval_lagrange (hat_x, hat_K, degree, sopt, inv_vdm, value);
238 Lambda = max(Lambda,
norm(value,1));
251 string command =
"gnuplot \"" + plot_name +
"\"";
252 if (verbose) clog <<
"! " << command << endl;
253 int status = system (command.c_str());
258 string command =
"/bin/rm -f " + filelist;
259 if (verbose) clog <<
"! " << command << endl;
260 int status = system (command.c_str());
273 plot <<
"set colors classic" << endl
274 <<
"set nokey" << endl
275 <<
"set noborder" << endl
276 <<
"set noxtics" << endl
277 <<
"set noytics" << endl
278 <<
"set noztics" << endl
281 plot <<
"set size ratio -1 # 2d equal scales" << endl;
283 plot <<
"set view 70, 140 # x to left, y to right, z to top" << endl;
284 plot <<
"set view equal xyz # 3d equal scales" << endl;
289 point_basic<T> xmin, xmax;
294 for (
size_type loc_isid = 0, loc_nsid = hat_K.
n_subgeo(1); loc_isid < loc_nsid; ++loc_isid) {
299 plot <<
"set arrow from ";
303 if (mu+1 != d2) plot <<
", ";
304 xmin[
mu] = min(xmin[mu],
T(a[mu]));
305 xmax[
mu] = max(xmax[mu],
T(a[mu]));
310 if (mu+1 != d2) plot <<
", ";
311 xmin[
mu] = min(xmin[mu],
T(b[mu]));
312 xmax[
mu] = max(xmax[mu],
T(b[mu]));
314 plot <<
" nohead lc 0 lw 1" << endl;
318 <<
"set xrange ["<<xmin[0]-
delta<<
":"<<xmax[0]+
delta<<
"]"<<endl;
320 plot <<
"set yrange ["<<xmin[1]-
delta<<
":"<<xmax[1]+
delta<<
"]"<<endl;
322 plot <<
"set yrange [-1:1]"<<endl;
325 plot <<
"set zrange ["<<xmin[2]-
delta<<
":"<<xmax[2]+
delta<<
"]"<<endl
326 <<
"set xyplane at "<< xmin[2]-
delta << endl;
332 point_basic<T> operator() (
const point_basic<T>& x)
const {
333 return is_continuous ? x : x -
a*(x-xc);
335 shrink_fun (
bool ic,
const point_basic<T>& xc1)
336 : is_continuous(ic), xc(xc1),
a(0.15) {}
346 bool verbose = iorheo::getverbose(os);
347 bool clean = iorheo::getclean(os);
348 bool execute = iorheo::getexecute(os);
350 _initialize_data_guard (hat_K);
351 const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& hat_x = hat_node(hat_K);
354 string basename =
"basis-" + name() +
"-" + hat_K.
name();
359 string gdat_name = basename +
".gdat";
360 ofstream gdat (gdat_name.c_str());
361 cerr <<
"! file \"" << gdat_name <<
"\" created" << endl;
362 filelist +=
" \"" + gdat_name +
"\"";
363 gdat << setprecision(std::numeric_limits<T>::digits10)
364 <<
"# basis " << name() << endl
365 <<
"# element " << hat_K.
name() << endl
366 <<
"# degree " << degree() << endl
367 <<
"# ndof " << ndof(hat_K) << endl
368 <<
"# raw_poly " << _sopt.get_raw_polynomial_name() << endl
369 <<
"# nodeset " << _sopt.get_node_name() << endl
370 <<
"# nnode " << hat_x.size() << endl
374 shrink_fun<T> shrink (option().is_continuous(), xc);
375 if (! option().is_trace_n()) {
377 for (
size_type loc_inod = first_inod_by_dimension(hat_K,dim), loc_nnod = first_inod_by_dimension(hat_K,dim+1); loc_inod < loc_nnod; ++loc_inod) {
380 gdat <<
" " << loc_inod << endl;
382 gdat << endl << endl;
386 size_type loc_inod = 0, loc_nnod = nnod(hat_K);
388 for (
size_type loc_isid = 0, loc_nsid = hat_K.
n_side(); loc_isid < loc_nsid; ++loc_isid) {
393 for (
size_type loc_jsidvert = 0; loc_jsidvert < loc_nsidvert; ++loc_jsidvert) {
395 xs += hat_K.
vertex (iloc);
397 xs /=
T(loc_nsidvert);
399 shrink_fun<T> shrink_s (option().is_continuous(), xs);
401 for (
size_type ins = 0; ins < nns; ++ins, ++loc_inod) {
402 check_macro (loc_inod < loc_nnod,
"invalid loc_inod");
403 point_basic<T> shrink_hat_x = shrink_s (shrink (hat_x [loc_inod]));
405 gdat <<
" " << loc_inod << endl;
412 string plot_name = basename +
".plot";
413 ofstream plot (plot_name.c_str());
414 cerr <<
"! file \"" << plot_name <<
"\" created" << endl;
415 filelist +=
" \"" + plot_name +
"\"";
416 plot <<
"gdat = \"" << gdat_name <<
"\"" << endl
417 <<
"set title \"" << name() <<
": nodes\"" << endl
419 put_hat_node_header<T> (plot, hat_K) ;
420 string color[4] = {
"#ff0000",
"#008800",
"#0000ff",
"#ff00ff"};
421 bool need_coma =
false;
423 plot << ((
d==3) ?
"s" :
"") <<
"plot \\" << endl;
424 if (! option().is_trace_n()) {
426 size_type nnod_dim = first_inod_by_dimension (hat_K, dim+1) - first_inod_by_dimension (hat_K, dim);
427 if (nnod_dim == 0)
continue;
428 plot <<
" " << (need_coma ?
", " :
" ") <<
"gdat \\" << endl
429 <<
" index " << index <<
" \\" << endl
430 <<
" w p ps 2 pt 7 lc rgb \"" << color[dim] <<
"\",\\" << endl
431 <<
" gdat \\" << endl
432 <<
" index " << index <<
" \\" << endl
433 <<
" u 1:2:3" << ((
d==3) ?
":4" :
"") <<
" w labels offset char 2 lc 0 \\" << endl
439 plot <<
" gdat \\" << endl
440 <<
" w p ps 2 pt 7 lc rgb \"" << color[0] <<
"\",\\" << endl
441 <<
" gdat \\" << endl
442 <<
" u 1:2:3" << ((
d==3) ?
":4" :
"") <<
" w labels offset char 2 lc 0" << endl
446 <<
"pause -1 \"<return>\"" << endl
451 string command =
"gnuplot \"" + plot_name +
"\"";
452 if (verbose) clog <<
"! " << command << endl;
453 int status = system (command.c_str());
458 string command =
"/bin/rm -f " + filelist;
459 if (verbose) clog <<
"! " << command << endl;
460 int status = system (command.c_str());
474 check_macro (is_nodal(),
"side-node visualization: only nodal basis supported");
475 Eigen::Matrix<size_type,Eigen::Dynamic,1> loc_inod;
476 local_idof_on_side (hat_K, sid, loc_inod);
478 bool verbose = iorheo::getverbose(os);
479 bool clean = iorheo::getclean(os);
480 bool execute = iorheo::getexecute(os);
482 _initialize_data_guard (hat_K);
483 const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& hat_x = hat_node (hat_K);
486 string basename =
"basis-" + name() +
"-" + hat_K.
name();
491 string gdat_name = basename +
".gdat";
492 ofstream gdat (gdat_name.c_str());
493 cerr <<
"! file \"" << gdat_name <<
"\" created" << endl;
494 filelist +=
" \"" + gdat_name +
"\"";
495 gdat << setprecision(std::numeric_limits<T>::digits10)
496 <<
"# basis " << name() << endl
497 <<
"# element " << hat_K.
name() << endl
498 <<
"# degree " << degree() << endl
499 <<
"# ndof " << ndof(hat_K) << endl
500 <<
"# raw_poly " << _sopt.get_raw_polynomial_name() << endl
501 <<
"# nodeset " << _sopt.get_node_name() << endl
502 <<
"# nnode " << hat_x.size() << endl
504 for (
size_type loc_sid_inod = 0, loc_sid_nnod = loc_inod.size(); loc_sid_inod < loc_sid_nnod; ++loc_sid_inod) {
505 hat_x [loc_inod[loc_sid_inod]].put (gdat, max(
d,
size_type(2)));
506 gdat <<
" " << loc_inod[loc_sid_inod] << endl;
511 string plot_name = basename +
".plot";
512 ofstream plot (plot_name.c_str());
513 cerr <<
"! file \"" << plot_name <<
"\" created" << endl;
514 filelist +=
" \"" + plot_name +
"\"";
515 plot <<
"gdat = \"" << gdat_name <<
"\"" << endl
516 <<
"set title \"" << name() <<
": nodes\"" << endl
518 put_hat_node_header<T> (plot, hat_K) ;
519 string color[4] = {
"#ff0000",
"#008800",
"#0000ff",
"#ff00ff"};
520 bool need_coma =
false;
522 plot << ((
d==3) ?
"s" :
"") <<
"plot \\" << endl
523 <<
" gdat \\" << endl
524 <<
" w p ps 2 pt 7 lc rgb \"" << color[
d] <<
"\",\\" << endl
525 <<
" gdat \\" << endl
526 <<
" u 1:2:3" << ((
d==3) ?
":4" :
"") <<
" w labels offset char 2 lc 0" << endl
527 <<
"pause -1 \"<return>\"" << endl
532 string command =
"gnuplot \"" + plot_name +
"\"";
533 if (verbose) clog <<
"! " << command << endl;
534 int status = system (command.c_str());
539 string command =
"/bin/rm -f " + filelist;
540 if (verbose) clog <<
"! " << command << endl;
541 int status = system (command.c_str());
547#define _RHEOLEF_instanciation(T) \
548template void basis_rep<Float>::put (ostream& os, reference_element hat_K) const; \
549template void basis_rep<Float>::put_scalar_valued (ostream& os, reference_element hat_K) const; \
550template void basis_rep<Float>::put_vector_valued (ostream& os, reference_element hat_K) const; \
551template void basis_rep<Float>::put_hat_node (ostream& os, reference_element hat_K) const; \
552template void basis_rep<Float>::put_hat_node_on_side ( \
553 std::ostream&, reference_element, const side_information_type&) const;
#define _RHEOLEF_instanciation(T, M, A)
basis - finite element method
field::size_type size_type
see the Float page for the full documentation
see the point page for the full documentation
reference_element::size_type size_type
void put_hat_node_on_side(std::ostream &os, reference_element hat_K, const side_information_type &sid) const
point_basic< T > virtual _RHEOLEF_compute_dofs("tensor", tensor_basic< T >) void put(std void put_scalar_valued(std::ostream &os, reference_element hat_K) const
void put_hat_node(std::ostream &os, reference_element hat_K) const
virtual void put_vector_valued(std::ostream &os, reference_element hat_K) const
std::ostream & put(std::ostream &s, int d=3) const
see the reference_element page for the full documentation
const point_basic< Float > & vertex(size_type iloc) const
static const variant_type H
static const variant_type q
static const variant_type e
reference_element side(size_type loc_isid) const
size_type dimension() const
size_type subgeo_size(size_type subgeo_dim, size_type loc_isid) const
static const variant_type p
variant_type variant() const
size_type subgeo_local_vertex(size_type subgeo_dim, size_type loc_isid, size_type loc_jsidvert) const
std::vector< int >::size_type size_type
size_type n_subgeo(size_type subgeo_dim) const
static const variant_type T
static const variant_type P
static const variant_type t
#define error_macro(message)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
Float delta(Float f, Float g)
void basis_on_pointset_evaluate(const Basis &b, const reference_element &hat_K, const Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_x, Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &vdm)
void put(const Basis &b, ostream &os, reference_element hat_K)
This file is part of Rheolef.
void reference_element_barycenter(reference_element hat_K, point_basic< T > &c)
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
void pointset_lagrange_equispaced(reference_element hat_K, size_t order_in, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_xnod, size_t internal=0)
space_constant::valued_type valued_tag() const