Flow123d
JS_before_hm-1804-gf2ad740aa
|
Go to the documentation of this file.
8 #ifndef SRC_FLOW_DARCY_FLOW_ASSEMBLY_HH_
9 #define SRC_FLOW_DARCY_FLOW_ASSEMBLY_HH_
53 template<
template<
int dim>
class Impl,
class Data>
55 return { std::make_shared<Impl<1> >(data),
56 std::make_shared<Impl<2> >(data),
57 std::make_shared<Impl<3> >(data) };
105 unsigned int nsides = dim+1;
109 for(
unsigned int i = 0; i < nsides; i++){
118 sp.submat(0, 0, nsides, nsides).ones();
125 sp.submat(0, nsides+1, nsides-1,
size()-1).diag().ones();
126 sp.submat(nsides+1, 0,
size()-1, nsides-1).diag().ones();
175 if (
ad_->balance !=
nullptr)
194 2*
ad_->conductivity.value( ele.
centre(), ele) *
195 arma::dot(
ad_->anisotropy.value( ele.
centre(), ele)*nv, nv) *
196 ad_->cross_section.value( neighb_side.
centre(), ele_higher ) *
197 ad_->cross_section.value( neighb_side.
centre(), ele_higher ) /
198 ad_->cross_section.value( ele.
centre(), ele ) *
226 const unsigned int nsides = ele->
n_sides();
229 unsigned int side_row, edge_row;
232 for (
unsigned int i = 0; i < nsides; i++) {
246 double cross_section =
ad_->cross_section.value(ele.
centre(), ele);
251 double bc_pressure =
ad_->bc_pressure.value(b_ele.
centre(), b_ele);
257 double bc_flux = -
ad_->bc_flux.value(b_ele.
centre(), b_ele);
258 double bc_pressure =
ad_->bc_pressure.value(b_ele.
centre(), b_ele);
259 double bc_sigma =
ad_->bc_robin_sigma.value(b_ele.
centre(), b_ele);
263 -b_ele.
measure() * bc_sigma * cross_section,
264 (bc_flux - bc_sigma * bc_pressure) * b_ele.
measure() * cross_section);
267 ad_->is_linear=
false;
270 char & switch_dirichlet =
ad_->bc_switch_dirichlet[loc_edge_idx];
271 double bc_pressure =
ad_->bc_switch_pressure.value(b_ele.
centre(), b_ele);
272 double bc_flux = -
ad_->bc_flux.value(b_ele.
centre(), b_ele);
273 double side_flux = bc_flux * b_ele.
measure() * cross_section;
276 if (switch_dirichlet) {
284 if ( solution_flux < side_flux) {
286 solution_flux = side_flux;
303 if ( solution_head > bc_pressure) {
305 solution_head = bc_pressure;
312 if (switch_dirichlet ||
ad_->force_no_neumann_bc ) {
322 ad_->is_linear=
false;
324 double bc_pressure =
ad_->bc_pressure.value(b_ele.
centre(), b_ele);
325 double bc_switch_pressure =
ad_->bc_switch_pressure.value(b_ele.
centre(), b_ele);
326 double bc_flux = -
ad_->bc_flux.value(b_ele.
centre(), b_ele);
327 double bc_sigma =
ad_->bc_robin_sigma.value(b_ele.
centre(), b_ele);
333 if (solution_head > bc_switch_pressure ||
ad_->force_no_neumann_bc) {
337 -b_ele.
measure() * bc_sigma * cross_section,
338 b_ele.
measure() * cross_section * (bc_flux - bc_sigma * bc_pressure) );
342 double bc_total_flux = bc_flux + bc_sigma*(bc_switch_pressure - bc_pressure);
348 THROW( ExcBCNotSupported() );
359 double cs =
ad_->cross_section.value(ele.
centre(), ele);
360 double conduct =
ad_->conductivity.value(ele.
centre(), ele);
361 double scale = 1 / cs /conduct;
376 for (
unsigned int k=0; k<qsize; k++)
377 for (
unsigned int i=0; i<ndofs; i++){
379 arma::dot(gravity_vec,velocity.value(i,k))
383 for (
unsigned int j=0; j<ndofs; j++){
385 arma::dot(velocity.value(i,k),
386 (
ad_->anisotropy.value(ele.
centre(), ele )).i()
387 * velocity.value(j,k))
405 for(
unsigned int i=0; i < ndofs; i++) {
406 double val_side = local_matrix(i,i);
407 double val_edge = -1./local_matrix(i,i);
411 static_cast<LinSys_BDDC*
>(
ad_->lin_sys)->diagonal_weights_set_value( side_row, val_side );
412 static_cast<LinSys_BDDC*
>(
ad_->lin_sys)->diagonal_weights_set_value( edge_row, val_edge );
421 for(
unsigned int side = 0; side <
loc_side_dofs.size(); side++){
436 auto ele = dh_cell.
elm();
455 higher_dim_dofs.resize(dh_neighb_cell.
n_dofs());
461 const unsigned int t = dh_neighb_cell.
n_dofs() - (dim+2) + neighb_side.side().side_idx();
474 static_cast<LinSys_BDDC*
>(
ad_->lin_sys)->diagonal_weights_set_value( ind, new_val );
483 unsigned int sidx = dh_side.side_idx();
484 if (dh_side.side().is_boundary()) {
485 ad_->balance->add_flux_values(
ad_->water_balance_idx, dh_side,
486 {local_dofs_[loc_side_dofs[sidx]]},
unsigned int n_neighs_vb() const
Return number of neighbours.
arma::Col< IntIdx > LocDofVec
Raviart-Thomas element of order 0.
Solver based on Multilevel BDDC - using corresponding class of OpenFTL package.
Definitions of Raviart-Thomas finite elements.
std::vector< unsigned int > loc_edge_dofs
const FEValuesViews::Vector< spacedim > & vector_view(unsigned int i) const
Accessor to vector values of multicomponent FE.
unsigned int n_points() const
Returns the number of quadrature points.
Side side() const
Return Side of given cell and side_idx.
void add_value(uint row, uint col, double mat_val, double rhs_val)
Adds a single entry into the local system.
virtual void fix_velocity(const DHCellAccessor &dh_cell)=0
arma::vec::fixed< spacedim > normal_vector(unsigned int point_no)
Returns the normal vector to a side at given quadrature point.
LocalSystem loc_system_vb_
double JxW(const unsigned int point_no)
Return the product of Jacobian determinant and the quadrature weight at given quadrature point.
bool is_boundary() const
Returns true for side on the boundary.
unsigned int get_dof_indices(std::vector< LongIdx > &indices) const
Fill vector of the global indices of dofs associated to the cell.
Class FEValues calculates finite element data on the actual cells such as shape function values,...
static constexpr bool value
Wrappers for linear systems based on MPIAIJ and MATIS format.
@ update_values
Shape function values.
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 initialize(Quadrature &_quadrature, FiniteElement< DIM > &_fe, UpdateFlags _flags)
Initialize structures and calculates cell-independent data.
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
std::shared_ptr< DarcyMH::EqData > AssemblyDataPtrMH
@ update_quadrature_points
Transformed quadrature points.
#define ASSERT_EQ_DBG(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
void assemble(const DHCellAccessor &dh_cell) override
unsigned int dim() const
Return dimension of element appropriate to cell.
@ update_normal_vectors
Normal vectors.
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
FEValues< 3 > fe_side_values_
void assemble_reconstruct(const DHCellAccessor &) override
virtual void assemble_reconstruct(const DHCellAccessor &dh_cell)=0
static MultidimAssembly create(Data data)
arma::vec3 centre() const
Side centre.
void set_dofs_and_bc(const DHCellAccessor &dh_cell)
QGauss velocity_interpolation_quad_
Definitions of basic Lagrangean finite elements with polynomial shape functions.
Side accessor allows to iterate over sides of DOF handler cell.
FEValues< 3 > velocity_interpolation_fv_
ArmaMat< double, N, M > mat
unsigned int n_dofs() const
Returns the number of shape functions.
Range< DHCellSide > side_range() const
Returns range of cell sides.
std::vector< LongIdx > global_dofs_
RangeConvert< DHNeighbSide, DHCellSide > neighb_sides() const
Returns range of neighbour cell of lower dimension corresponding to cell of higher dimension.
virtual void update_water_content(const DHCellAccessor &dh_cell)=0
Updates water content in Richards.
virtual void assemble(const DHCellAccessor &dh_cell)=0
Definitions of particular quadrature rules on simplices.
void assemble_sides(const DHCellAccessor &dh_cell)
Assembly explicit Schur complement for the given linear system. Provides method for resolution of the...
const ElementAccessor< 3 > elm() const
Return ElementAccessor to element of loc_ele_idx_.
std::shared_ptr< MortarAssemblyBase > mortar_assembly
unsigned int n_dofs() const
Return number of dofs on given cell.
ElementAccessor< 3 > element_accessor()
FE_P_disc< dim+1 > fe_p_disc_
void assembly_local_vb(ElementAccessor< 3 > ele, DHCellSide neighb_side)
Symmetric Gauss-Legendre quadrature formulae on simplices.
void assembly_dim_connections(const DHCellAccessor &dh_cell)
ElementAccessor< 3 > element() const
std::vector< unsigned int > dirichlet_edge
unsigned int n_sides() const
void eliminate_solution()
Cell accessor allow iterate over DOF handler cells.
void assemble_element(const DHCellAccessor &)
void reset()
Resets the matrix, RHS, dofs to zero and clears solution settings.
void update_water_content(const DHCellAccessor &) override
Updates water content in Richards.
std::vector< unsigned int > loc_side_dofs
mixed-hybrid model of linear Darcy flow, possibly unsteady.
@ update_JxW_values
Transformed quadrature weights.
#define ASSERT_LT_DBG(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
LocDofVec get_loc_dof_indices() const
Returns the local indices of dofs associated to the cell on the local process.
const arma::mat & get_matrix()
static unsigned int size()
unsigned int idx() const
Return local idx of element in boundary / bulk part of element vector.
void fix_velocity(const DHCellAccessor &dh_cell) override
void add_fluxes_in_balance_matrix(const DHCellAccessor &dh_cell)
void reinit(const ElementAccessor< spacedim > &cell)
Update cell-dependent data (gradients, Jacobians etc.)
void set_sparsity(const arma::umat &sp)
Sets the sparsity pattern for the local system.
SideIter side(const unsigned int loc_index)
double * get_solution_array()
NeighSideValues<(dim< 3) ? dim :2 > ngh_values_
DECLARE_EXCEPTION(ExcBCNotSupported,<< "BC type not supported.\n")
void assemble_sides_scale(const DHCellAccessor &dh_cell, double scale)
double measure() const
Computes the measure of the element.
arma::vec::fixed< spacedim > centre() const
Computes the barycenter.