17 #include <deal.II/base/quadrature_lib.h> 18 #include <deal.II/base/qprojector.h> 19 #include <deal.II/base/table.h> 20 #include <deal.II/grid/tria.h> 21 #include <deal.II/grid/tria_iterator.h> 22 #include <deal.II/dofs/dof_accessor.h> 23 #include <deal.II/fe/fe.h> 24 #include <deal.II/fe/mapping.h> 25 #include <deal.II/fe/fe_raviart_thomas.h> 26 #include <deal.II/fe/fe_values.h> 27 #include <deal.II/fe/fe_tools.h> 32 DEAL_II_NAMESPACE_OPEN
44 std::vector<bool>(dim,true)))
71 ref_case<RefinementCase<dim>::isotropic_refinement+1; ++ref_case)
75 for (
unsigned int i=0; i<nc; ++i)
76 this->
prolongation[ref_case-1][i].reinit (n_dofs, n_dofs);
82 for (
unsigned int i=0; i<GeometryInfo<dim>::max_children_per_face; ++i)
84 FETools::compute_face_embedding_matrices<dim,double>(*
this, face_embeddings, 0, 0);
87 unsigned int target_row=0;
88 for (
unsigned int d=0; d<GeometryInfo<dim>::max_children_per_face; ++d)
89 for (
unsigned int i=0; i<face_embeddings[d].m(); ++i)
91 for (
unsigned int j=0; j<face_embeddings[d].n(); ++j)
114 std::ostringstream namebuf;
115 namebuf <<
"FE_RaviartThomasNodal<" << dim <<
">(" << this->
degree-1 <<
")";
117 return namebuf.str();
143 unsigned int current = 0;
153 QGauss<dim-1> face_points (deg+1);
159 for (
unsigned int k=0;
163 ::DataSetDescriptor::face(0,
167 this->dofs_per_face));
169 current = this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;
179 for (
unsigned int d=0; d<dim; ++d)
184 ((d==1) ? low : high));
186 ((d==1) ? low : high),
187 ((d==2) ? low : high));
190 for (
unsigned int k=0; k<quadrature->
size(); ++k)
200 std::vector<unsigned int>
206 for (
unsigned int d=1; d<dim; ++d)
207 dofs_per_face *= deg+1;
213 std::vector<unsigned int> dpo(dim+1);
215 dpo[dim] = interior_dofs;
227 return std::vector<bool>();
238 for (
unsigned int d=2; d<dim; ++d)
239 dofs_per_face *= deg+1;
245 std::vector<bool> ret_val(dofs_per_cell,
false);
257 const unsigned int shape_index,
258 const unsigned int face_index)
const 268 const unsigned int support_face = shape_index / this->
degree;
288 std::vector<double> &nodal_values)
const 301 unsigned int fbase = 0;
303 for (; f<GeometryInfo<dim>::faces_per_cell;
314 const unsigned int istep = (this->
dofs_per_cell - fbase) / dim;
320 for (
unsigned int i=0; i<istep; ++i)
322 nodal_values[fbase+i] = support_point_values[fbase+i](f);
335 std::vector<double> &,
336 const std::vector<double> &)
const 345 std::vector<double> &local_dofs,
347 const unsigned int offset)
const 360 unsigned int fbase = 0;
362 for (; f<GeometryInfo<dim>::faces_per_cell;
373 const unsigned int istep = (this->
dofs_per_cell - fbase) / dim;
379 for (
unsigned int i=0; i<istep; ++i)
381 local_dofs[fbase+i] = values[fbase+i](offset+f);
393 std::vector<double> &local_dofs,
394 const VectorSlice<
const std::vector<std::vector<double> > > &values)
const 406 unsigned int fbase = 0;
408 for (; f<GeometryInfo<dim>::faces_per_cell;
418 const unsigned int istep = (this->
dofs_per_cell - fbase) / dim;
424 for (
unsigned int i=0; i<istep; ++i)
426 local_dofs[fbase+i] = values[f][fbase+i];
446 std::vector<std::pair<unsigned int, unsigned int> >
456 return std::vector<std::pair<unsigned int, unsigned int> > ();
460 return std::vector<std::pair<unsigned int, unsigned int> > ();
467 std::vector<std::pair<unsigned int, unsigned int> >
480 return std::vector<std::pair<unsigned int, unsigned int> >();
501 const unsigned int p = this->
degree-1;
502 const unsigned int q = fe_q_other->degree-1;
504 std::vector<std::pair<unsigned int, unsigned int> > identities;
507 for (
unsigned int i=0; i<p+1; ++i)
508 identities.push_back (std::make_pair(i,i));
510 else if (p%2==0 && q%2==0)
511 identities.push_back(std::make_pair(p/2,q/2));
518 return std::vector<std::pair<unsigned int, unsigned int> > ();
524 std::vector<std::pair<unsigned int, unsigned int> >
537 return std::vector<std::pair<unsigned int, unsigned int> >();
542 const unsigned int q = fe_q_other->dofs_per_quad;
544 std::vector<std::pair<unsigned int, unsigned int> > identities;
547 for (
unsigned int i=0; i<p; ++i)
548 identities.push_back (std::make_pair(i,i));
550 else if (p%2!=0 && q%2!=0)
551 identities.push_back(std::make_pair(p/2, q/2));
558 return std::vector<std::pair<unsigned int, unsigned int> > ();
571 if (this->degree < fe_q_other->
degree)
573 else if (this->degree == fe_q_other->degree)
662 double eps = 2e-13*this->
degree*(dim-1);
672 const Point<dim> &p = face_projection.point (i);
686 if ( std::fabs(matrix_entry - 1.0) < eps )
688 if ( std::fabs(matrix_entry) < eps )
691 interpolation_matrix(i,j) = matrix_entry;
705 sum += interpolation_matrix(j,i);
717 const unsigned int subface,
769 double eps = 2e-13*this->
degree*(dim-1);
780 const Point<dim> &p = subface_projection.point (i);
793 if ( std::fabs(matrix_entry - 1.0) < eps )
795 if ( std::fabs(matrix_entry) < eps )
798 interpolation_matrix(i,j) = matrix_entry;
812 sum += interpolation_matrix(j,i);
822 #include "fe_raviart_thomas_nodal.inst" 825 DEAL_II_NAMESPACE_CLOSE
virtual std::string get_name() const
std::vector< Point< dim > > generalized_support_points
FullMatrix< double > interface_constraints
const std::vector< Point< dim-1 > > & get_generalized_face_support_points() const
virtual FiniteElement< dim > * clone() const
FE_RaviartThomasNodal(const unsigned int p)
virtual void interpolate(std::vector< double > &local_dofs, const std::vector< double > &values) const 1
const unsigned int dofs_per_quad
const unsigned int degree
static unsigned int compute_n_pols(unsigned int degree)
#define AssertThrow(cond, exc)
static::ExceptionBase & ExcInterpolationNotImplemented()
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const Point< dim > & point(const unsigned int i) const
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
unsigned int size() const
static::ExceptionBase & ExcImpossibleInDim(int arg1)
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
std::vector< std::vector< FullMatrix< double > > > prolongation
void invert(const FullMatrix< number2 > &M)
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#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)
void initialize_support_points(const unsigned int rt_degree)
virtual std::string get_name() const =0
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)
virtual bool hp_constraints_are_implemented() const
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
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::ExceptionBase & ExcNotImplemented()
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
FullMatrix< double > inverse_node_matrix
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
virtual void convert_generalized_support_point_values_to_nodal_values(const std::vector< Vector< double > > &support_point_values, std::vector< double > &nodal_values) const
static::ExceptionBase & ExcInternalError()
static std::vector< bool > get_ria_vector(const unsigned int degree)