#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
#include <fe/fe_q.h>
-#include <fe/fe_values.h>
+#include <fe/hp_fe_values.h>
#include <fe/mapping_q.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
void solve ();
Triangulation<dim> triangulation;
- FE_Q<dim> fe;
- DoFHandler<dim> dof_handler;
- MappingQ<dim> mapping;
+ hp::FECollection<dim> fe;
+ hp::DoFHandler<dim> dof_handler;
+ hp::MappingCollection<dim> mapping;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
template <int dim>
LaplaceProblem<dim>::LaplaceProblem (const unsigned int mapping_degree) :
- fe (1),
+ fe (FE_Q<dim>(1)),
dof_handler (triangulation),
- mapping (mapping_degree)
+ mapping (MappingQ<dim>(mapping_degree))
{
deallog << "Using mapping with degree " << mapping_degree << ":"
<< std::endl
{
const unsigned int gauss_degree
- = std::max (static_cast<unsigned int>(std::ceil(1.*(mapping.get_degree()+1)/2)),
+ = std::max (static_cast<unsigned int>(std::ceil(1.*(static_cast<const MappingQ<dim>&>(mapping[0]).get_degree()+1)/2)),
2U);
MatrixTools::create_laplace_matrix (mapping, dof_handler,
QGauss<dim>(gauss_degree),
#include <grid/grid_refinement.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
-#include <fe/fe_values.h>
+#include <fe/hp_fe_values.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
public:
DGTransportEquation();
- void assemble_cell_term(const FEValues<dim>& fe_v,
+ void assemble_cell_term(const hp::FEValues<dim>& fe_v,
FullMatrix<double> &ui_vi_matrix,
Vector<double> &cell_vector) const;
- void assemble_boundary_term(const FEFaceValues<dim>& fe_v,
+ void assemble_boundary_term(const hp::FEFaceValues<dim>& fe_v,
FullMatrix<double> &ui_vi_matrix,
Vector<double> &cell_vector) const;
-
- void assemble_face_term1(const FEFaceValuesBase<dim>& fe_v,
- const FEFaceValuesBase<dim>& fe_v_neighbor,
+
+ template <class X, class Y>
+ void assemble_face_term1(const X& fe_v,
+ const Y& fe_v_neighbor,
FullMatrix<double> &ui_vi_matrix,
FullMatrix<double> &ue_vi_matrix) const;
- void assemble_face_term2(const FEFaceValuesBase<dim>& fe_v,
- const FEFaceValuesBase<dim>& fe_v_neighbor,
+ template <class X, class Y>
+ void assemble_face_term2(const X& fe_v,
+ const Y& fe_v_neighbor,
FullMatrix<double> &ui_vi_matrix,
FullMatrix<double> &ue_vi_matrix,
FullMatrix<double> &ui_ve_matrix,
template <int dim>
void DGTransportEquation<dim>::assemble_cell_term(
- const FEValues<dim> &fe_v,
+ const hp::FEValues<dim> &fe_v,
FullMatrix<double> &ui_vi_matrix,
Vector<double> &cell_vector) const
{
- const std::vector<double> &JxW = fe_v.get_JxW_values ();
+ const std::vector<double> &JxW = fe_v.get_present_fe_values().get_JxW_values ();
- std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
- std::vector<double> rhs (fe_v.n_quadrature_points);
+ std::vector<Point<dim> > beta (fe_v.get_present_fe_values().n_quadrature_points);
+ std::vector<double> rhs (fe_v.get_present_fe_values().n_quadrature_points);
- beta_function.value_list (fe_v.get_quadrature_points(), beta);
- rhs_function.value_list (fe_v.get_quadrature_points(), rhs);
+ beta_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), beta);
+ rhs_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), rhs);
- for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
- for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
+ for (unsigned int point=0; point<fe_v.get_present_fe_values().n_quadrature_points; ++point)
+ for (unsigned int i=0; i<fe_v.get_present_fe_values().dofs_per_cell; ++i)
{
- for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
- ui_vi_matrix(i,j) -= beta[point]*fe_v.shape_grad(i,point)*
- fe_v.shape_value(j,point) *
+ for (unsigned int j=0; j<fe_v.get_present_fe_values().dofs_per_cell; ++j)
+ ui_vi_matrix(i,j) -= beta[point]*fe_v.get_present_fe_values().shape_grad(i,point)*
+ fe_v.get_present_fe_values().shape_value(j,point) *
JxW[point];
- cell_vector(i) += rhs[point] * fe_v.shape_value(i,point) * JxW[point];
+ cell_vector(i) += rhs[point] * fe_v.get_present_fe_values().shape_value(i,point) * JxW[point];
}
}
template <int dim>
void DGTransportEquation<dim>::assemble_boundary_term(
- const FEFaceValues<dim>& fe_v,
+ const hp::FEFaceValues<dim>& fe_v,
FullMatrix<double> &ui_vi_matrix,
Vector<double> &cell_vector) const
{
- const std::vector<double> &JxW = fe_v.get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+ const std::vector<double> &JxW = fe_v.get_present_fe_values().get_JxW_values ();
+ const std::vector<Point<dim> > &normals = fe_v.get_present_fe_values().get_normal_vectors ();
- std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
- std::vector<double> g(fe_v.n_quadrature_points);
+ std::vector<Point<dim> > beta (fe_v.get_present_fe_values().n_quadrature_points);
+ std::vector<double> g(fe_v.get_present_fe_values().n_quadrature_points);
- beta_function.value_list (fe_v.get_quadrature_points(), beta);
- boundary_function.value_list (fe_v.get_quadrature_points(), g);
+ beta_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), beta);
+ boundary_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), g);
- for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
+ for (unsigned int point=0; point<fe_v.get_present_fe_values().n_quadrature_points; ++point)
{
const double beta_n=beta[point] * normals[point];
if (beta_n>0)
- for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
- for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+ for (unsigned int i=0; i<fe_v.get_present_fe_values().dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe_v.get_present_fe_values().dofs_per_cell; ++j)
ui_vi_matrix(i,j) += beta_n *
- fe_v.shape_value(j,point) *
- fe_v.shape_value(i,point) *
+ fe_v.get_present_fe_values().shape_value(j,point) *
+ fe_v.get_present_fe_values().shape_value(i,point) *
JxW[point];
else
- for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
+ for (unsigned int i=0; i<fe_v.get_present_fe_values().dofs_per_cell; ++i)
cell_vector(i) -= beta_n *
g[point] *
- fe_v.shape_value(i,point) *
+ fe_v.get_present_fe_values().shape_value(i,point) *
JxW[point];
}
}
template <int dim>
+template <class X, class Y>
void DGTransportEquation<dim>::assemble_face_term1(
- const FEFaceValuesBase<dim>& fe_v,
- const FEFaceValuesBase<dim>& fe_v_neighbor,
+ const X& fe_v,
+ const Y& fe_v_neighbor,
FullMatrix<double> &ui_vi_matrix,
FullMatrix<double> &ue_vi_matrix) const
{
- const std::vector<double> &JxW = fe_v.get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+ const std::vector<double> &JxW = fe_v.get_present_fe_values().get_JxW_values ();
+ const std::vector<Point<dim> > &normals = fe_v.get_present_fe_values().get_normal_vectors ();
- std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
- beta_function.value_list (fe_v.get_quadrature_points(), beta);
+ std::vector<Point<dim> > beta (fe_v.get_present_fe_values().n_quadrature_points);
+ beta_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), beta);
- for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
+ for (unsigned int point=0; point<fe_v.get_present_fe_values().n_quadrature_points; ++point)
{
const double beta_n=beta[point] * normals[point];
if (beta_n>0)
- for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
- for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+ for (unsigned int i=0; i<fe_v.get_present_fe_values().dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe_v.get_present_fe_values().dofs_per_cell; ++j)
ui_vi_matrix(i,j) += beta_n *
- fe_v.shape_value(j,point) *
- fe_v.shape_value(i,point) *
+ fe_v.get_present_fe_values().shape_value(j,point) *
+ fe_v.get_present_fe_values().shape_value(i,point) *
JxW[point];
else
- for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
- for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
+ for (unsigned int i=0; i<fe_v.get_present_fe_values().dofs_per_cell; ++i)
+ for (unsigned int k=0; k<fe_v_neighbor.get_present_fe_values().dofs_per_cell; ++k)
ue_vi_matrix(i,k) += beta_n *
- fe_v_neighbor.shape_value(k,point) *
- fe_v.shape_value(i,point) *
+ fe_v_neighbor.get_present_fe_values().shape_value(k,point) *
+ fe_v.get_present_fe_values().shape_value(i,point) *
JxW[point];
}
}
template <int dim>
+template <class X, class Y>
void DGTransportEquation<dim>::assemble_face_term2(
- const FEFaceValuesBase<dim>& fe_v,
- const FEFaceValuesBase<dim>& fe_v_neighbor,
+ const X& fe_v,
+ const Y& fe_v_neighbor,
FullMatrix<double> &ui_vi_matrix,
FullMatrix<double> &ue_vi_matrix,
FullMatrix<double> &ui_ve_matrix,
FullMatrix<double> &ue_ve_matrix) const
{
- const std::vector<double> &JxW = fe_v.get_JxW_values ();
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+ const std::vector<double> &JxW = fe_v.get_present_fe_values().get_JxW_values ();
+ const std::vector<Point<dim> > &normals = fe_v.get_present_fe_values().get_normal_vectors ();
- std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
+ std::vector<Point<dim> > beta (fe_v.get_present_fe_values().n_quadrature_points);
- beta_function.value_list (fe_v.get_quadrature_points(), beta);
+ beta_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), beta);
- for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
+ for (unsigned int point=0; point<fe_v.get_present_fe_values().n_quadrature_points; ++point)
{
const double beta_n=beta[point] * normals[point];
if (beta_n>0)
{
- for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
- for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+ for (unsigned int i=0; i<fe_v.get_present_fe_values().dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe_v.get_present_fe_values().dofs_per_cell; ++j)
ui_vi_matrix(i,j) += beta_n *
- fe_v.shape_value(j,point) *
- fe_v.shape_value(i,point) *
+ fe_v.get_present_fe_values().shape_value(j,point) *
+ fe_v.get_present_fe_values().shape_value(i,point) *
JxW[point];
- for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
- for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+ for (unsigned int k=0; k<fe_v_neighbor.get_present_fe_values().dofs_per_cell; ++k)
+ for (unsigned int j=0; j<fe_v.get_present_fe_values().dofs_per_cell; ++j)
ui_ve_matrix(k,j) -= beta_n *
- fe_v.shape_value(j,point) *
- fe_v_neighbor.shape_value(k,point) *
+ fe_v.get_present_fe_values().shape_value(j,point) *
+ fe_v_neighbor.get_present_fe_values().shape_value(k,point) *
JxW[point];
}
else
{
- for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
- for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
+ for (unsigned int i=0; i<fe_v.get_present_fe_values().dofs_per_cell; ++i)
+ for (unsigned int l=0; l<fe_v_neighbor.get_present_fe_values().dofs_per_cell; ++l)
ue_vi_matrix(i,l) += beta_n *
- fe_v_neighbor.shape_value(l,point) *
- fe_v.shape_value(i,point) *
+ fe_v_neighbor.get_present_fe_values().shape_value(l,point) *
+ fe_v.get_present_fe_values().shape_value(i,point) *
JxW[point];
- for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
- for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
+ for (unsigned int k=0; k<fe_v_neighbor.get_present_fe_values().dofs_per_cell; ++k)
+ for (unsigned int l=0; l<fe_v_neighbor.get_present_fe_values().dofs_per_cell; ++l)
ue_ve_matrix(k,l) -= beta_n *
- fe_v_neighbor.shape_value(l,point) *
- fe_v_neighbor.shape_value(k,point) *
+ fe_v_neighbor.get_present_fe_values().shape_value(l,point) *
+ fe_v_neighbor.get_present_fe_values().shape_value(k,point) *
JxW[point];
}
}
void output_results (const unsigned int cycle) const;
Triangulation<dim> triangulation;
- const MappingQ1<dim> mapping;
+ const hp::MappingCollection<dim> mapping;
- FE_DGQ<dim> fe;
- DoFHandler<dim> dof_handler;
+ hp::FECollection<dim> fe;
+ hp::DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
- const QGauss<dim> quadrature;
- const QGauss<dim-1> face_quadrature;
+ const hp::QCollection<dim> quadrature;
+ const hp::QCollection<dim-1> face_quadrature;
Vector<double> solution1;
Vector<double> solution2;
template <int dim>
DGMethod<dim>::DGMethod ()
:
- mapping (),
- fe (1),
+ mapping (MappingQ1<dim>()),
+ fe (FE_DGQ<dim>(1)),
dof_handler (triangulation),
- quadrature (4),
- face_quadrature (4),
+ quadrature (QGauss<dim>(4)),
+ face_quadrature (QGauss<dim-1>(4)),
dg ()
{}
sparsity_pattern.reinit (dof_handler.n_dofs(),
dof_handler.n_dofs(),
(GeometryInfo<dim>::faces_per_cell
- *GeometryInfo<dim>::subfaces_per_face+1)*fe.dofs_per_cell);
+ *GeometryInfo<dim>::subfaces_per_face+1)*fe[0].dofs_per_cell);
DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern);
template <int dim>
void DGMethod<dim>::assemble_system1 ()
{
- const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+ const unsigned int dofs_per_cell = dof_handler.get_fe()[0].dofs_per_cell;
std::vector<unsigned int> dofs (dofs_per_cell);
std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
const UpdateFlags neighbor_face_update_flags = update_values;
- FEValues<dim> fe_v (
+ hp::FEValues<dim> fe_v (
mapping, fe, quadrature, update_flags);
- FEFaceValues<dim> fe_v_face (
+ hp::FEFaceValues<dim> fe_v_face (
mapping, fe, face_quadrature, face_update_flags);
- FESubfaceValues<dim> fe_v_subface (
+ hp::FESubfaceValues<dim> fe_v_subface (
mapping, fe, face_quadrature, face_update_flags);
- FEFaceValues<dim> fe_v_face_neighbor (
+ hp::FEFaceValues<dim> fe_v_face_neighbor (
mapping, fe, face_quadrature, neighbor_face_update_flags);
- FESubfaceValues<dim> fe_v_subface_neighbor (
+ hp::FESubfaceValues<dim> fe_v_subface_neighbor (
mapping, fe, face_quadrature, neighbor_face_update_flags);
FullMatrix<double> ui_vi_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_vector (dofs_per_cell);
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
{
- typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
+ typename hp::DoFHandler<dim>::face_iterator face=cell->face(face_no);
ue_vi_matrix = 0;
}
else
{
- typename DoFHandler<dim>::cell_iterator neighbor=
+ typename hp::DoFHandler<dim>::cell_iterator neighbor=
cell->neighbor(face_no);;
if (face->has_children())
for (unsigned int subface_no=0;
subface_no<face->n_children(); ++subface_no)
{
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
neighbor_child
= cell->neighbor_child_on_subface (face_no, subface_no);
template <int dim>
void DGMethod<dim>::assemble_system2 ()
{
- const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+ const unsigned int dofs_per_cell = dof_handler.get_fe()[0].dofs_per_cell;
std::vector<unsigned int> dofs (dofs_per_cell);
std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
const UpdateFlags neighbor_face_update_flags = update_values;
- FEValues<dim> fe_v (
+ hp::FEValues<dim> fe_v (
mapping, fe, quadrature, update_flags);
- FEFaceValues<dim> fe_v_face (
+ hp::FEFaceValues<dim> fe_v_face (
mapping, fe, face_quadrature, face_update_flags);
- FESubfaceValues<dim> fe_v_subface (
+ hp::FESubfaceValues<dim> fe_v_subface (
mapping, fe, face_quadrature, face_update_flags);
- FEFaceValues<dim> fe_v_face_neighbor (
+ hp::FEFaceValues<dim> fe_v_face_neighbor (
mapping, fe, face_quadrature, neighbor_face_update_flags);
Vector<double> cell_vector (dofs_per_cell);
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (;cell!=endc; ++cell)
for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
{
- typename DoFHandler<dim>::face_iterator face=
+ typename hp::DoFHandler<dim>::face_iterator face=
cell->face(face_no);
if (face->at_boundary())
{
Assert (cell->neighbor(face_no).state() == IteratorState::valid,
ExcInternalError());
- typename DoFHandler<dim>::cell_iterator neighbor=
+ typename hp::DoFHandler<dim>::cell_iterator neighbor=
cell->neighbor(face_no);
if (face->has_children())
{
for (unsigned int subface_no=0;
subface_no<face->n_children(); ++subface_no)
{
- typename DoFHandler<dim>::cell_iterator neighbor_child
+ typename hp::DoFHandler<dim>::cell_iterator neighbor_child
= cell->neighbor_child_on_subface (face_no, subface_no);
Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
ExcInternalError());
PreconditionBlockSSOR<SparseMatrix<double> > preconditioner;
- preconditioner.initialize(system_matrix, fe.dofs_per_cell);
+ preconditioner.initialize(system_matrix, fe[0].dofs_per_cell);
solver.solve (system_matrix, solution, right_hand_side,
preconditioner);
solution2,
gradient_indicator);
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no)
deallog << "Writing solution to <" << filename << ">..."
<< std::endl << std::endl;
- DataOut<dim> data_out;
+ DataOut<dim,hp::DoFHandler<dim> > data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution2, "u");
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
#include <fe/fe_q.h>
-#include <fe/fe_values.h>
+#include <fe/hp_fe_values.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
#include <numerics/data_out.h>
void set_refinement_cycle (const unsigned int refinement_cycle);
- virtual void operator () (const DoFHandler<dim> &dof_handler,
+ virtual void operator () (const hp::DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const = 0;
protected:
unsigned int refinement_cycle;
PointValueEvaluation (const Point<dim> &evaluation_point,
TableHandler &results_table);
- virtual void operator () (const DoFHandler<dim> &dof_handler,
+ virtual void operator () (const hp::DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const;
DeclException1 (ExcEvaluationPointNotFound,
template <int dim>
void
PointValueEvaluation<dim>::
- operator () (const DoFHandler<dim> &dof_handler,
+ operator () (const hp::DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const
{
double point_value = 1e20;
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
bool evaluation_point_found = false;
++vertex)
if (cell->vertex(vertex) == evaluation_point)
{
- point_value = solution(cell->vertex_dof_index(vertex,0));
+ point_value = solution(cell->vertex_dof_index(vertex,0,
+ cell->active_fe_index()));
evaluation_point_found = true;
break;
SolutionOutput (const std::string &output_name_base,
const typename DataOut<dim>::OutputFormat output_format);
- virtual void operator () (const DoFHandler<dim> &dof_handler,
+ virtual void operator () (const hp::DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const;
private:
const std::string output_name_base;
template <int dim>
void
- SolutionOutput<dim>::operator () (const DoFHandler<dim> &dof_handler,
+ SolutionOutput<dim>::operator () (const hp::DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const
{
- DataOut<dim> data_out;
+ DataOut<dim,hp::DoFHandler<dim> > data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, "solution");
data_out.build_patches ();
{
public:
Solver (Triangulation<dim> &triangulation,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
const Function<dim> &boundary_values);
virtual
~Solver ();
n_dofs () const;
protected:
- const SmartPointer<const FiniteElement<dim> > fe;
- const SmartPointer<const Quadrature<dim> > quadrature;
- DoFHandler<dim> dof_handler;
+ const SmartPointer<const hp::FECollection<dim> > fe;
+ const SmartPointer<const hp::QCollection<dim> > quadrature;
+ hp::DoFHandler<dim> dof_handler;
Vector<double> solution;
const SmartPointer<const Function<dim> > boundary_values;
private:
struct LinearSystem
{
- LinearSystem (const DoFHandler<dim> &dof_handler);
+ LinearSystem (const hp::DoFHandler<dim> &dof_handler);
void solve (Vector<double> &solution) const;
void
assemble_matrix (LinearSystem &linear_system,
- const typename DoFHandler<dim>::active_cell_iterator &begin_cell,
- const typename DoFHandler<dim>::active_cell_iterator &end_cell,
+ const typename hp::DoFHandler<dim>::active_cell_iterator &begin_cell,
+ const typename hp::DoFHandler<dim>::active_cell_iterator &end_cell,
Threads::ThreadMutex &mutex) const;
};
template <int dim>
Solver<dim>::Solver (Triangulation<dim> &triangulation,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
const Function<dim> &boundary_values)
:
Base<dim> (triangulation),
Solver<dim>::assemble_linear_system (LinearSystem &linear_system)
{
typedef
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
active_cell_iterator;
const unsigned int n_threads = multithread_info.n_default_threads;
template <int dim>
void
Solver<dim>::assemble_matrix (LinearSystem &linear_system,
- const typename DoFHandler<dim>::active_cell_iterator &begin_cell,
- const typename DoFHandler<dim>::active_cell_iterator &end_cell,
+ const typename hp::DoFHandler<dim>::active_cell_iterator &begin_cell,
+ const typename hp::DoFHandler<dim>::active_cell_iterator &end_cell,
Threads::ThreadMutex &mutex) const
{
- FEValues<dim> fe_values (*fe, *quadrature,
+ hp::FEValues<dim> fe_values (*fe, *quadrature,
update_gradients | update_JxW_values);
- const unsigned int dofs_per_cell = fe->dofs_per_cell;
- const unsigned int n_q_points = quadrature->n_quadrature_points;
+ const unsigned int dofs_per_cell = (*fe)[0].dofs_per_cell;
+ const unsigned int n_q_points = (*quadrature)[0].n_quadrature_points;
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
- for (typename DoFHandler<dim>::active_cell_iterator cell=begin_cell;
+ for (typename hp::DoFHandler<dim>::active_cell_iterator cell=begin_cell;
cell!=end_cell; ++cell)
{
cell_matrix = 0;
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
- fe_values.shape_grad(j,q_point) *
- fe_values.JxW(q_point));
+ cell_matrix(i,j) += (fe_values.get_present_fe_values().shape_grad(i,q_point) *
+ fe_values.get_present_fe_values().shape_grad(j,q_point) *
+ fe_values.get_present_fe_values().JxW(q_point));
cell->get_dof_indices (local_dof_indices);
template <int dim>
Solver<dim>::LinearSystem::
- LinearSystem (const DoFHandler<dim> &dof_handler)
+ LinearSystem (const hp::DoFHandler<dim> &dof_handler)
{
hanging_node_constraints.clear ();
- void (*mhnc_p) (const DoFHandler<dim> &,
+ void (*mhnc_p) (const hp::DoFHandler<dim> &,
ConstraintMatrix &)
= &DoFTools::make_hanging_node_constraints;
{
public:
PrimalSolver (Triangulation<dim> &triangulation,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values);
protected:
template <int dim>
PrimalSolver<dim>::
PrimalSolver (Triangulation<dim> &triangulation,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values)
:
PrimalSolver<dim>::
assemble_rhs (Vector<double> &rhs) const
{
- FEValues<dim> fe_values (*this->fe, *this->quadrature,
+ hp::FEValues<dim> fe_values (*this->fe, *this->quadrature,
update_values | update_q_points |
update_JxW_values);
- const unsigned int dofs_per_cell = this->fe->dofs_per_cell;
- const unsigned int n_q_points = this->quadrature->n_quadrature_points;
+ const unsigned int dofs_per_cell = (*this->fe)[0].dofs_per_cell;
+ const unsigned int n_q_points = (*this->quadrature)[0].n_quadrature_points;
Vector<double> cell_rhs (dofs_per_cell);
std::vector<double> rhs_values (n_q_points);
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
cell = this->dof_handler.begin_active(),
endc = this->dof_handler.end();
for (; cell!=endc; ++cell)
{
cell_rhs = 0;
fe_values.reinit (cell);
- rhs_function->value_list (fe_values.get_quadrature_points(),
+ rhs_function->value_list (fe_values.get_present_fe_values().get_quadrature_points(),
rhs_values);
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
- cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+ cell_rhs(i) += (fe_values.get_present_fe_values().shape_value(i,q_point) *
rhs_values[q_point] *
- fe_values.JxW(q_point));
+ fe_values.get_present_fe_values().JxW(q_point));
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
public:
RefinementGlobal (Triangulation<dim> &coarse_grid,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values);
template <int dim>
RefinementGlobal<dim>::
RefinementGlobal (Triangulation<dim> &coarse_grid,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values)
:
{
public:
RefinementKelly (Triangulation<dim> &coarse_grid,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values);
template <int dim>
RefinementKelly<dim>::
RefinementKelly (Triangulation<dim> &coarse_grid,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values)
:
Triangulation<dim> triangulation;
GridGenerator::hyper_cube (triangulation, -1, 1);
triangulation.refine_global (2);
- const FE_Q<dim> fe(1);
- const QGauss<dim> quadrature(4);
+ const hp::FECollection<dim> fe(FE_Q<dim>(1));
+ const hp::QCollection<dim> quadrature(QGauss<dim>(4));
const RightHandSide<dim> rhs_function;
const Solution<dim> boundary_values;
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
#include <fe/fe_q.h>
-#include <fe/fe_values.h>
+#include <fe/hp_fe_values.h>
#include <fe/fe_tools.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
void set_refinement_cycle (const unsigned int refinement_cycle);
- virtual void operator () (const DoFHandler<dim> &dof_handler,
+ virtual void operator () (const hp::DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const = 0;
protected:
unsigned int refinement_cycle;
public:
PointValueEvaluation (const Point<dim> &evaluation_point);
- virtual void operator () (const DoFHandler<dim> &dof_handler,
+ virtual void operator () (const hp::DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const;
DeclException1 (ExcEvaluationPointNotFound,
template <int dim>
void
PointValueEvaluation<dim>::
- operator () (const DoFHandler<dim> &dof_handler,
+ operator () (const hp::DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const
{
double point_value = 1e20;
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
bool evaluation_point_found = false;
public:
PointXDerivativeEvaluation (const Point<dim> &evaluation_point);
- virtual void operator () (const DoFHandler<dim> &dof_handler,
+ virtual void operator () (const hp::DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const;
DeclException1 (ExcEvaluationPointNotFound,
template <int dim>
void
PointXDerivativeEvaluation<dim>::
- operator () (const DoFHandler<dim> &dof_handler,
+ operator () (const hp::DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const
{
double point_derivative = 0;
QTrapez<dim> vertex_quadrature;
- FEValues<dim> fe_values (dof_handler.get_fe(),
+ hp::FEValues<dim> fe_values (dof_handler.get_fe(),
vertex_quadrature,
update_gradients | update_q_points);
std::vector<Tensor<1,dim> >
solution_gradients (vertex_quadrature.n_quadrature_points);
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
unsigned int evaluation_point_hits = 0;
if (cell->vertex(vertex) == evaluation_point)
{
fe_values.reinit (cell);
- fe_values.get_function_grads (solution,
+ fe_values.get_present_fe_values().get_function_grads (solution,
solution_gradients);
unsigned int q_point = 0;
for (; q_point<solution_gradients.size(); ++q_point)
- if (fe_values.quadrature_point(q_point) ==
+ if (fe_values.get_present_fe_values().quadrature_point(q_point) ==
evaluation_point)
break;
public:
GridOutput (const std::string &output_name_base);
- virtual void operator () (const DoFHandler<dim> &dof_handler,
+ virtual void operator () (const hp::DoFHandler<dim> &dof_handler,
const Vector<double> &solution) const;
private:
const std::string output_name_base;
template <int dim>
void
- GridOutput<dim>::operator () (const DoFHandler<dim> &dof_handler,
+ GridOutput<dim>::operator () (const hp::DoFHandler<dim> &dof_handler,
const Vector<double> &/*solution*/) const
{
std::ostringstream filename;
{
public:
Solver (Triangulation<dim> &triangulation,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
+ const hp::QCollection<dim-1> &face_quadrature,
const Function<dim> &boundary_values);
virtual
~Solver ();
n_dofs () const;
protected:
- const SmartPointer<const FiniteElement<dim> > fe;
- const SmartPointer<const Quadrature<dim> > quadrature;
- const SmartPointer<const Quadrature<dim-1> > face_quadrature;
- DoFHandler<dim> dof_handler;
+ const SmartPointer<const hp::FECollection<dim> > fe;
+ const SmartPointer<const hp::QCollection<dim> > quadrature;
+ const SmartPointer<const hp::QCollection<dim-1> > face_quadrature;
+ hp::DoFHandler<dim> dof_handler;
Vector<double> solution;
const SmartPointer<const Function<dim> > boundary_values;
private:
struct LinearSystem
{
- LinearSystem (const DoFHandler<dim> &dof_handler);
+ LinearSystem (const hp::DoFHandler<dim> &dof_handler);
void solve (Vector<double> &solution) const;
void
assemble_matrix (LinearSystem &linear_system,
- const typename DoFHandler<dim>::active_cell_iterator &begin_cell,
- const typename DoFHandler<dim>::active_cell_iterator &end_cell,
+ const typename hp::DoFHandler<dim>::active_cell_iterator &begin_cell,
+ const typename hp::DoFHandler<dim>::active_cell_iterator &end_cell,
Threads::ThreadMutex &mutex) const;
};
template <int dim>
Solver<dim>::Solver (Triangulation<dim> &triangulation,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
+ const hp::QCollection<dim-1> &face_quadrature,
const Function<dim> &boundary_values)
:
Base<dim> (triangulation),
Solver<dim>::assemble_linear_system (LinearSystem &linear_system)
{
typedef
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
active_cell_iterator;
const unsigned int n_threads = multithread_info.n_default_threads;
template <int dim>
void
Solver<dim>::assemble_matrix (LinearSystem &linear_system,
- const typename DoFHandler<dim>::active_cell_iterator &begin_cell,
- const typename DoFHandler<dim>::active_cell_iterator &end_cell,
+ const typename hp::DoFHandler<dim>::active_cell_iterator &begin_cell,
+ const typename hp::DoFHandler<dim>::active_cell_iterator &end_cell,
Threads::ThreadMutex &mutex) const
{
- FEValues<dim> fe_values (*fe, *quadrature,
+ hp::FEValues<dim> fe_values (*fe, *quadrature,
update_gradients | update_JxW_values);
- const unsigned int dofs_per_cell = fe->dofs_per_cell;
- const unsigned int n_q_points = quadrature->n_quadrature_points;
+ const unsigned int dofs_per_cell = (*fe)[0].dofs_per_cell;
+ const unsigned int n_q_points = (*quadrature)[0].n_quadrature_points;
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
- for (typename DoFHandler<dim>::active_cell_iterator cell=begin_cell;
+ for (typename hp::DoFHandler<dim>::active_cell_iterator cell=begin_cell;
cell!=end_cell; ++cell)
{
cell_matrix = 0;
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
- fe_values.shape_grad(j,q_point) *
- fe_values.JxW(q_point));
+ cell_matrix(i,j) += (fe_values.get_present_fe_values().shape_grad(i,q_point) *
+ fe_values.get_present_fe_values().shape_grad(j,q_point) *
+ fe_values.get_present_fe_values().JxW(q_point));
cell->get_dof_indices (local_dof_indices);
template <int dim>
Solver<dim>::LinearSystem::
- LinearSystem (const DoFHandler<dim> &dof_handler)
+ LinearSystem (const hp::DoFHandler<dim> &dof_handler)
{
hanging_node_constraints.clear ();
- void (*mhnc_p) (const DoFHandler<dim> &,
+ void (*mhnc_p) (const hp::DoFHandler<dim> &,
ConstraintMatrix &)
= &DoFTools::make_hanging_node_constraints;
{
public:
PrimalSolver (Triangulation<dim> &triangulation,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
+ const hp::QCollection<dim-1> &face_quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values);
template <int dim>
PrimalSolver<dim>::
PrimalSolver (Triangulation<dim> &triangulation,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
+ const hp::QCollection<dim-1> &face_quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values)
:
void
PrimalSolver<dim>::output_solution () const
{
- DataOut<dim> data_out;
+ DataOut<dim,hp::DoFHandler<dim> > data_out;
data_out.attach_dof_handler (this->dof_handler);
data_out.add_data_vector (this->solution, "solution");
data_out.build_patches ();
PrimalSolver<dim>::
assemble_rhs (Vector<double> &rhs) const
{
- FEValues<dim> fe_values (*this->fe, *this->quadrature,
+ hp::FEValues<dim> fe_values (*this->fe, *this->quadrature,
update_values | update_q_points |
update_JxW_values);
- const unsigned int dofs_per_cell = this->fe->dofs_per_cell;
- const unsigned int n_q_points = this->quadrature->n_quadrature_points;
+ const unsigned int dofs_per_cell = (*this->fe)[0].dofs_per_cell;
+ const unsigned int n_q_points = (*this->quadrature)[0].n_quadrature_points;
Vector<double> cell_rhs (dofs_per_cell);
std::vector<double> rhs_values (n_q_points);
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
cell = this->dof_handler.begin_active(),
endc = this->dof_handler.end();
for (; cell!=endc; ++cell)
fe_values.reinit (cell);
- rhs_function->value_list (fe_values.get_quadrature_points(),
+ rhs_function->value_list (fe_values.get_present_fe_values().get_quadrature_points(),
rhs_values);
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
- cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+ cell_rhs(i) += (fe_values.get_present_fe_values().shape_value(i,q_point) *
rhs_values[q_point] *
- fe_values.JxW(q_point));
+ fe_values.get_present_fe_values().JxW(q_point));
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
public:
RefinementGlobal (Triangulation<dim> &coarse_grid,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
+ const hp::QCollection<dim-1> &face_quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values);
template <int dim>
RefinementGlobal<dim>::
RefinementGlobal (Triangulation<dim> &coarse_grid,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
+ const hp::QCollection<dim-1> &face_quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values)
:
{
public:
RefinementKelly (Triangulation<dim> &coarse_grid,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
+ const hp::QCollection<dim-1> &face_quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values);
template <int dim>
RefinementKelly<dim>::
RefinementKelly (Triangulation<dim> &coarse_grid,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
+ const hp::QCollection<dim-1> &face_quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values)
:
{
public:
RefinementWeightedKelly (Triangulation<dim> &coarse_grid,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
+ const hp::QCollection<dim-1> &face_quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values,
const Function<dim> &weighting_function);
template <int dim>
RefinementWeightedKelly<dim>::
RefinementWeightedKelly (Triangulation<dim> &coarse_grid,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
+ const hp::QCollection<dim-1> &face_quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values,
const Function<dim> &weighting_function)
{
Vector<float> estimated_error (this->triangulation->n_active_cells());
KellyErrorEstimator<dim>::estimate (this->dof_handler,
- *this->face_quadrature,
+ (*this->face_quadrature)[0],
typename FunctionMap<dim>::type(),
this->solution,
estimated_error);
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
cell = this->dof_handler.begin_active(),
endc = this->dof_handler.end();
for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index)
public:
virtual
void
- assemble_rhs (const DoFHandler<dim> &dof_handler,
+ assemble_rhs (const hp::DoFHandler<dim> &dof_handler,
Vector<double> &rhs) const = 0;
};
virtual
void
- assemble_rhs (const DoFHandler<dim> &dof_handler,
+ assemble_rhs (const hp::DoFHandler<dim> &dof_handler,
Vector<double> &rhs) const;
DeclException1 (ExcEvaluationPointNotFound,
template <int dim>
void
PointValueEvaluation<dim>::
- assemble_rhs (const DoFHandler<dim> &dof_handler,
+ assemble_rhs (const hp::DoFHandler<dim> &dof_handler,
Vector<double> &rhs) const
{
rhs.reinit (dof_handler.n_dofs());
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
virtual
void
- assemble_rhs (const DoFHandler<dim> &dof_handler,
+ assemble_rhs (const hp::DoFHandler<dim> &dof_handler,
Vector<double> &rhs) const;
DeclException1 (ExcEvaluationPointNotFound,
template <int dim>
void
PointXDerivativeEvaluation<dim>::
- assemble_rhs (const DoFHandler<dim> &dof_handler,
+ assemble_rhs (const hp::DoFHandler<dim> &dof_handler,
Vector<double> &rhs) const
{
rhs.reinit (dof_handler.n_dofs());
QGauss<dim> quadrature(4);
- FEValues<dim> fe_values (dof_handler.get_fe(), quadrature,
+ hp::FEValues<dim> fe_values (dof_handler.get_fe(), quadrature,
update_gradients |
update_q_points |
update_JxW_values);
- const unsigned int n_q_points = fe_values.n_quadrature_points;
+ const unsigned int n_q_points = fe_values.get_present_fe_values().n_quadrature_points;
const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
Vector<double> cell_rhs (dofs_per_cell);
double total_volume = 0;
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int i=0; i<dofs_per_cell; ++i)
- cell_rhs(i) += fe_values.shape_grad(i,q)[0] *
- fe_values.JxW (q);
- total_volume += fe_values.JxW (q);
+ cell_rhs(i) += fe_values.get_present_fe_values().shape_grad(i,q)[0] *
+ fe_values.get_present_fe_values().JxW (q);
+ total_volume += fe_values.get_present_fe_values().JxW (q);
};
cell->get_dof_indices (local_dof_indices);
{
public:
DualSolver (Triangulation<dim> &triangulation,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
+ const hp::QCollection<dim-1> &face_quadrature,
const DualFunctional::DualFunctionalBase<dim> &dual_functional);
virtual
template <int dim>
DualSolver<dim>::
DualSolver (Triangulation<dim> &triangulation,
- const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
+ const hp::QCollection<dim-1> &face_quadrature,
const DualFunctional::DualFunctionalBase<dim> &dual_functional)
:
Base<dim> (triangulation),
{
public:
WeightedResidual (Triangulation<dim> &coarse_grid,
- const FiniteElement<dim> &primal_fe,
- const FiniteElement<dim> &dual_fe,
- const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const hp::FECollection<dim> &primal_fe,
+ const hp::FECollection<dim> &dual_fe,
+ const hp::QCollection<dim> &quadrature,
+ const hp::QCollection<dim-1> &face_quadrature,
const Function<dim> &rhs_function,
const Function<dim> &boundary_values,
const DualFunctional::DualFunctionalBase<dim> &dual_functional);
void solve_dual_problem ();
typedef
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
active_cell_iterator;
typedef
- typename std::map<typename DoFHandler<dim>::face_iterator,double>
+ typename std::map<typename hp::DoFHandler<dim>::face_iterator,double>
FaceIntegrals;
struct CellData
{
- FEValues<dim> fe_values;
+ hp::FEValues<dim> fe_values;
const SmartPointer<const Function<dim> > right_hand_side;
std::vector<double> cell_residual;
std::vector<double> rhs_values;
std::vector<double> dual_weights;
typename std::vector<Tensor<2,dim> > cell_grad_grads;
- CellData (const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
+ CellData (const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
const Function<dim> &right_hand_side);
};
struct FaceData
{
- FEFaceValues<dim> fe_face_values_cell;
- FEFaceValues<dim> fe_face_values_neighbor;
- FESubfaceValues<dim> fe_subface_values_cell;
+ hp::FEFaceValues<dim> fe_face_values_cell;
+ hp::FEFaceValues<dim> fe_face_values_neighbor;
+ hp::FESubfaceValues<dim> fe_subface_values_cell;
std::vector<double> jump_residual;
std::vector<double> dual_weights;
typename std::vector<Tensor<1,dim> > cell_grads;
typename std::vector<Tensor<1,dim> > neighbor_grads;
- FaceData (const FiniteElement<dim> &fe,
- const Quadrature<dim-1> &face_quadrature);
+ FaceData (const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim-1> &face_quadrature);
};
template <int dim>
WeightedResidual<dim>::CellData::
- CellData (const FiniteElement<dim> &fe,
- const Quadrature<dim> &quadrature,
+ CellData (const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim> &quadrature,
const Function<dim> &right_hand_side)
:
fe_values (fe, quadrature,
right_hand_side (&right_hand_side)
{
const unsigned int n_q_points
- = quadrature.n_quadrature_points;
+ = quadrature[0].n_quadrature_points;
cell_residual.resize(n_q_points);
rhs_values.resize(n_q_points);
template <int dim>
WeightedResidual<dim>::FaceData::
- FaceData (const FiniteElement<dim> &fe,
- const Quadrature<dim-1> &face_quadrature)
+ FaceData (const hp::FECollection<dim> &fe,
+ const hp::QCollection<dim-1> &face_quadrature)
:
fe_face_values_cell (fe, face_quadrature,
update_values |
update_gradients)
{
const unsigned int n_face_q_points
- = face_quadrature.n_quadrature_points;
+ = face_quadrature[0].n_quadrature_points;
jump_residual.resize(n_face_q_points);
dual_weights.resize(n_face_q_points);
template <int dim>
WeightedResidual<dim>::
WeightedResidual (Triangulation<dim> &coarse_grid,
- const FiniteElement<dim> &primal_fe,
- const FiniteElement<dim> &dual_fe,
- const Quadrature<dim> &quadrature,
- const Quadrature<dim-1> &face_quadrature,
+ const hp::FECollection<dim> &primal_fe,
+ const hp::FECollection<dim> &dual_fe,
+ const hp::QCollection<dim> &quadrature,
+ const hp::QCollection<dim-1> &face_quadrature,
const Function<dim> &rhs_function,
const Function<dim> &bv,
const DualFunctional::DualFunctionalBase<dim> &dual_functional)
primal_hanging_node_constraints,
dual_solution);
- DataOut<dim> data_out;
+ DataOut<dim,hp::DoFHandler<dim> > data_out;
data_out.attach_dof_handler (primal_solver.dof_handler);
data_out.add_data_vector (primal_solver.solution,
{
cell_data.fe_values.reinit (cell);
cell_data.right_hand_side
- ->value_list (cell_data.fe_values.get_quadrature_points(),
+ ->value_list (cell_data.fe_values.get_present_fe_values().get_quadrature_points(),
cell_data.rhs_values);
- cell_data.fe_values.get_function_2nd_derivatives (primal_solution,
+ cell_data.fe_values.get_present_fe_values().get_function_2nd_derivatives (primal_solution,
cell_data.cell_grad_grads);
- cell_data.fe_values.get_function_values (dual_weights,
+ cell_data.fe_values.get_present_fe_values().get_function_values (dual_weights,
cell_data.dual_weights);
double sum = 0;
- for (unsigned int p=0; p<cell_data.fe_values.n_quadrature_points; ++p)
+ for (unsigned int p=0; p<cell_data.fe_values.get_present_fe_values().n_quadrature_points; ++p)
sum += ((cell_data.rhs_values[p]+trace(cell_data.cell_grad_grads[p])) *
cell_data.dual_weights[p] *
- cell_data.fe_values.JxW (p));
+ cell_data.fe_values.get_present_fe_values().JxW (p));
error_indicators(cell_index) += sum;
}
FaceIntegrals &face_integrals) const
{
const unsigned int
- n_q_points = face_data.fe_face_values_cell.n_quadrature_points;
+ n_q_points = face_data.fe_face_values_cell.get_present_fe_values().n_quadrature_points;
face_data.fe_face_values_cell.reinit (cell, face_no);
- face_data.fe_face_values_cell.get_function_grads (primal_solution,
+ face_data.fe_face_values_cell.get_present_fe_values().get_function_grads (primal_solution,
face_data.cell_grads);
Assert (cell->neighbor(face_no).state() == IteratorState::valid,
neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
const active_cell_iterator neighbor = cell->neighbor(face_no);
face_data.fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor);
- face_data.fe_face_values_neighbor.get_function_grads (primal_solution,
+ face_data.fe_face_values_neighbor.get_present_fe_values().get_function_grads (primal_solution,
face_data.neighbor_grads);
for (unsigned int p=0; p<n_q_points; ++p)
face_data.jump_residual[p]
= ((face_data.cell_grads[p] - face_data.neighbor_grads[p]) *
- face_data.fe_face_values_cell.normal_vector(p));
+ face_data.fe_face_values_cell.get_present_fe_values().normal_vector(p));
- face_data.fe_face_values_cell.get_function_values (dual_weights,
+ face_data.fe_face_values_cell.get_present_fe_values().get_function_values (dual_weights,
face_data.dual_weights);
double face_integral = 0;
for (unsigned int p=0; p<n_q_points; ++p)
face_integral += (face_data.jump_residual[p] *
face_data.dual_weights[p] *
- face_data.fe_face_values_cell.JxW(p));
+ face_data.fe_face_values_cell.get_present_fe_values().JxW(p));
Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
ExcInternalError());
FaceIntegrals &face_integrals) const
{
const unsigned int
- n_q_points = face_data.fe_face_values_cell.n_quadrature_points;
+ n_q_points = face_data.fe_face_values_cell.get_present_fe_values().n_quadrature_points;
- const typename DoFHandler<dim>::face_iterator
+ const typename hp::DoFHandler<dim>::face_iterator
face = cell->face(face_no);
- const typename DoFHandler<dim>::cell_iterator
+ const typename hp::DoFHandler<dim>::cell_iterator
neighbor = cell->neighbor(face_no);
Assert (neighbor.state() == IteratorState::valid,
ExcInternalError());
ExcInternalError());
face_data.fe_subface_values_cell.reinit (cell, face_no, subface_no);
- face_data.fe_subface_values_cell.get_function_grads (primal_solution,
+ face_data.fe_subface_values_cell.get_present_fe_values().get_function_grads (primal_solution,
face_data.cell_grads);
face_data.fe_face_values_neighbor.reinit (neighbor_child,
neighbor_neighbor);
- face_data.fe_face_values_neighbor.get_function_grads (primal_solution,
+ face_data.fe_face_values_neighbor.get_present_fe_values().get_function_grads (primal_solution,
face_data.neighbor_grads);
for (unsigned int p=0; p<n_q_points; ++p)
face_data.jump_residual[p]
= ((face_data.neighbor_grads[p] - face_data.cell_grads[p]) *
- face_data.fe_face_values_neighbor.normal_vector(p));
+ face_data.fe_face_values_neighbor.get_present_fe_values().normal_vector(p));
- face_data.fe_face_values_neighbor.get_function_values (dual_weights,
+ face_data.fe_face_values_neighbor.get_present_fe_values().get_function_values (dual_weights,
face_data.dual_weights);
double face_integral = 0;
for (unsigned int p=0; p<n_q_points; ++p)
face_integral += (face_data.jump_residual[p] *
face_data.dual_weights[p] *
- face_data.fe_face_values_neighbor.JxW(p));
+ face_data.fe_face_values_neighbor.get_present_fe_values().JxW(p));
face_integrals[neighbor_child->face(neighbor_neighbor)]
= face_integral;
};
triangulation (Triangulation<dim>::smoothing_on_refinement);
descriptor.data->create_coarse_grid (triangulation);
- const FE_Q<dim> primal_fe(descriptor.primal_fe_degree);
- const FE_Q<dim> dual_fe(descriptor.dual_fe_degree);
- const QGauss<dim> quadrature(descriptor.dual_fe_degree+1);
- const QGauss<dim-1> face_quadrature(descriptor.dual_fe_degree+1);
+ const hp::FECollection<dim> primal_fe(FE_Q<dim>(descriptor.primal_fe_degree));
+ const hp::FECollection<dim> dual_fe(FE_Q<dim>(descriptor.dual_fe_degree));
+ const hp::QCollection<dim> quadrature(QGauss<dim>(descriptor.dual_fe_degree+1));
+ const hp::QCollection<dim-1> face_quadrature(QGauss<dim-1>(descriptor.dual_fe_degree+1));
LaplaceSolver::Base<dim> * solver = 0;
switch (descriptor.refinement_criterion)
#include <dofs/dof_constraints.h>
#include <dofs/dof_tools.h>
#include <fe/fe_q.h>
-#include <fe/fe_values.h>
+#include <fe/hp_fe_values.h>
+#include <fe/q_collection.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
#include <numerics/data_out.h>
void output_results () const;
void refine_grid ();
- static double energy (const DoFHandler<dim> &dof_handler,
+ static double energy (const hp::DoFHandler<dim> &dof_handler,
const Vector<double> &function);
Triangulation<dim> triangulation;
- FE_Q<dim> fe;
- DoFHandler<dim> dof_handler;
+ hp::FECollection<dim> fe;
+ hp::DoFHandler<dim> dof_handler;
ConstraintMatrix hanging_node_constraints;
MinimizationProblem<dim>::MinimizationProblem (const unsigned int run_number)
:
run_number (run_number),
- fe (1),
+ fe (FE_Q<dim>(1)),
dof_handler (triangulation)
{}
InitializationValues(),
present_solution);
- DoFHandler<1>::cell_iterator cell;
+ hp::DoFHandler<1>::cell_iterator cell;
cell = dof_handler.begin(0);
while (cell->at_boundary(0) == false)
cell = cell->neighbor(0);
matrix.reinit (sparsity_pattern);
residual.reinit (dof_handler.n_dofs());
- QGauss<dim> quadrature_formula(4);
- FEValues<dim> fe_values (fe, quadrature_formula,
+ hp::QCollection<dim> quadrature_formula(QGauss<dim>(4));
+ hp::FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_q_points | update_JxW_values);
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
- const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
+ const unsigned int dofs_per_cell = fe[0].dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula[0].n_quadrature_points;
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
std::vector<double> local_solution_values (n_q_points);
std::vector<Tensor<1,dim> > local_solution_grads (n_q_points);
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
fe_values.reinit (cell);
- fe_values.get_function_values (present_solution,
+ fe_values.get_present_fe_values().get_function_values (present_solution,
local_solution_values);
- fe_values.get_function_grads (present_solution,
+ fe_values.get_present_fe_values().get_function_grads (present_solution,
local_solution_grads);
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
{
const double u = local_solution_values[q_point],
- x = fe_values.quadrature_point(q_point)(0);
+ x = fe_values.get_present_fe_values().quadrature_point(q_point)(0);
const double x_minus_u3 = (x-std::pow(u,3));
const Tensor<1,dim> u_prime = local_solution_grads[q_point];
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
cell_matrix(i,j)
- += (fe_values.shape_grad(i,q_point) *
- fe_values.shape_grad(j,q_point) *
+ += (fe_values.get_present_fe_values().shape_grad(i,q_point) *
+ fe_values.get_present_fe_values().shape_grad(j,q_point) *
cell->diameter() *
cell->diameter()
+
- fe_values.shape_value(i,q_point) *
- fe_values.shape_value(j,q_point)) *
- fe_values.JxW(q_point);
+ fe_values.get_present_fe_values().shape_value(i,q_point) *
+ fe_values.get_present_fe_values().shape_value(j,q_point)) *
+ fe_values.get_present_fe_values().JxW(q_point);
for (unsigned int i=0; i<dofs_per_cell; ++i)
cell_rhs(i) += -((6. * x_minus_u3 *
gradient_power (u_prime, 4) *
- fe_values.shape_value(i,q_point)
+ fe_values.get_present_fe_values().shape_value(i,q_point)
*
(x_minus_u3 *
(u_prime *
- fe_values.shape_grad(i,q_point))
+ fe_values.get_present_fe_values().shape_grad(i,q_point))
-
(u_prime*u_prime) * u * u *
- fe_values.shape_value(i,q_point))
+ fe_values.get_present_fe_values().shape_value(i,q_point))
)
*
- fe_values.JxW(q_point));
+ fe_values.get_present_fe_values().JxW(q_point));
}
cell->get_dof_indices (local_dof_indices);
void
MinimizationProblem<dim>::output_results () const
{
- DataOut<dim> data_out;
+ DataOut<dim,hp::DoFHandler<dim> > data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (present_solution, "solution");
data_out.build_patches ();
Vector<float> error_indicators (triangulation.n_active_cells());
- QTrapez<dim> quadrature;
- FEValues<dim> fe_values (fe, quadrature,
+ QTrapez<dim> q;
+ hp::QCollection<dim> quadrature(q);
+ hp::FEValues<dim> fe_values (fe, quadrature,
update_values | update_gradients |
update_second_derivatives |
update_q_points | update_JxW_values);
- FEValues<dim> neighbor_fe_values (fe, quadrature,
+ hp::FEValues<dim> neighbor_fe_values (fe, quadrature,
update_gradients);
- std::vector<double> local_values (quadrature.n_quadrature_points);
- std::vector<Tensor<1,dim> > local_gradients (quadrature.n_quadrature_points);
- std::vector<Tensor<2,dim> > local_2nd_derivs (quadrature.n_quadrature_points);
+ std::vector<double> local_values (quadrature[0].n_quadrature_points);
+ std::vector<Tensor<1,dim> > local_gradients (quadrature[0].n_quadrature_points);
+ std::vector<Tensor<2,dim> > local_2nd_derivs (quadrature[0].n_quadrature_points);
- DoFHandler<dim>::active_cell_iterator
+ hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active (),
endc = dof_handler.end ();
for (unsigned int cell_index = 0; cell!=endc; ++cell, ++cell_index)
{
fe_values.reinit (cell);
- fe_values.get_function_values (present_solution, local_values);
- fe_values.get_function_grads (present_solution, local_gradients);
- fe_values.get_function_2nd_derivatives (present_solution, local_2nd_derivs);
+ fe_values.get_present_fe_values().get_function_values (present_solution, local_values);
+ fe_values.get_present_fe_values().get_function_grads (present_solution, local_gradients);
+ fe_values.get_present_fe_values().get_function_2nd_derivatives (present_solution, local_2nd_derivs);
double cell_residual_norm = 0;
- for (unsigned int q=0; q<quadrature.n_quadrature_points; ++q)
+ for (unsigned int q=0; q<quadrature[0].n_quadrature_points; ++q)
{
- const double x = fe_values.quadrature_point(q)[0];
+ const double x = fe_values.get_present_fe_values().quadrature_point(q)[0];
const double u = local_values[q];
const double u_prime = local_gradients[q][0];
const double u_doubleprime = local_2nd_derivs[q][0][0];
2*u_prime*(1-3*u*u*u_prime)));
cell_residual_norm += (local_residual_value * local_residual_value *
- fe_values.JxW(q));
+ fe_values.get_present_fe_values().JxW(q));
}
error_indicators(cell_index) = cell_residual_norm *
cell->diameter() * cell->diameter();
- const double x_left = fe_values.quadrature_point(0)[0];
- const double x_right = fe_values.quadrature_point(1)[0];
+ const double x_left = fe_values.get_present_fe_values().quadrature_point(0)[0];
+ const double x_right = fe_values.get_present_fe_values().quadrature_point(1)[0];
Assert (x_left == cell->vertex(0)[0], ExcInternalError());
Assert (x_right == cell->vertex(1)[0], ExcInternalError());
if (cell->at_boundary(0) == false)
{
- DoFHandler<dim>::cell_iterator left_neighbor = cell->neighbor(0);
+ hp::DoFHandler<dim>::cell_iterator left_neighbor = cell->neighbor(0);
while (left_neighbor->has_children())
left_neighbor = left_neighbor->child(1);
neighbor_fe_values.reinit (left_neighbor);
- neighbor_fe_values.get_function_grads (present_solution, local_gradients);
+ neighbor_fe_values.get_present_fe_values().get_function_grads (present_solution, local_gradients);
const double neighbor_u_prime_left = local_gradients[1][0];
if (cell->at_boundary(1) == false)
{
- DoFHandler<dim>::cell_iterator right_neighbor = cell->neighbor(1);
+ hp::DoFHandler<dim>::cell_iterator right_neighbor = cell->neighbor(1);
while (right_neighbor->has_children())
right_neighbor = right_neighbor->child(0);
neighbor_fe_values.reinit (right_neighbor);
- neighbor_fe_values.get_function_grads (present_solution, local_gradients);
+ neighbor_fe_values.get_present_fe_values().get_function_grads (present_solution, local_gradients);
const double neighbor_u_prime_right = local_gradients[0][0];
template <int dim>
double
-MinimizationProblem<dim>::energy (const DoFHandler<dim> &dof_handler,
+MinimizationProblem<dim>::energy (const hp::DoFHandler<dim> &dof_handler,
const Vector<double> &function)
{
- QGauss<dim> quadrature_formula(4);
- FEValues<dim> fe_values (dof_handler.get_fe(), quadrature_formula,
+ hp::QCollection<dim> quadrature_formula(QGauss<dim>(4));
+ hp::FEValues<dim> fe_values (dof_handler.get_fe(), quadrature_formula,
update_values | update_gradients |
update_q_points | update_JxW_values);
- const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
+ const unsigned int n_q_points = quadrature_formula[0].n_quadrature_points;
std::vector<double> local_solution_values (n_q_points);
std::vector<Tensor<1,dim> > local_solution_grads (n_q_points);
double energy = 0.;
- typename DoFHandler<dim>::active_cell_iterator
+ typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
- fe_values.get_function_values (function,
+ fe_values.get_present_fe_values().get_function_values (function,
local_solution_values);
- fe_values.get_function_grads (function,
+ fe_values.get_present_fe_values().get_function_grads (function,
local_solution_grads);
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- energy += (std::pow (fe_values.quadrature_point(q_point)(0)
+ energy += (std::pow (fe_values.get_present_fe_values().quadrature_point(q_point)(0)
-
std::pow (local_solution_values[q_point], 3),
2) *
gradient_power (local_solution_grads[q_point], 6) *
- fe_values.JxW (q_point));
+ fe_values.get_present_fe_values().JxW (q_point));
}
return energy;