From 6ecaae0e619635842376ac79e4ea322f1742004d Mon Sep 17 00:00:00 2001 From: kanschat Date: Wed, 27 Oct 2010 00:11:18 +0000 Subject: [PATCH] add constraints in 3D... hopefully git-svn-id: https://svn.dealii.org/trunk@22512 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/source/fe/fe_nedelec.cc | 169 ++++++++++++++++++++++++++++++-- 1 file changed, 160 insertions(+), 9 deletions(-) diff --git a/deal.II/source/fe/fe_nedelec.cc b/deal.II/source/fe/fe_nedelec.cc index 0b58fd2d8a..5806d1af73 100644 --- a/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/source/fe/fe_nedelec.cc @@ -80,16 +80,20 @@ deg (p) 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); + + // In two space dimensions, things are easy, since the child + // faces do not share any degrees of freedom + Assert (dim<4, ExcNotImplemented()); 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) @@ -102,8 +106,155 @@ deg (p) ++target_row; } } + // Now the 2D case, where the children share degrees of freedom on + // edges. which means, we must go though all edges to reproduce the numbering. + // Lets pray that at least the subfaces are numbered as expected else - this->interface_constraints.reinit(0,0); + { + unsigned int target_row = 0; + // right edge of lower left child + for (unsigned int jj = 0; jj < this->dofs_per_line; ++jj) + { + const unsigned int j = this->fist_face_line_index + this->dofs_per_line + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[0] (j, k); + ++target_row; + } + // right edge of upper left child + for (unsigned int jj = 0; jj < this->dofs_per_line; ++jj) + { + const unsigned int j = this->fist_face_line_index + this->dofs_per_line + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[2] (j, k); + ++target_row; + } + // upper edge of lower left child + for (unsigned int jj = 0; jj < this->dofs_per_line; ++jj) + { + const unsigned int j = this->fist_face_line_index + 3*this->dofs_per_line + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[0] (j, k); + ++target_row; + } + // upper edge of lower right child + for (unsigned int jj = 0; jj < this->dofs_per_line; ++jj) + { + const unsigned int j = this->fist_face_line_index + 3*this->dofs_per_line + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[1] (j, k); + ++target_row; + } + // left edge of lower left child + for (unsigned int jj = 0; jj < this->dofs_per_line; ++jj) + { + const unsigned int j = this->fist_face_line_index + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[0] (j, k); + ++target_row; + } + // left edge of upper left child + for (unsigned int jj = 0; jj < this->dofs_per_line; ++jj) + { + const unsigned int j = this->fist_face_line_index + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[2] (j, k); + ++target_row; + } + // right edge of lower right child + for (unsigned int jj = 0; jj < this->dofs_per_line; ++jj) + { + const unsigned int j = this->fist_face_line_index + this->dofs_per_line + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[1] (j, k); + ++target_row; + } + // right edge of upper right child + for (unsigned int jj = 0; jj < this->dofs_per_line; ++jj) + { + const unsigned int j = this->fist_face_line_index + this->dofs_per_line + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[3] (j, k); + ++target_row; + } + // lower edge of lower left child + for (unsigned int jj = 0; jj < this->dofs_per_line; ++jj) + { + const unsigned int j = this->fist_face_line_index + 2*this->dofs_per_line + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[0] (j, k); + ++target_row; + } + // lower edge of lower right child + for (unsigned int jj = 0; jj < this->dofs_per_line; ++jj) + { + const unsigned int j = this->fist_face_line_index + 2*this->dofs_per_line + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[1] (j, k); + ++target_row; + } + // upper edge of upper left child + for (unsigned int jj = 0; jj < this->dofs_per_line; ++jj) + { + const unsigned int j = this->fist_face_line_index + 3*this->dofs_per_line + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[2] (j, k); + ++target_row; + } + // upper edge of upper right child + for (unsigned int jj = 0; jj < this->dofs_per_line; ++jj) + { + const unsigned int j = this->fist_face_line_index + 3*this->dofs_per_line + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[3] (j, k); + ++target_row; + } + // Now the interior dofs of each subface + for (unsigned int jj = 0; jj < this->dofs_per_quad; ++jj) + { + const unsigned int j = this->fist_face_quad_index + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[0] (j, k); + ++target_row; + } + for (unsigned int jj = 0; jj < this->dofs_per_quad; ++jj) + { + const unsigned int j = this->fist_face_quad_index + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[1] (j, k); + ++target_row; + } + for (unsigned int jj = 0; jj < this->dofs_per_quad; ++jj) + { + const unsigned int j = this->fist_face_quad_index + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[2] (j, k); + ++target_row; + } + for (unsigned int jj = 0; jj < this->dofs_per_quad; ++jj) + { + const unsigned int j = this->fist_face_quad_index + jj; + for (unsigned int k = 0; k < face_embeddings[0].n (); ++k) + this->interface_constraints (target_row, k) + = face_embeddings[3] (j, k); + ++target_row; + } + Assert(target_row == this->interface_constraints.m(), ExcInternalError()); + } } -- 2.39.5