]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
working on mg boundary
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 4 Jul 2006 05:49:15 +0000 (05:49 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 4 Jul 2006 05:49:15 +0000 (05:49 +0000)
git-svn-id: https://svn.dealii.org/trunk@13329 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/multigrid/mg_dof_tools.cc
tests/bits/Makefile
tests/bits/dof_tools_01c.cc
tests/dofs/Makefile [new file with mode: 0644]
tests/dofs/row_length_01.cc [new file with mode: 0644]
tests/multigrid/Makefile
tests/multigrid/boundary_01.cc [new file with mode: 0644]
tests/multigrid/boundary_01/.cvsignore [new file with mode: 0644]

index 9574b494bc4dc0893e96bc2a0ec114d14802ca70..124e7b3450a5436b199ede4cd52faefe7b6c752e 100644 (file)
@@ -1017,7 +1017,205 @@ MGTools::count_dofs_per_component (
 }
 
 
+#if deal_II_dimension == 1
+
+template <int dim>
+void
+MGTools::make_boundary_list(
+  const MGDoFHandler<dim>&,
+  const typename FunctionMap<dim>::type&,
+  std::vector<std::vector<unsigned int> >&,
+  const std::vector<bool>&)
+{
+  Assert(false, ExcNotImplemented());
+}
 
+
+#else
+
+template <int dim>
+void
+MGTools::make_boundary_list(
+  const MGDoFHandler<dim>& dof,
+  const typename FunctionMap<dim>::type& function_map,
+  std::vector<std::vector<unsigned int> >& boundary_indices,
+  const std::vector<bool>& component_mask_)
+{
+                                   // if for whatever reason we were
+                                   // passed an empty map, return
+                                   // immediately
+  if (function_map.size() == 0)
+    return;
+  
+  const unsigned int        n_components = DoFTools::n_components(dof);
+  const bool                fe_is_system = (n_components != 1);
+
+//    for (typename FunctionMap<dim>::type::const_iterator i=function_map.begin();
+//         i!=function_map.end(); ++i)
+//      Assert (n_components == i->second->n_components,
+//           ExcInvalidFE());
+
+                                   // set the component mask to either
+                                   // the original value or a vector
+                                   // of @p{true}s
+  const std::vector<bool> component_mask ((component_mask_.size() == 0) ?
+                                          std::vector<bool> (n_components, true) :
+                                          component_mask_);
+//   Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
+//        ExcComponentMismatch());
+
+                                   // field to store the indices
+  std::vector<unsigned int> face_dofs;
+  face_dofs.reserve (DoFTools::max_dofs_per_face(dof));
+  std::fill (face_dofs.begin (), face_dofs.end (), DoFHandler<dim>::invalid_dof_index);
+  
+  typename MGDoFHandler<dim>::cell_iterator cell = dof.begin(),
+                                           endc = dof.end();
+  for (; cell!=endc; ++cell)
+    for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
+         ++face_no)
+      {
+        const FiniteElement<dim> &fe = cell->get_fe();
+        const unsigned int level = cell->level();
+
+                                         // we can presently deal only with
+                                         // primitive elements for boundary
+                                         // values. this does not preclude
+                                         // us using non-primitive elements
+                                         // in components that we aren't
+                                         // interested in, however. make
+                                         // sure that all shape functions
+                                         // that are non-zero for the
+                                         // components we are interested in,
+                                         // are in fact primitive
+        for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+          {
+            const std::vector<bool> &nonzero_component_array
+              = cell->get_fe().get_nonzero_components (i);
+            for (unsigned int c=0; c<n_components; ++c)
+              if ((nonzero_component_array[c] == true)
+                  &&
+                  (component_mask[c] == true))
+                Assert (cell->get_fe().is_primitive (i),
+                        ExcMessage ("This function can only deal with requested boundary "
+                                    "values that correspond to primitive (scalar) base "
+                                    "elements"));
+          }
+
+        typename MGDoFHandler<dim>::face_iterator face = cell->face(face_no);
+        const unsigned char boundary_component = face->boundary_indicator();
+        if (function_map.find(boundary_component) != function_map.end())
+                                           // face is of the right component
+          {
+                                             // get indices, physical location and
+                                             // boundary values of dofs on this
+                                             // face
+            face_dofs.resize (fe.dofs_per_face);
+            face->get_mg_dof_indices (face_dofs);
+            if (fe_is_system)
+              {
+                                                 // enter those dofs
+                                                 // into the list that
+                                                 // match the
+                                                 // component
+                                                 // signature. avoid
+                                                 // the usual
+                                                 // complication that
+                                                 // we can't just use
+                                                 // *_system_to_component_index
+                                                 // for non-primitive
+                                                 // FEs
+                for (unsigned int i=0; i<face_dofs.size(); ++i)
+                  {
+                    unsigned int component;
+                    if (fe.is_primitive())
+                      component = fe.face_system_to_component_index(i).first;
+                    else
+                      {
+                                                         // non-primitive
+                                                         // case. make
+                                                         // sure that
+                                                         // this
+                                                         // particular
+                                                         // shape
+                                                         // function
+                                                         // _is_
+                                                         // primitive,
+                                                         // and get at
+                                                         // it's
+                                                         // component. use
+                                                         // usual
+                                                         // trick to
+                                                         // transfer
+                                                         // face dof
+                                                         // index to
+                                                         // cell dof
+
+                                                         // index
+                        const unsigned int cell_i
+                          = (dim == 1 ?
+                             i
+                             :
+                             (dim == 2 ?
+                              (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
+                              :
+                              (dim == 3 ?
+                               (i<4*fe.dofs_per_vertex ?
+                                i
+                                :
+                                (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
+                                 i+4*fe.dofs_per_vertex
+                                 :
+                                 i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
+                               :
+                               deal_II_numbers::invalid_unsigned_int)));
+                        Assert (cell_i < fe.dofs_per_cell, ExcInternalError());
+
+                                                         // make sure
+                                                         // that if
+                                                         // this is
+                                                         // not a
+                                                         // primitive
+                                                         // shape function,
+                                                         // then all
+                                                         // the
+                                                         // corresponding
+                                                         // components
+                                                         // in the
+                                                         // mask are
+                                                         // not set
+//                         if (!fe.is_primitive(cell_i))
+//                           for (unsigned int c=0; c<n_components; ++c)
+//                             if (fe.get_nonzero_components(cell_i)[c])
+//                               Assert (component_mask[c] == false,
+//                                       ExcFENotPrimitive());
+
+// let's pick the first of possibly more than one non-zero
+// components. if shape function is non-primitive, then we will ignore
+// the result in the following anyway, otherwise there's only one
+// non-zero component which we will use
+                        component = (std::find (fe.get_nonzero_components(cell_i).begin(),
+                                                fe.get_nonzero_components(cell_i).end(),
+                                                true)
+                                     -
+                                     fe.get_nonzero_components(cell_i).begin());
+                      }
+
+                    if (component_mask[component] == true)
+                      boundary_indices[level].push_back(face_dofs[i]);
+                  }
+              }
+            else
+             for (unsigned int i=0; i<face_dofs.size(); ++i)
+               boundary_indices[level].push_back(face_dofs[i]);
+          }
+      }
+}
+
+#endif
+
+
+  
 template<int dim, typename number>
 void
 MGTools::reinit_vector (const MGDoFHandler<dim> &mg_dof,
@@ -1218,3 +1416,9 @@ template void MGTools::count_dofs_per_component<deal_II_dimension> (
   const MGDoFHandler<deal_II_dimension>&,
   std::vector<std::vector<unsigned int> >&,
   std::vector<unsigned int>);
+
+template void MGTools::make_boundary_list(
+  const MGDoFHandler<deal_II_dimension>&,
+  const FunctionMap<deal_II_dimension>::type&,
+  std::vector<std::vector<unsigned int> >&,
+  const std::vector<bool>&);
index e0522f2222988866ddb610f3d89dd7274e87acee..52042bb4568cdd4e11059d8f71404edb0cc19acc 100644 (file)
@@ -32,6 +32,7 @@ tests_x = geometry_info_* \
          data_out_rotation_0* \
          data_out_stack_0* \
          dof_tools_[0-9]* \
+         dof_tools_row_[0-9] \
          fe_tools_[0-9]* \
          fe_tools_cpfqpm* \
           anna_? \
index 7a402f7ccd8651968c6d7e9c9bba2b416143d5bb..a6e3d1d6045b77eea092d134be1e5d097963da1f 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$ 
 //
-//    Copyright (C) 2003, 2004 by the deal.II authors
+//    Copyright (C) 2003, 2004, 2006 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -38,17 +38,22 @@ check_this (const DoFHandler<dim> &dof_handler)
   
                                    // create sparsity pattern
   const unsigned int n_components = dof_handler.get_fe().n_components();
-  BlockSparsityPattern sp (n_components,
-                           n_components);
-  std::vector<unsigned int> dofs_per_component(n_components);
-  DoFTools::count_dofs_per_component (dof_handler,
-                                      dofs_per_component);
-  for (unsigned int i=0; i<n_components; ++i)
-    for (unsigned int j=0; j<n_components; ++j)
-      sp.block(i,j).reinit(dofs_per_component[i],
-                           dofs_per_component[j],
-                           dof_handler.max_couplings_between_dofs());
-  sp.collect_sizes ();
+  const unsigned int n_blocks = dof_handler.get_fe().n_blocks();
+  
+  BlockSparsityPattern sp (n_blocks,
+                           n_blocks);
+  std::vector<unsigned int> dofs_per_component(n_blocks);
+  DoFTools::count_dofs_per_block (dof_handler,
+                                 dofs_per_component);
+  BlockIndices block_indices(dofs_per_component);
+  
+  std::vector<std::vector<unsigned int> >
+    row_length(n_blocks, std::vector<unsigned int>(dof_handler.n_dofs()));
+  Table<2,DoFTools::Coupling> couplings(n_components, n_components);
+  couplings.fill(DoFTools::always);
+  DoFTools::compute_row_length_vector(dof_handler, row_length,
+                                     couplings, couplings);
+  sp.reinit(block_indices, block_indices, row_length);
   
   DoFTools::make_sparsity_pattern (dof_handler, sp);
   sp.compress ();
diff --git a/tests/dofs/Makefile b/tests/dofs/Makefile
new file mode 100644 (file)
index 0000000..e169d15
--- /dev/null
@@ -0,0 +1,59 @@
+############################################################
+# Makefile,v 1.15 2002/06/13 12:51:13 hartmann Exp
+# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+############################################################
+
+############################################################
+# Include general settings for including DEAL libraries
+############################################################
+
+include ../Makefile.paths
+include $D/common/Make.global_options
+debug-mode = on
+
+libraries = $(lib-deal2-1d.g) \
+            $(lib-deal2-2d.g) \
+            $(lib-deal2-3d.g) \
+            $(lib-lac.g)  \
+            $(lib-base.g)
+
+default: run-tests
+
+############################################################
+
+tests_x = row_length_*
+
+# tests for the hp branch:
+#        fe_collection_*
+
+ifeq ($(USE_CONTRIB_PETSC),yes)
+  tests_x += petsc_*
+endif
+
+ifeq ($(USE_CONTRIB_METIS),yes)
+  tests_x += metis_*
+endif
+
+ifeq ($(USE_CONTRIB_HSL),yes)
+  tests_x += hsl_*
+endif
+
+ifeq ($(USE_CONTRIB_UMFPACK),yes)
+  tests_x += umfpack_*
+endif
+
+
+# from above list of regular expressions, generate the real set of
+# tests
+expand = $(shell echo $(addsuffix .cc,$(1)) \
+           | $(PERL) -pi -e 's/\.cc//g;')
+tests = $(call expand,$(tests_x))
+
+############################################################
+
+
+include ../Makefile.rules
+include Makefile.depend
+include Makefile.tests
+
+.PHONY: default
diff --git a/tests/dofs/row_length_01.cc b/tests/dofs/row_length_01.cc
new file mode 100644 (file)
index 0000000..9774b8e
--- /dev/null
@@ -0,0 +1,95 @@
+//----------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2003, 2004, 2005, 2006 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------------------------------------------------
+
+#include "../lib/dof_tools_frame.h"
+#include "../lib/test_grids.h"
+
+#include <lac/sparsity_pattern.h>
+#include <numeric>
+
+template <int dim>
+void
+check_this (const DoFHandler<dim> &dof)
+{
+  const unsigned int n = dof.n_dofs();
+  deallog << "dofs: " << n
+         << "\tcell: " << dof.get_fe().dofs_per_cell
+         << std::endl;
+  
+  std::vector<unsigned int> row_length(n);
+  DoFTools::compute_row_length_vector(dof, row_length);
+  SparsityPattern sparsity(n, n, row_length);
+  DoFTools::make_sparsity_pattern(dof, sparsity);
+  sparsity.compress();
+  unsigned int sum = std::accumulate(row_length.begin(), row_length.end(), 0U);
+  deallog << std::endl << "Entries estimated/actual " << sum
+         <<  '/' << sparsity.n_nonzero_elements() << std::endl;
+//  sparsity.clear();
+  output_vector(row_length);
+  
+  DoFTools::compute_row_length_vector(dof, row_length, DoFTools::always);
+  output_vector(row_length);
+}
+
+
+template <int dim>
+void check()
+{
+  Triangulation<dim> tr;
+  TestGrids::hypercube(tr);
+  deallog << "cube" << dim << std::endl;
+  check_grid(tr);
+  tr.refine_global(1);
+  deallog << "refined cube" << dim << std::endl;
+  check_grid(tr);  
+}
+
+
+int
+main()
+{
+  try
+    {
+      std::ofstream logfile("row_length_01/output");
+      logfile.precision (2);
+      deallog.attach(logfile);
+//      deallog.depth_console(0);
+      deallog.threshold_double(1.e-10);
+      check<2>();
+      return 0;
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+}
+
index 06559e9d0aa8d83cca0cd05ba3bc2f333d00a1f8..19074c3432a384ccefc936ddff427e757976d24b 100644 (file)
@@ -21,8 +21,9 @@ default: run-tests
 
 ############################################################
 
-tests_x = cycles transfer transfer_select transfer_block \
-       smoother_block count_* renumbering_*
+tests_x = cycles count_* boundary_* renumbering_* \
+       transfer transfer_select transfer_block \
+       smoother_block
 
 # from above list of regular expressions, generate the real set of
 # tests
diff --git a/tests/multigrid/boundary_01.cc b/tests/multigrid/boundary_01.cc
new file mode 100644 (file)
index 0000000..50693e9
--- /dev/null
@@ -0,0 +1,115 @@
+//----------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2006 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------------------------------------------------------
+
+// check MGTools::make_boundary_list
+
+#include <base/logstream.h>
+#include <base/function.h>
+#include <lac/vector.h>
+#include <lac/block_vector.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_generator.h>
+#include <dofs/function_map.h>
+#include <fe/fe_dgq.h>
+#include <fe/fe_q.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/fe_system.h>
+#include <multigrid/mg_dof_accessor.h>
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_tools.h>
+
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <algorithm>
+
+using namespace std;
+
+void log_vector (const std::vector<std::vector<unsigned int> >& count)
+{
+  for (unsigned int l=0;l<count.size();++l)
+    {
+      deallog << "Level " << l;
+      for (unsigned int c=0;c<count[l].size();++c)
+       deallog << '\t' << count[l][c];
+      deallog << std::endl;
+    }
+}
+
+
+template <int dim>
+void face_levels(const MGDoFHandler<dim>& dof)
+{
+  typename MGDoFHandler<dim>::face_iterator face;
+  const typename MGDoFHandler<dim>::face_iterator end = dof.end_face();
+  
+  for (face = dof.begin_face(); face != end; ++face)
+    {
+      deallog << "Face " << face->level();
+      for (unsigned int i=0;i<GeometryInfo<dim>::vertices_per_face;++i)
+       deallog << " v" << face->vertex(i);
+      deallog << " dofs ";
+      deallog << std::endl;
+    }
+}
+
+
+template <int dim>
+void check_fe(FiniteElement<dim>& fe)
+{
+  deallog << fe.get_name() << std::endl;
+  
+  Triangulation<dim> tr;
+  GridGenerator::hyper_cube(tr);
+  tr.refine_global(2);
+  ZeroFunction<dim> zero;
+  typename FunctionMap<dim>::type fmap;
+  fmap.insert(std::make_pair(0, &zero));
+  
+  MGDoFHandler<dim> mgdof(tr);
+  mgdof.distribute_dofs(fe);
+  face_levels(mgdof);
+  
+  std::vector<std::vector<unsigned int> > boundary_indices(tr.n_levels());
+  MGTools::make_boundary_list(mgdof, fmap, boundary_indices);
+  log_vector(boundary_indices);
+}
+
+
+template <int dim>
+void check()
+{
+  FE_Q<dim> q1(1);
+  FE_Q<dim> q2(2);
+//  FE_DGQ<dim> dq1(1);
+  
+  FESystem<dim> s1(q1, 2, q2,1);
+
+  check_fe(q1);
+  check_fe(q2);
+  check_fe(s1);
+}
+
+int main()
+{
+  std::ofstream logfile("boundary_01/output");
+  logfile.precision(3);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  check<2> ();
+  check<3> ();
+}
diff --git a/tests/multigrid/boundary_01/.cvsignore b/tests/multigrid/boundary_01/.cvsignore
new file mode 100644 (file)
index 0000000..d196dfb
--- /dev/null
@@ -0,0 +1 @@
+obj.* exe OK output

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.