Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
basis_fem_trace_n.cc
Go to the documentation of this file.
1
21// 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
32namespace rheolef {
33
34using namespace std;
35
36// =========================================================================
37// basis members
38// =========================================================================
39template<class T>
43template<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
59 std::string tr_name = base::standard_naming ("P", tr_n_degree, tr_n_opt);
61 base::_name = base::standard_naming (family_name(), family_index(), base::_sopt);
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}
71template<class T>
72void
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}
98template<class T>
99void
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}
127template<class T>
128const 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
137template<class T>
138void
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
159template<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}
169template<class T>
170void
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
187template<class T>
188void
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// -----------------------------------------------------------------------------
200template<class T>
201void
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()) {
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) \
319template class basis_fem_trace_n<T>;
320
322
323}// namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
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
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
bool is_continuous() const
Definition basis.h:240
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
#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)")
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
STL namespace.