Flow123d
JS_before_hm-1804-gf2ad740aa
|
Go to the documentation of this file.
18 #ifndef ASSEMBLY_DG_HH_
19 #define ASSEMBLY_DG_HH_
35 template <
unsigned int dim,
class Model>
42 static constexpr
const char *
name() {
return "MassAssemblyDG"; }
84 for (
unsigned int sbi=0; sbi<
eq_data_->n_substances(); ++sbi)
87 for (
unsigned int i=0; i<
ndofs_; i++)
89 for (
unsigned int j=0; j<
ndofs_; j++)
93 for (
auto p : this->
bulk_points(element_patch_idx) )
102 for (
unsigned int i=0; i<
ndofs_; i++)
107 for (
auto p : this->
bulk_points(element_patch_idx) )
136 shared_ptr<FiniteElement<dim>>
fe_;
153 template <
template<
IntDim...>
class DimAssembly>
162 template <
unsigned int dim,
class Model>
169 static constexpr
const char *
name() {
return "StiffnessAssemblyDG"; }
187 for (
auto a :
averages)
if (a !=
nullptr)
delete[] a;
188 for (
auto a :
waverages)
if (a !=
nullptr)
delete[] a;
189 for (
auto a :
jumps)
if (a !=
nullptr)
delete[] a;
234 for (
unsigned int s=0; s<2; s++)
246 if (!cell.
is_own())
return;
255 for (
unsigned int sbi=0; sbi<
eq_data_->n_substances(); sbi++)
257 for (
unsigned int i=0; i<
ndofs_; i++)
258 for (
unsigned int j=0; j<
ndofs_; j++)
262 for (
auto p : this->
bulk_points(element_patch_idx) )
264 for (
unsigned int i=0; i<
ndofs_; i++)
269 for (
unsigned int j=0; j<
ndofs_; j++)
284 ASSERT_EQ_DBG(cell_side.
dim(), dim).error(
"Dimension of element mismatch!");
295 for (
unsigned int sbi=0; sbi<
eq_data_->n_substances(); sbi++)
301 double side_flux = 0;
307 double transport_flux = side_flux/side.
measure();
311 unsigned int bc_type =
eq_fields_->bc_type[sbi](p_bdr);
322 transport_flux += gamma_l;
329 double flux_times_JxW;
346 for (
unsigned int i=0; i<
ndofs_; i++)
348 for (
unsigned int j=0; j<
ndofs_; j++)
373 double gamma_l, omega[2], transport_flux, delta[2], delta_sum;
374 double aniso1, aniso2;
375 double local_alpha=0.0;
379 auto dh_edge_cell =
eq_data_->
dh_->cell_accessor_from_element( edge_side.elem_idx() );
387 for (
unsigned int sbi=0; sbi<
eq_data_->n_substances(); sbi++)
390 double pflux = 0, nflux = 0;
400 fluxes[sid] /= edge_side.measure();
402 pflux += fluxes[sid];
404 nflux += fluxes[sid];
415 for (
unsigned int i=0; i<
fe_->n_dofs(); i++)
424 for(
DHCellSide edge_side1 : edge_side_range )
427 for(
DHCellSide edge_side2 : edge_side_range )
430 if (s2<=s1)
continue;
431 ASSERT(edge_side1.is_valid()).error(
"Invalid side of edge.");
437 if (fluxes[s2] > 0 && fluxes[s1] < 0)
438 transport_flux = fluxes[s1]*fabs(fluxes[s2]/pflux);
439 else if (fluxes[s2] < 0 && fluxes[s1] > 0)
440 transport_flux = fluxes[s1]*fabs(fluxes[s2]/nflux);
444 gamma_l = 0.5*fabs(transport_flux);
450 auto p2 = p1.point_on(edge_side2);
451 delta[0] += dot(
eq_fields_->diffusion_coef[sbi](p1)*normal_vector,normal_vector);
452 delta[1] += dot(
eq_fields_->diffusion_coef[sbi](p2)*normal_vector,normal_vector);
458 delta_sum = delta[0] + delta[1];
461 if (fabs(delta_sum) > 0)
463 omega[0] = delta[1]/delta_sum;
464 omega[1] = delta[0]/delta_sum;
465 double h = edge_side1.diameter();
468 gamma_l += local_alpha/h*aniso1*aniso2*(delta[0]*delta[1]/delta_sum);
471 for (
int i=0; i<2; i++) omega[i] = 0;
474 int sd[2];
bool is_side_own[2];
475 sd[0] = s1; is_side_own[0] = edge_side1.cell().is_own();
476 sd[1] = s2; is_side_own[1] = edge_side2.cell().is_own();
482 auto p2 = p1.point_on(edge_side2);
483 for (
unsigned int i=0; i<
fe_->n_dofs(); i++)
485 for (
int n=0; n<2; n++)
495 for (
int n=0; n<2; n++)
497 if (!is_side_own[n])
continue;
499 for (
int m=0; m<2; m++)
539 if (dim == 1)
return;
540 ASSERT_EQ_DBG(cell_lower_dim.
dim(), dim-1).error(
"Dimension of element mismatch!");
544 double comm_flux[2][2];
545 unsigned int n_dofs[2];
548 for(
unsigned int i=0; i<n_indices; ++i) {
552 n_dofs[0] =
fv_sb_[0]->n_dofs();
556 for(
unsigned int i=0; i<n_indices; ++i) {
560 n_dofs[1] =
fv_sb_[1]->n_dofs();
563 bool own_element_id[2];
564 own_element_id[0] = cell_lower_dim.
is_own();
565 own_element_id[1] = cell_higher_dim.
is_own();
568 for (
unsigned int sbi=0; sbi<
eq_data_->n_substances(); sbi++)
570 for (
unsigned int i=0; i<n_dofs[0]+n_dofs[1]; i++)
571 for (
unsigned int j=0; j<n_dofs[0]+n_dofs[1]; j++)
578 auto p_low = p_high.lower_dim(cell_lower_dim);
591 comm_flux[0][0] = (sigma-min(0.,transport_flux))*
fv_sb_[0]->JxW(k);
592 comm_flux[0][1] = -(sigma-min(0.,transport_flux))*
fv_sb_[0]->JxW(k);
593 comm_flux[1][0] = -(sigma+max(0.,transport_flux))*
fv_sb_[0]->JxW(k);
594 comm_flux[1][1] = (sigma+max(0.,transport_flux))*
fv_sb_[0]->JxW(k);
596 for (
int n=0; n<2; n++)
598 if (!own_element_id[n])
continue;
600 for (
unsigned int i=0; i<n_dofs[n]; i++)
601 for (
int m=0; m<2; m++)
602 for (
unsigned int j=0; j<n_dofs[m]; j++)
603 local_matrix_[(i+n*n_dofs[0])*(n_dofs[0]+n_dofs[1]) + m*n_dofs[0] + j] +=
614 shared_ptr<FiniteElement<dim>>
fe_;
641 template <
template<
IntDim...>
class DimAssembly>
650 template <
unsigned int dim,
class Model>
657 static constexpr
const char *
name() {
return "SourcesAssemblyDG"; }
701 for (
unsigned int sbi=0; sbi<
eq_data_->n_substances(); sbi++)
708 for (
auto p : this->
bulk_points(element_patch_idx) )
712 for (
unsigned int i=0; i<
ndofs_; i++)
718 for (
unsigned int i=0; i<
ndofs_; i++)
721 for (
auto p : this->
bulk_points(element_patch_idx) )
749 shared_ptr<FiniteElement<dim>>
fe_;
766 template <
template<
IntDim...>
class DimAssembly>
775 template <
unsigned int dim,
class Model>
782 static constexpr
const char *
name() {
return "BdrConditionAssemblyDG"; }
821 const unsigned int cond_idx = cell_side.
side().
cond_idx();
830 for (
unsigned int sbi=0; sbi<
eq_data_->n_substances(); sbi++)
836 double side_flux = 0;
842 double transport_flux = side_flux/cell_side.
measure();
846 unsigned int bc_type =
eq_fields_->bc_type[sbi](p_bdr);
852 auto p_bdr = p.point_bdr(bc_elm);
854 for (
unsigned int i=0; i<
ndofs_; i++)
858 for (
unsigned int i=0; i<
ndofs_; i++)
866 auto p_bdr = p.point_bdr(bc_elm);
869 for (
unsigned int i=0; i<
ndofs_; i++)
877 for (
unsigned int i=0; i<
ndofs_; i++)
886 for (
unsigned int i=0; i<
ndofs_; i++)
894 auto p_bdr = p.point_bdr(bc_elm);
897 for (
unsigned int i=0; i<
ndofs_; i++)
902 for (
unsigned int i=0; i<
ndofs_; i++)
906 auto p_bdr = p.point_bdr(bc_elm);
919 auto p_bdr = p.point_bdr(bc_elm);
922 for (
unsigned int i=0; i<
ndofs_; i++)
927 for (
unsigned int i=0; i<
ndofs_; i++)
931 auto p_bdr = p.point_bdr(bc_elm);
944 for (
unsigned int i=0; i<
ndofs_; i++)
971 shared_ptr<FiniteElement<dim>>
fe_;
988 template <
template<
IntDim...>
class DimAssembly>
997 template <
unsigned int dim,
class Model>
1004 static constexpr
const char *
name() {
return "InitProjectionAssemblyDG"; }
1042 for (
unsigned int sbi=0; sbi<
eq_data_->n_substances(); sbi++)
1044 for (
unsigned int i=0; i<
ndofs_; i++)
1047 for (
unsigned int j=0; j<
ndofs_; j++)
1052 for (
auto p : this->
bulk_points(element_patch_idx) )
1056 for (
unsigned int i=0; i<
ndofs_; i++)
1058 for (
unsigned int j=0; j<
ndofs_; j++)
1071 shared_ptr<FiniteElement<dim>>
fe_;
1087 template <
template<
IntDim...>
class DimAssembly>
1096 template <
unsigned int dim,
class Model>
1103 static constexpr
const char *
name() {
return "InitConditionAssemblyDG"; }
1112 for(
unsigned int i = 0; i<RefElement<dim>::n_nodes; i++)
1135 for (
unsigned int sbi=0; sbi<
eq_data_->n_substances(); sbi++)
1138 for (
auto p : this->
bulk_points(element_patch_idx) )
1140 double val =
eq_fields_->init_condition[sbi](p);
1157 template <
template<
IntDim...>
class DimAssembly>
TransportDG< Model >::EqFields EqFields
~InitConditionAssemblyDG()
Destructor.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
unsigned int dg_order
Polynomial order of finite elements.
std::string print_update_flags(UpdateFlags u) const
Print update flags to string format.
vector< PetscScalar > local_flux_balance_vector_
Auxiliary vector for set_boundary_conditions method.
FieldSet used_fields_
Sub field set contains fields used in calculation.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
vector< PetscScalar > local_matrix_
Auxiliary vector for assemble methods.
@ update_gradients
Shape function gradients.
TransportDG< Model >::EqFields EqFields
EqFields * eq_fields_
Data objects shared with TransportDG.
void begin() override
Implements AssemblyBase::begin.
~StiffnessAssemblyDG()
Destructor.
SourcesAssemblyDG(EqFields *eq_fields, EqData *eq_data)
Constructor.
Side side() const
Return Side of given cell and side_idx.
MultiField< 3, FieldValue< 3 >::Scalar > dg_penalty
Penalty enforcing inter-element continuity of solution (for each substance).
Discontinuous Galerkin method for equation of transport with dispersion.
unsigned int ndofs_
Number of dofs.
vector< LongIdx > dof_indices_
Vector of global DOF indices.
~InitProjectionAssemblyDG()
Destructor.
vector< PetscScalar > local_source_balance_rhs_
Auxiliary vector for set_sources method.
unsigned int max_edg_sides
Maximal number of edge sides (evaluate from dim 1,2,3)
Range< EdgePoint > edge_points(const DHCellSide &cell_side) const
Return EdgePoint range of appropriate dimension.
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.
~SourcesAssemblyDG()
Destructor.
Quadrature * quad_low_
Quadrature used in assembling methods (dim-1).
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
arma::vec::fixed< spacedim > normal_vector(unsigned int point_no)
Returns the normal vector to a side at given quadrature point.
EqFields * eq_fields_
Data objects shared with TransportDG.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
double JxW(const unsigned int point_no)
Return the product of Jacobian determinant and the quadrature weight at given quadrature point.
unsigned int ndofs_
Number of dofs.
FieldSet used_fields_
Sub field set contains fields used in calculation.
Directing class of FieldValueCache.
unsigned int get_dof_indices(std::vector< LongIdx > &indices) const
Fill vector of the global indices of dofs associated to the cell.
Class FEValues calculates finite element data on the actual cells such as shape function values,...
vector< PetscScalar > local_retardation_balance_vector_
Auxiliary vector for assemble mass matrix.
bool is_own() const
Return true if accessor represents own element (false for ghost element)
InitConditionAssemblyDG(EqFields *eq_fields, EqData *eq_data)
Constructor.
@ update_values
Shape function values.
unsigned int size() const
Returns number of quadrature points.
void initialize(Quadrature &_quadrature, FiniteElement< DIM > &_fe, UpdateFlags _flags)
Initialize structures and calculates cell-independent data.
static constexpr const char * name()
double measure() const
Calculate metrics of the side.
LinSys ** ls_dt
Linear algebra system for the time derivative (actually it is used only for handling the matrix struc...
@ update_quadrature_points
Transformed quadrature points.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
#define ASSERT_EQ_DBG(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
InitProjectionAssemblyDG(EqFields *eq_fields, EqData *eq_data)
Constructor.
static constexpr const char * name()
unsigned int dim() const
Return dimension of element appropriate to cell.
@ update_normal_vectors
Normal vectors.
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
unsigned int dim() const
Return dimension of element appropriate to the side.
static constexpr double almost_zero
void begin() override
Implements AssemblyBase::begin.
void set_DG_parameters_boundary(Side side, const int K_size, const std::vector< arma::mat33 > &K, const double flux, const arma::vec3 &normal_vector, const double alpha, double &gamma)
Sets up parameters of the DG method on a given boundary edge.
vector< LongIdx > dof_indices_
Vector of global DOF indices.
vector< FEValues< 3 > * > fv_sb_
Auxiliary vector, holds FEValues objects for assemble element-side.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
unsigned int ndofs_
Number of dofs.
PetscScalar local_flux_balance_rhs_
Auxiliary variable for set_boundary_conditions method.
std::vector< Vec > ret_vec
Auxiliary vectors for calculation of sources in balance due to retardation (e.g. sorption).
vector< PetscScalar > local_matrix_
Auxiliary vector for assemble methods.
StiffnessAssemblyDG(EqFields *eq_fields, EqData *eq_data)
Constructor.
Side accessor allows to iterate over sides of DOF handler cell.
int dg_variant
DG variant ((non-)symmetric/incomplete.
TransportDG< Model >::EqFields EqFields
LinSys ** ls
Linear algebra system for the transport equation.
void end() override
Implements AssemblyBase::end.
TransportDG< Model >::EqFields EqFields
static constexpr const char * name()
unsigned int cond_idx() const
Returns global index of the prescribed boundary condition.
std::shared_ptr< Balance > balance_
const DHCellAccessor & cell() const
Return DHCellAccessor appropriate to the side.
unsigned int ndofs_
Number of dofs.
FEValues< 3 > fe_values_vb_
FEValues of dim-1 object (of P disc finite element type)
void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side)
Assembles the fluxes between elements of different dimensions.
Definitions of particular quadrature rules on simplices.
EqFields * eq_fields_
Data objects shared with TransportDG.
Armor::Array< double >::ArrayMatSet set(uint i)
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
std::shared_ptr< DOFHandlerMultiDim > dh_
Object for distribution of dofs.
static constexpr const char * name()
TransportDG< Model >::EqData EqData
const TimeStep & step(int index=-1) const
unsigned int ndofs_
Number of dofs.
double weight(unsigned int i) const
Returns the ith weight.
void boundary_side_integral(DHCellSide cell_side)
Implements AssemblyBase::boundary_side_integral.
FEValues< 3 > fe_values_
FEValues of object (of P disc finite element type)
MassAssemblyDG(EqFields *eq_fields, EqData *eq_data)
Constructor.
FieldSet used_fields_
Sub field set contains fields used in calculation.
const ElementAccessor< 3 > elm() const
Return ElementAccessor to element of loc_ele_idx_.
void begin() override
Implements AssemblyBase::begin.
vector< vector< arma::mat33 > > dif_coef
Diffusion coefficients.
Container for various descendants of FieldCommonBase.
vector< PetscScalar > local_rhs_
Auxiliary vector for set_sources method.
ElementAccessor< 3 > element_accessor()
vector< LongIdx > dof_indices_
Vector of global DOF indices.
void end() override
Implements AssemblyBase::end.
virtual void rhs_set_values(int nrow, int *rows, double *vals)=0
vector< PetscScalar > local_matrix_
Auxiliary vector for assemble methods.
Abstract class for the description of a general finite element on a reference simplex in dim dimensio...
TransportDG< Model >::EqData EqData
MultiField< 3, FieldValue< 3 >::Scalar > fracture_sigma
Transition parameter for diffusive transfer on fractures (for each substance).
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
ElementAccessor< 3 > element() const
TransportDG< Model >::EqData EqData
Range< BulkPoint > bulk_points(unsigned int element_patch_idx) const
Return BulkPoint range of appropriate dimension.
EqFields * eq_fields_
Data objects shared with TransportDG.
unsigned int qsize_lower_dim_
Size of quadrature of dim-1.
FEValues< 3 > fe_values_
FEValues of object (of P disc finite element type)
static LocalPoint node_coords(unsigned int nid)
void boundary_side_integral(DHCellSide cell_side)
Assembles the fluxes on the boundary.
vector< PetscScalar > local_rhs_
Auxiliary vector for set_sources method.
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
Cell accessor allow iterate over DOF handler cells.
vector< double * > waverages
Auxiliary storage for weighted averages of shape functions.
vector< LongIdx > dof_indices_
Vector of global DOF indices.
void edge_integral(RangeConvert< DHEdgeSide, DHCellSide > edge_side_range)
Assembles the fluxes between sides of elements of the same dimension.
vector< LongIdx > dof_indices_
Vector of global DOF indices.
~MassAssemblyDG()
Destructor.
~BdrConditionAssemblyDG()
Destructor.
vector< PetscScalar > local_mass_balance_vector_
Same as previous.
vector< PetscScalar > local_source_balance_vector_
Auxiliary vector for set_sources method.
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
int active_integrals_
Holds mask of active integrals.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
TransportDG< Model >::EqFields EqFields
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.
FieldSet used_fields_
Sub field set contains fields used in calculation.
@ update_JxW_values
Transformed quadrature weights.
std::vector< std::vector< double > > gamma
Penalty parameters.
LocDofVec get_loc_dof_indices() const
Returns the local indices of dofs associated to the cell on the local process.
FEValues< 3 > fe_values_
FEValues of object (of P disc finite element type)
vector< vector< LongIdx > > side_dof_indices_
Vector of vectors of side DOF indices.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
unsigned int idx() const
Return local idx of element in boundary / bulk part of element vector.
UpdateFlags
Enum type UpdateFlags indicates which quantities are to be recomputed on each finite element cell.
vector< double * > averages
Auxiliary storage for averages of shape functions.
FEValues< 3 > fe_values_
FEValues of object (of P disc finite element type)
TransportDG< Model >::EqData EqData
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
shared_ptr< FiniteElement< dim-1 > > fe_low_
Finite element for the solution of the advection-diffusion equation (dim-1).
void reinit(const ElementAccessor< spacedim > &cell)
Update cell-dependent data (gradients, Jacobians etc.)
TransportDG< Model >::EqFields EqFields
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assembles the cell (volume) integral into the stiffness matrix.
FEValues< 3 > fe_values_side_
FEValues of object (of P disc finite element type)
@ update_side_JxW_values
Transformed quadrature weight for cell sides.
void end() override
Implements AssemblyBase::end.
std::vector< VectorMPI > output_vec
Vector of solution data.
vector< double * > jumps
Auxiliary storage for jumps of shape functions.
unsigned int index() const
static constexpr const char * name()
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.
FieldSet used_fields_
Sub field set contains fields used in calculation.
Range< CouplingPoint > coupling_points(const DHCellSide &cell_side) const
Return CouplingPoint range of appropriate dimension.
#define DebugOut()
Macro defining 'debug' record of log.
BdrConditionAssemblyDG(EqFields *eq_fields, EqData *eq_data)
Constructor.
IterConvert< ObjectIn, ObjectOut > begin()
Iterator to begin item of range.
Base class for quadrature rules on simplices in arbitrary dimensions.
vector< PetscScalar > local_rhs_
Auxiliary vector for set_sources method.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
EqFields * eq_fields_
Data objects shared with TransportDG.
FEValues< 3 > fe_values_side_
FEValues of object (of P disc finite element type)
FieldSet used_fields_
Sub field set contains fields used in calculation.
EqFields * eq_fields_
Data objects shared with TransportDG.
unsigned int bulk_idx() const
Returns index of the region in the bulk set.
vector< FEValues< 3 > > fe_values_vec_
Vector of FEValues of object (of P disc finite element types)
double elem_anisotropy(ElementAccessor< 3 > e) const
Compute and return anisotropy of given element.
Generic class of assemblation.
static constexpr const char * name()
TransportDG< Model >::EqData EqData
Range< BoundaryPoint > boundary_points(const DHCellSide &cell_side) const
Return BoundaryPoint range of appropriate dimension.
Discontinuous Lagrangean finite element on dim dimensional simplex.
ElementCacheMap * element_cache_map_
ElementCacheMap shared with GenericAssembly object.
vector< LongIdx > side_dof_indices_vb_
Vector of side DOF indices (assemble element-side fluxex)
TransportDG< Model >::EqData EqData
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
Quadrature * quad_
Quadrature used in assembling methods.
virtual void mat_set_values(int nrow, int *rows, int ncol, int *cols, double *vals)=0
unsigned int IntDim
A dimension index type.