Flow123d  master-f44eb46
fe_p.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file fe_p.cc
15  * @brief
16  */
17 
18 // !! implementation of specializations has to be i *.cc file to avoid multiple definition error during linking
19 #include "fe_p.hh"
20 #include "mesh/ref_element.hh"
21 
22 
23 
24 
25 PolynomialSpace::PolynomialSpace(unsigned int degree, unsigned int dim)
26  : FunctionSpace(dim, 1),
27  degree_(degree)
28 {
29 // computes powers of all monomials up to given @p degree
30 // the order is: 1,x,x^2, y, yx,y^2
31 //
32 // TODO: - check and possibly rewrite to be more clear (use sum_degree temporary
33 // - change order of monomials: 1, x, y, xy, x^2 , y^2 (increasing order)
34 // - allow Q polynomials: 1,x, y, xy
35 // - can use tensor products
36 
37  arma::uvec pows(dim);
38  pows.zeros();
39 
40  unsigned int degree_sum=0;
41  unsigned int i_dim;
42 
43  while (true) {
44  powers.push_back(pows);
45 
46  // increment pows
47  for(i_dim=0; i_dim < dim; i_dim++) {
48  if (degree_sum < degree) {
49  pows[i_dim]++;
50  degree_sum++;
51  break;
52  } else { // if degree_sum == degree, we find first non empty power, free it, and raise the next one
53  degree_sum-=pows[i_dim];
54  pows[i_dim]=0;
55  }
56  }
57  if (i_dim == dim) break; // just after pow == (0, 0, .., degree)
58  }
59 }
60 
61 
62 double PolynomialSpace::basis_value(unsigned int i,
63  const arma::vec &point,
64  unsigned int comp_index
65  ) const
66 {
67  ASSERT_EQ(comp_index, 0);
68  ASSERT_LE(i, powers.size());
69  ASSERT_EQ(point.size(), space_dim_);
70 
71  double v = 1;
72  for (unsigned int j=0; j<this->space_dim_; j++)
73  v *= pow(point[j], (int) powers[i][j]);
74  return v;
75 }
76 
77 
79  const arma::vec &p,
80  unsigned int comp_index
81  ) const
82 {
83  ASSERT_EQ(comp_index, 0);
84  ASSERT_LE(i, powers.size());
85 
86  arma::vec grad(this->space_dim_);
87 
88  for (unsigned int j=0; j<this->space_dim_; j++)
89  {
90  grad[j] = powers[i][j];
91  if (powers[i][j] == 0) continue;
92 
93  for (unsigned int k=0; k<this->space_dim_; k++)
94  {
95  grad[j] *= pow(p[k], (int) (k==j?powers[i][k]-1:powers[i][k]));
96  }
97  }
98  return grad;
99 }
100 
101 
102 
103 
104 
105 
106 
107 
108 
109 
110 
111 template<unsigned int dim>
113 {
114  if (degree_ == 0 || dim == 0)
115  {
116  // we define nodal dof:
117  // coords = barycentric coordinates of the support point,
118  // coefs = 1 (dof value = function value at the point)
119  arma::vec coords = arma::ones<arma::vec>(dim+1)/(dim+1);
120  this->dofs_.push_back(Dof(dim, 0, coords, { 1 }, Value));
121  }
122  else
123  {
124  // Create vector of barycentric coordinates.
125  // First we make vector uvbc which contains barycentric coordinates
126  // multiplied by degree_, so that its values are unsigned ints.
127  // Then by counting the nonzero barycentric coordinates we can decide
128  // whether the dof lies on node, line, triangle or tetrahedron.
130  arma::uvec ubc = arma::zeros<arma::uvec>(dim+1);
131  ubc[0] = degree_;
132  bool finish = false;
133  do {
134  uvbc.push_back(ubc);
135  if (ubc[0] > 0)
136  {
137  // by default increment the first coordinate
138  ubc[1] += 1;
139  ubc[0] -= 1;
140  }
141  else
142  {
143  // if sum of coordinates is maximal (0-th coordinate is zero)
144  // then find first nonzero coordinate,
145  // set it to zero, and increment the following coordinate.
146  unsigned int c = 1;
147  while (ubc[c] == 0) c++;
148  // if the first nonzero coordinate is the last but one, we reach the end
149  if (c == dim) finish = true;
150  else {
151  ubc[0] = ubc[c]-1;
152  ubc[c] = 0;
153  ubc[c+1] += 1;
154  }
155  }
156  } while (!finish);
157 
158  // define dofs
159  for (auto ubc : uvbc)
160  {
161  // we define nodal dof:
162  // coords = barycentric coordinates of the support point,
163  // coefs = 1 (dof value = function value at the point)
164 
165  // count nonzero coordinates in ubc: 1=>nodal dof, 2=>dof on line etc.
166  arma::uvec nonzeros = find(ubc);
167  // convert "unsigned barycentric coordinates" to real ones
168  arma::vec coords = arma::conv_to<arma::vec>::from(ubc);
169  coords /= degree_;
170 
171  // find index of n-face
172  std::pair<unsigned int, unsigned int> zeros = RefElement<dim>::zeros_positions(coords);
173  unsigned int n_face_idx = -1;
174  switch (dim-zeros.first) {
175  case 0:
176  n_face_idx = RefElement<dim>::template topology_idx<0>(zeros.second);
177  break;
178  case 1:
179  n_face_idx = RefElement<dim>::template topology_idx<1>(zeros.second);
180  break;
181  case 2:
182  n_face_idx = RefElement<dim>::template topology_idx<2>(zeros.second);
183  break;
184  case 3:
185  n_face_idx = RefElement<dim>::template topology_idx<3>(zeros.second);
186  break;
187  }
188  this->dofs_.push_back(Dof(nonzeros.size()-1, n_face_idx, coords, { 1 }, Value));
189  }
190  }
191 }
192 
193 
194 
195 
196 template<unsigned int dim>
197 FE_P<dim>::FE_P(unsigned int degree)
198  : FiniteElement<dim>(),
199  degree_(degree)
200 {
201  this->function_space_ = std::make_shared<PolynomialSpace>(degree,dim);
202 
203  init_dofs();
204 
205  this->setup_components();
206 
207  this->compute_node_matrix();
208 }
209 
210 
211 
212 
213 
214 
215 
216 
217 
218 
219 
220 
221 
222 
223 
224 template<unsigned int dim>
225 FE_P_disc<dim>::FE_P_disc(unsigned int degree)
226  : FE_P<dim>(degree)
227 {
228  // all dofs are "inside" the cell (not shared with neighbours)
229  for (unsigned int i=0; i<this->dofs_.size(); i++)
230  this->dofs_[i].dim = dim;
231 
232  this->setup_components();
233 
234  this->compute_node_matrix();
235 }
236 
237 
238 
239 
240 
241 
242 template<unsigned int dim>
244 : FiniteElement<dim>()
245 {
246  this->function_space_ = std::make_shared<PolynomialSpace>(1,dim);
247 
248  if (dim == 0)
249  {
250  this->dofs_.push_back(Dof(0, 0, { 1 }, { 1 }, Value));
251  }
252  else
253  {
254  arma::vec::fixed<dim> sp; // support point
255  for (unsigned int sid=0; sid<RefElement<dim>::n_sides; ++sid)
256  {
257  sp.fill(0);
258  for (unsigned int i=0; i<RefElement<dim>::n_nodes_per_side; ++i)
261  // barycentric coordinates
262  arma::vec::fixed<dim+1> bsp;
263  bsp.subvec(1,dim) = sp;
264  bsp[0] = 1. - arma::sum(sp);
265  this->dofs_.push_back(Dof(dim-1, sid, bsp, { 1 }, Value));
266  }
267  }
268 
269  this->setup_components();
270  this->compute_node_matrix();
271 }
272 
273 
274 
275 
276 template<unsigned int dim>
278 : FiniteElement<dim>()
279 {
280  this->function_space_ = std::make_shared<PolynomialSpace>(1,dim);
281 
282  if (dim == 0)
283  {
284  this->dofs_.push_back(Dof(0, 0, { 1 }, { 1 }, Value));
285  }
286  else
287  {
288  arma::vec::fixed<dim> sp; // support point
289  for (unsigned int sid=0; sid<RefElement<dim>::n_sides; ++sid)
290  {
291  sp.fill(0);
292  for (unsigned int i=0; i<RefElement<dim>::n_nodes_per_side; ++i)
295  // barycentric coordinates
296  arma::vec::fixed<dim+1> bsp;
297  bsp.subvec(1,dim) = sp;
298  bsp[0] = 1. - arma::sum(sp);
299  this->dofs_.push_back(Dof(dim, 0, bsp, { 1 }, Value));
300  }
301  }
302 
303  this->setup_components();
304  this->compute_node_matrix();
305 }
306 
307 
308 
309 
310 
311 template class FE_P<0>;
312 template class FE_P<1>;
313 template class FE_P<2>;
314 template class FE_P<3>;
315 
316 
317 template class FE_P_disc<0>;
318 template class FE_P_disc<1>;
319 template class FE_P_disc<2>;
320 template class FE_P_disc<3>;
321 
322 
323 template class FE_CR<0>;
324 template class FE_CR<1>;
325 template class FE_CR<2>;
326 template class FE_CR<3>;
327 
328 
329 template class FE_CR_disc<0>;
330 template class FE_CR_disc<1>;
331 template class FE_CR_disc<2>;
332 template class FE_CR_disc<3>;
333 
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
#define ASSERT_LE(a, b)
Definition of comparative assert macro (Less or Equal) only for debug mode.
Definition: asserts.hh:309
Discontinuos Crouzeix-Raviart finite element on dim dimensional simplex.
Definition: fe_p.hh:138
FE_CR_disc()
Definition: fe_p.cc:277
Crouzeix-Raviart finite element on dim dimensional simplex.
Definition: fe_p.hh:124
FE_CR()
Definition: fe_p.cc:243
Discontinuous Lagrangean finite element on dim dimensional simplex.
Definition: fe_p.hh:108
FE_P_disc(unsigned int degree)
Constructor.
Definition: fe_p.cc:225
Conforming Lagrangean finite element on dim dimensional simplex.
Definition: fe_p.hh:86
void init_dofs()
Definition: fe_p.cc:112
FE_P(unsigned int degree)
Constructor.
Definition: fe_p.cc:197
Abstract class for the description of a general finite element on a reference simplex in dim dimensio...
std::shared_ptr< FunctionSpace > function_space_
Function space defining the FE.
void setup_components()
Initialize vectors with information about components of basis functions.
virtual void compute_node_matrix()
Initializes the node_matrix for computing the coefficients of the shape functions in the raw basis of...
std::vector< Dof > dofs_
Set of degrees of freedom (functionals) defining the FE.
unsigned int space_dim_
Space dimension of function arguments (i.e. 1, 2 or 3).
unsigned int dim() const override
Dimension of function space (number of basis functions).
Definition: fe_p.hh:57
std::vector< arma::uvec > powers
Coefficients of basis functions.
Definition: fe_p.hh:70
PolynomialSpace(unsigned int degree, unsigned int dim)
Constructor.
Definition: fe_p.cc:25
double basis_value(unsigned int basis_index, const arma::vec &point, unsigned int comp_index=0) const override
Value of the i th basis function at point point.
Definition: fe_p.cc:62
const arma::vec basis_grad(unsigned int basis_index, const arma::vec &point, unsigned int comp_index=0) const override
Gradient of the i th basis function at point point.
Definition: fe_p.cc:78
static std::pair< unsigned int, unsigned int > zeros_positions(const BaryPoint &barycentric, double tolerance=std::numeric_limits< double >::epsilon() *2)
Definition: ref_element.cc:294
Definitions of basic Lagrangean finite elements with polynomial shape functions.
@ Value
ArmaVec< double, N > vec
Definition: armor.hh:885
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.