17 #include <deal.II/base/utilities.h> 18 #include <deal.II/base/quadrature.h> 19 #include <deal.II/base/quadrature_lib.h> 20 #include <deal.II/base/qprojector.h> 21 #include <deal.II/base/table.h> 22 #include <deal.II/grid/tria.h> 23 #include <deal.II/grid/tria_iterator.h> 24 #include <deal.II/dofs/dof_accessor.h> 25 #include <deal.II/fe/fe.h> 26 #include <deal.II/fe/mapping.h> 27 #include <deal.II/fe/fe_abf.h> 28 #include <deal.II/fe/fe_values.h> 29 #include <deal.II/fe/fe_tools.h> 39 DEAL_II_NAMESPACE_OPEN
53 std::vector<bool>(dim,true))),
88 std::vector<FullMatrix<double> >
96 unsigned int target_row=0;
97 for (
unsigned int d=0; d<face_embeddings.size(); ++d)
98 for (
unsigned int i=0; i<face_embeddings[d].m(); ++i)
100 for (
unsigned int j=0; j<face_embeddings[d].n(); ++j)
119 std::ostringstream namebuf;
121 namebuf <<
"FE_ABF<" << dim <<
">(" <<
rt_order <<
")";
123 return namebuf.str();
148 const unsigned int n_interior_points = cell_quadrature.
size();
150 unsigned int n_face_points = (dim>1) ? 1 : 0;
152 for (
unsigned int d=1; d<dim; ++d)
153 n_face_points *= deg+1;
156 + n_interior_points);
162 std::vector<AnisotropicPolynomials<dim>* > polynomials_abf(dim);
165 for (
unsigned int dd=0; dd<dim; ++dd)
167 std::vector<std::vector<Polynomials::Polynomial<double> > > poly(dim);
168 for (
unsigned int d=0; d<dim; ++d)
176 unsigned int current = 0;
180 QGauss<dim-1> face_points (deg+1);
189 for (
unsigned int k=0; k<n_face_points; ++k)
195 for (
unsigned int i=0; i<legendre.n(); ++i)
198 = face_points.weight(k)
199 * legendre.compute_value(i, face_points.point(k));
204 for (; current<GeometryInfo<dim>::faces_per_cell*n_face_points;
219 for (
unsigned int k=0; k < faces.
size(); ++k)
221 for (
unsigned int i=0; i<polynomials_abf[0]->n() * dim; ++i)
224 compute_value(i / dim, faces.
point(k)) * faces.
weight(k);
233 std::vector<AnisotropicPolynomials<dim>* > polynomials(dim);
235 for (
unsigned int dd=0; dd<dim; ++dd)
237 std::vector<std::vector<Polynomials::Polynomial<double> > > poly(dim);
238 for (
unsigned int d=0; d<dim; ++d)
247 for (
unsigned int k=0; k<cell_quadrature.
size(); ++k)
249 for (
unsigned int i=0; i<polynomials[0]->n(); ++i)
250 for (
unsigned int d=0; d<dim; ++d)
252 * polynomials[d]->compute_value(i,cell_quadrature.
point(k));
255 for (
unsigned int d=0; d<dim; ++d)
256 delete polynomials[d];
262 for (
unsigned int k=0; k<cell_quadrature.
size(); ++k)
270 polynomials_abf[0]->n() * dim, dim));
273 for (
unsigned int k=0; k<cell_quadrature.
size(); ++k)
275 for (
unsigned int i=0; i<polynomials_abf[0]->n() * dim; ++i)
277 poly_grad = polynomials_abf[i%dim]->compute_grad(i / dim,cell_quadrature.
point(k))
278 * cell_quadrature.
weight(k);
280 for (
unsigned int d=0; d<dim; ++d)
285 for (
unsigned int d=0; d<dim; ++d)
286 delete polynomials_abf[d];
311 for (
unsigned int i=0; i<GeometryInfo<dim>::max_children_per_cell; ++i)
317 const unsigned int n_face_points = q_base.size();
320 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
332 for (
unsigned int k=0; k<q_face.
size(); ++k)
334 cached_values_face(i,k)
338 for (
unsigned int sub=0; sub<GeometryInfo<dim>::max_children_per_face; ++sub)
346 const unsigned int child
362 for (
unsigned int k=0; k<n_face_points; ++k)
363 for (
unsigned int i_child = 0; i_child < this->
dofs_per_cell; ++i_child)
364 for (
unsigned int i_face = 0; i_face < this->
dofs_per_face; ++i_face)
371 this->
restriction[iso][child](face*this->dofs_per_face+i_face,
373 += Utilities::fixed_power<dim-1>(.5) * q_sub.
weight(k)
374 * cached_values_face(i_child, k)
387 std::vector<AnisotropicPolynomials<dim>* > polynomials(dim);
388 for (
unsigned int dd=0; dd<dim; ++dd)
390 std::vector<std::vector<Polynomials::Polynomial<double> > > poly(dim);
391 for (
unsigned int d=0; d<dim; ++d)
399 const unsigned int start_cell_dofs
406 for (
unsigned int k=0; k<q_cell.size(); ++k)
408 for (
unsigned int d=0; d<dim; ++d)
411 for (
unsigned int child=0; child<GeometryInfo<dim>::max_children_per_cell; ++child)
415 for (
unsigned int k=0; k<q_sub.
size(); ++k)
416 for (
unsigned int i_child = 0; i_child < this->
dofs_per_cell; ++i_child)
417 for (
unsigned int d=0; d<dim; ++d)
418 for (
unsigned int i_weight=0; i_weight<polynomials[d]->n(); ++i_weight)
420 this->
restriction[iso][child](start_cell_dofs+i_weight*dim+d,
423 * cached_values_cell(i_child, k, d)
424 * polynomials[d]->compute_value(i_weight, q_sub.
point(k));
428 for (
unsigned int d=0; d<dim; ++d)
429 delete polynomials[d];
435 std::vector<unsigned int>
441 return std::vector<unsigned int>();
450 for (
int d=0; d<dim-1; ++d)
451 dofs_per_face *= rt_order+1;
455 interior_dofs = dim*(rt_order+1)*dofs_per_face;
457 std::vector<unsigned int> dpo(dim+1);
459 dpo[dim] = interior_dofs;
471 const unsigned int face_index)
const 515 std::vector<double> &nodal_values)
const 524 std::fill(nodal_values.begin(), nodal_values.end(), 0.);
527 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
528 for (
unsigned int k=0; k<n_face_points; ++k)
540 for (
unsigned int d=0; d<dim; ++d)
541 nodal_values[start_cell_dofs+i*dim+d] +=
interior_weights(k,i,d) * support_point_values[k+start_cell_points][d];
548 for (
unsigned int d=0; d<dim; ++d)
549 nodal_values[start_abf_dofs+i] +=
interior_weights_abf(k,i,d) * support_point_values[k+start_cell_points][d];
552 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
555 for (
unsigned int fp=0; fp < n_face_points; ++fp)
567 if (std::fabs (nodal_values[start_abf_dofs+i]) < 1.0e-16)
568 nodal_values[start_abf_dofs+i] = 0.0;
576 std::vector<double> &,
577 const std::vector<double> &)
const 587 std::vector<double> &local_dofs,
589 const unsigned int offset)
const 598 std::fill(local_dofs.begin(), local_dofs.end(), 0.);
601 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
602 for (
unsigned int k=0; k<n_face_points; ++k)
614 for (
unsigned int d=0; d<dim; ++d)
615 local_dofs[start_cell_dofs+i*dim+d] +=
interior_weights(k,i,d) * values[k+start_cell_points](d+offset);
622 for (
unsigned int d=0; d<dim; ++d)
623 local_dofs[start_abf_dofs+i] +=
interior_weights_abf(k,i,d) * values[k+start_cell_points](d+offset);
626 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
629 for (
unsigned int fp=0; fp < n_face_points; ++fp)
641 if (std::fabs (local_dofs[start_abf_dofs+i]) < 1.0e-16)
642 local_dofs[start_abf_dofs+i] = 0.0;
650 std::vector<double> &local_dofs,
651 const VectorSlice<
const std::vector<std::vector<double> > > &values)
const 660 std::fill(local_dofs.begin(), local_dofs.end(), 0.);
663 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
664 for (
unsigned int k=0; k<n_face_points; ++k)
676 for (
unsigned int d=0; d<dim; ++d)
677 local_dofs[start_cell_dofs+i*dim+d] +=
interior_weights(k,i,d) * values[d][k+start_cell_points];
684 for (
unsigned int d=0; d<dim; ++d)
688 for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
691 for (
unsigned int fp=0; fp < n_face_points; ++fp)
703 if (std::fabs (local_dofs[start_abf_dofs+i]) < 1.0e-16)
704 local_dofs[start_abf_dofs+i] = 0.0;
720 #include "fe_abf.inst" 722 DEAL_II_NAMESPACE_CLOSE
void initialize_restriction()
static Quadrature< dim > project_to_child(const Quadrature< dim > &quadrature, const unsigned int child_no)
std::vector< std::vector< FullMatrix< double > > > restriction
std::vector< Point< dim > > generalized_support_points
Table< 3, double > interior_weights
FullMatrix< double > interface_constraints
Table< 2, double > boundary_weights
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
virtual void interpolate(std::vector< double > &local_dofs, const std::vector< double > &values) const 1
virtual std::size_t memory_consumption() const
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const Point< dim > & point(const unsigned int i) const
Table< 2, double > boundary_weights_abf
virtual void convert_generalized_support_point_values_to_nodal_values(const std::vector< Vector< double > > &support_point_values, std::vector< double > &nodal_values) const
unsigned int size() const
static::ExceptionBase & ExcImpossibleInDim(int arg1)
std::vector< std::vector< FullMatrix< double > > > prolongation
void invert(const FullMatrix< number2 > &M)
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim-1 > &face_refinement_case=RefinementCase< dim-1 >::isotropic_refinement)
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
Table< 3, double > interior_weights_abf
#define Assert(cond, exc)
static void project_to_subface(const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim > > &q_points, const RefinementCase< dim-1 > &ref_case=RefinementCase< dim-1 >::isotropic_refinement)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
double weight(const unsigned int i) const
static std::vector< Polynomial< number > > generate_complete_basis(const unsigned int degree)
unsigned int n_components() const
const unsigned int dofs_per_cell
std::vector< Point< dim-1 > > generalized_face_support_points
static Quadrature< dim > project_to_all_faces(const SubQuadrature &quadrature)
const unsigned int rt_order
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim > > &q_points)
const unsigned int dofs_per_face
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
static::ExceptionBase & ExcNotImplemented()
virtual std::string get_name() const
unsigned int size(const unsigned int i) const
void initialize_support_points(const unsigned int rt_degree)
FullMatrix< double > inverse_node_matrix
virtual FiniteElement< dim > * clone() const
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
static::ExceptionBase & ExcInternalError()
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
static DataSetDescriptor face(const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)