Flow123d  master-f44eb46
field_fe.hh
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 field_fe.hh
15  * @brief
16  */
17 
18 #ifndef FIELD_FE_HH_
19 #define FIELD_FE_HH_
20 
21 #include "petscmat.h"
22 #include "system/system.hh"
23 #include "system/index_types.hh"
25 #include "fields/field.hh"
26 #include "la/vector_mpi.hh"
27 #include "mesh/mesh.h"
28 #include "mesh/point.hh"
29 #include "mesh/bih_tree.hh"
30 #include "mesh/range_wrapper.hh"
31 #include "io/element_data_cache.hh"
32 #include "io/msh_basereader.hh"
33 #include "fem/fe_p.hh"
34 #include "fem/fe_system.hh"
35 #include "fem/dofhandler.hh"
36 #include "fem/finite_element.hh"
37 #include "fem/dh_cell_accessor.hh"
38 #include "fem/mapping_p1.hh"
39 #include "input/factory.hh"
40 
41 #include <memory>
42 
43 
44 
45 /**
46  * Class representing fields given by finite element approximation.
47  *
48  */
49 template <int spacedim, class Value>
50 class FieldFE : public FieldAlgorithmBase<spacedim, Value>
51 {
52 public:
56 
57  /**
58  * Possible interpolations of input data.
59  */
61  {
62  identic_msh, //!< identical mesh
63  equivalent_msh, //!< equivalent mesh (default value)
64  gauss_p0, //!< P0 interpolation (with the use of Gaussian distribution)
65  interp_p0 //!< P0 interpolation (with the use of calculation of intersections)
66  };
67 
68  /// Declaration of exception.
69  TYPEDEF_ERR_INFO( EI_Field, std::string);
70  TYPEDEF_ERR_INFO( EI_File, std::string);
71  TYPEDEF_ERR_INFO( EI_ElemIdx, unsigned int);
72  TYPEDEF_ERR_INFO( EI_Region, std::string);
73  DECLARE_EXCEPTION( ExcInvalidElemeDim, << "Dimension of element in target mesh must be 0, 1 or 2! elm.idx() = " << EI_ElemIdx::val << ".\n" );
74  DECLARE_INPUT_EXCEPTION( ExcUndefElementValue,
75  << "FieldFE " << EI_Field::qval << " on region " << EI_Region::qval << " have invalid value .\n"
76  << "Provided by file " << EI_File::qval << " at element ID " << EI_ElemIdx::val << ".\n"
77  << "Please specify in default_value key.\n");
78 
79 
80 
81  static const unsigned int undef_uint = -1;
82 
83  /**
84  * Default constructor, optionally we need number of components @p n_comp in the case of Vector valued fields.
85  */
86  FieldFE(unsigned int n_comp=0);
87 
88  /**
89  * Return Record for initialization of FieldFE that is derived from Abstract
90  */
91  static const Input::Type::Record & get_input_type();
92 
93  /**
94  * Return Input selection for discretization type (determines the section of VTK file).
95  */
97 
98  /**
99  * Return Input selection that allow to set interpolation of input data.
100  */
102 
103  /**
104  * Factory class (descendant of @p Field<...>::FactoryBase) that is necessary
105  * for setting pressure values are piezometric head values.
106  */
108  public:
109  /// Constructor.
110  NativeFactory(unsigned int index, std::shared_ptr<DOFHandlerMultiDim> conc_dof_handler, VectorMPI dof_vector = VectorMPI::sequential(0))
111  : index_(index),
112  conc_dof_handler_(conc_dof_handler),
113  dof_vector_(dof_vector)
114  {}
115 
116  typename Field<spacedim,Value>::FieldBasePtr create_field(Input::Record rec, const FieldCommon &field) override;
117 
118  unsigned int index_;
119  std::shared_ptr<DOFHandlerMultiDim> conc_dof_handler_;
121  };
122 
123 
124  /**
125  * Setter for the finite element data.
126  * @param dh Dof handler.
127  * @param dof_values Vector of dof values, optional (create own vector according to dofhandler).
128  * @param block_index Index of block (in FESystem) or '-1' for FieldFE under simple DOF handler.
129  * @return Data vector of dof values.
130  */
131  VectorMPI set_fe_data(std::shared_ptr<DOFHandlerMultiDim> dh, VectorMPI dof_values = VectorMPI::sequential(0),
132  unsigned int block_index = FieldFE<spacedim, Value>::undef_uint);
133 
134  /**
135  * Overload @p FieldAlgorithmBase::cache_update
136  */
138  ElementCacheMap &cache_map, unsigned int region_patch_idx) override;
139 
140  /**
141  * Overload @p FieldAlgorithmBase::cache_reinit
142  *
143  * Reinit fe_values_ data member.
144  */
145  void cache_reinit(const ElementCacheMap &cache_map) override;
146 
147  /**
148  * Initialization from the input interface.
149  */
150  virtual void init_from_input(const Input::Record &rec, const struct FieldAlgoBaseInitData& init_data);
151 
152 
153  /**
154  * Update time and possibly update data from GMSH file.
155  */
156  bool set_time(const TimeStep &time) override;
157 
158 
159  /**
160  * Set target mesh.
161  */
162  void set_mesh(const Mesh *mesh) override;
163 
164 
165  /**
166  * Copy data vector to given output ElementDataCache
167  */
168  void native_data_to_cache(ElementDataCache<double> &output_data_cache);
169 
170 
171  /**
172  * Return size of vector of data stored in Field
173  */
174  unsigned int data_size() const;
175 
176  inline std::shared_ptr<DOFHandlerMultiDim> get_dofhandler() const {
177  return dh_;
178  }
179 
180  inline VectorMPI& vec() {
181  return data_vec_;
182  }
183 
184  inline const VectorMPI& vec() const {
185  return data_vec_;
186  }
187 
188  /// Call begin scatter functions (local to ghost) on data vector
190 
191  /// Call end scatter functions (local to ghost) on data vector
193 
194  /// Destructor.
195  virtual ~FieldFE();
196 
197 private:
198  /**
199  * Helper class holds specific data of field evaluation.
200  */
201  class FEItem {
202  public:
203  unsigned int comp_index_;
204  unsigned int range_begin_;
205  unsigned int range_end_;
206  };
207 
208  /**
209  * Helper class holds data of invalid values of all regions.
210  *
211  * If region contains invalid element value (typically 'not a number') is_invalid_ flag is set to true
212  * and other information (element id an value) are stored. Check of invalid values is performed
213  * during processing data of reader cache and possible exception is thrown only if FieldFE is defined
214  * on appropriate region.
215  *
216  * Data of all regions are stored in vector of RegionValueErr instances.
217  */
219  public:
220  /// Default constructor, sets valid region
222 
223  /// Constructor, sets invalid region, element and value specification
224  RegionValueErr(const std::string &region_name, unsigned int elm_id, double value) {
225  is_invalid_ = true;
226  region_name_ = region_name;
227  elm_id_ = elm_id;
228  value_ = value;
229  }
230 
232  std::string region_name_;
233  unsigned int elm_id_;
234  double value_;
235  };
236 
237  /// Create DofHandler object
238  void make_dof_handler(const MeshBase *mesh);
239 
240  /// Interpolate data (use Gaussian distribution) over all elements of target mesh.
241  void interpolate_gauss();
242 
243  /// Interpolate data (use intersection library) over all elements of target mesh.
245 
246 // /// Calculate native data over all elements of target mesh.
247 // void calculate_native_values(ElementDataCache<double>::CacheData data_cache);
248 //
249 // /// Calculate data of equivalent_mesh interpolation on input over all elements of target mesh.
250 // void calculate_equivalent_values(ElementDataCache<double>::CacheData data_cache);
251 
252  /// Calculate data of equivalent_mesh interpolation or native data on input over all elements of target mesh.
254 
255  /// Initialize FEValues object of given dimension.
256  template <unsigned int dim>
257  Quadrature init_quad(std::shared_ptr<EvalPoints> eval_points);
258 
260  unsigned int i_dof, unsigned int i_qp)
261  {
263  for (unsigned int c=0; c<Value::NRows_*Value::NCols_; ++c)
264  v(c/spacedim,c%spacedim) = fe_values_[dim].shape_value_component(i_dof, i_qp, c);
265  if (Value::NRows_ == Value::NCols_)
266  return v;
267  else
268  return v.t();
269  }
270 
271  template<unsigned int dim>
272  void fill_fe_system_data(unsigned int block_index) {
273  auto fe_system_ptr = std::dynamic_pointer_cast<FESystem<dim>>( dh_->ds()->fe()[Dim<dim>{}] );
274  ASSERT(fe_system_ptr != nullptr).error("Wrong type, must be FESystem!\n");
275  this->fe_item_[dim].comp_index_ = fe_system_ptr->function_space()->dof_indices()[block_index].component_offset;
276  this->fe_item_[dim].range_begin_ = fe_system_ptr->fe_dofs(block_index)[0];
277  this->fe_item_[dim].range_end_ = this->fe_item_[dim].range_begin_ + fe_system_ptr->fe()[block_index]->n_dofs();
278  }
279 
280 
281  template<unsigned int dim>
282  void fill_fe_item() {
283  this->fe_item_[dim].comp_index_ = 0;
284  this->fe_item_[dim].range_begin_ = 0;
285  this->fe_item_[dim].range_end_ = dh_->ds()->fe()[Dim<dim>{}]->n_dofs();
286  }
287 
288  /**
289  * Method computes value of given input cache element.
290  *
291  * If computed value is invalid (e.g. NaN value) sets the data specifying error value.
292  * @param i_cache_el Index of element of input ElementDataCache
293  * @param elm_idx Idx of element of computational mesh
294  * @param region_name Region of computational mesh
295  * @param actual_compute_region_error Data object holding data of region with invalid value.
296  */
297  double get_scaled_value(int i_cache_el, unsigned int elm_idx, const std::string &region_name, RegionValueErr &actual_compute_region_error);
298 
299  /**
300  * Helper method. Compute real coordinates and weights (use QGauss) of given element.
301  *
302  * Method is needs in Gauss interpolation.
303  */
304  template<int elemdim>
305  unsigned int compute_fe_quadrature(std::vector<arma::vec::fixed<3>> & q_points, std::vector<double> & q_weights,
306  const ElementAccessor<spacedim> &elm, unsigned int order=3)
307  {
308  static_assert(elemdim <= spacedim, "Dimension of element must be less equal than spacedim.");
309  static const double weight_coefs[] = { 1., 1., 2., 6. };
310 
311  QGauss qgauss(elemdim, order);
313 
314  for(unsigned i=0; i<qgauss.size(); ++i) {
315  q_weights[i] = qgauss.weight(i)*weight_coefs[elemdim];
317  }
318 
319  return qgauss.size();
320  }
321 
322 
323  /// DOF handler object
324  std::shared_ptr<DOFHandlerMultiDim> dh_;
325  /// Store data of Field
327 
328  /// mesh reader file
330 
331  /// field name read from input
332  std::string field_name_;
333 
334  /// Specify section where to find the field data in input mesh file
336 
337  /// Specify type of FE data interpolation
339 
340  /// Field flags.
342 
343  /// Default value of element if not set in mesh data file
345 
346  /// Accessor to Input::Record
348 
349  /// Is set in set_mesh method. Value true means, that we accept only boundary element accessors in the @p value method.
351 
352  /// List of FEValues objects of dimensions 0,1,2,3 used for value calculation
354 
355  /// Maps element indices from computational mesh to the source (data).
356  std::shared_ptr<EquivalentMeshMap> source_target_mesh_elm_map_;
357 
358  /// Holds specific data of field evaluation over all dimensions.
359  std::array<FEItem, 4> fe_item_;
361 
362  /// Set holds data of valid / invalid element values on all regions
364 
365  /// Input ElementDataCache is stored in set_time and used in all evaluation and interpolation methods.
367 
368  /// Registrar of class to factory
369  static const int registrar;
370 };
371 
372 
373 /** Create FieldFE from dof handler */
374 template <int spacedim, class Value>
375 std::shared_ptr<FieldFE<spacedim, Value> > create_field_fe(std::shared_ptr<DOFHandlerMultiDim> dh,
376  VectorMPI *vec = nullptr,
377  unsigned int block_index = FieldFE<spacedim, Value>::undef_uint)
378 {
379  // Construct FieldFE
380  std::shared_ptr< FieldFE<spacedim, Value> > field_ptr = std::make_shared< FieldFE<spacedim, Value> >();
381  if (vec == nullptr)
382  field_ptr->set_fe_data( dh, dh->create_vector(), block_index );
383  else
384  field_ptr->set_fe_data( dh, *vec, block_index );
385 
386  return field_ptr;
387 }
388 
389 
390 /** Create FieldFE with parallel VectorMPI from finite element */
391 template <int spacedim, class Value>
392 std::shared_ptr<FieldFE<spacedim, Value> > create_field_fe(Mesh & mesh, const MixedPtr<FiniteElement> &fe)
393 {
394  // Prepare DOF handler
395  std::shared_ptr<DOFHandlerMultiDim> dh_par = std::make_shared<DOFHandlerMultiDim>(mesh);
396  std::shared_ptr<DiscreteSpace> ds = std::make_shared<EqualOrderDiscreteSpace>( &mesh, fe);
397  dh_par->distribute_dofs(ds);
398 
399  return create_field_fe<spacedim,Value>(dh_par);
400 }
401 
402 
403 #endif /* FIELD_FE_HH_ */
#define ASSERT(expr)
Definition: asserts.hh:351
Directing class of FieldValueCache.
unsigned int n_comp() const
Space< spacedim >::Point Point
Common abstract parent of all Field<...> classes.
Definition: field_common.hh:77
unsigned int range_end_
Definition: field_fe.hh:205
unsigned int range_begin_
Definition: field_fe.hh:204
unsigned int comp_index_
Definition: field_fe.hh:203
Field< spacedim, Value >::FieldBasePtr create_field(Input::Record rec, const FieldCommon &field) override
Definition: field_fe.cc:144
std::shared_ptr< DOFHandlerMultiDim > conc_dof_handler_
Definition: field_fe.hh:119
NativeFactory(unsigned int index, std::shared_ptr< DOFHandlerMultiDim > conc_dof_handler, VectorMPI dof_vector=VectorMPI::sequential(0))
Constructor.
Definition: field_fe.hh:110
unsigned int index_
Definition: field_fe.hh:118
RegionValueErr(const std::string &region_name, unsigned int elm_id, double value)
Constructor, sets invalid region, element and value specification.
Definition: field_fe.hh:224
std::string region_name_
Definition: field_fe.hh:232
unsigned int elm_id_
Definition: field_fe.hh:233
RegionValueErr()
Default constructor, sets valid region.
Definition: field_fe.hh:221
FieldAlgorithmBase< spacedim, Value >::Point Point
Definition: field_fe.hh:53
DECLARE_EXCEPTION(ExcInvalidElemeDim,<< "Dimension of element in target mesh must be 0, 1 or 2! elm.idx() = "<< EI_ElemIdx::val<< ".\n")
void cache_reinit(const ElementCacheMap &cache_map) override
Definition: field_fe.cc:262
Armor::ArmaMat< typename Value::element_type, Value::NRows_, Value::NCols_ > handle_fe_shape(unsigned int dim, unsigned int i_dof, unsigned int i_qp)
Definition: field_fe.hh:259
double get_scaled_value(int i_cache_el, unsigned int elm_idx, const std::string &region_name, RegionValueErr &actual_compute_region_error)
Definition: field_fe.cc:797
MixedPtr< FiniteElement > fe_
Definition: field_fe.hh:360
FieldAlgorithmBase< spacedim, Value > FactoryBaseType
Definition: field_fe.hh:54
VectorMPI data_vec_
Store data of Field.
Definition: field_fe.hh:326
static const Input::Type::Selection & get_interp_selection_input_type()
Definition: field_fe.cc:105
unsigned int data_size() const
Definition: field_fe.cc:776
std::vector< RegionValueErr > region_value_err_
Set holds data of valid / invalid element values on all regions.
Definition: field_fe.hh:363
void local_to_ghost_data_scatter_begin()
Call begin scatter functions (local to ghost) on data vector.
Definition: field_fe.cc:783
TYPEDEF_ERR_INFO(EI_File, std::string)
void calculate_element_values()
Calculate data of equivalent_mesh interpolation or native data on input over all elements of target m...
Definition: field_fe.cc:635
virtual void init_from_input(const Input::Record &rec, const struct FieldAlgoBaseInitData &init_data)
Definition: field_fe.cc:285
virtual ~FieldFE()
Destructor.
Definition: field_fe.cc:833
Field< spacedim, Value >::FactoryBase FieldFactoryBaseType
Definition: field_fe.hh:55
DataInterpolation
Definition: field_fe.hh:61
@ interp_p0
P0 interpolation (with the use of calculation of intersections)
Definition: field_fe.hh:65
@ equivalent_msh
equivalent mesh (default value)
Definition: field_fe.hh:63
@ gauss_p0
P0 interpolation (with the use of Gaussian distribution)
Definition: field_fe.hh:64
@ identic_msh
identical mesh
Definition: field_fe.hh:62
void local_to_ghost_data_scatter_end()
Call end scatter functions (local to ghost) on data vector.
Definition: field_fe.cc:790
static const Input::Type::Record & get_input_type()
Implementation.
Definition: field_fe.cc:60
unsigned int compute_fe_quadrature(std::vector< arma::vec::fixed< 3 >> &q_points, std::vector< double > &q_weights, const ElementAccessor< spacedim > &elm, unsigned int order=3)
Definition: field_fe.hh:305
std::vector< FEValues< spacedim > > fe_values_
List of FEValues objects of dimensions 0,1,2,3 used for value calculation.
Definition: field_fe.hh:353
void interpolate_gauss()
Interpolate data (use Gaussian distribution) over all elements of target mesh.
Definition: field_fe.cc:436
TYPEDEF_ERR_INFO(EI_Field, std::string)
Declaration of exception.
void make_dof_handler(const MeshBase *mesh)
Create DofHandler object.
Definition: field_fe.cc:345
Quadrature init_quad(std::shared_ptr< EvalPoints > eval_points)
Initialize FEValues object of given dimension.
Definition: field_fe.cc:275
OutputTime::DiscreteSpace discretization_
Specify section where to find the field data in input mesh file.
Definition: field_fe.hh:335
std::shared_ptr< DOFHandlerMultiDim > get_dofhandler() const
Definition: field_fe.hh:176
TYPEDEF_ERR_INFO(EI_Region, std::string)
FilePath reader_file_
mesh reader file
Definition: field_fe.hh:329
VectorMPI & vec()
Definition: field_fe.hh:180
std::array< FEItem, 4 > fe_item_
Holds specific data of field evaluation over all dimensions.
Definition: field_fe.hh:359
TYPEDEF_ERR_INFO(EI_ElemIdx, unsigned int)
bool set_time(const TimeStep &time) override
Definition: field_fe.cc:385
static const int registrar
Registrar of class to factory.
Definition: field_fe.hh:369
void interpolate_intersection()
Interpolate data (use intersection library) over all elements of target mesh.
Definition: field_fe.cc:532
VectorMPI set_fe_data(std::shared_ptr< DOFHandlerMultiDim > dh, VectorMPI dof_values=VectorMPI::sequential(0), unsigned int block_index=FieldFE< spacedim, Value >::undef_uint)
Definition: field_fe.cc:175
double default_value_
Default value of element if not set in mesh data file.
Definition: field_fe.hh:344
DECLARE_INPUT_EXCEPTION(ExcUndefElementValue,<< "FieldFE "<< EI_Field::qval<< " on region "<< EI_Region::qval<< " have invalid value .\n"<< "Provided by file "<< EI_File::qval<< " at element ID "<< EI_ElemIdx::val<< ".\n"<< "Please specify in default_value key.\n")
const VectorMPI & vec() const
Definition: field_fe.hh:184
void cache_update(FieldValueCache< typename Value::element_type > &data_cache, ElementCacheMap &cache_map, unsigned int region_patch_idx) override
Definition: field_fe.cc:215
bool boundary_domain_
Is set in set_mesh method. Value true means, that we accept only boundary element accessors in the va...
Definition: field_fe.hh:350
DataInterpolation interpolation_
Specify type of FE data interpolation.
Definition: field_fe.hh:338
Input::Record in_rec_
Accessor to Input::Record.
Definition: field_fe.hh:347
FieldFlag::Flags flags_
Field flags.
Definition: field_fe.hh:341
void native_data_to_cache(ElementDataCache< double > &output_data_cache)
Definition: field_fe.cc:757
std::string field_name_
field name read from input
Definition: field_fe.hh:332
FieldFE(unsigned int n_comp=0)
Definition: field_fe.cc:134
std::shared_ptr< DOFHandlerMultiDim > dh_
DOF handler object.
Definition: field_fe.hh:324
void fill_fe_item()
Definition: field_fe.hh:282
std::shared_ptr< EquivalentMeshMap > source_target_mesh_elm_map_
Maps element indices from computational mesh to the source (data).
Definition: field_fe.hh:356
void fill_fe_system_data(unsigned int block_index)
Definition: field_fe.hh:272
static const unsigned int undef_uint
Definition: field_fe.hh:81
static const Input::Type::Selection & get_disc_selection_input_type()
Definition: field_fe.cc:93
ElementDataCache< double >::CacheData input_data_cache_
Input ElementDataCache is stored in set_time and used in all evaluation and interpolation methods.
Definition: field_fe.hh:366
void set_mesh(const Mesh *mesh) override
Definition: field_fe.cc:309
std::shared_ptr< FieldBaseType > FieldBasePtr
Definition: field.hh:96
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
Record type proxy class.
Definition: type_record.hh:182
Template for classes storing finite set of named values.
static RealPoint project_unit_to_real(const BaryPoint &point, const ElementMap &map)
Definition: mapping_p1.cc:85
static ElementMap element_map(ElementAccessor< 3 > elm)
Definition: mapping_p1.cc:48
Base class for Mesh and BCMesh.
Definition: mesh.h:96
Definition: mesh.h:362
Symmetric Gauss-Legendre quadrature formulae on simplices.
Base class for quadrature rules on simplices in arbitrary dimensions.
Definition: quadrature.hh:48
unsigned int size() const
Returns number of quadrature points.
Definition: quadrature.hh:86
Armor::ArmaVec< double, point_dim > point(unsigned int i) const
Returns the ith quadrature point.
Definition: quadrature.hh:151
double weight(unsigned int i) const
Returns the ith weight.
Definition: quadrature.hh:103
Representation of one time step..
static VectorMPI sequential(unsigned int size)
Definition: vector_mpi.cc:44
Declaration of class which handles the ordering of degrees of freedom (dof) and mappings between loca...
Definitions of basic Lagrangean finite elements with polynomial shape functions.
Class FESystem for compound finite elements.
std::shared_ptr< FieldFE< spacedim, Value > > create_field_fe(std::shared_ptr< DOFHandlerMultiDim > dh, VectorMPI *vec=nullptr, unsigned int block_index=FieldFE< spacedim, Value >::undef_uint)
Definition: field_fe.hh:375
Abstract class for description of finite elements.
static constexpr bool value
Definition: json.hpp:87
Class MappingP1 implements the affine transformation of the unit cell onto the actual cell.
typename arma::Mat< Type >::template fixed< nr, nc > ArmaMat
Definition: armor.hh:502
ArmaMat< double, N, M > mat
Definition: armor.hh:888
ArmaVec< double, N > vec
Definition: armor.hh:885
Implementation of range helper class.
Definition: mixed.hh:25
Helper struct stores data for initizalize descentants of FieldAlgorithmBase.