}
+#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,
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>&);
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_? \
// $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
// 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 ();
--- /dev/null
+############################################################
+# 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
--- /dev/null
+//----------------------------------------------------------------------
+// $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;
+ };
+}
+
############################################################
-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
--- /dev/null
+//----------------------------------------------------------------------------
+// $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> ();
+}
--- /dev/null
+obj.* exe OK output