From 5912d53fa57b7f54b4b0f7b2f56ae0923dd33db5 Mon Sep 17 00:00:00 2001 From: hartmann Date: Mon, 27 Jun 2005 07:24:59 +0000 Subject: [PATCH] Simplify implementation of FE_Q<3>::initialize_constraints by using the QProjector::project_to_subface function. git-svn-id: https://svn.dealii.org/trunk@10946 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/fe_q.cc | 148 ++++++++++++++---------------- 1 file changed, 71 insertions(+), 77 deletions(-) diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 146bebb623..45a67ad71c 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -11,6 +11,8 @@ // //--------------------------------------------------------------------------- +#include +#include #include #include @@ -950,86 +952,78 @@ FE_Q<3>::initialize_constraints () // Generate destination points. std::vector > constraint_points; - const std::vector > &un_supp - = this->get_unit_face_support_points (); - const unsigned int pnts = un_supp.size (); // Add midpoint constraint_points.push_back (Point (0.5, 0.5)); - // Add midpoints of lines of "mother-face" - for (unsigned int face = 0; - face < GeometryInfo::subfaces_per_face; ++face) - { - Point pnt_temp = un_supp[(face + 1) % 4]; - pnt_temp *= 0.5; - pnt_temp += (GeometryInfo::unit_cell_vertex (face) * 0.5); - constraint_points.push_back (pnt_temp); - } - - // Add nodes of lines interior in the + // Add midpoints of lines of // "mother-face" - for (unsigned int face = 0; - face < GeometryInfo::subfaces_per_face; ++face) - { - const unsigned int line_offset - = 4 + ((face + 1) % 4) * (this->degree-1); - for (unsigned int line_dof = 0; line_dof < this->degree-1; ++line_dof) - { - Point pnt_temp = un_supp[line_offset + line_dof]; - pnt_temp *= 0.5; - pnt_temp += (GeometryInfo::unit_cell_vertex (face) * 0.5); - constraint_points.push_back (pnt_temp); - } - } + constraint_points.push_back (Point (0.5, 0)); + constraint_points.push_back (Point (1, 0.5)); + constraint_points.push_back (Point (0.5, 1)); + constraint_points.push_back (Point (0, 0.5)); - // DoFs on bordering lines - for (unsigned int line = 0; - line < GeometryInfo::lines_per_face; ++line) + if (this->degree>1) { - // This face index runs through the two - // faces, which share the line "line" - // with the coarse element. - for (unsigned int face = 0; face < 2; ++face) - { - const unsigned int line_offset = 4 + (line * (this->degree-1)); - - // Line 2 and 3 have a different - // ordering - const unsigned int offset - = ((line < 2) ? - ((line + face) % 4) : - ((line + 1 - face) % 4)); - - for (unsigned int line_dof = 0; - line_dof < this->degree-1; ++line_dof) - { - Point pnt_temp = un_supp[line_offset + line_dof]; - pnt_temp *= 0.5; - pnt_temp += (GeometryInfo::unit_cell_vertex (offset) * 0.5); - constraint_points.push_back (pnt_temp); - } - } - } - - // Create constraints for interior nodes - const unsigned int dofs_per_face = (this->degree-1) * (this->degree-1); - for (unsigned int face = 0; - face < GeometryInfo::subfaces_per_face; ++face) - { - const unsigned int face_offset = 4 + (4 * (this->degree-1)); - for (unsigned int face_dof = 0; face_dof < dofs_per_face; ++face_dof) - { - Point pnt_temp = un_supp[face_offset + face_dof]; - pnt_temp *= 0.5; - pnt_temp += (GeometryInfo::unit_cell_vertex (face) * 0.5); - constraint_points.push_back (pnt_temp); - } + const unsigned int n=this->degree-1; + const double step=1./this->degree; + vector > line_support_points(n); + for (unsigned int i=0; i qline(line_support_points); + + // auxiliary points in 2d + vector > p_line(n); + + // Add nodes of lines interior + // in the "mother-face" + + // line 5: use line 15 + QProjector::project_to_subface(qline, 3, 0, p_line); + for (unsigned int i=0; i (0.5, 0)); + // line 6: use line 10 + QProjector::project_to_subface(qline, 0, 1, p_line); + for (unsigned int i=0; i (0, 0.5)); + // line 7: use line 16 + QProjector::project_to_subface(qline, 3, 1, p_line); + for (unsigned int i=0; i (0.5, 0)); + // line 8: use line 9 + QProjector::project_to_subface(qline, 0, 0, p_line); + for (unsigned int i=0; i (0, 0.5)); + + // DoFs on bordering lines + // lines 9-16 + for (unsigned int face=0; face::faces_per_cell; ++face) + for (unsigned int subface=0; + subface::subfaces_per_face; ++subface) + { + QProjector::project_to_subface(qline, face, subface, p_line); + constraint_points.insert(constraint_points.end(), + p_line.begin(), p_line.end()); + } + + // Create constraints for + // interior nodes + vector > inner_points(n*n); + for (unsigned int i=0, iy=1; iy<=n; ++iy) + for (unsigned int ix=1; ix<=n; ++ix) + inner_points[i++] = Point (ix*step, iy*step); + + for (unsigned int child=0; + child::children_per_cell; ++child) + for (unsigned int i=0; i::child_to_cell_coordinates(inner_points[i], child)); } // Now construct relation between // destination (child) and source (mother) // dofs. + const unsigned int pnts=(this->degree+1)*(this->degree+1); std::vector v; for (unsigned int i=0;i<=this->degree;++i) v.push_back(Polynomials::LagrangeEquidistant(this->degree,i)); @@ -1039,7 +1033,7 @@ FE_Q<3>::initialize_constraints () this->interface_constraints .TableBase<2,double>::reinit (this->interface_constraints_size()); - for (unsigned int j = 0; j < constraint_points.size(); ++j) + for (unsigned int i=0; idegree * 2); bool mirror[dim - 1]; @@ -1062,10 +1056,10 @@ FE_Q<3>::initialize_constraints () // will be scaled back. The equal // treatment of all coordinates should // eliminate any FP errors. - for (unsigned int k = 0; k < dim - 1; ++k) + for (unsigned int k=0; k (constraint_points[j](k) * interval + 0.25); + static_cast (constraint_points[i](k) * interval + 0.25); constraint_point(k) = 1.*coord_int / interval; // The following lines of code @@ -1112,11 +1106,11 @@ FE_Q<3>::initialize_constraints () constraint_point(k) = 1.0 - constraint_point(k); } - for (unsigned i = 0; i < pnts; ++i) + for (unsigned j=0; jdegree + 1), - face_index_map[i] / (this->degree + 1) }; + = { face_index_map[j] % (this->degree + 1), + face_index_map[j] / (this->degree + 1) }; for (unsigned int k = 0; k<2; ++k) if (mirror[k]) @@ -1125,7 +1119,7 @@ FE_Q<3>::initialize_constraints () const unsigned int new_index = indices[1] * (this->degree + 1) + indices[0]; - this->interface_constraints(j,i) = + this->interface_constraints(i,j) = face_polynomials.compute_value (new_index, constraint_point); // if the value is small up @@ -1138,8 +1132,8 @@ FE_Q<3>::initialize_constraints () // of other DoFs a // constrained DoF would // couple to) - if (std::fabs(this->interface_constraints(j,i)) < 1e-14) - this->interface_constraints(j,i) = 0; + if (std::fabs(this->interface_constraints(i,j)) < 1e-14) + this->interface_constraints(i,j) = 0; } } } -- 2.39.5