Flow123d  master-1fea4ce
assembly_lmh.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_lmh.hh
15  * @brief
16  */
17 
18 #ifndef ASSEMBLY_LMH_HH_
19 #define ASSEMBLY_LMH_HH_
20 
23 #include "flow/darcy_flow_lmh.hh"
25 
26 #include "system/index_types.hh"
27 #include "mesh/mesh.h"
28 #include "mesh/accessors.hh"
29 #include "mesh/neighbours.h"
30 #include "fem/fe_p.hh"
31 #include "fem/fe_values.hh"
32 #include "fem/fe_rt.hh"
33 #include "fem/fe_values_views.hh"
34 #include "fem/fe_system.hh"
36 
37 #include "la/linsys_PETSC.hh"
38 #include "la/schur.hh"
39 #include "la/local_system.hh"
40 #include "la/local_constraint.hh"
41 
42 #include "coupling/balance.hh"
43 
44 
45 template <unsigned int dim>
47 {
48 public:
49  typedef typename DarcyLMH::EqFields EqFields;
50  typedef typename DarcyLMH::EqData EqData;
51 
52  static constexpr const char * name() { return "ReadInitCondAssemblyLMH"; }
53 
54  /// Constructor.
55  ReadInitCondAssemblyLMH(EqFields *eq_fields, EqData *eq_data)
56  : AssemblyBase<dim>(0), eq_fields_(eq_fields), eq_data_(eq_data) {
59  }
60 
61  /// Destructor.
63 
64  /// Initialize auxiliary vectors and other data members
65  void initialize(ElementCacheMap *element_cache_map) {
66  this->element_cache_map_ = element_cache_map;
67  }
68 
69 
70  /// Assemble integral over element
71  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
72  {
73  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
74 
76  ASSERT(l_indices_.n_elem == cell.elm().element()->n_sides());
77 
78  // set initial condition
79  auto p = *( this->bulk_points(element_patch_idx).begin() );
80  double init_value = eq_fields_->init_pressure(p);
81 
82  for (unsigned int i=0; i<cell.elm()->n_sides(); i++) {
83  double init_value_on_edge = init_value / cell.elm().side(i)->edge().n_sides();
84  eq_data_->p_edge_solution.add(l_indices_[i], init_value_on_edge);
85  }
86  }
87 
88  /// Implements @p AssemblyBase::end.
89  void end() override
90  {
95  }
96 
97 
98 
99 protected:
100  /// Sub field set contains fields used in calculation.
102 
103  /// Data objects shared with Flow equation
106 
107  /// Vector of pre-computed local DOF indices
109 
110  template < template<IntDim...> class DimAssembly>
111  friend class GenericAssembly;
112 
113 };
114 
115 template <unsigned int dim>
117 {
118 public:
119  typedef typename DarcyLMH::EqFields EqFields;
120  typedef typename DarcyLMH::EqData EqData;
121 
122  DECLARE_EXCEPTION( ExcBCNotSupported, << "BC type not supported.\n" );
123 
124  static constexpr const char * name() { return "MHMatrixAssemblyLMH"; }
125 
126  /// Constructor.
127  MHMatrixAssemblyLMH(EqFields *eq_fields, EqData *eq_data)
128  : AssemblyBase<dim>(0), eq_fields_(eq_fields), eq_data_(eq_data), quad_rt_(dim, 2) {
133  this->used_fields_ += eq_fields_->sigma;
138  this->used_fields_ += eq_fields_->bc_type;
140  this->used_fields_ += eq_fields_->bc_flux;
143  }
144 
145  /// Destructor.
146  virtual ~MHMatrixAssemblyLMH() {}
147 
148  /// Initialize auxiliary vectors and other data members
149  void initialize(ElementCacheMap *element_cache_map) {
150  //this->balance_ = eq_data_->balance_;
151  this->element_cache_map_ = element_cache_map;
152 
153  fe_ = std::make_shared< FE_P_disc<dim> >(0);
155 
157 
158  // local numbering of dofs for MH system
159  // note: this shortcut supposes that the fe_system is the same on all elements
160  // TODO the function should be DiscreteSpace.fe(ElementAccessor)
161  // and correct form fe(ad_->dh_->own_range().begin()->elm()) (see documentation of DiscreteSpace::fe)
162  auto fe = eq_data_->dh_->ds()->fe()[Dim<dim>{}];
163  FESystem<dim>* fe_system = dynamic_cast<FESystem<dim>*>(fe.get());
164  eq_data_->loc_side_dofs[dim-1] = fe_system->fe_dofs(0);
165  eq_data_->loc_ele_dof[dim-1] = fe_system->fe_dofs(1)[0];
166  eq_data_->loc_edge_dofs[dim-1] = fe_system->fe_dofs(2);
167 
168  // reconstructed vector for the velocity and pressure
169  // length = local Schur complement offset in LocalSystem
170  eq_data_->schur_offset_[dim-1] = eq_data_->loc_edge_dofs[dim-1][0];
172  }
173 
174 
175  /**
176  * Integral over element.
177  *
178  * Override in descendants.
179  */
180  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
181  {
182  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
183 
184  // evaluation point
185  auto p = *( this->bulk_points(element_patch_idx).begin() );
186  bulk_local_idx_ = cell.local_idx();
187 
188  this->asm_sides(cell, p, eq_fields_->conductivity(p));
189  this->asm_element();
190  this->asm_source_term_darcy(cell, p);
191  }
192 
193 
194  /**
195  * Assembles between boundary element and corresponding side on bulk element.
196  *
197  * Override in descendants.
198  */
199  inline void boundary_side_integral(DHCellSide cell_side)
200  {
201  ASSERT_EQ(cell_side.dim(), dim).error("Dimension of element mismatch!");
202  if (!cell_side.cell().is_own()) return;
203 
204  auto p_side = *( this->boundary_points(cell_side).begin() );
205  auto p_bdr = p_side.point_bdr(cell_side.cond().element_accessor() );
206  ElementAccessor<3> b_ele = cell_side.side().cond().element_accessor(); // ??
207 
208  precompute_boundary_side(cell_side, p_side, p_bdr);
209 
211  this->use_dirichlet_switch(cell_side, b_ele, p_bdr);
212  }
213 
214  this->boundary_side_integral_in(cell_side, b_ele, p_bdr);
215  }
216 
217 
218  /**
219  * Assembles the fluxes between elements of different dimensions.
220  *
221  * Common in all descendants.
222  */
223  inline void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side) {
224  if (dim == 1) return;
225  ASSERT_EQ(cell_lower_dim.dim(), dim-1).error("Dimension of element mismatch!");
226 
227  unsigned int neigh_idx = ngh_idx(cell_lower_dim, neighb_side); // TODO use better evaluation of neighbour_idx
228  unsigned int loc_dof_higher = (2*(cell_lower_dim.dim()+1) + 1) + neigh_idx; // loc dof of higher ele edge
229  bulk_local_idx_ = cell_lower_dim.local_idx();
230 
231  // Evaluation points
232  auto p_high = *( this->coupling_points(neighb_side).begin() );
233  auto p_low = p_high.lower_dim(cell_lower_dim);
234 
235  fe_values_side_.reinit(neighb_side.side());
237 
238  ngh_value_ = eq_fields_->sigma(p_low) *
239  2*eq_fields_->conductivity(p_low) *
240  arma::dot(eq_fields_->anisotropy(p_low)*nv_, nv_) *
241  eq_fields_->cross_section(p_high) * // cross-section of higher dim. (2d)
242  eq_fields_->cross_section(p_high) /
243  eq_fields_->cross_section(p_low) * // crossection of lower dim.
244  neighb_side.measure();
245 
247  eq_data_->loc_system_[bulk_local_idx_].add_value(eq_data_->loc_ele_dof[dim-2], loc_dof_higher, ngh_value_);
248  eq_data_->loc_system_[bulk_local_idx_].add_value(loc_dof_higher, eq_data_->loc_ele_dof[dim-2], ngh_value_);
249  eq_data_->loc_system_[bulk_local_idx_].add_value(loc_dof_higher, loc_dof_higher, -ngh_value_);
250 
251 // // update matrix for weights in BDDCML
252 // if ( typeid(*eq_data_->lin_sys) == typeid(LinSys_BDDC) ) {
253 // int ind = eq_data_->loc_system_[bulk_local_idx_].row_dofs[p];
254 // // there is -value on diagonal in block C!
255 // static_cast<LinSys_BDDC*>(eq_data_->lin_sys)->diagonal_weights_set_value( ind, -value );
256 // }
257  }
258 
259 
260  /// Implements @p AssemblyBase::begin.
261  void begin() override
262  {
263  this->begin_mh_matrix();
264  }
265 
266 
267  /// Implements @p AssemblyBase::end.
268  void end() override
269  {
270  this->end_mh_matrix();
271  }
272 
273 protected:
274  /// Common code of begin method of MH matrix assembly (Darcy and Richards)
276  {
277  // DebugOut() << "assembly_mh_matrix \n";
278  // set auxiliary flag for switchting Dirichlet like BC
280 
281  eq_data_->reset();
282 
283  eq_data_->balance_->start_flux_assembly(eq_data_->water_balance_idx);
284  eq_data_->balance_->start_source_assembly(eq_data_->water_balance_idx);
285  eq_data_->balance_->start_mass_assembly(eq_data_->water_balance_idx);
286 
287  this->set_dofs();
288  }
289 
290 
291  /// Common code of end method of MH matrix assembly (Darcy and Richards)
293  {
294  std::sort(eq_data_->loc_constraint_.begin(), eq_data_->loc_constraint_.end());
295  eq_data_->loc_constraint_.emplace_back(uint(-1), 0, 0.0); // add constraint of invalid element to end of vector
296  unsigned int i_constr = 0;
297  for ( DHCellAccessor dh_cell : eq_data_->dh_cr_->own_range() ) {
298  this->set_loc_schur(dh_cell);
299  while (eq_data_->loc_constraint_[i_constr].i_element() == dh_cell.local_idx()) {
301  i_constr++;
302  }
303  eq_data_->loc_system_[dh_cell.local_idx()].compute_schur_complement(
304  eq_data_->schur_offset_[dh_cell.dim()-1], this->loc_schur_, true);
305 
306  // for seepage BC, save local system
307  if (eq_data_->save_local_system_[dh_cell.local_idx()])
308  eq_data_->seepage_bc_systems[dh_cell.elm_idx()] = eq_data_->loc_system_[dh_cell.local_idx()];
309 
311  eq_data_->lin_sys_schur->set_local_system(this->loc_schur_, eq_data_->dh_cr_->get_local_to_global_map());
312 
313  // TODO:
314  // if (mortar_assembly)
315  // mortar_assembly->assembly(dh_cell);
316  }
317 
318  eq_data_->balance_->finish_mass_assembly(eq_data_->water_balance_idx);
319  eq_data_->balance_->finish_source_assembly(eq_data_->water_balance_idx);
320  eq_data_->balance_->finish_flux_assembly(eq_data_->water_balance_idx);
321 
322  eq_data_->loc_constraint_.clear();
323  }
324 
325 
326  /// Common code of begin method of Reconstruct Schur assembly (Darcy and Richards)
328  {
332  }
333 
334 
335  /// Common code of end method of Reconstruct Schur assembly (Darcy and Richards)
337  {
338  for ( DHCellAccessor dh_cell : this->eq_data_->dh_cr_->own_range() ) {
339  this->bulk_local_idx_ = dh_cell.local_idx();
340  this->set_loc_schur(dh_cell);
341  arma::vec schur_solution = this->eq_data_->p_edge_solution.get_subvec(this->loc_schur_.row_dofs);
342  // reconstruct the velocity and pressure
343  this->eq_data_->loc_system_[this->bulk_local_idx_].reconstruct_solution_schur(this->eq_data_->schur_offset_[dh_cell.dim()-1],
344  schur_solution, this->reconstructed_solution_);
345 
347 
348  this->eq_data_->full_solution.set_subvec(this->eq_data_->loc_system_[this->bulk_local_idx_].row_dofs.head(
349  this->eq_data_->schur_offset_[dh_cell.dim()-1]), this->reconstructed_solution_);
350  this->eq_data_->full_solution.set_subvec(this->eq_data_->loc_system_[this->bulk_local_idx_].row_dofs.tail(
351  this->loc_schur_.row_dofs.n_elem), schur_solution);
352  }
353 
356 
357  eq_data_->loc_constraint_.clear();
358  }
359 
360  /// Part of cell_integral method, common in all descendants
361  inline void asm_sides(const DHCellAccessor& cell, BulkPoint &p, double conductivity)
362  {
363  auto ele = cell.elm();
364  double scale_sides = 1 / eq_fields_->cross_section(p) / conductivity;
365 
366  fe_values_.reinit(ele);
367  auto velocity = fe_values_.vector_view(0);
368 
369  for (unsigned int k=0; k<fe_values_.n_points(); k++)
370  for (unsigned int i=0; i<fe_values_.n_dofs(); i++){
371  double rhs_val = arma::dot(eq_data_->gravity_vec_, velocity.value(i,k))
372  * fe_values_.JxW(k);
373  eq_data_->loc_system_[bulk_local_idx_].add_value(i, rhs_val);
374 
375  for (unsigned int j=0; j<fe_values_.n_dofs(); j++){
376  double mat_val =
377  arma::dot( velocity.value(i,k), //TODO: compute anisotropy before
378  (eq_fields_->anisotropy(p)).i() * velocity.value(j,k)
379  )
380  * scale_sides * fe_values_.JxW(k);
381 
382  eq_data_->loc_system_[bulk_local_idx_].add_value(i, j, mat_val);
383  }
384  }
385 
386  // assemble matrix for weights in BDDCML
387  // approximation to diagonal of
388  // S = -C - B*inv(A)*B'
389  // as
390  // diag(S) ~ - diag(C) - 1./diag(A)
391  // the weights form a partition of unity to average a discontinuous solution from neighbouring subdomains
392  // to a continuous one
393  // it is important to scale the effect - if conductivity is low for one subdomain and high for the other,
394  // trust more the one with low conductivity - it will be closer to the truth than an arithmetic average
395 // if ( typeid(*ad_->lin_sys) == typeid(LinSys_BDDC) ) {
396 // const arma::mat& local_matrix = loc_system_.get_matrix();
397 // for(unsigned int i=0; i < ndofs; i++) {
398 // double val_side = local_matrix(i,i);
399 // double val_edge = -1./local_matrix(i,i);
400 //
401 // unsigned int side_row = loc_system_.row_dofs[loc_side_dofs[i]];
402 // unsigned int edge_row = loc_system_.row_dofs[loc_edge_dofs[i]];
403 // static_cast<LinSys_BDDC*>(ad_->lin_sys)->diagonal_weights_set_value( side_row, val_side );
404 // static_cast<LinSys_BDDC*>(ad_->lin_sys)->diagonal_weights_set_value( edge_row, val_edge );
405 // }
406 // }
407  }
408 
409  /// Part of cell_integral method, common in all descendants
410  inline void asm_element()
411  {
412  // set block B, B': element-side, side-element
413 
414  for(unsigned int side = 0; side < eq_data_->loc_side_dofs[dim-1].size(); side++){
415  eq_data_->loc_system_[bulk_local_idx_].add_value(eq_data_->loc_ele_dof[dim-1], eq_data_->loc_side_dofs[dim-1][side], -1.0);
416  eq_data_->loc_system_[bulk_local_idx_].add_value(eq_data_->loc_side_dofs[dim-1][side], eq_data_->loc_ele_dof[dim-1], -1.0);
417  }
418 
419 // if ( typeid(*ad_->lin_sys) == typeid(LinSys_BDDC) ) {
420 // double val_ele = 1.;
421 // static_cast<LinSys_BDDC*>(ad_->lin_sys)->
422 // diagonal_weights_set_value( loc_system_.row_dofs[loc_ele_dof], val_ele );
423 // }
424  }
425 
426  /// Part of cell_integral method, specialized in Darcy equation
427  inline void asm_source_term_darcy(const DHCellAccessor& cell, BulkPoint &p)
428  {
429  // compute lumped source
430  uint n_sides = cell.elm()->n_sides();
431  coef_ = (1.0 / n_sides) * cell.elm().measure() * eq_fields_->cross_section(p);
433 
434  // in unsteady, compute time term
435  time_term_diag_ = 0.0;
436  time_term_ = 0.0;
437  time_term_rhs_ = 0.0;
438 
440  {
442  }
443 
444  const DHCellAccessor cr_cell = cell.cell_with_other_dh(eq_data_->dh_cr_.get());
445  auto loc_dof_vec = cr_cell.get_loc_dof_indices();
446 
447  for (unsigned int i=0; i<n_sides; i++)
448  {
450  {
453 
454  eq_data_->balance_->add_mass_values(eq_data_->water_balance_idx, cell,
455  {eq_data_->loc_system_[bulk_local_idx_].row_dofs[eq_data_->loc_edge_dofs[dim-1][i]]}, {time_term_}, 0);
456  }
457  else
458  {
459  // Add zeros explicitely to keep the sparsity pattern.
460  // Otherwise Petsc would compress out the zeros in FinishAssembly.
461  eq_data_->balance_->add_mass_values(eq_data_->water_balance_idx, cell,
462  {eq_data_->loc_system_[bulk_local_idx_].row_dofs[eq_data_->loc_edge_dofs[dim-1][i]]}, {0}, 0);
463  }
464 
468 
469  eq_data_->balance_->add_source_values(eq_data_->water_balance_idx, cell.elm().region().bulk_idx(),
470  {eq_data_->loc_system_[bulk_local_idx_].row_dofs[eq_data_->loc_edge_dofs[dim-1][i]]}, {0}, {source_term_});
471  }
472  }
473 
474  /// Precompute values used in boundary side integral on given DHCellSide
475  inline void precompute_boundary_side(DHCellSide &cell_side, BoundaryPoint &p_side, BulkPoint &p_bdr)
476  {
477  bulk_local_idx_ = cell_side.cell().local_idx();
478 
479  cross_section_ = eq_fields_->cross_section(p_side);
480  type_ = (DarcyLMH::EqFields::BC_Type)eq_fields_->bc_type(p_bdr);
481 
482  sidx_ = cell_side.side_idx();
483  side_row_ = eq_data_->loc_side_dofs[dim-1][sidx_]; //local
484  edge_row_ = eq_data_->loc_edge_dofs[dim-1][sidx_]; //local
485  }
486 
487  /// Update BC switch dirichlet in MH matrix assembly if BC type is seepage
488  inline void use_dirichlet_switch(DHCellSide &cell_side, const ElementAccessor<3> &b_ele, BulkPoint &p_bdr)
489  {
490  char & switch_dirichlet = eq_data_->bc_switch_dirichlet[b_ele.idx()];
491  if (switch_dirichlet) {
492  // check and possibly switch to flux BC
493  // The switch raise error on the corresponding edge row.
494  // Magnitude of the error is abs(solution_flux - side_flux).
495 
496  // try reconstructing local system for seepage BC
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_;
499 
500  if ( reconstr_solution[side_row_] < side_flux) {
501  //DebugOut().fmt("x: {}, to neum, p: {} f: {} -> f: {}\n", b_ele.centre()[0], bc_pressure, reconstructed_solution_[side_row_], side_flux);
502  switch_dirichlet = 0;
503  }
504  } else {
505  // check and possibly switch to pressure BC
506  // TODO: What is the appropriate DOF in not local?
507  // The switch raise error on the corresponding side row.
508  // Magnitude of the error is abs(solution_head - bc_pressure)
509  // Since usually K is very large, this error would be much
510  // higher then error caused by the inverse switch, this
511  // cause that a solution with the flux violating the
512  // flux inequality leading may be accepted, while the error
513  // in pressure inequality is always satisfied.
514 
515  const DHCellAccessor cr_cell = cell_side.cell().cell_with_other_dh(eq_data_->dh_cr_.get());
516  auto loc_dof_vec = cr_cell.get_loc_dof_indices();
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);
519 
520  if ( solution_head > bc_pressure) {
521  //DebugOut().fmt("x: {}, to dirich, p: {} -> p: {} f: {}\n",b_ele.centre()[0], solution_head, bc_pressure, bc_flux);
522  switch_dirichlet=1;
523  }
524  }
525  }
526 
527  /**
528  * Common code of boundary_side_integral.
529  */
530  inline void boundary_side_integral_in(DHCellSide cell_side, const ElementAccessor<3> &b_ele, BulkPoint &p_bdr)
531  {
532  if ( type_ == DarcyLMH::EqFields::none) {
533  // homogeneous neumann
534  } else if ( type_ == DarcyLMH::EqFields::dirichlet ) {
535  eq_data_->loc_constraint_.emplace_back( bulk_local_idx_, sidx_, eq_fields_->bc_pressure(p_bdr) );
536 
537  } else if ( type_ == DarcyLMH::EqFields::total_flux) {
538  // internally we work with outward flux
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_);
542  } else if (type_==DarcyLMH::EqFields::seepage) {
543  eq_data_->is_linear=false;
544 
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_;
547 
548  eq_data_->save_local_system_[bulk_local_idx_] = (bool) eq_data_->bc_switch_dirichlet[b_ele.idx()];
549 
550  // ** Apply BCUpdate BC type. **
551  // Force Dirichlet type during the first iteration of the unsteady case.
552  if (eq_data_->save_local_system_[bulk_local_idx_] || eq_data_->force_no_neumann_bc ) {
553  //DebugOut().fmt("x: {}, dirich, bcp: {}\n", b_ele.centre()[0], bc_pressure);
554  eq_data_->loc_constraint_.emplace_back(bulk_local_idx_, sidx_, bc_pressure);
555  } else {
556  //DebugOut()("x: {}, neuman, q: {} bcq: {}\n", b_ele.centre()[0], side_flux, bc_flux);
557  eq_data_->loc_system_[bulk_local_idx_].add_value(edge_row_, side_flux);
558  }
559 
560  } else if (type_==DarcyLMH::EqFields::river) {
561  eq_data_->is_linear=false;
562 
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);
567 
568  const DHCellAccessor cr_cell = cell_side.cell().cell_with_other_dh(eq_data_->dh_cr_.get());
569  auto loc_dof_vec = cr_cell.get_loc_dof_indices();
570  double solution_head = eq_data_->p_edge_solution.get(loc_dof_vec[sidx_]);
571 
572  // Force Robin type during the first iteration of the unsteady case.
573  if (solution_head > bc_switch_pressure || eq_data_->force_no_neumann_bc) {
574  // Robin BC
575  //DebugOut().fmt("x: {}, robin, bcp: {}\n", b_ele.centre()[0], bc_pressure);
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) );
579  } else {
580  // Neumann BC
581  //DebugOut().fmt("x: {}, neuman, q: {} bcq: {}\n", b_ele.centre()[0], bc_switch_pressure, bc_pressure);
582  double bc_total_flux = bc_flux + bc_sigma*(bc_switch_pressure - bc_pressure);
583 
584  eq_data_->loc_system_[bulk_local_idx_].add_value(edge_row_, bc_total_flux * b_ele.measure() * cross_section_);
585  }
586  }
587  else {
588  THROW( ExcBCNotSupported() );
589  }
590 
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_]]},
593  {1}, 0);
594  }
595 
596  /// Precompute loc_system and loc_schur data members.
597  void set_dofs() {
598  unsigned int size, loc_size, elm_dim;;
599  for ( DHCellAccessor dh_cell : eq_data_->dh_->local_range() ) {
600  const ElementAccessor<3> ele = dh_cell.elm();
601 
602  elm_dim = ele.dim();
603  size = (elm_dim+1) + 1 + (elm_dim+1); // = n_sides + 1 + n_sides
604  loc_size = size + ele->n_neighs_vb();
605  LocDofVec dofs(loc_size);
606 
607  // add full vec indices
608  dofs.head(dh_cell.n_dofs()) = dh_cell.get_loc_dof_indices();
609 
610  if(ele->n_neighs_vb() != 0)
611  {
612  //D, E',E block: compatible connections: element-edge
613  unsigned int i = 0;
614 
615  for ( DHCellSide neighb_side : dh_cell.neighb_sides() ) {
616 
617  // read neighbor dofs (dh dofhandler)
618  // neighbor cell owning neighb_side
619  DHCellAccessor dh_neighb_cell = neighb_side.cell();
620 
621  // read neighbor dofs (dh_cr dofhandler)
622  // neighbor cell owning neighb_side
623  DHCellAccessor dh_neighb_cr_cell = dh_neighb_cell.cell_with_other_dh(eq_data_->dh_cr_.get());
624 
625  // local index of pedge dof in local system
626  const unsigned int p = size+i;
627  // local index of pedge dof on neighboring cell
628  const unsigned int t = dh_neighb_cell.n_dofs() - dh_neighb_cr_cell.n_dofs() + neighb_side.side().side_idx();
629  dofs[p] = dh_neighb_cell.get_loc_dof_indices()[t];
630  i++;
631  }
632  }
633  eq_data_->loc_system_[dh_cell.local_idx()].reset(dofs, dofs);
634 
635  // Set side-edge (flux-lambda) terms
636  for (DHCellSide dh_side : dh_cell.side_range()) {
637  unsigned int sidx = dh_side.side_idx();
638  // side-edge (flux-lambda) terms
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);
641  }
642  }
643  }
644 
645 
646  /// Precompute loc_schur data member of given cell.
647  void set_loc_schur(const DHCellAccessor dh_cr_cell) {
648  const ElementAccessor<3> ele = dh_cr_cell.elm();
649 
650  unsigned int loc_size_schur = ele->n_sides() + ele->n_neighs_vb();
651  LocDofVec dofs_schur(loc_size_schur);
652 
653  // add schur vec indices
654  dofs_schur.head(dh_cr_cell.n_dofs()) = dh_cr_cell.get_loc_dof_indices();
655 
656  if(ele->n_neighs_vb() != 0)
657  {
658  //D, E',E block: compatible connections: element-edge
659  unsigned int i = 0;
660 
661  for ( DHCellSide neighb_side : dh_cr_cell.neighb_sides() ) {
662 
663  // read neighbor dofs (dh dofhandler)
664  // neighbor cell owning neighb_side
665  DHCellAccessor dh_neighb_cr_cell = neighb_side.cell();
666 
667  // local index of pedge dof in local schur system
668  const unsigned int tt = dh_cr_cell.n_dofs()+i;
669  dofs_schur[tt] = dh_neighb_cr_cell.get_loc_dof_indices()[neighb_side.side().side_idx()];
670  i++;
671  }
672  }
673 
674 
675  loc_schur_.reset(dofs_schur, dofs_schur);
676  }
677 
678 
679  /// Temporary method find neighbour index in higher-dim cell
680  inline unsigned int ngh_idx(DHCellAccessor &dh_cell, DHCellSide &neighb_side) {
681  for (uint n_i=0; n_i<dh_cell.elm()->n_neighs_vb(); ++n_i) {
682  auto side = dh_cell.elm()->neigh_vb[n_i]->side();
683  if ( (side->elem_idx() == neighb_side.elem_idx()) && (side->side_idx() == neighb_side.side_idx()) ) return n_i;
684  }
685  ASSERT_PERMANENT(false)(dh_cell.elm_idx())(neighb_side.side_idx()).error("Side is not a part of neighbour!\n");
686  return 0;
687  }
688 
689 
690  /** Loads the local system from a map: element index -> LocalSystem,
691  * if it exits, or if the full solution is not yet reconstructed,
692  * and reconstructs the full solution on the element.
693  * Currently used only for seepage BC.
694  */
696  {
697  // possibly find the corresponding local system
698  auto ls = eq_data_->seepage_bc_systems.find(dh_cell.elm_idx());
699  if (ls != eq_data_->seepage_bc_systems.end())
700  {
701  const DHCellAccessor cr_cell = dh_cell.cell_with_other_dh(eq_data_->dh_cr_.get());
702  auto loc_dof_vec = cr_cell.get_loc_dof_indices();
703  arma::vec schur_solution = eq_data_->p_edge_solution.get_subvec(loc_dof_vec);
704  // reconstruct the velocity and pressure
705  ls->second.reconstruct_solution_schur(eq_data_->schur_offset_[dim-1], schur_solution, reconstructed_solution_);
706 
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_);
710 
711  eq_data_->bc_fluxes_reconstruted[bulk_local_idx_] = true;
712  }
713 
714  return reconstructed_solution_;
715  }
716 
717 
718  /**
719  * Precompute edge_scale and edge_source_term.
720  *
721  * This method must be calls in methods postprocess_velocity_darcy and postprocess_velocity_richards
722  */
724  {
725  edge_scale_ = dh_cell.elm().measure()
726  * eq_fields_->cross_section(p)
727  / dh_cell.elm()->n_sides();
728 
729  edge_source_term_ = edge_scale_ *
730  ( eq_fields_->water_source_density(p)
731  + eq_fields_->extra_source(p));
732  }
733 
734 
735  /// Postprocess velocity during loading of local system and after calculating of cell integral.
736  void postprocess_velocity_darcy(const DHCellAccessor& dh_cell, BulkPoint &p, arma::vec& solution)
737  {
738  postprocess_velocity(dh_cell, p);
739 
740  time_term_ = 0.0;
741  const DHCellAccessor cr_cell = dh_cell.cell_with_other_dh(eq_data_->dh_cr_.get());
742  auto loc_dof_vec = cr_cell.get_loc_dof_indices();
743 
744  for (unsigned int i=0; i<dh_cell.elm()->n_sides(); i++) {
745 
746  if( ! eq_data_->use_steady_assembly_)
747  {
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);
752  }
753  solution[eq_data_->loc_side_dofs[dim-1][i]] += edge_source_term_ - time_term_;
754  }
755  }
756 
757 
758  /// Sub field set contains fields used in calculation.
760 
761  /// Data objects shared with DarcyFlow
764 
765  /// Assembly volume integrals
769 
770  shared_ptr<FiniteElement<dim>> fe_; ///< Finite element for the solution of the advection-diffusion equation.
771  FEValues<3> fe_values_side_; ///< FEValues of object (of P disc finite element type)
772 
773  /// Vector for reconstructed solution (velocity and pressure on element) from Schur complement.
775 
776  double coef_, source_term_; ///< Variables used in compute lumped source.
777  double time_term_diag_, time_term_, time_term_rhs_; ///< Variables used in compute time term (unsteady)
778  double cross_section_; ///< Precomputed cross-section value
779  unsigned int bulk_local_idx_; ///< Local idx of bulk element
780  unsigned int sidx_, side_row_, edge_row_; ///< Helper indices in boundary assembly
781  DarcyLMH::EqFields::BC_Type type_; ///< Type of boundary condition
782  arma::vec3 nv_; ///< Normal vector
783  double ngh_value_; ///< Precomputed ngh value
784  double edge_scale_, edge_source_term_; ///< Precomputed values in postprocess_velocity
785 
787 
788  template < template<IntDim...> class DimAssembly>
789  friend class GenericAssembly;
790 
791 };
792 
793 template <unsigned int dim>
795 {
796 public:
797  typedef typename DarcyLMH::EqFields EqFields;
798  typedef typename DarcyLMH::EqData EqData;
799 
800  static constexpr const char * name() { return "ReconstructSchurAssemblyLMH"; }
801 
802  /// Constructor.
804  : MHMatrixAssemblyLMH<dim>(eq_fields, eq_data) {}
805 
806  /// Integral over element.
807  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
808  {
809  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
810 
811  // evaluation point
812  auto p = *( this->bulk_points(element_patch_idx).begin() );
813  this->bulk_local_idx_ = cell.local_idx();
814 
815  { // postprocess the velocity
816  this->eq_data_->postprocess_solution_[this->bulk_local_idx_].zeros(this->eq_data_->schur_offset_[dim-1]);
817  this->postprocess_velocity_darcy(cell, p, this->eq_data_->postprocess_solution_[this->bulk_local_idx_]);
818  }
819  }
820 
821 
822  /// Assembles between boundary element and corresponding side on bulk element.
824  {}
825 
826  inline void dimjoin_intergral(FMT_UNUSED DHCellAccessor cell_lower_dim, FMT_UNUSED DHCellSide neighb_side)
827  {}
828 
829 
830  /// Implements @p AssemblyBase::begin.
831  void begin() override
832  {
833  this->begin_reconstruct_schur();
834  }
835 
836 
837  /// Implements @p AssemblyBase::end.
838  void end() override
839  {
840  this->end_reconstruct_schur();
841  }
842 
843 };
844 
845 #endif /* ASSEMBLY_LMH_HH_ */
846 
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
#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.
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...
Definition: eval_subset.hh:209
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_.
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.
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.
double measure() const
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_
arma::vec3 gravity_vec_
std::map< LongIdx, LocalSystem > seepage_bc_systems
VectorMPI full_solution
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_
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.
Definition: accessors.hh:329
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...
Definition: accessors.hh:215
unsigned int dim() const
Definition: accessors.hh:188
const Element * element() const
Definition: accessors.hh:193
Directing class of FieldValueCache.
Neighbour ** neigh_vb
Definition: elements.h:83
unsigned int n_neighs_vb() const
Return number of neighbours.
Definition: elements.h:65
unsigned int n_sides() const
Definition: elements.h:131
Compound finite element on dim dimensional simplex.
Definition: fe_system.hh:102
std::vector< unsigned int > fe_dofs(unsigned int fe_index)
Return dof indices belonging to given sub-FE.
Definition: fe_system.cc:267
double JxW(const unsigned int point_no)
Return the product of Jacobian determinant and the quadrature weight at given quadrature point.
Definition: fe_values.hh:223
const FEValuesViews::Vector< spacedim > & vector_view(unsigned int i) const
Accessor to vector values of multicomponent FE.
Definition: fe_values.hh:277
void initialize(Quadrature &_quadrature, FiniteElement< DIM > &_fe, UpdateFlags _flags)
Initialize structures and calculates cell-independent data.
Definition: fe_values.cc:137
void reinit(const ElementAccessor< spacedim > &cell)
Update cell-dependent data (gradients, Jacobians etc.)
Definition: fe_values.cc:537
unsigned int n_points() const
Returns the number of quadrature points.
Definition: fe_values.hh:296
arma::vec::fixed< spacedim > normal_vector(unsigned int point_no)
Returns the normal vector to a side at given quadrature point.
Definition: fe_values.hh:257
unsigned int n_dofs() const
Returns the number of shape functions.
Definition: fe_values.hh:302
Raviart-Thomas element of order 0.
Definition: fe_rt.hh:61
Container for various descendants of FieldCommonBase.
Definition: field_set.hh:159
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.
Definition: local_system.cc:82
void eliminate_solution()
LocDofVec row_dofs
Definition: local_system.hh:54
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")
FEValues< 3 > fe_values_
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)
LocalSystem loc_schur_
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::EqData EqData
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)
SideIter side() const
Definition: neighbours.h:147
Symmetric Gauss-Legendre quadrature formulae on simplices.
FieldSet used_fields_
Sub field set contains fields used in calculation.
static constexpr const char * name()
Definition: assembly_lmh.hh:52
DarcyLMH::EqData EqData
Definition: assembly_lmh.hh:50
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
Definition: assembly_lmh.hh:65
ReadInitCondAssemblyLMH(EqFields *eq_fields, EqData *eq_data)
Constructor.
Definition: assembly_lmh.hh:55
void end() override
Implements AssemblyBase::end.
Definition: assembly_lmh.hh:89
DarcyLMH::EqFields EqFields
Definition: assembly_lmh.hh:49
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
Definition: assembly_lmh.hh:71
LocDofVec l_indices_
Vector of pre-computed local DOF indices.
virtual ~ReadInitCondAssemblyLMH()
Destructor.
Definition: assembly_lmh.hh:62
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
Boundary cond() const
Edge edge() const
Returns pointer to the edge connected to the side.
arma::vec get_subvec(const LocDofVec &loc_indices)
Definition: vector_mpi.cc:110
double get(unsigned int pos) const
Return value on given position.
Definition: vector_mpi.hh:105
void local_to_ghost_end()
local_to_ghost_{begin,end} updates the ghost values on neighbouring processors from local values
Definition: vector_mpi.hh:160
void set_subvec(const LocDofVec &loc_indices, const arma::vec &values)
Definition: vector_mpi.cc:138
void add(unsigned int pos, double val)
Add value to item on given position.
Definition: vector_mpi.hh:125
void local_to_ghost_begin()
local_to_ghost_{begin,end} updates the ghost values on neighbouring processors from local values
Definition: vector_mpi.hh:156
void zero_entries()
Definition: vector_mpi.hh:82
void ghost_to_local_end()
ghost_to_local_{begin,end} updates the local values by adding ghost values from neighbouring processo...
Definition: vector_mpi.hh:168
void ghost_to_local_begin()
ghost_to_local_{begin,end} updates the local values by adding ghost values from neighbouring processo...
Definition: vector_mpi.hh:164
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,...
@ coupling
@ boundary
@ bulk
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
arma::Col< IntIdx > LocDofVec
Definition: index_types.hh:28
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
unsigned int uint
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19
ArmaVec< double, N > vec
Definition: armor.hh:933
#define FMT_UNUSED
Definition: posix.h:75
Definitions of particular quadrature rules on simplices.
Assembly explicit Schur complement for the given linear system. Provides method for resolution of the...
Definition: mixed.hh:25
@ update_values
Shape function values.
Definition: update_flags.hh:87
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_quadrature_points
Transformed quadrature points.