]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Use functionality to evaluate finite element functions and fix bugs in error estimator.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 29 Apr 1998 12:09:38 +0000 (12:09 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 29 Apr 1998 12:09:38 +0000 (12:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@224 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/numerics/assembler.cc
deal.II/deal.II/source/numerics/base.cc
deal.II/deal.II/source/numerics/error_estimator.cc

index a88b02e5217bc2a551ca94cc4d4c5a7df68b6143..475545b464b0bfaca294f1fe6dad3c36396867be 100644 (file)
@@ -101,7 +101,8 @@ void Assembler<dim>::assemble (const Equation<dim> &equation) {
                                   // re-init fe values for this cell
   fe_values.reinit (DoFHandler<dim>::cell_iterator (tria,
                                                    present_level,
-                                                   present_index),
+                                                   present_index,
+                                                   dof_handler),
                    fe,
                    boundary);
   const unsigned int n_dofs = dof_handler->get_selected_fe().total_dofs;
index da58fa61af4cf96eaa28293297314af3b2435279..b411e92bcd9e01686b1578843c81f0b16920d977 100644 (file)
@@ -211,12 +211,7 @@ void ProblemBase<dim>::integrate_difference (const Function<dim>      &exact_sol
                                             // reference_function(x_j)-\psi_j
                                             // and assign that to the vector
                                             // \psi.
-           const unsigned int n_dofs = fe.total_dofs;
            const unsigned int n_q_points = q.n_quadrature_points;
-           const dFMatrix & shape_values = fe_values.get_shape_values();
-           vector<double>   dof_values(fe.total_dofs, 0);
-           cell->get_dof_values (solution, dof_values);
-           
            vector<double>   psi;
 
                                             // in praxi: first compute
@@ -225,9 +220,17 @@ void ProblemBase<dim>::integrate_difference (const Function<dim>      &exact_sol
                                       psi);
                                             // then subtract finite element
                                             // solution
-           for (unsigned int j=0; j<n_q_points; ++j) 
-             for (unsigned int i=0; i<n_dofs; ++i)
-               psi[j] -= dof_values[i]*shape_values(i,j);
+           if (true) 
+             {
+               vector<double> function_values (n_q_points, 0);
+               fe_values.get_function_values (solution, function_values);
+
+               transform (psi.begin(), psi.end(),
+                          function_values.begin(),
+                          psi.begin(),
+                          minus<double>());
+             };
+           
 
                                             // for L1_norm and Linfty_norm:
                                             // take absolute
@@ -296,12 +299,7 @@ void ProblemBase<dim>::integrate_difference (const Function<dim>      &exact_sol
 
                                             // same procedure as above, but now
                                             // psi is a vector of gradients
-           const unsigned int n_dofs = fe.total_dofs;
            const unsigned int n_q_points = q.n_quadrature_points;
-           const vector<vector<Point<dim> > > & shape_grads = fe_values.get_shape_grads();
-           vector<double>   dof_values(fe.total_dofs, 0);
-           cell->get_dof_values (solution, dof_values);
-           
            vector<Point<dim> >   psi;
 
                                             // in praxi: first compute
@@ -311,10 +309,16 @@ void ProblemBase<dim>::integrate_difference (const Function<dim>      &exact_sol
            
                                             // then subtract finite element
                                             // solution
-           for (unsigned int j=0; j<n_q_points; ++j) 
-             for (unsigned int i=0; i<n_dofs; ++i)
-               psi[j] -= dof_values[i]*shape_grads[i][j];
+           if (true) 
+             {
+               vector<Point<dim> > function_grads (n_q_points, Point<dim>());
+               fe_values.get_function_grads (solution, function_grads);
 
+               transform (psi.begin(), psi.end(),
+                          function_grads.begin(),
+                          psi.begin(),
+                          minus<Point<dim> >());
+             };
                                             // take square of integrand
            vector<double> psi_square (psi.size(), 0.0);
            for (unsigned int i=0; i<n_q_points; ++i)
index 26f568a07c7bacace1221b49a8a88c67f26afdf7..24ba965881a255869434892f7047df23b6538677 100644 (file)
@@ -51,8 +51,6 @@ void KellyErrorEstimator<dim>::estimate_error (const DoFHandler<dim>    &dof,
 
                                   // number of integration points per face
   const unsigned int n_q_points = quadrature.n_quadrature_points;
-                                  // number of dofs per cell
-  const unsigned int n_dofs = fe.total_dofs;
   
                                   // make up a fe face values object for the
                                   // restriction of the finite element function
@@ -98,27 +96,8 @@ void KellyErrorEstimator<dim>::estimate_error (const DoFHandler<dim>    &dof,
                                         // let psi be a short name for
                                         // [grad u_h]
        vector<Point<dim> > psi(n_q_points);
+       fe_face_values_cell.get_function_grads (solution, psi);
        
-                                        // get a list of the values of
-                                        // the solution for the ansatz
-                                        // functions on this cell
-       vector<double>   dof_values(fe.total_dofs, 0);
-       cell->get_dof_values (solution, dof_values);
-
-                                        // get a list of the gradients of
-                                        // the ansatz functions on this
-                                        // cell at the quadrature points
-       const vector<vector<Point<dim> > > &shape_grads(fe_face_values_cell.
-                                                       get_shape_grads());
-
-                                        // compute the gradients of the solution
-                                        // at the quadrature points by summing
-                                        // over the ansatz functions.
-       for (unsigned int j=0; j<n_q_points; ++j) 
-         for (unsigned int i=0; i<n_dofs; ++i)
-           psi[j] += dof_values[i]*shape_grads[i][j];
-
-
        
                                         // now compute over the other side of
                                         // the face
@@ -147,21 +126,17 @@ void KellyErrorEstimator<dim>::estimate_error (const DoFHandler<dim>    &dof,
            fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor,
                                            fe, boundary);
 
-                                            // get a list of the values of
-                                            // the solution for the ansatz
-                                            // functions on this neighbor
-           neighbor->get_dof_values (solution, dof_values);
-                                            // get a list of the gradients of the
-                                            // 
-           const vector<vector<Point<dim> > > &neighbor_grads (fe_face_values_cell.
-                                                               get_shape_grads());
-                                            // subtract the gradients of the
-                                            // solution on the neigbor cell
-                                            // at the quadrature points from
-                                            // those of the present cell
-           for (unsigned int j=0; j<n_q_points; ++j) 
-             for (unsigned int i=0; i<n_dofs; ++i)
-               psi[j] -= dof_values[i]*neighbor_grads[i][j];       
+                                            // get a list of the gradients of
+                                            // the finite element solution
+                                            // restricted to the neighbor cell
+           vector<Point<dim> > neighbor_psi (n_q_points);
+           fe_face_values_neighbor.get_function_grads (solution, neighbor_psi);
+
+                                            // compute the jump in the gradients
+           transform (psi.begin(), psi.end(),
+                      neighbor_psi.begin(),
+                      psi.begin(),
+                      minus<Point<dim> >());
          };
 
 
@@ -188,6 +163,7 @@ void KellyErrorEstimator<dim>::estimate_error (const DoFHandler<dim>    &dof,
        vector<double> phi(n_q_points,0);
        const vector<Point<dim> > &normal_vectors(fe_face_values_cell.
                                                  get_normal_vectors());
+       
        for (unsigned int point=0; point<n_q_points; ++point)
          phi[point] = psi[point]*normal_vectors[point];
        

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.