Flow123d  master-f44eb46
assembly_dg.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_dg.hh
15  * @brief
16  */
17 
18 #ifndef ASSEMBLY_DG_HH_
19 #define ASSEMBLY_DG_HH_
20 
24 #include "fem/fe_p.hh"
25 #include "fem/fe_values.hh"
27 #include "coupling/balance.hh"
29 
30 
31 
32 /**
33  * Auxiliary container class for Finite element and related objects of given dimension.
34  */
35 template <unsigned int dim, class Model>
36 class MassAssemblyDG : public AssemblyBase<dim>
37 {
38 public:
41 
42  static constexpr const char * name() { return "MassAssemblyDG"; }
43 
44  /// Constructor.
45  MassAssemblyDG(EqFields *eq_fields, EqData *eq_data)
46  : AssemblyBase<dim>(eq_data->dg_order), eq_fields_(eq_fields), eq_data_(eq_data) {
48  this->used_fields_ += eq_fields_->mass_matrix_coef;
49  this->used_fields_ += eq_fields_->retardation_coef;
50  }
51 
52  /// Destructor.
54 
55  /// Initialize auxiliary vectors and other data members
56  void initialize(ElementCacheMap *element_cache_map) {
57  this->element_cache_map_ = element_cache_map;
58 
59  fe_ = std::make_shared< FE_P_disc<dim> >(eq_data_->dg_order);
61  fe_values_.initialize(*this->quad_, *fe_, u);
62  if (dim==1) // print to log only one time
63  DebugOut() << "List of MassAssembly FEValues updates flags: " << this->print_update_flags(u);
64  ndofs_ = fe_->n_dofs();
65  dof_indices_.resize(ndofs_);
66  local_matrix_.resize(4*ndofs_*ndofs_);
69 
70  }
71 
72 
73  /// Assemble integral over element
74  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
75  {
76  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
77 
78  ElementAccessor<3> elm = cell.elm();
79  unsigned int k;
80 
81  fe_values_.reinit(elm);
83 
84  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
85  {
86  // assemble the local mass matrix
87  for (unsigned int i=0; i<ndofs_; i++)
88  {
89  for (unsigned int j=0; j<ndofs_; j++)
90  {
91  local_matrix_[i*ndofs_+j] = 0;
92  k=0;
93  for (auto p : this->bulk_points(element_patch_idx) )
94  {
95  local_matrix_[i*ndofs_+j] += (eq_fields_->mass_matrix_coef(p)+eq_fields_->retardation_coef[sbi](p)) *
97  k++;
98  }
99  }
100  }
101 
102  for (unsigned int i=0; i<ndofs_; i++)
103  {
106  k=0;
107  for (auto p : this->bulk_points(element_patch_idx) )
108  {
109  local_mass_balance_vector_[i] += eq_fields_->mass_matrix_coef(p)*fe_values_.shape_value(i,k)*fe_values_.JxW(k);
110  local_retardation_balance_vector_[i] -= eq_fields_->retardation_coef[sbi](p)*fe_values_.shape_value(i,k)*fe_values_.JxW(k);
111  k++;
112  }
113  }
114 
115  eq_data_->balance_->add_mass_values(eq_data_->subst_idx()[sbi], cell, cell.get_loc_dof_indices(),
117 
119  VecSetValues(eq_data_->ret_vec[sbi], ndofs_, &(dof_indices_[0]), &(local_retardation_balance_vector_[0]), ADD_VALUES);
120  }
121  }
122 
123  /// Implements @p AssemblyBase::begin.
124  void begin() override
125  {
126  eq_data_->balance_->start_mass_assembly( eq_data_->subst_idx() );
127  }
128 
129  /// Implements @p AssemblyBase::end.
130  void end() override
131  {
132  eq_data_->balance_->finish_mass_assembly( eq_data_->subst_idx() );
133  }
134 
135  private:
136  shared_ptr<FiniteElement<dim>> fe_; ///< Finite element for the solution of the advection-diffusion equation.
137 
138  /// Data objects shared with TransportDG
141 
142  /// Sub field set contains fields used in calculation.
144 
145  unsigned int ndofs_; ///< Number of dofs
146  FEValues<3> fe_values_; ///< FEValues of object (of P disc finite element type)
147 
148  vector<LongIdx> dof_indices_; ///< Vector of global DOF indices
149  vector<PetscScalar> local_matrix_; ///< Auxiliary vector for assemble methods
150  vector<PetscScalar> local_retardation_balance_vector_; ///< Auxiliary vector for assemble mass matrix.
152 
153  template < template<IntDim...> class DimAssembly>
154  friend class GenericAssembly;
155 
156 };
157 
158 
159 // return the ratio of longest and shortest edge
161 {
162  double h_max = 0, h_min = numeric_limits<double>::infinity();
163  for (unsigned int i=0; i<e->n_nodes(); i++)
164  for (unsigned int j=i+1; j<e->n_nodes(); j++)
165  {
166  double dist = arma::norm(*e.node(i) - *e.node(j));
167  h_max = max(h_max, dist);
168  h_min = min(h_min, dist);
169  }
170  return h_max/h_min;
171 }
172 
173 /**
174  * @brief Computes average normal diffusivity over a set of points
175  *
176  * @param diff_coef Diffusion tensor.
177  * @param pts Points.
178  * @param nv Normal vector.
179  * @return double
180  */
182 {
183  double delta = 0;
184  unsigned int n = 0;
185  for (auto p : pts )
186  {
187  delta += dot(diff_coef(p)*nv, nv);
188  n++;
189  }
190  return n == 0 ? 0 : (delta/n);
191 }
192 
193 
194 /**
195  * @brief Computes the penalty parameter of the DG method on a given boundary edge.
196  *
197  * Assumption is that the edge consists of only 1 side.
198  * @param side The boundary side.
199  * @param diff_delta Average normal dispersivity K*n.n computed by diffusion_delta()
200  * @param ad_vector Advection vector.
201  * @param alpha Penalty parameter that influences the continuity
202  * of the solution (large value=more continuity).
203  */
205  const double &diff_delta,
206  const double flux,
207  const double alpha)
208 {
209  return 0.5*fabs(flux) + alpha/side.diameter()*diff_delta*elem_anisotropy(side.element());
210 }
211 
212 
213 /**
214  * @brief Computes advective flux.
215  *
216  * @param advection_coef Advection vector.
217  * @param pts Quadrature points.
218  * @param fv FE values object.
219  * @return double
220  */
221 template <class PointType>
223 {
224  double side_flux = 0;
225  unsigned int k=0;
226  for (auto p : pts) {
227  side_flux += arma::dot(advection_coef(p), fv.normal_vector(k))*fv.JxW(k);
228  k++;
229  }
230  return side_flux;
231 }
232 
233 
234 /**
235  * Auxiliary container class for Finite element and related objects of given dimension.
236  */
237 template <unsigned int dim, class Model>
239 {
240 public:
243 
244  static constexpr const char * name() { return "StiffnessAssemblyDG"; }
245 
246  /// Constructor.
247  StiffnessAssemblyDG(EqFields *eq_fields, EqData *eq_data)
248  : AssemblyBase<dim>(eq_data->dg_order), eq_fields_(eq_fields), eq_data_(eq_data) {
250  this->used_fields_ += eq_fields_->advection_coef;
251  this->used_fields_ += eq_fields_->diffusion_coef;
252  this->used_fields_ += eq_fields_->cross_section;
254  this->used_fields_ += eq_fields_->sources_sigma_out;
256  this->used_fields_ += eq_fields_->bc_type;
257  this->used_fields_ += eq_fields_->bc_robin_sigma;
258  }
259 
260  /// Destructor.
262  for (auto a : averages) if (a != nullptr) delete[] a;
263  for (auto a : waverages) if (a != nullptr) delete[] a;
264  for (auto a : jumps) if (a != nullptr) delete[] a;
265  }
266 
267  /// Initialize auxiliary vectors and other data members
268  void initialize(ElementCacheMap *element_cache_map) {
269  this->element_cache_map_ = element_cache_map;
270 
271  fe_ = std::make_shared< FE_P_disc<dim> >(eq_data_->dg_order);
272  fe_low_ = std::make_shared< FE_P_disc<dim-1> >(eq_data_->dg_order);
275  fe_values_.initialize(*this->quad_, *fe_, u);
276  if (dim>1) {
277  fe_values_vb_.initialize(*this->quad_low_, *fe_low_, u);
278  }
279  fe_values_side_.initialize(*this->quad_low_, *fe_, u_side);
280  if (dim==1) { // print to log only one time
281  DebugOut() << "List of StiffnessAssemblyDG FEValues (cell) updates flags: " << this->print_update_flags(u);
282  DebugOut() << "List of StiffnessAssemblyDG FEValues (side) updates flags: " << this->print_update_flags(u_side);
283  }
284  ndofs_ = fe_->n_dofs();
285  qsize_lower_dim_ = this->quad_low_->size();
286  dof_indices_.resize(ndofs_);
287  side_dof_indices_vb_.resize(2*ndofs_);
288  local_matrix_.resize(4*ndofs_*ndofs_);
289 
291  for (unsigned int sid=0; sid<eq_data_->max_edg_sides; sid++)
292  {
294  fe_values_vec_[sid].initialize(*this->quad_low_, *fe_,
296  }
297 
298  // index 0 = element with lower dimension,
299  // index 1 = side of element with higher dimension
300  fv_sb_.resize(2);
301  fv_sb_[0] = &fe_values_vb_;
302  fv_sb_[1] = &fe_values_side_;
303 
305  for (unsigned int s=0; s<eq_data_->max_edg_sides; s++)
306  averages[s] = new double[qsize_lower_dim_*fe_->n_dofs()];
307  waverages.resize(2);
308  jumps.resize(2);
309  for (unsigned int s=0; s<2; s++)
310  {
311  waverages[s] = new double[qsize_lower_dim_*fe_->n_dofs()];
312  jumps[s] = new double[qsize_lower_dim_*fe_->n_dofs()];
313  }
314  }
315 
316 
317  /// Assembles the cell (volume) integral into the stiffness matrix.
318  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
319  {
320  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
321  if (!cell.is_own()) return;
322 
323  ElementAccessor<3> elm = cell.elm();
324 
325  fe_values_.reinit(elm);
327  unsigned int k;
328 
329  // assemble the local stiffness matrix
330  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
331  {
332  for (unsigned int i=0; i<ndofs_; i++)
333  for (unsigned int j=0; j<ndofs_; j++)
334  local_matrix_[i*ndofs_+j] = 0;
335 
336  k=0;
337  for (auto p : this->bulk_points(element_patch_idx) )
338  {
339  for (unsigned int i=0; i<ndofs_; i++)
340  {
341  arma::vec3 Kt_grad_i = eq_fields_->diffusion_coef[sbi](p).t()*fe_values_.shape_grad(i,k);
342  double ad_dot_grad_i = arma::dot(eq_fields_->advection_coef[sbi](p), fe_values_.shape_grad(i,k));
343 
344  for (unsigned int j=0; j<ndofs_; j++)
345  local_matrix_[i*ndofs_+j] += (arma::dot(Kt_grad_i, fe_values_.shape_grad(j,k))
346  -fe_values_.shape_value(j,k)*ad_dot_grad_i
347  +eq_fields_->sources_sigma_out[sbi](p)*fe_values_.shape_value(j,k)*fe_values_.shape_value(i,k))*fe_values_.JxW(k);
348  }
349  k++;
350  }
352  }
353  }
354 
355 
356  /// Assembles the fluxes on the boundary.
357  inline void boundary_side_integral(DHCellSide cell_side)
358  {
359  ASSERT_EQ(cell_side.dim(), dim).error("Dimension of element mismatch!");
360  if (!cell_side.cell().is_own()) return;
361 
362  Side side = cell_side.side();
363  const DHCellAccessor &cell = cell_side.cell();
364 
366  fe_values_side_.reinit(side);
367  unsigned int k;
368  double gamma_l;
369 
370  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
371  {
372  std::fill(local_matrix_.begin(), local_matrix_.end(), 0);
373 
374  double side_flux = advective_flux(eq_fields_->advection_coef[sbi], this->boundary_points(cell_side), fe_values_side_);
375  double transport_flux = side_flux/side.measure();
376 
377  // On Neumann boundaries we have only term from integrating by parts the advective term,
378  // on Dirichlet boundaries we additionally apply the penalty which enforces the prescribed value.
379  auto p_side = *( this->boundary_points(cell_side).begin() );
380  auto p_bdr = p_side.point_bdr( side.cond().element_accessor() );
381  unsigned int bc_type = eq_fields_->bc_type[sbi](p_bdr);
383  {
384  // set up the parameters for DG method
385  gamma_l = DG_penalty_boundary(side,
386  diffusion_delta(eq_fields_->diffusion_coef[sbi], this->boundary_points(cell_side), fe_values_side_.normal_vector(0)),
387  transport_flux,
388  eq_fields_->dg_penalty[sbi](p_side));
389  transport_flux += gamma_l;
390  }
391 
392  // fluxes and penalty
393  k=0;
394  for (auto p : this->boundary_points(cell_side) )
395  {
396  double flux_times_JxW;
398  {
399  //sigma_ corresponds to robin_sigma
400  auto p_bdr = p.point_bdr(side.cond().element_accessor());
401  flux_times_JxW = eq_fields_->cross_section(p)*eq_fields_->bc_robin_sigma[sbi](p_bdr)*fe_values_side_.JxW(k);
402  }
404  {
405  auto p_bdr = p.point_bdr(side.cond().element_accessor());
406  flux_times_JxW = (transport_flux + eq_fields_->cross_section(p)*eq_fields_->bc_robin_sigma[sbi](p_bdr))*fe_values_side_.JxW(k);
407  }
408  else if (bc_type == AdvectionDiffusionModel::abc_inflow && side_flux < 0)
409  flux_times_JxW = 0;
410  else
411  flux_times_JxW = transport_flux*fe_values_side_.JxW(k);
412 
413  for (unsigned int i=0; i<ndofs_; i++)
414  {
415  for (unsigned int j=0; j<ndofs_; j++)
416  {
417  // flux due to advection and penalty
419 
420  // flux due to diffusion (only on dirichlet and inflow boundary)
422  local_matrix_[i*ndofs_+j] -= (arma::dot(eq_fields_->diffusion_coef[sbi](p)*fe_values_side_.shape_grad(j,k),fe_values_side_.normal_vector(k))*fe_values_side_.shape_value(i,k)
424  )*fe_values_side_.JxW(k);
425  }
426  }
427  k++;
428  }
429 
431  }
432  }
433 
434 
435  /// Assembles the fluxes between sides of elements of the same dimension.
437  ASSERT_EQ(edge_side_range.begin()->element().dim(), dim).error("Dimension of element mismatch!");
438 
439  unsigned int k;
440  double gamma_l, omega[2], transport_flux, delta[2], delta_sum;
441  double aniso1, aniso2;
442  double local_alpha=0.0;
443  int sid=0, s1, s2;
444  for( DHCellSide edge_side : edge_side_range )
445  {
446  auto dh_edge_cell = eq_data_->dh_->cell_accessor_from_element( edge_side.elem_idx() );
447  dh_edge_cell.get_dof_indices(side_dof_indices_[sid]);
448  fe_values_vec_[sid].reinit(edge_side.side());
449  ++sid;
450  }
451  arma::vec3 normal_vector = fe_values_vec_[0].normal_vector(0);
452 
453  // fluxes and penalty
454  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
455  {
456  vector<double> fluxes(edge_side_range.begin()->n_edge_sides());
457  double pflux = 0, nflux = 0; // calculate the total in- and out-flux through the edge
458  sid=0;
459  for( DHCellSide edge_side : edge_side_range )
460  {
461  fluxes[sid] = advective_flux(eq_fields_->advection_coef[sbi], this->edge_points(edge_side), fe_values_vec_[sid]) / edge_side.measure();
462  if (fluxes[sid] > 0)
463  pflux += fluxes[sid];
464  else
465  nflux += fluxes[sid];
466  ++sid;
467  }
468 
469  // precompute averages of shape functions over pairs of sides
470  s1=0;
471  for (DHCellSide edge_side : edge_side_range)
472  {
473  (void)edge_side;
474  for (unsigned int k=0; k<qsize_lower_dim_; k++)
475  {
476  for (unsigned int i=0; i<fe_->n_dofs(); i++)
477  averages[s1][k*fe_->n_dofs()+i] = fe_values_vec_[s1].shape_value(i,k)*0.5;
478  }
479  s1++;
480  }
481 
482 
483 
484  s1=0;
485  for( DHCellSide edge_side1 : edge_side_range )
486  {
487  s2=-1; // need increment at begin of loop (see conditionally 'continue' directions)
488  for( DHCellSide edge_side2 : edge_side_range )
489  {
490  s2++;
491  if (s2<=s1) continue;
492  ASSERT(edge_side1.is_valid()).error("Invalid side of edge.");
493 
494  arma::vec3 nv = fe_values_vec_[s1].normal_vector(0);
495 
496  // set up the parameters for DG method
497  // calculate the flux from edge_side1 to edge_side2
498  if (fluxes[s2] > 0 && fluxes[s1] < 0)
499  transport_flux = fluxes[s1]*fabs(fluxes[s2]/pflux);
500  else if (fluxes[s2] < 0 && fluxes[s1] > 0)
501  transport_flux = fluxes[s1]*fabs(fluxes[s2]/nflux);
502  else
503  transport_flux = 0;
504 
505  gamma_l = 0.5*fabs(transport_flux);
506 
507  delta[0] = 0;
508  delta[1] = 0;
509  for (auto p1 : this->edge_points(edge_side1) )
510  {
511  auto p2 = p1.point_on(edge_side2);
512  delta[0] += dot(eq_fields_->diffusion_coef[sbi](p1)*normal_vector,normal_vector);
513  delta[1] += dot(eq_fields_->diffusion_coef[sbi](p2)*normal_vector,normal_vector);
514  local_alpha = max(eq_fields_->dg_penalty[sbi](p1), eq_fields_->dg_penalty[sbi](p2));
515  }
516  delta[0] /= qsize_lower_dim_;
517  delta[1] /= qsize_lower_dim_;
518 
519  delta_sum = delta[0] + delta[1];
520 
521 // if (delta_sum > numeric_limits<double>::epsilon())
522  if (fabs(delta_sum) > 0)
523  {
524  omega[0] = delta[1]/delta_sum;
525  omega[1] = delta[0]/delta_sum;
526  double h = edge_side1.diameter();
527  aniso1 = elem_anisotropy(edge_side1.element());
528  aniso2 = elem_anisotropy(edge_side2.element());
529  gamma_l += local_alpha/h*aniso1*aniso2*(delta[0]*delta[1]/delta_sum);
530  }
531  else
532  for (int i=0; i<2; i++) omega[i] = 0;
533  // end of set up the parameters for DG method
534 
535  int sd[2]; bool is_side_own[2];
536  sd[0] = s1; is_side_own[0] = edge_side1.cell().is_own();
537  sd[1] = s2; is_side_own[1] = edge_side2.cell().is_own();
538 
539  // precompute jumps and weighted averages of shape functions over the pair of sides (s1,s2)
540  k=0;
541  for (auto p1 : this->edge_points(edge_side1) )
542  {
543  auto p2 = p1.point_on(edge_side2);
544  for (unsigned int i=0; i<fe_->n_dofs(); i++)
545  {
546  for (int n=0; n<2; n++)
547  {
548  jumps[n][k*fe_->n_dofs()+i] = (n==0)*fe_values_vec_[s1].shape_value(i,k) - (n==1)*fe_values_vec_[s2].shape_value(i,k);
549  waverages[n][k*fe_->n_dofs()+i] = arma::dot(eq_fields_->diffusion_coef[sbi]( (n==0 ? p1 : p2) )*fe_values_vec_[sd[n]].shape_grad(i,k),nv)*omega[n];
550  }
551  }
552  k++;
553  }
554 
555  // For selected pair of elements:
556  for (int n=0; n<2; n++)
557  {
558  if (!is_side_own[n]) continue;
559 
560  for (int m=0; m<2; m++)
561  {
562  for (unsigned int i=0; i<fe_values_vec_[sd[n]].n_dofs(); i++)
563  for (unsigned int j=0; j<fe_values_vec_[sd[m]].n_dofs(); j++)
564  local_matrix_[i*fe_values_vec_[sd[m]].n_dofs()+j] = 0;
565 
566  for (k=0; k<this->quad_low_->size(); ++k)
567  {
568  for (unsigned int i=0; i<fe_values_vec_[sd[n]].n_dofs(); i++)
569  {
570  for (unsigned int j=0; j<fe_values_vec_[sd[m]].n_dofs(); j++)
571  {
572  int index = i*fe_values_vec_[sd[m]].n_dofs()+j;
573 
574  local_matrix_[index] += (
575  // flux due to transport (applied on interior edges) (average times jump)
576  transport_flux*jumps[n][k*fe_->n_dofs()+i]*averages[sd[m]][k*fe_->n_dofs()+j]
577 
578  // penalty enforcing continuity across edges (applied on interior and Dirichlet edges) (jump times jump)
579  + gamma_l*jumps[n][k*fe_->n_dofs()+i]*jumps[m][k*fe_->n_dofs()+j]
580 
581  // terms due to diffusion
582  - jumps[n][k*fe_->n_dofs()+i]*waverages[m][k*fe_->n_dofs()+j]
583  - eq_data_->dg_variant*waverages[n][k*fe_->n_dofs()+i]*jumps[m][k*fe_->n_dofs()+j]
585  }
586  }
587  }
588  eq_data_->ls[sbi]->mat_set_values(fe_values_vec_[sd[n]].n_dofs(), &(side_dof_indices_[sd[n]][0]), fe_values_vec_[sd[m]].n_dofs(), &(side_dof_indices_[sd[m]][0]), &(local_matrix_[0]));
589  }
590  }
591  }
592  s1++;
593  }
594  }
595  }
596 
597 
598  /// Assembles the fluxes between elements of different dimensions.
599  inline void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side) {
600  if (dim == 1) return;
601  ASSERT_EQ(cell_lower_dim.dim(), dim-1).error("Dimension of element mismatch!");
602 
603  // Note: use data members csection_ and velocity_ for appropriate quantities of lower dim element
604 
605  double comm_flux[2][2];
606  unsigned int n_dofs[2];
607  ElementAccessor<3> elm_lower_dim = cell_lower_dim.elm();
608  unsigned int n_indices = cell_lower_dim.get_dof_indices(dof_indices_);
609  for(unsigned int i=0; i<n_indices; ++i) {
611  }
612  fe_values_vb_.reinit(elm_lower_dim);
613  n_dofs[0] = fv_sb_[0]->n_dofs();
614 
615  DHCellAccessor cell_higher_dim = eq_data_->dh_->cell_accessor_from_element( neighb_side.element().idx() );
616  n_indices = cell_higher_dim.get_dof_indices(dof_indices_);
617  for(unsigned int i=0; i<n_indices; ++i) {
618  side_dof_indices_vb_[i+n_dofs[0]] = dof_indices_[i];
619  }
620  fe_values_side_.reinit(neighb_side.side());
621  n_dofs[1] = fv_sb_[1]->n_dofs();
622 
623  // Testing element if they belong to local partition.
624  bool own_element_id[2];
625  own_element_id[0] = cell_lower_dim.is_own();
626  own_element_id[1] = cell_higher_dim.is_own();
627 
628  unsigned int k;
629  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++) // Optimize: SWAP LOOPS
630  {
631  for (unsigned int i=0; i<n_dofs[0]+n_dofs[1]; i++)
632  for (unsigned int j=0; j<n_dofs[0]+n_dofs[1]; j++)
633  local_matrix_[i*(n_dofs[0]+n_dofs[1])+j] = 0;
634 
635  // set transmission conditions
636  k=0;
637  for (auto p_high : this->coupling_points(neighb_side) )
638  {
639  auto p_low = p_high.lower_dim(cell_lower_dim);
640  // The communication flux has two parts:
641  // - "diffusive" term containing sigma
642  // - "advective" term representing usual upwind
643  //
644  // The calculation differs from the reference manual, since ad_coef and dif_coef have different meaning
645  // than b and A in the manual.
646  // In calculation of sigma there appears one more csection_lower in the denominator.
647  double sigma = eq_fields_->fracture_sigma[sbi](p_low)*arma::dot(eq_fields_->diffusion_coef[sbi](p_low)*fe_values_side_.normal_vector(k),fe_values_side_.normal_vector(k))*
648  2*eq_fields_->cross_section(p_high)*eq_fields_->cross_section(p_high)/(eq_fields_->cross_section(p_low)*eq_fields_->cross_section(p_low));
649 
650  double transport_flux = arma::dot(eq_fields_->advection_coef[sbi](p_high), fe_values_side_.normal_vector(k));
651 
652  comm_flux[0][0] = (sigma-min(0.,transport_flux))*fv_sb_[0]->JxW(k);
653  comm_flux[0][1] = -(sigma-min(0.,transport_flux))*fv_sb_[0]->JxW(k);
654  comm_flux[1][0] = -(sigma+max(0.,transport_flux))*fv_sb_[0]->JxW(k);
655  comm_flux[1][1] = (sigma+max(0.,transport_flux))*fv_sb_[0]->JxW(k);
656 
657  for (int n=0; n<2; n++)
658  {
659  if (!own_element_id[n]) continue;
660 
661  for (unsigned int i=0; i<n_dofs[n]; i++)
662  for (int m=0; m<2; m++)
663  for (unsigned int j=0; j<n_dofs[m]; j++)
664  local_matrix_[(i+n*n_dofs[0])*(n_dofs[0]+n_dofs[1]) + m*n_dofs[0] + j] +=
665  comm_flux[m][n]*fv_sb_[m]->shape_value(j,k)*fv_sb_[n]->shape_value(i,k) + LocalSystem::almost_zero;
666  }
667  k++;
668  }
669  eq_data_->ls[sbi]->mat_set_values(n_dofs[0]+n_dofs[1], &(side_dof_indices_vb_[0]), n_dofs[0]+n_dofs[1], &(side_dof_indices_vb_[0]), &(local_matrix_[0]));
670  }
671  }
672 
673 
674 private:
675  shared_ptr<FiniteElement<dim>> fe_; ///< Finite element for the solution of the advection-diffusion equation.
676  shared_ptr<FiniteElement<dim-1>> fe_low_; ///< Finite element for the solution of the advection-diffusion equation (dim-1).
677 
678  /// Data objects shared with TransportDG
681 
682  /// Sub field set contains fields used in calculation.
684 
685  unsigned int ndofs_; ///< Number of dofs
686  unsigned int qsize_lower_dim_; ///< Size of quadrature of dim-1
687  FEValues<3> fe_values_; ///< FEValues of object (of P disc finite element type)
688  FEValues<3> fe_values_vb_; ///< FEValues of dim-1 object (of P disc finite element type)
689  FEValues<3> fe_values_side_; ///< FEValues of object (of P disc finite element type)
690  vector<FEValues<3>> fe_values_vec_; ///< Vector of FEValues of object (of P disc finite element types)
691  vector<FEValues<3>*> fv_sb_; ///< Auxiliary vector, holds FEValues objects for assemble element-side
692 
693  vector<LongIdx> dof_indices_; ///< Vector of global DOF indices
694  vector< vector<LongIdx> > side_dof_indices_; ///< Vector of vectors of side DOF indices
695  vector<LongIdx> side_dof_indices_vb_; ///< Vector of side DOF indices (assemble element-side fluxex)
696  vector<PetscScalar> local_matrix_; ///< Auxiliary vector for assemble methods
697 
698  vector<double*> averages; ///< Auxiliary storage for averages of shape functions.
699  vector<double*> waverages; ///< Auxiliary storage for weighted averages of shape functions.
700  vector<double*> jumps; ///< Auxiliary storage for jumps of shape functions.
701 
702  template < template<IntDim...> class DimAssembly>
703  friend class GenericAssembly;
704 
705 };
706 
707 
708 /**
709  * Auxiliary container class for Finite element and related objects of given dimension.
710  */
711 template <unsigned int dim, class Model>
712 class SourcesAssemblyDG : public AssemblyBase<dim>
713 {
714 public:
717 
718  static constexpr const char * name() { return "SourcesAssemblyDG"; }
719 
720  /// Constructor.
721  SourcesAssemblyDG(EqFields *eq_fields, EqData *eq_data)
722  : AssemblyBase<dim>(eq_data->dg_order), eq_fields_(eq_fields), eq_data_(eq_data) {
724  this->used_fields_ += eq_fields_->sources_density_out;
725  this->used_fields_ += eq_fields_->sources_conc_out;
726  this->used_fields_ += eq_fields_->sources_sigma_out;
727  }
728 
729  /// Destructor.
731 
732  /// Initialize auxiliary vectors and other data members
733  void initialize(ElementCacheMap *element_cache_map) {
734  this->element_cache_map_ = element_cache_map;
735 
736  fe_ = std::make_shared< FE_P_disc<dim> >(eq_data_->dg_order);
738  fe_values_.initialize(*this->quad_, *fe_, u);
739  if (dim==1) // print to log only one time
740  DebugOut() << "List of SourcesAssemblyDG FEValues updates flags: " << this->print_update_flags(u);
741  ndofs_ = fe_->n_dofs();
742  dof_indices_.resize(ndofs_);
743  local_rhs_.resize(ndofs_);
746  }
747 
748 
749  /// Assemble integral over element
750  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
751  {
752  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
753 
754  ElementAccessor<3> elm = cell.elm();
755  unsigned int k;
756  double source;
757 
758  fe_values_.reinit(elm);
760 
761  // assemble the local stiffness matrix
762  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
763  {
764  fill_n( &(local_rhs_[0]), ndofs_, 0 );
767 
768  k=0;
769  for (auto p : this->bulk_points(element_patch_idx) )
770  {
771  source = (eq_fields_->sources_density_out[sbi](p) + eq_fields_->sources_conc_out[sbi](p)*eq_fields_->sources_sigma_out[sbi](p))*fe_values_.JxW(k);
772 
773  for (unsigned int i=0; i<ndofs_; i++)
774  local_rhs_[i] += source*fe_values_.shape_value(i,k);
775  k++;
776  }
777  eq_data_->ls[sbi]->rhs_set_values(ndofs_, &(dof_indices_[0]), &(local_rhs_[0]));
778 
779  for (unsigned int i=0; i<ndofs_; i++)
780  {
781  k=0;
782  for (auto p : this->bulk_points(element_patch_idx) )
783  {
784  local_source_balance_vector_[i] -= eq_fields_->sources_sigma_out[sbi](p)*fe_values_.shape_value(i,k)*fe_values_.JxW(k);
785  k++;
786  }
787 
789  }
790  eq_data_->balance_->add_source_values(eq_data_->subst_idx()[sbi], elm.region().bulk_idx(),
791  cell.get_loc_dof_indices(),
793  }
794  }
795 
796  /// Implements @p AssemblyBase::begin.
797  void begin() override
798  {
799  eq_data_->balance_->start_source_assembly( eq_data_->subst_idx() );
800  }
801 
802  /// Implements @p AssemblyBase::end.
803  void end() override
804  {
805  eq_data_->balance_->finish_source_assembly( eq_data_->subst_idx() );
806  }
807 
808 
809  private:
810  shared_ptr<FiniteElement<dim>> fe_; ///< Finite element for the solution of the advection-diffusion equation.
811 
812  /// Data objects shared with TransportDG
815 
816  /// Sub field set contains fields used in calculation.
818 
819  unsigned int ndofs_; ///< Number of dofs
820  FEValues<3> fe_values_; ///< FEValues of object (of P disc finite element type)
821 
822  vector<LongIdx> dof_indices_; ///< Vector of global DOF indices
823  vector<PetscScalar> local_rhs_; ///< Auxiliary vector for set_sources method.
824  vector<PetscScalar> local_source_balance_vector_; ///< Auxiliary vector for set_sources method.
825  vector<PetscScalar> local_source_balance_rhs_; ///< Auxiliary vector for set_sources method.
826 
827  template < template<IntDim...> class DimAssembly>
828  friend class GenericAssembly;
829 
830 };
831 
832 
833 /**
834  * Assembles the r.h.s. components corresponding to the Dirichlet boundary conditions..
835  */
836 template <unsigned int dim, class Model>
838 {
839 public:
842 
843  static constexpr const char * name() { return "BdrConditionAssemblyDG"; }
844 
845  /// Constructor.
846  BdrConditionAssemblyDG(EqFields *eq_fields, EqData *eq_data)
847  : AssemblyBase<dim>(eq_data->dg_order), eq_fields_(eq_fields), eq_data_(eq_data) {
849  this->used_fields_ += eq_fields_->advection_coef;
850  this->used_fields_ += eq_fields_->diffusion_coef;
851  this->used_fields_ += eq_fields_->cross_section;
852  this->used_fields_ += eq_fields_->bc_type;
853  this->used_fields_ += eq_fields_->bc_dirichlet_value;
854  this->used_fields_ += eq_fields_->bc_robin_sigma;
855  this->used_fields_ += eq_fields_->bc_flux;
856  }
857 
858  /// Destructor.
860 
861  /// Initialize auxiliary vectors and other data members
862  void initialize(ElementCacheMap *element_cache_map) {
863  this->element_cache_map_ = element_cache_map;
864 
865  fe_ = std::make_shared< FE_P_disc<dim> >(eq_data_->dg_order);
868  if (dim==1) // print to log only one time
869  DebugOut() << "List of BdrConditionAssemblyDG FEValues updates flags: " << this->print_update_flags(u);
870  ndofs_ = fe_->n_dofs();
871  dof_indices_.resize(ndofs_);
872  local_rhs_.resize(ndofs_);
874  }
875 
876 
877  /// Implements @p AssemblyBase::boundary_side_integral.
878  inline void boundary_side_integral(DHCellSide cell_side)
879  {
880  unsigned int k;
881 
882  ElementAccessor<3> bc_elm = cell_side.cond().element_accessor();
883 
884  fe_values_side_.reinit(cell_side.side());
885 
886  const DHCellAccessor &cell = cell_side.cell();
888 
889  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
890  {
891  fill_n(&(local_rhs_[0]), ndofs_, 0);
894 
895  double side_flux = advective_flux(eq_fields_->advection_coef[sbi], this->boundary_points(cell_side), fe_values_side_);
896  double transport_flux = side_flux/cell_side.measure();
897 
898  auto p_side = *( this->boundary_points(cell_side)).begin();
899  auto p_bdr = p_side.point_bdr(cell_side.cond().element_accessor() );
900  unsigned int bc_type = eq_fields_->bc_type[sbi](p_bdr);
901  if (bc_type == AdvectionDiffusionModel::abc_inflow && side_flux < 0)
902  {
903  k=0;
904  for (auto p : this->boundary_points(cell_side) )
905  {
906  auto p_bdr = p.point_bdr(bc_elm);
907  double bc_term = -transport_flux*eq_fields_->bc_dirichlet_value[sbi](p_bdr)*fe_values_side_.JxW(k);
908  for (unsigned int i=0; i<ndofs_; i++)
909  local_rhs_[i] += bc_term*fe_values_side_.shape_value(i,k);
910  k++;
911  }
912  for (unsigned int i=0; i<ndofs_; i++)
914  }
915  else if (bc_type == AdvectionDiffusionModel::abc_dirichlet)
916  {
917  double side_flux = advective_flux(eq_fields_->advection_coef[sbi], this->boundary_points(cell_side), fe_values_side_);
918  double transport_flux = side_flux/cell_side.measure();
919 
920  double gamma_l = DG_penalty_boundary(cell_side.side(),
921  diffusion_delta(eq_fields_->diffusion_coef[sbi], this->boundary_points(cell_side), fe_values_side_.normal_vector(0)),
922  transport_flux,
923  eq_fields_->dg_penalty[sbi](p_bdr));
924 
925  k=0;
926  for (auto p : this->boundary_points(cell_side) )
927  {
928  auto p_bdr = p.point_bdr(bc_elm);
929  double bc_term = gamma_l*eq_fields_->bc_dirichlet_value[sbi](p_bdr)*fe_values_side_.JxW(k);
930  arma::vec3 bc_grad = -eq_fields_->bc_dirichlet_value[sbi](p_bdr)*fe_values_side_.JxW(k)*eq_data_->dg_variant*(arma::trans(eq_fields_->diffusion_coef[sbi](p))*fe_values_side_.normal_vector(k));
931  for (unsigned int i=0; i<ndofs_; i++)
932  local_rhs_[i] += bc_term*fe_values_side_.shape_value(i,k)
933  + arma::dot(bc_grad,fe_values_side_.shape_grad(i,k));
934  k++;
935  }
936  k=0;
937  for (auto p : this->boundary_points(cell_side) )
938  {
939  for (unsigned int i=0; i<ndofs_; i++)
940  {
941  local_flux_balance_vector_[i] += (arma::dot(eq_fields_->advection_coef[sbi](p), fe_values_side_.normal_vector(k))*fe_values_side_.shape_value(i,k)
942  - arma::dot(eq_fields_->diffusion_coef[sbi](p)*fe_values_side_.shape_grad(i,k),fe_values_side_.normal_vector(k))
943  + gamma_l*fe_values_side_.shape_value(i,k))*fe_values_side_.JxW(k);
944  }
945  k++;
946  }
947  if (eq_data_->time_->step().index() > 0)
948  for (unsigned int i=0; i<ndofs_; i++)
950  }
951  else if (bc_type == AdvectionDiffusionModel::abc_total_flux)
952  {
953  k=0;
954  for (auto p : this->boundary_points(cell_side) )
955  {
956  auto p_bdr = p.point_bdr(bc_elm);
957  double bc_term = eq_fields_->cross_section(p) * (eq_fields_->bc_robin_sigma[sbi](p_bdr)*eq_fields_->bc_dirichlet_value[sbi](p_bdr) +
958  eq_fields_->bc_flux[sbi](p_bdr)) * fe_values_side_.JxW(k);
959  for (unsigned int i=0; i<ndofs_; i++)
960  local_rhs_[i] += bc_term*fe_values_side_.shape_value(i,k);
961  k++;
962  }
963 
964  for (unsigned int i=0; i<ndofs_; i++)
965  {
966  k=0;
967  for (auto p : this->boundary_points(cell_side) ) {
968  auto p_bdr = p.point_bdr(bc_elm);
969  local_flux_balance_vector_[i] += eq_fields_->cross_section(p) * eq_fields_->bc_robin_sigma[sbi](p_bdr) *
971  k++;
972  }
974  }
975  }
977  {
978  k=0;
979  for (auto p : this->boundary_points(cell_side) )
980  {
981  auto p_bdr = p.point_bdr(bc_elm);
982  double bc_term = eq_fields_->cross_section(p) * (eq_fields_->bc_robin_sigma[sbi](p_bdr)*eq_fields_->bc_dirichlet_value[sbi](p_bdr) +
983  eq_fields_->bc_flux[sbi](p_bdr)) * fe_values_side_.JxW(k);
984  for (unsigned int i=0; i<ndofs_; i++)
985  local_rhs_[i] += bc_term*fe_values_side_.shape_value(i,k);
986  k++;
987  }
988 
989  for (unsigned int i=0; i<ndofs_; i++)
990  {
991  k=0;
992  for (auto p : this->boundary_points(cell_side) ) {
993  auto p_bdr = p.point_bdr(bc_elm);
994  local_flux_balance_vector_[i] += eq_fields_->cross_section(p)*(arma::dot(eq_fields_->advection_coef[sbi](p), fe_values_side_.normal_vector(k)) +
995  eq_fields_->bc_robin_sigma[sbi](p_bdr))*fe_values_side_.JxW(k)*fe_values_side_.shape_value(i,k);
996  k++;
997  }
999  }
1000  }
1001  else if (bc_type == AdvectionDiffusionModel::abc_inflow && side_flux >= 0)
1002  {
1003  k=0;
1004  for (auto p : this->boundary_points(cell_side) )
1005  {
1006  for (unsigned int i=0; i<ndofs_; i++)
1008  k++;
1009  }
1010  }
1011  eq_data_->ls[sbi]->rhs_set_values(ndofs_, &(dof_indices_[0]), &(local_rhs_[0]));
1012 
1013  eq_data_->balance_->add_flux_values(eq_data_->subst_idx()[sbi], cell_side,
1014  cell.get_loc_dof_indices(),
1016  }
1017  }
1018 
1019  /// Implements @p AssemblyBase::begin.
1020  void begin() override
1021  {
1022  eq_data_->balance_->start_flux_assembly( eq_data_->subst_idx() );
1023  }
1024 
1025  /// Implements @p AssemblyBase::end.
1026  void end() override
1027  {
1028  eq_data_->balance_->finish_flux_assembly( eq_data_->subst_idx() );
1029  }
1030 
1031 
1032  private:
1033  shared_ptr<FiniteElement<dim>> fe_; ///< Finite element for the solution of the advection-diffusion equation.
1034 
1035  /// Data objects shared with TransportDG
1038 
1039  /// Sub field set contains fields used in calculation.
1041 
1042  unsigned int ndofs_; ///< Number of dofs
1043  FEValues<3> fe_values_side_; ///< FEValues of object (of P disc finite element type)
1044 
1045  vector<LongIdx> dof_indices_; ///< Vector of global DOF indices
1046  vector<PetscScalar> local_rhs_; ///< Auxiliary vector for set_sources method.
1047  vector<PetscScalar> local_flux_balance_vector_; ///< Auxiliary vector for set_boundary_conditions method.
1048  PetscScalar local_flux_balance_rhs_; ///< Auxiliary variable for set_boundary_conditions method.
1049 
1050  template < template<IntDim...> class DimAssembly>
1051  friend class GenericAssembly;
1052 
1053 };
1054 
1055 
1056 /**
1057  * Auxiliary container class sets the initial condition.
1058  */
1059 template <unsigned int dim, class Model>
1061 {
1062 public:
1065 
1066  static constexpr const char * name() { return "InitProjectionAssemblyDG"; }
1067 
1068  /// Constructor.
1070  : AssemblyBase<dim>(eq_data->dg_order), eq_fields_(eq_fields), eq_data_(eq_data) {
1072  this->used_fields_ += eq_fields_->init_condition;
1073  }
1074 
1075  /// Destructor.
1077 
1078  /// Initialize auxiliary vectors and other data members
1079  void initialize(ElementCacheMap *element_cache_map) {
1080  this->element_cache_map_ = element_cache_map;
1081 
1082  fe_ = std::make_shared< FE_P_disc<dim> >(eq_data_->dg_order);
1084  fe_values_.initialize(*this->quad_, *fe_, u);
1085  // if (dim==1) // print to log only one time
1086  // DebugOut() << "List of InitProjectionAssemblyDG FEValues updates flags: " << this->print_update_flags(u);
1087  ndofs_ = fe_->n_dofs();
1088  dof_indices_.resize(ndofs_);
1089  local_matrix_.resize(4*ndofs_*ndofs_);
1090  local_rhs_.resize(ndofs_);
1091  }
1092 
1093 
1094  /// Assemble integral over element
1095  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
1096  {
1097  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
1098 
1099  unsigned int k;
1100  ElementAccessor<3> elem = cell.elm();
1102  fe_values_.reinit(elem);
1103 
1104  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
1105  {
1106  for (unsigned int i=0; i<ndofs_; i++)
1107  {
1108  local_rhs_[i] = 0;
1109  for (unsigned int j=0; j<ndofs_; j++)
1110  local_matrix_[i*ndofs_+j] = 0;
1111  }
1112 
1113  k=0;
1114  for (auto p : this->bulk_points(element_patch_idx) )
1115  {
1116  double rhs_term = eq_fields_->init_condition[sbi](p)*fe_values_.JxW(k);
1117 
1118  for (unsigned int i=0; i<ndofs_; i++)
1119  {
1120  for (unsigned int j=0; j<ndofs_; j++)
1122 
1123  local_rhs_[i] += fe_values_.shape_value(i,k)*rhs_term;
1124  }
1125  k++;
1126  }
1127  eq_data_->ls[sbi]->set_values(ndofs_, &(dof_indices_[0]), ndofs_, &(dof_indices_[0]), &(local_matrix_[0]), &(local_rhs_[0]));
1128  }
1129  }
1130 
1131 
1132  private:
1133  shared_ptr<FiniteElement<dim>> fe_; ///< Finite element for the solution of the advection-diffusion equation.
1134 
1135  /// Data objects shared with TransportDG
1138 
1139  /// Sub field set contains fields used in calculation.
1141 
1142  unsigned int ndofs_; ///< Number of dofs
1143  FEValues<3> fe_values_; ///< FEValues of object (of P disc finite element type)
1144 
1145  vector<LongIdx> dof_indices_; ///< Vector of global DOF indices
1146  vector<PetscScalar> local_matrix_; ///< Auxiliary vector for assemble methods
1147  vector<PetscScalar> local_rhs_; ///< Auxiliary vector for set_sources method.
1148 
1149  template < template<IntDim...> class DimAssembly>
1150  friend class GenericAssembly;
1151 
1152 };
1153 
1154 
1155 /**
1156  * Auxiliary container class sets the initial condition.
1157  */
1158 template <unsigned int dim, class Model>
1160 {
1161 public:
1164 
1165  static constexpr const char * name() { return "InitConditionAssemblyDG"; }
1166 
1167  /// Constructor.
1169  : AssemblyBase<dim>(), eq_fields_(eq_fields), eq_data_(eq_data) {
1171  this->used_fields_ += eq_fields_->init_condition;
1172 
1173  this->quad_ = new Quadrature(dim, RefElement<dim>::n_nodes);
1174  for(unsigned int i = 0; i<RefElement<dim>::n_nodes; i++)
1175  {
1176  this->quad_->weight(i) = 1.0;
1177  this->quad_->set(i) = RefElement<dim>::node_coords(i);
1178  }
1179  }
1180 
1181  /// Destructor.
1183 
1184  /// Initialize auxiliary vectors and other data members
1185  void initialize(ElementCacheMap *element_cache_map) {
1186  this->element_cache_map_ = element_cache_map;
1187  }
1188 
1189 
1190  /// Assemble integral over element
1191  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
1192  {
1193  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
1194 
1195  unsigned int k;
1196 
1197  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
1198  {
1199  k=0;
1200  for (auto p : this->bulk_points(element_patch_idx) )
1201  {
1202  double val = eq_fields_->init_condition[sbi](p);
1203 
1204  eq_data_->output_vec[sbi].set(cell.get_loc_dof_indices()[k], val);
1205  k++;
1206  }
1207  }
1208  }
1209 
1210 
1211  private:
1212  /// Data objects shared with TransportDG
1215 
1216  /// Sub field set contains fields used in calculation.
1218 
1219  template < template<IntDim...> class DimAssembly>
1220  friend class GenericAssembly;
1221 
1222 };
1223 
1224 #endif /* ASSEMBLY_DG_HH_ */
1225 
double elem_anisotropy(ElementAccessor< 3 > e)
Definition: assembly_dg.hh:160
double DG_penalty_boundary(Side side, const double &diff_delta, const double flux, const double alpha)
Computes the penalty parameter of the DG method on a given boundary edge.
Definition: assembly_dg.hh:204
double diffusion_delta(Field< 3, FieldValue< 3 >::TensorFixed > &diff_coef, Range< BoundaryPoint > pts, const arma::vec3 &nv)
Computes average normal diffusivity over a set of points.
Definition: assembly_dg.hh:181
double advective_flux(Field< 3, FieldValue< 3 >::VectorFixed > &advection_coef, Range< PointType > pts, FEValues< 3 > &fv)
Computes advective flux.
Definition: assembly_dg.hh:222
#define ASSERT(expr)
Definition: asserts.hh:351
#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.
Range< EdgePoint > edge_points(const DHCellSide &cell_side) const
Return EdgePoint range of appropriate dimension.
Quadrature * quad_
Quadrature used in assembling methods.
std::string print_update_flags(UpdateFlags u) const
Print update flags to string format.
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.
static constexpr const char * name()
Definition: assembly_dg.hh:843
FieldSet used_fields_
Sub field set contains fields used in calculation.
vector< PetscScalar > local_flux_balance_vector_
Auxiliary vector for set_boundary_conditions method.
void boundary_side_integral(DHCellSide cell_side)
Implements AssemblyBase::boundary_side_integral.
Definition: assembly_dg.hh:878
FEValues< 3 > fe_values_side_
FEValues of object (of P disc finite element type)
TransportDG< Model >::EqData EqData
Definition: assembly_dg.hh:841
TransportDG< Model >::EqFields EqFields
Definition: assembly_dg.hh:840
void begin() override
Implements AssemblyBase::begin.
~BdrConditionAssemblyDG()
Destructor.
Definition: assembly_dg.hh:859
PetscScalar local_flux_balance_rhs_
Auxiliary variable for set_boundary_conditions method.
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
Definition: assembly_dg.hh:862
vector< LongIdx > dof_indices_
Vector of global DOF indices.
EqFields * eq_fields_
Data objects shared with TransportDG.
void end() override
Implements AssemblyBase::end.
vector< PetscScalar > local_rhs_
Auxiliary vector for set_sources method.
unsigned int ndofs_
Number of dofs.
BdrConditionAssemblyDG(EqFields *eq_fields, EqData *eq_data)
Constructor.
Definition: assembly_dg.hh:846
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_.
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.
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.
Side accessor allows to iterate over sides of DOF handler cell.
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
ElementAccessor< 3 > element() const
NodeAccessor< 3 > node(unsigned int ni) const
Definition: accessors.hh:230
Region region() const
Definition: accessors.hh:198
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
Definition: accessors.hh:215
Directing class of FieldValueCache.
unsigned int n_nodes() const
Definition: elements.h:125
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
double shape_value(const unsigned int function_no, const unsigned int point_no) const
Return the value of the function_no-th shape function at the point_no-th quadrature point.
Definition: fe_values.hh:149
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
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
arma::vec::fixed< spacedim > shape_grad(const unsigned int function_no, const unsigned int point_no) const
Return the gradient of the function_no-th shape function at the point_no-th quadrature point.
Definition: fe_values.hh:164
Discontinuous Lagrangean finite element on dim dimensional simplex.
Definition: fe_p.hh:108
Container for various descendants of FieldCommonBase.
Definition: field_set.hh:159
Class template representing a field with values dependent on: point, element, and region.
Definition: field.hh:92
Abstract class for the description of a general finite element on a reference simplex in dim dimensio...
Generic class of assemblation.
FieldSet used_fields_
Sub field set contains fields used in calculation.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
static constexpr const char * name()
~InitConditionAssemblyDG()
Destructor.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
InitConditionAssemblyDG(EqFields *eq_fields, EqData *eq_data)
Constructor.
EqFields * eq_fields_
Data objects shared with TransportDG.
TransportDG< Model >::EqData EqData
TransportDG< Model >::EqFields EqFields
EqFields * eq_fields_
Data objects shared with TransportDG.
vector< PetscScalar > local_rhs_
Auxiliary vector for set_sources method.
vector< LongIdx > dof_indices_
Vector of global DOF indices.
TransportDG< Model >::EqData EqData
FieldSet used_fields_
Sub field set contains fields used in calculation.
unsigned int ndofs_
Number of dofs.
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
FEValues< 3 > fe_values_
FEValues of object (of P disc finite element type)
static constexpr const char * name()
vector< PetscScalar > local_matrix_
Auxiliary vector for assemble methods.
InitProjectionAssemblyDG(EqFields *eq_fields, EqData *eq_data)
Constructor.
~InitProjectionAssemblyDG()
Destructor.
TransportDG< Model >::EqFields EqFields
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
virtual void rhs_set_values(int nrow, int *rows, double *vals)=0
virtual void mat_set_values(int nrow, int *rows, int ncol, int *cols, double *vals)=0
void set_values(int nrow, int *rows, int ncol, int *cols, PetscScalar *mat_vals, PetscScalar *rhs_vals)
Set values in the system matrix and values in the right-hand side vector on corresponding rows.
Definition: linsys.hh:386
static constexpr double almost_zero
vector< LongIdx > dof_indices_
Vector of global DOF indices.
Definition: assembly_dg.hh:148
vector< PetscScalar > local_mass_balance_vector_
Same as previous.
Definition: assembly_dg.hh:151
TransportDG< Model >::EqFields EqFields
Definition: assembly_dg.hh:39
FieldSet used_fields_
Sub field set contains fields used in calculation.
Definition: assembly_dg.hh:143
FEValues< 3 > fe_values_
FEValues of object (of P disc finite element type)
Definition: assembly_dg.hh:146
vector< PetscScalar > local_matrix_
Auxiliary vector for assemble methods.
Definition: assembly_dg.hh:149
void end() override
Implements AssemblyBase::end.
Definition: assembly_dg.hh:130
~MassAssemblyDG()
Destructor.
Definition: assembly_dg.hh:53
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
Definition: assembly_dg.hh:136
EqData * eq_data_
Definition: assembly_dg.hh:140
void begin() override
Implements AssemblyBase::begin.
Definition: assembly_dg.hh:124
static constexpr const char * name()
Definition: assembly_dg.hh:42
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
Definition: assembly_dg.hh:56
unsigned int ndofs_
Number of dofs.
Definition: assembly_dg.hh:145
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
Definition: assembly_dg.hh:74
TransportDG< Model >::EqData EqData
Definition: assembly_dg.hh:40
EqFields * eq_fields_
Data objects shared with TransportDG.
Definition: assembly_dg.hh:139
vector< PetscScalar > local_retardation_balance_vector_
Auxiliary vector for assemble mass matrix.
Definition: assembly_dg.hh:150
MassAssemblyDG(EqFields *eq_fields, EqData *eq_data)
Constructor.
Definition: assembly_dg.hh:45
Base class for quadrature rules on simplices in arbitrary dimensions.
Definition: quadrature.hh:48
unsigned int size() const
Returns number of quadrature points.
Definition: quadrature.hh:86
Armor::Array< double >::ArrayMatSet set(uint i)
Definition: quadrature.hh:93
double weight(unsigned int i) const
Returns the ith weight.
Definition: quadrature.hh:103
IterConvert< ObjectIn, ObjectOut > begin()
Iterator to begin item of range.
Range helper class.
static LocalPoint node_coords(unsigned int nid)
Definition: ref_element.hh:619
unsigned int bulk_idx() const
Returns index of the region in the bulk set.
Definition: region.hh:90
double measure() const
Calculate metrics of the side.
Definition: accessors.cc:27
Boundary cond() const
ElementAccessor< 3 > element() const
Returns iterator to the element of the side.
double diameter() const
Calculate the side diameter.
Definition: accessors.cc:132
vector< PetscScalar > local_rhs_
Auxiliary vector for set_sources method.
Definition: assembly_dg.hh:823
TransportDG< Model >::EqData EqData
Definition: assembly_dg.hh:716
TransportDG< Model >::EqFields EqFields
Definition: assembly_dg.hh:715
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
Definition: assembly_dg.hh:750
~SourcesAssemblyDG()
Destructor.
Definition: assembly_dg.hh:730
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
Definition: assembly_dg.hh:810
vector< PetscScalar > local_source_balance_vector_
Auxiliary vector for set_sources method.
Definition: assembly_dg.hh:824
void end() override
Implements AssemblyBase::end.
Definition: assembly_dg.hh:803
FEValues< 3 > fe_values_
FEValues of object (of P disc finite element type)
Definition: assembly_dg.hh:820
static constexpr const char * name()
Definition: assembly_dg.hh:718
EqFields * eq_fields_
Data objects shared with TransportDG.
Definition: assembly_dg.hh:813
vector< PetscScalar > local_source_balance_rhs_
Auxiliary vector for set_sources method.
Definition: assembly_dg.hh:825
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
Definition: assembly_dg.hh:733
FieldSet used_fields_
Sub field set contains fields used in calculation.
Definition: assembly_dg.hh:817
unsigned int ndofs_
Number of dofs.
Definition: assembly_dg.hh:819
vector< LongIdx > dof_indices_
Vector of global DOF indices.
Definition: assembly_dg.hh:822
void begin() override
Implements AssemblyBase::begin.
Definition: assembly_dg.hh:797
SourcesAssemblyDG(EqFields *eq_fields, EqData *eq_data)
Constructor.
Definition: assembly_dg.hh:721
TransportDG< Model >::EqData EqData
Definition: assembly_dg.hh:242
void edge_integral(RangeConvert< DHEdgeSide, DHCellSide > edge_side_range)
Assembles the fluxes between sides of elements of the same dimension.
Definition: assembly_dg.hh:436
StiffnessAssemblyDG(EqFields *eq_fields, EqData *eq_data)
Constructor.
Definition: assembly_dg.hh:247
FieldSet used_fields_
Sub field set contains fields used in calculation.
Definition: assembly_dg.hh:683
unsigned int ndofs_
Number of dofs.
Definition: assembly_dg.hh:685
vector< vector< LongIdx > > side_dof_indices_
Vector of vectors of side DOF indices.
Definition: assembly_dg.hh:694
shared_ptr< FiniteElement< dim-1 > > fe_low_
Finite element for the solution of the advection-diffusion equation (dim-1).
Definition: assembly_dg.hh:676
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
Definition: assembly_dg.hh:675
void boundary_side_integral(DHCellSide cell_side)
Assembles the fluxes on the boundary.
Definition: assembly_dg.hh:357
vector< double * > jumps
Auxiliary storage for jumps of shape functions.
Definition: assembly_dg.hh:700
vector< double * > waverages
Auxiliary storage for weighted averages of shape functions.
Definition: assembly_dg.hh:699
vector< FEValues< 3 > * > fv_sb_
Auxiliary vector, holds FEValues objects for assemble element-side.
Definition: assembly_dg.hh:691
void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side)
Assembles the fluxes between elements of different dimensions.
Definition: assembly_dg.hh:599
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
Definition: assembly_dg.hh:268
vector< LongIdx > side_dof_indices_vb_
Vector of side DOF indices (assemble element-side fluxex)
Definition: assembly_dg.hh:695
FEValues< 3 > fe_values_side_
FEValues of object (of P disc finite element type)
Definition: assembly_dg.hh:689
vector< FEValues< 3 > > fe_values_vec_
Vector of FEValues of object (of P disc finite element types)
Definition: assembly_dg.hh:690
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assembles the cell (volume) integral into the stiffness matrix.
Definition: assembly_dg.hh:318
TransportDG< Model >::EqFields EqFields
Definition: assembly_dg.hh:241
vector< double * > averages
Auxiliary storage for averages of shape functions.
Definition: assembly_dg.hh:698
vector< PetscScalar > local_matrix_
Auxiliary vector for assemble methods.
Definition: assembly_dg.hh:696
~StiffnessAssemblyDG()
Destructor.
Definition: assembly_dg.hh:261
static constexpr const char * name()
Definition: assembly_dg.hh:244
vector< LongIdx > dof_indices_
Vector of global DOF indices.
Definition: assembly_dg.hh:693
FEValues< 3 > fe_values_vb_
FEValues of dim-1 object (of P disc finite element type)
Definition: assembly_dg.hh:688
FEValues< 3 > fe_values_
FEValues of object (of P disc finite element type)
Definition: assembly_dg.hh:687
EqFields * eq_fields_
Data objects shared with TransportDG.
Definition: assembly_dg.hh:679
unsigned int qsize_lower_dim_
Size of quadrature of dim-1.
Definition: assembly_dg.hh:686
const TimeStep & step(int index=-1) const
unsigned int index() const
std::shared_ptr< DOFHandlerMultiDim > dh_
Object for distribution of dofs.
std::vector< VectorMPI > output_vec
Vector of solution data.
LinSys ** ls_dt
Linear algebra system for the time derivative (actually it is used only for handling the matrix struc...
int dg_variant
DG variant ((non-)symmetric/incomplete.
TimeGovernor * time_
std::shared_ptr< Balance > balance_
unsigned int dg_order
Polynomial order of finite elements.
LinSys ** ls
Linear algebra system for the transport equation.
std::vector< Vec > ret_vec
Auxiliary vectors for calculation of sources in balance due to retardation (e.g. sorption).
unsigned int max_edg_sides
Maximal number of edge sides (evaluate from dim 1,2,3)
MultiField< 3, FieldValue< 3 >::Scalar > dg_penalty
Penalty enforcing inter-element continuity of solution (for each substance).
MultiField< 3, FieldValue< 3 >::Scalar > fracture_sigma
Transition parameter for diffusive transfer on fractures (for each substance).
Definitions of basic Lagrangean finite elements with polynomial shape functions.
Class FEValues calculates finite element data on the actual cells such as shape function values,...
@ coupling
@ boundary
@ bulk
@ edge
#define DebugOut()
Macro defining 'debug' record of log.
Definition: logger.hh:284
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19
Definitions of particular quadrature rules on simplices.
Discontinuous Galerkin method for equation of transport with dispersion.
UpdateFlags
Enum type UpdateFlags indicates which quantities are to be recomputed on each finite element cell.
Definition: update_flags.hh:68
@ update_values
Shape function values.
Definition: update_flags.hh:87
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_side_JxW_values
Transformed quadrature weight for cell sides.
@ update_gradients
Shape function gradients.
Definition: update_flags.hh:95
@ update_quadrature_points
Transformed quadrature points.