quad.push_back (quad_temp);
quadrature->add_quadrature (*quad_temp);
- //TODO: This is currently a design error in the hp Framework.
- // While the idea of using a QuadratureCollection for the
- // element stiffness matrices is straightforward, it is
- // not complete for the face terms. At these places, we can
- // have different polynomial degrees on the two faces. Therefore
- // a quadrature rule, which is sufficient for the highest polynomial
- // degree on both faces, would be required. This will certainly need
- // some changes in the hpFEValues classes. But for a first example,
- // simply using a uniformly high quadrature rule will do the job
- // although is obviously not very efficient, as even linear polynomials
- // will be integrated with 7th order Gauss formulas.
- face_quad_temp = new QGauss<dim-1> (7); // (i+2);
+ // face_quad_temp = new QGauss<dim-1> (7);
+ face_quad_temp = new QGauss<dim-1> (i+2);
face_quad.push_back (face_quad_temp);
face_quadrature->add_quadrature (*face_quad_temp);
}
neighbor->child(GeometryInfo<dim>::
child_cell_on_face(neighbor2,subface_no));
+ // an additional speciality
+ // for the hp method appears
+ // on the faces. To get an
+ // efficient assembly, the
+ // lowest order but
+ // sufficient quadrature
+ // rule should be used. Hence
+ // the face quadrature rule of the
+ // higher order element
+ // will be used.
+ const unsigned int use_fe_index =
+ neighbor_child->active_fe_index () > cell->active_fe_index () ?
+ neighbor_child->active_fe_index () : cell->active_fe_index ();
+
+
// As these are
// quite
// complicated
// of the
// neighboring
// child cell.
- fe_v_subface_x.reinit (cell, face_no, subface_no);
- fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2);
+ fe_v_subface_x.reinit (cell, face_no, subface_no, use_fe_index);
+ fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2, use_fe_index);
dg.assemble_face_term1(fe_v_subface_x.get_present_fe_values (),
fe_v_face_neighbor_x.get_present_fe_values (),
// cell.
const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
+ // Like before. Use
+ // quadrature rule
+ // of higher order
+ // cell.
+ const unsigned int use_fe_index =
+ neighbor->active_fe_index () > cell->active_fe_index () ?
+ neighbor->active_fe_index () : cell->active_fe_index ();
+
// We reinit
// the
// ``FEFaceValues''
// the
// corresponding
// face terms.
- fe_v_face_x.reinit (cell, face_no);
- fe_v_face_neighbor_x.reinit (neighbor, neighbor2);
+ fe_v_face_x.reinit (cell, face_no, use_fe_index);
+ fe_v_face_neighbor_x.reinit (neighbor, neighbor2, use_fe_index);
dg.assemble_face_term1(fe_v_face_x.get_present_fe_values (),
fe_v_face_neighbor_x.get_present_fe_values (),
->child(GeometryInfo<dim>::child_cell_on_face(
face_no,neighbor_subface_no)) == cell, ExcInternalError());
+
+ // Like before. Use
+ // quadrature rule
+ // of higher order
+ // cell.
+ const unsigned int use_fe_index =
+ neighbor->active_fe_index () > cell->active_fe_index () ?
+ neighbor->active_fe_index () : cell->active_fe_index ();
+
// Reinit the
// appropriate
// ``FEFaceValues''
// and assemble
// the face
// terms.
- fe_v_face_x.reinit (cell, face_no);
+ fe_v_face_x.reinit (cell, face_no, use_fe_index);
fe_v_subface_neighbor_x.reinit (neighbor, neighbor_face_no,
- neighbor_subface_no);
+ neighbor_subface_no, use_fe_index);
dg.assemble_face_term1(fe_v_face_x.get_present_fe_values (),
fe_v_subface_neighbor_x.get_present_fe_values (),
neighbor2,subface_no));
const unsigned int dofs_on_neighbor_child = neighbor_child->get_fe().dofs_per_cell;
+ const unsigned int use_fe_index =
+ neighbor_child->active_fe_index () > cell->active_fe_index () ?
+ neighbor_child->active_fe_index () : cell->active_fe_index ();
+
Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
ExcInternalError());
Assert (!neighbor_child->has_children(), ExcInternalError());
u_vn_matrix = 0;
un_vn_matrix = 0;
- fe_v_subface_x.reinit (cell, face_no, subface_no);
- fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2);
+ fe_v_subface_x.reinit (cell, face_no, subface_no, use_fe_index);
+ fe_v_face_neighbor_x.reinit (neighbor_child, neighbor2, use_fe_index);
dg.assemble_face_term2(fe_v_subface_x.get_present_fe_values (),
fe_v_face_neighbor_x.get_present_fe_values (),
{
const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
+ const unsigned int use_fe_index =
+ neighbor->active_fe_index () > cell->active_fe_index () ?
+ neighbor->active_fe_index () : cell->active_fe_index ();
+
un_v_matrix = 0;
u_vn_matrix = 0;
un_vn_matrix = 0;
- fe_v_face_x.reinit (cell, face_no);
- fe_v_face_neighbor_x.reinit (neighbor, neighbor2);
+ fe_v_face_x.reinit (cell, face_no, use_fe_index);
+ fe_v_face_neighbor_x.reinit (neighbor, neighbor2, use_fe_index);
dg.assemble_face_term2(fe_v_face_x.get_present_fe_values (),
fe_v_face_neighbor_x.get_present_fe_values (),
{
try
{
- DGMethod<2> dgmethod;
+ DGMethod<3> dgmethod;
dgmethod.run ();
}
catch (std::exception &exc)