]> https://gitweb.dealii.org/ - dealii.git/commitdiff
push forward one coordinate at a time for jacobian derivatives
authorMaien Hamed <tomaien@hotmail.com>
Mon, 5 Oct 2015 16:24:39 +0000 (18:24 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 6 Oct 2015 06:47:32 +0000 (08:47 +0200)
source/fe/mapping_q_generic.cc

index 3a1bbda100e756e562ebf8674410272cbde39dd9..f85c3f77eec0ad4b98beb82a22510a043417d987 100644 (file)
@@ -2263,6 +2263,7 @@ namespace internal
 
           if (cell_similarity != CellSimilarity::translation)
             {
+              double tmp[spacedim][spacedim][spacedim];
               for (unsigned int point=0; point<n_q_points; ++point)
                 {
                   const Tensor<2,dim> *second =
@@ -2282,23 +2283,34 @@ namespace internal
                               *
                               data.mapping_support_points[k][i]);
 
-                  // pushing forward the derivative coordinates
+                  // first push forward the j-components
                   for (unsigned int i=0; i<spacedim; ++i)
                     for (unsigned int j=0; j<spacedim; ++j)
-                      for (unsigned int l=0; l<spacedim; ++l)
-                        {
-                          jacobian_pushed_forward_grads[point][i][j][l] = result[i][0][0] *
-                                                                          data.covariant[point][j][0] *
-                                                                          data.covariant[point][l][0];
-                          for (unsigned int jr=0; jr<dim; ++jr)
-                            {
-                              const unsigned int lr_start = jr==0? 1:0;
-                              for (unsigned int lr=lr_start; lr<dim; ++lr)
-                                jacobian_pushed_forward_grads[point][i][j][l] += result[i][jr][lr] *
-                                                                                 data.covariant[point][j][jr] *
-                                                                                 data.covariant[point][l][lr];
-                            }
-                        }
+                        for (unsigned int l=0; l<dim; ++l)
+                          {
+                            tmp[i][j][l] = result[i][0][l] *
+                                           data.covariant[point][j][0];
+                            for (unsigned int jr=1; jr<dim; ++jr)
+                              {
+                                tmp[i][j][l] += result[i][jr][l] *
+                                                data.covariant[point][j][jr];
+                              }
+                          }
+
+                  // now, pushing forward the l-components
+                  for (unsigned int i=0; i<spacedim; ++i)
+                    for (unsigned int j=0; j<spacedim; ++j)
+                        for (unsigned int l=0; l<spacedim; ++l)
+                          {
+                            jacobian_pushed_forward_grads[point][i][j][l] = tmp[i][j][0] *
+                                                                            data.covariant[point][l][0];
+                            for (unsigned int lr=1; lr<dim; ++lr)
+                              {
+                                jacobian_pushed_forward_grads[point][i][j][l] += tmp[i][j][lr] *
+                                                                                data.covariant[point][l][lr];
+                              }
+
+                          }
                 }
             }
         }
@@ -2376,6 +2388,7 @@ namespace internal
 
           if (cell_similarity != CellSimilarity::translation)
             {
+              double tmp[spacedim][spacedim][spacedim][spacedim];
               for (unsigned int point=0; point<n_q_points; ++point)
                 {
                   const Tensor<3,dim> *third =
@@ -2397,27 +2410,49 @@ namespace internal
                                 *
                                 data.mapping_support_points[k][i]);
 
+                  // push forward the j-coordinate
+                  for (unsigned int i=0; i<spacedim; ++i)
+                    for (unsigned int j=0; j<spacedim; ++j)
+                      for (unsigned int l=0; l<dim; ++l)
+                        for (unsigned int m=0; m<dim; ++m)
+                          {
+                            jacobian_pushed_forward_2nd_derivatives[point][i][j][l][m]
+                              = result[i][0][l][m]*
+                                data.covariant[point][j][0];
+                            for (unsigned int jr=1; jr<dim; ++jr)
+                              jacobian_pushed_forward_2nd_derivatives[point][i][j][l][m]
+                              += result[i][jr][l][m]*
+                                 data.covariant[point][j][jr];
+                          }
+
+                  // push forward the l-coordinate
+                  for (unsigned int i=0; i<spacedim; ++i)
+                    for (unsigned int j=0; j<spacedim; ++j)
+                      for (unsigned int l=0; l<spacedim; ++l)
+                        for (unsigned int m=0; m<dim; ++m)
+                          {
+                            tmp[i][j][l][m]
+                              = jacobian_pushed_forward_2nd_derivatives[point][i][j][0][m]*
+                                data.covariant[point][l][0];
+                            for (unsigned int lr=1; lr<dim; ++lr)
+                              tmp[i][j][l][m]
+                              += jacobian_pushed_forward_2nd_derivatives[point][i][j][lr][m]*
+                                 data.covariant[point][l][lr];
+                          }
+
+                  // push forward the m-coordinate
                   for (unsigned int i=0; i<spacedim; ++i)
                     for (unsigned int j=0; j<spacedim; ++j)
                       for (unsigned int l=0; l<spacedim; ++l)
                         for (unsigned int m=0; m<spacedim; ++m)
                           {
                             jacobian_pushed_forward_2nd_derivatives[point][i][j][l][m]
-                              = result[i][0][0][0]*
-                                data.covariant[point][j][0] *
-                                data.covariant[point][l][0] *
+                              = tmp[i][j][l][0]*
                                 data.covariant[point][m][0];
-                            for (unsigned int jr=0; jr<dim; ++jr)
-                              for (unsigned int lr=0; lr<dim; ++lr)
-                                {
-                                  const unsigned int mr0 = (jr+lr == 0)? 1:0;
-                                  for (unsigned int mr=mr0; mr<dim; ++mr)
-                                    jacobian_pushed_forward_2nd_derivatives[point][i][j][l][m]
-                                    += result[i][jr][lr][mr] *
-                                       data.covariant[point][j][jr] *
-                                       data.covariant[point][l][lr] *
-                                       data.covariant[point][m][mr];
-                                }
+                            for (unsigned int mr=1; mr<dim; ++mr)
+                              jacobian_pushed_forward_2nd_derivatives[point][i][j][l][m]
+                              += tmp[i][j][l][mr]*
+                                 data.covariant[point][m][mr];
                           }
                 }
             }
@@ -2498,6 +2533,7 @@ namespace internal
 
           if (cell_similarity != CellSimilarity::translation)
             {
+              double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
               for (unsigned int point=0; point<n_q_points; ++point)
                 {
                   const Tensor<4,dim> *fourth =
@@ -2521,6 +2557,53 @@ namespace internal
                                   *
                                   data.mapping_support_points[k][i]);
 
+                  // push-forward the j-coordinate
+                  for (unsigned int i=0; i<spacedim; ++i)
+                    for (unsigned int j=0; j<spacedim; ++j)
+                      for (unsigned int l=0; l<dim; ++l)
+                        for (unsigned int m=0; m<dim; ++m)
+                          for (unsigned int n=0; n<dim; ++n)
+                            {
+                              tmp[i][j][l][m][n] = result[i][0][l][m][n] *
+                                                   data.covariant[point][j][0];
+                              for (unsigned int jr=1; jr<dim; ++jr)
+                                tmp[i][j][l][m][n] += result[i][jr][l][m][n] *
+                                                      data.covariant[point][j][jr];
+                            }
+
+                  // push-forward the l-coordinate
+                  for (unsigned int i=0; i<spacedim; ++i)
+                    for (unsigned int j=0; j<spacedim; ++j)
+                      for (unsigned int l=0; l<spacedim; ++l)
+                        for (unsigned int m=0; m<dim; ++m)
+                          for (unsigned int n=0; n<dim; ++n)
+                            {
+                              jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
+                                = tmp[i][j][0][m][n] *
+                                  data.covariant[point][l][0];
+                              for (unsigned int lr=1; lr<dim; ++lr)
+                                jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
+                                += tmp[i][j][lr][m][n] *
+                                   data.covariant[point][l][lr];
+                            }
+
+                  // push-forward the m-coordinate
+                  for (unsigned int i=0; i<spacedim; ++i)
+                    for (unsigned int j=0; j<spacedim; ++j)
+                      for (unsigned int l=0; l<spacedim; ++l)
+                        for (unsigned int m=0; m<spacedim; ++m)
+                          for (unsigned int n=0; n<dim; ++n)
+                            {
+                              tmp[i][j][l][m][n]
+                                = jacobian_pushed_forward_3rd_derivatives[point][i][j][l][0][n] *
+                                  data.covariant[point][m][0];
+                              for (unsigned int mr=1; mr<dim; ++mr)
+                                tmp[i][j][l][m][n]
+                                += jacobian_pushed_forward_3rd_derivatives[point][i][j][l][mr][n] *
+                                   data.covariant[point][m][mr];
+                            }
+
+                  // push-forward the n-coordinate
                   for (unsigned int i=0; i<spacedim; ++i)
                     for (unsigned int j=0; j<spacedim; ++j)
                       for (unsigned int l=0; l<spacedim; ++l)
@@ -2528,24 +2611,12 @@ namespace internal
                           for (unsigned int n=0; n<spacedim; ++n)
                             {
                               jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
-                                = result[i][0][0][0][0] *
-                                  data.covariant[point][j][0] *
-                                  data.covariant[point][l][0] *
-                                  data.covariant[point][m][0] *
+                                = tmp[i][j][l][m][0] *
                                   data.covariant[point][n][0];
-                              for (unsigned int jr=0; jr<dim; ++jr)
-                                for (unsigned int lr=0; lr<dim; ++lr)
-                                  for (unsigned int mr=0; mr<dim; ++mr)
-                                    {
-                                      const unsigned int nr0 = (jr+lr+mr==0)? 1:0;
-                                      for (unsigned int nr=nr0; nr<dim; ++nr)
-                                        jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
-                                        += result[i][jr][lr][mr][nr] *
-                                           data.covariant[point][j][jr] *
-                                           data.covariant[point][l][lr] *
-                                           data.covariant[point][m][mr] *
-                                           data.covariant[point][n][nr];
-                                    }
+                              for (unsigned int nr=1; nr<dim; ++nr)
+                                jacobian_pushed_forward_3rd_derivatives[point][i][j][l][m][n]
+                                += tmp[i][j][l][m][nr] *
+                                   data.covariant[point][n][nr];
                             }
                 }
             }

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.