Flow123d  master-f44eb46
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  ASSERT_PERMANENT(false).error("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( comp < n_components() );
117  ASSERT( 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( comp < n_components() );
132  ASSERT( 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(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_PERMANENT(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 
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
double evaluate(const FS &function_space, unsigned int basis_idx) const
Evaulate dof for basis function of given function space.
Compound finite element on dim dimensional simplex.
Definition: fe_system.hh:102
std::vector< unsigned int > get_vector_components() const
Definition: fe_system.hh:134
std::vector< unsigned int > get_scalar_components() const
Definition: fe_system.hh:131
std::vector< unsigned int > get_tensor_components() const
Definition: fe_system.hh:137
Abstract class for the description of a general finite element on a reference simplex in dim dimensio...
unsigned int n_space_components(unsigned int spacedim)
Number of components of FE in a mapped space with dimension spacedim.
virtual UpdateFlags update_each(UpdateFlags flags)
Decides which additional quantities have to be computed for each cell.
virtual std::vector< arma::vec::fixed< dim+1 > > dof_points() const
void init(bool primitive=true, FEType type=FEScalar)
Clears all internal structures.
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...
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...
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...
FiniteElement()
Constructor.
Declaration of class which handles the ordering of degrees of freedom (dof) and mappings between loca...
Class FESystem for compound finite elements.
Class FEValues calculates finite element data on the actual cells such as shape function values,...
Abstract class for description of finite elements.
FEType
@ FEScalar
@ FEMixedSystem
@ FEVectorPiola
@ FETensor
@ FEVectorContravariant
@ FEVector
@ Value
static constexpr bool value
Definition: json.hpp:87
ArmaMat< double, N, M > mat
Definition: armor.hh:888
ArmaVec< double, N > vec
Definition: armor.hh:885
Basic definitions of numerical quadrature rules.
UpdateFlags
Enum type UpdateFlags indicates which quantities are to be recomputed on each finite element cell.
Definition: update_flags.hh:68
@ update_volume_elements
Determinant of the Jacobian.
@ update_values
Shape function values.
Definition: update_flags.hh:87
@ update_jacobians
Volume element.
@ update_inverse_jacobians
Volume element.
@ update_gradients
Shape function gradients.
Definition: update_flags.hh:95