18 #ifndef ASSEMBLY_LMH_HH_
19 #define ASSEMBLY_LMH_HH_
45 template <
unsigned int dim>
52 static constexpr
const char *
name() {
return "ReadInitCondAssemblyLMH"; }
73 ASSERT_EQ(cell.
dim(), dim).error(
"Dimension of element mismatch!");
79 auto p = *( this->
bulk_points(element_patch_idx).begin() );
82 for (
unsigned int i=0; i<cell.
elm()->
n_sides(); i++) {
110 template <
template<
IntDim...>
class DimAssembly>
115 template <
unsigned int dim>
124 static constexpr
const char *
name() {
return "MHMatrixAssemblyLMH"; }
153 fe_ = std::make_shared< FE_P_disc<dim> >(0);
182 ASSERT_EQ(cell.
dim(), dim).error(
"Dimension of element mismatch!");
185 auto p = *( this->
bulk_points(element_patch_idx).begin() );
201 ASSERT_EQ(cell_side.
dim(), dim).error(
"Dimension of element mismatch!");
224 if (dim == 1)
return;
225 ASSERT_EQ(cell_lower_dim.
dim(), dim-1).error(
"Dimension of element mismatch!");
227 unsigned int neigh_idx =
ngh_idx(cell_lower_dim, neighb_side);
228 unsigned int loc_dof_higher = (2*(cell_lower_dim.
dim()+1) + 1) + neigh_idx;
233 auto p_low = p_high.lower_dim(cell_lower_dim);
296 unsigned int i_constr = 0;
344 schur_solution, this->reconstructed_solution_);
349 this->eq_data_->schur_offset_[dh_cell.dim()-1]), this->reconstructed_solution_);
351 this->loc_schur_.row_dofs.n_elem), schur_solution);
363 auto ele = cell.
elm();
377 arma::dot( velocity.value(i,k),
447 for (
unsigned int i=0; i<n_sides; i++)
479 cross_section_ = eq_fields_->cross_section(p_side);
483 side_row_ = eq_data_->loc_side_dofs[dim-1][sidx_];
484 edge_row_ = eq_data_->loc_edge_dofs[dim-1][sidx_];
490 char & switch_dirichlet = eq_data_->bc_switch_dirichlet[b_ele.
idx()];
491 if (switch_dirichlet) {
497 auto reconstr_solution = this->load_local_system(cell_side.
cell());
498 double side_flux = -eq_fields_->bc_flux(p_bdr) * b_ele.
measure() * cross_section_;
500 if ( reconstr_solution[side_row_] < side_flux) {
502 switch_dirichlet = 0;
517 double solution_head = eq_data_->p_edge_solution.get(loc_dof_vec[sidx_]);
518 double bc_pressure = eq_fields_->bc_switch_pressure(p_bdr);
520 if ( solution_head > bc_pressure) {
535 eq_data_->loc_constraint_.emplace_back( bulk_local_idx_, sidx_, eq_fields_->bc_pressure(p_bdr) );
539 eq_data_->loc_system_[bulk_local_idx_].add_value(edge_row_, edge_row_,
540 -b_ele.
measure() * eq_fields_->bc_robin_sigma(p_bdr) * cross_section_,
541 (-eq_fields_->bc_flux(p_bdr) - eq_fields_->bc_robin_sigma(p_bdr) * eq_fields_->bc_pressure(p_bdr)) * b_ele.
measure() * cross_section_);
543 eq_data_->is_linear=
false;
545 double bc_pressure = eq_fields_->bc_switch_pressure(p_bdr);
546 double side_flux = -eq_fields_->bc_flux(p_bdr) * b_ele.
measure() * cross_section_;
548 eq_data_->save_local_system_[bulk_local_idx_] = (bool) eq_data_->bc_switch_dirichlet[b_ele.
idx()];
552 if (eq_data_->save_local_system_[bulk_local_idx_] || eq_data_->force_no_neumann_bc ) {
554 eq_data_->loc_constraint_.emplace_back(bulk_local_idx_, sidx_, bc_pressure);
557 eq_data_->loc_system_[bulk_local_idx_].add_value(edge_row_, side_flux);
561 eq_data_->is_linear=
false;
563 double bc_pressure = eq_fields_->bc_pressure(p_bdr);
564 double bc_switch_pressure = eq_fields_->bc_switch_pressure(p_bdr);
565 double bc_flux = -eq_fields_->bc_flux(p_bdr);
566 double bc_sigma = eq_fields_->bc_robin_sigma(p_bdr);
570 double solution_head = eq_data_->p_edge_solution.get(loc_dof_vec[sidx_]);
573 if (solution_head > bc_switch_pressure || eq_data_->force_no_neumann_bc) {
576 eq_data_->loc_system_[bulk_local_idx_].add_value(edge_row_, edge_row_,
577 -b_ele.
measure() * bc_sigma * cross_section_,
578 b_ele.
measure() * cross_section_ * (bc_flux - bc_sigma * bc_pressure) );
582 double bc_total_flux = bc_flux + bc_sigma*(bc_switch_pressure - bc_pressure);
584 eq_data_->loc_system_[bulk_local_idx_].add_value(edge_row_, bc_total_flux * b_ele.
measure() * cross_section_);
588 THROW( ExcBCNotSupported() );
591 eq_data_->balance_->add_flux_values(eq_data_->water_balance_idx, cell_side,
592 {eq_data_->loc_system_[bulk_local_idx_].row_dofs[eq_data_->loc_side_dofs[dim-1][sidx_]]},
598 unsigned int size, loc_size, elm_dim;;
603 size = (elm_dim+1) + 1 + (elm_dim+1);
608 dofs.head(dh_cell.n_dofs()) = dh_cell.get_loc_dof_indices();
615 for (
DHCellSide neighb_side : dh_cell.neighb_sides() ) {
626 const unsigned int p = size+i;
628 const unsigned int t = dh_neighb_cell.
n_dofs() - dh_neighb_cr_cell.
n_dofs() + neighb_side.side().side_idx();
633 eq_data_->loc_system_[dh_cell.local_idx()].reset(dofs, dofs);
636 for (
DHCellSide dh_side : dh_cell.side_range()) {
637 unsigned int sidx = dh_side.side_idx();
639 eq_data_->loc_system_[dh_cell.local_idx()].add_value(eq_data_->loc_side_dofs[elm_dim-1][sidx], eq_data_->loc_edge_dofs[elm_dim-1][sidx], 1.0);
640 eq_data_->loc_system_[dh_cell.local_idx()].add_value(eq_data_->loc_edge_dofs[elm_dim-1][sidx], eq_data_->loc_side_dofs[elm_dim-1][sidx], 1.0);
668 const unsigned int tt = dh_cr_cell.
n_dofs()+i;
675 loc_schur_.reset(dofs_schur, dofs_schur);
681 for (
uint n_i=0; n_i<dh_cell.
elm()->n_neighs_vb(); ++n_i) {
683 if ( (side->elem_idx() == neighb_side.
elem_idx()) && (side->side_idx() == neighb_side.
side_idx()) )
return n_i;
698 auto ls = eq_data_->seepage_bc_systems.find(dh_cell.
elm_idx());
699 if (ls != eq_data_->seepage_bc_systems.end())
703 arma::vec schur_solution = eq_data_->p_edge_solution.get_subvec(loc_dof_vec);
705 ls->second.reconstruct_solution_schur(eq_data_->schur_offset_[dim-1], schur_solution, reconstructed_solution_);
707 unsigned int pos_in_cache = this->element_cache_map_->position_in_cache(dh_cell.
elm_idx());
708 auto p = *( this->bulk_points(pos_in_cache).begin() );
709 postprocess_velocity_darcy(dh_cell, p, reconstructed_solution_);
711 eq_data_->bc_fluxes_reconstruted[bulk_local_idx_] =
true;
714 return reconstructed_solution_;
726 * eq_fields_->cross_section(p)
729 edge_source_term_ = edge_scale_ *
730 ( eq_fields_->water_source_density(p)
731 + eq_fields_->extra_source(p));
738 postprocess_velocity(dh_cell, p);
744 for (
unsigned int i=0; i<dh_cell.
elm()->n_sides(); i++) {
746 if( ! eq_data_->use_steady_assembly_)
748 double new_pressure = eq_data_->p_edge_solution.get(loc_dof_vec[i]);
749 double old_pressure = eq_data_->p_edge_solution_previous_time.get(loc_dof_vec[i]);
750 time_term_ = edge_scale_ * (eq_fields_->storativity(p) + eq_fields_->extra_storativity(p))
751 / eq_data_->time_step_ * (new_pressure - old_pressure);
753 solution[eq_data_->loc_side_dofs[dim-1][i]] += edge_source_term_ - time_term_;
770 shared_ptr<FiniteElement<dim>>
fe_;
788 template <
template<
IntDim...>
class DimAssembly>
793 template <
unsigned int dim>
800 static constexpr
const char *
name() {
return "ReconstructSchurAssemblyLMH"; }
809 ASSERT_EQ(cell.
dim(), dim).error(
"Dimension of element mismatch!");
812 auto p = *( this->
bulk_points(element_patch_idx).begin() );
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
#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.
Quadrature * quad_low_
Quadrature used in assembling methods (dim-1).
Range< CouplingPoint > coupling_points(const DHCellSide &cell_side) const
Return CouplingPoint 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.
Point accessor allow iterate over quadrature points of given side defined in local element coordinate...
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_.
unsigned int n_dofs() const
Return number of dofs on given cell.
RangeConvert< DHNeighbSide, DHCellSide > neighb_sides() const
Returns range of neighbour cell of lower dimension corresponding to cell of higher dimension.
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 elm_idx() const
Return serial idx to element of loc_ele_idx_.
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.
unsigned int elem_idx() const
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.
unsigned int side_idx() const
std::array< unsigned int, 3 > schur_offset_
Index offset in the local system for the Schur complement (of dim = 1,2,3).
unsigned int nonlinear_iteration_
std::map< LongIdx, LocalSystem > seepage_bc_systems
std::vector< bool > save_local_system_
Flag for saving the local system. Currently used only in case of seepage BC.
void reset()
Reset data members.
std::vector< arma::vec > postprocess_solution_
bool use_steady_assembly_
std::shared_ptr< SubDOFHandlerMultiDim > dh_cr_
DOF handler represents DOFs of edges.
bool force_no_neumann_bc
auxiliary flag for switchting Dirichlet like BC
VectorMPI p_edge_solution_previous_time
VectorMPI p_edge_solution
std::shared_ptr< DOFHandlerMultiDim > dh_
full DOF handler represents DOFs of sides, elements and edges
std::vector< LocalConstraint > loc_constraint_
std::vector< LocalSystem > loc_system_
Following data members are stored in vectors, one item for every cell.
std::shared_ptr< LinSys > lin_sys_schur
std::array< unsigned int, 3 > loc_ele_dof
std::array< std::vector< unsigned int >, 3 > loc_edge_dofs
std::shared_ptr< Balance > balance_
Shared Balance object.
std::array< std::vector< unsigned int >, 3 > loc_side_dofs
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 > init_pressure
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
unsigned int n_sides() const
Returns number of sides aligned with the edge.
double measure() const
Computes the measure of the element.
SideIter side(const unsigned int loc_index)
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
const Element * element() const
Directing class of FieldValueCache.
unsigned int n_neighs_vb() const
Return number of neighbours.
unsigned int n_sides() const
Compound finite element on dim dimensional simplex.
std::vector< unsigned int > fe_dofs(unsigned int fe_index)
Return dof indices belonging to given sub-FE.
double JxW(const unsigned int point_no)
Return the product of Jacobian determinant and the quadrature weight at given quadrature point.
const FEValuesViews::Vector< spacedim > & vector_view(unsigned int i) const
Accessor to vector values of multicomponent FE.
void initialize(Quadrature &_quadrature, FiniteElement< DIM > &_fe, UpdateFlags _flags)
Initialize structures and calculates cell-independent data.
void reinit(const ElementAccessor< spacedim > &cell)
Update cell-dependent data (gradients, Jacobians etc.)
unsigned int n_points() const
Returns the number of quadrature points.
arma::vec::fixed< spacedim > normal_vector(unsigned int point_no)
Returns the normal vector to a side at given quadrature point.
unsigned int n_dofs() const
Returns the number of shape functions.
Raviart-Thomas element of order 0.
Container for various descendants of FieldCommonBase.
Generic class of assemblation.
void set_solution(uint loc_dof, double solution, double diag=0.0)
Set the position and value of known solution. E.g. Dirichlet boundary condition.
void eliminate_solution()
virtual ~MHMatrixAssemblyLMH()
Destructor.
void end_mh_matrix()
Common code of end method of MH matrix assembly (Darcy and Richards)
arma::vec3 nv_
Normal vector.
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 boundary_side_integral(DHCellSide cell_side)
DECLARE_EXCEPTION(ExcBCNotSupported,<< "BC type not supported.\n")
unsigned int ngh_idx(DHCellAccessor &dh_cell, DHCellSide &neighb_side)
Temporary method find neighbour index in higher-dim cell.
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side)
arma::vec reconstructed_solution_
Vector for reconstructed solution (velocity and pressure on element) from Schur complement.
void asm_sides(const DHCellAccessor &cell, BulkPoint &p, double conductivity)
Part of cell_integral method, common in all descendants.
FEValues< 3 > fe_values_side_
FEValues of object (of P disc finite element type)
void set_dofs()
Precompute loc_system and loc_schur data members.
double ngh_value_
Precomputed ngh value.
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)
FE_RT0< dim > fe_rt_
Assembly volume integrals.
void begin() override
Implements AssemblyBase::begin.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
void precompute_boundary_side(DHCellSide &cell_side, BoundaryPoint &p_side, BulkPoint &p_bdr)
Precompute values used in boundary side integral on given DHCellSide.
EqFields * eq_fields_
Data objects shared with DarcyFlow.
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)
MHMatrixAssemblyLMH(EqFields *eq_fields, EqData *eq_data)
Constructor.
DarcyLMH::EqFields::BC_Type type_
Type of boundary condition.
double time_term_rhs_
Variables used in compute time term (unsteady)
void postprocess_velocity_darcy(const DHCellAccessor &dh_cell, BulkPoint &p, arma::vec &solution)
Postprocess velocity during loading of local system and after calculating of cell integral.
void set_loc_schur(const DHCellAccessor dh_cr_cell)
Precompute loc_schur data member of given cell.
double cross_section_
Precomputed cross-section value.
arma::vec load_local_system(const DHCellAccessor &dh_cell)
static constexpr const char * name()
void begin_reconstruct_schur()
Common code of begin method of Reconstruct Schur assembly (Darcy and Richards)
unsigned int edge_row_
Helper indices in boundary assembly.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
DarcyLMH::EqFields EqFields
void end() override
Implements AssemblyBase::end.
void asm_source_term_darcy(const DHCellAccessor &cell, BulkPoint &p)
Part of cell_integral method, specialized in Darcy equation.
double source_term_
Variables used in compute lumped source.
void postprocess_velocity(const DHCellAccessor &dh_cell, BulkPoint &p)
Symmetric Gauss-Legendre quadrature formulae on simplices.
FieldSet used_fields_
Sub field set contains fields used in calculation.
static constexpr const char * name()
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
ReadInitCondAssemblyLMH(EqFields *eq_fields, EqData *eq_data)
Constructor.
void end() override
Implements AssemblyBase::end.
DarcyLMH::EqFields EqFields
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
LocDofVec l_indices_
Vector of pre-computed local DOF indices.
virtual ~ReadInitCondAssemblyLMH()
Destructor.
EqFields * eq_fields_
Data objects shared with Flow equation.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Integral over element.
static constexpr const char * name()
void begin() override
Implements AssemblyBase::begin.
void dimjoin_intergral(FMT_UNUSED DHCellAccessor cell_lower_dim, FMT_UNUSED DHCellSide neighb_side)
Assembles the fluxes between elements of different dimensions.
void boundary_side_integral(FMT_UNUSED DHCellSide cell_side)
Assembles between boundary element and corresponding side on bulk element.
ReconstructSchurAssemblyLMH(EqFields *eq_fields, EqData *eq_data)
Constructor.
void end() override
Implements AssemblyBase::end.
DarcyLMH::EqFields EqFields
Edge edge() const
Returns pointer to the edge connected to the side.
arma::vec get_subvec(const LocDofVec &loc_indices)
double get(unsigned int pos) const
Return value on given position.
void local_to_ghost_end()
local_to_ghost_{begin,end} updates the ghost values on neighbouring processors from local values
void set_subvec(const LocDofVec &loc_indices, const arma::vec &values)
void add(unsigned int pos, double val)
Add value to item on given position.
void local_to_ghost_begin()
local_to_ghost_{begin,end} updates the ghost values on neighbouring processors from local values
void ghost_to_local_end()
ghost_to_local_{begin,end} updates the local values by adding ghost values from neighbouring processo...
void ghost_to_local_begin()
ghost_to_local_{begin,end} updates the local values by adding ghost values from neighbouring processo...
Lumped mixed-hybrid model of linear Darcy flow, possibly unsteady.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
Definitions of Raviart-Thomas finite elements.
Class FESystem for compound finite elements.
Class FEValues calculates finite element data on the actual cells such as shape function values,...
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
arma::Col< IntIdx > LocDofVec
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
unsigned int IntDim
A dimension index type.
Definitions of particular quadrature rules on simplices.
Assembly explicit Schur complement for the given linear system. Provides method for resolution of the...
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_quadrature_points
Transformed quadrature points.