Flow123d  master-f44eb46
assembly_richards.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_richards.hh
15  * @brief
16  */
17 
18 #ifndef ASSEMBLY_RICHARDS_HH_
19 #define ASSEMBLY_RICHARDS_HH_
20 
23 #include "flow/assembly_lmh.hh"
24 #include "flow/soil_models.hh"
26 
27 
28 template <unsigned int dim>
30 {
31 public:
32  typedef typename RichardsLMH::EqFields EqFields;
33  typedef typename RichardsLMH::EqData EqData;
34 
35  static constexpr const char * name() { return "InitCondPostprocessAssembly"; }
36 
37  /// Constructor.
39  : AssemblyBase<dim>(0), eq_fields_(eq_fields), eq_data_(eq_data) {
41  this->used_fields_ += this->eq_fields_->storativity;
47  this->used_fields_ += this->eq_fields_->conductivity;
48  this->used_fields_ += this->eq_fields_->cross_section;
49  }
50 
51  /// Destructor.
53 
54  /// Initialize auxiliary vectors and other data members
55  void initialize(ElementCacheMap *element_cache_map) {
56  //this->balance_ = eq_data_->balance_;
57  this->element_cache_map_ = element_cache_map;
58  }
59 
60  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx) {
61  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
62 
65  const DHCellAccessor dh_cell = cell.cell_with_other_dh(this->eq_data_->dh_.get());
66 
67  auto p = *( this->bulk_points(element_patch_idx).begin() );
68  bool genuchten_on = reset_soil_model(cell, p);
70  + this->eq_fields_->extra_storativity(p);
71  VectorMPI water_content_vec = this->eq_fields_->water_content_ptr->vec();
72  double diagonal_coef = cell.elm().measure() * eq_fields_->cross_section(p) / cell.elm()->n_sides();
73 
74  for (unsigned int i=0; i<cell.elm()->n_sides(); i++) {
75  const int local_side = cr_disc_dofs_[i];
76  capacity = 0;
77  water_content = 0;
79 
80  if (genuchten_on) {
81  fadbad::B<double> x_phead(phead);
82  fadbad::B<double> evaluated( this->eq_data_->soil_model_->water_content_diff(x_phead) );
83  evaluated.diff(0,1);
84  water_content = evaluated.val();
85  capacity = x_phead.d(0);
86  }
88  water_content_vec.set( cr_disc_dofs_[i], water_content + storativity_ * phead);
89 
90  this->eq_data_->balance_->add_mass_values(eq_data_->water_balance_idx, dh_cell, {local_side},
91  {0.0}, diagonal_coef*(water_content + storativity_ * phead) );
92  }
93  }
94 
95  /// Implements @p AssemblyBase::begin.
96  void begin() override
97  {
98  this->eq_data_->balance_->start_mass_assembly(this->eq_data_->water_balance_idx);
99  }
100 
101 
102  /// Implements @p AssemblyBase::end.
103  void end() override
104  {
105  this->eq_data_->balance_->finish_mass_assembly(this->eq_data_->water_balance_idx);
106  }
107 
108 
109 private:
110  bool reset_soil_model(const DHCellAccessor& cell, BulkPoint &p) {
111  bool genuchten_on = (this->eq_fields_->genuchten_p_head_scale.field_result({cell.elm().region()}) != result_zeros);
112  if (genuchten_on) {
113  SoilData soil_data;
114  soil_data.n = this->eq_fields_->genuchten_n_exponent(p);
115  soil_data.alpha = this->eq_fields_->genuchten_p_head_scale(p);
116  soil_data.Qr = this->eq_fields_->water_content_residual(p);
117  soil_data.Qs = this->eq_fields_->water_content_saturated(p);
118  soil_data.Ks = this->eq_fields_->conductivity(p);
119  //soil_data.cut_fraction = 0.99; // set by model
120 
121  this->eq_data_->soil_model_->reset(soil_data);
122  }
123  return genuchten_on;
124  }
125 
126 
127  /// Sub field set contains fields used in calculation.
129 
130  /// Data objects shared with ConvectionTransport
133 
134  LocDofVec cr_disc_dofs_; ///< Vector of local DOF indices pre-computed on different DOF handlers
135  LocDofVec edge_indices_; ///< Dofs of discontinuous fields on element edges.
136  double storativity_;
138 
139  template < template<IntDim...> class DimAssembly>
140  friend class GenericAssembly;
141 };
142 
143 
144 template <unsigned int dim>
146 {
147 public:
148  typedef typename RichardsLMH::EqFields EqFields;
149  typedef typename RichardsLMH::EqData EqData;
150 
151  static constexpr const char * name() { return "MHMatrixAssemblyRichards"; }
152 
154  : MHMatrixAssemblyLMH<dim>(eq_fields, eq_data), eq_fields_(eq_fields), eq_data_(eq_data) {
159  this->used_fields_ += eq_fields_->sigma;
168  this->used_fields_ += eq_fields_->bc_type;
170  this->used_fields_ += eq_fields_->bc_flux;
174  }
175 
176  /// Destructor.
178 
179  /// Initialize auxiliary vectors and other data members
180  void initialize(ElementCacheMap *element_cache_map) {
181  //this->balance_ = eq_data_->balance_;
182  MHMatrixAssemblyLMH<dim>::initialize(element_cache_map);
183  }
184 
185 
186  /// Integral over element.
187  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
188  {
189  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
190 
191  // evaluation point
192  auto p = *( this->bulk_points(element_patch_idx).begin() );
193  this->bulk_local_idx_ = cell.local_idx();
194 
195  this->asm_sides(cell, p, this->compute_conductivity(cell, p));
196  this->asm_element();
197  this->asm_source_term_richards(cell, p);
198  }
199 
200 
201  /// Assembles between boundary element and corresponding side on bulk element.
202  inline void boundary_side_integral(DHCellSide cell_side)
203  {
204  ASSERT_EQ(cell_side.dim(), dim).error("Dimension of element mismatch!");
205  if (!cell_side.cell().is_own()) return;
206 
207  auto p_side = *( this->boundary_points(cell_side).begin() );
208  auto p_bdr = p_side.point_bdr(cell_side.cond().element_accessor() );
209  ElementAccessor<3> b_ele = cell_side.side().cond().element_accessor(); // ??
210 
211  this->precompute_boundary_side(cell_side, p_side, p_bdr);
212 
213  if (this->type_==DarcyLMH::EqFields::seepage) {
214  this->use_dirichlet_switch(cell_side, b_ele, p_bdr);
215  }
216 
217  this->boundary_side_integral_in(cell_side, b_ele, p_bdr);
218  }
219 
220 
221  /// Implements @p AssemblyBase::begin.
222  void begin() override
223  {
224  this->begin_mh_matrix();
225  }
226 
227 
228  /// Implements @p AssemblyBase::end.
229  void end() override
230  {
231  this->end_mh_matrix();
232  }
233 
234 
235 protected:
236  /// Part of cell_integral method, specialized in Richards equation
237  inline void asm_source_term_richards(const DHCellAccessor& cell, BulkPoint &p)
238  {
239  update_water_content(cell, p);
240  const ElementAccessor<3> ele = cell.elm();
241 
242  // set lumped source
243  diagonal_coef_ = ele.measure() * eq_fields_->cross_section(p) / ele->n_sides();
245 
246  VectorMPI water_content_vec = eq_fields_->water_content_ptr->vec();
247 
248  const DHCellAccessor cr_cell = cell.cell_with_other_dh(eq_data_->dh_cr_.get());
249  auto loc_dof_vec = cr_cell.get_loc_dof_indices();
250 
251  for (unsigned int i=0; i<ele->n_sides(); i++)
252  {
253 
254  const int local_side = cr_disc_dofs_[i];
255  /*if (this->dirichlet_edge[i] == 0)*/ { //TODO Fix condition evaluating dirichlet_edge
256 
257  water_content_diff_ = -water_content_vec.get(local_side) + eq_data_->water_content_previous_time.get(local_side);
259 
260  /*
261  DebugOut().fmt("w diff: {:g} mass: {:g} w prev: {:g} w time: {:g} c: {:g} p: {:g} z: {:g}",
262  water_content_diff,
263  mass_diagonal * eq_data_->p_edge_solution[local_edge],
264  -eq_data_->water_content_previous_it[local_side],
265  eq_data_->water_content_previous_time[local_side],
266  capacity,
267  eq_data_->p_edge_solution[local_edge],
268  ele.centre()[0] );
269  */
270 
271 
274 
275  /*
276  DBGCOUT(<< "source [" << loc_system_.row_dofs[this->loc_edge_dofs[i]] << ", " << loc_system_.row_dofs[this->loc_edge_dofs[i]]
277  << "] mat: " << -mass_diagonal/this->eq_data_->time_step_
278  << " rhs: " << -source_diagonal_ - mass_rhs
279  << "\n");
280  */
281  eq_data_->loc_system_[cell.local_idx()].add_value(eq_data_->loc_edge_dofs[dim-1][i], eq_data_->loc_edge_dofs[dim-1][i],
284  }
285 
286  eq_data_->balance_->add_mass_values(eq_data_->water_balance_idx, cell, {local_side},
287  {0.0}, diagonal_coef_*water_content_vec.get(local_side));
288  eq_data_->balance_->add_source_values(eq_data_->water_balance_idx, ele.region().bulk_idx(),
289  {this->eq_data_->loc_system_[cell.local_idx()].row_dofs[eq_data_->loc_edge_dofs[dim-1][i]]},
290  {0},{source_diagonal_});
291  }
292  }
293 
294  bool reset_soil_model(const DHCellAccessor& cell, BulkPoint &p) {
295  bool genuchten_on = (this->eq_fields_->genuchten_p_head_scale.field_result({cell.elm().region()}) != result_zeros);
296  if (genuchten_on) {
297  SoilData soil_data;
298  soil_data.n = this->eq_fields_->genuchten_n_exponent(p);
299  soil_data.alpha = this->eq_fields_->genuchten_p_head_scale(p);
300  soil_data.Qr = this->eq_fields_->water_content_residual(p);
301  soil_data.Qs = this->eq_fields_->water_content_saturated(p);
302  soil_data.Ks = this->eq_fields_->conductivity(p);
303  //soil_data.cut_fraction = 0.99; // set by model
304 
305  this->eq_data_->soil_model_->reset(soil_data);
306  }
307  return genuchten_on;
308  }
309 
310 
312  edge_indices_ = cell.cell_with_other_dh(this->eq_data_->dh_cr_.get()).get_loc_dof_indices();
313  cr_disc_dofs_ = cell.cell_with_other_dh(this->eq_data_->dh_cr_disc_.get()).get_loc_dof_indices();
314 
315  bool genuchten_on = reset_soil_model(cell, p);
316  storativity_ = this->eq_fields_->storativity(p)
317  + this->eq_fields_->extra_storativity(p);
318  VectorMPI water_content_vec = this->eq_fields_->water_content_ptr->vec();
319 
320  for (unsigned int i=0; i<cell.elm()->n_sides(); i++) {
321  capacity = 0;
322  water_content = 0;
323  phead = this->eq_data_->p_edge_solution.get( edge_indices_[i] );
324 
325  if (genuchten_on) {
326  fadbad::B<double> x_phead(phead);
327  fadbad::B<double> evaluated( this->eq_data_->soil_model_->water_content_diff(x_phead) );
328  evaluated.diff(0,1);
329  water_content = evaluated.val();
330  capacity = x_phead.d(0);
331  }
332  this->eq_data_->capacity.set( cr_disc_dofs_[i], capacity + storativity_ );
333  water_content_vec.set( cr_disc_dofs_[i], water_content + storativity_ * phead);
334  }
335  }
336 
337  /// Precompute conductivity on bulk point.
339  {
340  bool genuchten_on = reset_soil_model(cell, p);
341 
342  double conductivity = 0;
343  if (genuchten_on) {
344  const DHCellAccessor cr_cell = cell.cell_with_other_dh(eq_data_->dh_cr_.get());
345  auto loc_dof_vec = cr_cell.get_loc_dof_indices();
346 
347  for (unsigned int i=0; i<cell.elm()->n_sides(); i++)
348  {
349  double phead = eq_data_->p_edge_solution.get( loc_dof_vec[i] );
350  conductivity += eq_data_->soil_model_->conductivity(phead);
351  }
352  conductivity /= cell.elm()->n_sides();
353  }
354  else {
355  conductivity = eq_fields_->conductivity(p);
356  }
357  return conductivity;
358  }
359 
360 
361  /// Postprocess velocity after calculating of cell integral.
363  {
364  this->postprocess_velocity(dh_cell, p);
365 
366  this->update_water_content(dh_cell, p);
367 
368  VectorMPI water_content_vec = eq_fields_->water_content_ptr->vec();
369 
370  for (unsigned int i=0; i<dh_cell.elm()->n_sides(); i++) {
371  water_content = water_content_vec.get( this->cr_disc_dofs_[i] );
372  water_content_previous_time = eq_data_->water_content_previous_time.get( this->cr_disc_dofs_[i] );
373 
374  solution[eq_data_->loc_side_dofs[dim-1][i]]
375  += this->edge_source_term_ - this->edge_scale_ * (water_content - water_content_previous_time) / eq_data_->time_step_;
376  }
377 
378  IntIdx p_dof = dh_cell.cell_with_other_dh(eq_data_->dh_p_.get()).get_loc_dof_indices()(0);
379  eq_fields_->conductivity_ptr->vec().set( p_dof, compute_conductivity(dh_cell, p) );
380  }
381 
382 
383  /// Data objects shared with ConvectionTransport
386 
387  LocDofVec cr_disc_dofs_; ///< Vector of local DOF indices pre-computed on different DOF handlers
388  LocDofVec edge_indices_; ///< Dofs of discontinuous fields on element edges.
389  double storativity_;
390  double capacity, phead;
391  double water_content, water_content_previous_time;
392  double diagonal_coef_, source_diagonal_;
393  double water_content_diff_, mass_diagonal_, mass_rhs_;
394 
395  template < template<IntDim...> class DimAssembly>
396  friend class GenericAssembly;
397 };
398 
399 
400 template <unsigned int dim>
402 {
403 public:
404  typedef typename RichardsLMH::EqFields EqFields;
405  typedef typename RichardsLMH::EqData EqData;
406 
407  static constexpr const char * name() { return "ReconstructSchurAssemblyRichards"; }
408 
410  : MHMatrixAssemblyRichards<dim>(eq_fields, eq_data) {
411  }
412 
413  /// Integral over element.
414  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
415  {
416  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
417 
418  // evaluation point
419  auto p = *( this->bulk_points(element_patch_idx).begin() );
420  this->bulk_local_idx_ = cell.local_idx();
421 
422  { // postprocess the velocity
423  this->eq_data_->postprocess_solution_[this->bulk_local_idx_].zeros(this->eq_data_->schur_offset_[dim-1]);
424  this->postprocess_velocity_richards(cell, p, this->eq_data_->postprocess_solution_[this->bulk_local_idx_]);
425  }
426  }
427 
428 
429  /// Assembles between boundary element and corresponding side on bulk element.
431  {}
432 
433  inline void dimjoin_intergral(FMT_UNUSED DHCellAccessor cell_lower_dim, FMT_UNUSED DHCellSide neighb_side)
434  {}
435 
436 
437  /// Implements @p AssemblyBase::begin.
438  void begin() override
439  {
440  this->begin_reconstruct_schur();
441  }
442 
443 
444  /// Implements @p AssemblyBase::end.
445  void end() override
446  {
447  this->end_reconstruct_schur();
448  }
449 protected:
450  template < template<IntDim...> class DimAssembly>
451  friend class GenericAssembly;
452 };
453 
454 
455 
456 #endif /* ASSEMBLY_RICHARDS_HH_ */
457 
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
Range< BulkPoint > bulk_points(unsigned int element_patch_idx) const
Return BulkPoint range of appropriate dimension.
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()
Base point accessor class.
Definition: eval_subset.hh:55
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.
DHCellAccessor cell_with_other_dh(const DOFHandlerMultiDim *dh) const
Create new accessor with same local idx and given DOF handler. Actual and given DOF handler must be c...
unsigned int dim() const
Return dimension of element appropriate to cell.
unsigned int local_idx() const
Return local index to element (index of DOF handler).
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.
unsigned int dim() const
Return dimension of element appropriate to the side.
std::array< unsigned int, 3 > schur_offset_
Index offset in the local system for the Schur complement (of dim = 1,2,3).
std::vector< arma::vec > postprocess_solution_
std::shared_ptr< SubDOFHandlerMultiDim > dh_cr_
DOF handler represents DOFs of edges.
std::shared_ptr< DOFHandlerMultiDim > dh_cr_disc_
DOF handler represents DOFs of sides.
VectorMPI p_edge_solution
std::shared_ptr< DOFHandlerMultiDim > dh_
full DOF handler represents DOFs of sides, elements and edges
std::vector< LocalSystem > loc_system_
Following data members are stored in vectors, one item for every cell.
std::array< std::vector< unsigned int >, 3 > loc_edge_dofs
std::shared_ptr< Balance > balance_
Shared Balance object.
Field< 3, FieldValue< 3 >::Scalar > water_source_density
Field< 3, FieldValue< 3 >::Scalar > extra_storativity
BCField< 3, FieldValue< 3 >::Enum > bc_type
Field< 3, FieldValue< 3 >::Scalar > sigma
Field< 3, FieldValue< 3 >::Scalar > storativity
BCField< 3, FieldValue< 3 >::Scalar > bc_pressure
BCField< 3, FieldValue< 3 >::Scalar > bc_flux
Field< 3, FieldValue< 3 >::TensorFixed > anisotropy
Field< 3, FieldValue< 3 >::Scalar > cross_section
Field< 3, FieldValue< 3 >::Scalar > conductivity
BCField< 3, FieldValue< 3 >::Scalar > bc_robin_sigma
Field< 3, FieldValue< 3 >::Scalar > extra_source
Externally added storativity.
BCField< 3, FieldValue< 3 >::Scalar > bc_switch_pressure
double measure() const
Computes the measure of the element.
Region region() const
Definition: accessors.hh:198
Directing class of FieldValueCache.
unsigned int n_sides() const
Definition: elements.h:131
Container for various descendants of FieldCommonBase.
Definition: field_set.hh:159
FieldResult field_result(RegionSet region_set) const override
Indicates special field states.
Definition: field.impl.hh:391
Generic class of assemblation.
bool reset_soil_model(const DHCellAccessor &cell, BulkPoint &p)
EqFields * eq_fields_
Data objects shared with ConvectionTransport.
static constexpr const char * name()
void end() override
Implements AssemblyBase::end.
InitCondPostprocessAssembly(EqFields *eq_fields, EqData *eq_data)
Constructor.
void begin() override
Implements AssemblyBase::begin.
RichardsLMH::EqFields EqFields
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
LocDofVec cr_disc_dofs_
Vector of local DOF indices pre-computed on different DOF handlers.
FieldSet used_fields_
Sub field set contains fields used in calculation.
LocDofVec edge_indices_
Dofs of discontinuous fields on element edges.
void end_mh_matrix()
Common code of end method of MH matrix assembly (Darcy and Richards)
FieldSet used_fields_
Sub field set contains fields used in calculation.
void use_dirichlet_switch(DHCellSide &cell_side, const ElementAccessor< 3 > &b_ele, BulkPoint &p_bdr)
Update BC switch dirichlet in MH matrix assembly if BC type is seepage.
void asm_sides(const DHCellAccessor &cell, BulkPoint &p, double conductivity)
Part of cell_integral method, common in all descendants.
void begin_mh_matrix()
Common code of begin method of MH matrix assembly (Darcy and Richards)
void boundary_side_integral_in(DHCellSide cell_side, const ElementAccessor< 3 > &b_ele, BulkPoint &p_bdr)
void precompute_boundary_side(DHCellSide &cell_side, BoundaryPoint &p_side, BulkPoint &p_bdr)
Precompute values used in boundary side integral on given DHCellSide.
void asm_element()
Part of cell_integral method, common in all descendants.
unsigned int bulk_local_idx_
Local idx of bulk element.
void end_reconstruct_schur()
Common code of end method of Reconstruct Schur assembly (Darcy and Richards)
DarcyLMH::EqFields::BC_Type type_
Type of boundary condition.
void begin_reconstruct_schur()
Common code of begin method of Reconstruct Schur assembly (Darcy and Richards)
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
LocDofVec cr_disc_dofs_
Vector of local DOF indices pre-computed on different DOF handlers.
~MHMatrixAssemblyRichards()
Destructor.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
void boundary_side_integral(DHCellSide cell_side)
Assembles between boundary element and corresponding side on bulk element.
RichardsLMH::EqFields EqFields
void update_water_content(const DHCellAccessor &cell, BulkPoint &p)
MHMatrixAssemblyRichards(EqFields *eq_fields, EqData *eq_data)
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Integral over element.
LocDofVec edge_indices_
Dofs of discontinuous fields on element edges.
void begin() override
Implements AssemblyBase::begin.
void end() override
Implements AssemblyBase::end.
void postprocess_velocity_richards(const DHCellAccessor &dh_cell, BulkPoint &p, arma::vec &solution)
Postprocess velocity after calculating of cell integral.
bool reset_soil_model(const DHCellAccessor &cell, BulkPoint &p)
RichardsLMH::EqData EqData
static constexpr const char * name()
void asm_source_term_richards(const DHCellAccessor &cell, BulkPoint &p)
Part of cell_integral method, specialized in Richards equation.
double compute_conductivity(const DHCellAccessor &cell, BulkPoint &p)
Precompute conductivity on bulk point.
EqFields * eq_fields_
Data objects shared with ConvectionTransport.
void end() override
Implements AssemblyBase::end.
void begin() override
Implements AssemblyBase::begin.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Integral over element.
void boundary_side_integral(FMT_UNUSED DHCellSide cell_side)
Assembles between boundary element and corresponding side on bulk element.
void dimjoin_intergral(FMT_UNUSED DHCellAccessor cell_lower_dim, FMT_UNUSED DHCellSide neighb_side)
Assembles the fluxes between elements of different dimensions.
static constexpr const char * name()
ReconstructSchurAssemblyRichards(EqFields *eq_fields, EqData *eq_data)
unsigned int bulk_idx() const
Returns index of the region in the bulk set.
Definition: region.hh:90
std::shared_ptr< SoilModelBase > soil_model_
Definition: richards_lmh.hh:98
VectorMPI water_content_previous_time
Definition: richards_lmh.hh:95
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > water_content_ptr
Definition: richards_lmh.hh:82
Field< 3, FieldValue< 3 >::Scalar > genuchten_p_head_scale
Definition: richards_lmh.hh:77
Field< 3, FieldValue< 3 >::Scalar > water_content_saturated
Definition: richards_lmh.hh:75
Field< 3, FieldValue< 3 >::Scalar > water_content_residual
Definition: richards_lmh.hh:76
Field< 3, FieldValue< 3 >::Scalar > genuchten_n_exponent
Definition: richards_lmh.hh:78
Boundary cond() const
double get(unsigned int pos) const
Return value on given position.
Definition: vector_mpi.hh:105
void set(unsigned int pos, double val)
Set value on given position.
Definition: vector_mpi.hh:111
@ result_zeros
@ coupling
@ boundary
@ bulk
int IntIdx
Definition: index_types.hh:25
arma::Col< IntIdx > LocDofVec
Definition: index_types.hh:28
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19
ArmaVec< double, N > vec
Definition: armor.hh:885
#define FMT_UNUSED
Definition: posix.h:75
double Ks
Definition: soil_models.hh:45
double Qr
Definition: soil_models.hh:42
double n
Definition: soil_models.hh:40
double alpha
Definition: soil_models.hh:41
double Qs
Definition: soil_models.hh:43