Rheolef  7.1
an efficient C++ finite element environment
equispaced.icc
Go to the documentation of this file.
1 #ifndef _RHEOLEF_EQUISPACED_ICC
2 #define _RHEOLEF_EQUISPACED_ICC
23 //
24 // compute the equispaced point set on the reference element
25 //
26 // author: Pierre.Saramito@imag.fr
27 //
28 // date: 2 september 2017
29 //
30 #include "rheolef/reference_element.h"
32 #include "rheolef/compiler_eigen.h"
33 
34 namespace rheolef {
35 
36 // The optional argument "interior" specifies
37 // how many points from the boundary to omit.
38 // For example, on an "edge" reference element with order=2 and interior=0,
39 // this function will return the vertices and midpoint,
40 // but with interior=1, it will only return the midpoint.
41 
42 template<class T>
43 void
45  reference_element hat_K,
46  size_t order_in,
47  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& hat_xnod,
48  size_t internal = 0)
49 {
50 trace_macro("hat_K="<<hat_K.name()<<": order_in="<<order_in<<", internal="<<internal);
52  typedef point_basic<size_type> ilat;
53  if (order_in == 0) {
54  hat_xnod.resize(1);
55  reference_element_barycenter (hat_K, hat_xnod[0]);
56 trace_macro("barycenter: hat_xnod [0]="<<hat_xnod[0]);
57  return;
58  }
59  size_t d = hat_K.dimension();
60  // when only internal nodes: compute the effective "order" and the lattice size
61  size_t order = 0, size = 0;
62  switch (hat_K.variant()) {
64  order = order_in;
65  size = reference_element::n_node(hat_K.variant(), order);
66  break;
70  order = order_in + (d+1)*internal;
71  size = (order >= (d+1)*internal) ? reference_element::n_node(hat_K.variant(), order-(d+1)*internal) : 0;
72  break;
75  order = order_in + 2*internal;
76  size = (order >= 2*internal) ? reference_element::n_node(hat_K.variant(), order-2*internal) : 0;
77  break;
79  check_macro (internal == 0, "internal: not yet");
80  order = order_in;
81  size = reference_element::n_node(hat_K.variant(), order);
82  break;
83  default:
84  error_macro ("unexpected element type `"<<hat_K.name()<<"'");
85  }
86  size_t first = internal;
87  size_t last = (order >= internal) ? order-internal : 0;
88 trace_macro("hat_K="<<hat_K.name()<<": order="<<order<<", first="<<first<<", last="<<last<<", size="<<size);
89  hat_xnod.resize (size);
90  switch (hat_K.variant()) {
91  case reference_element::p: {
92  hat_xnod.resize (1);
93  hat_xnod [0] = point_basic<T> (T(0));
94  break;
95  }
96  case reference_element::e: {
97  for (size_type i = first; i <= last; i++) {
98  size_type loc_idof = reference_element_e::ilat2loc_inod (order-(d+1)*internal, ilat(i-internal));
99  hat_xnod [loc_idof] = point_basic<T> ((T(int(i)))/T(int(order)));
100 trace_macro("e::lattice("<<i<<"): hat_xnod ["<<loc_idof<<"]="<<hat_xnod [loc_idof]);
101  }
102  break;
103  }
104  case reference_element::t: {
105  for (size_type j = first; j <= last; j++) {
106  for (size_type i = first; i+j <= last; i++) {
107  size_type loc_idof = reference_element_t::ilat2loc_inod (order-(d+1)*internal, ilat(i-internal,j-internal));
108  hat_xnod [loc_idof] = point_basic<T> ((T(int(i)))/T(int(order)), T(int(j))/T(int(order)));
109 trace_macro("t::lattice("<<i<<","<<j<<"): hat_xnod ["<<loc_idof<<"]="<<hat_xnod [loc_idof]);
110  }}
111  break;
112  }
113  case reference_element::q: {
114  for (size_type j = first; j <= last; j++) {
115  for (size_type i = first; i <= last; i++) {
116  size_type loc_idof = reference_element_q::ilat2loc_inod (order-2*internal, ilat(i-internal,j-internal));
117  hat_xnod [loc_idof] = point_basic<T> (-1+2*(T(int(i)))/T(int(order)), -1+2*T(int(j))/T(int(order)));
118 trace_macro("q::lattice("<<i<<","<<j<<"): hat_xnod ["<<loc_idof<<"]="<<hat_xnod [loc_idof]);
119  }}
120  break;
121  }
122  case reference_element::T: {
123  check_macro (internal == 0, "internal: not yet");
124  for (size_type k = first; k <= last; k++) {
125  for (size_type j = first; j+k <= last; j++) {
126  for (size_type i = first; i+j+k <= last; i++) {
127  size_type loc_idof = reference_element_T::ilat2loc_inod (order, ilat(i,j,k));
128  hat_xnod [loc_idof] = point_basic<T> ((T(int(i)))/T(int(order)), T(int(j))/T(int(order)), T(int(k))/T(int(order)));
129  }}}
130  break;
131  }
132  case reference_element::P: {
133  check_macro (internal == 0, "internal: not yet");
134  for (size_type k = first; k <= last; k++) {
135  for (size_type j = first; j <= last; j++) {
136  for (size_type i = first; i+j <= last; i++) {
137  size_type loc_idof = reference_element_P::ilat2loc_inod (order, ilat(i,j,k));
138  hat_xnod [loc_idof] = point_basic<T> ((T(int(i)))/T(int(order)), T(int(j))/T(int(order)), -1+2*T(int(k))/T(int(order)));
139  }}}
140  break;
141  }
142  case reference_element::H: {
143  check_macro (internal == 0, "internal: not yet");
144  for (size_type k = first; k <= last; k++) {
145  for (size_type j = first; j <= last; j++) {
146  for (size_type i = first; i <= last; i++) {
147  size_type loc_idof = reference_element_H::ilat2loc_inod (order, ilat(i,j,k));
148  hat_xnod [loc_idof] = point_basic<T> (-1+2*(T(int(i)))/T(int(order)), -1+2*T(int(j))/T(int(order)), -1+2*T(int(k))/T(int(order)));
149  }}}
150  break;
151  }
152  default: error_macro ("unexpected element type `"<<hat_K.name()<<"'");
153  }
154 }
155 
156 } // namespace rheolef
157 #endif // _RHEOLEF_EQUISPACED_ICC
field::size_type size_type
Definition: branch.cc:425
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
see the reference_element page for the full documentation
static const variant_type H
static const variant_type q
static const variant_type e
static const variant_type p
variant_type variant() const
std::vector< int >::size_type size_type
static size_type n_node(variant_type variant, size_type order)
static const variant_type T
static const variant_type P
static const variant_type t
size_t size_type
Definition: basis_get.cc:76
point_basic< T >
Definition: piola_fem.h:135
Expr1::float_type T
Definition: field_expr.h:261
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.
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)
Definition: equispaced.icc:44
void reference_element_barycenter(reference_element hat_K, point_basic< T > &c)