Flow123d  JS_before_hm-1621-g63a12c7
finite_element.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 finite_element.cc
15  * @brief Abstract class for description of finite elements.
16  * @author Jan Stebel
17  */
18 
19 #include "system/system.hh"
20 #include "quadrature/quadrature.hh"
21 #include "fem/dofhandler.hh"
22 #include "fem/finite_element.hh"
23 #include "fem/fe_values.hh"
24 #include "fem/fe_system.hh"
25 
26 
27 
28 using namespace std;
29 
30 
31 
32 
33 
34 template<class FS> double Dof::evaluate(const FS &function_space,
35  unsigned int basis_idx) const
36 {
37  // Check that FS is derived from FunctionSpace.
38  static_assert(std::is_base_of<FunctionSpace, FS>::value, "FS must be derived from FunctionSpace.");
39 
40  // We cannot evaluate dof on dim-dimensional n-face if the function space lies on lower-dimensional n-face.
41  ASSERT(function_space.space_dim()+1 == coords.size());
42 
43  switch (type)
44  {
45  case Value:
46  {
47  // evaluate basis function and return the linear combination of components
48  arma::vec vec_value(function_space.n_components());
49 
50  // drop off the 0-th barycentric coordinate
51  // in case of space_dim=0 subvec does not work
52  arma::vec f_coords(function_space.space_dim());
53  if (function_space.space_dim() > 0)
54  f_coords = coords.subvec(1,coords.size()-1);
55  for (unsigned int c=0; c<function_space.n_components(); c++)
56  vec_value[c] = function_space.basis_value(basis_idx, f_coords, c);
57  return dot(coefs, vec_value);
58  break;
59  }
60 
61  default:
62  OLD_ASSERT(false, "Dof evaluation not implemented for this type.");
63  }
64  return 0;
65 }
66 
67 
68 
69 
70 
71 
72 
73 template<unsigned int dim>
75  : function_space_(nullptr)
76 {
77  init();
78 }
79 
80 template<unsigned int dim>
81 void FiniteElement<dim>::init(bool primitive, FEType type)
82 {
83  dofs_.clear();
84  is_primitive_ = primitive;
85  type_ = type;
86 }
87 
88 
89 template<unsigned int dim>
91 {
92  component_indices_.resize(dofs_.size(), 0);
93  nonzero_components_.resize(dofs_.size(), { true });
94 }
95 
96 
97 template<unsigned int dim> inline
99 {
100  arma::mat M(dofs_.size(), dofs_.size());
101 
102  for (unsigned int i = 0; i < dofs_.size(); i++)
103  for (unsigned int j = 0; j < dofs_.size(); j++) {
104  M(j, i) = dofs_[i].evaluate(*function_space_, j);
105 
106  }
107  node_matrix = arma::inv(M);
108 }
109 
110 
111 template<unsigned int dim>
112 double FiniteElement<dim>::shape_value(const unsigned int i,
113  const arma::vec::fixed<dim> &p,
114  const unsigned int comp) const
115 {
116  ASSERT_DBG( comp < n_components() );
117  ASSERT_DBG( i < dofs_.size()).error("Index of basis function is out of range.");
118 
119  double value = 0;
120  for (unsigned int j=0; j<function_space_->dim(); j++)
121  value += function_space_->basis_value(j, p, comp) * node_matrix(i,j);
122 
123  return value;
124 }
125 
126 template<unsigned int dim>
127 arma::vec::fixed<dim> FiniteElement<dim>::shape_grad(const unsigned int i,
128  const arma::vec::fixed<dim> &p,
129  const unsigned int comp) const
130 {
131  ASSERT_DBG( comp < n_components() );
132  ASSERT_DBG( i < dofs_.size()).error("Index of basis function is out of range.");
133 
134  // uword is a typedef for an unsigned integer type; it is used for matrix indices as well as all internal counters and loops
135  // sword is a typedef for a signed integer type
136  arma::vec grad((arma::uword)dim);
137  grad.zeros();
138  for (unsigned int j=0; j<function_space_->dim(); j++)
139  grad += function_space_->basis_grad(j, p, comp) * node_matrix(i,j);
140 
141  return grad;
142 }
143 
144 
145 template<unsigned int dim> inline
147 {
148  UpdateFlags f = flags;
149 
150  switch (type_)
151  {
152  case FEScalar:
153  case FEVector:
154  case FETensor:
155  if (flags & update_gradients)
157  break;
159  if (flags & update_values)
160  f |= update_jacobians;
161  if (flags & update_gradients)
163  break;
164  case FEVectorPiola:
165  if (flags & update_values)
167  if (flags & update_gradients)
169  break;
170  default:;
171  }
172 
173  return f;
174 }
175 
176 
177 template<unsigned int dim>
178 unsigned int FiniteElement<dim>::n_space_components(unsigned int spacedim)
179 {
180  switch (type_) {
181  case FEScalar:
182  return 1;
183  break;
184  case FEVector:
186  case FEVectorPiola:
187  return spacedim;
188  break;
189  case FETensor:
190  return spacedim*spacedim;
191  break;
192  case FEMixedSystem:
193  const FESystem<dim> *fe_sys = dynamic_cast<const FESystem<dim>*>(this);
194  ASSERT_DBG(fe_sys != nullptr).error("Mixed system must be represented by FESystem.");
195  return fe_sys->get_scalar_components().size()
196  +fe_sys->get_vector_components().size()*spacedim
197  +fe_sys->get_tensor_components().size()*spacedim*spacedim;
198  break;
199  }
200 
201  // should be never reached
202  ASSERT(0).error("Unknown type of FiniteElement.");
203  return 0;
204 }
205 
206 
207 template<unsigned int dim>
210  points.resize(0);
211  for(auto dof : this->dofs_)
212  points.push_back(dof.coords);
213  return points;
214 }
215 
216 
217 
218 
219 
220 template class FiniteElement<0>;
221 template class FiniteElement<1>;
222 template class FiniteElement<2>;
223 template class FiniteElement<3>;
224 
225 
FiniteElement()
Constructor.
Shape function values.
Definition: update_flags.hh:87
UpdateFlags
Enum type UpdateFlags indicates which quantities are to be recomputed on each finite element cell...
Definition: update_flags.hh:67
Declaration of class which handles the ordering of degrees of freedom (dof) and mappings between loca...
std::vector< unsigned int > get_tensor_components() const
Definition: fe_system.hh:137
arma::vec coords
Barycentric coordinates.
ArmaVec< double, N > vec
Definition: armor.hh:885
arma::mat node_matrix
Matrix that determines the coefficients of the raw basis functions from the values at the support poi...
Determinant of the Jacobian.
bool is_primitive_
Primitive FE is using componentwise shape functions, i.e. only one component is nonzero for each shap...
arma::vec::fixed< dim > shape_grad(const unsigned int i, const arma::vec::fixed< dim > &p, const unsigned int comp=0) const
Calculates the comp-th component of the gradient of the i-th shape function at the point p on the ref...
std::vector< std::vector< bool > > nonzero_components_
Footprints of nonzero components of shape functions.
Class FEValues calculates finite element data on the actual cells such as shape function values...
Class FESystem for compound finite elements.
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library...
Definition: asserts.hh:347
FEType
std::vector< unsigned int > get_scalar_components() const
Definition: fe_system.hh:131
static constexpr bool value
Definition: json.hpp:87
Volume element.
#define OLD_ASSERT(...)
Definition: global_defs.h:108
ArmaMat< double, N, M > mat
Definition: armor.hh:888
std::vector< unsigned int > get_vector_components() const
Definition: fe_system.hh:134
double shape_value(const unsigned int i, const arma::vec::fixed< dim > &p, const unsigned int comp=0) const
Calculates the value of the comp-th component of the i-th shape function at the point p on the refere...
Compound finite element on dim dimensional simplex.
Definition: fe_system.hh:101
std::vector< Dof > dofs_
Set of degrees of freedom (functionals) defining the FE.
virtual UpdateFlags update_each(UpdateFlags flags)
Decides which additional quantities have to be computed for each cell.
unsigned int n_components() const
Returns numer of components of the basis function.
virtual void compute_node_matrix()
Initializes the node_matrix for computing the coefficients of the shape functions in the raw basis of...
Basic definitions of numerical quadrature rules.
Shape function gradients.
Definition: update_flags.hh:95
void init(bool primitive=true, FEType type=FEScalar)
Clears all internal structures.
unsigned int n_space_components(unsigned int spacedim)
Number of components of FE in a mapped space with dimension spacedim.
std::shared_ptr< FunctionSpace > function_space_
Function space defining the FE.
const Dof & dof(unsigned int i) const
Returns i -th degree of freedom.
void setup_components()
Initialize vectors with information about components of basis functions.
double evaluate(const FS &function_space, unsigned int basis_idx) const
Evaulate dof for basis function of given function space.
virtual std::vector< arma::vec::fixed< dim+1 > > dof_points() const
#define ASSERT_DBG(expr)
Abstract class for description of finite elements.
std::vector< unsigned int > component_indices_
Indices of nonzero components of shape functions (for primitive FE).
FEType type_
Type of FiniteElement.