Flow123d  release_3.0.0-973-g92f55e826
darcy_flow_assembly.hh
Go to the documentation of this file.
1 /*
2  * darcy_flow_assembly.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 "mesh/long_idx.hh"
12 #include "mesh/mesh.h"
13 #include "mesh/accessors.hh"
14 #include "mesh/neighbours.h"
15 #include "fem/mapping_p1.hh"
16 #include "fem/fe_p.hh"
17 #include "fem/fe_values.hh"
18 #include "fem/fe_rt.hh"
19 #include "fem/fe_values_views.hh"
20 #include "quadrature/quadrature_lib.hh"
21 #include "flow/mh_dofhandler.hh"
22 
23 #include "la/linsys.hh"
24 #include "la/linsys_PETSC.hh"
25 #include "la/linsys_BDDC.hh"
26 #include "la/schur.hh"
27 
28 #include "la/local_system.hh"
29 
30 #include "coupling/balance.hh"
31 #include "flow/mortar_assembly.hh"
32 
33 
34 class AssemblyBase
35 {
36 public:
37  typedef std::shared_ptr<DarcyMH::EqData> AssemblyDataPtr;
38  typedef std::vector<std::shared_ptr<AssemblyBase> > MultidimAssembly;
39 
40  virtual ~AssemblyBase() {}
41 
42  /**
43  * Generic creator of multidimensional assembly, i.e. vector of
44  * particular assembly objects.
45  */
46  template< template<int dim> class Impl >
47  static MultidimAssembly create(typename Impl<1>::AssemblyDataPtr data) {
48  return { std::make_shared<Impl<1> >(data),
49  std::make_shared<Impl<2> >(data),
50  std::make_shared<Impl<3> >(data) };
51 
52  }
53 
54 // virtual LocalSystem & get_local_system() = 0;
55  virtual void fix_velocity(LocalElementAccessorBase<3> ele_ac) = 0;
56  virtual void assemble(LocalElementAccessorBase<3> ele_ac) = 0;
57 
58  // assembly compatible neighbourings
59  virtual void assembly_local_vb(ElementAccessor<3> ele, Neighbour *ngh) = 0;
60 
61  // compute velocity value in the barycenter
62  // TODO: implement and use general interpolations between discrete spaces
63  virtual arma::vec3 make_element_vector(ElementAccessor<3> ele) = 0;
64 
65  virtual void update_water_content(LocalElementAccessorBase<3> ele)
66  {}
67 
68 protected:
69 
70  virtual void assemble_sides(LocalElementAccessorBase<3> ele) =0;
71 
72  virtual void assemble_source_term(LocalElementAccessorBase<3> ele)
73  {}
74 };
75 
76 
77 
78 template <int dim>
79 class NeighSideValues {
80 private:
81  // assembly face integrals (BC)
82  MappingP1<dim+1,3> side_map_;
83  QGauss<dim> side_quad_;
84  FE_P_disc<dim+1> fe_p_disc_;
85 public:
86  NeighSideValues<dim>()
87  : side_quad_(1),
88  fe_p_disc_(0),
89  fe_side_values_(side_map_, side_quad_, fe_p_disc_, update_normal_vectors)
90  {}
91  FESideValues<dim+1,3> fe_side_values_;
92 
93 };
94 
95 
96 
97 template<int dim>
98 class AssemblyMH : public AssemblyBase
99 {
100 public:
101  AssemblyMH<dim>(AssemblyDataPtr data)
102  : quad_(3),
103  fe_values_(map_, quad_, fe_rt_,
104  update_values | update_gradients | update_JxW_values | update_quadrature_points),
105 
106  velocity_interpolation_quad_(0), // veloctiy values in barycenter
107  velocity_interpolation_fv_(map_,velocity_interpolation_quad_, fe_rt_, update_values | update_quadrature_points),
108 
109  ad_(data),
110  loc_system_(size(), size()),
111  loc_system_vb_(2,2)
112  {
113  // local numbering of dofs for MH system
114  unsigned int nsides = dim+1;
115  loc_side_dofs.resize(nsides);
116  loc_ele_dof = nsides;
117  loc_edge_dofs.resize(nsides);
118  for(unsigned int i = 0; i < nsides; i++){
119  loc_side_dofs[i] = i;
120  loc_edge_dofs[i] = nsides + i + 1;
121  }
122  //DebugOut() << print_var(this) << print_var(side_quad_.size());
123 
124  // create local sparsity pattern
125  arma::umat sp(size(),size());
126  sp.zeros();
127  sp.submat(0, 0, nsides, nsides).ones();
128  sp.diag().ones();
129  // armadillo 8.4.3 bug with negative sub diagonal index
130  // sp.diag(nsides+1).ones();
131  // sp.diag(-nsides-1).ones();
132  // sp.print();
133 
134  sp.submat(0, nsides+1, nsides-1, size()-1).diag().ones();
135  sp.submat(nsides+1, 0, size()-1, nsides-1).diag().ones();
136 
137  loc_system_.set_sparsity(sp);
138 
139  // local system 2x2 for vb neighbourings is full matrix
140  // this matrix cannot be influenced by any BC (no elimination can take place)
141  sp.ones(2,2);
142  loc_system_vb_.set_sparsity(sp);
143 
144  if (ad_->mortar_method_ == DarcyMH::MortarP0) {
145  mortar_assembly = std::make_shared<P0_CouplingAssembler>(ad_);
146  } else if (ad_->mortar_method_ == DarcyMH::MortarP1) {
147  mortar_assembly = std::make_shared<P1_CouplingAssembler>(ad_);
148  }
149 
150  }
151 
152 
153  ~AssemblyMH<dim>() override
154  {}
155 
156 // LocalSystem& get_local_system() override
157 // { return loc_system_;}
158 
159  void fix_velocity(LocalElementAccessorBase<3> ele_ac) override
160  {
161  if (mortar_assembly)
162  mortar_assembly->fix_velocity(ele_ac);
163  }
164 
165  void assemble(LocalElementAccessorBase<3> ele_ac) override
166  {
167  ASSERT_EQ_DBG(ele_ac.dim(), dim);
168  loc_system_.reset();
169 
170  set_dofs_and_bc(ele_ac);
171 
172  assemble_sides(ele_ac);
173  assemble_element(ele_ac);
174  assemble_source_term(ele_ac);
175 
176  ad_->lin_sys->set_local_system(loc_system_);
177 
178  assembly_dim_connections(ele_ac);
179 
180  if (ad_->balance != nullptr)
181  add_fluxes_in_balance_matrix(ele_ac);
182 
183  if (mortar_assembly)
184  mortar_assembly->assembly(ele_ac);
185  }
186 
187  void assembly_local_vb(ElementAccessor<3> ele, Neighbour *ngh) override
188  {
189  ASSERT_LT_DBG(ele->dim(), 3);
190  //DebugOut() << "alv " << print_var(this);
191  //START_TIMER("Assembly<dim>::assembly_local_vb");
192  // compute normal vector to side
193  arma::vec3 nv;
194  ElementAccessor<3> ele_higher = ad_->mesh->element_accessor( ngh->side()->element().idx() );
195  ngh_values_.fe_side_values_.reinit(ele_higher, ngh->side()->side_idx());
196  nv = ngh_values_.fe_side_values_.normal_vector(0);
197 
198  double value = ad_->sigma.value( ele.centre(), ele) *
199  2*ad_->conductivity.value( ele.centre(), ele) *
200  arma::dot(ad_->anisotropy.value( ele.centre(), ele)*nv, nv) *
201  ad_->cross_section.value( ngh->side()->centre(), ele_higher ) * // cross-section of higher dim. (2d)
202  ad_->cross_section.value( ngh->side()->centre(), ele_higher ) /
203  ad_->cross_section.value( ele.centre(), ele ) * // crossection of lower dim.
204  ngh->side()->measure();
205 
206  loc_system_vb_.add_value(0,0, -value);
207  loc_system_vb_.add_value(0,1, value);
208  loc_system_vb_.add_value(1,0, value);
209  loc_system_vb_.add_value(1,1, -value);
210  }
211 
212 
213  arma::vec3 make_element_vector(ElementAccessor<3> ele) override
214  {
215  //START_TIMER("Assembly<dim>::make_element_vector");
216  arma::vec3 flux_in_center;
217  flux_in_center.zeros();
218 
219  velocity_interpolation_fv_.reinit(ele);
220  for (unsigned int li = 0; li < ele->n_sides(); li++) {
221  flux_in_center += ad_->mh_dh->side_flux( *(ele.side( li ) ) )
222  * velocity_interpolation_fv_.vector_view(0).value(li,0);
223  }
224 
225  flux_in_center /= ad_->cross_section.value(ele.centre(), ele );
226  return flux_in_center;
227  }
228 
229 protected:
230  static const unsigned int size()
231  {
232  // dofs: velocity, pressure, edge pressure
233  return RefElement<dim>::n_sides + 1 + RefElement<dim>::n_sides;
234  }
235 
236  void set_dofs_and_bc(LocalElementAccessorBase<3> ele_ac){
237 
238  ASSERT_DBG(ele_ac.dim() == dim);
239 
240  //set global dof for element (pressure)
241  loc_system_.row_dofs[loc_ele_dof] = loc_system_.col_dofs[loc_ele_dof] = ele_ac.ele_row();
242 
243  //shortcuts
244  const unsigned int nsides = ele_ac.n_sides();
245  LinSys *ls = ad_->lin_sys;
246 
247  Boundary *bcd;
248  unsigned int side_row, edge_row;
249 
250  dirichlet_edge.resize(nsides);
251  for (unsigned int i = 0; i < nsides; i++) {
252 
253  side_row = loc_side_dofs[i]; //local
254  edge_row = loc_edge_dofs[i]; //local
255  loc_system_.row_dofs[side_row] = loc_system_.col_dofs[side_row] = ele_ac.side_row(i); //global
256  loc_system_.row_dofs[edge_row] = loc_system_.col_dofs[edge_row] = ele_ac.edge_row(i); //global
257 
258  bcd = ele_ac.side(i)->cond();
259  dirichlet_edge[i] = 0;
260  if (bcd) {
261  ElementAccessor<3> b_ele = bcd->element_accessor();
262  DarcyMH::EqData::BC_Type type = (DarcyMH::EqData::BC_Type)ad_->bc_type.value(b_ele.centre(), b_ele);
263 
264  double cross_section = ad_->cross_section.value(ele_ac.centre(), ele_ac.element_accessor());
265 
266  if ( type == DarcyMH::EqData::none) {
267  // homogeneous neumann
268  } else if ( type == DarcyMH::EqData::dirichlet ) {
269  double bc_pressure = ad_->bc_pressure.value(b_ele.centre(), b_ele);
270  loc_system_.set_solution(loc_edge_dofs[i],bc_pressure,-1);
271  dirichlet_edge[i] = 1;
272 
273  } else if ( type == DarcyMH::EqData::total_flux) {
274  // internally we work with outward flux
275  double bc_flux = -ad_->bc_flux.value(b_ele.centre(), b_ele);
276  double bc_pressure = ad_->bc_pressure.value(b_ele.centre(), b_ele);
277  double bc_sigma = ad_->bc_robin_sigma.value(b_ele.centre(), b_ele);
278 
279  dirichlet_edge[i] = 2; // to be skipped in LMH source assembly
280  loc_system_.add_value(edge_row, edge_row,
281  -b_ele.measure() * bc_sigma * cross_section,
282  (bc_flux - bc_sigma * bc_pressure) * b_ele.measure() * cross_section);
283  }
284  else if (type==DarcyMH::EqData::seepage) {
285  ad_->is_linear=false;
286 
287  unsigned int loc_edge_idx = bcd->bc_ele_idx_;
288  char & switch_dirichlet = ad_->bc_switch_dirichlet[loc_edge_idx];
289  double bc_pressure = ad_->bc_switch_pressure.value(b_ele.centre(), b_ele);
290  double bc_flux = -ad_->bc_flux.value(b_ele.centre(), b_ele);
291  double side_flux = bc_flux * b_ele.measure() * cross_section;
292 
293  // ** Update BC type. **
294  if (switch_dirichlet) {
295  // check and possibly switch to flux BC
296  // The switch raise error on the corresponding edge row.
297  // Magnitude of the error is abs(solution_flux - side_flux).
298  ASSERT_DBG(ad_->mh_dh->rows_ds->is_local(ele_ac.side_row(i)))(ele_ac.side_row(i));
299  unsigned int loc_side_row = ele_ac.side_local_row(i);
300  double & solution_flux = ls->get_solution_array()[loc_side_row];
301 
302  if ( solution_flux < side_flux) {
303  //DebugOut().fmt("x: {}, to neum, p: {} f: {} -> f: {}\n", b_ele.centre()[0], bc_pressure, solution_flux, side_flux);
304  solution_flux = side_flux;
305  switch_dirichlet=0;
306  }
307  } else {
308  // check and possibly switch to pressure BC
309  // TODO: What is the appropriate DOF in not local?
310  // The switch raise error on the corresponding side row.
311  // Magnitude of the error is abs(solution_head - bc_pressure)
312  // Since usually K is very large, this error would be much
313  // higher then error caused by the inverse switch, this
314  // cause that a solution with the flux violating the
315  // flux inequality leading may be accepted, while the error
316  // in pressure inequality is always satisfied.
317  ASSERT_DBG(ad_->mh_dh->rows_ds->is_local(ele_ac.edge_row(i)))(ele_ac.edge_row(i));
318  unsigned int loc_edge_row = ele_ac.edge_local_row(i);
319  double & solution_head = ls->get_solution_array()[loc_edge_row];
320 
321  if ( solution_head > bc_pressure) {
322  //DebugOut().fmt("x: {}, to dirich, p: {} -> p: {} f: {}\n",b_ele.centre()[0], solution_head, bc_pressure, bc_flux);
323  solution_head = bc_pressure;
324  switch_dirichlet=1;
325  }
326  }
327 
328  // ** Apply BCUpdate BC type. **
329  // Force Dirichlet type during the first iteration of the unsteady case.
330  if (switch_dirichlet || ad_->force_bc_switch ) {
331  //DebugOut().fmt("x: {}, dirich, bcp: {}\n", b_ele.centre()[0], bc_pressure);
332  loc_system_.set_solution(loc_edge_dofs[i],bc_pressure, -1);
333  dirichlet_edge[i] = 1;
334  } else {
335  //DebugOut()("x: {}, neuman, q: {} bcq: {}\n", b_ele.centre()[0], side_flux, bc_flux);
336  loc_system_.add_value(edge_row, side_flux);
337  }
338 
339  } else if (type==DarcyMH::EqData::river) {
340  ad_->is_linear=false;
341 
342  double bc_pressure = ad_->bc_pressure.value(b_ele.centre(), b_ele);
343  double bc_switch_pressure = ad_->bc_switch_pressure.value(b_ele.centre(), b_ele);
344  double bc_flux = -ad_->bc_flux.value(b_ele.centre(), b_ele);
345  double bc_sigma = ad_->bc_robin_sigma.value(b_ele.centre(), b_ele);
346  ASSERT_DBG(ad_->mh_dh->rows_ds->is_local(ele_ac.edge_row(i)))(ele_ac.edge_row(i));
347  unsigned int loc_edge_row = ele_ac.edge_local_row(i);
348  double & solution_head = ls->get_solution_array()[loc_edge_row];
349 
350  // Force Robin type during the first iteration of the unsteady case.
351  if (solution_head > bc_switch_pressure || ad_->force_bc_switch) {
352  // Robin BC
353  //DebugOut().fmt("x: {}, robin, bcp: {}\n", b_ele.centre()[0], bc_pressure);
354  loc_system_.add_value(edge_row, edge_row,
355  -b_ele.measure() * bc_sigma * cross_section,
356  b_ele.measure() * cross_section * (bc_flux - bc_sigma * bc_pressure) );
357  } else {
358  // Neumann BC
359  //DebugOut().fmt("x: {}, neuman, q: {} bcq: {}\n", b_ele.centre()[0], bc_switch_pressure, bc_pressure);
360  double bc_total_flux = bc_flux + bc_sigma*(bc_switch_pressure - bc_pressure);
361 
362  loc_system_.add_value(edge_row, bc_total_flux * b_ele.measure() * cross_section);
363  }
364  }
365  else {
366  xprintf(UsrErr, "BC type not supported.\n");
367  }
368  }
369  loc_system_.add_value(side_row, edge_row, 1.0);
370  loc_system_.add_value(edge_row, side_row, 1.0);
371  }
372 
373 // DBGCOUT(<< "ele " << ele_ac.ele_global_idx() << ": ");
374 // for (unsigned int i = 0; i < nsides; i++) cout << loc_system_.row_dofs[loc_side_dofs[i]] << " ";
375 // cout << loc_system_.row_dofs[loc_ele_dof] << " ";
376 // for (unsigned int i = 0; i < nsides; i++) cout << loc_system_.row_dofs[loc_edge_dofs[i]] << " ";
377 // cout << "\n";
378  }
379 
380  void assemble_sides(LocalElementAccessorBase<3> ele_ac) override
381  {
382  double cs = ad_->cross_section.value(ele_ac.centre(), ele_ac.element_accessor());
383  double conduct = ad_->conductivity.value(ele_ac.centre(), ele_ac.element_accessor());
384  double scale = 1 / cs /conduct;
385 
386  assemble_sides_scale(ele_ac, scale);
387  }
388 
389  void assemble_sides_scale(LocalElementAccessorBase<3> ele_ac, double scale)
390  {
391  arma::vec3 &gravity_vec = ad_->gravity_vec_;
392 
393  ElementAccessor<3> ele =ele_ac.element_accessor();
394  fe_values_.reinit(ele);
395  unsigned int ndofs = fe_values_.get_fe()->n_dofs();
396  unsigned int qsize = fe_values_.get_quadrature()->size();
397  auto velocity = fe_values_.vector_view(0);
398 
399  for (unsigned int k=0; k<qsize; k++)
400  for (unsigned int i=0; i<ndofs; i++){
401  double rhs_val =
402  arma::dot(gravity_vec,velocity.value(i,k))
403  * fe_values_.JxW(k);
404  loc_system_.add_value(i, rhs_val);
405 
406  for (unsigned int j=0; j<ndofs; j++){
407  double mat_val =
408  arma::dot(velocity.value(i,k), //TODO: compute anisotropy before
409  (ad_->anisotropy.value(ele_ac.centre(), ele_ac.element_accessor() )).i()
410  * velocity.value(j,k))
411  * scale * fe_values_.JxW(k);
412 
413  loc_system_.add_value(i, j, mat_val);
414  }
415  }
416 
417  // assemble matrix for weights in BDDCML
418  // approximation to diagonal of
419  // S = -C - B*inv(A)*B'
420  // as
421  // diag(S) ~ - diag(C) - 1./diag(A)
422  // the weights form a partition of unity to average a discontinuous solution from neighbouring subdomains
423  // to a continuous one
424  // it is important to scale the effect - if conductivity is low for one subdomain and high for the other,
425  // trust more the one with low conductivity - it will be closer to the truth than an arithmetic average
426  if ( typeid(*ad_->lin_sys) == typeid(LinSys_BDDC) ) {
427  const arma::mat& local_matrix = loc_system_.get_matrix();
428  for(unsigned int i=0; i < ndofs; i++) {
429  double val_side = local_matrix(i,i);
430  double val_edge = -1./local_matrix(i,i);
431 
432  unsigned int side_row = loc_system_.row_dofs[loc_side_dofs[i]];
433  unsigned int edge_row = loc_system_.row_dofs[loc_edge_dofs[i]];
434  static_cast<LinSys_BDDC*>(ad_->lin_sys)->diagonal_weights_set_value( side_row, val_side );
435  static_cast<LinSys_BDDC*>(ad_->lin_sys)->diagonal_weights_set_value( edge_row, val_edge );
436  }
437  }
438  }
439 
440 
441  void assemble_element(LocalElementAccessorBase<3> ele_ac){
442  // set block B, B': element-side, side-element
443 
444  for(unsigned int side = 0; side < loc_side_dofs.size(); side++){
445  loc_system_.add_value(loc_ele_dof, loc_side_dofs[side], -1.0);
446  loc_system_.add_value(loc_side_dofs[side], loc_ele_dof, -1.0);
447  }
448 
449  if ( typeid(*ad_->lin_sys) == typeid(LinSys_BDDC) ) {
450  double val_ele = 1.;
451  static_cast<LinSys_BDDC*>(ad_->lin_sys)->
452  diagonal_weights_set_value( loc_system_.row_dofs[loc_ele_dof], val_ele );
453  }
454  }
455 
456 
457  void assembly_dim_connections(LocalElementAccessorBase<3> ele_ac){
458  //D, E',E block: compatible connections: element-edge
459  ElementAccessor<3> ele = ele_ac.element_accessor();
460 
461  // no Neighbours => nothing to asssemble here
462  if(ele->n_neighs_vb() == 0) return;
463 
464  int ele_row = ele_ac.ele_row();
465  Neighbour *ngh;
466 
467  //DebugOut() << "adc " << print_var(this) << print_var(side_quad_.size());
468  for (unsigned int i = 0; i < ele->n_neighs_vb(); i++) {
469  // every compatible connection adds a 2x2 matrix involving
470  // current element pressure and a connected edge pressure
471  ngh = ele_ac.element_accessor()->neigh_vb[i];
472  loc_system_vb_.reset();
473  loc_system_vb_.row_dofs[0] = loc_system_vb_.col_dofs[0] = ele_row;
474  loc_system_vb_.row_dofs[1] = loc_system_vb_.col_dofs[1] = ad_->mh_dh->row_4_edge[ ngh->edge_idx() ];
475 
476  assembly_local_vb(ele, ngh);
477 
478  ad_->lin_sys->set_local_system(loc_system_vb_);
479 
480  // update matrix for weights in BDDCML
481  if ( typeid(*ad_->lin_sys) == typeid(LinSys_BDDC) ) {
482  int ind = loc_system_vb_.row_dofs[1];
483  // there is -value on diagonal in block C!
484  double new_val = loc_system_vb_.get_matrix()(0,0);
485  static_cast<LinSys_BDDC*>(ad_->lin_sys)->diagonal_weights_set_value( ind, new_val );
486  }
487  }
488  }
489 
490  void add_fluxes_in_balance_matrix(LocalElementAccessorBase<3> ele_ac){
491 
492  for (unsigned int i = 0; i < ele_ac.n_sides(); i++) {
493  Boundary* bcd = ele_ac.side(i)->cond();
494 
495  if (bcd) {
496  /*
497  DebugOut().fmt("add_flux: {} {} {}\n",
498  ad_->mh_dh->el_ds->myp(),
499  ele_ac.ele_global_idx(),
500  ele_ac.side_row(i));
501  */
502  ad_->balance->add_flux_matrix_values(ad_->water_balance_idx, ele_ac.side(i),
503  {(LongIdx)(ele_ac.side_row(i))}, {1});
504  }
505  }
506  }
507 
508 
509 
510  // assembly volume integrals
511  FE_RT0<dim> fe_rt_;
512  MappingP1<dim,3> map_;
513  QGauss<dim> quad_;
514  FEValues<dim,3> fe_values_;
515 
516  NeighSideValues<dim<3?dim:2> ngh_values_;
517 
518  // Interpolation of velocity into barycenters
519  QGauss<dim> velocity_interpolation_quad_;
520  FEValues<dim,3> velocity_interpolation_fv_;
521 
522  // data shared by assemblers of different dimension
523  AssemblyDataPtr ad_;
524  std::vector<unsigned int> dirichlet_edge;
525 
526  LocalSystem loc_system_;
527  LocalSystem loc_system_vb_;
528  std::vector<unsigned int> loc_side_dofs;
529  std::vector<unsigned int> loc_edge_dofs;
530  unsigned int loc_ele_dof;
531 
532  std::shared_ptr<MortarAssemblyBase> mortar_assembly;
533 };
534 
535 
536 #endif /* SRC_FLOW_DARCY_FLOW_ASSEMBLY_HH_ */
AssemblyMH::assemble_sides_scale
void assemble_sides_scale(LocalElementAccessorBase< 3 > ele_ac, double scale)
Definition: darcy_flow_assembly.hh:389
Element::n_neighs_vb
unsigned int n_neighs_vb() const
Return number of neighbours.
Definition: elements.h:70
LocalElementAccessorBase::edge_local_row
uint edge_local_row(uint i)
Definition: mh_dofhandler.hh:176
FE_RT0
Raviart-Thomas element of order 0.
Definition: fe_rt.hh:60
AssemblyBase::update_water_content
virtual void update_water_content(LocalElementAccessorBase< 3 > ele)
Definition: darcy_flow_assembly.hh:65
DarcyMH::EqData::BC_Type
BC_Type
Definition: darcy_flow_mh.hh:155
linsys_BDDC.hh
Solver based on Multilevel BDDC - using corresponding class of OpenFTL package.
Boundary
Definition: boundaries.h:48
AssemblyMH::map_
MappingP1< dim, 3 > map_
Definition: darcy_flow_assembly.hh:512
update_gradients
@ update_gradients
Shape function gradients.
Definition: update_flags.hh:95
AssemblyMH::make_element_vector
arma::vec3 make_element_vector(ElementAccessor< 3 > ele) override
Definition: darcy_flow_assembly.hh:213
fe_rt.hh
Definitions of Raviart-Thomas finite elements.
LocalElementAccessorBase::n_sides
uint n_sides()
Definition: mh_dofhandler.hh:128
LinSys
Definition: la_linsys_new.hh:169
fe_values_views.hh
NeighSideValues::side_quad_
QGauss< dim > side_quad_
Definition: darcy_flow_assembly.hh:83
LocalSystem
Definition: local_system.hh:27
UsrErr
@ UsrErr
Definition: system.hh:64
AssemblyMH::assemble_element
void assemble_element(LocalElementAccessorBase< 3 > ele_ac)
Definition: darcy_flow_assembly.hh:441
RefElement
Definition: ref_element.hh:162
DarcyMH::EqData::dirichlet
@ dirichlet
Definition: darcy_flow_mh.hh:157
mh_dofhandler.hh
neighbours.h
AssemblyMH::fe_rt_
FE_RT0< dim > fe_rt_
Definition: darcy_flow_assembly.hh:511
FEValues::reinit
void reinit(ElementAccessor< 3 > &cell)
Update cell-dependent data (gradients, Jacobians etc.)
Definition: fe_values.cc:528
AssemblyBase::assemble_sides
virtual void assemble_sides(LocalElementAccessorBase< 3 > ele)=0
AssemblyMH
Definition: darcy_flow_assembly.hh:98
Element::dim
unsigned int dim() const
Definition: elements.h:124
AssemblyBase::MultidimAssembly
std::vector< std::shared_ptr< AssemblyBase > > MultidimAssembly
Definition: darcy_flow_assembly.hh:38
fe_values.hh
Class FEValues calculates finite element data on the actual cells such as shape function values,...
Neighbour::side
SideIter side()
Definition: neighbours.h:147
value
static constexpr bool value
Definition: json.hpp:87
ASSERT_DBG
#define ASSERT_DBG(expr)
Definition: asserts.hh:349
linsys.hh
Wrappers for linear systems based on MPIAIJ and MATIS format.
AssemblyMH::fix_velocity
void fix_velocity(LocalElementAccessorBase< 3 > ele_ac) override
Definition: darcy_flow_assembly.hh:159
update_values
@ update_values
Shape function values.
Definition: update_flags.hh:87
AssemblyMH::assemble
void assemble(LocalElementAccessorBase< 3 > ele_ac) override
Definition: darcy_flow_assembly.hh:165
AssemblyBase::assembly_local_vb
virtual void assembly_local_vb(ElementAccessor< 3 > ele, Neighbour *ngh)=0
std::vector
Definition: doxy_dummy_defs.hh:7
Side::measure
double measure() const
Definition: sides.cc:30
ElementAccessor< 3 >
AssemblyMH::size
static const unsigned int size()
Definition: darcy_flow_assembly.hh:230
LocalElementAccessorBase< 3 >
DarcyMH::EqData::total_flux
@ total_flux
Definition: darcy_flow_mh.hh:158
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: update_flags.hh:102
arma::vec3
Definition: doxy_dummy_defs.hh:17
ASSERT_EQ_DBG
#define ASSERT_EQ_DBG(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:331
AssemblyBase::assemble
virtual void assemble(LocalElementAccessorBase< 3 > ele_ac)=0
update_normal_vectors
@ update_normal_vectors
Normal vectors.
Definition: update_flags.hh:124
linsys_PETSC.hh
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
AssemblyBase::assemble_source_term
virtual void assemble_source_term(LocalElementAccessorBase< 3 > ele)
Definition: darcy_flow_assembly.hh:72
Neighbour
Definition: neighbours.h:119
FEValuesBase::get_fe
FiniteElement< dim > * get_fe() const
Returns the finite element in use.
Definition: fe_values.hh:432
FEValuesBase::get_quadrature
Quadrature< dim > * get_quadrature() const
Returns the quadrature in use.
Definition: fe_values.hh:424
Side::cond
Boundary * cond() const
Definition: side_impl.hh:71
AssemblyMH::add_fluxes_in_balance_matrix
void add_fluxes_in_balance_matrix(LocalElementAccessorBase< 3 > ele_ac)
Definition: darcy_flow_assembly.hh:490
FEValues< dim, 3 >
fe_p.hh
Definitions of basic Lagrangean finite elements with polynomial shape functions.
LocalElementAccessorBase::centre
const arma::vec3 centre() const
Definition: mh_dofhandler.hh:136
mortar_assembly.hh
DarcyMH::EqData::none
@ none
Definition: darcy_flow_mh.hh:156
AssemblyBase::fix_velocity
virtual void fix_velocity(LocalElementAccessorBase< 3 > ele_ac)=0
DarcyMH::EqData::river
@ river
Definition: darcy_flow_mh.hh:160
AssemblyDataPtr
std::shared_ptr< DarcyMH::EqData > AssemblyDataPtr
Definition: mortar_assembly.hh:19
accessors.hh
LocalElementAccessorBase::ele_row
uint ele_row()
Definition: mh_dofhandler.hh:156
AssemblyMH::assemble_sides
void assemble_sides(LocalElementAccessorBase< 3 > ele_ac) override
Definition: darcy_flow_assembly.hh:380
quadrature_lib.hh
Definitions of particular quadrature rules on simplices.
xprintf
#define xprintf(...)
Definition: system.hh:92
schur.hh
Assembly explicit Schur complement for the given linear system. Provides method for resolution of the...
AssemblyMH::assembly_dim_connections
void assembly_dim_connections(LocalElementAccessorBase< 3 > ele_ac)
Definition: darcy_flow_assembly.hh:457
FEValuesBase::JxW
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:336
Neighbour::edge_idx
unsigned int edge_idx()
Definition: neighbours.h:153
LinSys_BDDC
Definition: linsys_BDDC.hh:41
local_system.hh
FEValuesBase::vector_view
const FEValuesViews::Vector< dim, spacedim > & vector_view(unsigned int i) const
Accessor to vector values of multicomponent FE.
Definition: fe_values.hh:388
Side::element
ElementAccessor< 3 > element() const
Definition: side_impl.hh:53
AssemblyMH::fe_values_
FEValues< dim, 3 > fe_values_
Definition: darcy_flow_assembly.hh:514
mesh.h
Boundary::element_accessor
ElementAccessor< 3 > element_accessor()
Definition: boundaries.cc:48
NeighSideValues::side_map_
MappingP1< dim+1, 3 > side_map_
Definition: darcy_flow_assembly.hh:82
NeighSideValues::fe_p_disc_
FE_P_disc< dim+1 > fe_p_disc_
Definition: darcy_flow_assembly.hh:84
AssemblyMH::quad_
QGauss< dim > quad_
Definition: darcy_flow_assembly.hh:513
Side::centre
arma::vec3 centre() const
Definition: sides.cc:121
Side::side_idx
unsigned int side_idx() const
Definition: side_impl.hh:82
AssemblyMH::set_dofs_and_bc
void set_dofs_and_bc(LocalElementAccessorBase< 3 > ele_ac)
Definition: darcy_flow_assembly.hh:236
QGauss
Symmetric Gauss-Legendre quadrature formulae on simplices.
Definition: quadrature_lib.hh:33
Element::neigh_vb
Neighbour ** neigh_vb
Definition: elements.h:87
LocalElementAccessorBase::side_row
uint side_row(uint i)
Definition: mh_dofhandler.hh:197
LocalElementAccessorBase::edge_row
uint edge_row(uint i)
Definition: mh_dofhandler.hh:172
Element::n_sides
unsigned int n_sides() const
Definition: elements.h:135
LocalElementAccessorBase::element_accessor
ElementAccessor< 3 > element_accessor()
Definition: mh_dofhandler.hh:132
MappingP1< dim+1, 3 >
LocalElementAccessorBase::dim
uint dim()
Definition: mh_dofhandler.hh:124
AssemblyBase::~AssemblyBase
virtual ~AssemblyBase()
Definition: darcy_flow_assembly.hh:40
DarcyMH::MortarP0
@ MortarP0
Definition: darcy_flow_mh.hh:142
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: update_flags.hh:114
LocalElementAccessorBase::side_local_row
uint side_local_row(uint i)
Definition: mh_dofhandler.hh:201
NeighSideValues
Definition: darcy_flow_assembly.hh:79
ASSERT_LT_DBG
#define ASSERT_LT_DBG(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:299
AssemblyBase::create
static MultidimAssembly create(typename Impl< 1 >::AssemblyDataPtr data)
Definition: darcy_flow_assembly.hh:47
ElementAccessor::idx
unsigned int idx() const
Return local idx of element in boundary / bulk part of element vector.
Definition: accessors.hh:111
AssemblyBase
Definition: darcy_flow_assembly.hh:34
long_idx.hh
mapping_p1.hh
Class MappingP1 implements the affine transformation of the unit cell onto the actual cell.
balance.hh
AssemblyMH::assembly_local_vb
void assembly_local_vb(ElementAccessor< 3 > ele, Neighbour *ngh) override
Definition: darcy_flow_assembly.hh:187
AssemblyBase::AssemblyDataPtr
std::shared_ptr< DarcyMH::EqData > AssemblyDataPtr
Definition: darcy_flow_assembly.hh:37
FESideValues< dim+1, 3 >
DarcyMH::EqData::seepage
@ seepage
Definition: darcy_flow_mh.hh:159
ElementAccessor::side
SideIter side(const unsigned int loc_index)
Definition: accessors.hh:137
DarcyMH::MortarP1
@ MortarP1
Definition: darcy_flow_mh.hh:143
LinSys::get_solution_array
double * get_solution_array()
Definition: linsys.hh:325
AssemblyBase::make_element_vector
virtual arma::vec3 make_element_vector(ElementAccessor< 3 > ele)=0
FE_P_disc< dim+1 >
ElementAccessor::measure
double measure() const
Computes the measure of the element.
Definition: accessors.hh:254
LocalElementAccessorBase::side
SideIter side(uint i)
Definition: mh_dofhandler.hh:185
ElementAccessor::centre
arma::vec::fixed< spacedim > centre() const
Computes the barycenter.
Definition: accessors.hh:285
Boundary::bc_ele_idx_
unsigned int bc_ele_idx_
Definition: boundaries.h:79
NeighSideValues::fe_side_values_
FESideValues< dim+1, 3 > fe_side_values_
Definition: darcy_flow_assembly.hh:91
Generated by   doxygen 1.8.17