The thing is that the FE_DGQ is discontinuous anyway, so there can't be any constraints to begin with.
<h3>Specific improvements</h3>
<ol>
+ <li> Fixed: Trying to have FE_Q(p) and FE_DGQ(r) elements next to each
+ other in an hp::DoFHandler object led to assertions saying that these two
+ elements don't know how to compute interface constraints where such
+ elements touch. This has now been fixed: FE_DGQ is a discontinuous element,
+ so there cannot be any interface constraints at all.
+ <br>
+ (Wolfgang Bangerth, 2014/09/26)
+ </li>
+
<li> Fixed: The function TrilinosWrappers::VectorBase::sadd(double factor,
VectorBase &v) erroneously added factor*v instead of scaling the calling
vector by factor. This is now fixed.
FE_DGQ<dim, spacedim>::
hp_vertex_dof_identities (const FiniteElement<dim, spacedim> &fe_other) const
{
- // there are no such constraints for DGQ
- // elements at all
- if (dynamic_cast<const FE_DGQ<dim, spacedim>*>(&fe_other) != 0)
- return
- std::vector<std::pair<unsigned int, unsigned int> > ();
- else
- {
- Assert (false, ExcNotImplemented());
- return std::vector<std::pair<unsigned int, unsigned int> > ();
- }
+ // this element is discontinuous, so by definition there can
+ // be no identities between its dofs and those of any neighbor
+ // (of whichever type the neighbor may be -- after all, we have
+ // no face dofs on this side to begin with)
+ return
+ std::vector<std::pair<unsigned int, unsigned int> > ();
}
FE_DGQ<dim, spacedim>::
hp_line_dof_identities (const FiniteElement<dim, spacedim> &fe_other) const
{
- // there are no such constraints for DGQ
- // elements at all
- if (dynamic_cast<const FE_DGQ<dim, spacedim>*>(&fe_other) != 0)
- return
- std::vector<std::pair<unsigned int, unsigned int> > ();
- else
- {
- Assert (false, ExcNotImplemented());
- return std::vector<std::pair<unsigned int, unsigned int> > ();
- }
+ // this element is discontinuous, so by definition there can
+ // be no identities between its dofs and those of any neighbor
+ // (of whichever type the neighbor may be -- after all, we have
+ // no face dofs on this side to begin with)
+ return
+ std::vector<std::pair<unsigned int, unsigned int> > ();
}
FE_DGQ<dim, spacedim>::
hp_quad_dof_identities (const FiniteElement<dim, spacedim> &fe_other) const
{
- // there are no such constraints for DGQ
- // elements at all
- if (dynamic_cast<const FE_DGQ<dim, spacedim>*>(&fe_other) != 0)
- return
- std::vector<std::pair<unsigned int, unsigned int> > ();
- else
- {
- Assert (false, ExcNotImplemented());
- return std::vector<std::pair<unsigned int, unsigned int> > ();
- }
+ // this element is discontinuous, so by definition there can
+ // be no identities between its dofs and those of any neighbor
+ // (of whichever type the neighbor may be -- after all, we have
+ // no face dofs on this side to begin with)
+ return
+ std::vector<std::pair<unsigned int, unsigned int> > ();
}
FiniteElementDomination::Domination
FE_DGQ<dim, spacedim>::compare_for_face_domination (const FiniteElement<dim, spacedim> &fe_other) const
{
- // check whether both are discontinuous
- // elements, see the description of
- // FiniteElementDomination::Domination
- if (dynamic_cast<const FE_DGQ<dim, spacedim>*>(&fe_other) != 0)
- return FiniteElementDomination::no_requirements;
-
- Assert (false, ExcNotImplemented());
- return FiniteElementDomination::neither_element_dominates;
+ // this is a discontinuous element, so by definition there will
+ // be no constraints wherever this element comes together with
+ // any other kind of element
+ return FiniteElementDomination::no_requirements;
}
// equivalencies to be recorded
return std::vector<std::pair<unsigned int, unsigned int> > ();
}
+ else if (fe_other.dofs_per_face == 0)
+ {
+ // if the other element has no elements on faces at all,
+ // then it would be impossible to enforce any kind of
+ // continuity even if we knew exactly what kind of element
+ // we have -- simply because the other element declares
+ // that it is discontinuous because it has no DoFs on
+ // its faces. in that case, just state that we have no
+ // constraints to declare
+ return std::vector<std::pair<unsigned int, unsigned int> > ();
+ }
else
{
Assert (false, ExcNotImplemented());
// equivalencies to be recorded
return std::vector<std::pair<unsigned int, unsigned int> > ();
}
+ else if (fe_other.dofs_per_face == 0)
+ {
+ // if the other element has no elements on faces at all,
+ // then it would be impossible to enforce any kind of
+ // continuity even if we knew exactly what kind of element
+ // we have -- simply because the other element declares
+ // that it is discontinuous because it has no DoFs on
+ // its faces. in that case, just state that we have no
+ // constraints to declare
+ return std::vector<std::pair<unsigned int, unsigned int> > ();
+ }
else
{
Assert (false, ExcNotImplemented());
// equivalencies to be recorded
return std::vector<std::pair<unsigned int, unsigned int> > ();
}
+ else if (fe_other.dofs_per_face == 0)
+ {
+ // if the other element has no elements on faces at all,
+ // then it would be impossible to enforce any kind of
+ // continuity even if we knew exactly what kind of element
+ // we have -- simply because the other element declares
+ // that it is discontinuous because it has no DoFs on
+ // its faces. in that case, just state that we have no
+ // constraints to declare
+ return std::vector<std::pair<unsigned int, unsigned int> > ();
+ }
else
{
Assert (false, ExcNotImplemented());
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test by Korosh Taebi: check that we can gave Q(p) and DGQ(r) in the same
+// mesh
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_dgq.h>
+
+
+namespace Step27
+{
+ using namespace dealii;
+
+ template <int dim>
+ class MixedFECollection
+ {
+ public:
+ MixedFECollection ();
+ ~MixedFECollection ();
+
+ void run ();
+
+ private:
+
+ Triangulation<dim> triangulation;
+ hp::DoFHandler<dim> dof_handler;
+ hp::FECollection<dim> fe_collection;
+
+ };
+
+ template <int dim>
+ MixedFECollection<dim>::MixedFECollection ()
+ :
+ dof_handler (triangulation)
+ {
+
+ }
+
+ template <int dim>
+ MixedFECollection<dim>::~MixedFECollection ()
+ {
+ dof_handler.clear ();
+ }
+
+ template <int dim>
+ void MixedFECollection<dim>::run ()
+ {
+ // add two a CG and a DG finite element object to fe_collection
+ fe_collection.push_back (FE_Q<dim>(1));
+ fe_collection.push_back (FE_DGQ<dim>(1));
+ deallog << " fe_collection size = " << fe_collection.size() << std::endl;
+
+ // produce a simple grid with 4 cells
+ GridGenerator::hyper_cube (triangulation, 0, 1);
+ triangulation.refine_global (1);
+
+ // looping over all cells and assigning the FE_DG object to the first cell
+ typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (unsigned int counter = 0; cell!=endc; ++cell, counter ++)
+ if (counter == 0)
+ {
+ cell->set_active_fe_index (cell->active_fe_index() + 1);
+ }
+
+ dof_handler.distribute_dofs (fe_collection);
+
+ deallog << " Number of active cells: " << triangulation.n_active_cells() << std::endl
+ << " Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl;
+ }
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("output");
+ logfile.precision(2);
+
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ Step27::MixedFECollection<1>().run();
+ Step27::MixedFECollection<2>().run();
+ Step27::MixedFECollection<3>().run();
+
+ deallog << "OK" << std::endl;
+}
--- /dev/null
+
+DEAL:: fe_collection size = 2
+DEAL:: Number of active cells: 2
+DEAL:: Number of degrees of freedom: 4
+DEAL:: fe_collection size = 2
+DEAL:: Number of active cells: 4
+DEAL:: Number of degrees of freedom: 12
+DEAL:: fe_collection size = 2
+DEAL:: Number of active cells: 8
+DEAL:: Number of degrees of freedom: 34
+DEAL::OK
projects / dealii.git / commitdiff