]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implemented support for vector valued FEs in some methods for
authoroliver <oliver@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 3 Jun 2005 15:02:14 +0000 (15:02 +0000)
committeroliver <oliver@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 3 Jun 2005 15:02:14 +0000 (15:02 +0000)
constructing RHS vectors and matrices.

git-svn-id: https://svn.dealii.org/trunk@10831 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/numerics/matrices.cc
deal.II/deal.II/source/numerics/vectors.cc

index 2e5bba699b18799044598c90640efce773215f8e..5660fb286aee4b55b72e65686ba87d2fb39aa1d8 100644 (file)
@@ -144,6 +144,7 @@ void MatrixCreator::create_mass_matrix_1 (const Mapping<dim>       &mapping,
        {
          if (coefficient->n_components==1)
            {
+             // Version for variable coefficient with 1 component
              coefficient->value_list (fe_values.get_quadrature_points(),
                                       coefficient_values);
              for (unsigned int point=0; point<n_q_points; ++point)
@@ -169,54 +170,122 @@ void MatrixCreator::create_mass_matrix_1 (const Mapping<dim>       &mapping,
            {
              coefficient->vector_value_list (fe_values.get_quadrature_points(),
                                              coefficient_vector_values);
+             if (fe.is_primitive ())
+               {
+                 // Version for variable coefficient with multiple components
+                 for (unsigned int point=0; point<n_q_points; ++point)
+                   {
+                     const double weight = fe_values.JxW(point);
+                     for (unsigned int i=0; i<dofs_per_cell; ++i)
+                       {
+                         const double v = fe_values.shape_value(i,point);
+                         const unsigned int component_i=
+                           fe.system_to_component_index(i).first;
+                         for (unsigned int j=0; j<dofs_per_cell; ++j)
+                           {
+                             const double u = fe_values.shape_value(j,point);
+                             if ((n_components==1) ||
+                                 (fe.system_to_component_index(j).first == component_i))
+                               cell_matrix(i,j) +=
+                                 (u * v * weight *
+                                  coefficient_vector_values[point](component_i));
+                           }
+                       }
+                   }
+               }
+             else
+               {
+                 // Version for variable coefficient with multiple components and
+                 // vector values FE.
+                 for (unsigned int point=0; point<n_q_points; ++point)
+                   {
+                     const double weight = fe_values.JxW(point);
+                     for (unsigned int i=0; i<dofs_per_cell; ++i)
+                       for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+                         if (fe.get_nonzero_components(i)[comp_i])
+                           {
+                             const double v = fe_values.shape_value_component(i,point,comp_i);
+                             for (unsigned int j=0; j<dofs_per_cell; ++j)
+                               for (unsigned int comp_j = 0; comp_j < n_components; ++comp_j)
+                                 if (fe.get_nonzero_components(j)[comp_j])
+                                   {
+                                     const double u = fe_values.shape_value_component(j,point,comp_j);
+                                     if (comp_i == comp_j)
+                                       cell_matrix(i,j) +=
+                                         (u * v * weight *
+                                          coefficient_vector_values[point](comp_i));
+                                   }
+                           }
+                   }
+               }
+           }
+       }
+      else
+       {
+         if (fe.is_primitive ())
+           {
+             // Version for primitive FEs
              for (unsigned int point=0; point<n_q_points; ++point)
                {
                  const double weight = fe_values.JxW(point);
                  for (unsigned int i=0; i<dofs_per_cell; ++i)
                    {
-                     const double v = fe_values.shape_value(i,point);
-                     const unsigned int component_i=
-                       fe.system_to_component_index(i).first;
+                     const double v = fe_values.shape_value(i,point); 
                      for (unsigned int j=0; j<dofs_per_cell; ++j)
                        {
                          const double u = fe_values.shape_value(j,point);
                          if ((n_components==1) ||
-                             (fe.system_to_component_index(j).first == component_i))
-                           cell_matrix(i,j) +=
-                             (u * v * weight *
-                              coefficient_vector_values[point](component_i));
+                             (fe.system_to_component_index(i).first ==
+                              fe.system_to_component_index(j).first))
+                           cell_matrix(i,j) += (u * v * weight);
                        }
                    }
                }
            }
+         else
+           {
+             // Version for vector valued FEs
+             for (unsigned int point=0; point<n_q_points; ++point)
+               {
+                 const double weight = fe_values.JxW(point);
+                 for (unsigned int i=0; i<dofs_per_cell; ++i)
+                   for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+                     if (fe.get_nonzero_components(i)[comp_i])
+                       {
+                         const double v = fe_values.shape_value_component(i,point,comp_i); 
+                         for (unsigned int j=0; j<dofs_per_cell; ++j)
+                           for (unsigned int comp_j = 0; comp_j < n_components; ++comp_j)
+                             if (fe.get_nonzero_components(j)[comp_j])
+                               {
+                                 const double u = fe_values.shape_value_component(j,point,comp_j);
+                                 if (comp_i == comp_j)
+                                   cell_matrix(i,j) += (u * v * weight);
+                               }
+                       }
+               }
+           }
        }
-      else
-       for (unsigned int point=0; point<n_q_points; ++point)
-         {
-           const double weight = fe_values.JxW(point);
-           for (unsigned int i=0; i<dofs_per_cell; ++i)
-             {
-               const double v = fe_values.shape_value(i,point); 
-               for (unsigned int j=0; j<dofs_per_cell; ++j)
-                 {
-                   const double u = fe_values.shape_value(j,point);
-                   if ((n_components==1) ||
-                       (fe.system_to_component_index(i).first ==
-                        fe.system_to_component_index(j).first))
-                     cell_matrix(i,j) += (u * v * weight);
-                 }
-             }
-         }
                                       // transfer everything into the
                                       // global object
       mutex.acquire ();
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       for (unsigned int j=0; j<dofs_per_cell; ++j)
-         if ((n_components==1) ||
-             (fe.system_to_component_index(i).first ==
-              fe.system_to_component_index(j).first))
-           matrix.add (dof_indices[i], dof_indices[j],
-                       cell_matrix(i,j));
+      if (fe.is_primitive ())
+       {
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           for (unsigned int j=0; j<dofs_per_cell; ++j)
+             if ((n_components==1) ||
+                 (fe.system_to_component_index(i).first ==
+                  fe.system_to_component_index(j).first))
+               matrix.add (dof_indices[i], dof_indices[j],
+                           cell_matrix(i,j));
+       }
+      else
+       {
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           for (unsigned int j=0; j<dofs_per_cell; ++j)
+             matrix.add (dof_indices[i], dof_indices[j],
+                         cell_matrix(i,j));
+       }
+
       mutex.release ();
     };
 }
index 1a5775defb11e2a7e617f62557a27403566cd78b..2bbad617ac37e6b54752c3aff7abb1a994504dfa 100644 (file)
@@ -493,29 +493,60 @@ void VectorTools::create_right_hand_side (const Mapping<dim>    &mapping,
     {
       std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
       
-      for (; cell!=endc; ++cell) 
+      // Use the faster code if the FiniteElement is primitive
+      if (fe.is_primitive ())
        {
-         fe_values.reinit(cell);
-         
-         const std::vector<double> &weights   = fe_values.get_JxW_values ();
-         rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
-         
-         cell_vector = 0;
-         for (unsigned int point=0; point<n_q_points; ++point)
-           for (unsigned int i=0; i<dofs_per_cell; ++i)
-             {
-               const unsigned int component
-                 = fe.system_to_component_index(i).first;
-               
-               cell_vector(i) += rhs_values[point](component) *
-                                 fe_values.shape_value(i,point) *
-                                 weights[point];
-             }
-         
-         cell->get_dof_indices (dofs);
-         
-         for (unsigned int i=0; i<dofs_per_cell; ++i)
-           rhs_vector(dofs[i]) += cell_vector(i);
+         for (; cell!=endc; ++cell) 
+           {
+             fe_values.reinit(cell);
+             
+             const std::vector<double> &weights   = fe_values.get_JxW_values ();
+             rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
+             
+             cell_vector = 0;
+             for (unsigned int point=0; point<n_q_points; ++point)
+               for (unsigned int i=0; i<dofs_per_cell; ++i)
+                 {
+                   const unsigned int component
+                     = fe.system_to_component_index(i).first;
+                   
+                   cell_vector(i) += rhs_values[point](component) *
+                                     fe_values.shape_value(i,point) *
+                                     weights[point];
+                 }
+             
+             cell->get_dof_indices (dofs);
+             
+             for (unsigned int i=0; i<dofs_per_cell; ++i)
+               rhs_vector(dofs[i]) += cell_vector(i);
+           }
+       }
+      else
+       // Otherwise do it the way proposed for vector valued elements
+       {
+         for (; cell!=endc; ++cell) 
+           {
+             fe_values.reinit(cell);
+             
+             const std::vector<double> &weights   = fe_values.get_JxW_values ();
+             rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
+             
+             cell_vector = 0;
+             for (unsigned int point=0; point<n_q_points; ++point)
+               for (unsigned int i=0; i<dofs_per_cell; ++i)
+                 for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+                   if (fe.get_nonzero_components(i)[comp_i])
+                     {
+                       cell_vector(i) += rhs_values[point](comp_i) *
+                                         fe_values.shape_value_component(i,point,comp_i) *
+                                         weights[point];
+                     }
+             
+             cell->get_dof_indices (dofs);
+             
+             for (unsigned int i=0; i<dofs_per_cell; ++i)
+               rhs_vector(dofs[i]) += cell_vector(i);
+           }
        }
     }
 }
@@ -610,22 +641,40 @@ VectorTools::create_boundary_right_hand_side (const Mapping<dim>      &mapping,
               boundary_indicators.end()))
            {
              fe_values.reinit(cell, face);
-         
+             
              const std::vector<double> &weights   = fe_values.get_JxW_values ();
              rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
-         
+             
              cell_vector = 0;
-             for (unsigned int point=0; point<n_q_points; ++point)
-               for (unsigned int i=0; i<dofs_per_cell; ++i)
-                 {
-                   const unsigned int component
-                     = fe.system_to_component_index(i).first;
-                   
-                   cell_vector(i) += rhs_values[point](component) *
-                                     fe_values.shape_value(i,point) *
-                                     weights[point];
-                 }
              
+             // Use the faster code if the FiniteElement is primitive
+             if (fe.is_primitive ())
+               {                 
+                 for (unsigned int point=0; point<n_q_points; ++point)
+                   for (unsigned int i=0; i<dofs_per_cell; ++i)
+                     {
+                       const unsigned int component
+                         = fe.system_to_component_index(i).first;
+                       
+                       cell_vector(i) += rhs_values[point](component) *
+                                         fe_values.shape_value(i,point) *
+                                         weights[point];
+                     }
+               }
+             else
+               {
+                 // And the full featured code, if vector valued FEs are used
+                 for (unsigned int point=0; point<n_q_points; ++point)
+                   for (unsigned int i=0; i<dofs_per_cell; ++i)
+                     for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+                       if (fe.get_nonzero_components(i)[comp_i])
+                         {
+                           cell_vector(i) += rhs_values[point](comp_i) *
+                                             fe_values.shape_value_component(i,point,comp_i) *
+                                             weights[point];
+                         }
+               }
+                 
              cell->get_dof_indices (dofs);
              
              for (unsigned int i=0; i<dofs_per_cell; ++i)

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.