]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Clean up code, fix comments by David. 5410/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 15 Nov 2017 15:23:09 +0000 (16:23 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 16 Nov 2017 07:04:12 +0000 (08:04 +0100)
include/deal.II/grid/manifold_lib.h
source/grid/manifold_lib.cc

index fa58b747e488e407c17c6f0c4646a8de9e7ad927..e4434e8121894797ffe572aa6ffd84c6bc12f891 100644 (file)
@@ -805,11 +805,12 @@ private:
    * Gradient of the push_forward method.
    *
    * @note This internal function is not compatible with the
-   * ChartManifold::pull_back() function because the given class represents an
-   * atlas of charts, not a single chart. Furthermore, this private function
-   * also requires the user to provide the result of the push_forward() call
-   * on the chart point for the use case of this function, namely inside a
-   * Newton iteration where the gradient is computed by finite differences.
+   * ChartManifold::push_forward_gradient() function because the given class
+   * represents an atlas of charts, not a single chart. Furthermore, this
+   * private function also requires the user to provide the result of the
+   * push_forward() call on the chart point for the single use case of this
+   * function, namely inside a Newton iteration where the gradient is computed
+   * by finite differences.
    */
   DerivativeForm<1,dim,spacedim>
   push_forward_gradient(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
index f7b6eea6cc89b58e7b8f98369796a96f7f3679c4..e4d0b77e7c5cd4a167f5dac4d1b85a643b61af1b 100644 (file)
@@ -941,32 +941,31 @@ namespace
     const std::array<Point<spacedim>, 4> vertices
     {{cell.vertex(0), cell.vertex(1), cell.vertex(2), cell.vertex(3)}};
 
-    // store the components of the linear shape functions because we need them
-    // repeatedly
-    double linear_shapes[4];
-    for (unsigned int d=0; d<2; ++d)
+    // this evaluates all bilinear shape functions because we need them
+    // repeatedly. we will update this values in the complicated case with
+    // curved lines below
+    std::array<double,4> weights_vertices
+    {
       {
-        linear_shapes[2*d]   = 1.-chart_point[d];
-        linear_shapes[2*d+1] = chart_point[d];
+        (1.-chart_point[0]) *(1.-chart_point[1]),
+        chart_point[0] *(1.-chart_point[1]),
+        (1.-chart_point[0]) *chart_point[1],
+        chart_point[0] *chart_point[1]
       }
+    };
 
     Point<spacedim> new_point;
     if (cell_is_flat)
-      for (unsigned int i1=0, v=0; i1<2; ++i1)
-        for (unsigned int i0=0; i0<2; ++i0, ++v)
-          new_point += linear_shapes[2+i1] * linear_shapes[i0] * vertices[v];
+      for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
+        new_point += weights_vertices[v] * vertices[v];
     else
       {
-        // We subtract the contribution of the vertices (second line in formula).
-        // If a line applies the same manifold as the cell, we also subtract a
-        // weighted part of the vertex, so accumulate the final weight of the
-        // the vertices while going through the faces (this is a bit artificial
-        // in 2D but it becomes clear in 3D where we avoid looking at the faces'
-        // orientation and other complications).
-        double weights_vertices[GeometryInfo<2>::vertices_per_cell];
-        for (unsigned int i1=0, v=0; i1<2; ++i1)
-          for (unsigned int i0=0; i0<2; ++i0, ++v)
-            weights_vertices[v] = -linear_shapes[2+i1] * linear_shapes[i0];
+        // The second line in the formula tells us to subtract the
+        // contribution of the vertices.  If a line employs the same manifold
+        // as the cell, we can merge the weights of the line with the weights
+        // of the vertex with a negative sign while going through the faces
+        // (this is a bit artificial in 2D but it becomes clear in 3D where we
+        // avoid looking at the faces' orientation and other complications).
 
         // add the contribution from the lines around the cell (first line in
         // formula)
@@ -981,16 +980,18 @@ namespace
             const double my_weight = line%2 ? chart_point[line/2] : 1-chart_point[line/2];
             const double line_point = chart_point[1-line/2];
 
-            // Same manifold or invalid id which will go back to the same class ->
-            // adds to the vertices
+            // Same manifold or invalid id which will go back to the same
+            // class -> contribution should be added for the final point,
+            // which means that we subtract the current weight from the
+            // negative weight applied to the vertex
             const types::manifold_id line_manifold_id = cell.line(line)->manifold_id();
             if (line_manifold_id == my_manifold_id ||
                 line_manifold_id == numbers::invalid_manifold_id)
               {
                 weights_vertices[GeometryInfo<2>::line_to_cell_vertices(line,0)]
-                += my_weight * (1.-line_point);
+                -= my_weight * (1.-line_point);
                 weights_vertices[GeometryInfo<2>::line_to_cell_vertices(line,1)]
-                += my_weight * line_point;
+                -= my_weight * line_point;
               }
             else
               {
@@ -1006,7 +1007,7 @@ namespace
 
         // subtract contribution from the vertices (second line in formula)
         for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
-          new_point += weights_vertices[v] * vertices[v];
+          new_point -= weights_vertices[v] * vertices[v];
       }
 
     return new_point;
@@ -1072,22 +1073,25 @@ namespace
         linear_shapes[2*d+1] = chart_point[d];
       }
 
+    // wrap linear shape functions around for access in face loop
+    for (unsigned int d=6; d<10; ++d)
+      linear_shapes[d] = linear_shapes[d-6];
+
+    std::array<double,8> weights_vertices;
+    for (unsigned int i2=0, v=0; i2<2; ++i2)
+      for (unsigned int i1=0; i1<2; ++i1)
+        for (unsigned int i0=0; i0<2; ++i0, ++v)
+          weights_vertices[v] = (linear_shapes[4+i2] * linear_shapes[2+i1]) *
+                                linear_shapes[i0];
+
     Point<spacedim> new_point;
     if (cell_is_flat)
-      for (unsigned int i2=0, v=0; i2<2; ++i2)
-        for (unsigned int i1=0; i1<2; ++i1)
-          for (unsigned int i0=0; i0<2; ++i0, ++v)
-            new_point += (linear_shapes[4+i2] * linear_shapes[2+i1]) *
-                         linear_shapes[i0] * vertices[v];
+      for (unsigned int v=0; v<8; ++v)
+        new_point += weights_vertices[v] * vertices[v];
     else
       {
-        // identify the weights for the vertices and lines to be accumulated
-        double weights_vertices[GeometryInfo<3>::vertices_per_cell];
-        for (unsigned int i2=0, v=0; i2<2; ++i2)
-          for (unsigned int i1=0; i1<2; ++i1)
-            for (unsigned int i0=0; i0<2; ++i0, ++v)
-              weights_vertices[v] = (linear_shapes[4+i2] * linear_shapes[2+i1]) *
-                                    linear_shapes[i0];
+        // identify the weights for the lines to be accumulated (vertex
+        // weights are set outside and coincide with the flat manifold case)
 
         double weights_lines[GeometryInfo<3>::lines_per_cell];
         for (unsigned int line=0; line<GeometryInfo<3>::lines_per_cell; ++line)
@@ -1100,10 +1104,6 @@ namespace
         const auto weights_view = make_array_view(weights.begin(), weights.end());
         const auto points_view = make_array_view(points.begin(), points.end());
 
-        // wrap linear shape functions around for access in face loop
-        for (unsigned int d=6; d<10; ++d)
-          linear_shapes[d] = linear_shapes[d-6];
-
         for (unsigned int face=0; face<GeometryInfo<3>::faces_per_cell; ++face)
           {
             const double my_weight = linear_shapes[face];
@@ -1124,6 +1124,10 @@ namespace
                     weights_lines[face_to_cell_lines_3d[face][line]] +=
                       my_weight * line_weight;
                   }
+                // as to the indices inside linear_shapes: we use the index
+                // wrapped around at 2*d, ensuring the correct orientation of
+                // the face's coordinate system with respect to the
+                // lexicographic indices
                 weights_vertices[face_to_cell_vertices_3d[face][0]] -=
                   linear_shapes[face_even+2]*(linear_shapes[face_even+4]*my_weight);
                 weights_vertices[face_to_cell_vertices_3d[face][1]] -=
@@ -1254,7 +1258,7 @@ TransfiniteInterpolationManifold<dim,spacedim>
 {
   Point<dim> outside;
   for (unsigned int d=0; d<dim; ++d)
-    outside[d] = 20;
+    outside[d] = 20.;
 
   // project the user-given input to unit cell
   Point<dim> chart_point = GeometryInfo<dim>::project_to_unit_cell(initial_guess);
@@ -1295,8 +1299,11 @@ TransfiniteInterpolationManifold<dim,spacedim>
       // convergence
       if (i%8 == 0)
         {
-          // if the determinant is zero, the mapping is not invertible and we are
-          // outside the valid chart region
+          // if the determinant is zero or negative, the mapping is either not
+          // invertible or already has inverted and we are outside the valid
+          // chart region. Note that the Jacobian here represents the
+          // derivative of the forward map and should have a positive
+          // determinant since we use properly oriented meshes.
           DerivativeForm<1,dim,spacedim> grad
             = push_forward_gradient(cell, chart_point,
                                     Point<spacedim>(point-residual));
@@ -1481,7 +1488,7 @@ TransfiniteInterpolationManifold<dim, spacedim>
           else if (surrounding_points.size() == 8 && i > 4)
             {
               const Point<dim> guess = chart_points[i-4] +
-                                       Point<dim>(chart_points[4] - chart_points[0]);
+                                       (chart_points[4] - chart_points[0]);
               chart_points[i] = pull_back(cell, surrounding_points[i], guess);
             }
           else

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.