From 30e9e9e5b4c2132402b6b90c094bd0e214222708 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Wed, 6 Apr 2016 19:14:55 +0200 Subject: [PATCH] Doc changes and minor adjustments. --- source/fe/mapping_manifold.cc | 54 +++++++++++++++++++++-------------- 1 file changed, 32 insertions(+), 22 deletions(-) diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index 082879cd22..0137b9e2b9 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -119,8 +119,9 @@ initialize_face (const UpdateFlags update_flags, const unsigned int nfaces = GeometryInfo::faces_per_cell; unit_tangentials.resize (nfaces*(dim-1), std::vector > (n_original_q_points)); - switch(dim) { - case 2: + switch (dim) + { + case 2: { // ensure a counterclockwise // orientation of tangentials @@ -134,7 +135,7 @@ initialize_face (const UpdateFlags update_flags, } break; } - case 3: + case 3: { for (unsigned int i=0; i(data.vertices, data.cell_manifold_quadrature_weights[point+data_set])); - // Always get the maximum length from the point to the - // boundary of the reference element, to compute the - // tangent vectors from the Manifold object + // To compute the Jacobian, we choose dim points aligned + // with with the dim reference axes, which are still in the + // given cell, and ask for the tangent vector in these + // directions. Choosing the points is somewhat arbitrary, + // so we try to be smart and we pick points which are + // on the opposite quadrant w.r.t. the evaluation + // point. for (unsigned int i=0; i ei = Point::unit_vector(i); - const double ai = ei*p; - Assert(ai >=0, ExcInternalError("Was expecting a quadrature point " - "inside the unit reference element.")); - const Point np(ai > .5 ? p-ai *ei : p+(1-ai)*ei); + const double pi = p[i]; + Assert(pi >=0 && pi <= 1.0, + ExcInternalError("Was expecting a quadrature point " + "inside the unit reference element.")); + const Point np(pi > .5 ? p-pi *ei : p+(1-pi)*ei); - // In the lenghts, we store also the direction sign, + // In the length L, we store also the direction sign, // which is positive, if the coordinate is < .5, - double L = ai > .5 ? -ai: 1-ai; + double L = pi > .5 ? -pi: 1-pi; // Get the weights to compute the np point in real space for (unsigned int j=0; j::vertices_per_cell; ++j) data.vertex_weights[j] = GeometryInfo::d_linear_shape_function(np, j); - Point NP=data.manifold-> - get_new_point(Quadrature(data.vertices, - data.vertex_weights)); + const Point NP= + data.manifold->get_new_point(Quadrature(data.vertices, + data.vertex_weights)); Tensor<1,spacedim> T = data.manifold->get_tangent_vector(P, NP); @@ -563,7 +569,7 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, G[i][j] = DX_t[i] * DX_t[j]; output_data.JxW_values[point] - = sqrt(determinant(G)) * weights[point]; + = std::sqrt(determinant(G)) * weights[point]; if (cell_similarity == CellSimilarity::inverted_translation) { @@ -573,12 +579,16 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, } else { - const unsigned int codim = spacedim-dim; - (void)codim; - if (update_flags & update_normal_vectors) { - Assert( codim==1 , ExcMessage("There is no cell normal in codim 2.")); + Assert(spacedim == dim+1, + ExcMessage("There is no (unique) cell normal for " + + Utilities::int_to_string(dim) + + "-dimensional cells in " + + Utilities::int_to_string(spacedim) + + "-dimensional space. This only works if the " + "space dimension is one greater than the " + "dimensionality of the mesh cells.")); if (dim==1) output_data.normal_vectors[point] = -- 2.39.5