Flow123d  master-f44eb46
assembly_hm.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 assembly_hm.hh
15  * @brief
16  */
17 
18 #ifndef ASSEMBLY_HM_HH_
19 #define ASSEMBLY_HM_HH_
20 
23 #include "coupling/hm_iterative.hh"
24 #include "fem/fe_p.hh"
25 #include "fem/fe_values.hh"
27 #include "coupling/balance.hh"
29 
30 
31 /**
32  * Auxiliary container class for Finite element and related objects of given dimension.
33  */
34 template <unsigned int dim>
36 {
37 public:
38  typedef typename HM_Iterative::EqFields EqFields;
39  typedef typename HM_Iterative::EqData EqData;
40  typedef typename DarcyLMH::EqFields FlowEqFields;
41 
42  static constexpr const char * name() { return "FlowPotentialAssemblyHM"; }
43 
44  /// Constructor.
45  FlowPotentialAssemblyHM(EqFields *eq_fields, EqData *eq_data)
46  : AssemblyBase<dim>(1), eq_fields_(eq_fields), eq_data_(eq_data) {
48  this->used_fields_ += eq_fields_->alpha;
51  this->used_fields_ += eq_data_->flow_->eq_fields().bc_type;
52  this->used_fields_ += eq_data_->flow_->eq_fields().bc_pressure;
53  }
54 
55  /// Destructor.
57 
58  /// Initialize auxiliary vectors and other data members
59  void initialize(ElementCacheMap *element_cache_map) {
60  this->element_cache_map_ = element_cache_map;
61 
62  shared_ptr<FE_P<dim>> fe_p = std::make_shared< FE_P<dim> >(1);
63  shared_ptr<FiniteElement<dim>> fe_ = std::make_shared<FESystem<dim>>(fe_p, FEVector, 3);
65 
66  dof_indices_.resize(fe_->n_dofs());
67 
69  }
70 
71 
72  /// Assemble integral over element
73  inline void boundary_side_integral(DHCellSide dh_side)
74  {
76 
77  double ref_pot = 0;
78  if (dh_side.side().is_boundary())
79  {
80  auto p_side = *this->boundary_points(dh_side).begin();
81  auto p_bdr = p_side.point_bdr( dh_side.cond().element_accessor() );
82  unsigned int flow_bc_type = eq_data_->flow_->eq_fields().bc_type(p_bdr);
83  if (flow_bc_type == DarcyLMH::EqFields::dirichlet || flow_bc_type == DarcyLMH::EqFields::total_flux)
84  {
85  unsigned int k=0;
86  for ( auto p : this->boundary_points(dh_side) )
87  {
88  // The reference potential is applied only on dirichlet and total_flux b.c.,
89  // i.e. where only mechanical traction is prescribed.
90  auto p_bdr = p.point_bdr(dh_side.cond().element_accessor());
91  double alpha = eq_fields_->alpha(p);
92  double density = eq_fields_->density(p);
93  double gravity = eq_fields_->gravity(p);
94  double bc_pressure = eq_data_->flow_->eq_fields().bc_pressure(p_bdr);
95  ref_pot += -alpha*density*gravity*bc_pressure * fe_values_side_.JxW(k) / dh_side.measure();
96  ++k;
97  }
98  }
99  }
100  ref_potential_vec_.set(dof_indices_[dh_side.side_idx()], ref_pot);
101  }
102 
103 
104 
105 private:
106 
107  EqFields *eq_fields_; ///< Fields shared with HM_Iterative
108  EqData *eq_data_; ///< Data objects shared with HM_Iterative
109 
110  /// Sub field set contains fields used in calculation.
112 
113  FEValues<3> fe_values_side_; ///< FEValues of side object
114  LocDofVec dof_indices_; ///< Vector of global DOF indices
115  VectorMPI ref_potential_vec_; ///< Vector of dofs of field ref_potential
116 
117  template < template<IntDim...> class DimAssembly>
118  friend class GenericAssembly;
119 
120 };
121 
122 
123 template <unsigned int dim>
124 class ResidualAssemblyHM : public AssemblyBase<dim>
125 {
126 public:
128  typedef typename HM_Iterative::EqData EqData;
129 
130  static constexpr const char * name() { return "ResidualAssemblyHM"; }
131 
132  /// Constructor.
133  ResidualAssemblyHM(EqFields *eq_fields, EqData *eq_data)
134  : AssemblyBase<dim>(1), eq_fields_(eq_fields), eq_data_(eq_data) {
136  this->used_fields_ += eq_data_->flow_->eq_fields().field_ele_pressure;
138  }
139 
140  /// Destructor.
142 
143  /// Initialize auxiliary vectors and other data members
144  void initialize(ElementCacheMap *element_cache_map) {
145  this->element_cache_map_ = element_cache_map;
146  shared_ptr<FE_P<dim>> fe_ = std::make_shared< FE_P<dim> >(0);
148  }
149 
150 
151  /// Assemble integral over element
152  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
153  {
154  if (cell.dim() != dim) return;
155  if (!cell.is_own()) return;
156 
157  fe_values_.reinit(cell.elm());
158 
159  // compute pressure error
160  unsigned int k=0;
161  for (auto p : this->bulk_points(element_patch_idx) )
162  {
163  double new_p = eq_data_->flow_->eq_fields().field_ele_pressure(p);
164  double old_p = eq_fields_->old_iter_pressure(p);
165  eq_data_->p_dif2 += (new_p - old_p)*(new_p - old_p) * fe_values_.JxW(k);
166  eq_data_->p_norm2 += old_p*old_p * fe_values_.JxW(k);
167  ++k;
168  }
169  }
170 
171 
172 
173 private:
174  /// Data objects shared with Elasticity
177 
178  /// Sub field set contains fields used in calculation.
180 
181  FEValues<3> fe_values_; ///< FEValues of cell object
182 
183  template < template<IntDim...> class DimAssembly>
184  friend class GenericAssembly;
185 
186 };
187 
188 
189 #endif /* ASSEMBLY_HM_HH_ */
190 
Range< BulkPoint > bulk_points(unsigned int element_patch_idx) const
Return BulkPoint range of appropriate dimension.
Quadrature * quad_
Quadrature used in assembling methods.
Quadrature * quad_low_
Quadrature used in assembling methods (dim-1).
int active_integrals_
Holds mask of active integrals.
Range< BoundaryPoint > boundary_points(const DHCellSide &cell_side) const
Return BoundaryPoint range of appropriate dimension.
ElementCacheMap * element_cache_map_
ElementCacheMap shared with GenericAssembly object.
ElementAccessor< 3 > element_accessor()
Cell accessor allow iterate over DOF handler cells.
const ElementAccessor< 3 > elm() const
Return ElementAccessor to element of loc_ele_idx_.
bool is_own() const
Return true if accessor represents own element (false for ghost element)
LocDofVec get_loc_dof_indices() const
Returns the local indices of dofs associated to the cell on the local process.
unsigned int dim() const
Return dimension of element appropriate to cell.
Side accessor allows to iterate over sides of DOF handler cell.
Side side() const
Return Side of given cell and side_idx.
Boundary cond() const
const DHCellAccessor & cell() const
Return DHCellAccessor appropriate to the side.
double measure() const
unsigned int side_idx() const
Directing class of FieldValueCache.
double JxW(const unsigned int point_no)
Return the product of Jacobian determinant and the quadrature weight at given quadrature point.
Definition: fe_values.hh:223
void initialize(Quadrature &_quadrature, FiniteElement< DIM > &_fe, UpdateFlags _flags)
Initialize structures and calculates cell-independent data.
Definition: fe_values.cc:137
void reinit(const ElementAccessor< spacedim > &cell)
Update cell-dependent data (gradients, Jacobians etc.)
Definition: fe_values.cc:537
Container for various descendants of FieldCommonBase.
Definition: field_set.hh:159
FEValues< 3 > fe_values_side_
FEValues of side object.
Definition: assembly_hm.hh:113
HM_Iterative::EqFields EqFields
Definition: assembly_hm.hh:38
~FlowPotentialAssemblyHM()
Destructor.
Definition: assembly_hm.hh:56
EqData * eq_data_
Data objects shared with HM_Iterative.
Definition: assembly_hm.hh:108
FieldSet used_fields_
Sub field set contains fields used in calculation.
Definition: assembly_hm.hh:111
VectorMPI ref_potential_vec_
Vector of dofs of field ref_potential.
Definition: assembly_hm.hh:115
EqFields * eq_fields_
Fields shared with HM_Iterative.
Definition: assembly_hm.hh:107
LocDofVec dof_indices_
Vector of global DOF indices.
Definition: assembly_hm.hh:114
static constexpr const char * name()
Definition: assembly_hm.hh:42
HM_Iterative::EqData EqData
Definition: assembly_hm.hh:39
void boundary_side_integral(DHCellSide dh_side)
Assemble integral over element.
Definition: assembly_hm.hh:73
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
Definition: assembly_hm.hh:59
FlowPotentialAssemblyHM(EqFields *eq_fields, EqData *eq_data)
Constructor.
Definition: assembly_hm.hh:45
DarcyLMH::EqFields FlowEqFields
Definition: assembly_hm.hh:40
Generic class of assemblation.
double p_dif2
Squared norm of pressure difference in two subsequent iterations.
std::shared_ptr< DarcyLMH > flow_
steady or unsteady water flow simulator based on MH scheme
double p_norm2
Squared pressure norm in the last iteration.
Field< 3, FieldValue< 3 >::Scalar > density
Density of fluid.
Field< 3, FieldValue< 3 >::Scalar > old_iter_pressure
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > ref_potential_ptr_
FieldFE for pressure_potential field.
Field< 3, FieldValue< 3 >::Scalar > alpha
Biot coefficient.
Field< 3, FieldValue< 3 >::Scalar > gravity
Standard gravity.
~ResidualAssemblyHM()
Destructor.
Definition: assembly_hm.hh:141
ResidualAssemblyHM(EqFields *eq_fields, EqData *eq_data)
Constructor.
Definition: assembly_hm.hh:133
HM_Iterative::EqFields EqFields
Definition: assembly_hm.hh:127
HM_Iterative::EqData EqData
Definition: assembly_hm.hh:128
FieldSet used_fields_
Sub field set contains fields used in calculation.
Definition: assembly_hm.hh:179
EqFields * eq_fields_
Data objects shared with Elasticity.
Definition: assembly_hm.hh:175
static constexpr const char * name()
Definition: assembly_hm.hh:130
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
Definition: assembly_hm.hh:144
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
Definition: assembly_hm.hh:152
FEValues< 3 > fe_values_
FEValues of cell object.
Definition: assembly_hm.hh:181
bool is_boundary() const
Returns true for side on the boundary.
void set(unsigned int pos, double val)
Set value on given position.
Definition: vector_mpi.hh:111
Definitions of basic Lagrangean finite elements with polynomial shape functions.
Class FEValues calculates finite element data on the actual cells such as shape function values,...
@ FEVector
@ boundary
@ bulk
arma::Col< IntIdx > LocDofVec
Definition: index_types.hh:28
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19
Definitions of particular quadrature rules on simplices.
@ update_JxW_values
Transformed quadrature weights.
@ update_side_JxW_values
Transformed quadrature weight for cell sides.