Rheolef  7.2
an efficient C++ finite element environment
basis_fem_trace_n.cc
Go to the documentation of this file.
1 // See trace-n-fig.fig
22 //
23 #include "basis_fem_trace_n.h"
24 #include "rheolef/rheostream.h"
25 #include "rheolef/iorheo.h"
26 
27 #include "equispaced.icc"
28 #include "warburton.icc"
29 #include "fekete.icc"
30 #include "eigen_util.h"
31 
32 namespace rheolef {
33 
34 using namespace std;
35 
36 // =========================================================================
37 // basis members
38 // =========================================================================
39 template<class T>
41 {
42 }
43 template<class T>
45  : basis_rep<T>(vec_basis.option()),
46  _vec_basis (vec_basis),
47  _tr_n_basis (),
48  _hat_node(),
49  _first_sid_inod(),
50  _first_sid_idof(),
51  _sid_value()
52 {
54  "unexpected "<< _vec_basis.valued() <<"-valued basis \""<<_vec_basis.name()
55  <<"\": a vector-valued basis is expected for the trace_n(.) operator");
56  base::_sopt.set_trace_n (true);
57  size_type tr_n_degree = _vec_basis.family_index(); // RTkd.degree=k+1 but RTkd.family_index=k: it fully contains (Pkd)^d
58  basis_option tr_n_opt = _vec_basis.option();
59  std::string tr_name = base::standard_naming ("P", tr_n_degree, tr_n_opt);
60  _tr_n_basis = basis_basic<T> (tr_name);
62  base::_piola_fem = _tr_n_basis.get_piola_fem();
64  if (base::option().is_continuous()) {
65  // "trace_n(RTkd)" : discontinuous and multi-valued along sides connections (edges in 3d, vertives in 2d)
66  // "trace_n(RTk)" : could be continuous, but not yet supported
67  warning_macro("normal trace of continuous approximation not yet supported");
68  error_macro ("HINT: replace \""<<_vec_basis.name()<<"\" by its dicontinuous version");
69  }
70 }
71 template<class T>
72 void
74 {
75  size_type k = degree();
76  for (size_type d = 0; d < 4; ++d) {
77  base::_ndof_on_subgeo_internal [d].fill (0);
78  base::_nnod_on_subgeo_internal [d].fill (0);
79  if (d == 0) continue;
82  subgeo_variant++) {
83  // customized for _vec_basis="RTkd" on subgeo-variant : all dofs in "_vec_basis" are merged on (d-1) = dim(subgeo_variant)
84  // => all dofs are associated to the "d" dimension in the present element
85  // => cannot handle yet the case _vec_basis="RTkd" here : it requires "_vec_basis.ndof_on_subgeo_internal(d,subgeo_variant)
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);
88  }
89  }
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);
97 }
98 template<class T>
99 void
101 {
102  size_type k = degree();
103  size_type d = tilde_K.dimension();
104  size_type variant = tilde_K.variant();
105 
106  // nodes:
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);
111  size_type loc_first_sid_inod = 0;
112  for (size_type loc_isid = 0, loc_nsid = tilde_K.n_side(); loc_isid < loc_nsid; ++loc_isid) {
113  side_information_type sid (tilde_K, loc_isid);
114  reference_element hat_S = sid.hat;
115  const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& side_hat_node = _tr_n_basis.hat_node (hat_S);
116  size_type loc_sid_nnod = side_hat_node.size();
117  for (size_type loc_sid_inod = 0; loc_sid_inod < loc_sid_nnod; ++loc_sid_inod) {
118  // piola transform side nodes from hat_S to partial tilde_K
119  size_type loc_inod = loc_first_sid_inod + loc_sid_inod;
120  tilde_node[loc_inod] = reference_element_face_transformation (tilde_K, sid, side_hat_node[loc_sid_inod]);
121  }
122  loc_first_sid_inod += loc_sid_nnod;
123  _first_sid_idof [variant][loc_isid+1] = loc_first_sid_inod;
124  }
125  _first_sid_inod [variant] = _first_sid_idof [variant];
126 }
127 template<class T>
128 const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&
130 {
131  base::_initialize_data_guard (hat_K);
132  return _hat_node [hat_K.variant()];
133 }
134 // evaluation on a side of all related basis functions at hat_x in hat_S
135 // NOTE: by convention, the value array has a full size for all sides dofs
136 // and is filled by zero on all others sides
137 template<class T>
138 void
140  reference_element tilde_K,
141  const side_information_type& sid,
142  const point_basic<T>& hat_x,
143  Eigen::Matrix<T,Eigen::Dynamic,1>& value) const
144 {
145  base::_initialize_data_guard (tilde_K);
146  size_type variant = tilde_K.variant();
147  size_type loc_ndof = base::ndof (tilde_K);
148  value.resize (loc_ndof);
149  value.fill(T(0));
150  size_type loc_first_sid_idof = _first_sid_inod [variant][sid.loc_isid];
151  reference_element hat_S = sid.hat;
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];
156  }
157 }
158 // extract local dof-indexes on a side
159 template<class T>
162  reference_element tilde_K,
163  const side_information_type& sid) const
164 {
165  base::_initialize_data_guard (tilde_K);
166  return _first_sid_inod [tilde_K.variant()][sid.loc_isid+1]
167  - _first_sid_inod [tilde_K.variant()][sid.loc_isid];
168 }
169 template<class T>
170 void
172  reference_element tilde_K,
173  const side_information_type& sid,
174  Eigen::Matrix<size_type,Eigen::Dynamic,1>& loc_idof) const
175 {
176  base::_initialize_data_guard (tilde_K);
177  size_type variant = tilde_K.variant();
178  size_type first_loc_inod = _first_sid_inod [variant][sid.loc_isid];
179  size_type loc_sid_nnod = _first_sid_inod [variant][sid.loc_isid+1]
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;
184  }
185 }
186 // dofs for a scalar-valued function
187 template<class T>
188 void
190  reference_element tilde_K,
191  const Eigen::Matrix<T,Eigen::Dynamic,1>& f_xnod,
192  Eigen::Matrix<T,Eigen::Dynamic,1>& dof) const
193 {
194  check_macro (is_nodal(), "_compute_dofs: only supported for nodal basis");
195  dof = f_xnod;
196 }
197 // -----------------------------------------------------------------------------
198 // visu: gnuplot in elevation for 2d (t,q)
199 // -----------------------------------------------------------------------------
200 template<class T>
201 void
203 {
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; // default value
209 
210  size_type loc_ndof = base::ndof (tilde_K);
211  string basename = "basis-" + base::name() + "-" + tilde_K.name();
212  string filelist;
213  // --------------------
214  // .plot
215  // --------------------
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 + "\"";
221  size_type d = tilde_K.dimension();
222  check_macro (d==2, "unsupported dimension " << d);
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
228  ;
229  if (tilde_K.variant() == reference_element::t) {
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
236  ;
237  } else {
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
245  ;
246  }
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
253  << " }" << endl
254  << "}" << endl
255  << "pause -1 \"<return>\"" << endl
256  ;
257  plot.close();
258  // --------------------
259  // .gdat
260  // --------------------
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
268  ;
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) {
271  side_information_type sid (tilde_K, loc_isid);
272  reference_element hat_S = sid.hat;
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
275  << "# x y";
276  for (size_type loc_sid_idof = 0; loc_sid_idof < loc_sid_ndof; ++loc_sid_idof) {
277  gdat << " L" << loc_sid_idof+1;
278  }
279  gdat << endl;
280  switch (hat_S.variant()) {
281  case reference_element::e: {
282  for (size_type i = 0; i <= nsub; i++) {
283  point_basic<T> sid_hat_x (T(int(i))/T(int(nsub)));
284  evaluate_on_side (tilde_K, sid, sid_hat_x, value);
285  point_basic<T> tilde_x = reference_element_face_transformation (tilde_K, sid, sid_hat_x);
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];
290  }
291  gdat << endl;
292  }
293  gdat << endl << endl;
294  break;
295  }
296  default:
297  error_macro ("unexpected side type `"<<hat_S.name()<<"'");
298  }
299  }
300  // -----------
301  if (execute) {
302  // -----------
303  string command = "gnuplot \"" + plot_name + "\"";
304  if (verbose) clog << "! " << command << endl;
305  int status = system (command.c_str());
306  }
307  // -----------
308  if (clean) {
309  // -----------
310  string command = "/bin/rm -f " + filelist;
311  if (verbose) clog << "! " << command << endl;
312  int status = system (command.c_str());
313  }
314 }
315 // ----------------------------------------------------------------------------
316 // instanciation in library
317 // ----------------------------------------------------------------------------
318 #define _RHEOLEF_instanciation(T) \
319 template class basis_fem_trace_n<T>;
320 
322 
323 }// namespace rheolef
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
void local_idof_on_side(reference_element hat_K, const side_information_type &sid, Eigen::Matrix< size_type, Eigen::Dynamic, 1 > &loc_idof) 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)
void _initialize_data(reference_element hat_K) const
see the basis_option page for the full documentation
Definition: basis_option.h:93
void set_trace_n(bool r=true)
Definition: basis_option.h:281
bool is_continuous() const
Definition: basis.h:240
reference_element::size_type size_type
Definition: basis.h:214
piola_fem< T > _piola_fem
Definition: basis.h:424
basis_option _sopt
Definition: basis.h:423
const basis_option & option() const
Definition: basis.h:238
static std::string standard_naming(std::string family_name, size_t degree, const basis_option &sopt)
Definition: basis_rep.cc:44
std::string _name
Definition: basis.h:422
see the reference_element page for the full documentation
static const variant_type e
static variant_type last_variant_by_dimension(size_type dim)
static variant_type first_variant_by_dimension(size_type dim)
variant_type variant() const
static const variant_type t
rheolef::std value
#define error_macro(message)
Definition: dis_macros.h:49
#define warning_macro(message)
Definition: dis_macros.h:53
Expr1::float_type T
Definition: field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
string command
Definition: mkgeo_ball.sh:136
size_type ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
This file is part of Rheolef.
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
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
point_basic< T > reference_element_face_transformation(reference_element tilde_K, const side_information_type &sid, const point_basic< T > &sid_hat_x)