Flow123d  master-f44eb46
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/accessors.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(fs_vector.size() > 0);
40  unsigned int space_dim = fs_vector[0]->space_dim();
41  for (auto fs : fs_vector)
42  ASSERT(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 double FESystemFunctionSpace::basis_value(unsigned int i,
68  const arma::vec &p,
69  unsigned int comp) const
70 {
71  ASSERT(i < dim_).error("Index of basis function is out of range.");
72  ASSERT(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 < static_cast<int>(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(i < dim_).error("Index of basis function is out of range.");
88  ASSERT(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 < static_cast<int>(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  ASSERT_EQ(fe->n_components(), 1).error("FEVectorContravariant and FEVectorPiola can only be created from scalar FE.");
106  .error("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  .error("This constructor can be used only for FEVector, FETensor or FEMixedSystem.");
119  ASSERT(fe->n_components() == 1 || t == FEType::FEMixedSystem)
120  .error("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  ASSERT_PERMANENT(false).error("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 (uint 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 template<unsigned int dim>
235  points.resize(0);
236  for (unsigned int i = 0; i < fe_.size(); i++)
237  {
238  auto fe_points = fe_[i]->dof_points();
239  points.insert(points.end(), fe_points.begin(), fe_points.end());
240  }
241  return points;
242 }
243 
244 
245 template<unsigned int dim>
247 {
248  // form the node_matrix of the FESystem as block diagonal matrix
249  // composed of node_matrices of each base FE class
250 
251  this->node_matrix.resize(this->dofs_.size(), this->dofs_.size());
252 
253  unsigned int offset = 0;
254  for (unsigned int i = 0; i < fe_.size(); i++)
255  {
256  this->node_matrix.submat(offset, offset, offset+fe_[i]->n_dofs()-1, offset+fe_[i]->n_dofs()-1)
257  = fe_[i]->node_matrix;
258 
259  offset += fe_[i]->n_dofs();
260  }
261 }
262 
263 
264 
265 
266 template<unsigned int dim>
268 {
269  auto fs = dynamic_cast<FESystemFunctionSpace *>(this->function_space_.get());
270  std::vector<unsigned int> fe_dof_indices;
271 
272  for (unsigned int i=0; i<fs->dof_indices().size(); i++)
273  {
274  if (fs->dof_indices()[i].fe_index == fe_index)
275  fe_dof_indices.push_back(i);
276  }
277 
278  return fe_dof_indices;
279 }
280 
281 
282 
283 
284 
285 
286 
287 
288 template class FESystem<0>;
289 template class FESystem<1>;
290 template class FESystem<2>;
291 template class FESystem<3>;
292 
293 
294 
#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
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
FESystemFunctionSpace(const std::vector< std::shared_ptr< FunctionSpace > > &fs_vector)
Constructor.
Definition: fe_system.cc:48
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
std::vector< std::shared_ptr< FunctionSpace > > fs_
Function spaces that are put together.
Definition: fe_system.hh:78
unsigned int dim_
Number of basis functions.
Definition: fe_system.hh:89
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
std::vector< DofComponentData > dof_indices_
Indices of basis functions relative to the sub-spaces.
Definition: fe_system.hh:86
Compound finite element on dim dimensional simplex.
Definition: fe_system.hh:102
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:246
std::vector< unsigned int > fe_dofs(unsigned int fe_index)
Return dof indices belonging to given sub-FE.
Definition: fe_system.cc:267
virtual std::vector< arma::vec::fixed< dim+1 > > dof_points() const
Definition: fe_system.cc:233
UpdateFlags update_each(UpdateFlags flags) override
Decides which additional quantities have to be computed for each cell.
Definition: fe_system.cc:222
FESystem(std::shared_ptr< FiniteElement< dim > > fe, FEType t)
Constructor for FEVectorContravariant and FEVectorPiola.
Definition: fe_system.cc:102
void initialize()
Initialization of the internal structures from the vector of base FE.
Definition: fe_system.cc:140
Abstract class for the description of a general finite element on a reference simplex in dim dimensio...
void init(bool primitive=true, FEType type=FEScalar)
Clears all internal structures.
unsigned int n_components() const
Getter for number of components.
unsigned int space_dim() const
Getter for space dimension.
unsigned int count_components(const vector< shared_ptr< FunctionSpace > > &fs_vector)
Definition: fe_system.cc:27
unsigned int check_spacedim(const vector< shared_ptr< FunctionSpace > > &fs_vector)
Definition: fe_system.cc:37
Class FESystem for compound finite elements.
FEType
@ FEScalar
@ FEMixedSystem
@ FEVectorPiola
@ FETensor
@ FEVectorContravariant
@ FEVector
Global macros to enhance readability and debugging, general constants.
unsigned int uint
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