21 #ifndef PATCH_FE_VALUES_HH_
22 #define PATCH_FE_VALUES_HH_
51 template <
class ValueType>
72 template <
class ValueType>
97 template <
class ValueType>
116 unsigned int begin,
unsigned int begin_side,
unsigned int n_dofs_bulk,
unsigned int n_dofs_side,
unsigned int join_idx)
176 template<
unsigned int dim>
185 template<
unsigned int FE_dim>
188 if (fe_sys !=
nullptr) {
189 return fe_sys->
fe()[component_idx];
191 ASSERT_EQ(component_idx, 0).warning(
"Non-zero component_idx can only be used for FESystem.");
203 template<
unsigned int FE_dim>
207 arma::mat shape_values(fe->n_dofs(), fe->n_components());
208 for (
unsigned int i=0; i<q->
size(); i++)
210 for (
unsigned int j=0; j<fe->n_dofs(); j++)
212 for (
unsigned int c=0; c<fe->n_components(); c++)
213 shape_values(j,c) = fe->shape_value(j, q->
point<FE_dim>(i), c);
215 ref_shape_vals[i][j] = trans(shape_values.row(j));
219 return ref_shape_vals;
230 template<
unsigned int FE_dim>
234 arma::mat shape_values(fe->n_dofs(), fe->n_components());
236 for (
unsigned int sid=0; sid<FE_dim+1; sid++) {
238 for (
unsigned int i=0; i<quad.size(); i++)
240 for (
unsigned int j=0; j<fe->n_dofs(); j++)
242 for (
unsigned int c=0; c<fe->n_components(); c++) {
243 shape_values(j,c) = fe->shape_value(j, quad.template point<FE_dim>(i), c);
246 ref_shape_vals[sid][i][j] = trans(shape_values.row(j));
251 return ref_shape_vals;
256 template<
unsigned int dim>
304 auto fe_component = this->
fe_comp(
fe_, component_idx);
305 ASSERT_EQ(fe_component->fe_type(),
FEType::FEScalar).error(
"Type of FiniteElement of scalar_shape accessor must be FEScalar!\n");
311 for (
unsigned int j = 0; j < fe_component->n_dofs(); j++) {
312 shape_values[i][j] = ref_shape_vals[i][j][0];
319 uint begin = scalar_shape_bulk_op.result_row();
338 auto fe_component = this->
fe_comp(
fe_, component_idx);
339 ASSERT_EQ(fe_component->fe_type(),
FEType::FEScalar).error(
"Type of FiniteElement of grad_scalar_shape accessor must be FEScalar!\n");
345 bulk_reinit::ptop_scalar_shape_grads<dim>(operations, op_results, ref_shape_grads, scalar_shape_grads_op_idx);
348 uint begin = grad_scalar_shape_bulk_op.result_row();
369 for (
unsigned int i_pt=0; i_pt<q->
size(); i_pt++)
371 for (
unsigned int i_dof=0; i_dof<fe->n_dofs(); i_dof++)
374 for (
unsigned int c=0; c<fe->n_components(); c++)
375 grad.col(c) += fe->shape_grad(i_dof, q->
point<dim>(i_pt), c);
377 ref_shape_grads[i_pt][i_dof] = grad;
381 return ref_shape_grads;
385 std::shared_ptr< FiniteElement<dim> >
fe_;
389 template<
unsigned int dim>
435 auto fe_component = this->
fe_comp(
fe_, component_idx);
436 ASSERT_EQ(fe_component->fe_type(),
FEType::FEScalar).error(
"Type of FiniteElement of scalar_shape accessor must be FEScalar!\n");
444 for (
unsigned int s=0; s<dim+1; ++s)
446 for (
unsigned int j = 0; j < fe_component->n_dofs(); j++) {
447 shape_values[s][i][j] = ref_shape_vals[s][i][j][0];
454 uint begin = scalar_shape_bulk_op.result_row();
462 auto fe_component = this->
fe_comp(
fe_, component_idx);
463 ASSERT_EQ(fe_component->fe_type(),
FEType::FEScalar).error(
"Type of FiniteElement of grad_scalar_shape accessor must be FEScalar!\n");
469 side_reinit::ptop_scalar_shape_grads<dim>(operations, op_results, el_table, ref_shape_grads, scalar_shape_grads_op_idx);
472 uint begin = grad_scalar_shape_side_op.result_row();
491 for (
unsigned int sid=0; sid<dim+1; sid++) {
493 for (
unsigned int i_pt=0; i_pt<quad.size(); i_pt++)
495 for (
unsigned int i_dof=0; i_dof<fe->n_dofs(); i_dof++)
498 for (
unsigned int c=0; c<fe->n_components(); c++)
499 grad.col(c) += fe->shape_grad(i_dof, quad.template point<dim>(i_pt), c);
501 ref_shape_grads[sid][i_pt][i_dof] = grad;
506 return ref_shape_grads;
510 std::shared_ptr< FiniteElement<dim> >
fe_;
514 template<
unsigned int dim>
531 ASSERT_EQ(fe_component_low->fe_type(),
FEType::FEScalar).error(
"Type of FiniteElement of scalar_shape accessor must be FEScalar!\n");
536 for (
unsigned int j = 0; j < fe_component_low->n_dofs(); j++) {
537 shape_values_bulk[i][j] = ref_shape_vals_bulk[i][j][0];
544 uint begin_bulk = grad_scalar_shape_bulk_op.result_row();
548 ASSERT_EQ(fe_component_high->fe_type(),
FEType::FEScalar).error(
"Type of FiniteElement of scalar_shape accessor must be FEScalar!\n");
555 for (
unsigned int s=0; s<dim+1; ++s)
557 for (
unsigned int j = 0; j < fe_component_high->n_dofs(); j++) {
558 shape_values_side[s][i][j] = ref_shape_vals_side[s][i][j][0];
565 uint begin_side = grad_scalar_shape_side_op.result_row();
568 begin_bulk, begin_side, fe_component_low->n_dofs(), fe_component_high->n_dofs(), 0) );
569 unsigned int end_idx = fe_component_low->n_dofs() + fe_component_high->n_dofs();
571 begin_bulk, begin_side, fe_component_low->n_dofs(), fe_component_high->n_dofs(), end_idx) );
600 template<
unsigned int spacedim = 3>
647 used_quads_[0] =
false; used_quads_[1] =
false;
651 : asm_arena_(1024 * 1024, 256),
652 patch_arena_(1024 * 1024, 256),
661 used_quads_[0] =
false; used_quads_[1] =
false;
671 if (is_bulk)
return patch_point_vals_bulk_[dim-1].get_quadrature();
672 else return patch_point_vals_side_[dim-1].get_quadrature();
682 template<
unsigned int DIM>
685 if ( _quadrature.
dim() == DIM ) {
686 used_quads_[0] =
true;
687 patch_point_vals_bulk_[DIM-1].initialize(3);
689 used_quads_[1] =
true;
690 patch_point_vals_side_[DIM-1].initialize(4);
697 for (
unsigned int i=0; i<3; ++i) {
698 if (used_quads_[0]) patch_point_vals_bulk_[i].reset();
699 if (used_quads_[1]) patch_point_vals_side_[i].reset();
706 for (
unsigned int i=0; i<3; ++i) {
707 if (used_quads_[0]) patch_point_vals_bulk_[i].reinit_patch();
708 if (used_quads_[1]) patch_point_vals_side_[i].reinit_patch();
715 template<
unsigned int dim>
717 ASSERT((dim>=0) && (dim<=3))(dim).error(
"Dimension must be 0, 1, 2 or 3.");
722 template<
unsigned int dim>
724 ASSERT((dim>0) && (dim<=3))(dim).error(
"Dimension must be 1, 2 or 3.");
729 template<
unsigned int dim>
731 ASSERT((dim>0) && (dim<=3))(dim).error(
"Dimension must be 1, 2 or 3.");
736 template<
unsigned int dim>
739 return JoinValues<dim>(&patch_point_vals_bulk_[dim-2], &patch_point_vals_side_[dim-1], fe_);
746 for (
uint i=0; i<3; ++i) {
747 if (used_quads_[0]) patch_point_vals_bulk_[i].resize_tables(table_sizes.
elem_sizes_[0][i], table_sizes.
point_sizes_[0][i]);
748 if (used_quads_[1]) patch_point_vals_side_[i].resize_tables(table_sizes.
elem_sizes_[1][i], table_sizes.
point_sizes_[1][i]);
755 switch (cell.
dim()) {
758 return patch_point_vals_bulk_[0].register_element(coords, element_patch_idx);
762 return patch_point_vals_bulk_[1].register_element(coords, element_patch_idx);
766 return patch_point_vals_bulk_[2].register_element(coords, element_patch_idx);
778 for (
unsigned int n=0; n<cell_side.
dim(); n++)
779 for (
unsigned int c=0; c<spacedim; c++)
780 side_coords(c,n) = (*cell_side.
side().
node(n))[c];
784 switch (cell.
dim()) {
787 return patch_point_vals_side_[0].register_side(elm_coords, side_coords);
791 return patch_point_vals_side_[1].register_side(elm_coords, side_coords);
795 return patch_point_vals_side_[2].register_side(elm_coords, side_coords);
806 return patch_point_vals_bulk_[cell.
dim()-1].register_bulk_point(elem_table_row, value_patch_idx, cell.
elm_idx());
811 return patch_point_vals_side_[cell_side.
dim()-1].register_side_point(elem_table_row, value_patch_idx, cell_side.
elem_idx(), cell_side.
side_idx());
816 stream << endl <<
"Table of patch FE data:" << endl;
817 for (
uint i=0; i<3; ++i) {
818 stream << std::setfill(
'-') << setw(100) <<
"" << endl;
819 stream <<
"Bulk, dimension " << (i+1) << endl;
820 patch_point_vals_bulk_[i].print_data_tables(stream, points, ints);
823 for (
uint i=0; i<3; ++i) {
824 stream << std::setfill(
'-') << setw(100) <<
"" << endl;
825 stream <<
"Side, dimension " << (i+1) << endl;
826 patch_point_vals_side_[i].print_data_tables(stream, points, ints);
828 stream << std::setfill(
'=') << setw(100) <<
"" << endl;
833 stream << endl <<
"Table of patch FE operations:" << endl;
834 for (
uint i=0; i<3; ++i) {
835 stream << std::setfill(
'-') << setw(100) <<
"" << endl;
836 stream <<
"Bulk, dimension " << (i+1) <<
", n_rows " << patch_point_vals_bulk_[i].n_rows() << endl;
837 patch_point_vals_bulk_[i].print_operations(stream, 0);
839 for (
uint i=0; i<3; ++i) {
840 stream << std::setfill(
'-') << setw(100) <<
"" << endl;
841 stream <<
"Side, dimension " << (i+1) <<
", n_rows " << patch_point_vals_side_[i].n_rows() << endl;
842 patch_point_vals_side_[i].print_operations(stream, 1);
844 stream << std::setfill(
'=') << setw(100) <<
"" << endl;
856 template <
class ValueType>
858 template <
class ValueType>
860 template <
class ValueType>
865 template <
class ValueType>
868 return patch_point_vals_.scalar_value(op_idx_, value_cache_idx);
874 return patch_point_vals_.vector_value(op_idx_, value_cache_idx);
880 return patch_point_vals_.tensor_value(op_idx_, value_cache_idx);
883 template <
class ValueType>
886 return patch_point_vals_.scalar_val(
begin_, value_cache_idx);
892 return patch_point_vals_.vector_val(
begin_, value_cache_idx);
898 return patch_point_vals_.tensor_val(
begin_, value_cache_idx);
901 template <
class ValueType>
904 return patch_point_vals_.scalar_val(
begin_+shape_idx, value_cache_idx);
910 return patch_point_vals_.vector_val(
begin_+3*shape_idx, value_cache_idx);
915 Tensor tens; tens.zeros();
919 template <
class ValueType>
922 return patch_point_vals_.scalar_val(
begin_+shape_idx, value_cache_idx);
928 return patch_point_vals_.vector_val(
begin_+3*shape_idx, value_cache_idx);
933 Tensor tens; tens.zeros();
938 template <
class ValueType>
950 Vector vect; vect.zeros();
956 Tensor tens; tens.zeros();
960 template <
class ValueType>
972 Vector vect; vect.zeros();
978 Tensor tens; tens.zeros();
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
std::vector< std::vector< std::vector< arma::vec > > > ref_shape_values_side(Quadrature *q, std::shared_ptr< FiniteElement< FE_dim >> fe)
Precomputed values of basis functions at the side quadrature points.
std::shared_ptr< FiniteElement< FE_dim > > fe_comp(std::shared_ptr< FiniteElement< FE_dim > > fe, uint component_idx)
Return FiniteElement of component_idx for FESystem or fe for other types.
std::vector< std::vector< arma::vec > > ref_shape_values_bulk(Quadrature *q, std::shared_ptr< FiniteElement< FE_dim >> fe)
Precomputed values of basis functions at the bulk quadrature points.
Base point accessor class.
const ElementCacheMap * elm_cache_map() const
unsigned int elem_patch_idx() const
unsigned int eval_point_idx() const
Return index in EvalPoints object.
ElQ< Scalar > JxW()
Register the product of Jacobian determinant and the quadrature weight at bulk quadrature points.
PatchPointValues< 3 > & patch_point_vals_
ElQ< Vector > coords()
Create bulk accessor of coords entity.
std::vector< std::vector< arma::mat > > ref_shape_gradients(std::shared_ptr< FiniteElement< dim >> fe)
Precomputed gradients of basis functions at the quadrature points.
ElQ< Scalar > determinant()
Create bulk accessor of jac determinant entity.
std::shared_ptr< FiniteElement< dim > > fe_
FeQ< Vector > grad_scalar_shape(uint component_idx=0)
Return the value of the function_no-th gradient shape function at the p bulk quadrature point.
FeQ< Scalar > scalar_shape(uint component_idx=0)
Return the value of the function_no-th shape function at the p bulk quadrature point.
BulkValues(PatchPointValues< 3 > &patch_point_vals, MixedPtr< FiniteElement > fe)
Constructor.
Cell accessor allow iterate over DOF handler cells.
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_.
ElementAccessor< 3 > elm() const
Return ElementAccessor to element of loc_ele_idx_.
Side accessor allows to iterate over sides of DOF handler cell.
unsigned int elem_idx() const
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.
unsigned int side_idx() const
ElQ(PatchPointValues< 3 > &patch_point_vals, unsigned int begin, unsigned int op_idx)
Constructor.
unsigned int op_idx_
Index of the first component of the bulk Quantity. Size is given by ValueType.
ElQ()=delete
Forbidden default constructor.
PatchPointValues< 3 > & patch_point_vals_
Reference to PatchPointValues.
ValueType operator()(FMT_UNUSED const SidePoint &point)
ValueType operator()(FMT_UNUSED const BulkPoint &point)
int element_eval_point(unsigned int i_elem_in_cache, unsigned int i_eval_point) const
Compound finite element on dim dimensional simplex.
const std::vector< std::shared_ptr< FiniteElement< dim > > > & fe() const
Bulk data specialization, order of item in operations_ vector corresponds to the BulkOps enum.
unsigned int begin_
Index of the first component of the Quantity. Size is given by ValueType.
PatchPointValues< 3 > & patch_point_vals_
Reference to PatchPointValues.
unsigned int n_dofs_
Number of DOFs.
ValueType operator()(FMT_UNUSED unsigned int shape_idx, FMT_UNUSED const BulkPoint &point)
FeQ(PatchPointValues< 3 > &patch_point_vals, unsigned int begin, unsigned int n_dofs)
FeQ()=delete
Forbidden default constructor.
ValueType operator()(FMT_UNUSED unsigned int shape_idx, FMT_UNUSED const SidePoint &point)
Bulk Side specialization, order of item in operations_ vector corresponds to the SideOps enum.
Abstract class for the description of a general finite element on a reference simplex in dim dimensio...
unsigned int n_dofs_both() const
unsigned int begin_side_
Index of the first component of the side Quantity. Size is given by ValueType.
unsigned int local_idx() const
Return local index of DOF (on low / high-dim) - should be private method.
bool operator==(const JoinShapeAccessor< ValueType > &other) const
Comparison of accessors.
unsigned int begin_
Index of the first component of the bulk Quantity. Size is given by ValueType.
PatchPointValues< 3 > * patch_point_vals_bulk_
Pointer to bulk PatchPointValues.
JoinShapeAccessor(PatchPointValues< 3 > *patch_point_vals_bulk, PatchPointValues< 3 > *patch_point_vals_side, unsigned int begin, unsigned int begin_side, unsigned int n_dofs_bulk, unsigned int n_dofs_side, unsigned int join_idx)
unsigned int n_dofs_high() const
JoinShapeAccessor()
Default constructor.
void inc()
Iterates to next item.
unsigned int join_idx_
Index of processed DOF.
unsigned int n_dofs_high_
Number of DOFs on high-dim element.
unsigned int join_idx() const
Return global index of DOF.
unsigned int n_dofs_low() const
ValueType operator()(const BulkPoint &point)
unsigned int n_dofs_low_
Number of DOFs on low-dim element.
PatchPointValues< 3 > * patch_point_vals_side_
Pointer to side PatchPointValues.
JoinValues(FMT_UNUSED PatchPointValues< 3 > *patch_point_vals_bulk, FMT_UNUSED PatchPointValues< 3 > *patch_point_vals_side, FMT_UNUSED MixedPtr< FiniteElement > fe)
Constructor.
Range< JoinShapeAccessor< Scalar > > scalar_join_shape(FMT_UNUSED uint component_idx=0)
std::shared_ptr< FiniteElement< dim > > fe_high_dim_
JoinValues(PatchPointValues< 3 > *patch_point_vals_bulk, PatchPointValues< 3 > *patch_point_vals_side, MixedPtr< FiniteElement > fe)
Constructor.
PatchPointValues< 3 > * patch_point_vals_bulk_
Range< JoinShapeAccessor< Scalar > > scalar_join_shape(uint component_idx=0)
std::shared_ptr< FiniteElement< dim-1 > > fe_low_dim_
PatchPointValues< 3 > * patch_point_vals_side_
static ElementMap element_map(ElementAccessor< 3 > elm)
uint register_element(DHCellAccessor cell, uint element_patch_idx)
Register element to patch_point_vals_ table by dimension of element.
~PatchFEValues()
Destructor.
std::array< FeBulk::PatchPointValues< spacedim >, 3 > patch_point_vals_bulk_
Sub objects of bulk data of dimensions 1,2,3.
uint register_bulk_point(DHCellAccessor cell, uint elem_table_row, uint value_patch_idx)
Register bulk point to patch_point_vals_ table by dimension of element.
uint register_side(DHCellSide cell_side)
Register side to patch_point_vals_ table by dimension of side.
void print_data_tables(ostream &stream, bool points, bool ints, bool only_bulk=true) const
Temporary development method.
JoinValues< dim > join_values()
Return JoinValue object of dimension given by template parameter.
Quadrature * get_quadrature(uint dim, bool is_bulk) const
Return bulk or side quadrature of given dimension.
BulkValues< dim > bulk_values()
Return BulkValue object of dimension given by template parameter.
void initialize(Quadrature &_quadrature)
Initialize structures and calculates cell-independent data.
void print_operations(ostream &stream) const
Temporary development method.
std::array< FeSide::PatchPointValues< spacedim >, 3 > patch_point_vals_side_
Sub objects of side data of dimensions 1,2,3.
MixedPtr< FiniteElement > fe_
Mixed of shared pointers of FiniteElement object.
SideValues< dim > side_values()
Return SideValue object of dimension given by template parameter.
AssemblyArena patch_arena_
void resize_tables(TableSizes table_sizes)
Resize tables of patch_point_vals_.
PatchFEValues(unsigned int quad_order, MixedPtr< FiniteElement > fe)
void reinit_patch()
Reinit data.
unsigned int n_dofs() const
Returns the number of shape functions.
uint register_side_point(DHCellSide cell_side, uint elem_table_row, uint value_patch_idx)
Register side point to patch_point_vals_ table by dimension of side.
void reset()
Reset PatchpointValues structures.
ElOp< spacedim > & make_fe_op(std::initializer_list< uint > shape, ReinitFunction reinit_f, std::vector< uint > input_ops_vec, uint n_dofs, OpSizeType size_type=pointOp)
std::vector< ElOp< spacedim > > operations_
Vector of all defined operations.
Quadrature * get_quadrature() const
Getter for quadrature.
uint dim() const
Getter for dim_.
Scalar scalar_val(uint result_row, uint point_idx) const
Base class for quadrature rules on simplices in arbitrary dimensions.
Quadrature make_from_side(unsigned int sid) const
unsigned int size() const
Returns number of quadrature points.
Armor::ArmaVec< double, point_dim > point(unsigned int i) const
Returns the ith quadrature point.
General point a+ side_begin_ + ccessor allow iterate over quadrature points of given side defined in ...
unsigned int eval_point_idx() const
Return index in EvalPoints object.
SideValues(PatchPointValues< 3 > &patch_point_vals, MixedPtr< FiniteElement > fe)
Constructor.
ElQ< Vector > normal_vector()
Register the normal vector to a side at side quadrature points.
std::shared_ptr< FiniteElement< dim > > fe_
ElQ< Scalar > JxW()
Same as BulkValues::JxW but register at side quadrature points.
ElQ< Vector > coords()
Create side accessor of coords entity.
FeQ< Scalar > scalar_shape(uint component_idx=0)
Same as BulkValues::scalar_shape but register at side quadrature points.
FeQ< Vector > grad_scalar_shape(uint component_idx=0)
Same as BulkValues::grad_scalar_shape but register at side quadrature points.
std::vector< std::vector< std::vector< arma::mat > > > ref_shape_gradients(std::shared_ptr< FiniteElement< dim >> fe)
Precomputed gradients of basis functions at the quadrature points.
PatchPointValues< 3 > & patch_point_vals_
ElQ< Scalar > determinant()
Create bulk accessor of jac determinant entity.
NodeAccessor< 3 > node(unsigned int i) const
Returns node for given local index i on the side.
Class ElementValues calculates data related to transformation of reference cell to actual cell (Jacob...
Class FESystem for compound finite elements.
Class FEValues calculates finite element data on the actual cells such as shape function values,...
Iter< Object > make_iter(Object obj)
Class MappingP1 implements the affine transformation of the unit cell onto the actual cell.
ArmaMat< double, N, M > mat
@ opCoords
operations evaluated on quadrature points
@ opInvJac
inverse Jacobian
@ opJxW
JxW value of quadrature point.
@ opNormalVec
normal vector of quadrature point
std::vector< std::array< uint, 3 > > DimPointTable
Holds triplet (dim; bulk/side; idx of point in subtable)
Store finite element data on the actual patch such as shape function values, gradients,...
Definitions of particular quadrature rules on simplices.
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.
Struct for pre-computing number of elements, sides, bulk points and side points on each dimension.
std::vector< std::vector< uint > > point_sizes_
void reset()
Set all values to zero.
std::vector< std::vector< uint > > elem_sizes_
static void ptop_scalar_shape(std::vector< ElOp< 3 >> &operations, TableDbl &op_results, std::vector< std::vector< double > > shape_values, uint scalar_shape_op_idx)
static void ptop_scalar_shape(std::vector< ElOp< 3 >> &operations, TableDbl &op_results, TableInt &el_table, std::vector< std::vector< std::vector< double > > > shape_values, uint scalar_shape_op_idx)
Enum type UpdateFlags indicates which quantities are to be recomputed on each finite element cell.