From 2531320bc18909618b25c14a2266631c42d31831 Mon Sep 17 00:00:00 2001 From: kanschat Date: Sun, 24 Oct 2010 18:27:59 +0000 Subject: [PATCH] fix interface constraints in 2D and disable in 3D git-svn-id: https://svn.dealii.org/trunk@22460 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_nedelec.h | 10 ++++ deal.II/deal.II/source/fe/fe_nedelec.cc | 66 ++++++++++++++++--------- 2 files changed, 53 insertions(+), 23 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_nedelec.h b/deal.II/deal.II/include/fe/fe_nedelec.h index ffcf7adbf9..db49d30724 100644 --- a/deal.II/deal.II/include/fe/fe_nedelec.h +++ b/deal.II/deal.II/include/fe/fe_nedelec.h @@ -33,6 +33,11 @@ template class MappingQ; /*@{*/ /** + * @warning Several aspects of the implementation are + * experimental. For the moment, it is safe to use the element on + * globally refined meshes with consistent orientation of faces. See + * the todo entries below for more detailed caveats. + * * Implementation of Nédélec elements, conforming with the * space Hcurl. These elements generate vector fields with * tangential components continuous between mesh cells. @@ -59,6 +64,11 @@ template class MappingQ; * This class is not implemented for the codimension one case * (spacedim != dim). * + * @todo The constraint matrices for hanging nodes are only + * implemented in 2D. Currently, the 3D version will run without an + * exception being triggered, but results at refinement edges will be + * wrong. + * * @todo Even if this element is implemented for two and three space * dimensions, the definition of the node values relies on * consistently oriented faces in 3D. Therefore, care should be taken diff --git a/deal.II/deal.II/source/fe/fe_nedelec.cc b/deal.II/deal.II/source/fe/fe_nedelec.cc index 802dc6a05f..0b58fd2d8a 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec.cc @@ -62,28 +62,48 @@ deg (p) // Fill prolongation matrices with embedding operators FETools::compute_embedding_matrices (*this, this->prolongation); initialize_restriction (); + + switch (dim) + { + case 2: + this->interface_constraints.reinit ( + GeometryInfo::max_children_per_face * this->dofs_per_face, + this->dofs_per_face); + break; + case 3: + this->interface_constraints.reinit ( + GeometryInfo::max_children_per_face * this->dofs_per_face + - 4 * this->dofs_per_line, + this->dofs_per_face); + break; + default: + Assert(false, ExcNotImplemented()); + break; + } - FullMatrix face_embeddings[GeometryInfo::max_children_per_face]; - - for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) - face_embeddings[i].reinit (this->dofs_per_face, this->dofs_per_face); - - FETools::compute_face_embedding_matrices - (*this, face_embeddings, 0, 0); - this->interface_constraints.reinit ((1 << (dim - 1)) * this->dofs_per_face, - this->dofs_per_face); - - unsigned int target_row = 0; - - for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) - for (unsigned int j = 0; j < face_embeddings[i].m (); ++j) - { - for (unsigned int k = 0; k < face_embeddings[i].n (); ++k) - this->interface_constraints (target_row, k) - = face_embeddings[i] (j, k); - - ++target_row; - } + if (dim==2) + { + FullMatrix face_embeddings[GeometryInfo::max_children_per_face]; + + for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) + face_embeddings[i].reinit (this->dofs_per_face, this->dofs_per_face); + + FETools::compute_face_embedding_matrices + (*this, face_embeddings, 0, 0); + unsigned int target_row = 0; + + for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) + for (unsigned int j = 0; j < face_embeddings[i].m (); ++j) + { + for (unsigned int k = 0; k < face_embeddings[i].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[i] (j, k); + + ++target_row; + } + } + else + this->interface_constraints.reinit(0,0); } @@ -100,7 +120,7 @@ FE_Nedelec::get_name () const // have to be kept in synch std::ostringstream namebuf; - namebuf << "FE_Nedelec<" << dim << ">(" << deg << ")"; + namebuf << "FE_Nedelec<" << dim << ">(" << this->tensor_degree()-1 << ")"; return namebuf.str(); } @@ -1304,7 +1324,7 @@ template bool FE_Nedelec::hp_constraints_are_implemented () const { - return true; + return dim != 2; } template -- 2.39.5