From 87b82457f0aef0f1bd4c3b9b49a34462c08bc66c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 26 Sep 2014 10:21:10 -0500 Subject: [PATCH] Allow Q(p) and DGQ(r) elements to touch each other in hp::DoFHandlers. The thing is that the FE_DGQ is discontinuous anyway, so there can't be any constraints to begin with. --- doc/news/changes.h | 9 ++++ source/fe/fe_dgq.cc | 60 ++++++++------------- source/fe/fe_q_base.cc | 33 ++++++++++++ tests/hp/q_vs_dgq.cc | 112 +++++++++++++++++++++++++++++++++++++++ tests/hp/q_vs_dgq.output | 11 ++++ 5 files changed, 187 insertions(+), 38 deletions(-) create mode 100644 tests/hp/q_vs_dgq.cc create mode 100644 tests/hp/q_vs_dgq.output diff --git a/doc/news/changes.h b/doc/news/changes.h index b3829b40ab..d7d5a6b361 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -295,6 +295,15 @@ inconvenience this causes.

Specific improvements

    +
  1. 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. +
    + (Wolfgang Bangerth, 2014/09/26) +
    2. +
  2. 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. diff --git a/source/fe/fe_dgq.cc b/source/fe/fe_dgq.cc index a61ecfb8ae..b00486470a 100644 --- a/source/fe/fe_dgq.cc +++ b/source/fe/fe_dgq.cc @@ -617,16 +617,12 @@ std::vector > FE_DGQ:: hp_vertex_dof_identities (const FiniteElement &fe_other) const { - // there are no such constraints for DGQ - // elements at all - if (dynamic_cast*>(&fe_other) != 0) - return - std::vector > (); - else - { - Assert (false, ExcNotImplemented()); - return std::vector > (); - } + // 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 > (); } @@ -636,16 +632,12 @@ std::vector > FE_DGQ:: hp_line_dof_identities (const FiniteElement &fe_other) const { - // there are no such constraints for DGQ - // elements at all - if (dynamic_cast*>(&fe_other) != 0) - return - std::vector > (); - else - { - Assert (false, ExcNotImplemented()); - return std::vector > (); - } + // 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 > (); } @@ -655,16 +647,12 @@ std::vector > FE_DGQ:: hp_quad_dof_identities (const FiniteElement &fe_other) const { - // there are no such constraints for DGQ - // elements at all - if (dynamic_cast*>(&fe_other) != 0) - return - std::vector > (); - else - { - Assert (false, ExcNotImplemented()); - return std::vector > (); - } + // 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 > (); } @@ -673,14 +661,10 @@ template FiniteElementDomination::Domination FE_DGQ::compare_for_face_domination (const FiniteElement &fe_other) const { - // check whether both are discontinuous - // elements, see the description of - // FiniteElementDomination::Domination - if (dynamic_cast*>(&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; } diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index 99c65476d8..b384c1a3a1 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -683,6 +683,17 @@ hp_vertex_dof_identities (const FiniteElement &fe_other) const // equivalencies to be recorded return std::vector > (); } + 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 > (); + } else { Assert (false, ExcNotImplemented()); @@ -734,6 +745,17 @@ hp_line_dof_identities (const FiniteElement &fe_other) const // equivalencies to be recorded return std::vector > (); } + 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 > (); + } else { Assert (false, ExcNotImplemented()); @@ -789,6 +811,17 @@ hp_quad_dof_identities (const FiniteElement &fe_other) cons // equivalencies to be recorded return std::vector > (); } + 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 > (); + } else { Assert (false, ExcNotImplemented()); diff --git a/tests/hp/q_vs_dgq.cc b/tests/hp/q_vs_dgq.cc new file mode 100644 index 0000000000..abfdf5a161 --- /dev/null +++ b/tests/hp/q_vs_dgq.cc @@ -0,0 +1,112 @@ +// --------------------------------------------------------------------- +// +// 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 + +#include +#include +#include +#include +#include + + +namespace Step27 +{ + using namespace dealii; + + template + class MixedFECollection + { + public: + MixedFECollection (); + ~MixedFECollection (); + + void run (); + + private: + + Triangulation triangulation; + hp::DoFHandler dof_handler; + hp::FECollection fe_collection; + + }; + + template + MixedFECollection::MixedFECollection () + : + dof_handler (triangulation) + { + + } + + template + MixedFECollection::~MixedFECollection () + { + dof_handler.clear (); + } + + template + void MixedFECollection::run () + { + // add two a CG and a DG finite element object to fe_collection + fe_collection.push_back (FE_Q(1)); + fe_collection.push_back (FE_DGQ(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::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; +} diff --git a/tests/hp/q_vs_dgq.output b/tests/hp/q_vs_dgq.output new file mode 100644 index 0000000000..16ded5a89e --- /dev/null +++ b/tests/hp/q_vs_dgq.output @@ -0,0 +1,11 @@ + +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 -- 2.39.5