18 #ifndef ASSEMBLY_RICHARDS_HH_
19 #define ASSEMBLY_RICHARDS_HH_
28 template <
unsigned int dim>
35 static constexpr
const char *
name() {
return "InitCondPostprocessAssembly"; }
61 ASSERT_EQ(cell.
dim(), dim).error(
"Dimension of element mismatch!");
67 auto p = *( this->
bulk_points(element_patch_idx).begin() );
74 for (
unsigned int i=0; i<cell.
elm()->n_sides(); i++) {
81 fadbad::B<double> x_phead(
phead);
139 template <
template<
IntDim...>
class DimAssembly>
144 template <
unsigned int dim>
151 static constexpr
const char *
name() {
return "MHMatrixAssemblyRichards"; }
189 ASSERT_EQ(cell.
dim(), dim).error(
"Dimension of element mismatch!");
192 auto p = *( this->
bulk_points(element_patch_idx).begin() );
204 ASSERT_EQ(cell_side.
dim(), dim).error(
"Dimension of element mismatch!");
251 for (
unsigned int i=0; i<ele->
n_sides(); i++)
295 bool genuchten_on = (this->eq_fields_->genuchten_p_head_scale.field_result({cell.
elm().region()}) !=
result_zeros);
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);
305 this->eq_data_->soil_model_->reset(soil_data);
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();
320 for (
unsigned int i=0; i<cell.
elm()->n_sides(); i++) {
323 phead = this->eq_data_->p_edge_solution.
get( edge_indices_[i] );
326 fadbad::B<double> x_phead(phead);
327 fadbad::B<double> evaluated( this->eq_data_->soil_model_->water_content_diff(x_phead) );
329 water_content = evaluated.val();
330 capacity = x_phead.d(0);
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);
340 bool genuchten_on = reset_soil_model(cell, p);
342 double conductivity = 0;
347 for (
unsigned int i=0; i<cell.
elm()->n_sides(); i++)
349 double phead = eq_data_->p_edge_solution.get( loc_dof_vec[i] );
350 conductivity += eq_data_->soil_model_->conductivity(phead);
355 conductivity = eq_fields_->conductivity(p);
364 this->postprocess_velocity(dh_cell, p);
366 this->update_water_content(dh_cell, p);
368 VectorMPI water_content_vec = eq_fields_->water_content_ptr->vec();
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] );
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_;
379 eq_fields_->conductivity_ptr->vec().set( p_dof, compute_conductivity(dh_cell, p) );
395 template <
template<
IntDim...>
class DimAssembly>
400 template <
unsigned int dim>
407 static constexpr
const char *
name() {
return "ReconstructSchurAssemblyRichards"; }
416 ASSERT_EQ(cell.
dim(), dim).error(
"Dimension of element mismatch!");
419 auto p = *( this->
bulk_points(element_patch_idx).begin() );
450 template <
template<
IntDim...>
class DimAssembly>
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
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.
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.
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.
Directing class of FieldValueCache.
unsigned int n_sides() const
Container for various descendants of FieldCommonBase.
FieldResult field_result(RegionSet region_set) const override
Indicates special field states.
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::EqData EqData
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.
~InitCondPostprocessAssembly()
Destructor.
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.
double water_content_diff_
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.
RichardsLMH::EqData EqData
static constexpr const char * name()
RichardsLMH::EqFields EqFields
ReconstructSchurAssemblyRichards(EqFields *eq_fields, EqData *eq_data)
unsigned int bulk_idx() const
Returns index of the region in the bulk set.
std::shared_ptr< SoilModelBase > soil_model_
VectorMPI water_content_previous_time
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > water_content_ptr
Field< 3, FieldValue< 3 >::Scalar > genuchten_p_head_scale
Field< 3, FieldValue< 3 >::Scalar > water_content_saturated
Field< 3, FieldValue< 3 >::Scalar > water_content_residual
Field< 3, FieldValue< 3 >::Scalar > genuchten_n_exponent
double get(unsigned int pos) const
Return value on given position.
void set(unsigned int pos, double val)
Set value on given position.
arma::Col< IntIdx > LocDofVec
unsigned int IntDim
A dimension index type.