]> https://gitweb.dealii.org/ - dealii.git/commitdiff
deprecate functions, simplify code, cleanup
authorTimo Heister <timo.heister@gmail.com>
Tue, 2 Aug 2016 19:32:14 +0000 (21:32 +0200)
committerTimo Heister <timo.heister@gmail.com>
Tue, 9 Aug 2016 13:43:04 +0000 (09:43 -0400)
- cleanup
- update tests
- deprecate some functions
- simplify code

include/deal.II/integrators/divergence.h
include/deal.II/integrators/grad_div.h
tests/integrators/cochain_01.cc
tests/integrators/grad_div_01.cc

index 3dedf483cb55baa166c7fba50e414c7b98fb733e..791b671973afa13b28a5afec132cf2554cfad941 100644 (file)
@@ -24,6 +24,7 @@
 #include <deal.II/fe/mapping.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/meshworker/dof_info.h>
+#include <deal.II/integrators/grad_div.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -39,35 +40,6 @@ namespace LocalIntegrators
    */
   namespace Divergence
   {
-    /**
-     * @deprecated Use LocalIntegrators::GradDiv::
-     * Auxiliary function. Computes the grad-div-operator from a set of
-     * Hessians.
-     *
-     * @note The third tensor argument is not used in two dimensions and can
-     * for instance duplicate one of the previous.
-     *
-     * @author Guido Kanschat
-     * @date 2011
-     */
-    template <int dim>
-    Tensor<1,dim>
-    grad_div(
-      const Tensor<2,dim> &h0,
-      const Tensor<2,dim> &h1,
-      const Tensor<2,dim> &h2)
-    {
-      Tensor<1,dim> result;
-      for (unsigned int d=0; d<dim; ++d)
-        {
-          result[d] += h0[d][0];
-          if (dim >=2) result[d] += h1[d][1];
-          if (dim >=3) result[d] += h2[d][2];
-        }
-      return result;
-    }
-
-
     /**
      * Cell matrix for divergence. The derivative is on the trial function.
      * \f[ \int_Z v\nabla \cdot \mathbf u \,dx \f] This is the strong
@@ -452,77 +424,41 @@ namespace LocalIntegrators
     }
 
     /**
-     * @deprecated Use LocalIntegrators::GradDiv::cell_matrix()
-     * instead.
-     *
-     * @author Guido Kanschat
-     * @date 2011
+     * @deprecated Use LocalIntegrators::GradDiv::cell_matrix() instead.
      */
     template <int dim>
     void grad_div_matrix (
       FullMatrix<double> &M,
       const FEValuesBase<dim> &fe,
-      double factor = 1.) DEAL_II_DEPRECATED
-    {
-      const unsigned int n_dofs = fe.dofs_per_cell;
-
-      AssertDimension(fe.get_fe().n_components(), dim);
-      AssertDimension(M.m(), n_dofs);
-      AssertDimension(M.n(), n_dofs);
-
-      for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
-        {
-          const double dx = factor * fe.JxW(k);
-          for (unsigned int i=0; i<n_dofs; ++i)
-            for (unsigned int j=0; j<n_dofs; ++j)
-              {
-                double dv = 0.;
-                double du = 0.;
-                for (unsigned int d=0; d<dim; ++d)
-                  {
-                    dv += fe.shape_grad_component(i,k,d)[d];
-                    du += fe.shape_grad_component(j,k,d)[d];
-                  }
+      const double factor = 1.) DEAL_II_DEPRECATED;
 
-                M(i,j) += dx * du * dv;
-              }
-        }
+    template <int dim>
+    void grad_div_matrix (
+      FullMatrix<double> &M,
+      const FEValuesBase<dim> &fe,
+      const double factor)
+    {
+      GradDiv::cell_matrix(M, fe, factor);
     }
 
     /**
      * @deprecated Use LocalIntegrators::GradDiv::cell_residual() instead.
-     *
-     * @author Guido Kanschat
-     * @date 2014
      */
     template <int dim, typename number>
     void grad_div_residual (
       Vector<number> &result,
       const FEValuesBase<dim> &fetest,
       const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &input,
-      const double factor = 1.) DEAL_II_DEPRECATED
-    {
-      const unsigned int n_dofs = fetest.dofs_per_cell;
-
-      AssertDimension(fetest.get_fe().n_components(), dim);
-      AssertVectorVectorDimension(input, dim, fetest.n_quadrature_points);
-
-      for (unsigned int k=0; k<fetest.n_quadrature_points; ++k)
-        {
-          const double dx = factor * fetest.JxW(k);
-          for (unsigned int i=0; i<n_dofs; ++i)
-            {
-              double dv = 0.;
-              double du = 0.;
-              for (unsigned int d=0; d<dim; ++d)
-                {
-                  dv += fetest.shape_grad_component(i,k,d)[d];
-                  du += input[d][k][d];
-                }
+      const double factor = 1.) DEAL_II_DEPRECATED;
 
-              result(i) += dx * du * dv;
-            }
-        }
+    template <int dim, typename number>
+    void grad_div_residual (
+      Vector<number> &result,
+      const FEValuesBase<dim> &fetest,
+      const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &input,
+      const double factor)
+    {
+      GradDiv::cell_residual(result, fetest, input, factor);
     }
 
     /**
index 2b2fa9714e03befefa255b2d18ac2364a52602f9..58b8868fe890472ad51d50aba06f25fa417e957a 100644 (file)
@@ -39,33 +39,6 @@ namespace LocalIntegrators
    */
   namespace GradDiv
   {
-    /**
-     * Auxiliary function. Computes the grad-div-operator from a set of
-     * Hessians.
-     *
-     * @note The third tensor argument is not used in two dimensions and can
-     * for instance duplicate one of the previous.
-     *
-     * @author Guido Kanschat
-     * @date 2011
-     */
-    template <int dim>
-    Tensor<1,dim>
-    grad_div(
-      const Tensor<2,dim> &h0,
-      const Tensor<2,dim> &h1,
-      const Tensor<2,dim> &h2)
-    {
-      Tensor<1,dim> result;
-      for (unsigned int d=0; d<dim; ++d)
-        {
-          result[d] += h0[d][0];
-          if (dim >=2) result[d] += h1[d][1];
-          if (dim >=3) result[d] += h2[d][2];
-        }
-      return result;
-    }
-
     /**
      * The weak form of the grad-div operator penalizing volume changes
      * @f[
@@ -93,15 +66,10 @@ namespace LocalIntegrators
           for (unsigned int i=0; i<n_dofs; ++i)
             for (unsigned int j=0; j<n_dofs; ++j)
               {
-                double dv = 0.;
-                double du = 0.;
-                for (unsigned int d=0; d<dim; ++d)
-                  {
-                    dv += fe.shape_grad_component(i,k,d)[d];
-                    du += fe.shape_grad_component(j,k,d)[d];
-                  }
-               
-                M(i,j) += dx * du * dv;
+                const double divu = fe[FEValuesExtractors::Vector(0)].divergence(j,k);
+                const double divv = fe[FEValuesExtractors::Vector(0)].divergence(i,k);
+
+                M(i,j) += dx * divu * divv;
               }
         }
     }
@@ -123,28 +91,25 @@ namespace LocalIntegrators
       const double factor = 1.)
     {
       const unsigned int n_dofs = fetest.dofs_per_cell;
-      
+
       AssertDimension(fetest.get_fe().n_components(), dim);
       AssertVectorVectorDimension(input, dim, fetest.n_quadrature_points);
-      
+
       for (unsigned int k=0; k<fetest.n_quadrature_points; ++k)
         {
           const double dx = factor * fetest.JxW(k);
           for (unsigned int i=0; i<n_dofs; ++i)
             {
-              double dv = 0.;
+              const double divv = fetest[FEValuesExtractors::Vector(0)].divergence(i,k);
               double du = 0.;
               for (unsigned int d=0; d<dim; ++d)
-                {
-                  dv += fetest.shape_grad_component(i,k,d)[d];
-                  du += input[d][k][d];
-                }
-             
-              result(i) += dx * du * dv;
+                du += input[d][k][d];
+
+              result(i) += dx * du * divv;
             }
         }
     }
-    
+
     /**
      * The matrix for the weak boundary condition of Nitsche type for linear elasticity:
      * @f[
@@ -172,13 +137,13 @@ namespace LocalIntegrators
           for (unsigned int i=0; i<n_dofs; ++i)
             for (unsigned int j=0; j<n_dofs; ++j)
               {
-                double un = 0., vn = 0., divu = 0., divv = 0.;
+                const double divu = fe[FEValuesExtractors::Vector(0)].divergence(j,k);
+                const double divv = fe[FEValuesExtractors::Vector(0)].divergence(i,k);
+                double un = 0., vn = 0.;
                 for (unsigned int d=0; d<dim; ++d)
                   {
                     un += fe.shape_value_component(j,k,d) * n[d];
                     vn += fe.shape_value_component(i,k,d) * n[d];
-                    divu += fe.shape_grad_component(j,k,d)[d];
-                    divv += fe.shape_grad_component(i,k,d)[d];
                   }
 
                 M(i,j) += dx * 2. * penalty * un * vn;
@@ -217,7 +182,7 @@ namespace LocalIntegrators
     {
       const unsigned int n_dofs = fe.dofs_per_cell;
       AssertDimension(fe.get_fe().n_components(), dim)
-       AssertVectorVectorDimension(input, dim, fe.n_quadrature_points);
+      AssertVectorVectorDimension(input, dim, fe.n_quadrature_points);
       AssertVectorVectorDimension(Dinput, dim, fe.n_quadrature_points);
       AssertVectorVectorDimension(data, dim, fe.n_quadrature_points);
 
@@ -225,30 +190,27 @@ namespace LocalIntegrators
         {
           const double dx = factor * fe.JxW(k);
           const Tensor<1,dim> n = fe.normal_vector(k);
-         
-         double umgn = 0.;
-         double divu = 0.;
-         for (unsigned int d=0; d<dim; ++d)
-           {
-             umgn += (input[d][k] - data[d][k]) * n[d];
-             divu += Dinput[d][k][d];
-           }
-         
+
+          double umgn = 0.;
+          double divu = 0.;
+          for (unsigned int d=0; d<dim; ++d)
+            {
+              umgn += (input[d][k] - data[d][k]) * n[d];
+              divu += Dinput[d][k][d];
+            }
+
           for (unsigned int i=0; i<n_dofs; ++i)
-           {
-             double vn = 0.;
-             double divv = 0.;
-             for (unsigned int d=0; d<dim; ++d)
-               {
-                 vn += fe.shape_value_component(i,k,d) * n[d];
-                 divv += fe.shape_grad_component(i,k,d)[d];
-               }
-             
-             result(i) += dx*(2.*penalty*umgn*vn - divv*umgn - divu*vn);
-           }
+            {
+              double vn = 0.;
+              const double divv = fe[FEValuesExtractors::Vector(0)].divergence(i,k);
+              for (unsigned int d=0; d<dim; ++d)
+                vn += fe.shape_value_component(i,k,d) * n[d];
+
+              result(i) += dx*(2.*penalty*umgn*vn - divv*umgn - divu*vn);
+            }
         }
     }
-    
+
     /**
      * The interior penalty flux for the grad-div operator. See
      * ip_residual() for details.
@@ -288,34 +250,30 @@ namespace LocalIntegrators
           const double dx = fe1.JxW(k);
           const Tensor<1,dim> n = fe1.normal_vector(k);
           for (unsigned int i=0; i<n_dofs; ++i)
-           for (unsigned int j=0; j<n_dofs; ++j)
-             {
-               double uni = 0.;
-               double une = 0.;
-               double vni = 0.;
-               double vne = 0.;
-               double divui = 0.;
-               double divue = 0.;
-               double divvi = 0.;
-               double divve = 0.;
-               
-               for (unsigned int d=0; d<dim; ++d)
-                 {
-                   uni += fe1.shape_value_component(j,k,d) * n[d];
-                   une += fe2.shape_value_component(j,k,d) * n[d];
-                   vni += fe1.shape_value_component(i,k,d) * n[d];
-                   vne += fe2.shape_value_component(i,k,d) * n[d];
-                   divui += fe1.shape_grad_component(j,k,d)[d];
-                   divue += fe2.shape_grad_component(j,k,d)[d];
-                   divvi += fe1.shape_grad_component(i,k,d)[d];
-                   divve += fe2.shape_grad_component(i,k,d)[d];
-                 }
-               M11(i,j) += dx*(-.5*fi*divvi*uni-.5*fi*divui*vni+f*penalty*uni*vni);
-               M12(i,j) += dx*( .5*fi*divvi*une-.5*fe*divue*vni-f*penalty*vni*une);
-               M21(i,j) += dx*(-.5*fe*divve*uni+.5*fi*divui*vne-f*penalty*uni*vne);
-               M22(i,j) += dx*( .5*fe*divve*une+.5*fe*divue*vne+f*penalty*une*vne);
-             }
-       }
+            for (unsigned int j=0; j<n_dofs; ++j)
+              {
+                double uni = 0.;
+                double une = 0.;
+                double vni = 0.;
+                double vne = 0.;
+                const double divui = fe1[FEValuesExtractors::Vector(0)].divergence(j,k);
+                const double divue = fe2[FEValuesExtractors::Vector(0)].divergence(j,k);
+                const double divvi = fe1[FEValuesExtractors::Vector(0)].divergence(i,k);
+                const double divve = fe2[FEValuesExtractors::Vector(0)].divergence(i,k);
+
+                for (unsigned int d=0; d<dim; ++d)
+                  {
+                    uni += fe1.shape_value_component(j,k,d) * n[d];
+                    une += fe2.shape_value_component(j,k,d) * n[d];
+                    vni += fe1.shape_value_component(i,k,d) * n[d];
+                    vne += fe2.shape_value_component(i,k,d) * n[d];
+                  }
+                M11(i,j) += dx*(-.5*fi*divvi*uni-.5*fi*divui*vni+f*penalty*uni*vni);
+                M12(i,j) += dx*( .5*fi*divvi*une-.5*fe*divue*vni-f*penalty*vni*une);
+                M21(i,j) += dx*(-.5*fe*divve*uni+.5*fi*divui*vne-f*penalty*uni*vne);
+                M22(i,j) += dx*( .5*fe*divve*une+.5*fe*divue*vne+f*penalty*une*vne);
+              }
+        }
     }
 
     /**
@@ -364,37 +322,35 @@ namespace LocalIntegrators
         {
           const double dx = fe1.JxW(k);
           const Tensor<1,dim> n = fe1.normal_vector(k);
-         double uni = 0.;
-         double une = 0.;
-         double divui = 0.;
-         double divue = 0.;
-         for (unsigned int d=0; d<dim; ++d)
-           {
-             uni += input1[d][k]*n[d];
-             une += input2[d][k]*n[d];
-             divui += Dinput1[d][k][d];
-             divue += Dinput2[d][k][d];              
-           }
-         
+          double uni = 0.;
+          double une = 0.;
+          double divui = 0.;
+          double divue = 0.;
+          for (unsigned int d=0; d<dim; ++d)
+            {
+              uni += input1[d][k]*n[d];
+              une += input2[d][k]*n[d];
+              divui += Dinput1[d][k][d];
+              divue += Dinput2[d][k][d];
+            }
+
           for (unsigned int i=0; i<n1; ++i)
-           {
-             double vni = 0.;
-             double vne = 0.;
-             double divvi = 0.;
-             double divve = 0.;
-             for (unsigned int d=0; d<dim; ++d)
-               {
-                 vni += fe1.shape_value_component(i,k,d)*n[d];
-                 vne += fe2.shape_value_component(i,k,d)*n[d];
-                 divvi += fe1.shape_grad_component(i,k,d)[d];
-                 divve += fe2.shape_grad_component(i,k,d)[d];        
-               }
-             
-             result1(i) += dx*(-.5*fi*divvi*uni-.5*fi*divui*vni+penalty*uni*vni);
-             result1(i) += dx*( .5*fi*divvi*une-.5*fe*divue*vni-penalty*vni*une);
-             result2(i) += dx*(-.5*fe*divve*uni+.5*fi*divui*vne-penalty*uni*vne);
-             result2(i) += dx*( .5*fe*divve*une+.5*fe*divue*vne+penalty*une*vne);
-           }
+            {
+              double vni = 0.;
+              double vne = 0.;
+              const double divvi = fe1[FEValuesExtractors::Vector(0)].divergence(i,k);
+              const double divve = fe2[FEValuesExtractors::Vector(0)].divergence(i,k);
+              for (unsigned int d=0; d<dim; ++d)
+                {
+                  vni += fe1.shape_value_component(i,k,d)*n[d];
+                  vne += fe2.shape_value_component(i,k,d)*n[d];
+                }
+
+              result1(i) += dx*(-.5*fi*divvi*uni-.5*fi*divui*vni+penalty*uni*vni);
+              result1(i) += dx*( .5*fi*divvi*une-.5*fe*divue*vni-penalty*vni*une);
+              result2(i) += dx*(-.5*fe*divve*uni+.5*fi*divui*vne-penalty*uni*vne);
+              result2(i) += dx*( .5*fe*divve*une+.5*fe*divue*vne+penalty*une*vne);
+            }
         }
     }
   }
index 65bb5b8387593ab0bed5509497aeba0fd3c4aebf..671747cec754ef305e96c4021e80dd2c50efc8ab 100644 (file)
@@ -47,6 +47,7 @@
 
 #include <deal.II/integrators/l2.h>
 #include <deal.II/integrators/divergence.h>
+#include <deal.II/integrators/grad_div.h>
 #include <deal.II/integrators/laplace.h>
 #include <deal.II/integrators/maxwell.h>
 #include <deal.II/integrators/elasticity.h>
@@ -77,7 +78,7 @@ void cell_matrix(
     {
       ++de;
       L2::mass_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de));
-      Divergence::grad_div_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de));
+      GradDiv::cell_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de));
       Divergence::cell_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de), info.fe_values(de+1));
     }
 
index 2bcc6db40b5f9b49d2796009b79f2ff82f57511c..23773a42b182b4a84794fb9d2065a516f14d73cb 100644 (file)
@@ -14,7 +14,7 @@
 // ---------------------------------------------------------------------
 
 
-// Test the functions in integrators/elasticity.h
+// Test the functions in integrators/grad_div.h
 // Output matrices and assert consistency of residuals
 
 #include "../tests.h"

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.