18 #ifndef ASSEMBLY_ELASTICITY_HH_
19 #define ASSEMBLY_ELASTICITY_HH_
34 template <
unsigned int dim>
41 static constexpr
const char *
name() {
return "StiffnessAssemblyElasticity"; }
63 shared_ptr<FE_P<dim>> fe_p = std::make_shared< FE_P<dim> >(1);
64 shared_ptr<
FE_P<dim-1>> fe_p_low = std::make_shared<
FE_P<dim-1> >(1);
65 fe_ = std::make_shared<FESystem<dim>>(fe_p,
FEVector, 3);
83 for (
uint m=0; m<2; ++m) {
85 for (
uint n=0; n<2; ++n)
97 if (cell.
dim() != dim)
return;
105 for (
unsigned int i=0; i<
n_dofs_; i++)
106 for (
unsigned int j=0; j<
n_dofs_; j++)
110 for (
auto p : this->
bulk_points(element_patch_idx) )
112 for (
unsigned int i=0; i<
n_dofs_; i++)
114 for (
unsigned int j=0; j<
n_dofs_; j++)
128 ASSERT_EQ(cell_side.
dim(), dim).error(
"Dimension of element mismatch!");
136 for (
unsigned int i=0; i<
n_dofs_; i++)
137 for (
unsigned int j=0; j<
n_dofs_; j++)
143 double side_measure = cell_side.
measure();
148 for (
unsigned int i=0; i<
n_dofs_; i++)
149 for (
unsigned int j=0; j<
n_dofs_; j++)
159 for (
unsigned int i=0; i<
n_dofs_; i++)
160 for (
unsigned int j=0; j<
n_dofs_; j++)
174 if (dim == 1)
return;
175 ASSERT_EQ(cell_lower_dim.
dim(), dim-1).error(
"Dimension of element mismatch!");
186 bool own_element_id[2];
187 own_element_id[0] = cell_lower_dim.
is_own();
188 own_element_id[1] = cell_higher_dim.
is_own();
190 for (
unsigned int n=0; n<2; ++n)
191 for (
unsigned int i=0; i<
n_dofs_; i++)
192 for (
unsigned int m=0; m<2; ++m)
193 for (
unsigned int j=0; j<
n_dofs_; j++)
200 auto p_low = p_high.lower_dim(cell_lower_dim);
203 for (
int n=0; n<2; n++)
205 if (!own_element_id[n])
continue;
213 for (
int m=0; m<2; m++)
218 double divuft = (m==0) ? arma::trace(guft) : 0;
236 for (
unsigned int n=0; n<2; ++n)
237 for (
unsigned int m=0; m<2; ++m)
252 shared_ptr<FiniteElement<dim>>
fe_;
277 template <
template<
IntDim...>
class DimAssembly>
283 template <
unsigned int dim>
290 static constexpr
const char *
name() {
return "RhsAssemblyElasticity"; }
317 shared_ptr<FE_P<dim>> fe_p = std::make_shared< FE_P<dim> >(1);
318 shared_ptr<
FE_P<dim-1>> fe_p_low = std::make_shared<
FE_P<dim-1> >(1);
319 fe_ = std::make_shared<FESystem<dim>>(fe_p,
FEVector, 3);
349 if (cell.
dim() != dim)
return;
350 if (!cell.
is_own())
return;
364 for (
auto p : this->
bulk_points(element_patch_idx) )
366 for (
unsigned int i=0; i<
n_dofs_; i++)
390 ASSERT_EQ(cell_side.
dim(), dim).error(
"Dimension of element mismatch!");
411 for (
unsigned int i=0; i<
n_dofs_; i++)
422 double side_measure = cell_side.
measure();
426 for (
unsigned int i=0; i<
n_dofs_; i++)
435 double side_measure = cell_side.
measure();
439 for (
unsigned int i=0; i<
n_dofs_; i++)
452 for (
unsigned int i=0; i<
n_dofs_; i++)
464 for (
unsigned int i=0; i<
n_dofs_; i++)
484 if (dim == 1)
return;
485 ASSERT_EQ(cell_lower_dim.
dim(), dim-1).error(
"Dimension of element mismatch!");
496 bool own_element_id[2];
497 own_element_id[0] = cell_lower_dim.
is_own();
498 own_element_id[1] = cell_higher_dim.
is_own();
500 for (
unsigned int n=0; n<2; ++n)
501 for (
unsigned int i=0; i<
n_dofs_; i++)
508 auto p_low = p_high.lower_dim(cell_lower_dim);
511 for (
int n=0; n<2; n++)
513 if (!own_element_id[n])
continue;
527 for (
unsigned int n=0; n<2; ++n)
534 shared_ptr<FiniteElement<dim>>
fe_;
562 template <
template<
IntDim...>
class DimAssembly>
567 template <
unsigned int dim>
574 static constexpr
const char *
name() {
return "OutpuFieldsAssemblyElasticity"; }
594 shared_ptr<FE_P<dim>> fe_p = std::make_shared< FE_P<dim> >(1);
595 fe_ = std::make_shared<FESystem<dim>>(fe_p,
FEVector, 3);
617 if (cell.
dim() != dim)
return;
618 if (!cell.
is_own())
return;
622 auto elm = cell.
elm();
629 auto p = *( this->
bulk_points(element_patch_idx).begin() );
633 for (
unsigned int i=0; i<
n_dofs_; i++)
640 arma::mat33 stress_dev = stress - arma::trace(stress)/3*arma::eye(3,3);
641 double von_mises_stress = sqrt(1.5*arma::dot(stress_dev, stress_dev));
642 double mean_stress = arma::trace(stress) / 3;
645 for (
unsigned int i=0; i<3; i++)
646 for (
unsigned int j=0; j<3; j++)
657 if (dim == 1)
return;
658 ASSERT_EQ(cell_lower_dim.
dim(), dim-1).error(
"Dimension of element mismatch!");
670 auto p_low = p_high.lower_dim(cell_lower_dim);
672 for (
unsigned int i=0; i<
n_dofs_; i++)
681 for (
unsigned int i=0; i<3; i++)
682 for (
unsigned int j=0; j<3; j++)
691 shared_ptr<FiniteElement<dim>>
fe_;
721 template <
template<
IntDim...>
class DimAssembly>
731 template <
unsigned int dim>
738 static constexpr
const char *
name() {
return "ConstraintAssemblyElasticity"; }
755 shared_ptr<FE_P<dim>> fe_p = std::make_shared< FE_P<dim> >(1);
756 fe_ = std::make_shared<FESystem<dim>>(fe_p,
FEVector, 3);
769 if (dim == 1)
return;
770 if (!cell_lower_dim.
is_own())
return;
772 ASSERT_EQ(cell_lower_dim.
dim(), dim-1).error(
"Dimension of element mismatch!");
779 for (
unsigned int i=0; i<
n_dofs_; i++)
786 double local_vector = 0;
789 auto p_low = p_high.lower_dim(cell_lower_dim);
794 for (
unsigned int i=0; i<
n_dofs_; i++)
812 shared_ptr<FiniteElement<dim>>
fe_;
829 template <
template<
IntDim...>
class DimAssembly>
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Range< BulkPoint > bulk_points(unsigned int element_patch_idx) const
Return BulkPoint range of appropriate dimension.
Quadrature * quad_
Quadrature used in assembling methods.
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.
ElementAccessor< 3 > element_accessor()
~ConstraintAssemblyElasticity()
Destructor.
FEValues< 3 > fe_values_side_
FEValues of side object.
EqFields * eq_fields_
Data objects shared with Elasticity.
void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side)
Assembles between elements of different dimensions.
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
static constexpr const char * name()
ConstraintAssemblyElasticity(EqFields *eq_fields, EqData *eq_data)
Constructor.
vector< vector< LongIdx > > side_dof_indices_
2 items vector of DOF indices in neighbour calculation.
vector< PetscScalar > local_matrix_
Auxiliary vector for assemble methods.
unsigned int n_dofs_
Number of dofs.
vector< LongIdx > dof_indices_
Vector of global DOF indices.
const FEValuesViews::Vector< 3 > * vec_view_side_
Vector view in boundary / neighbour calculation.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
FieldSet used_fields_
Sub field set contains fields used in calculation.
Elasticity::EqFields EqFields
Elasticity::EqData EqData
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.
DHCellAccessor cell_with_other_dh(const DOFHandlerMultiDim *dh) const
Create new accessor with same local idx and given DOF handler. Actual and given DOF handler must be c...
unsigned int 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.
unsigned int elm_idx() const
Return serial idx to element of loc_ele_idx_.
Side accessor allows to iterate over sides of DOF handler cell.
Side side() const
Return Side of given cell and side_idx.
const DHCellAccessor & cell() const
Return DHCellAccessor appropriate to the side.
unsigned int dim() const
Return dimension of element appropriate to the side.
ElementAccessor< 3 > element() const
LinSys * ls
Linear algebraic system.
std::shared_ptr< DOFHandlerMultiDim > dh_scalar_
std::shared_ptr< DOFHandlerMultiDim > dh_tensor_
std::shared_ptr< DOFHandlerMultiDim > dh_
Objects for distribution of dofs.
std::map< LongIdx, LongIdx > constraint_idx
Field< 3, FieldValue< 3 >::Scalar > potential_load
Potential of an additional (external) load.
@ bc_type_displacement_normal
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::TensorFixed > > output_stress_ptr
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > output_von_mises_stress_ptr
Field< 3, FieldValue< 3 >::Scalar > cross_section_min
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > output_cross_section_ptr
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::VectorFixed > > output_field_ptr
BCField< 3, FieldValue< 3 >::Enum > bc_type
Field< 3, FieldValue< 3 >::Scalar > lame_lambda
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > output_div_ptr
BCField< 3, FieldValue< 3 >::VectorFixed > bc_displacement
Field< 3, FieldValue< 3 >::TensorFixed > initial_stress
BCField< 3, FieldValue< 3 >::VectorFixed > bc_traction
Field< 3, FieldValue< 3 >::Scalar > ref_potential_load
Potential of reference external load on boundary. TODO: Switch to BCField when possible.
Field< 3, FieldValue< 3 >::VectorFixed > load
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > output_mean_stress_ptr
Field< 3, FieldValue< 3 >::Scalar > fracture_sigma
Transition parameter for diffusive transfer on fractures.
Field< 3, FieldValue< 3 >::Scalar > dirichlet_penalty
Field< 3, FieldValue< 3 >::Scalar > cross_section
Pointer to DarcyFlow field cross_section.
BCField< 3, FieldValue< 3 >::TensorFixed > bc_stress
Field< 3, FieldValue< 3 >::Scalar > lame_mu
double measure() const
Computes the measure of the element.
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
Directing class of FieldValueCache.
unsigned int n_neighs_vb() const
Return number of neighbours.
Compound finite element on dim dimensional simplex.
arma::mat::fixed< spacedim, spacedim > grad(unsigned int function_no, unsigned int point_no) const
Return gradient of vector-valued shape function.
arma::mat::fixed< spacedim, spacedim > sym_grad(unsigned int function_no, unsigned int point_no) const
Return symmetric gradient of vector-valued shape function.
arma::vec::fixed< spacedim > value(unsigned int function_no, unsigned int point_no) const
Return value of vector-valued shape function.
double divergence(unsigned int function_no, unsigned int point_no) const
Return divergence of vector-valued shape function.
double JxW(const unsigned int point_no)
Return the product of Jacobian determinant and the quadrature weight at given quadrature point.
const FEValuesViews::Vector< spacedim > & vector_view(unsigned int i) const
Accessor to vector values of multicomponent FE.
void initialize(Quadrature &_quadrature, FiniteElement< DIM > &_fe, UpdateFlags _flags)
Initialize structures and calculates cell-independent data.
void reinit(const ElementAccessor< spacedim > &cell)
Update cell-dependent data (gradients, Jacobians etc.)
arma::vec::fixed< spacedim > normal_vector(unsigned int point_no)
Returns the normal vector to a side at given quadrature point.
Conforming Lagrangean finite element on dim dimensional simplex.
Container for various descendants of FieldCommonBase.
Abstract class for the description of a general finite element on a reference simplex in dim dimensio...
Generic class of assemblation.
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
Elasticity::EqData EqData
LocDofVec dof_indices_scalar_
Vector of local DOF indices of scalar fields.
VectorMPI output_vec_
Data vectors of output fields (FieldFE).
unsigned int n_dofs_
Number of dofs.
FieldSet used_fields_
Sub field set contains fields used in calculation.
FEValues< 3 > fv_
FEValues of cell object (FESystem of P disc finite element type)
static constexpr const char * name()
LocDofVec dof_indices_
Vector of local DOF indices of vector fields.
VectorMPI output_div_vec_
VectorMPI output_stress_vec_
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
VectorMPI output_von_mises_stress_vec_
LocDofVec dof_indices_tensor_
Vector of local DOF indices of tensor fields.
const FEValuesViews::Vector< 3 > * vec_view_side_
Vector view in neighbour calculation.
EqFields * eq_fields_
Data objects shared with Elasticity.
arma::mat33 normal_stress_
Holds constributions of normal stress.
FEValues< 3 > fsv_
FEValues of side (neighbour integral) object.
Elasticity::EqFields EqFields
VectorMPI output_cross_sec_vec_
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side)
Assembles between elements of different dimensions.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
VectorMPI output_mean_stress_vec_
double normal_displacement_
Holds constributions of normal displacement.
OutpuFieldsAssemblyElasticity(EqFields *eq_fields, EqData *eq_data)
Constructor.
const FEValuesViews::Vector< 3 > * vec_view_
Vector view in cell integral calculation.
~OutpuFieldsAssemblyElasticity()
Destructor.
FEValues< 3 > fe_values_sub_
FEValues of lower dimension cell object.
void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side)
Assembles between elements of different dimensions.
vector< LongIdx > dof_indices_
Vector of global DOF indices.
~RhsAssemblyElasticity()
Destructor.
shared_ptr< FiniteElement< dim-1 > > fe_low_
Finite element for the solution of the advection-diffusion equation (dim-1).
vector< vector< LongIdx > > side_dof_indices_
2 items vector of DOF indices in neighbour calculation.
const FEValuesViews::Vector< 3 > * vec_view_bdr_
Vector view in boundary calculation.
vector< vector< PetscScalar > > local_rhs_ngh_
Auxiliary vectors for assemble ngh integral.
FEValues< 3 > fe_values_bdr_side_
FEValues of side (boundary integral) object.
const FEValuesViews::Vector< 3 > * vec_view_sub_
Vector view of low dim element in neighbour calculation.
EqFields * eq_fields_
Data objects shared with Elasticity.
Elasticity::EqFields EqFields
const FEValuesViews::Vector< 3 > * vec_view_
Vector view in cell integral calculation.
FEValues< 3 > fe_values_side_
FEValues of side (neighbour integral) object.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
std::vector< unsigned int > n_dofs_ngh_
Number of dofs on lower and higher dimension element (vector of 2 items)
Elasticity::EqData EqData
FEValues< 3 > fe_values_
FEValues of cell object (FESystem of P disc finite element type)
RhsAssemblyElasticity(EqFields *eq_fields, EqData *eq_data)
Constructor.
unsigned int n_dofs_sub_
Number of dofs (on lower dim element)
vector< PetscScalar > local_rhs_
Auxiliary vector for assemble methods.
unsigned int n_dofs_
Number of dofs.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
static constexpr const char * name()
const FEValuesViews::Vector< 3 > * vec_view_side_
Vector view in neighbour calculation.
void boundary_side_integral(DHCellSide cell_side)
Assembles boundary integral.
FieldSet used_fields_
Sub field set contains fields used in calculation.
FEValues< 3 > fe_values_side_
FEValues of side object.
void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side)
Assembles between elements of different dimensions.
arma::mat33 mat_t(const arma::mat33 &m, const arma::vec3 &n)
FieldSet used_fields_
Sub field set contains fields used in calculation.
Elasticity::EqFields EqFields
const FEValuesViews::Vector< 3 > * vec_view_
Vector view in cell integral calculation.
Elasticity::EqData EqData
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
static constexpr const char * name()
const FEValuesViews::Vector< 3 > * vec_view_side_
Vector view in boundary / neighbour calculation.
shared_ptr< FiniteElement< dim-1 > > fe_low_
Finite element for the solution of the advection-diffusion equation (dim-1).
vector< PetscScalar > local_matrix_
Auxiliary vector for assemble methods.
~StiffnessAssemblyElasticity()
Destructor.
vector< vector< LongIdx > > side_dof_indices_
2 items vector of DOF indices in neighbour calculation.
vector< LongIdx > dof_indices_
Vector of global DOF indices.
EqFields * eq_fields_
Data objects shared with Elasticity.
FEValues< 3 > fe_values_sub_
FEValues of lower dimension cell object.
const FEValuesViews::Vector< 3 > * vec_view_sub_
Vector view of low dim element in neighbour calculation.
unsigned int n_dofs_sub_
Number of dofs (on lower dim element)
unsigned int n_dofs_
Number of dofs.
FEValues< 3 > fe_values_
FEValues of cell object (FESystem of P disc finite element type)
vector< vector< vector< PetscScalar > > > local_matrix_ngh_
Auxiliary vectors for assemble ngh integral.
std::vector< unsigned int > n_dofs_ngh_
Number of dofs on lower and higher dimension element (vector of 2 items)
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.
StiffnessAssemblyElasticity(EqFields *eq_fields, EqData *eq_data)
Constructor.
void boundary_side_integral(DHCellSide cell_side)
Assembles boundary integral.
double get(unsigned int pos) const
Return value on given position.
void set(unsigned int pos, double val)
Set value on given position.
void add(unsigned int pos, double val)
Add value to item on given position.
FEM for linear elasticity.
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,...
arma::Col< IntIdx > LocDofVec
unsigned int IntDim
A dimension index type.
Definitions of particular quadrature rules on simplices.
@ update_values
Shape function values.
@ 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.
@ update_quadrature_points
Transformed quadrature points.