24#include "rheolef/rheostream.h"
25#include "rheolef/iorheo.h"
46 _vec_basis (vec_basis),
55 <<
"\": a vector-valued basis is expected for the trace_n(.) operator");
56 base::_sopt.set_trace_n (
true);
67 warning_macro(
"normal trace of continuous approximation not yet supported");
77 base::_ndof_on_subgeo_internal [
d].fill (0);
78 base::_nnod_on_subgeo_internal [
d].fill (0);
86 base::_ndof_on_subgeo_internal [
d][subgeo_variant] = _tr_n_basis.ndof_on_subgeo (
d-1, subgeo_variant);
87 base::_nnod_on_subgeo_internal [
d][subgeo_variant] = _tr_n_basis.nnod_on_subgeo (
d-1, subgeo_variant);
90 bool is_continuous =
false;
91 base::_helper_make_discontinuous_ndof_on_subgeo (is_continuous, base::_ndof_on_subgeo_internal, base::_ndof_on_subgeo);
92 base::_helper_make_discontinuous_ndof_on_subgeo (is_continuous, base::_nnod_on_subgeo_internal, base::_nnod_on_subgeo);
93 base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (base::_ndof_on_subgeo_internal, base::_first_idof_by_dimension_internal);
94 base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (base::_ndof_on_subgeo, base::_first_idof_by_dimension);
95 base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (base::_nnod_on_subgeo_internal, base::_first_inod_by_dimension_internal);
96 base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (base::_nnod_on_subgeo, base::_first_inod_by_dimension);
107 _first_sid_idof [variant].fill (0);
108 size_type loc_nnod = base::_first_inod_by_dimension [variant][
d+1];
109 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& tilde_node = _hat_node [variant];
110 tilde_node.resize (loc_nnod);
112 for (
size_type loc_isid = 0, loc_nsid = tilde_K.
n_side(); loc_isid < loc_nsid; ++loc_isid) {
115 const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& side_hat_node = _tr_n_basis.hat_node (hat_S);
117 for (
size_type loc_sid_inod = 0; loc_sid_inod < loc_sid_nnod; ++loc_sid_inod) {
119 size_type loc_inod = loc_first_sid_inod + loc_sid_inod;
122 loc_first_sid_inod += loc_sid_nnod;
123 _first_sid_idof [variant][loc_isid+1] = loc_first_sid_inod;
125 _first_sid_inod [variant] = _first_sid_idof [variant];
128const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&
131 base::_initialize_data_guard (hat_K);
132 return _hat_node [hat_K.
variant()];
143 Eigen::Matrix<T,Eigen::Dynamic,1>& value)
const
145 base::_initialize_data_guard (tilde_K);
147 size_type loc_ndof = base::ndof (tilde_K);
148 value.resize (loc_ndof);
152 _tr_n_basis.evaluate (hat_S, hat_x, _sid_value);
153 for (
size_type loc_sid_idof = 0, loc_sid_ndof = _sid_value.size(); loc_sid_idof < loc_sid_ndof; ++loc_sid_idof) {
154 size_type loc_idof = loc_first_sid_idof + loc_sid_idof;
155 value [loc_idof] = _sid_value [loc_sid_idof];
165 base::_initialize_data_guard (tilde_K);
174 Eigen::Matrix<size_type,Eigen::Dynamic,1>& loc_idof)
const
176 base::_initialize_data_guard (tilde_K);
180 - _first_sid_inod [variant][sid.
loc_isid];
181 loc_idof.resize (loc_sid_nnod);
182 for (
size_type loc_sid_inod = 0; loc_sid_inod < loc_sid_nnod; ++loc_sid_inod) {
183 loc_idof [loc_sid_inod] = first_loc_inod + loc_sid_inod;
191 const Eigen::Matrix<T,Eigen::Dynamic,1>& f_xnod,
192 Eigen::Matrix<T,Eigen::Dynamic,1>& dof)
const
194 check_macro (is_nodal(),
"_compute_dofs: only supported for nodal basis");
204 bool verbose = iorheo::getverbose(os);
205 bool clean = iorheo::getclean(os);
206 bool execute = iorheo::getexecute(os);
207 size_type nsub = iorheo::getsubdivide(os);
208 if (nsub <= 1) nsub = 20;
210 size_type loc_ndof = base::ndof (tilde_K);
211 string basename =
"basis-" + base::name() +
"-" + tilde_K.
name();
216 string plot_name = basename +
".plot";
217 string gdat_name = basename +
".gdat";
218 ofstream plot (plot_name.c_str());
219 cerr <<
"! file \"" << plot_name <<
"\" created" << endl;
220 filelist +=
" \"" + plot_name +
"\"";
223 plot <<
"gdat = \"" << gdat_name <<
"\"" << endl
224 <<
"set colors classic" << endl
225 <<
"set size square" << endl
226 <<
"pause_duration = 1.5" << endl
227 <<
"set view 55, 145" << endl
230 plot <<
"set xrange [-0.1:1.1]" << endl
231 <<
"set yrange [-0.1:1.1]" << endl
232 <<
"set zrange [-1.1:1.1]" << endl
233 <<
"set arrow from 0,0,0 to 1,0,0 nohead lc 0" << endl
234 <<
"set arrow from 1,0,0 to 0,1,0 nohead lc 0" << endl
235 <<
"set arrow from 0,1,0 to 0,0,0 nohead lc 0" << endl
238 plot <<
"set xrange [-1.1:1.1]" << endl
239 <<
"set yrange [-1.1:1.1]" << endl
240 <<
"set zrange [-1.1:1.1]" << endl
241 <<
"set arrow from -1,-1,0 to 1,-1,0 nohead lc 0" << endl
242 <<
"set arrow from 1,-1,0 to 1, 1,0 nohead lc 0" << endl
243 <<
"set arrow from 1, 1,0 to -1, 1,0 nohead lc 0" << endl
244 <<
"set arrow from -1, 1,0 to -1,-1,0 nohead lc 0" << endl
247 plot <<
"do for [loc_isid=0:"<<tilde_K.
n_side()-1<<
"] {" << endl
248 <<
" loc_sid_ndof=" << degree()+1 << endl
249 <<
" do for [i=0:loc_sid_ndof-1] {" << endl
250 <<
" ip3=i+3" << endl
251 <<
" splot gdat index loc_isid u 1:2:ip3 t sprintf(\"side %d: L%d\",loc_isid,i) w l lc 1" << endl
252 <<
" pause pause_duration" << endl
255 <<
"pause -1 \"<return>\"" << endl
261 ofstream gdat (gdat_name.c_str());
262 cerr <<
"! file \"" << gdat_name <<
"\" created" << endl;
263 filelist +=
" \"" + gdat_name +
"\"";
264 gdat << setprecision(std::numeric_limits<T>::digits10)
265 <<
"# basis " << base::name() << endl
266 <<
"# element " << tilde_K.
name() << endl
267 <<
"# degree " << degree() << endl
269 Eigen::Matrix<T,Eigen::Dynamic,1> value;
270 for (
size_type loc_isid = 0, loc_nsid = tilde_K.
n_side(); loc_isid < loc_nsid; ++loc_isid) {
273 size_type loc_sid_ndof = first_sid_idof (tilde_K,loc_isid+1) - first_sid_idof (tilde_K,loc_isid);
274 gdat <<
"# side " << loc_isid <<
": size " << loc_sid_ndof << endl
276 for (
size_type loc_sid_idof = 0; loc_sid_idof < loc_sid_ndof; ++loc_sid_idof) {
277 gdat <<
" L" << loc_sid_idof+1;
286 gdat << tilde_x[0] <<
" " << tilde_x[1];
287 for (
size_type loc_idof = first_sid_idof (tilde_K,loc_isid), last_loc_idof = first_sid_idof (tilde_K,loc_isid+1);
288 loc_idof < last_loc_idof; ++loc_idof) {
289 gdat <<
" " << value[loc_idof];
293 gdat << endl << endl;
303 string command =
"gnuplot \"" + plot_name +
"\"";
304 if (verbose) clog <<
"! " << command << endl;
305 int status = system (command.c_str());
310 string command =
"/bin/rm -f " + filelist;
311 if (verbose) clog <<
"! " << command << endl;
312 int status = system (command.c_str());
318#define _RHEOLEF_instanciation(T) \
319template class basis_fem_trace_n<T>;
#define _RHEOLEF_instanciation(T, M, A)
see the Float page for the full documentation
size_type local_ndof_on_side(reference_element hat_K, const side_information_type &sid) const
void evaluate_on_side(reference_element hat_K, const side_information_type &sid, const point_basic< T > &hat_x, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value) const
reference_element::size_type size_type
basis_basic< T > _tr_n_basis
void local_idof_on_side(reference_element hat_K, const side_information_type &sid, Eigen::Matrix< size_type, Eigen::Dynamic, 1 > &loc_idof) const
void _initialize_cstor_sizes() const
virtual void put_scalar_valued(std::ostream &os, reference_element hat_K) const
void _compute_dofs(reference_element hat_K, const Eigen::Matrix< T, Eigen::Dynamic, 1 > &f_xnod, Eigen::Matrix< T, Eigen::Dynamic, 1 > &dof) const
size_type family_index() const
const Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > & hat_node(reference_element hat_K) const
std::string family_name() const
basis_fem_trace_n(const basis_basic< T > &vec_basis)
basis_basic< T > _vec_basis
void _initialize_data(reference_element hat_K) const
see the basis_option page for the full documentation
bool is_continuous() const
see the reference_element page for the full documentation
static const variant_type e
static variant_type last_variant_by_dimension(size_type dim)
size_type dimension() const
static variant_type first_variant_by_dimension(size_type dim)
variant_type variant() const
static const variant_type t
#define error_macro(message)
#define warning_macro(message)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
point_basic< T > reference_element_face_transformation(reference_element tilde_K, const side_information_type &sid, const point_basic< T > &sid_hat_x)
void evaluate_on_side(const geo_basic< float_type, M > &omega_K, const geo_element &K, const side_information_type &sid, Eigen::Matrix< Result, Eigen::Dynamic, 1 > &value) const