Flow123d  master-f44eb46
assembly_mh_old.hh
Go to the documentation of this file.
1 /*
2  * assembly_mh_old.hh
3  *
4  * Created on: Apr 21, 2016
5  * Author: jb
6  */
7 
8 #ifndef SRC_FLOW_DARCY_FLOW_ASSEMBLY_HH_
9 #define SRC_FLOW_DARCY_FLOW_ASSEMBLY_HH_
10 
11 #include "system/index_types.hh"
12 #include "mesh/mesh.h"
13 #include "mesh/accessors.hh"
14 #include "mesh/neighbours.h"
15 #include "fem/fe_p.hh"
16 #include "fem/fe_values.hh"
17 #include "fem/fe_rt.hh"
18 #include "fem/fe_values_views.hh"
20 
21 #include "la/linsys.hh"
22 #include "la/linsys_PETSC.hh"
23 #include "la/linsys_BDDC.hh"
24 #include "la/schur.hh"
25 
26 #include "la/local_system.hh"
27 
28 #include "coupling/balance.hh"
29 #include "flow/darcy_flow_mh.hh"
30 #include "flow/mortar_assembly.hh"
31 
32 
33 /** Common abstract class for the assembly routines in Darcy flow.
34  * Is implemented in DarcyMH, DarcyLMH and RichardsLMH assembly classes,
35  * which are independent of each other.
36  */
38 {
39 public:
40  DECLARE_EXCEPTION( ExcBCNotSupported, << "BC type not supported.\n" );
41 
42  virtual void fix_velocity(const DHCellAccessor& dh_cell) = 0;
43  virtual void assemble(const DHCellAccessor& dh_cell) = 0;
44  virtual void assemble_reconstruct(const DHCellAccessor& dh_cell) = 0;
45 
46  /// Updates water content in Richards.
47  virtual void update_water_content(const DHCellAccessor& dh_cell) = 0;
48 
49  /**
50  * Generic creator of multidimensional assembly, i.e. vector of
51  * particular assembly objects.
52  */
53  template< template<int dim> class Impl, class Fields, class Data>
54  static MultidimAssembly create(Fields eq_fields, Data eq_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) };
58  }
59 
60  virtual ~AssemblyFlowBase() {}
61 };
62 
63 
64 template <int dim>
66 private:
67  // assembly face integrals (BC)
70 public:
72  : side_quad_(dim, 1),
73  fe_p_disc_(0)
74  {
76  }
78 
79 };
80 
81 
82 /** MH version of Darcy flow assembly. It is supposed not to be improved anymore,
83  * however it is kept functioning aside of the LMH lumped version until
84  * the LMH version is stable and optimized.
85  */
86 template<int dim>
88 {
89 public:
90  typedef std::shared_ptr<DarcyMH::EqFields> AssemblyFieldsPtrMH;
91  typedef std::shared_ptr<DarcyMH::EqData> AssemblyDataPtrMH;
92 
94  : quad_(dim, 2),
95  velocity_interpolation_quad_(dim, 0), // veloctiy values in barycenter
96  af_(eq_fields),
97  ad_(eq_data),
98  loc_system_(size(), size()),
99  loc_system_vb_(2,2)
100 
101  {
105 
106  // local numbering of dofs for MH system
107  unsigned int nsides = dim+1;
108  loc_side_dofs.resize(nsides);
109  loc_ele_dof = nsides;
110  loc_edge_dofs.resize(nsides);
111  for(unsigned int i = 0; i < nsides; i++){
112  loc_side_dofs[i] = i;
113  loc_edge_dofs[i] = nsides + i + 1;
114  }
115  //DebugOut() << print_var(this) << print_var(side_quad_.size());
116 
117  // create local sparsity pattern
118  arma::umat sp(size(),size());
119  sp.zeros();
120  sp.submat(0, 0, nsides, nsides).ones();
121  sp.diag().ones();
122  // armadillo 8.4.3 bug with negative sub diagonal index
123  // sp.diag(nsides+1).ones();
124  // sp.diag(-nsides-1).ones();
125  // sp.print();
126 
127  sp.submat(0, nsides+1, nsides-1, size()-1).diag().ones();
128  sp.submat(nsides+1, 0, size()-1, nsides-1).diag().ones();
129 
131 
132  // local system 2x2 for vb neighbourings is full matrix
133  // this matrix cannot be influenced by any BC (no elimination can take place)
134  sp.ones(2,2);
136 
137  if (ad_->mortar_method_ == DarcyFlowInterface::MortarP0) {
138  mortar_assembly = std::make_shared<P0_CouplingAssembler>(af_, ad_);
139  } else if (ad_->mortar_method_ == DarcyFlowInterface::MortarP1) {
140  mortar_assembly = std::make_shared<P1_CouplingAssembler>(af_, ad_);
141  }
142 
143  }
144 
145  void assemble_reconstruct(const DHCellAccessor& dh_cell) override
146  {
147  ASSERT_EQ(dh_cell.dim(), dim);
148  loc_system_.reset();
149 
150  // assemble_local_system(dh_cell); //do not switch dirichlet in seepage when reconstructing
151  set_dofs_and_bc(dh_cell, false);
152  // assemble_bc(dh_cell, use_dirichlet_switch);
153  assemble_sides(dh_cell);
154  assemble_element(dh_cell);
155  // assemble_source_term(dh_cell);
156  // assembly_dim_connections(dh_cell);
157 
158  // TODO:
159  // if (mortar_assembly)
160  // mortar_assembly->assembly(ele_ac);
161  // if (mortar_assembly)
162  // mortar_assembly->fix_velocity(ele_ac);
163 
164  // arma::vec schur_solution = ad_->p_edge_solution.get_subvec(loc_schur_.row_dofs);
165  unsigned int nsides = dim+1;
166  unsigned int schur_offset = loc_edge_dofs[0];
167  LocDofVec loc_dofs = dh_cell.get_loc_dof_indices();
168  // get edge pressure
169  arma::vec schur_solution = ad_->full_solution.get_subvec(loc_dofs.subvec(nsides+1, size()-1));
170 
171  // prepare vector for velocity and ele pressure
172  arma::vec reconstructed_solution;
173  reconstructed_solution.zeros(schur_offset);
174 
175  // reconstruct the velocity and pressure
176  loc_system_.reconstruct_solution_schur(schur_offset, schur_solution, reconstructed_solution);
177 
178  // postprocess the velocity
179  // postprocess_velocity(dh_cell, reconstructed_solution);
180 
181  ad_->full_solution.set_subvec(loc_dofs.head(schur_offset), reconstructed_solution);
182  // ad_->full_solution.set_subvec(loc_system_.row_dofs.tail(loc_schur_.row_dofs.n_elem), schur_solution);
183  }
184 
185  void update_water_content(const DHCellAccessor&) override
186  {};
187 
188  ~AssemblyMH<dim>() override
189  {}
190 
191 // LocalSystem& get_local_system() override
192 // { return loc_system_;}
193 
194  void fix_velocity(const DHCellAccessor& dh_cell) override
195  {
196  if (mortar_assembly)
197  mortar_assembly->fix_velocity(dh_cell);
198  }
199 
200  void assemble(const DHCellAccessor& dh_cell) override
201  {
202  ASSERT_EQ(dh_cell.dim(), dim);
203  loc_system_.reset();
204 
205  set_dofs_and_bc(dh_cell, true);
206 
207  assemble_sides(dh_cell);
208  assemble_element(dh_cell);
209 
211  ad_->lin_sys->set_local_system(loc_system_);
212 
213  assembly_dim_connections(dh_cell);
214 
215  if (ad_->balance != nullptr)
217 
218  if (mortar_assembly)
219  mortar_assembly->assembly(dh_cell);
220  }
221 
222  void assembly_local_vb(ElementAccessor<3> ele, DHCellSide neighb_side) //override
223  {
224  ASSERT_LT(ele->dim(), 3);
225  //DebugOut() << "alv " << print_var(this);
226  //START_TIMER("Assembly<dim>::assembly_local_vb");
227  // compute normal vector to side
228  arma::vec3 nv;
229  ElementAccessor<3> ele_higher = ad_->mesh->element_accessor( neighb_side.element().idx() );
230  ngh_values_.fe_side_values_.reinit(neighb_side.side());
231  nv = ngh_values_.fe_side_values_.normal_vector(0);
232 
233  double value = af_->sigma.value( ele.centre(), ele) *
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 ) * // cross-section of higher dim. (2d)
237  af_->cross_section.value( neighb_side.centre(), ele_higher ) /
238  af_->cross_section.value( ele.centre(), ele ) * // crossection of lower dim.
239  neighb_side.measure();
240 
245  }
246 
247 protected:
248  static unsigned int size()
249  {
250  // dofs: velocity, pressure, edge pressure
252  }
253 
254  void set_dofs_and_bc(const DHCellAccessor& dh_cell, bool use_dirichlet_switch){
255 
256  local_dofs_ = dh_cell.get_loc_dof_indices();
257  global_dofs_.resize(dh_cell.n_dofs());
258  dh_cell.get_dof_indices(global_dofs_);
259 
260  const ElementAccessor<3> ele = dh_cell.elm();
261 
262  //set global dof for element (pressure)
264 
265  //shortcuts
266  const unsigned int nsides = ele->n_sides();
267  LinSys *ls = ad_->lin_sys;
268 
269  unsigned int side_row, edge_row;
270 
271  dirichlet_edge.resize(nsides);
272  for (unsigned int i = 0; i < nsides; i++) {
273 
274  side_row = loc_side_dofs[i]; //local
275  edge_row = loc_edge_dofs[i]; //local
276  loc_system_.row_dofs[side_row] = loc_system_.col_dofs[side_row] = global_dofs_[side_row]; //global
277  loc_system_.row_dofs[edge_row] = loc_system_.col_dofs[edge_row] = global_dofs_[edge_row]; //global
278 
279  dirichlet_edge[i] = 0;
280  Side side = *dh_cell.elm().side(i);
281  if (side.is_boundary()) {
282  Boundary bcd = side.cond();
283  ElementAccessor<3> b_ele = bcd.element_accessor();
284  DarcyMH::EqFields::BC_Type type = (DarcyMH::EqFields::BC_Type)af_->bc_type.value(b_ele.centre(), b_ele);
285 
286  double cross_section = af_->cross_section.value(ele.centre(), ele);
287 
288  if ( type == DarcyMH::EqFields::none) {
289  // homogeneous neumann
290  } else if ( type == DarcyMH::EqFields::dirichlet ) {
291  double bc_pressure = af_->bc_pressure.value(b_ele.centre(), b_ele);
292  loc_system_.set_solution(loc_edge_dofs[i],bc_pressure,-1);
293  dirichlet_edge[i] = 1;
294 
295  } else if ( type == DarcyMH::EqFields::total_flux) {
296  // internally we work with outward flux
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);
300 
301  dirichlet_edge[i] = 2; // to be skipped in LMH source assembly
302  loc_system_.add_value(edge_row, edge_row,
303  -b_ele.measure() * bc_sigma * cross_section,
304  (bc_flux - bc_sigma * bc_pressure) * b_ele.measure() * cross_section);
305  }
306  else if (type == DarcyMH::EqFields::seepage) {
307  ad_->is_linear=false;
308 
309  unsigned int loc_edge_idx = bcd.bc_ele_idx();
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;
314 
315  // ** Update BC type. **
316  if (use_dirichlet_switch) {
317  if (switch_dirichlet) {
318  // check and possibly switch to flux BC
319  // The switch raise error on the corresponding edge row.
320  // Magnitude of the error is abs(solution_flux - side_flux).
321  ASSERT(ad_->dh_->distr()->is_local(global_dofs_[side_row]))(global_dofs_[side_row]);
322  unsigned int loc_side_row = local_dofs_[side_row];
323  double & solution_flux = ls->get_solution_array()[loc_side_row];
324 
325  if ( solution_flux < side_flux) {
326  //DebugOut().fmt("x: {}, to neum, p: {} f: {} -> f: {}\n", b_ele.centre()[0], bc_pressure, solution_flux, side_flux);
327  solution_flux = side_flux;
328  switch_dirichlet=0;
329  }
330  } else {
331  // check and possibly switch to pressure BC
332  // TODO: What is the appropriate DOF in not local?
333  // The switch raise error on the corresponding side row.
334  // Magnitude of the error is abs(solution_head - bc_pressure)
335  // Since usually K is very large, this error would be much
336  // higher then error caused by the inverse switch, this
337  // cause that a solution with the flux violating the
338  // flux inequality leading may be accepted, while the error
339  // in pressure inequality is always satisfied.
340  ASSERT(ad_->dh_->distr()->is_local(global_dofs_[edge_row]))(global_dofs_[edge_row]);
341  unsigned int loc_edge_row = local_dofs_[edge_row];
342  double & solution_head = ls->get_solution_array()[loc_edge_row];
343 
344  if ( solution_head > bc_pressure) {
345  //DebugOut().fmt("x: {}, to dirich, p: {} -> p: {} f: {}\n",b_ele.centre()[0], solution_head, bc_pressure, bc_flux);
346  solution_head = bc_pressure;
347  switch_dirichlet=1;
348  }
349  }
350  }
351 
352  // ** Apply BCUpdate BC type. **
353  // Force Dirichlet type during the first iteration of the unsteady case.
354  if (switch_dirichlet || ad_->force_no_neumann_bc ) {
355  //DebugOut().fmt("x: {}, dirich, bcp: {}\n", b_ele.centre()[0], bc_pressure);
356  loc_system_.set_solution(loc_edge_dofs[i],bc_pressure, -1);
357  dirichlet_edge[i] = 1;
358  } else {
359  //DebugOut()("x: {}, neuman, q: {} bcq: {}\n", b_ele.centre()[0], side_flux, bc_flux);
360  loc_system_.add_value(edge_row, side_flux);
361  }
362 
363  } else if (type == DarcyMH::EqFields::river) {
364  ad_->is_linear=false;
365 
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);
370  ASSERT(ad_->dh_->distr()->is_local(global_dofs_[edge_row]))(global_dofs_[edge_row]);
371  unsigned int loc_edge_row = local_dofs_[edge_row];
372  double & solution_head = ls->get_solution_array()[loc_edge_row];
373 
374  // Force Robin type during the first iteration of the unsteady case.
375  if (solution_head > bc_switch_pressure || ad_->force_no_neumann_bc) {
376  // Robin BC
377  //DebugOut().fmt("x: {}, robin, bcp: {}\n", b_ele.centre()[0], bc_pressure);
378  loc_system_.add_value(edge_row, edge_row,
379  -b_ele.measure() * bc_sigma * cross_section,
380  b_ele.measure() * cross_section * (bc_flux - bc_sigma * bc_pressure) );
381  } else {
382  // Neumann BC
383  //DebugOut().fmt("x: {}, neuman, q: {} bcq: {}\n", b_ele.centre()[0], bc_switch_pressure, bc_pressure);
384  double bc_total_flux = bc_flux + bc_sigma*(bc_switch_pressure - bc_pressure);
385 
386  loc_system_.add_value(edge_row, bc_total_flux * b_ele.measure() * cross_section);
387  }
388  }
389  else {
390  THROW( ExcBCNotSupported() );
391  }
392  }
393  loc_system_.add_value(side_row, edge_row, 1.0);
394  loc_system_.add_value(edge_row, side_row, 1.0);
395  }
396  }
397 
398  void assemble_sides(const DHCellAccessor& dh_cell)
399  {
400  const ElementAccessor<3> ele = dh_cell.elm();
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;
404 
405  assemble_sides_scale(dh_cell, scale);
406  }
407 
408  void assemble_sides_scale(const DHCellAccessor& dh_cell, double scale)
409  {
410  arma::vec3 &gravity_vec = ad_->gravity_vec_;
411  const ElementAccessor<3> ele = dh_cell.elm();
412 
413  fe_values_.reinit(ele);
414  unsigned int ndofs = fe_values_.n_dofs();
415  unsigned int qsize = fe_values_.n_points();
416  auto velocity = fe_values_.vector_view(0);
417 
418  for (unsigned int k=0; k<qsize; k++)
419  for (unsigned int i=0; i<ndofs; i++){
420  double rhs_val =
421  arma::dot(gravity_vec,velocity.value(i,k))
422  * fe_values_.JxW(k);
423  loc_system_.add_value(i, rhs_val);
424 
425  for (unsigned int j=0; j<ndofs; j++){
426  double mat_val =
427  arma::dot(velocity.value(i,k), //TODO: compute anisotropy before
428  (af_->anisotropy.value(ele.centre(), ele )).i()
429  * velocity.value(j,k))
430  * scale * fe_values_.JxW(k);
431 
432  loc_system_.add_value(i, j, mat_val);
433  }
434  }
435 
436  // assemble matrix for weights in BDDCML
437  // approximation to diagonal of
438  // S = -C - B*inv(A)*B'
439  // as
440  // diag(S) ~ - diag(C) - 1./diag(A)
441  // the weights form a partition of unity to average a discontinuous solution from neighbouring subdomains
442  // to a continuous one
443  // it is important to scale the effect - if conductivity is low for one subdomain and high for the other,
444  // trust more the one with low conductivity - it will be closer to the truth than an arithmetic average
445  if ( typeid(*ad_->lin_sys) == typeid(LinSys_BDDC) ) {
446  const arma::mat& local_matrix = loc_system_.get_matrix();
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);
450 
451  unsigned int side_row = loc_system_.row_dofs[loc_side_dofs[i]];
452  unsigned int edge_row = loc_system_.row_dofs[loc_edge_dofs[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 );
455  }
456  }
457  }
458 
459 
461  // set block B, B': element-side, side-element
462 
463  for(unsigned int side = 0; side < loc_side_dofs.size(); side++){
466  }
467 
468  if ( typeid(*ad_->lin_sys) == typeid(LinSys_BDDC) ) {
469  double val_ele = 1.;
470  static_cast<LinSys_BDDC*>(ad_->lin_sys)->
471  diagonal_weights_set_value( loc_system_.row_dofs[loc_ele_dof], val_ele );
472  }
473  }
474 
475 
477  //D, E',E block: compatible connections: element-edge
478  auto ele = dh_cell.elm(); //ElementAccessor<3>
479 
480  // no Neighbours => nothing to asssemble here
481  if(dh_cell.elm()->n_neighs_vb() == 0) return;
482 
483  std::vector<LongIdx> higher_dim_dofs;
484  //DebugOut() << "adc " << print_var(this) << print_var(side_quad_.size());
485  unsigned int i = 0;
486  for ( DHCellSide neighb_side : dh_cell.neighb_sides() ) {
487  // every compatible connection adds a 2x2 matrix involving
488  // current element pressure and a connected edge pressure
489 
492 
493  // // read neighbor dofs (dh dofhandler)
494  // // neighbor cell owning neighb_side
495  DHCellAccessor dh_neighb_cell = neighb_side.cell();
496 
497  higher_dim_dofs.resize(dh_neighb_cell.n_dofs());
498  dh_neighb_cell.get_dof_indices(higher_dim_dofs);
499 
500  // local index of pedge dof on neighboring cell
501  // (dim+2) is number of edges of higher dim element
502  // TODO: replace with DHCell getter when available for FESystem component
503  const unsigned int t = dh_neighb_cell.n_dofs() - (dim+2) + neighb_side.side().side_idx();
504  loc_system_vb_.row_dofs[1] = loc_system_vb_.col_dofs[1] = higher_dim_dofs[t];
505 
506  assembly_local_vb(ele, neighb_side);
507 
509  ad_->lin_sys->set_local_system(loc_system_vb_);
510 
511  // update matrix for weights in BDDCML
512  if ( typeid(*ad_->lin_sys) == typeid(LinSys_BDDC) ) {
513  int ind = loc_system_vb_.row_dofs[1];
514  // there is -value on diagonal in block C!
515  double new_val = loc_system_vb_.get_matrix()(0,0);
516  static_cast<LinSys_BDDC*>(ad_->lin_sys)->diagonal_weights_set_value( ind, new_val );
517  }
518  ++i;
519  }
520  }
521 
523 
524  for(DHCellSide dh_side : dh_cell.side_range()){
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]]},
529  {1}, 0);
530  }
531  }
532  }
533 
534 
535 
536  // assembly volume integrals
540 
541  NeighSideValues< (dim<3) ? dim : 2> ngh_values_;
542 
543  // Interpolation of velocity into barycenters
546 
547  // data shared by assemblers of different dimension
551 
556  unsigned int loc_ele_dof;
557 
558  std::shared_ptr<MortarAssemblyBase> mortar_assembly;
559 
562 };
563 
564 
565 #endif /* SRC_FLOW_DARCY_FLOW_ASSEMBLY_HH_ */
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:301
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
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)
FEValues< 3 > fe_values_
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()
LocDofVec local_dofs_
std::vector< unsigned int > loc_side_dofs
unsigned int loc_ele_dof
FEValues< 3 > velocity_interpolation_fv_
std::shared_ptr< MortarAssemblyBase > mortar_assembly
LocalSystem loc_system_
void add_fluxes_in_balance_matrix(const DHCellAccessor &dh_cell)
std::vector< unsigned int > loc_edge_dofs
FE_RT0< dim > fe_rt_
AssemblyFieldsPtrMH af_
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
AssemblyDataPtrMH ad_
uint bc_ele_idx()
Definition: accessors.hh:374
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.
double measure() const
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...
Definition: accessors.hh:215
arma::vec::fixed< spacedim > centre() const
Computes the barycenter.
unsigned int n_neighs_vb() const
Return number of neighbours.
Definition: elements.h:65
unsigned int n_sides() const
Definition: elements.h:131
unsigned int dim() const
Definition: elements.h:120
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
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
double * get_solution_array()
Definition: linsys.hh:325
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.
Definition: local_system.cc:82
LocDofVec col_dofs
Definition: local_system.hh:54
void eliminate_solution()
const arma::mat & get_matrix()
Definition: local_system.hh:82
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.
LocDofVec row_dofs
Definition: local_system.hh:54
void reset()
Resets the matrix, RHS, dofs to zero and clears solution settings.
Definition: local_system.cc:46
FE_P_disc< dim+1 > fe_p_disc_
FEValues< 3 > fe_side_values_
Symmetric Gauss-Legendre quadrature formulae on simplices.
Boundary cond() const
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.
Definition: exceptions.hh:53
arma::Col< IntIdx > LocDofVec
Definition: index_types.hh:28
static constexpr bool value
Definition: json.hpp:87
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
Definition: armor.hh:888
ArmaVec< double, N > vec
Definition: armor.hh:885
BasicData Data
Definition: format.h:848
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.
Definition: update_flags.hh:87
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_quadrature_points
Transformed quadrature points.