8 #ifndef SRC_FLOW_DARCY_FLOW_ASSEMBLY_HH_
9 #define SRC_FLOW_DARCY_FLOW_ASSEMBLY_HH_
29 #include "flow/darcy_flow_mh.hh"
30 #include "flow/mortar_assembly.hh"
53 template<
template<
int dim>
class Impl,
class Fields,
class Data>
55 return { std::make_shared<Impl<1> >(eq_fields, eq_data),
56 std::make_shared<Impl<2> >(eq_fields, eq_data),
57 std::make_shared<Impl<3> >(eq_fields, eq_data) };
107 unsigned int nsides = dim+1;
111 for(
unsigned int i = 0; i < nsides; i++){
120 sp.submat(0, 0, nsides, nsides).ones();
127 sp.submat(0, nsides+1, nsides-1,
size()-1).diag().ones();
128 sp.submat(nsides+1, 0,
size()-1, nsides-1).diag().ones();
165 unsigned int nsides = dim+1;
169 arma::vec schur_solution =
ad_->full_solution.get_subvec(loc_dofs.subvec(nsides+1,
size()-1));
173 reconstructed_solution.zeros(schur_offset);
181 ad_->full_solution.set_subvec(loc_dofs.head(schur_offset), reconstructed_solution);
215 if (
ad_->balance !=
nullptr)
234 2*
af_->conductivity.value( ele.
centre(), ele) *
235 arma::dot(
af_->anisotropy.value( ele.
centre(), ele)*nv, nv) *
236 af_->cross_section.value( neighb_side.
centre(), ele_higher ) *
237 af_->cross_section.value( neighb_side.
centre(), ele_higher ) /
238 af_->cross_section.value( ele.
centre(), ele ) *
266 const unsigned int nsides = ele->
n_sides();
269 unsigned int side_row, edge_row;
272 for (
unsigned int i = 0; i < nsides; i++) {
286 double cross_section =
af_->cross_section.value(ele.
centre(), ele);
291 double bc_pressure =
af_->bc_pressure.value(b_ele.
centre(), b_ele);
297 double bc_flux = -
af_->bc_flux.value(b_ele.
centre(), b_ele);
298 double bc_pressure =
af_->bc_pressure.value(b_ele.
centre(), b_ele);
299 double bc_sigma =
af_->bc_robin_sigma.value(b_ele.
centre(), b_ele);
303 -b_ele.
measure() * bc_sigma * cross_section,
304 (bc_flux - bc_sigma * bc_pressure) * b_ele.
measure() * cross_section);
307 ad_->is_linear=
false;
310 char & switch_dirichlet =
ad_->bc_switch_dirichlet[loc_edge_idx];
311 double bc_pressure =
af_->bc_switch_pressure.value(b_ele.
centre(), b_ele);
312 double bc_flux = -
af_->bc_flux.value(b_ele.
centre(), b_ele);
313 double side_flux = bc_flux * b_ele.
measure() * cross_section;
316 if (use_dirichlet_switch) {
317 if (switch_dirichlet) {
325 if ( solution_flux < side_flux) {
327 solution_flux = side_flux;
344 if ( solution_head > bc_pressure) {
346 solution_head = bc_pressure;
354 if (switch_dirichlet ||
ad_->force_no_neumann_bc ) {
364 ad_->is_linear=
false;
366 double bc_pressure =
af_->bc_pressure.value(b_ele.
centre(), b_ele);
367 double bc_switch_pressure =
af_->bc_switch_pressure.value(b_ele.
centre(), b_ele);
368 double bc_flux = -
af_->bc_flux.value(b_ele.
centre(), b_ele);
369 double bc_sigma =
af_->bc_robin_sigma.value(b_ele.
centre(), b_ele);
375 if (solution_head > bc_switch_pressure ||
ad_->force_no_neumann_bc) {
379 -b_ele.
measure() * bc_sigma * cross_section,
380 b_ele.
measure() * cross_section * (bc_flux - bc_sigma * bc_pressure) );
384 double bc_total_flux = bc_flux + bc_sigma*(bc_switch_pressure - bc_pressure);
390 THROW( ExcBCNotSupported() );
401 double cs =
af_->cross_section.value(ele.
centre(), ele);
402 double conduct =
af_->conductivity.value(ele.
centre(), ele);
403 double scale = 1 / cs /conduct;
418 for (
unsigned int k=0; k<qsize; k++)
419 for (
unsigned int i=0; i<ndofs; i++){
421 arma::dot(gravity_vec,velocity.value(i,k))
425 for (
unsigned int j=0; j<ndofs; j++){
427 arma::dot(velocity.value(i,k),
428 (
af_->anisotropy.value(ele.
centre(), ele )).i()
429 * velocity.value(j,k))
447 for(
unsigned int i=0; i < ndofs; i++) {
448 double val_side = local_matrix(i,i);
449 double val_edge = -1./local_matrix(i,i);
453 static_cast<LinSys_BDDC*
>(
ad_->lin_sys)->diagonal_weights_set_value( side_row, val_side );
454 static_cast<LinSys_BDDC*
>(
ad_->lin_sys)->diagonal_weights_set_value( edge_row, val_edge );
463 for(
unsigned int side = 0; side <
loc_side_dofs.size(); side++){
478 auto ele = dh_cell.
elm();
497 higher_dim_dofs.resize(dh_neighb_cell.
n_dofs());
503 const unsigned int t = dh_neighb_cell.
n_dofs() - (dim+2) + neighb_side.side().side_idx();
516 static_cast<LinSys_BDDC*
>(
ad_->lin_sys)->diagonal_weights_set_value( ind, new_val );
525 unsigned int sidx = dh_side.side_idx();
526 if (dh_side.side().is_boundary()) {
527 ad_->balance->add_flux_values(
ad_->water_balance_idx, dh_side,
528 {local_dofs_[loc_side_dofs[sidx]]},
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
virtual void assemble_reconstruct(const DHCellAccessor &dh_cell)=0
static MultidimAssembly create(Fields eq_fields, Data eq_data)
DECLARE_EXCEPTION(ExcBCNotSupported,<< "BC type not supported.\n")
virtual void fix_velocity(const DHCellAccessor &dh_cell)=0
virtual void update_water_content(const DHCellAccessor &dh_cell)=0
Updates water content in Richards.
virtual void assemble(const DHCellAccessor &dh_cell)=0
virtual ~AssemblyFlowBase()
void update_water_content(const DHCellAccessor &) override
Updates water content in Richards.
void assembly_dim_connections(const DHCellAccessor &dh_cell)
std::vector< unsigned int > dirichlet_edge
std::shared_ptr< DarcyMH::EqData > AssemblyDataPtrMH
void set_dofs_and_bc(const DHCellAccessor &dh_cell, bool use_dirichlet_switch)
void assemble_element(const DHCellAccessor &)
LocalSystem loc_system_vb_
void assemble_reconstruct(const DHCellAccessor &dh_cell) override
void fix_velocity(const DHCellAccessor &dh_cell) override
NeighSideValues<(dim< 3) ? dim :2 > ngh_values_
static unsigned int size()
std::vector< unsigned int > loc_side_dofs
FEValues< 3 > velocity_interpolation_fv_
std::shared_ptr< MortarAssemblyBase > mortar_assembly
void add_fluxes_in_balance_matrix(const DHCellAccessor &dh_cell)
std::vector< unsigned int > loc_edge_dofs
void assemble_sides(const DHCellAccessor &dh_cell)
QGauss velocity_interpolation_quad_
std::shared_ptr< DarcyMH::EqFields > AssemblyFieldsPtrMH
void assembly_local_vb(ElementAccessor< 3 > ele, DHCellSide neighb_side)
void assemble_sides_scale(const DHCellAccessor &dh_cell, double scale)
std::vector< LongIdx > global_dofs_
void assemble(const DHCellAccessor &dh_cell) override
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_.
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.
LocDofVec get_loc_dof_indices() const
Returns the local indices of dofs associated to the cell on the local process.
unsigned int get_dof_indices(std::vector< LongIdx > &indices) const
Fill vector of the global indices of dofs associated to the cell.
unsigned int dim() const
Return dimension of element appropriate to cell.
Range< DHCellSide > side_range() const
Returns range of cell sides.
Side accessor allows to iterate over sides of DOF handler cell.
Side side() const
Return Side of given cell and side_idx.
ElementAccessor< 3 > element() const
arma::vec3 centre() const
Side centre.
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...
arma::vec::fixed< spacedim > centre() const
Computes the barycenter.
unsigned int n_neighs_vb() const
Return number of neighbours.
unsigned int n_sides() const
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.
unsigned int n_dofs() const
Returns the number of shape functions.
Raviart-Thomas element of order 0.
double * get_solution_array()
void set_sparsity(const arma::umat &sp)
Sets the sparsity pattern for the local system.
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()
const arma::mat & get_matrix()
void reconstruct_solution_schur(uint offset, const arma::vec &schur_solution, arma::vec &reconstructed_solution) const
Reconstructs the solution from the Schur complement solution: x = invA*b - invA * Bt * schur_solution...
void add_value(uint row, uint col, double mat_val, double rhs_val)
Adds a single entry into the local system.
void reset()
Resets the matrix, RHS, dofs to zero and clears solution settings.
FE_P_disc< dim+1 > fe_p_disc_
FEValues< 3 > fe_side_values_
Symmetric Gauss-Legendre quadrature formulae on simplices.
bool is_boundary() const
Returns true for side on the boundary.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
Definitions of Raviart-Thomas 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
static constexpr bool value
Wrappers for linear systems based on MPIAIJ and MATIS format.
Solver based on Multilevel BDDC - using corresponding class of OpenFTL package.
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
ArmaMat< double, N, M > mat
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.