Rheolef  7.1
an efficient C++ finite element environment
basis_fem_sides.cc
Go to the documentation of this file.
1 // There are four possibilities :
22 // 1) space Mh (omega["sides"], "P1d");
23 // => constructed side by side
24 // with local discontinuous connections between sides of an element
25 // along vertices (2d) or edges (3d)
26 // and global continuity for uni-valued approx along side
27 // status : supported in Rheolef
28 // used by HDG for "lambda"
29 // 2) space Mh (omega["sides"], "P1");
30 // => constructed side by side
31 // with local continuous connections between sides of an element
32 // and global continuity for uni-valued approx along side
33 // status : *NOT* supported in Rheolef
34 // 3) space Mhs (omega, "P1d[sides]");
35 // => constructed element by element, defined on each element boundary
36 // with local discontinuous connections between sides of an element
37 // along vertices (2d) or edges (3d)
38 // and global continuity for uni-valued approx along side
39 // status : supported in Rheolef
40 // used by HDG-POST for "lambda_s"
41 // 4) space Mhs (omega, "P1[sides,global_discontinuous]");
42 // => constructed element by element, defined on each element boundary
43 // with local continuous connections between sides of an element
44 // along vertices (2d) or edges (3d)
45 // and global discontinuity for uni-valued approx along side
46 // status : *NOT* supported in Rheolef
47 //
48 // The case 3 is the aim of this file: "Pkd[sides]"
49 //
50 #include "basis_fem_sides.h"
51 #include "rheolef/rheostream.h"
52 #include "rheolef/iorheo.h"
53 
54 #include "equispaced.icc"
55 #include "warburton.icc"
56 #include "fekete.icc"
57 #include "eigen_util.h"
58 
59 namespace rheolef {
60 
61 using namespace std;
62 
63 // =========================================================================
64 // basis members
65 // =========================================================================
66 template<class T>
68 {
69 }
70 template<class T>
72  : basis_rep<T>(side_basis.option()),
73  _side_basis(side_basis),
74  _hat_node(),
75  _first_sid_inod(),
76  _first_sid_idof(),
77  _sid_value()
78 {
82  base::_piola_fem = _side_basis.get_piola_fem();
83  if (base::option().is_continuous()) {
84  // "Pkd[sices]" : discontinuous and multi-valued along sides connections (edges in 3d, vertives in 2d)
85  // "Pk[sides]" : could be continuous, but not yet supported
86  warning_macro("continuous approximation restricted to sides not yet supported");
87  error_macro ("HINT: replace \""<<base::name()<<"\" by \"P"<<degree()<<"d[sides]\"");
88  }
89 }
90 template<class T>
91 void
93 {
94  size_type k = degree();
95  for (size_type d = 0; d < 4; ++d) {
96  base::_ndof_on_subgeo [d].fill (0);
97  base::_nnod_on_subgeo [d].fill (0);
98  if (d == 0) continue;
101  subgeo_variant++) {
102  // customized for _side_basis="Pkd" on subgeo-variant : all dofs in "_side_basis" are merged on (d-1) = dim(subgeo_variant)
103  // => all dofs are associated to the "d" dimension in the present element
104  // => cannot handle yet the case _side_basis="RTkd" here : it requires "_side_basis.ndof_on_subgeo_continuous(d,subgeo_variant)
105  base::_ndof_on_subgeo [d][subgeo_variant] = _side_basis.ndof_on_subgeo (d-1, subgeo_variant);
106  base::_nnod_on_subgeo [d][subgeo_variant] = _side_basis.nnod_on_subgeo (d-1, subgeo_variant);
107  }
108  }
109  bool is_continuous = false;
110  base::_helper_discontinuous_ndof_on_subgeo_inplace_change (is_continuous, base::_ndof_on_subgeo);
111  base::_helper_discontinuous_ndof_on_subgeo_inplace_change (is_continuous, base::_nnod_on_subgeo);
112  base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (base::_ndof_on_subgeo, base::_first_idof_by_dimension);
113  base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (base::_nnod_on_subgeo, base::_first_inod_by_dimension);
114 }
115 template<class T>
116 void
118 {
119  size_type k = degree();
120  size_type d = tilde_K.dimension();
121  size_type variant = tilde_K.variant();
122 
123  // nodes:
124  _first_sid_idof [variant].fill (0);
125  size_type loc_nnod = base::_first_inod_by_dimension [variant][d+1];
126  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& tilde_node = _hat_node [variant];
127  tilde_node.resize (loc_nnod);
128  size_type loc_first_sid_inod = 0;
129  for (size_type loc_isid = 0, loc_nsid = tilde_K.n_side(); loc_isid < loc_nsid; ++loc_isid) {
130  side_information_type sid (tilde_K, loc_isid);
131  reference_element hat_S = sid.hat;
132  const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& side_hat_node = _side_basis.hat_node (hat_S);
133  size_type loc_sid_nnod = side_hat_node.size();
134  for (size_type loc_sid_inod = 0; loc_sid_inod < loc_sid_nnod; ++loc_sid_inod) {
135  // piola transform side nodes from hat_S to partial tilde_K
136  size_type loc_inod = loc_first_sid_inod + loc_sid_inod;
137  tilde_node[loc_inod] = reference_element_face_transformation (tilde_K, sid, side_hat_node[loc_sid_inod]);
138  }
139  loc_first_sid_inod += loc_sid_nnod;
140  _first_sid_idof [variant][loc_isid+1] = loc_first_sid_inod;
141  }
142  _first_sid_inod [variant] = _first_sid_idof [variant];
143 }
144 template<class T>
145 const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&
147 {
148  base::_initialize_data_guard (hat_K);
149  return _hat_node [hat_K.variant()];
150 }
151 // evaluation on a side of all related basis functions at hat_x in hat_S
152 // NOTE: by convention, the value array has a full size for all sides dofs
153 // and is filled by zero on all others sides
154 template<class T>
155 void
157  reference_element tilde_K,
158  const side_information_type& sid,
159  const point_basic<T>& hat_x,
160  Eigen::Matrix<T,Eigen::Dynamic,1>& value) const
161 {
162  base::_initialize_data_guard (tilde_K);
163  size_type variant = tilde_K.variant();
164  size_type loc_ndof = base::ndof (tilde_K);
165  value.resize (loc_ndof);
166  value.fill(T(0));
167  size_type loc_first_sid_idof = _first_sid_inod [variant][sid.loc_isid];
168  reference_element hat_S = sid.hat;
169  _side_basis.evaluate (hat_S, hat_x, _sid_value);
170  for (size_type loc_sid_idof = 0, loc_sid_ndof = _sid_value.size(); loc_sid_idof < loc_sid_ndof; ++loc_sid_idof) {
171  size_type loc_idof = loc_first_sid_idof + loc_sid_idof;
172  value [loc_idof] = _sid_value [loc_sid_idof];
173  }
174 }
175 // extract local dof-indexes on a side
176 template<class T>
179  reference_element tilde_K,
180  const side_information_type& sid) const
181 {
182  base::_initialize_data_guard (tilde_K);
183  return _first_sid_inod [tilde_K.variant()][sid.loc_isid+1]
184  - _first_sid_inod [tilde_K.variant()][sid.loc_isid];
185 }
186 template<class T>
187 void
189  reference_element tilde_K,
190  const side_information_type& sid,
191  Eigen::Matrix<size_type,Eigen::Dynamic,1>& loc_idof) const
192 {
193  base::_initialize_data_guard (tilde_K);
194  size_type variant = tilde_K.variant();
195  size_type first_loc_inod = _first_sid_inod [variant][sid.loc_isid];
196  size_type loc_sid_nnod = _first_sid_inod [variant][sid.loc_isid+1]
197  - _first_sid_inod [variant][sid.loc_isid];
198  loc_idof.resize (loc_sid_nnod);
199  for (size_type loc_sid_inod = 0; loc_sid_inod < loc_sid_nnod; ++loc_sid_inod) {
200  loc_idof [loc_sid_inod] = first_loc_inod + loc_sid_inod;
201  }
202 }
203 // dofs for a scalar-valued function
204 template<class T>
205 void
207  reference_element tilde_K,
208  const Eigen::Matrix<T,Eigen::Dynamic,1>& f_xnod,
209  Eigen::Matrix<T,Eigen::Dynamic,1>& dof) const
210 {
211  check_macro (is_nodal(), "_compute_dofs: only supported for nodal basis");
212  dof = f_xnod;
213 }
214 // -----------------------------------------------------------------------------
215 // visu: gnuplot in elevation for 2d (t,q)
216 // -----------------------------------------------------------------------------
217 template<class T>
218 void
220 {
221  bool verbose = iorheo::getverbose(os);
222  bool clean = iorheo::getclean(os);
223  bool execute = iorheo::getexecute(os);
224  size_type nsub = iorheo::getsubdivide(os);
225  if (nsub <= 1) nsub = 20; // default value
226 
227  size_type loc_ndof = base::ndof (tilde_K);
228  string basename = "basis-" + base::name() + "-" + tilde_K.name();
229  string filelist;
230  // --------------------
231  // .plot
232  // --------------------
233  string plot_name = basename + ".plot";
234  string gdat_name = basename + ".gdat";
235  ofstream plot (plot_name.c_str());
236  cerr << "! file \"" << plot_name << "\" created" << endl;
237  filelist += " " + plot_name;
238  size_type d = tilde_K.dimension();
239  check_macro (d==2, "unsupported dimension " << d);
240  plot << "gdat = \"" << gdat_name << "\"" << endl
241  << "set colors classic" << endl
242  << "set size square" << endl
243  << "pause_duration = 1.5" << endl
244  << "set view 55, 145" << endl
245  ;
246  if (tilde_K.variant() == reference_element::t) {
247  plot << "set xrange [-0.1:1.1]" << endl
248  << "set yrange [-0.1:1.1]" << endl
249  << "set zrange [-1.1:1.1]" << endl
250  << "set arrow from 0,0,0 to 1,0,0 nohead lc 0" << endl
251  << "set arrow from 1,0,0 to 0,1,0 nohead lc 0" << endl
252  << "set arrow from 0,1,0 to 0,0,0 nohead lc 0" << endl
253  ;
254  } else {
255  plot << "set xrange [-1.1:1.1]" << endl
256  << "set yrange [-1.1:1.1]" << endl
257  << "set zrange [-1.1:1.1]" << endl
258  << "set arrow from -1,-1,0 to 1,-1,0 nohead lc 0" << endl
259  << "set arrow from 1,-1,0 to 1, 1,0 nohead lc 0" << endl
260  << "set arrow from 1, 1,0 to -1, 1,0 nohead lc 0" << endl
261  << "set arrow from -1, 1,0 to -1,-1,0 nohead lc 0" << endl
262  ;
263  }
264  plot << "do for [loc_isid=0:"<<tilde_K.n_side()-1<<"] {" << endl
265  << " loc_sid_ndof=" << degree()+1 << endl
266  << " do for [i=0:loc_sid_ndof-1] {" << endl
267  << " ip3=i+3" << endl
268  << " splot gdat index loc_isid u 1:2:ip3 t sprintf(\"side %d: L%d\",loc_isid,i) w l lc 1" << endl
269  << " pause pause_duration" << endl
270  << " }" << endl
271  << "}" << endl
272  << "pause -1 \"<return>\"" << endl
273  ;
274  plot.close();
275  // --------------------
276  // .gdat
277  // --------------------
278  ofstream gdat (gdat_name.c_str());
279  cerr << "! file \"" << gdat_name << "\" created" << endl;
280  filelist += " " + gdat_name;
281  gdat << setprecision(std::numeric_limits<T>::digits10)
282  << "# basis " << base::name() << endl
283  << "# element " << tilde_K.name() << endl
284  << "# degree " << degree() << endl
285  ;
286  Eigen::Matrix<T,Eigen::Dynamic,1> value;
287  for (size_type loc_isid = 0, loc_nsid = tilde_K.n_side(); loc_isid < loc_nsid; ++loc_isid) {
288  side_information_type sid (tilde_K, loc_isid);
289  reference_element hat_S = sid.hat;
290  size_type loc_sid_ndof = first_sid_idof (tilde_K,loc_isid+1) - first_sid_idof (tilde_K,loc_isid);
291  gdat << "# side " << loc_isid << ": size " << loc_sid_ndof << endl
292  << "# x y";
293  for (size_type loc_sid_idof = 0; loc_sid_idof < loc_sid_ndof; ++loc_sid_idof) {
294  gdat << " L" << loc_sid_idof+1;
295  }
296  gdat << endl;
297  switch (hat_S.variant()) {
298  case reference_element::e: {
299  for (size_type i = 0; i <= nsub; i++) {
300  point_basic<T> sid_hat_x (T(int(i))/T(int(nsub)));
301  evaluate_on_side (tilde_K, sid, sid_hat_x, value);
302  point_basic<T> tilde_x = reference_element_face_transformation (tilde_K, sid, sid_hat_x);
303  gdat << tilde_x[0] << " " << tilde_x[1];
304  for (size_type loc_idof = first_sid_idof (tilde_K,loc_isid), last_loc_idof = first_sid_idof (tilde_K,loc_isid+1);
305  loc_idof < last_loc_idof; ++loc_idof) {
306  gdat << " " << value[loc_idof];
307  }
308  gdat << endl;
309  }
310  gdat << endl << endl;
311  break;
312  }
313  default:
314  error_macro ("unexpected side type `"<<hat_S.name()<<"'");
315  }
316  }
317  // -----------
318  if (execute) {
319  // -----------
320  string command = "gnuplot " + plot_name;
321  if (verbose) clog << "! " << command << endl;
322  int status = system (command.c_str());
323  }
324  // -----------
325  if (clean) {
326  // -----------
327  string command = "/bin/rm -f " + filelist;
328  if (verbose) clog << "! " << command << endl;
329  int status = system (command.c_str());
330  }
331 #ifdef TODO
332 #endif // TODO
333 }
334 // ----------------------------------------------------------------------------
335 // instanciation in library
336 // ----------------------------------------------------------------------------
337 #define _RHEOLEF_instanciation(T) \
338 template class basis_fem_sides<T>;
339 
341 
342 }// namespace rheolef
rheolef::space_numbering::ndof
size_type ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
Definition: space_numbering.cc:28
rheolef::reference_element::last_variant_by_dimension
static variant_type last_variant_by_dimension(size_type dim)
Definition: reference_element.h:150
rheolef::reference_element::n_side
size_type n_side() const
Definition: reference_element.h:104
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
warburton.icc
rheolef::basis_fem_sides::family_index
size_type family_index() const
Definition: basis_fem_sides.h:65
rheolef::reference_element::size
size_type size() const
Definition: reference_element.h:102
eigen_util.h
warning_macro
#define warning_macro(message)
Definition: dis_macros.h:53
rheolef::basis_fem_sides::_initialize_cstor_sizes
void _initialize_cstor_sizes() const
Definition: basis_fem_sides.cc:92
rheolef::point_basic
Definition: point.h:87
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::basis_fem_sides::degree
size_type degree() const
Definition: basis_fem_sides.h:66
rheolef::basis_fem_sides::put_scalar_valued
virtual void put_scalar_valued(std::ostream &os, reference_element hat_K) const
Definition: basis_fem_sides.cc:219
rheolef::reference_element_face_transformation
point_basic< T > reference_element_face_transformation(reference_element tilde_K, const side_information_type &sid, const point_basic< T > &sid_hat_x)
Definition: reference_element_face_transformation.cc:33
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::value
rheolef::std value
rheolef::basis_fem_sides::_compute_dofs
void _compute_dofs(reference_element hat_K, const Eigen::Matrix< T, Eigen::Dynamic, 1 > &f_xnod, Eigen::Matrix< T, Eigen::Dynamic, 1 > &dof) const
Definition: basis_fem_sides.cc:206
rheolef::basis_option::set_restricted_to_sides
void set_restricted_to_sides(bool r=true)
Definition: basis_option.h:281
rheolef::basis_fem_sides::hat_node
const Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > & hat_node(reference_element hat_K) const
Definition: basis_fem_sides.cc:146
rheolef::basis_rep::size_type
reference_element::size_type size_type
Definition: basis.h:214
rheolef::reference_element::first_variant_by_dimension
static variant_type first_variant_by_dimension(size_type dim)
Definition: reference_element.h:148
rheolef::basis_fem_sides::family_name
std::string family_name() const
Definition: basis_fem_sides.h:64
fekete.icc
rheolef::evaluate_on_side
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
Definition: field_expr_recursive.h:894
rheolef::side_information_type
Definition: reference_element_face_transformation.h:37
rheolef::basis_basic
Definition: basis.h:206
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
mkgeo_ball.variant
variant
Definition: mkgeo_ball.sh:149
rheolef::basis_rep::_piola_fem
piola_fem< T > _piola_fem
Definition: basis.h:394
rheolef::reference_element::name
char name() const
Definition: reference_element.h:100
basis_fem_sides.h
mkgeo_ball.clean
clean
Definition: mkgeo_ball.sh:335
rheolef::reference_element::variant
variant_type variant() const
Definition: reference_element.h:99
rheolef::basis_rep::_sopt
basis_option _sopt
Definition: basis.h:393
rheolef::side_information_type::loc_isid
size_t loc_isid
Definition: reference_element_face_transformation.h:38
rheolef::basis_fem_sides::evaluate_on_side
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
Definition: basis_fem_sides.cc:156
mkgeo_couette.basename
basename
Definition: mkgeo_couette.sh:73
rheolef::basis_rep::_name
std::string _name
Definition: basis.h:392
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::basis_fem_sides::~basis_fem_sides
~basis_fem_sides()
Definition: basis_fem_sides.cc:67
Float
see the Float page for the full documentation
rheolef::basis_fem_sides::_initialize_data
void _initialize_data(reference_element hat_K) const
Definition: basis_fem_sides.cc:117
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
rheolef::basis_fem_sides::size_type
reference_element::size_type size_type
Definition: basis_fem_sides.h:54
rheolef::reference_element::dimension
size_type dimension() const
Definition: reference_element.h:101
mkgeo_ball.verbose
verbose
Definition: mkgeo_ball.sh:133
equispaced.icc
rheolef::basis_rep::standard_naming
static std::string standard_naming(std::string family_name, size_t degree, const basis_option &sopt)
Definition: basis_rep.cc:44
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
rheolef::basis_fem_sides::_side_basis
basis_basic< T > _side_basis
Definition: basis_fem_sides.h:108
rheolef::basis_fem_sides::local_idof_on_side
void local_idof_on_side(reference_element hat_K, const side_information_type &sid, Eigen::Matrix< size_type, Eigen::Dynamic, 1 > &loc_idof) const
Definition: basis_fem_sides.cc:188
rheolef::basis_rep::option
const basis_option & option() const
Definition: basis.h:238
rheolef::side_information_type::hat
reference_element hat
Definition: reference_element_face_transformation.h:43
rheolef::basis_fem_sides::local_ndof_on_side
size_type local_ndof_on_side(reference_element hat_K, const side_information_type &sid) const
Definition: basis_fem_sides.cc:178
mkgeo_contraction.status
status
Definition: mkgeo_contraction.sh:290
mkgeo_ball.command
command
Definition: mkgeo_ball.sh:136
rheolef::basis_rep::name
std::string name() const
Definition: basis.h:228
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133
rheolef::basis_rep::is_continuous
bool is_continuous() const
Definition: basis.h:240
rheolef::basis_rep
Definition: basis.h:209
T
Expr1::float_type T
Definition: field_expr.h:218
rheolef::basis_fem_sides::basis_fem_sides
basis_fem_sides(const basis_basic< T > &side_basis)
Definition: basis_fem_sides.cc:71