]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Move the stuff that computes no-flux constraints into the library proper.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 24 Jan 2008 23:28:37 +0000 (23:28 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 24 Jan 2008 23:28:37 +0000 (23:28 +0000)
git-svn-id: https://svn.dealii.org/trunk@15679 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-22/step-22.cc

index c635db8287b80b8fe1c43f5d4fb3775951df1e3b..fd58cac8eaffc63d43ca65d76f5aec81a081cfc2 100644 (file)
@@ -231,479 +231,6 @@ RightHandSide<dim>::vector_value (const Point<dim> &p,
 
 
 
-namespace internals
-{
-  namespace VectorTools
-  {
-                                    /**
-                                     * A structure that stores the dim DoF
-                                     * indices that correspond to a
-                                     * vector-valued quantity at a single
-                                     * support point.
-                                     */
-    template <int dim>
-    struct VectorDoFTuple
-    {
-       unsigned int dof_indices[dim];
-
-       bool operator < (const VectorDoFTuple<dim> &other) const
-         {
-           for (unsigned int i=0; i<dim; ++i)
-             if (dof_indices[i] < other.dof_indices[i])
-               return true;
-             else
-               if (dof_indices[i] > other.dof_indices[i])
-                 return false;
-           return false;
-         }
-
-       bool operator == (const VectorDoFTuple<dim> &other) const
-         {
-           for (unsigned int i=0; i<dim; ++i)
-             if (dof_indices[i] != other.dof_indices[i])
-               return false;
-
-           return true;
-         }
-
-       bool operator != (const VectorDoFTuple<dim> &other) const
-         {
-           return ! (*this == other);
-         }
-    };
-  }
-}
-
-
-
-template <int dim, template <int> class DH>
-void
-compute_no_normal_flux_constraints (const DH<dim>         &dof_handler,
-                                   const unsigned int     first_vector_component,
-                                   const std::set<unsigned char> &boundary_ids,
-                                   ConstraintMatrix      &constraints,
-                                   const Mapping<dim>    &mapping = StaticMappingQ1<dim>::mapping)
-{
-  Assert (dim > 1,
-         ExcMessage ("This function is not useful in 1d because it amounts "
-                     "to imposing Dirichlet values on the vector-valued "
-                     "quantity."));
-         
-  const FiniteElement<dim> &fe = dof_handler.get_fe();
-
-  std::vector<unsigned int> face_dofs (fe.dofs_per_face);
-  std::vector<Point<dim> >  dof_locations  (fe.dofs_per_face);
-
-                                  // have a map that stores normal vectors
-                                  // for each vector-dof tuple we want to
-                                  // constrain. since we can get at the same
-                                  // vector dof tuple more than once (for
-                                  // example if it is located at a vertex
-                                  // that we visit from all adjacent cells),
-                                  // we will want to average later on the
-                                  // normal vectors computed on different
-                                  // cells as described in the documentation
-                                  // of this function. however, we can only
-                                  // average if the contributions came from
-                                  // different cells, whereas we want to
-                                  // constrain twice or more in case the
-                                  // contributions came from different faces
-                                  // of the same cell. consequently, we also
-                                  // have to store which cell a normal vector
-                                  // was computed on
-  typedef 
-    std::multimap< ::internals::VectorTools::VectorDoFTuple<dim>,
-    std::pair<Tensor<1,dim>, typename DH<dim>::active_cell_iterator> >
-    DoFToNormalsMap;
-
-  DoFToNormalsMap dof_to_normals_map;
-
-                                  // now loop over all cells and all faces
-  typename DH<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-  for (; cell!=endc; ++cell)
-    for (unsigned int face_no=0; face_no < GeometryInfo<dim>::faces_per_cell;
-        ++face_no)
-      if (boundary_ids.find(cell->face(face_no)->boundary_indicator())
-         != boundary_ids.end())
-       {
-         typename DH<dim>::face_iterator face = cell->face(face_no);
-
-                                          // get the indices of the
-                                          // dofs on this cell...
-         face->get_dof_indices (face_dofs, cell->active_fe_index());
-
-                                          // ...and the normal
-                                          // vectors at the locations
-                                          // where they are defined:
-         const std::vector<Point<dim-1> > &
-           unit_support_points = fe.get_unit_face_support_points();  
-         Quadrature<dim-1> aux_quad (unit_support_points);
-         FEFaceValues<dim> fe_values (mapping, fe, aux_quad,
-                                      update_normal_vectors);
-         fe_values.reinit(cell, face_no);
-
-                                          // then identify which of
-                                          // them correspond to the
-                                          // selected set of vector
-                                          // components
-         for (unsigned int i=0; i<face_dofs.size(); ++i)
-           if (fe.face_system_to_component_index(i).first ==
-               first_vector_component)
-             {
-                                                // find corresponding other
-                                                // components of vector
-               ::internals::VectorTools::VectorDoFTuple<dim> vector_dofs;
-               vector_dofs.dof_indices[0] = face_dofs[i];
-               
-               for (unsigned int k=0; k<fe.dofs_per_face; ++k)
-                 if ((k != i)
-                     &&
-                     (unit_support_points[k] == unit_support_points[i])
-                     &&
-                     (fe.face_system_to_component_index(k).first >=
-                     first_vector_component)
-                     &&
-                     (fe.face_system_to_component_index(k).first <
-                      first_vector_component + dim))
-                   vector_dofs.dof_indices[fe.face_system_to_component_index(k).first]
-                     = face_dofs[k];
-
-                                                // and enter the
-                                                // (dofs,(normal_vector,cell))
-                                                // entry into the map
-               dof_to_normals_map
-                 .insert (std::make_pair (vector_dofs,
-                                          std::make_pair (fe_values.normal_vector(i),
-                                                          cell)));
-             }
-       }
-
-                                  // Now do something with the
-                                  // collected information. To this
-                                  // end, loop through all sets of
-                                  // pairs (dofs,normal_vector) and
-                                  // identify which entries belong to
-                                  // the same set of dofs and then do
-                                  // as described in the
-                                  // documentation, i.e. either
-                                  // average the normal vector or
-                                  // don't for this particular set of
-                                  // dofs
-  typename DoFToNormalsMap::const_iterator
-    p = dof_to_normals_map.begin();
-  
-  while (p != dof_to_normals_map.end())
-    {
-                                      // first find the range of entries in
-                                      // the multimap that corresponds to the
-                                      // same vector-dof tuple. as usual, we
-                                      // define the range half-open. the
-                                      // first entry of course is 'p'
-      typename DoFToNormalsMap::const_iterator same_dof_range[2]
-       = { p };
-      for (++p; p != dof_to_normals_map.end(); ++p)
-       if (p->first != same_dof_range[0]->first)
-         {
-           same_dof_range[1] = p;
-           break;
-         }
-      if (p == dof_to_normals_map.end())
-       same_dof_range[1] = dof_to_normals_map.end();
-
-                                      // now compute the reverse mapping: for
-                                      // each of the cells that contributed
-                                      // to the current set of vector dofs,
-                                      // add up the normal vectors. the
-                                      // values of the map are pairs of
-                                      // normal vectors and number of cells
-                                      // that have contributed
-      typedef
-       std::map
-       <typename DH<dim>::active_cell_iterator,
-       std::pair<Tensor<1,dim>, unsigned int> >
-       CellToNormalsMap;
-
-      CellToNormalsMap cell_to_normals_map;
-      for (typename DoFToNormalsMap::const_iterator
-            q = same_dof_range[0];
-          q != same_dof_range[1]; ++q)
-       if (cell_to_normals_map.find (q->second.second)
-           == cell_to_normals_map.end())
-         cell_to_normals_map[q->second.second]
-           = std::make_pair (q->second.first, 1U);
-       else
-         {
-           const Tensor<1,dim> old_normal
-             = cell_to_normals_map[q->second.second].first;
-           const unsigned int old_count
-             = cell_to_normals_map[q->second.second].second;
-
-           Assert (old_count > 0, ExcInternalError());
-
-                                            // in the same entry,
-                                            // store again the now
-                                            // averaged normal vector
-                                            // and the new count
-           cell_to_normals_map[q->second.second]
-             = std::make_pair ((old_normal * old_count + q->second.first) / (old_count + 1),
-                               old_count + 1);
-         }
-
-      Assert (cell_to_normals_map.size() >= 1, ExcInternalError());
-
-                                      // count the maximum number of
-                                      // contributions from each cell
-      unsigned int max_n_contributions_per_cell = 1;
-      for (typename CellToNormalsMap::const_iterator
-            x = cell_to_normals_map.begin();
-          x != cell_to_normals_map.end(); ++x)
-       max_n_contributions_per_cell
-         = std::max (max_n_contributions_per_cell,
-                     x->second.second);
-
-                                      // verify that each cell can have only
-                                      // contributed at most dim times, since
-                                      // that is the maximum number of faces
-                                      // that come together at a single place
-      Assert (max_n_contributions_per_cell <= dim, ExcInternalError());
-
-      switch (max_n_contributions_per_cell)
-       {
-                                          // first deal with the case that a
-                                          // number of cells all have
-                                          // registered that they have a
-                                          // normal vector defined at the
-                                          // location of a given vector dof,
-                                          // and that each of them have
-                                          // encountered this vector dof
-                                          // exactly once while looping over
-                                          // all their faces. as stated in
-                                          // the documentation, this is the
-                                          // case where we want to simply
-                                          // average over all normal vectors
-         case 1:
-         {
-           
-                                            // compute the average
-                                            // normal vector from all
-                                            // the ones that have the
-                                            // same set of dofs. we
-                                            // could add them up and
-                                            // divide them by the
-                                            // number of additions,
-                                            // or simply normalize
-                                            // them right away since
-                                            // we want them to have
-                                            // unit length anyway
-           Tensor<1,dim> normal;
-           for (typename CellToNormalsMap::const_iterator
-                  x = cell_to_normals_map.begin();
-                x != cell_to_normals_map.end(); ++x)
-             normal += x->second.first;
-           normal /= normal.norm();
-
-           const ::internals::VectorTools::VectorDoFTuple<dim> &
-               dof_indices = same_dof_range[0]->first;
-           
-           if (cell_to_normals_map.size() == 2)
-             {
-               std::cout << "XX " << same_dof_range[0]->first.dof_indices[0]
-                         << ' ' <<  same_dof_range[0]->first.dof_indices[1]
-                         << std::endl;
-               std::cout << "   " << cell_to_normals_map.begin()->first
-                         << ' '   << cell_to_normals_map.begin()->second.first
-                         << std::endl;
-               std::cout << "   " << (++cell_to_normals_map.begin())->first
-                         << ' '   << (++cell_to_normals_map.begin())->second.first
-                         << std::endl;
-               std::cout << "   " << normal << std::endl;
-             }
-           if (cell_to_normals_map.size() == 1)
-             {
-               std::cout << "YY " << same_dof_range[0]->first.dof_indices[0]
-                         << ' ' <<  same_dof_range[0]->first.dof_indices[1]
-                         << std::endl;
-               std::cout << "   " << cell_to_normals_map.begin()->first
-                         << ' '   << cell_to_normals_map.begin()->second.first
-                         << std::endl;
-               std::cout << "   " << normal << std::endl;
-             }
-           
-           
-                                            // then construct constraints
-                                            // from this. choose the DoF that
-                                            // has the largest component in
-                                            // the normal vector as the one
-                                            // to be constrained as this
-                                            // makes the process stable in
-                                            // cases where the normal vector
-                                            // has the form n=(1,0) or
-                                            // n=(0,1)
-           switch (dim)
-             {
-               case 2:
-               {
-                 if (std::fabs(normal[0]) > std::fabs(normal[1]))
-                   {
-                     constraints.add_line (dof_indices.dof_indices[0]);
-                     constraints.add_entry (dof_indices.dof_indices[0],
-                                            dof_indices.dof_indices[1],
-                                            -normal[1]/normal[0]);
-                   }
-                 else
-                   {
-                     constraints.add_line (dof_indices.dof_indices[1]);
-                     constraints.add_entry (dof_indices.dof_indices[1],
-                                            dof_indices.dof_indices[0],
-                                            -normal[0]/normal[1]);
-                   }
-                 break;
-               }
-
-               case 3:
-               {
-                 if ((std::fabs(normal[0]) >= std::fabs(normal[1]))
-                     &&
-                     (std::fabs(normal[0]) >= std::fabs(normal[2])))
-                   {
-                     constraints.add_line (dof_indices.dof_indices[0]);
-                     constraints.add_entry (dof_indices.dof_indices[0],
-                                            dof_indices.dof_indices[1],
-                                            -normal[1]/normal[0]);
-                     constraints.add_entry (dof_indices.dof_indices[0],
-                                            dof_indices.dof_indices[2],
-                                            -normal[2]/normal[0]);
-                   }
-                 else
-                   if ((std::fabs(normal[1]) >= std::fabs(normal[0]))
-                       &&
-                       (std::fabs(normal[1]) >= std::fabs(normal[2])))
-                     {
-                       constraints.add_line (dof_indices.dof_indices[1]);
-                       constraints.add_entry (dof_indices.dof_indices[1],
-                                              dof_indices.dof_indices[0],
-                                              -normal[0]/normal[1]);
-                       constraints.add_entry (dof_indices.dof_indices[1],
-                                              dof_indices.dof_indices[2],
-                                              -normal[2]/normal[1]);
-                     }
-                   else
-                     {
-                       constraints.add_line (dof_indices.dof_indices[2]);
-                       constraints.add_entry (dof_indices.dof_indices[2],
-                                              dof_indices.dof_indices[0],
-                                              -normal[0]/normal[2]);
-                       constraints.add_entry (dof_indices.dof_indices[2],
-                                              dof_indices.dof_indices[1],
-                                              -normal[1]/normal[2]);
-                     }
-           
-                 break;
-               }
-
-               default:
-                     Assert (false, ExcNotImplemented());
-             }
-
-           break;
-         }
-
-
-                                          // this is the slightly
-                                          // more complicated case
-                                          // that a single cell has
-                                          // contributed with exactly
-                                          // DIM normal vectors to
-                                          // the same set of vector
-                                          // dofs. this is what
-                                          // happens in a corner in
-                                          // 2d and 3d (but not on an
-                                          // edge in 3d, where we
-                                          // have only 2, i.e. <DIM,
-                                          // contributions. Here we
-                                          // do not want to average
-                                          // the normal
-                                          // vectors. Since we have
-                                          // DIM contributions, let's
-                                          // assume (and verify) that
-                                          // they are in fact all
-                                          // linearly independent; in
-                                          // that case, all vector
-                                          // components are
-                                          // constrained and we need
-                                          // to set them to zero
-         case dim:
-         {
-                                            // assert that indeed
-                                            // only a single cell has
-                                            // contributed
-           Assert (cell_to_normals_map.size() == 1,
-                   ExcInternalError());
-
-                                            // check linear
-                                            // independence by
-                                            // computing the
-                                            // determinant of the
-                                            // matrix created from
-                                            // all the normal
-                                            // vectors. if they are
-                                            // linearly independent,
-                                            // then the determinant
-                                            // is nonzero. if they
-                                            // are orthogonal, then
-                                            // the matrix is in fact
-                                            // equal to 1 (since they
-                                            // are all unit vectors);
-                                            // make sure the
-                                            // determinant is larger
-                                            // than 1e-3 to avoid
-                                            // cases where cells are
-                                            // degenerate
-           {
-             Tensor<2,dim> t;
-
-             typename DoFToNormalsMap::const_iterator x = same_dof_range[0];
-             for (unsigned int i=0; i<dim; ++i, ++x)
-               for (unsigned int j=0; j<dim; ++j)
-                 t[i][j] = x->second.first[j];
-             
-             Assert (std::fabs(determinant (t)) > 1e-3,
-                     ExcMessage("Found a set of normal vectors that are nearly collinear."));
-           }
-           
-                                            // so all components of
-                                            // this vector dof are
-                                            // constrained. enter
-                                            // this into the
-                                            // constraint matrix
-           for (unsigned int i=0; i<dim; ++i)
-             {
-               constraints.add_line (same_dof_range[0]->first.dof_indices[i]);
-                                                // no add_entries here
-             }
-           
-           std::cout << "DIM contributions at "
-                     << same_dof_range[0]->first.dof_indices[0]
-                     << ' '
-                     << same_dof_range[0]->first.dof_indices[1]
-                     << std::endl;
-           break;
-         }
-         
-
-                                          // this is the case of an
-                                          // edge contribution in 3d
-         default:
-               Assert (dim >= 3, ExcInternalError());
-               Assert (false, ExcNotImplemented());
-       }
-    }
-}
-
-
-
 
 
 
@@ -862,9 +389,9 @@ void InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double>       &dst,
     {
       cg.solve (*matrix, dst, src, preconditioner);
     }
-  catch (...)
+  catch (std::exception &e)
     {
-      Assert (false, ExcInternalError());
+      Assert (false, ExcMessage(e.what()));
     }
 }
 
@@ -939,8 +466,9 @@ void BoussinesqFlowProblem<dim>::setup_dofs (const bool setup_matrices)
                                           hanging_node_constraints);
   std::set<unsigned char> no_normal_flux_boundaries;
   no_normal_flux_boundaries.insert (0);
-  compute_no_normal_flux_constraints (dof_handler, 0, no_normal_flux_boundaries,
-                                     hanging_node_constraints);
+  VectorTools::compute_no_normal_flux_constraints (dof_handler, 0,
+                                                  no_normal_flux_boundaries,
+                                                  hanging_node_constraints);
   hanging_node_constraints.close ();
 
   std::vector<unsigned int> dofs_per_component (dim+2);
@@ -1727,9 +1255,14 @@ void BoussinesqFlowProblem<dim>::run ()
     {
       case 2:
       {
+//     GridGenerator::hyper_ball (triangulation);
+
+//     static const HyperBallBoundary<dim> boundary;
+//     triangulation.set_boundary (0, boundary);
+
        GridGenerator::hyper_cube (triangulation);
        
-       triangulation.refine_global (5);
+       triangulation.refine_global (6);
 
        break;
       }

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.