Flow123d  release_3.0.0-1113-g35b167a
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 const double PolynomialSpace::basis_value(unsigned int i,
63  const arma::vec &point,
64  unsigned int comp_index
65  ) const
66 {
67  ASSERT(comp_index == 0);
68  OLD_ASSERT(i<=powers.size(), "Index of basis function is out of range.");
69  ASSERT(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 
75  return v;
76 }
77 
78 
79 const arma::vec PolynomialSpace::basis_grad(unsigned int i,
80  const arma::vec &p,
81  unsigned int comp_index
82  ) const
83 {
84  ASSERT(comp_index == 0);
85  OLD_ASSERT(i<=powers.size(), "Index of basis function is out of range.");
86 
87  arma::vec grad(this->space_dim_);
88 
89  for (unsigned int j=0; j<this->space_dim_; j++)
90  {
91  grad[j] = powers[i][j];
92  if (powers[i][j] == 0) continue;
93 
94  for (unsigned int k=0; k<this->space_dim_; k++)
95  {
96  grad[j] *= pow(p[k], (int) (k==j?powers[i][k]-1:powers[i][k]));
97  }
98  }
99  return grad;
100 }
101 
102 
103 
104 
105 
106 
107 
108 
109 
110 
111 
112 template<unsigned int dim>
114 {
115  if (degree_ == 0 || dim == 0)
116  {
117  // we define nodal dof:
118  // coords = barycentric coordinates of the support point,
119  // coefs = 1 (dof value = function value at the point)
120  arma::vec coords = arma::ones<arma::vec>(dim+1)/(dim+1);
121  this->dofs_.push_back(Dof(dim, 0, coords, { 1 }, Value));
122  }
123  else
124  {
125  // Create vector of barycentric coordinates.
126  // First we make vector uvbc which contains barycentric coordinates
127  // multiplied by degree_, so that its values are unsigned ints.
128  // Then by counting the nonzero barycentric coordinates we can decide
129  // whether the dof lies on node, line, triangle or tetrahedron.
131  arma::uvec ubc = arma::zeros<arma::uvec>(dim+1);
132  ubc[0] = degree_;
133  bool finish = false;
134  do {
135  uvbc.push_back(ubc);
136  if (ubc[0] > 0)
137  {
138  // by default increment the first coordinate
139  ubc[1] += 1;
140  ubc[0] -= 1;
141  }
142  else
143  {
144  // if sum of coordinates is maximal (0-th coordinate is zero)
145  // then find first nonzero coordinate,
146  // set it to zero, and increment the following coordinate.
147  unsigned int c = 1;
148  while (ubc[c] == 0) c++;
149  // if the first nonzero coordinate is the last but one, we reach the end
150  if (c == dim) finish = true;
151  else {
152  ubc[0] = ubc[c]-1;
153  ubc[c] = 0;
154  ubc[c+1] += 1;
155  }
156  }
157  } while (!finish);
158 
159  // define dofs
160  for (auto ubc : uvbc)
161  {
162  // we define nodal dof:
163  // coords = barycentric coordinates of the support point,
164  // coefs = 1 (dof value = function value at the point)
165 
166  // count nonzero coordinates in ubc: 1=>nodal dof, 2=>dof on line etc.
167  arma::uvec nonzeros = find(ubc);
168  // convert "unsigned barycentric coordinates" to real ones
169  arma::vec coords = arma::conv_to<arma::vec>::from(ubc);
170  coords /= degree_;
171 
172  // find index of n-face
173  std::pair<unsigned int, unsigned int> zeros = RefElement<dim>::zeros_positions(coords);
174  unsigned int n_face_idx;
175  switch (dim-zeros.first) {
176  case 0:
177  n_face_idx = RefElement<dim>::template topology_idx<0>(zeros.second);
178  break;
179  case 1:
180  n_face_idx = RefElement<dim>::template topology_idx<1>(zeros.second);
181  break;
182  case 2:
183  n_face_idx = RefElement<dim>::template topology_idx<2>(zeros.second);
184  break;
185  case 3:
186  n_face_idx = RefElement<dim>::template topology_idx<3>(zeros.second);
187  break;
188  }
189  this->dofs_.push_back(Dof(nonzeros.size()-1, n_face_idx, coords, { 1 }, Value));
190  }
191  }
192 }
193 
194 
195 
196 
197 template<unsigned int dim>
198 FE_P<dim>::FE_P(unsigned int degree)
199  : FiniteElement<dim>(),
200  degree_(degree)
201 {
202  this->function_space_ = std::make_shared<PolynomialSpace>(degree,dim);
203 
204  init_dofs();
205 
206  this->setup_components();
207 
208  this->compute_node_matrix();
209 }
210 
211 
212 
213 
214 
215 
216 
217 
218 
219 
220 
221 
222 
223 
224 
225 template<unsigned int dim>
226 FE_P_disc<dim>::FE_P_disc(unsigned int degree)
227  : FE_P<dim>(degree)
228 {
229  // all dofs are "inside" the cell (not shared with neighbours)
230  for (unsigned int i=0; i<this->dofs_.size(); i++)
231  this->dofs_[i].dim = dim;
232 
233  this->setup_components();
234 
235  this->compute_node_matrix();
236 }
237 
238 
239 
240 
241 
242 
243 template<unsigned int dim>
245 : FiniteElement<dim>()
246 {
247  this->function_space_ = std::make_shared<PolynomialSpace>(1,dim);
248 
249  if (dim == 0)
250  {
251  this->dofs_.push_back(Dof(0, 0, { 1 }, { 1 }, Value));
252  }
253  else
254  {
255  arma::vec::fixed<dim> sp; // support point
256  for (unsigned int sid=0; sid<RefElement<dim>::n_sides; ++sid)
257  {
258  sp.fill(0);
259  for (unsigned int i=0; i<RefElement<dim>::n_nodes_per_side; ++i)
262  // barycentric coordinates
263  arma::vec::fixed<dim+1> bsp;
264  bsp.subvec(1,dim) = sp;
265  bsp[0] = 1. - arma::sum(sp);
266  this->dofs_.push_back(Dof(dim-1, sid, bsp, { 1 }, Value));
267  }
268  }
269 
270  this->setup_components();
271  this->compute_node_matrix();
272 }
273 
274 
275 
276 
277 template<unsigned int dim>
279 : FiniteElement<dim>()
280 {
281  this->function_space_ = std::make_shared<PolynomialSpace>(1,dim);
282 
283  if (dim == 0)
284  {
285  this->dofs_.push_back(Dof(0, 0, { 1 }, { 1 }, Value));
286  }
287  else
288  {
289  arma::vec::fixed<dim> sp; // support point
290  for (unsigned int sid=0; sid<RefElement<dim>::n_sides; ++sid)
291  {
292  sp.fill(0);
293  for (unsigned int i=0; i<RefElement<dim>::n_nodes_per_side; ++i)
296  // barycentric coordinates
297  arma::vec::fixed<dim+1> bsp;
298  bsp.subvec(1,dim) = sp;
299  bsp[0] = 1. - arma::sum(sp);
300  this->dofs_.push_back(Dof(dim, 0, bsp, { 1 }, Value));
301  }
302  }
303 
304  this->setup_components();
305  this->compute_node_matrix();
306 }
307 
308 
309 
310 
311 
312 template class FE_P<0>;
313 template class FE_P<1>;
314 template class FE_P<2>;
315 template class FE_P<3>;
316 
317 
318 template class FE_P_disc<0>;
319 template class FE_P_disc<1>;
320 template class FE_P_disc<2>;
321 template class FE_P_disc<3>;
322 
323 
324 template class FE_CR<0>;
325 template class FE_CR<1>;
326 template class FE_CR<2>;
327 template class FE_CR<3>;
328 
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:79
PolynomialSpace(unsigned int degree, unsigned int dim)
Constructor.
Definition: fe_p.cc:25
const 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
FE_P_disc(unsigned int degree)
Constructor.
Definition: fe_p.cc:226
static std::pair< unsigned int, unsigned int > zeros_positions(const BaryPoint &barycentric, double tolerance=std::numeric_limits< double >::epsilon()*2)
Definition: ref_element.cc:393
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library...
Definition: asserts.hh:346
std::vector< arma::uvec > powers
Coefficients of basis functions.
Definition: fe_p.hh:73
#define OLD_ASSERT(...)
Definition: global_defs.h:131
FE_CR()
Definition: fe_p.cc:244
std::vector< Dof > dofs_
Set of degrees of freedom (functionals) defining the FE.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
FE_P(unsigned int degree)
Constructor.
Definition: fe_p.cc:198
void init_dofs()
Definition: fe_p.cc:113
virtual void compute_node_matrix()
Initializes the node_matrix for computing the coefficients of the shape functions in the raw basis of...
std::shared_ptr< FunctionSpace > function_space_
Function space defining the FE.
void setup_components()
Initialize vectors with information about components of basis functions.
Conforming Lagrangean finite element on dim dimensional simplex.
Definition: fe_p.hh:88
const unsigned int dim() const override
Dimension of function space (number of basis functions).
Definition: fe_p.hh:60
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc...
unsigned int space_dim_
Space dimension of function arguments (i.e. 1, 2 or 3).
Abstract class for the description of a general finite element on a reference simplex in dim dimensio...
Crouzeix-Raviart finite element on dim dimensional simplex.
Definition: fe_p.hh:126
FE_CR_disc()
Definition: fe_p.cc:278
const unsigned int degree_
Max. degree of polynomials.
Definition: fe_p.hh:65