Flow123d  release_3.0.0-879-gd0f69e9
fe_system.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_system.cc
15  * @brief Class FESystem for compound finite elements.
16  * @author Jan Stebel
17  */
18 
19 #include "mesh/side_impl.hh"
20 #include "fem/fe_system.hh"
21 #include "system/global_defs.h"
22 #include "quadrature/quadrature.hh"
23 
24 using namespace std;
25 
26 
27 unsigned int count_components(const vector<shared_ptr<FunctionSpace> > &fs_vector)
28 {
29  unsigned int n_comp = 0;
30  for (auto fs : fs_vector)
31  n_comp += fs->n_components();
32 
33  return n_comp;
34 }
35 
36 
37 unsigned int check_spacedim(const vector<shared_ptr<FunctionSpace> > &fs_vector)
38 {
39  ASSERT_DBG(fs_vector.size() > 0);
40  unsigned int space_dim = fs_vector[0]->space_dim();
41  for (auto fs : fs_vector)
42  ASSERT_DBG(fs->space_dim() == space_dim).error("FunctionSpace space_dim mismatch.");
43 
44  return space_dim;
45 }
46 
47 
48 FESystemFunctionSpace::FESystemFunctionSpace(const vector<shared_ptr<FunctionSpace> > &fs_vector)
49  : FunctionSpace(check_spacedim(fs_vector), count_components(fs_vector)),
50  fs_(fs_vector)
51 {
52  dim_ = 0;
53  unsigned int fe_index = 0;
54  unsigned int comp_offset = 0;
55  for (auto fs : fs_vector)
56  {
57  for (unsigned int i=0; i<fs->dim(); i++)
58  dof_indices_.push_back(DofComponentData(fe_index, i, comp_offset));
59 
60  dim_ += fs->dim();
61  fe_index++;
62  comp_offset += fs->n_components();
63  }
64 }
65 
66 
67 const double FESystemFunctionSpace::basis_value(unsigned int i,
68  const arma::vec &p,
69  unsigned int comp) const
70 {
71  ASSERT_DBG(i < dim_).error("Index of basis function is out of range.");
72  ASSERT_DBG(comp < n_components()).error("Index of component is out of range.");
73 
74  // component index in the base FE
75  int l_comp = comp-dof_indices_[i].component_offset;
76  if (l_comp >= 0 && l_comp < fs_[dof_indices_[i].fe_index]->n_components())
77  return fs_[dof_indices_[i].fe_index]->basis_value(dof_indices_[i].basis_index, p, l_comp);
78  else
79  return 0;
80 }
81 
82 
83 const arma::vec FESystemFunctionSpace::basis_grad(const unsigned int i,
84  const arma::vec &p,
85  const unsigned int comp) const
86 {
87  ASSERT_DBG(i < dim_).error("Index of basis function is out of range.");
88  ASSERT_DBG(comp < n_components()).error("Index of component is out of range.");
89 
90  // component index in the base FE
91  int l_comp = comp-dof_indices_[i].component_offset;
92  if (l_comp >= 0 && l_comp < fs_[dof_indices_[i].fe_index]->n_components())
93  return fs_[dof_indices_[i].fe_index]->basis_grad(dof_indices_[i].basis_index, p, l_comp);
94  else
95  return arma::zeros(space_dim());
96 }
97 
98 
99 
100 
101 template<unsigned int dim>
103 {
104  OLD_ASSERT(fe->n_components() == 1, "FEVectorContravariant and FEVectorPiola can only be created from scalar FE.");
106  "This constructor can be used only for FEVectorContravariant or FEVectorPiola.");
107 
108  FiniteElement<dim>::init(false, t);
110  initialize();
111 }
112 
113 
114 template<unsigned int dim>
115 FESystem<dim>::FESystem(const std::shared_ptr<FiniteElement<dim> > &fe, FEType t, unsigned int n)
116 {
118  "This constructor can be used only for FEVector, FETensor or FEMixedSystem.");
119  OLD_ASSERT(fe->n_components() == 1 || t == FEType::FEMixedSystem,
120  "FEVector and FETensor can only be created from scalar FE.");
121 
122  FiniteElement<dim>::init(false, t);
124  initialize();
125 }
126 
127 
128 template<unsigned int dim>
130 {
132  for (std::shared_ptr<FiniteElement<dim> > fe_object : fe)
133  fe_.push_back(fe_object);
134  initialize();
135 }
136 
137 
138 
139 template<unsigned int dim>
141 {
142  unsigned int fe_index = 0;
143  unsigned int comp_offset = 0;
145  // for each base FE add components, support points, and other
146  // information to the system
147  for (auto fe : fe_)
148  {
149  switch (fe->type_)
150  {
151  case FEType::FEScalar:
152  scalar_components_.push_back(comp_offset);
153  break;
154  case FEType::FEVector:
157  vector_components_.push_back(comp_offset);
158  break;
159  case FEType::FETensor:
160  tensor_components_.push_back(comp_offset);
161  break;
162  default:
163  OLD_ASSERT(false, "Not implemented.");
164  break;
165  }
166 
167  fe_index++;
168  comp_offset += fe->n_components();
169  fs_vector.push_back(shared_ptr<FunctionSpace>(fe->function_space_));
170  }
171 
172  this->function_space_ = make_shared<FESystemFunctionSpace>(fs_vector);
173 
174  double dof_index = 0;
175  comp_offset = 0;
176  for (auto fe : fe_)
177  {
178  for (unsigned int i=0; i<fe->n_dofs(); i++)
179  {
180  arma::vec coefs(this->function_space_->n_components());
181  coefs.subvec(comp_offset, comp_offset+fe->dof(i).coefs.size()-1) = fe->dof(i).coefs;
182  this->dofs_.push_back(Dof(fe->dof(i).dim, fe->dof(i).n_face_idx, fe->dof(i).coords, coefs, fe->dof(i).type));
183  dof_index++;
184  }
185  comp_offset += fe->n_components();
186  }
187 
188 // if (this->is_primitive_)
189 // {
190 // double dof_index = 0;
191 // // add index of nonzero component for each dof in FESystem
192 // for (auto fe : fe_)
193 // {
194 // for (int i=0; i<fe->n_dofs(); ++i)
195 // this->component_indices_.push_back(fe_dof_indices_[dof_index++].fe_index);
196 // }
197 // } else {
198  dof_index = 0;
199  comp_offset = 0;
200  // add footprint of nonzero components for each dof in FESystem
201  for (auto fe : fe_)
202  {
203  for (int i=0; i<fe->n_dofs(); ++i)
204  {
205  std::vector<bool> nonzeros(this->function_space_->n_components(), false);
206  for (unsigned int c=0; c<fe->n_components(); c++)\
207  nonzeros[comp_offset+c] = fe->get_nonzero_components(i)[c];
208  this->nonzero_components_.push_back(nonzeros);
209  dof_index++;
210  }
211  comp_offset += fe->n_components();
212  }
213 // }
214 
215  compute_node_matrix();
216 }
217 
218 
219 
220 
221 template<unsigned int dim> inline
223 {
224  UpdateFlags f = flags;
225 
226  for (auto fe : fe_)
227  f |= fe->update_each(flags);
228 
229  return f;
230 }
231 
232 
233 template<unsigned int dim>
235 {
236  // form the node_matrix of the FESystem as block diagonal matrix
237  // composed of node_matrices of each base FE class
238 
239  this->node_matrix.resize(this->dofs_.size(), this->dofs_.size());
240 
241  unsigned int offset = 0;
242  for (unsigned int i = 0; i < fe_.size(); i++)
243  {
244  this->node_matrix.submat(offset, offset, offset+fe_[i]->n_dofs()-1, offset+fe_[i]->n_dofs()-1)
245  = fe_[i]->node_matrix;
246 
247  offset += fe_[i]->n_dofs();
248  }
249 }
250 
251 
252 
253 
254 
255 
256 
257 
258 template class FESystem<0>;
259 template class FESystem<1>;
260 template class FESystem<2>;
261 template class FESystem<3>;
262 
263 
264 
const unsigned int space_dim() const
Getter for space dimension.
void initialize()
Initialization of the internal structures from the vector of base FE.
Definition: fe_system.cc:140
UpdateFlags
Enum type UpdateFlags indicates which quantities are to be recomputed on each finite element cell...
Definition: update_flags.hh:67
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_system.cc:67
const unsigned int dim() const override
Dimension of function space (number of basis functions).
Definition: fe_system.hh:74
std::vector< DofComponentData > dof_indices_
Indices of basis functions relative to the sub-spaces.
Definition: fe_system.hh:84
unsigned int dim_
Number of basis functions.
Definition: fe_system.hh:87
unsigned int check_spacedim(const vector< shared_ptr< FunctionSpace > > &fs_vector)
Definition: fe_system.cc:37
Class FESystem for compound finite elements.
FEType
std::vector< std::shared_ptr< FunctionSpace > > fs_
Function spaces that are put together.
Definition: fe_system.hh:76
#define OLD_ASSERT(...)
Definition: global_defs.h:131
Global macros to enhance readability and debugging, general constants.
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_system.cc:83
Compound finite element on dim dimensional simplex.
Definition: fe_system.hh:99
Basic definitions of numerical quadrature rules.
void init(bool primitive=true, FEType type=FEScalar)
Clears all internal structures.
const unsigned int n_components() const
Getter for number of components.
UpdateFlags update_each(UpdateFlags flags) override
Decides which additional quantities have to be computed for each cell.
Definition: fe_system.cc:222
unsigned int count_components(const vector< shared_ptr< FunctionSpace > > &fs_vector)
Definition: fe_system.cc:27
#define ASSERT_DBG(expr)
Definition: asserts.hh:349
Abstract class for the description of a general finite element on a reference simplex in dim dimensio...
void compute_node_matrix() override
Initializes the node_matrix for computing the coefficients of the shape functions in the raw basis of...
Definition: fe_system.cc:234
FESystemFunctionSpace(const std::vector< std::shared_ptr< FunctionSpace > > &fs_vector)
Constructor.
Definition: fe_system.cc:48
FESystem(std::shared_ptr< FiniteElement< dim > > fe, FEType t)
Constructor for FEVectorContravariant and FEVectorPiola.
Definition: fe_system.cc:102