Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
34namespace 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
42template<class T>
43void
46 size_t order_in,
47 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& hat_xnod,
48 size_t internal = 0)
49{
50trace_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]);
56trace_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;
88trace_macro("hat_K="<<hat_K.name()<<": order="<<order<<", first="<<first<<", last="<<last<<", size="<<size);
89 hat_xnod.resize (size);
90 switch (hat_K.variant()) {
92 hat_xnod.resize (1);
93 hat_xnod [0] = point_basic<T> (T(0));
94 break;
95 }
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)));
100trace_macro("e::lattice("<<i<<"): hat_xnod ["<<loc_idof<<"]="<<hat_xnod [loc_idof]);
101 }
102 break;
103 }
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)));
109trace_macro("t::lattice("<<i<<","<<j<<"): hat_xnod ["<<loc_idof<<"]="<<hat_xnod [loc_idof]);
110 }}
111 break;
112 }
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)));
118trace_macro("q::lattice("<<i<<","<<j<<"): hat_xnod ["<<loc_idof<<"]="<<hat_xnod [loc_idof]);
119 }}
120 break;
121 }
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 }
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 }
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:430
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
#define trace_macro(message)
Definition dis_macros.h:111
#define error_macro(message)
Definition dis_macros.h:49
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.
void reference_element_barycenter(reference_element hat_K, point_basic< T > &c)
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)