From: hartmann Date: Mon, 27 Jun 2005 07:34:54 +0000 (+0000) Subject: Total reimplemenation of the FE_Q<2>::initialize_constraints function to be in line... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5c16d7c29da21e4f35973aa72dfcc741ebe3a399;p=dealii-svn.git Total reimplemenation of the FE_Q<2>::initialize_constraints function to be in line with the implementation of the FE_Q<3>::initialize_constraints function. Doc updates. git-svn-id: https://svn.dealii.org/trunk@10947 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 45a67ad71c..1a115b3177 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -596,44 +596,23 @@ FE_Q<2>::initialize_constraints () // traces of the shape functions is // an element of P_{k} (in 2d), or // Q_{k} (in 3d), where k is the - // degree of the element - // - // from this, we interpolate - // between mother and cell - // face. for the general case, this - // may be a little complicated if - // we don't use Lagrange - // interpolation polynomials, since - // then we can't just use point - // interpolation. what we do - // instead is to evaluate at a - // number of points and then invert - // the interpolation matrix. here, - // for the FE_Q elements, we - // actually do have Lagrange - // polynomials, but we still follow - // the general scheme since this - // code here is the master copy for - // what we use in other elements as - // well. however, there are places - // where we make use of the fact - // that we have Lagrange - // interpolation polynomials. - - // mathematically speaking, the - // interpolation process works in - // the following way: on each - // subface, we want that finite - // element solututions from both - // sides coincide. i.e. if a and b - // are expansion coefficients for - // the shape functions from both - // sides, we seek a relation - // between x and y such that - // sum_i a_i phi^c_i(x) - // == sum_j b_j phi_j(x) + // degree of the element. from + // this, we interpolate between + // mother and cell face. + + // the interpolation process works + // as followings: on each subface, + // we want that finite element + // solutions from both sides + // coincide. i.e. if a and b are + // expansion coefficients for the + // shape functions from both sides, + // we seek a relation between x and + // y such that + // sum_j a_j phi^c_j(x) + // == sum_j b_j phi_j(x) // for all points x on the - // interface. here, phi^c_i are the + // interface. here, phi^c_j are the // shape functions on the small // cell on one side of the face, // and phi_j those on the big cell @@ -645,229 +624,94 @@ FE_Q<2>::initialize_constraints () // need n evaluation points, and we // choose them equidistantly. // - // what one then gets is a matrix - // system - // a A == b B + // we obtain the matrix system + // A a == B b // where - // A_ij = phi^c_i(x_j) - // B_ij = phi_i(x_j) - // and the relation we are looking for - // is - // a = (A^T)^-1 B^T b + // A_ij = phi^c_j(x_i) + // B_ij = phi_j(x_i) + // and the relation we are looking + // for is + // a = A^-1 B b // - // below, we build up these - // matrices, but rather than - // transposing them after the - // fact, we do so while building - // them. A will be - // subface_interpolation, B will be - // face_interpolation. note that we - // build up these matrices for all - // faces at once, rather than - // considering them separately. the - // reason is that we finally will - // want to have them in this order - // anyway, as this is the format we - // need inside deal.II - TensorProductPolynomials - face_polynomials (Polynomials::LagrangeEquidistant:: - generate_complete_basis (this->degree)); - Assert (face_polynomials.n() == this->dofs_per_face, ExcInternalError()); - - const unsigned int n_small_functions = this->interface_constraints_size()[0]; - - FullMatrix face_interpolation (n_small_functions, this->dofs_per_face); - FullMatrix subface_interpolation (n_small_functions, n_small_functions); - - for (unsigned int i=0; i > constraint_points; + // Add midpoint + constraint_points.push_back (Point (0.5)); + + if (this->degree>1) { - // generate a quadrature point - // xi. it is actually not so - // important where this point - // lies, as long as we make - // sure that they are not - // equal. however, we will want - // them to be the (equidistant) - // Lagrange points, since then - // the subface_interpolation - // matrix has a most positive - // property: it is a - // permutation of the identity - // matrix. so create an - // equidistant mesh of points - // in the interior of the face - // (in 2d). for 3d, things are - // somewhat more convoluted as - // usual, since the new (child) - // shape functions are not only - // located in the interior of - // the face, but also on the - // edges, with the exception of - // the four vertices of the - // face. the function we call - // takes care of all this - const Point p_face - = FE_Q_Helper::generate_face_unit_point (i, n_small_functions); - - // evaluate the big face - // shape function at this - // point. note that the - // numbering of our shape - // functions is different - // from that of the - // polynomial, which orders - // them in the order of - // interpolation points. - // - // face_index_map_inverse will - // get us over this little - // conversion - for (unsigned int j=0; jdofs_per_face; ++j) - { - face_interpolation(i,j) - = face_polynomials.compute_value(face_index_map[j], p_face); - // if the value is small up - // to round-off, then - // simply set it to zero to - // avoid unwanted fill-in - // of the constraint - // matrices (which would - // then increase the number - // of other DoFs a - // constrained DoF would - // couple to) - if (std::fabs(face_interpolation(i,j)) < 1e-14) - face_interpolation(i,j) = 0; - } - - // then evaluate all the - // small shape functions at - // this point. - for (unsigned int j=0; j::child_cell_from_point (p_face); - - // then check whether small - // shape function number j - // is nonzero on this - // face. as usual with our - // numbering of shape - // functions in constraint - // matrices, this is messy, - // so have a function that - // does this for us - // - // if not active, then the - // entry in the matrix will - // remain zero, and we - // simply go on with the - // next entry - if (! - FE_Q_Helper:: - constraint_function_is_active_on_child (j, subface, *this)) - continue; - - // otherwise: compute the - // coordinates of this - // evaluation point on - // the small face - const Point p_subface - = GeometryInfo::cell_to_child_coordinates (p_face, subface); - - // then get the index of - // small shape function j - // on this subface. again, - // divert to a function - // that is specialized for - // this - const unsigned int local_j - = FE_Q_Helper::constraint_get_local_j (j, subface, *this); - - // so evaluate this shape - // function there. now, - // since we have been - // careful with our choice - // of evaluation points, - // this is not actually - // necessary: the values of - // the small shape - // functions at these - // points should be either - // zero, and we can - // precompute which they - // are. However, we double - // check just to be sure we - // didn't do something - // stupid... - // - // (we could just set the - // evaluated value, but - // we'd end up with a lot - // of almost-zero entries, - // which will then carry - // over to the final - // result. this clutters up - // the constraint matrices, - // which we want to keep as - // small as possible.) - if (FE_Q_Helper::constraint_is_support_point (i, j, subface, *this)) - subface_interpolation(i, j) = 1.; - else - subface_interpolation(i, j) = 0.; - Assert (std::fabs (subface_interpolation(i, j) - - face_polynomials.compute_value(local_j, p_subface)) - < 1e-12, - ExcInternalError()); - } + const unsigned int n=this->degree-1; + const double step=1./this->degree; + // subface 0 + for (unsigned int i=1; i<=n; ++i) + constraint_points.push_back ( + GeometryInfo::child_to_cell_coordinates(Point(i*step),0)); + // subface 1 + for (unsigned int i=1; i<=n; ++i) + constraint_points.push_back ( + GeometryInfo::child_to_cell_coordinates(Point(i*step),1)); } - // what we now want to do is to - // compute - // (subface_intp)^-1 face_intp - // which should give us the - // desired hanging node constraints. - // rather than actually doing this, - // we note that we have constructed - // subface_interpolation to be a - // permutation of the unit matrix. - // rather than doing a gauss jordan - // inversion, we note that the - // inverse is actually given by the - // transpose of the matrix. This has - // the additional benefit of being - // more stable and in particular of - // not adding almost-zeros + // Now construct relation between + // destination (child) and source (mother) + // dofs. + const std::vector > polynomials= + Polynomials::LagrangeEquidistant::generate_complete_basis(this->degree); + this->interface_constraints .TableBase<2,double>::reinit (this->interface_constraints_size()); - subface_interpolation.Tmmult (this->interface_constraints, - face_interpolation); + + for (unsigned int i=0; idegree+1; ++j) + { + this->interface_constraints(i,j) = + polynomials[face_index_map[j]].value (constraint_points[i](0)); + + // if the value is small up + // to round-off, then + // simply set it to zero to + // avoid unwanted fill-in + // of the constraint + // matrices (which would + // then increase the number + // of other DoFs a + // constrained DoF would + // couple to) + if (std::fabs(this->interface_constraints(i,j)) < 1e-14) + this->interface_constraints(i,j) = 0; + } } #endif @@ -879,78 +723,25 @@ FE_Q<3>::initialize_constraints () { const unsigned int dim = 3; - // This algorithm for the automatic - // generation of the constraint matrices is - // different from the one implemented for - // the 2D elements. Hence it is only suited - // for standard Finite Elements with a - // Lagrangian basis. This algorithm - // consists of two parts. In the first - // part, the coordinates of the hanging - // nodes on the master element will be - // determined. These points are constructed - // in a special order (as described in the - // fe_base.h file for the class - // FiniteElementBase). First the hanging - // node in the mid of the coarser element - // is considered: - // - // Coarse Fine - // +-----+ +--+--+ - // | | | | | - // | * | +--+--+ - // | | | | | - // +-----+ +--+--+ - // - // Then the coordinates of the hanging - // nodes at the midpoint of the outline of - // the coarse element follow: - // - // Coarse Fine - // +--*--+ +--+--+ - // | | | | | - // * * +--+--+ - // | | | | | - // +--*--+ +--+--+ - // - // For Q1 that was it. But for higher order - // elements some more constraints are - // required. Hanging nodes on the lines - // which are inside the coarse element: - // - // Coarse Fine - // +-----+ +--+--+ - // | * | | | | - // | * * | +--+--+ - // | * | | | | - // +-----+ +--+--+ - // - // Hanging nodes on the outside lines: - // - // Coarse Fine - // +-*-*-+ +--+--+ - // * * | | | - // | | +--+--+ - // * * | | | - // +-*-*-+ +--+--+ - // - // And finally the interior nodes: - // - // Coarse Fine - // +-----+ +--+--+ - // | * * | | | | - // | | +--+--+ - // | * * | | | | - // +-----+ +--+--+ - // - // Once these points are known, it is - // pretty easy to get the contribution of - // each node on the coarse face to the - // value at the hanging nodes. This task - // is accomplished in the second part of - // the algorithm - - // Generate destination points. + // For a detailed documentation of + // the interpolation see the + // FE_Q<2>::initialize_constraints + // function. + + // In the following the points x_i + // are constructed in the order as + // described in the documentation + // of the FiniteElementBase class + // (fe_base.h), i.e. + // *--13--3--14--* + // | | | + // 16 20 7 19 12 + // | | | + // 4--8---0--6---2 + // | | | + // 15 17 5 18 11 + // | | | + // *--9---1--10--* std::vector > constraint_points; // Add midpoint